上海土壤水文分组高精度栅格制作与SCS-CN径流模拟应用

简介：上海市土壤水文分组（HSG）高精度栅格数据，是基于USDA曲线数（CN）方法估算降雨径流的关键基础数据，适用于SWAT水文模型建模、城市排水规划及流域水文分析。数据源自HYSOGs250m与soilGrids250m系统，按省级行政区裁剪为上海市范围，WGS84地理坐标系，空间分辨率约250m，提供A、B、C、D四类径流潜力等级，并包含湿润土壤双重HSG分类信息。压缩包大小3.46MB，文件总数0，类型明细暂无数据，可用GIS软件直接加载分析。目前已有212人学习，适合水文、环境、GIS领域的研究者与工程师，可直接用于计算径流曲线数、划分水文响应单元及开展区域水资源评估。

1. 上海市土壤水文分组高精度栅格：给产流模型一张能直接查表的底图

做上海市域内涝模拟或海绵城市绩效评估的人，迟早会被同一个问题卡住：SCS-CN 查表需要土壤水文分组（HSG，A/B/C/D 四档），而手头常见的是二普扫描图、HWSD 的 1 km 网格或 SoilGrids 的 250 m 插值。上海市土壤水文分组高精度栅格，就是把上海地表土壤按入渗能力归入四档，做成覆盖全市的 10 m 或 30 m 整幅栅格，让 SWMM、InfoWorks ICM、HEC-HMS 在每个像元上直接完成 HSG 到 CN 的映射。难在上海几乎每一项地表条件都在抬升难度：高差几米、地下水位贴着地表、硬化面把土体切得支离破碎，而 1:50 万土壤图投影到 10 m 栅格，本身就是一种"假精细"。下面按判定、制作、验证、对接四步把链路拆开。

2. A 到 D 四档怎么定：USDA 水文土壤分组的判定参数与上海土壤归属逻辑

2.1 水文土壤分组不是土壤分类，而是入渗能力分级

中国土壤发生分类里的水稻土、潮土、滨海盐土回答的是"土壤怎么来的"，USDA 水文土壤分组回答的是另一件事：长时间湿润、无冻土、无地表封盖条件下，土壤在强降雨里能吞下多少水。NRCS 在 National Engineering Handbook Part 630 中按入渗能力把土壤分成 A、B、C、D 四组：A 组通常为砂质或砾质，入渗极快；B 组中等；C 组偏慢；D 组极慢。D 组有一条补充定义对上海特别致命：即使质地不粘，只要下挖 60 cm 内遇到持续地下水位，或 50 cm 内遇到致密不透水层，也判为 D。上海大部分地区地下水位在 0.5～1.5 m，雨季贴着地表，这一条直接锁定了大片 D 组。

2.2 先看质地，再看 Ks，最后用排水条件兜底

判定流程我习惯按三步走。第一步查土种的质地剖面，按 USDA 质地三角图映射出候选组；第二步用饱和导水率 Ks 复核，实验室定水头法、现场双环入渗仪都可以；第三步用地下水埋深和不透水层深度做一票否决。三者对应关系如下：

表里的 Ks 取的是 NRCS 常用参考区间。实际判定时看的是整层剖面中最差的那个层次，不是表层——上海水稻土的耕作层透水性尚可，真正决定分组的是犁底层。

2.3 上海主要土种的归属判断

按土类分布看：南汇、奉贤沿海的滨海盐土，粉砂质粘壤土加常年高水位，归 D；青浦、松江、金山的青紫泥、青黄土等水稻土，犁底层紧实、下有潜育层，归 D；嘉定、宝山及黄浦江以西的灰潮土，粉砂壤土到粘壤土、水位 1～2 m，多为 C；崇明东滩、长兴岛新近沉积的砂质潮土，局部到 B，砂性强的可到 A。老城区内的回填土和建筑渣土扰动层，属性无法从二普图获得，单独编码为 9 类"扰动土"，不参与常规查表。具体归属如下：

2.4 从粒径或质地名称自动落组：字典映射就够

