尧图网站建设 
简介这份黑龙江省土壤水文分组高精度栅格数据集源自HYSOGs250m经省级裁剪配准面向水文建模、水土保持与SWAT模型应用人群。数据采用WGS84坐标系提供约250米分辨率的A/B/C/D四类土壤径流潜力栅格同时标注浅地下水位引发的双重HSG可直观用于降雨径流潜力分区与CN值估算。包体共7个文件以tif栅格为主配套tfw坐标参数、dbf属性表、xml元数据、pdf说明文档与cpg编码满足主流GIS平台便捷加载压缩包仅4.28MB轻量易用。已有85人学习浏览。借助该数据集可直接获得黑龙江省统一的土壤水文分组成果省去全球数据下载、拼接与裁剪流程。数据附带的属性表与元数据便于批量读取和制图适合SWAT建模前期参数准备、省域径流潜力对比及教学科研中的快速演示。1. 黑龙江省HSG栅格从默认土壤到250米径流潜力分级做SWAT模型的人大多在土壤数据上栽过跟头。默认的1千米土壤图在黑龙江这种地形破碎、三江平原湿地连片的地方跑出来的CN值不是偏大就是偏小率定过程很痛苦。HYSOGs250m把USDA曲线数方法依赖的土壤水文分组HSG放到了250米分辨率上每个像元给出A、B、C、D四档径流潜力湿润区还标了双重HSG。这个栅格基于SoilGrids250m派生已经按黑龙江省界裁剪好坐标系是WGS84适合直接进入SWAT的径流模拟流程。对于做水文建模和GIS数据预处理的人来说关键是先读懂分类逻辑再把它转成SWAT认识的格式。2. HYSOGs250m的底账A/B/C/D分类、双重HSG与栅格文件结构2.1 四档HSG如何从SoilGrids250m派生HYSOGs250m没有直接钻探土壤而是用FAO的SoilGrids250m系统给出的土壤质地等级和到基岩的深度在1/480十进制度网格上推演出一个“排水性能”分级。SCS-CN方法里土壤被分成A、B、C、D四档代表降雨入渗和径流潜力的顺序。A类对应砂质或砾质土壤导水性强降雨下去很快B类是砂壤土下渗中等C类偏黏壤土下渗受限D类则是黏土或浅层不透水基本靠地表产流。这里有一组常用的最小下渗率参考值虽然不是HYSOGs250m直接输出的字段但能帮你理解分级边界HSG典型土壤条件最小下渗率mm/h径流潜力A砂土、砾质土深厚且排水好7.6以上低B砂壤土土层较深3.8–7.6中等偏低C黏壤土、细砂壤土下渗受限1.3–3.8中等偏高D黏土、浅层土壤或有地下水位小于1.3高这个表是SCS经典公式的参数。你在SWAT里填CN值时如果CN表里给出A/B/C/D四列用的就是这套逻辑。HYSOGs250m把SoilGrids250m的连续属性离散成了这四档所以它本质上是一个分类图而不是属性图。不要把栅格值当作可平滑处理的连续数值做重采样必须用最近邻。2.2 双重HSG是怎么回事USDA的HSG分类里还有一个特殊情况如果地表以下60厘米范围内存在地下水位土壤不管表层是砂质还是黏质在雨季都表现出高径流潜力。HYSOGs250m会为这类区域赋予双重HSG比如A/D、B/D或C/D。前者是排水良好的干燥状态后者是饱和状态。黑龙江省的三江平原、挠力河流域和松嫩低洼地区经常出现这种双代码。SWAT标准的usersoil表里只有一列HYDGRP装不下双代码。我一般的做法是在雨季模拟时取保守值比如把A/D当作D使用如果在做季节率定也可以利用双代码做月尺度切换。要注意的是只看字母A就把它当成低产流在三江平原会明显低估洪峰。2.3 解压后每个文件都有用解压后拿到一堆文件初看好像冗余其实每一项都有固定角色尤其是.vat.dbf经常被忽略很多人的做法是直接读tif然后靠猜这很容易在双代码上翻车。下面先列最常见的几个文件文件作用HYSOGs250m_黑龙江省.tif主栅格像元值对应HSG代码.tfw世界文件记录影像左上角坐标和像元尺寸.aux.xmlESRI辅助元数据记录统计信息和色彩映射.vat.dbf栅格值属性表把栅格Value翻译成HSG字母.vat.cpgVAT表的字符编码文件土壤水文分组.pdf数据说明文档.vat.dbf是理解栅格值的关键。可以用GDAL自带的ogrinfo直接看ogrinfo -al HYSOGs250m_黑龙江省.tif.vat.dbf输出里通常有Value和HSG两个字段Value是tif里存的整数HSG是对应的字母或双字母代码。先跑这条命令确认代码含义再继续处理不要凭经验猜。否则后面统计面积或重分类时把A/D当成某个数字就会出偏差。如果PDF说明文档和VAT表对不上以VAT表为准因为VAT是随栅格一起发布的PDF可能描述的是整个中国范围而黑龙江省版已经做过裁剪和重编码。另外比起全国版栅格这个版本经过了省级裁剪文件体积更小但边缘像元仍可能覆盖邻省或省界外的无意义区域。后续处理建议用黑龙江行政边界再做一次裁剪避免统计面积时混入界外像元。3. GDALPython处理黑龙江省HSG栅格投影、统计与重映射现在进入实操。处理这类区域裁剪后的HSG栅格我一般分成四步看元数据、读VAT、投影后做面积统计、重新编码输出。3.1 打开栅格前先做一次体检拿到HYSOGs250m_黑龙江省.tif先跑gdalinfo HYSOGs250m_黑龙江省.tif重点看四行Size、Origin、Pixel Size和Coordinate System。这个数据集是WGS84地理坐标系像素分辨率是0.00208333度约250米。NoData也要记下来通常可能是-9999或0后面统计时要过滤。再看旁边的.tfw世界文件。一个标准的TFW文件有6行数字前两行是像元在X和Y方向上的尺寸中间两行是旋转项正常为0第五、六行是左上角像元中心坐标。用记事本打开就能确认行列范围是否覆盖黑龙江边界尤其是跨到俄罗斯一侧的区域不需要额外做范围裁切。3.2 用VAT把Value翻译成HSG代码我建议一上来就把VAT读进内存生成一个Value到HSG字母的映射。这样后续所有统计都能直接打印字母。from osgeo import ogr vat ogr.Open(HYSOGs250m_黑龙江省.tif.vat.dbf) layer vat.GetLayer(0) mapping {} while True: feat layer.GetNextFeature() if feat is None: break value feat.GetField(Value) hsg feat.GetField(HSG) mapping[int(value)] hsg print(f{value} - {hsg})这段代码遍历VAT表的每一行取出Value和HSG两个字段放进dict。后续对栅格数组做np.unique后可以用这个dict把整数值翻译成HSG字母。注意字段名不一定叫HSG有的版本叫Label或Class可以先打印layer.GetLayerDefn()确认字段列表。3.3 投影后做面积统计WGS84地理坐标系下每个像元在不同纬度对应的真实面积不一样。黑龙江省跨北纬43度到53度如果直接用250乘以250统计面积误差在5%以上。面积统计前我先转到Albers等积投影。gdalwarp -t_srs projaea lat_145 lat_260 lat_025 lon_0105 datumWGS84 \ -tr 250 250 -r near -overwrite \ HYSOGs250m_黑龙江省.tif HYSOGs250m_hlj_aea.tif该命令中-t_srs指定Albers等积圆锥投影两条标准纬线选45和60度覆盖黑龙江主体-tr 250 250把输出像元严格设为250米-r near用最近邻插值保证分类值不产生新层级。投影完成后栅格每个像元面积都是62500平方米统计才有意义。Python统计各HSG面积from osgeo import gdal import numpy as np ds gdal.Open(HYSOGs250m_hlj_aea.tif) arr ds.GetRasterBand(1).ReadAsArray() arr arr.astype(np.int16) # 过滤NoData按原始Value统计 valid_mask (arr 0) values, counts np.unique(arr[valid_mask], return_countsTrue) pixel_area_km2 250 * 250 / 1e6 for v, c in zip(values, counts): hsg mapping.get(int(v), unknown) print(f{hsg}: {c * pixel_area_km2:.2f} km2)astype(np.int16)是为了防止某些配置下栅格以无符号整型读入负号和NoData变成大数valid_mask过滤掉小于等于0的像元。如果NoData是-9999arr 0同样可以把它滤掉。统计结果出来后你会看到黑龙江各HSG类的大致比例通常C和D占大头这与黑土和草甸土分布基本吻合。3.4 重分类并导出SWAT可用的GRID统计完之后通常要重编码。比如把双代码A/D统一赋成4生成一个只含1、2、3、4的四类栅格。这样后面给SWAT的土壤图时字段更简洁。reclass np.zeros_like(arr) for v, h in mapping.items(): if D in h: reclass[arr v] 4 elif h C: reclass[arr v] 3 elif h B: reclass[arr v] 2 elif h A: reclass[arr v] 1这段代码把任何包含D的双代码都归到第4类其余A/B/C对应1/2/3。如果你做湿润季节模拟这是最稳妥的保守处理。保存重分类结果时要继承原始投影和仿射参数driver gdal.GetDriverByName(GTiff) out_ds driver.Create(reclass.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Int16) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(reclass) out_band.SetNoDataValue(0) out_ds.FlushCache()最后用gdal_translate转成文本GridSWAT可以直接读取gdal_translate -of AAIGrid -ot Int16 -a_nodata 0 reclass.tif hsg_reclass.asc这样HSG空间分布图就准备好进入SWAT流程了。注意AAIGrid文件头里的ncols和nrows必须和DEM一致如果DEM分辨率不一样需要先重采样到同一格网。4. 在SWAT模型里把HSG变成CN值从土壤属性到HRU4.1 SWAT中的HYDGRP与usersoil表SWAT不是直接用栅格里的HSG参与产流计算的。它把HSG存在土壤数据库的HYDGRP字段里土地利用和土壤类型组合成HRU后再通过CN2值进入径流曲线数计算。所以在SWAT里看到的是“每一类土壤对应一个A/B/C/D”不是一个逐像元字段。如果你的模型已经有黑龙江的土壤类型图可以把HSG栅格作为辅助图层和土壤类型做交集统计得到每个土壤类型的主要HSG。如果没有高精度土壤类型图直接把第3章生成的1/2/3/4栅格当作四类土壤图也能跑通但需要在usersoil表里补齐每类的物理属性包括土层深度、砂粒、黏粒、有机碳和容重。HYSOGs250m本身不提供这些属性所以别指望只复制一个HSG字母就能完成SWAT土壤数据库配置。4.2 按土地利用和HSG查CN值CN值的确定方法在USDA TR-55里有标准查找表。对黑龙江这种垦殖率高的区域最常见的组合是农田、草地和林地。下表给出一组常用参考值土地利用/覆盖HSG AHSG BHSG CHSG D普通行播作物耕地67788589牧场/草地39617480混合林地30557077不透水商业区89929495这些数值是湿润吸水条件下的CN2来自TR-55常用查算表实际项目要根据耕作方式、植被覆盖率做调整。比如免耕地的CN值比普通翻耕地低3到5个点三江平原水田的CN值又有另一套取值需要结合当地实验数据。另外CN2会随土壤含水量分三档调整SWAT会自己计算CN1和CN3所以HSG只负责提供基准档。4.3 生成SWAT可直接读入的土壤数据行如果你采用“四类HSG当土壤图”的简化方案可以写一个Python脚本按HSG类别生成usersoil表需要的一些字段soil_rows [] for hsg_code in [1, 2, 3, 4]: if hsg_code 1: hydgrp A elif hsg_code 2: hydgrp B elif hsg_code 3: hydgrp C else: hydgrp D soil_rows.append({ SNAM: fHSG_{hydgrp}, NLAYERS: 2, HYDGRP: hydgrp, SOL_Z1: 300, SOL_Z2: 1200, SOL_BD1: 1.4, SOL_AWC1: 0.1, SOL_K1: 100.0 if hydgrp A else 10.0 if hydgrp B else 3.0 if hydgrp C else 0.3, })这只是一个最小示例真正进SWAT还需要填入每个土层的黏粒、砂粒、有机碳和电导率。我的做法是先用HSG分组再把原SoilGrids250m对应位置的平均属性填进去。这样既保留了250米空间差异又满足了SWAT对土壤物理参数的要求。土层厚度和饱和导水率不能直接从HSG映射必须参照同区域的土壤调查数据否则CN值对、产流却对不上。执行SWAT的HRU分析时建议把土壤面积阈值调低。因为四类HSG在黑龙江省的空间分布极不均匀三江平原大部分是D类或CD类如果把土壤阈值设成20%小面积的B类林地黄土可能被过滤导致HRU类型缺失。常见做法是先用5%以下阈值生成HRU再在率定时合并相近水文响应单元。5. 验证黑龙江省径流模拟时容易忽略的HSG细节5.1 双HSG和洼地地下水黑龙江省的HSG分布不能只看A/D比例。三江平原、挠力河、乌裕尔河一带春季土壤解冻和冻层顶托会形成临时地下水地表60厘米内出现水位的概率极高。此时双代码里的D才是实际状态。我建议你在雨季模拟前把重分类结果中任意含D的像元挑出来叠加到DEM上看是否落在洼地、平坦耕地或沼泽区。如果吻合度高说明这个保守处理是合理的。5.2 用实测流量算纳什系数来判断HSG是否合适一个快速验证HSG是否低估或高估产流的方法把SWAT的模拟日径流与水文站实测日径流放在一起算Nash-Sutcliffe系数NSE。import numpy as np obs np.array(obs_values) # 实测日径流单位m3/s sim np.array(sim_values) # 模拟日径流单位m3/s nse 1 - np.sum((obs - sim) ** 2) / np.sum((obs - np.mean(obs)) ** 2) print(fNSE {nse:.2f})NSE接近1说明模拟精度很好接近0说明模型只比把平均值当预测好一点。如果NSE为负优先检查HSG双代码处理尤其要看你是否把A/D都按A类处理了。把含D的双代码统一改为D后CN2值会上升6到10个点洪峰响应会有明显改善。5.3 检查投影和NoData最后提交前再跑一次gdalinfo HYSOGs250m_hlj_aea.tif | grep -E Corner|NoData看输出里四个角点坐标是否覆盖黑龙江全域NoData是否和重分类后的0一致。如果边界出现空隙多半是在gdalwarp时没有保留全栅格范围用-crop_to_cutline配合省级矢量边界再切一次即可。本文还有配套的精品资源点击获取 黑龙江省HSG栅格处理与SWAT模型径流模拟实战
简介这份黑龙江省土壤水文分组高精度栅格数据集源自HYSOGs250m经省级裁剪配准面向水文建模、水土保持与SWAT模型应用人群。数据采用WGS84坐标系提供约250米分辨率的A/B/C/D四类土壤径流潜力栅格同时标注浅地下水位引发的双重HSG可直观用于降雨径流潜力分区与CN值估算。包体共7个文件以tif栅格为主配套tfw坐标参数、dbf属性表、xml元数据、pdf说明文档与cpg编码满足主流GIS平台便捷加载压缩包仅4.28MB轻量易用。已有85人学习浏览。借助该数据集可直接获得黑龙江省统一的土壤水文分组成果省去全球数据下载、拼接与裁剪流程。数据附带的属性表与元数据便于批量读取和制图适合SWAT建模前期参数准备、省域径流潜力对比及教学科研中的快速演示。1. 黑龙江省HSG栅格从默认土壤到250米径流潜力分级做SWAT模型的人大多在土壤数据上栽过跟头。默认的1千米土壤图在黑龙江这种地形破碎、三江平原湿地连片的地方跑出来的CN值不是偏大就是偏小率定过程很痛苦。HYSOGs250m把USDA曲线数方法依赖的土壤水文分组HSG放到了250米分辨率上每个像元给出A、B、C、D四档径流潜力湿润区还标了双重HSG。这个栅格基于SoilGrids250m派生已经按黑龙江省界裁剪好坐标系是WGS84适合直接进入SWAT的径流模拟流程。对于做水文建模和GIS数据预处理的人来说关键是先读懂分类逻辑再把它转成SWAT认识的格式。2. HYSOGs250m的底账A/B/C/D分类、双重HSG与栅格文件结构2.1 四档HSG如何从SoilGrids250m派生HYSOGs250m没有直接钻探土壤而是用FAO的SoilGrids250m系统给出的土壤质地等级和到基岩的深度在1/480十进制度网格上推演出一个“排水性能”分级。SCS-CN方法里土壤被分成A、B、C、D四档代表降雨入渗和径流潜力的顺序。A类对应砂质或砾质土壤导水性强降雨下去很快B类是砂壤土下渗中等C类偏黏壤土下渗受限D类则是黏土或浅层不透水基本靠地表产流。这里有一组常用的最小下渗率参考值虽然不是HYSOGs250m直接输出的字段但能帮你理解分级边界HSG典型土壤条件最小下渗率mm/h径流潜力A砂土、砾质土深厚且排水好7.6以上低B砂壤土土层较深3.8–7.6中等偏低C黏壤土、细砂壤土下渗受限1.3–3.8中等偏高D黏土、浅层土壤或有地下水位小于1.3高这个表是SCS经典公式的参数。你在SWAT里填CN值时如果CN表里给出A/B/C/D四列用的就是这套逻辑。HYSOGs250m把SoilGrids250m的连续属性离散成了这四档所以它本质上是一个分类图而不是属性图。不要把栅格值当作可平滑处理的连续数值做重采样必须用最近邻。2.2 双重HSG是怎么回事USDA的HSG分类里还有一个特殊情况如果地表以下60厘米范围内存在地下水位土壤不管表层是砂质还是黏质在雨季都表现出高径流潜力。HYSOGs250m会为这类区域赋予双重HSG比如A/D、B/D或C/D。前者是排水良好的干燥状态后者是饱和状态。黑龙江省的三江平原、挠力河流域和松嫩低洼地区经常出现这种双代码。SWAT标准的usersoil表里只有一列HYDGRP装不下双代码。我一般的做法是在雨季模拟时取保守值比如把A/D当作D使用如果在做季节率定也可以利用双代码做月尺度切换。要注意的是只看字母A就把它当成低产流在三江平原会明显低估洪峰。2.3 解压后每个文件都有用解压后拿到一堆文件初看好像冗余其实每一项都有固定角色尤其是.vat.dbf经常被忽略很多人的做法是直接读tif然后靠猜这很容易在双代码上翻车。下面先列最常见的几个文件文件作用HYSOGs250m_黑龙江省.tif主栅格像元值对应HSG代码.tfw世界文件记录影像左上角坐标和像元尺寸.aux.xmlESRI辅助元数据记录统计信息和色彩映射.vat.dbf栅格值属性表把栅格Value翻译成HSG字母.vat.cpgVAT表的字符编码文件土壤水文分组.pdf数据说明文档.vat.dbf是理解栅格值的关键。可以用GDAL自带的ogrinfo直接看ogrinfo -al HYSOGs250m_黑龙江省.tif.vat.dbf输出里通常有Value和HSG两个字段Value是tif里存的整数HSG是对应的字母或双字母代码。先跑这条命令确认代码含义再继续处理不要凭经验猜。否则后面统计面积或重分类时把A/D当成某个数字就会出偏差。如果PDF说明文档和VAT表对不上以VAT表为准因为VAT是随栅格一起发布的PDF可能描述的是整个中国范围而黑龙江省版已经做过裁剪和重编码。另外比起全国版栅格这个版本经过了省级裁剪文件体积更小但边缘像元仍可能覆盖邻省或省界外的无意义区域。后续处理建议用黑龙江行政边界再做一次裁剪避免统计面积时混入界外像元。3. GDALPython处理黑龙江省HSG栅格投影、统计与重映射现在进入实操。处理这类区域裁剪后的HSG栅格我一般分成四步看元数据、读VAT、投影后做面积统计、重新编码输出。3.1 打开栅格前先做一次体检拿到HYSOGs250m_黑龙江省.tif先跑gdalinfo HYSOGs250m_黑龙江省.tif重点看四行Size、Origin、Pixel Size和Coordinate System。这个数据集是WGS84地理坐标系像素分辨率是0.00208333度约250米。NoData也要记下来通常可能是-9999或0后面统计时要过滤。再看旁边的.tfw世界文件。一个标准的TFW文件有6行数字前两行是像元在X和Y方向上的尺寸中间两行是旋转项正常为0第五、六行是左上角像元中心坐标。用记事本打开就能确认行列范围是否覆盖黑龙江边界尤其是跨到俄罗斯一侧的区域不需要额外做范围裁切。3.2 用VAT把Value翻译成HSG代码我建议一上来就把VAT读进内存生成一个Value到HSG字母的映射。这样后续所有统计都能直接打印字母。from osgeo import ogr vat ogr.Open(HYSOGs250m_黑龙江省.tif.vat.dbf) layer vat.GetLayer(0) mapping {} while True: feat layer.GetNextFeature() if feat is None: break value feat.GetField(Value) hsg feat.GetField(HSG) mapping[int(value)] hsg print(f{value} - {hsg})这段代码遍历VAT表的每一行取出Value和HSG两个字段放进dict。后续对栅格数组做np.unique后可以用这个dict把整数值翻译成HSG字母。注意字段名不一定叫HSG有的版本叫Label或Class可以先打印layer.GetLayerDefn()确认字段列表。3.3 投影后做面积统计WGS84地理坐标系下每个像元在不同纬度对应的真实面积不一样。黑龙江省跨北纬43度到53度如果直接用250乘以250统计面积误差在5%以上。面积统计前我先转到Albers等积投影。gdalwarp -t_srs projaea lat_145 lat_260 lat_025 lon_0105 datumWGS84 \ -tr 250 250 -r near -overwrite \ HYSOGs250m_黑龙江省.tif HYSOGs250m_hlj_aea.tif该命令中-t_srs指定Albers等积圆锥投影两条标准纬线选45和60度覆盖黑龙江主体-tr 250 250把输出像元严格设为250米-r near用最近邻插值保证分类值不产生新层级。投影完成后栅格每个像元面积都是62500平方米统计才有意义。Python统计各HSG面积from osgeo import gdal import numpy as np ds gdal.Open(HYSOGs250m_hlj_aea.tif) arr ds.GetRasterBand(1).ReadAsArray() arr arr.astype(np.int16) # 过滤NoData按原始Value统计 valid_mask (arr 0) values, counts np.unique(arr[valid_mask], return_countsTrue) pixel_area_km2 250 * 250 / 1e6 for v, c in zip(values, counts): hsg mapping.get(int(v), unknown) print(f{hsg}: {c * pixel_area_km2:.2f} km2)astype(np.int16)是为了防止某些配置下栅格以无符号整型读入负号和NoData变成大数valid_mask过滤掉小于等于0的像元。如果NoData是-9999arr 0同样可以把它滤掉。统计结果出来后你会看到黑龙江各HSG类的大致比例通常C和D占大头这与黑土和草甸土分布基本吻合。3.4 重分类并导出SWAT可用的GRID统计完之后通常要重编码。比如把双代码A/D统一赋成4生成一个只含1、2、3、4的四类栅格。这样后面给SWAT的土壤图时字段更简洁。reclass np.zeros_like(arr) for v, h in mapping.items(): if D in h: reclass[arr v] 4 elif h C: reclass[arr v] 3 elif h B: reclass[arr v] 2 elif h A: reclass[arr v] 1这段代码把任何包含D的双代码都归到第4类其余A/B/C对应1/2/3。如果你做湿润季节模拟这是最稳妥的保守处理。保存重分类结果时要继承原始投影和仿射参数driver gdal.GetDriverByName(GTiff) out_ds driver.Create(reclass.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Int16) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(reclass) out_band.SetNoDataValue(0) out_ds.FlushCache()最后用gdal_translate转成文本GridSWAT可以直接读取gdal_translate -of AAIGrid -ot Int16 -a_nodata 0 reclass.tif hsg_reclass.asc这样HSG空间分布图就准备好进入SWAT流程了。注意AAIGrid文件头里的ncols和nrows必须和DEM一致如果DEM分辨率不一样需要先重采样到同一格网。4. 在SWAT模型里把HSG变成CN值从土壤属性到HRU4.1 SWAT中的HYDGRP与usersoil表SWAT不是直接用栅格里的HSG参与产流计算的。它把HSG存在土壤数据库的HYDGRP字段里土地利用和土壤类型组合成HRU后再通过CN2值进入径流曲线数计算。所以在SWAT里看到的是“每一类土壤对应一个A/B/C/D”不是一个逐像元字段。如果你的模型已经有黑龙江的土壤类型图可以把HSG栅格作为辅助图层和土壤类型做交集统计得到每个土壤类型的主要HSG。如果没有高精度土壤类型图直接把第3章生成的1/2/3/4栅格当作四类土壤图也能跑通但需要在usersoil表里补齐每类的物理属性包括土层深度、砂粒、黏粒、有机碳和容重。HYSOGs250m本身不提供这些属性所以别指望只复制一个HSG字母就能完成SWAT土壤数据库配置。4.2 按土地利用和HSG查CN值CN值的确定方法在USDA TR-55里有标准查找表。对黑龙江这种垦殖率高的区域最常见的组合是农田、草地和林地。下表给出一组常用参考值土地利用/覆盖HSG AHSG BHSG CHSG D普通行播作物耕地67788589牧场/草地39617480混合林地30557077不透水商业区89929495这些数值是湿润吸水条件下的CN2来自TR-55常用查算表实际项目要根据耕作方式、植被覆盖率做调整。比如免耕地的CN值比普通翻耕地低3到5个点三江平原水田的CN值又有另一套取值需要结合当地实验数据。另外CN2会随土壤含水量分三档调整SWAT会自己计算CN1和CN3所以HSG只负责提供基准档。4.3 生成SWAT可直接读入的土壤数据行如果你采用“四类HSG当土壤图”的简化方案可以写一个Python脚本按HSG类别生成usersoil表需要的一些字段soil_rows [] for hsg_code in [1, 2, 3, 4]: if hsg_code 1: hydgrp A elif hsg_code 2: hydgrp B elif hsg_code 3: hydgrp C else: hydgrp D soil_rows.append({ SNAM: fHSG_{hydgrp}, NLAYERS: 2, HYDGRP: hydgrp, SOL_Z1: 300, SOL_Z2: 1200, SOL_BD1: 1.4, SOL_AWC1: 0.1, SOL_K1: 100.0 if hydgrp A else 10.0 if hydgrp B else 3.0 if hydgrp C else 0.3, })这只是一个最小示例真正进SWAT还需要填入每个土层的黏粒、砂粒、有机碳和电导率。我的做法是先用HSG分组再把原SoilGrids250m对应位置的平均属性填进去。这样既保留了250米空间差异又满足了SWAT对土壤物理参数的要求。土层厚度和饱和导水率不能直接从HSG映射必须参照同区域的土壤调查数据否则CN值对、产流却对不上。执行SWAT的HRU分析时建议把土壤面积阈值调低。因为四类HSG在黑龙江省的空间分布极不均匀三江平原大部分是D类或CD类如果把土壤阈值设成20%小面积的B类林地黄土可能被过滤导致HRU类型缺失。常见做法是先用5%以下阈值生成HRU再在率定时合并相近水文响应单元。5. 验证黑龙江省径流模拟时容易忽略的HSG细节5.1 双HSG和洼地地下水黑龙江省的HSG分布不能只看A/D比例。三江平原、挠力河、乌裕尔河一带春季土壤解冻和冻层顶托会形成临时地下水地表60厘米内出现水位的概率极高。此时双代码里的D才是实际状态。我建议你在雨季模拟前把重分类结果中任意含D的像元挑出来叠加到DEM上看是否落在洼地、平坦耕地或沼泽区。如果吻合度高说明这个保守处理是合理的。5.2 用实测流量算纳什系数来判断HSG是否合适一个快速验证HSG是否低估或高估产流的方法把SWAT的模拟日径流与水文站实测日径流放在一起算Nash-Sutcliffe系数NSE。import numpy as np obs np.array(obs_values) # 实测日径流单位m3/s sim np.array(sim_values) # 模拟日径流单位m3/s nse 1 - np.sum((obs - sim) ** 2) / np.sum((obs - np.mean(obs)) ** 2) print(fNSE {nse:.2f})NSE接近1说明模拟精度很好接近0说明模型只比把平均值当预测好一点。如果NSE为负优先检查HSG双代码处理尤其要看你是否把A/D都按A类处理了。把含D的双代码统一改为D后CN2值会上升6到10个点洪峰响应会有明显改善。5.3 检查投影和NoData最后提交前再跑一次gdalinfo HYSOGs250m_hlj_aea.tif | grep -E Corner|NoData看输出里四个角点坐标是否覆盖黑龙江全域NoData是否和重分类后的0一致。如果边界出现空隙多半是在gdalwarp时没有保留全栅格范围用-crop_to_cutline配合省级矢量边界再切一次即可。本文还有配套的精品资源点击获取