四川高分辨率水文土壤组HSG栅格数据解析与SWAT建模实践

简介：四川省土壤水文分组高精度栅格数据集专为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/hr；B 类为中等偏低的径流潜力，砂壤土或壤土，饱和导水率 0.15 到 0.40 in/hr；C 类为中等偏高的径流潜力，砂质黏壤土或黏壤土，导水率 0.05 到 0.15 in/hr；D 类为高径流潜力，黏土或膨胀性土壤，导水率小于 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 或 Int16，NoData 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_counts=True) for v, c in zip(vals, counts): print(f'HSG {v}: {c} pixels, {c * 250 * 250 / 1e6:.2f} km²')

这段代码的作用是统计每个 HSG 级别的像素数量和对应的面积，便于快速判断分类结果是否符合区域认知。如果 D 类面积超过 50%，就要回去检查是否把 NoData 值当成了有效值。np.unique会把 NoData 也统计进来，所以需要在读取时就设置masked=True。

3.2 重投影与重采样策略

HYSOGs250m 的原始坐标系是 WGS84 地理坐标系，但 SWAT 建模通常使用投影坐标系。四川地处中纬度，推荐使用 Albers 等积圆锥投影（中央经线 105°E，标准纬线 25°N 和 47°N），这是中国区域制图的标准方案。

gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84" \ -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, resolution=250) profile = src.profile.copy() profile.update(crs=dst_crs, transform=transform, width=width, height=height) with rasterio.open('sichuan_hsg_aea_250m.tif', 'w', **profile) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=dst_crs, resampling=Resampling.nearest)

calculate_default_transform会根据目标分辨率和范围自动计算输出栅格的行列数。注意如果源栅格是 1/480 度的地理坐标，resolution=250对 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.tif

gdalbuildvrt不复制数据，而是生成一个虚拟镶嵌文件，最后用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 土壤库的转换逻辑如下：

1. 在 ArcGIS 中把 HSG 栅格转为矢量多边形。
2. 与 SWAT 的 HRU 边界做空间连接，提取每个 HRU 中的优势 HSG 类别。
3. 在土壤属性表里为每种 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', encoding='gbk', 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 --outfile=chengdu_plain_hsg.tif \ --calc="A * (A > 0)" --quiet

这不是一个精确的按区域裁剪，而是通过--calc表达式过滤 NoData 的值。精确的成都平原范围裁剪需要准备矢量边界文件，然后用gdalwarp -cutline完成。

5. 潮湿土壤双重 HSG 的处理细节与常见误用排查

5.1 双重 HSG 在数据中如何识别

官方 HYSOGs250m 数据中，双重 HSG 通常被编码为组合值，例如 5 表示 A/D，6 表示 B/D，7 表示 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 特征对照

以四川省几种典型地貌为例：

提示：如果成渝城市群的大部分区域模拟结果都出现 C 类或 D 类，不要急着怀疑数据有问题，先核对原始 soilGrids250m 的质地分类。

5.3 高频错误与快速排查方法

最常见的一个错误是：直接用重分类工具把无双值关系的栅格值转换为文本 HSG 代码时，代码写反了。1=A 和 4=A 是完全相反的结果，前者代表低径流，后者代表高径流。检查方式是在成都平原区域做一次局部验证。

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])

如果非常确定该坐标位于成都平原（壤土为主），输出值却是 4（D 类），那就要检查原始数据在裁剪阶段是否发生了值域的翻转或偏移。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 建模中把它用对，能把率定初期的不确定性挪走一块，后续调参也就更有方向感。

本文还有配套的精品资源，点击获取