如果数据源给的是粘粒、砂粒、粉粒百分数，先按 USDA 三角图定质地类，再按字典映射落到四档：

TEXTURE_HSG = { "sand": "A", "loamy sand": "A", "sandy loam": "B", "loam": "B", "silt loam": "B", "silt": "B", "sandy clay loam": "C", "clay loam": "C", "silty clay loam": "C", "sandy clay": "D", "silty clay": "D", "clay": "D", } def texture_to_hsg(texture: str) -> str: return TEXTURE_HSG.get(texture.strip().lower(), None)

映射本身很朴素，功夫在前处理：把《上海土壤》里每个土属的典型质地剖面整理成 CSV，逐层读入，取剖面中最差一层的质地做最终归属。老资料里"中壤""重壤"这类卡庆斯基制名称，要先转换成 USDA 制的砂粒、粉粒、粘粒百分数再进三角图，两套制式的质地边界并不重合。粒径数据充足时，用 Saxton-Rawls pedotransfer 函数把粒径、容重、有机质换算成 Ks，再按 2.2 节的阈值分组，比纯质地映射更贴实测，但该函数标定数据主要来自美国土壤，对上海高粉粒、高盐基的冲积土外推时，建议留 20% 左右的余量。

注意：城市绿地里的"客土"（绿化工程换填的种植土），资料上常标为砂壤土，实际压实后入渗率往往跌到 C 组，遇到这种情况以现场探测为准。

3. 高精度栅格怎么造：数据源分工、投影统一与 rasterio 栅格化参数

3.1 先搞清楚"高精度"到底来自哪里

1:50 万上海市第二次土壤普查数字化图的精度在制图综合层面，最小图斑对应的地面范围远大于 10 m 像元；HWSD 是 30 弧秒网格，SoilGrids 是 250 m。直接拿任何一个单一来源做 10 m 输出，都只是把粗边界切细，叫"高精度"名不副实。常见做法是把来源按职责拆开，让各自干各自擅长的事：

这样组合出来的栅格，边界精度来自土壤图，内部细分来自地形与用地信息，才算真正对得起"高精度"三个字。

3.2 投影与像元对齐：统一到本地投影坐标系

上海市常用 CGCS2000 基准、3° 带高斯-克吕格投影，中央经线取 121.5°E，单位为米。生产上第一步是把所有矢量重投影到与模板栅格一致的坐标系，模板一般取土地利用栅格或 DEM。像元对齐有三条硬约束：像元尺寸取 10 m 的整数倍；像元原点（左上角）与模板完全一致；行列数不能四舍五入，否则整幅栅格会沿东北方向漂移半个像元，与路网叠合时错位肉眼可见。最省事的做法是让模板栅格先定义好 transform，后续所有栅格化和重采样都以它为唯一基准。

3.3 rasterize 参数与输出写法

import geopandas as gpd import numpy as np import rasterio from rasterio.features import rasterize # hsg_code: 1=A, 2=B, 3=C, 4=D, 9=扰动土, -9999=无数据 gdf = gpd.read_file("sh_soil_50w_hsg.shp") with rasterio.open("sh_lu_10m_template.tif") as ref: gdf = gdf.to_crs(ref.crs) # 统一坐标系 transform = ref.transform out_shape = (ref.height, ref.width) shapes = [(geom, int(code)) for geom, code in zip(gdf.geometry, gdf["hsg_code"])] raster = rasterize( shapes, out_shape=out_shape, transform=transform, fill=-9999, dtype=np.int16, all_touched=False, ) profile = ref.profile.copy() profile.update(dtype=np.int16, count=1, nodata=-9999, compress="deflate") with rasterio.open("sh_hsg_10m_raw.tif", "w", **profile) as dst: dst.write(raster, 1)

rasterize 的三个参数值得单独说。all_touched=False 表示只有像元中心落在图斑内才赋值，避免沿边界多占一圈像元；对 10 m 栅格和细碎图斑，先设 False，若边界处出现大量空洞再改 True。dtype 用 int16 不是因为分组只有 4 个值，而是要为 9（扰动土）和 -9999（无数据）留编码空间，同时保持与土地利用栅格一致的整型语义，后续 CN 查表时布尔索引的速度也更接近 C 级。fill=-9999 先铺底，栅格化后未覆盖区域全部是 -9999，后面统一用上海市域掩膜裁掉；不加掩膜的话，沿海潮间带和江面里的沉积物会被当成真实土壤算进 D 组。

