简介四川省土壤水文分组高精度栅格数据集专为SWAT水文建模、降雨径流估算及土壤入渗特性分析而设计提供基于USDA曲线数CN方法的HSG分类结果。数据源采用HYSOGs250m方案依据FAO soilGrids250m提供的土壤质地等级与基岩深度生成空间分辨率约250米栅格将土壤划分为A、B、C、D四个等级分别对应低、中低、中高及高径流潜力并对60厘米深度内存在地下水位的潮湿土壤附双重HSG标识便于模型捕捉特殊水文行为、减少预处理工作量。该成果已按四川省界裁剪坐标系为WGS84可直接用于SWAT模型参数输入也可在ArcGIS、QGIS等平台中完成流域汇水分析与水土保持评价对复杂地形下的产流潜力判断有明显帮助。压缩包体积约4.11MB内容以栅格数据文件为主整体轻量易用。目前已有211人学习/下载适合水文水利、地理信息及环境科学等方向的中高级研究人员作为基础数据使用。1. 从 CN 数到栅格为什么四川需要一套独立的水文土壤分组做 SWAT 或 HEC-HMS 的人都知道曲线数Curve Number里最容易被“凑合用”的输入就是水文土壤组HSG。大多数项目直接套用美国本土的 STATSGO 派生数据或者用土壤质地粗略一估结果率定的时候怎么调参数都不对。这套 HYSOGs250m 数据的关键在于它把 FAO SoilGrids250m 的土壤质地和基岩深度按 USDA 的标准重分类成了 A、B、C、D 四类径流潜力等级并且在 250 米分辨率上覆盖四川全省。四川的情况比较特殊。成都平原的冲积土、川西高原的高山草甸土、盆周山地的紫色土物理性质差异极大更要命的是若尔盖一带存在大量地表 60 厘米内就有地下水位的潮湿土壤这类土壤无论质地如何都应该归为高径流潜力的双重 HSG 等级。这套数据已经按省级行政区裁剪过坐标系是 WGS84 地理坐标系拿过来可以直接进 ArcGIS 或 QGIS 做 SWAT 建模不用自己拼图、重投影或者用行政边界去裁剪。下面先把分类原理拆开再讲实际操作中的重分类、镶嵌和坑。2. HSG 分类原理与 HYSOGs250m 的派生链路2.1 USDA 四类 HSG 的判定标准USDA-NRCS 的 HSG 分类不是直接测出来的而是根据土壤的入渗率、导水率和地下水位深度综合判定的。标准定义是A 类为低径流潜力砂质或砾质土壤饱和导水率大于 0.40 in/hrB 类为中等偏低的径流潜力砂壤土或壤土饱和导水率 0.15 到 0.40 in/hrC 类为中等偏高的径流潜力砂质黏壤土或黏壤土导水率 0.05 到 0.15 in/hrD 类为高径流潜力黏土或膨胀性土壤导水率小于 0.05 in/hr。实际应用中还有一个容易忽略的规则即便土壤质地属于 A 类或 B 类如果地下水位在地表 60 厘米以内也要强制归为 D 类这就是所谓的双重 HSG 或潮湿土壤等级。这套 HYSOGs250m 数据正是考虑了这一点所以才在分类逻辑里引入了到基岩深度和潜水位因素。2.2 从 soilGrids250m 到 HSG 重分类的具体规则HYSOGs250m 的原始派生逻辑来自 Ross 等人发表的全球水文土壤组数据集。它基于 soilGrids250m 提供的 0-30 cm 和 30-60 cm 两层的 USDA 土壤质地类别以及到基岩的深度bedrock depth执行如下规则任意一层质地为砂土Sand或壤砂土Loamy Sand且基岩深度大于 100 cm则归为 A 类。任意一层质地为砂壤土Sandy Loam或壤土Loam且不满足 A 类条件则归为 B 类。任意一层质地为砂质黏壤土Sandy Clay Loam、粉砂壤土Silt Loam或黏壤土Clay Loam则归为 C 类。任意一层质地为砂质黏土Sandy Clay、粉砂黏土Silty Clay或黏土Clay则归为 D 类。任一层出现有机土壤或基岩深度小于 60 cm直接归为 D 类。地表 60 cm 内检测到潜水位不论质地如何均赋予双重 HSG代码中单独标识为 A/D 或 B/D 等。这套数据在四川的应用中有一个值得注意的点川西高原的土壤基岩深度普遍较浅很多区域直接命中“基岩深度小于 60 cm 归 D 类”这条规则导致 D 类面积占比看起来偏高。实际上这是符合 USDA 定义的因为浅薄土层确实几乎没有蓄水能力降雨会快速形成地表径流。2.3 WGS84 坐标系下 250 米栅格的实际含义数据的地理分辨率为 1/480 十进制度也就是 0.002083 度对应赤道附近大约 250 米。但需要注意在高纬度地区东西方向的实际距离会变短四川的纬度范围大致在北纬 26° 到 34° 之间在这个区间内 0.002083 度经度对应约 180 到 220 米纬度方向始终约 232 米。这意味着 SWAT 建模时如果 HRU 划分阈值设置得过细会出现同一个 HRU 边界横跨多个 HSG 栅格单元的情况。常见的做法是在 SWAT 的土壤库中直接使用 HSG 代码因为这个栅格数据本身不提供完整的土壤物理属性容重、有机碳、饱和导水率等那些属性仍然要来自 HWSD 或者本地的土壤普查数据。提示这套栅格数据的值域是 1 到 4分别对应 A、B、C、D。双重 HSG 在部分版本中有单独编码拿到数据后第一步就是检查属性表和唯一值不要默认只有四个值。3. 四川省 HSG 栅格的预处理与完整复现流程3.1 数据检查与属性确认拿到四川省的 HSG 栅格后第一步不要急着做分析先用 GDAL 检查基本属性gdalinfo sichuan_hsg.tif重点看Type是否为 Byte 或 Int16NoData Value是否为 -9999 或 0以及Coordinate Reference System是否为 WGS84。如果源数据是从全球镶嵌图中按省界裁剪出来的边缘区域可能出现 NoData 环绕这属于正常现象不需要惊慌。再检查唯一值分布gdallinfo -stats sichuan_hsg.tif或者在 Python 中用 rasterio 直接读取并计数import rasterio import numpy as np with rasterio.open(sichuan_hsg.tif) as src: data src.read(1) vals, counts np.unique(data, return_countsTrue) for v, c in zip(vals, counts): print(fHSG {v}: {c} pixels, {c * 250 * 250 / 1e6:.2f} km²)这段代码的作用是统计每个 HSG 级别的像素数量和对应的面积便于快速判断分类结果是否符合区域认知。如果 D 类面积超过 50%就要回去检查是否把 NoData 值当成了有效值。np.unique会把 NoData 也统计进来所以需要在读取时就设置maskedTrue。3.2 重投影与重采样策略HYSOGs250m 的原始坐标系是 WGS84 地理坐标系但 SWAT 建模通常使用投影坐标系。四川地处中纬度推荐使用 Albers 等积圆锥投影中央经线 105°E标准纬线 25°N 和 47°N这是中国区域制图的标准方案。gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 \ -tr 250 250 -r near \ sichuan_hsg.tif sichuan_hsg_aea_250m.tif这里的-r near是分类栅格重采样的关键参数。HSG 是类别数据不是连续变量使用双线性或三次卷积会插值出诸如 2.5 之类的无效值。虽然 gdalwarp 会在输出时对非整数做取整但在类别边界上会产生非语义性的混合分类宁可损失少许空间精度也要保证类别纯净。import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:3415 with rasterio.open(sichuan_hsg.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, resolution250) profile src.profile.copy() profile.update(crsdst_crs, transformtransform, widthwidth, heightheight) with rasterio.open(sichuan_hsg_aea_250m.tif, w, **profile) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest)calculate_default_transform会根据目标分辨率和范围自动计算输出栅格的行列数。注意如果源栅格是 1/480 度的地理坐标resolution250对 Albers 投影来说才是真正的 250 米。代码中所有重采样方式都是Resampling.nearest与命令行版本的-r near保持一致的语义。3.3 镶嵌多省数据或拼接相邻图幅四川省面积约 48.6 万平方公里在全中国范围数据集中通常会被拆成多块。拿到手后若发现四川被分成多个文件需要先检查重叠区域的分类是否一致再执行镶嵌。gdalbuildvrt sichuan_hsg.vrt sichuan_part1.tif sichuan_part2.tif sichuan_part3.tif gdal_translate sichuan_hsg.vrt sichuan_hsg_merged.tifgdalbuildvrt不复制数据而是生成一个虚拟镶嵌文件最后用gdal_translate落盘。重叠区域如果分类值一致VRT 默认取第一个文件的像素如果分类值不一致说明裁剪边界附近存在分类冲突需要回到原始全球数据核查。注意不要使用gdal_merge.py做拼接。该工具遇到带有 NoData 的区域时可能会把 NoData 源覆盖有效像元产生不规则的洞。VRT 方案在语义上更安全。4. SWAT 建模中 HSG 栅格的实际集成与参数映射4.1 SWAT 土壤库中 HSG 字段的填写方式SWAT 模型的usersoil数据库中有HYDGRP字段取值范围就是 A、B、C、D。这个字段直接参与了 SCS-CN 方程中 CN2 的查表计算。SWAT 自带的cn2查表是基于 HSG 和土地利用类型的组合确定 CN 值的所以 HSG 的错误分类会直接传导至径流模拟结果。从栅格到 SWAT 土壤库的转换逻辑如下在 ArcGIS 中把 HSG 栅格转为矢量多边形。与 SWAT 的 HRU 边界做空间连接提取每个 HRU 中的优势 HSG 类别。在土壤属性表里为每种 HSG 分配对应的HYDGRP字符值。这里有一个关键问题SWAT 支持的最细土壤单位是单个 HRU但每个 HRU 的面积通常远大于 250 米栅格单元一个 HRU 内部完全可能同时包含 B 类和 C 类。常见做法是用面积占比最大的类别或者按照 CN 值加权平均后反查 HSG。后者更精确但需要先算出各 CN 值再聚合。4.2 用 Python 批量生成 SWAT 土壤库 HSG 字段import pandas as pd soil_db pd.read_csv(usersoil.csv, encodinggbk, dtype{MUID: str}) def map_hsg(value): mapping {1: A, 2: B, 3: C, 4: D} return mapping.get(value, D) soil_db[HYDGRP] soil_db[S5ID].apply( lambda x: map_hsg(int(x.split(_)[-1])))这个片段演示了如何从土壤编号中解析 HSG 值并将其映射到 SWAT 的HYDGRP字段。实际项目里S5ID可能不包含 HSG 信息这时需要通过空间连接的结果构建映射表。由于 SWAT 的土壤数据库对字符串大小写敏感统一用大写字母否则后续 HRU 划分时会报错找不到对应的 HSG。4.3 在 QGIS 中快速预览 HSG 分布QGIS 加载 HSG 栅格后用栅格计算器或直接设置Pseudo-color渲染方案即可。建议的颜色方案为A 类用浅蓝低径流B 类用浅绿中低径流C 类用橙色中高径流D 类用红色高径流。这样可以在 30 秒内识别出成都平原和若尔盖湿地的 HSG 差异对检查数据合理性非常直观。# 输出成都平原的 HSG 面积统计 gdal_calc.py -A sichuan_hsg.tif --outfilechengdu_plain_hsg.tif \ --calcA * (A 0) --quiet这不是一个精确的按区域裁剪而是通过--calc表达式过滤 NoData 的值。精确的成都平原范围裁剪需要准备矢量边界文件然后用gdalwarp -cutline完成。5. 潮湿土壤双重 HSG 的处理细节与常见误用排查5.1 双重 HSG 在数据中如何识别官方 HYSOGs250m 数据中双重 HSG 通常被编码为组合值例如 5 表示 A/D6 表示 B/D7 表示 C/D。但四川省裁剪版可能只保留标准值 1-4因为裁剪过程重新映射了属性表。务必做唯一值检查gdalinfo -hist sichuan_hsg.tif | tail -50如果Histogram中出现 5、6、7 或超出 1-4 范围的数值说明双重 HSG 信息被保留了下来。SWAT 本身不接受组合 HSG 输入必须把 A/D 等组合值映射为单一 HSG。两条经验规则对于常年积水的湿地、沼泽、水稻田映射为 D 类因为 SWAT 的 CN 查表在 D 类下更接近实际积水产流行为。对于表层干燥但 60 cm 内有地下水位的区域映射为 C 类避免 D 类导致径流过度高估。5.2 四川典型区域 HSG 特征对照以四川省几种典型地貌为例区域主要 HSG成因说明成都平原B 类为主局部 A 类岷江冲积物质地以壤土和砂壤土为主排水良好川中丘陵C 类为主紫色页岩风化形成的紫色土黏粒含量偏高若尔盖湿地D 类双重 HSG地表 60 cm 内潜水位高无论质地一律高径流川西高山峡谷D 类浅基岩区域基岩深度小于 60 cm土层极薄提示如果成渝城市群的大部分区域模拟结果都出现 C 类或 D 类不要急着怀疑数据有问题先核对原始 soilGrids250m 的质地分类。5.3 高频错误与快速排查方法最常见的一个错误是直接用重分类工具把无双值关系的栅格值转换为文本 HSG 代码时代码写反了。1A 和 4A 是完全相反的结果前者代表低径流后者代表高径流。检查方式是在成都平原区域做一次局部验证。import rasterio with rasterio.open(sichuan_hsg_aea_250m.tif) as src: # 成都平原约在 103.8E, 30.7N row, col src.index(103.8, 30.7) print(src.read(1)[row, col])如果非常确定该坐标位于成都平原壤土为主输出值却是 4D 类那就要检查原始数据在裁剪阶段是否发生了值域的翻转或偏移。src.index根据地理坐标计算行列号注意 WGS84 下经度在前纬度在后。另一个高频问题出在投影变换后栅格的像元值全部变成 NoData。原因通常是gdalwarp输出范围计算错误源栅格的部分像元落在了目标坐标系的有效区域之外或者-te参数设置不当。解决方法是删掉-te参数让 gdalwarp 自动计算边界或提供略大于四川省边界的范围。6. 基于 HSG 栅格验证 CN 数与径流结果的实用技巧6.1 用历史降雨事件反向验证 HSG 分类如果项目区有实测的降雨-径流序列可以用一个非常轻量级的验证方法对一次孤立降雨事件使用 SCS-CN 方程反算 CN 值。$$Q \frac{(P - 0.2S)^2}{P 0.8S}, \quad S \frac{25400}{CN} - 254$$将实测降雨量 P 和直接径流深 Q 代入公式反解 CN 值。然后将反算的 CN 值对照 SCS 标准 CN 查算表判断当前 HSG 分类是否合理。如果实测反算 CN 落在 A 类区间但你的模型用的是 C 类说明 HSG 分类过重土壤实际产流能力比预期低。这个方法不适合溶蚀性基岩地区喀斯特地貌因为地下漏失会使地表径流偏小反算出的 CN 值系统性偏低。四川兴文、叙永一带的喀斯特区域如果使用此方法需要先判断是否存在明显的落水洞和地下河。6.2 使用 Python 实现单事件 CN 反算def cn_from_runoff(P, Q): if Q 0 or P 0: return None S 0.2 * P 0.8 * Q * P / (P - Q 1e-6) - 0.8 * Q CN 25400 / (S 254) return max(0.01, min(100, CN)) # 示例某次降雨 42mm实测径流深 8mm print(cn_from_runoff(42, 8))这段代码先由 P 和 Q 反解 S再代入 CN 换算公式。注意公式中隐含了初损 Ia 0.2S 的假设这是 SCS-CN 模型的标准配置。代码中的1e-6是为了防止除零。如果反算结果大于 100 或小于 0说明实测数据可能不符合 SCS-CN 的适用条件应该剔除该场次。6.3 多栅格交叉对比提升用数据可信度HYSOGs250m 毕竟只有 250 米分辨率在四川盆周山区一个栅格单元可能横跨河谷和坡地HSG 分类存在混合语义。如果项目精度要求高于 250 米建议与以下数据做交叉验证HWSD v2.0 提供的土壤质地和有效土层深度1 公里分辨率中国 1:100 万土壤图数字化属性矢量若项目区有当地土壤普查的典型剖面数据交叉验证的逻辑很简单在同一个空间位置上如果 HYSOGs250m 判定为 A 类砂土、深基岩但 HWSD 显示该位置是黏土且浅层有岩石出露那就要做一次人工核查。gdal_calc.py可以快速生成两个数据源的差值栅格配合 QGIS 的底图进行目视判读。整体来看这套四川省 HSG 数据最值得依赖的部分是对潮湿土壤和浅基岩区域的处理逻辑这是许多国产土壤数据缺失的维度。在 SWAT 建模中把它用对能把率定初期的不确定性挪走一块后续调参也就更有方向感。本文还有配套的精品资源点击获取