3.4 边界修正与孔洞清理：gdal_sieve 的用法

原始栅格化结果有两个典型毛病：河岸边土壤图与水体相交，产生锯齿状伪像元；细碎图斑在 all_touched=False 下变成孔洞。处理分两步。先做水体掩膜，把黄浦江、长江口、淀山湖等水体多边形栅格化为 -9999；再做连通域清理，用 GDAL 自带的 sieve 工具去掉小于阈值的碎斑：

gdal_sieve.py -st 20 -8 sh_hsg_10m_raw.tif sh_hsg_10m_sieve.tif

-st 20 表示小于 20 个像元的连通域并入周边最大类，-8 表示八邻域连通判定。阈值按碎斑分布调：上海市区地块被道路切得很碎，20 太小会保留地块内部的图斑噪声，我一般给到 30～50。最后把城市建设用地边界内的像元用土地利用图覆写为 9 类扰动土，避免老图斑在建成区"复活"。sieve 只改空间连续性，不改变大类内部的真实过渡，对 10 m 栅格是安全的。

提示：如果下游是 SWMM 这类子汇水区模型，栅格分辨率取子汇水区平均尺寸的 1/4 到 1/2 即可，过细只增加重采样耗时，不会提高 CN 精度。

4. 上海土壤水文分组栅格的精度验证：样点布局、混淆矩阵与边界效应排查

4.1 样点布局：分层抽样比均匀撒点可靠

验证样本按 HSG 和地貌单元双重分层抽取。每个 HSG 组不少于 20 个点，全市 80～120 个点足以支撑四类混淆矩阵。按二项分布估算，期望总体精度 80%、允许误差 10%、置信水平 95% 时最少约 62 个点，和这个范围吻合。点位约束有三条：距图斑边界至少一个像元；避开道路硬化面；避开已竣工基坑。现场每点挖 60 cm 剖面读质地与层次，再用双环入渗仪测稳定入渗率；剖面判分组，入渗率用于解释偏差来源而不是直接改判——单点入渗受压实、根系和裂隙影响很大，一次读数不能代表整个图斑。

4.2 混淆矩阵、总体精度与 Kappa

把实测编码与栅格提取编码对齐后，用 sklearn 一次算出三个指标：

import numpy as np from sklearn.metrics import confusion_matrix, cohen_kappa_score obs = np.array([1, 1, 2, 2, 3, 3, 3, 4, 4, 4]) # 实测 HSG pred = np.array([1, 2, 2, 2, 3, 4, 3, 4, 4, 4]) # 栅格提取 cm = confusion_matrix(obs, pred, labels=[1, 2, 3, 4]) oa = np.trace(cm) / cm.sum() kappa = cohen_kappa_score(obs, pred) producer = np.diag(cm) / cm.sum(axis=0) # 生产者精度 user = np.diag(cm) / cm.sum(axis=1) # 用户精度

一组合格的 90 点验证结果如下：

总体精度 84.4%，Kappa 0.79。这个数放在 1:50 万源图衍生的 10 m 产品里是合格的；真正的诊断信息在混淆结构上：B 与 C 互相错分占了大头，说明质地过渡带被硬切成两档的边界误差，远大于"整块分错类"的误差。Kappa 低于 0.8 时不要急着改栅格，先回到 2.3 节的剖面最差层判定逻辑，复查原始属性表。

4.3 三类最容易出错的边界

河岸带是重灾区。黄浦江和吴淞江沿岸的砂质透镜体多是古河道摆动留下的，土壤图上没有，要用 DEM 的河漫滩相对高度和影像纹理补。做法是取距水系 100 m 缓冲带，叠加相对高程小于 2 m 的区域，把原判 C 的像元复核为 B。第二类是田埂和大棚菜地，犁底层和压实层让剖面整体下移半档，这类地类面积逐年缩小，但残余图斑对 CN 的影响不小。第三类是建设用地边缘，回填土与原地表犬牙交错，栅格上表现为 9 类像元包围 D 组的椒盐噪声。直接做 majority filter 会抹掉真实信息，我一般保留原值，在属性里加 flag 字段，模型侧对 9 类单独设 CN，而不是硬并进 C 或 D。

5. 把 HSG 栅格接进上海 SCS-CN 径流模拟：CN 查表合成、AMC 修正与反推验证

5.1 TR-55 的 HSG-CN 查表关系

HSG 单独不产生水文意义，意义在 SCS-CN 的查表环节。TR-55 给每个土地利用类别配了四组 CN 默认值，上海的常见类别如下：

这张表的工程潜台词是：土地利用对 CN 的影响远大于 HSG 分组。上海市区以高密度建成区为主，HSG 的误差只有在绿地、农田这类透水地类上才会放大，所以验证资源应当优先投到郊野和蓝绿空间。

5.2 栅格级 CN 合成的 numpy 写法

把 HSG 栅格与同期土地利用栅格对齐到同一像元网格后，按查表逐项赋值：

import numpy as np import rasterio CN_TABLE = { (1, 1): 39, (2, 1): 61, (3, 1): 74, (4, 1): 80, # 绿地 (1, 2): 89, (2, 2): 92, (3, 2): 94, (4, 2): 95, # 商服/工业 (1, 3): 98, (2, 3): 98, (3, 3): 98, (4, 3): 98, # 不透水 (1, 4): 67, (2, 4): 75, (3, 4): 83, (4, 4): 87, # 旱作农田 } with rasterio.open("sh_hsg_10m.tif") as hs, \ rasterio.open("sh_lu_10m.tif") as lu: hsg = hs.read(1) land = lu.read(1) cn = np.full(hsg.shape, -9999, dtype=np.float32) for (h, l), val in CN_TABLE.items(): cn[(hsg == h) & (land == l)] = val profile = hs.profile.copy() profile.update(dtype=np.float32, nodata=-9999, compress="deflate") with rasterio.open("sh_cn_10m.tif", "w", **profile) as dst: dst.write(cn, 1)

这里有一个覆盖顺序的坑：字典遍历的先后决定同名像元谁胜出。上面这张表各键互斥所以安全；若引入水域、湿地等多条件叠加的类别，必须把优先级最高的类别放最后赋值，或用 np.select 按条件数组顺序求值。

5.3 三个必调参数与一个反推验证技巧

第一个必调参数是前期土壤湿度条件 AMC。上海梅雨和台风季前期土壤接近饱和，AMC-III 的 CN 比 AMC-II 高 5～15，只出一版栅格会系统性低估产流。常见做法是准备两版 CN 栅格，用前 5 日降雨量判断当次取值，AMC-II 到 AMC-III 的换算用 SWAT 近似式 CN3 = CN2 * exp(0.00673 * (100 - CN2))。第二个是初损率 λ。TR-55 默认 0.2，上海城镇流域的短历时降雨用 0.05 与实测过程线拟合更好，S 的换算按 Hawkins 的 S05 = 1.33 * S02^1.15（英寸制）重算，不能直接套 0.2 的 S 值。第三个是聚合重采样。10 m 栅格聚合到子汇水区时只能取众数或面积加权平均，双线性插值会把四类整数插出 3.7 这样的伪类。

注意：凡 CN_TABLE 未覆盖的像元（nodata、9 类扰动土）最后都会落在 -9999 上，模型侧要按无效值剔除，不要用 0 兜底，0 会被当成 CN=100 参与计算。

验证技巧：挑上海近年来 3～5 场有完整雨量与出口流量记录的短历时降雨，反推流域平均 CN，再与同流域栅格 CN 的面积加权均值对比。偏差在 ±5 以内说明 HSG 栅格与土地利用编码的配合是自洽的；若偏差方向一致且普遍偏大，优先检查土地利用的 CN 取值和不透水率假设，而不是回头改土壤分组。

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