ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

90米土壤可蚀性因子K数据:从原理到RUSLE/CSLE建模实操指南

90米土壤可蚀性因子K数据:从原理到RUSLE/CSLE建模实操指南 做水土流失调查的同行应该都有过这种经历模型参数表里明明列着“土壤可蚀性因子K”可真到给自己项目区填一个像样的K值时翻文献、查数据库、找土壤图忙活一上午也拿不准。K因子不是简单查一张表就能得到的参数它由砂粒、粉粒、黏粒含量、有机碳、土壤结构、渗透性共同制约而不同来源的数据差异还不小。最近我试用了“中国90米分辨率可蚀性因子K数据”整体用下来比以往常见的1公里产品顺手不少尤其在中小流域尺度做RUSLE、CSLE建模时省去了大量自算K值的功夫。这篇文章打算从K因子的物理含义、这套数据的价值、算法原理、实操流程、常见坑五个角度拆开讲清楚适合正在做土壤侵蚀建模、水土保持规划、生态评价的科研人员和基层技术员参考。1. 先弄明白土壤可蚀性因子K是什么为什么重要1.1 一个定义背后藏着“标准小区”的故事K因子在通用土壤流失方程USLE中的官方定义是在标准小区上单位降雨侵蚀力所引起的土壤流失量。标准小区的条件很苛刻——坡长22.1米坡度9%顺坡耕作地面裸露且无覆盖每年翻耕。之所以限定这么死板的条件是为了把“土壤本身抵抗侵蚀的能力”单独拎出来排除坡长、坡度、植被、措施的影响。你可以把K因子理解成土壤的“体质报告”同样是下一场暴雨有的土壤表面很快形成结皮径流带着泥沙哗哗往下冲有的土壤颗粒粗、透水快雨水直接渗下去流失量就小。K值就是把这种差异量化成数字。通常来说粉粒含量高的土壤K值偏高因为粉粒既不象砂粒那样粗重也不象黏粒那样容易团聚最容易被径流带走砂质土虽然颗粒松但渗透性强K值反而中等高有机质、结构好的土壤K值偏低因为团聚体稳定不易分散。土壤侵蚀模型发展到今天USLE、RUSLE、RUSLE2、CSLE换了多少代K因子始终是必须的乘法因子。它的地位就像炒菜里的盐——比例不大但没有它整道菜就不成立。在CSLE方程里K因子同样放在公式的乘法项中和降雨侵蚀力R、地形因子LS、生物措施因子B、工程措施因子E、耕作措施因子T一起共同决定土壤流失量A。1.2 K因子在模型链条里的位置以最常用的RUSLE为例年土壤流失量A R × K × L × S × C × P。这里每个因子代表一类控制手段R是老天爷给的降雨能量K是土壤本身的脆弱程度L和S是地形地貌C是植被覆盖P是水土保持措施。六个因子相乘任何一个为0流失量就是0任何一个翻倍流失量也翻倍。K因子作为“本底脆弱性”不随季节、作物、管理措施变化所以它是区域土壤侵蚀评价里最稳定的输入之一。做项目时我常跟人打比方R因子是“有多大的雨在冲”K因子是“脚下的土有多不禁冲”。同一个降雨条件下同样坡度同样植被K值0.05的土和K值0.02的土流失量能差出两倍多。因此K因子的精度直接决定整体预报结果的可靠程度。很多初学者把精力全放在调C因子和P因子上却忽视K因子本身可能带来的误差这其实是本末倒置。1.3 量纲与单位陷阱K因子的单位是模型使用中最容易翻车的细节。国际制单位是t·ha·h/(ha·MJ·mm)美制单位是ton·ac·h/(ac·ft-tonf·in)。两者数值相差约13.17倍国际制数值约为美制数值乘以0.1317。国内大多数文献和指南要求使用国际制单位数值通常在0.005到0.09之间美制单位则在0.01到0.07之间。单看区间存在重叠不仔细辨别很容易弄混。判断一个K因子栅格用的是哪个单位最直接的办法是看数值分布直方图如果大量像元集中在0.01到0.05可能是美制如果集中在0.005到0.03更可能是国际制。另有一个笨办法——查数据元数据或随附说明文件正规数据一定写清楚单位。万一单位搞错后续侵蚀量计算结果就会系统性偏差十倍以上这种错误在论文里是致命的。2. 这套90米数据“值”在哪2.1 从1公里到90米意味着什么过去很多全国尺度的K因子产品空间分辨率为1公里甚至10公里做省级、国家级宏观评估没问题但一旦把研究区缩小到某个县、某个小流域1公里栅格往往连一条沟谷都分不清更谈不上反映不同土种的空间镶嵌。90米分辨率相当于把全国划分成约3弧秒的网格这个尺度能较好匹配SRTM等90米DEM产品也和当前多数区域水文模型的网格设置兼容。90米能带来的实际差异我在一个西南丘陵区的项目里体会很深。1公里产品在该区域几乎全部是均一数值完全看不出坡脚冲积土和坡顶残积土的差别换成90米数据后河谷两侧的冲积土、低丘顶部的风化壳、山腰的黏土母质层都能区分出来K值从0.015到0.045的渐变过程清清楚楚。这种空间细节对判断侵蚀热点区非常重要——同一面坡上不同土层的可蚀性差异可能比不同坡面之间的差异还大。2.2 面向模型应用的数据设计思路这套90米数据在设计上有几个值得点赞的地方。首先是坐标系和网格定义清晰数据处理时能直接和主流DEM产品对齐不需要反复投影重采样其次是像元值已经是国际制单位且经过掩膜处理水体、建筑区、裸岩等非土壤区域被设为无效值避免了后续计算里数学上有意义但实际无意义的数字干扰。更重要的是它是从土壤属性数据出发经过K因子计算公式逐像元计算得到的而不是从某张小比例尺土壤图直接矢量化转栅格。这意味着栅格之间的过渡是渐变的更接近自然土壤属性的连续分布特征用在模型里不会出现明显的“块状伪影”。2.3 适用场景与边界这套数据适合的场景可以列得很具体县域或小流域尺度的水土流失动态监测、全国水土保持区划的辅助分析、生态修复工程的选址评估、土壤侵蚀敏感性评价、高校和研究机构的模型教学。凡是研究区在几百平方公里到几万平方公里之间、模型网格在30米到250米之间这套90米数据基本都能直接或经重采样后使用。但也要说清楚边界。它不适合田块尺度的精准水土保持设计——你不可能用它指导一个具体梯田地块的施工因为90米网格内往往包含多种土地利用和土壤类型的混合。同样如果研究区地形极其破碎、属于典型的山区细碎地貌90米分辨率也可能忽略微小地形对土壤属性的再分布影响。遇到这种情况建议以这套数据为背景叠加野外采样点实测K值做局部校正。3. K因子是怎么算出来的从土壤属性到栅格3.1 两条主流算法路线K因子的计算方法大致分两派一是经典诺模图法二是在此基础上发展的公式法。诺模图法由Wischmeier等人提出需要四个输入土壤质地等级、有机质含量等级、土壤结构等级、土壤渗透性等级。操作方式是把这些等级值放到专门的诺模图上连线交点就是K值。优点是直观、考虑了土壤结构和渗透性的综合影响缺点是高度依赖专家经验把连续属性离散成等级空间化困难在大范围制图时不现实。公式法以EPIC公式最为常用。它把K因子表达为砂粒含量SAN、粉粒含量SIL、黏粒含量CLA和有机碳含量C的函数K 0.1317 × [0.2 0.3·exp(-0.0256·SAN·(1 - SIL/100))] × [SIL/(CLA SIL)]^0.3 × [1.0 - 0.25·C/(C exp(3.72 - 2.95·C))] × [1.0 - 0.7·SN1/(SN1 exp(-5.51 22.9·SN1))]其中SN1 1 - SAN/100SAN、SIL、CLA、C均以百分数表示。等式最前面的0.1317是把美制结果转换到国际制的系数如果只想要美制单位可去掉。这个公式看起来吓人实际计算用GIS栅格计算器几分钟就完成而且输入参数只有四个百分数特别适合大范围自动化制图。3.2 输入数据从哪里来做全国尺度的K因子数据最核心的输入是土壤质地和有机碳分布。国内主流做法是基于全国土壤调查形成的土种志、剖面数据和土壤类型图把每个土种典型剖面的质地数据、有机碳数据整理成属性表再通过空间连接和插值方法推广到整个面状分布区域。这个过程的核心难点在于属性数据是“点”而我们需要“面”点与面之间的对应关系靠土壤类型图来桥接。土壤质地数据到90米网格的插值不是简单的反距离权重而要先根据成土母质、地形部位、地貌类型建立分区控制再在分区内部做平滑。否则很容易出现“同一个土种内部K值完全一样边界处突然跳变”的阶梯效应。有机碳含量则往往做深度加权到表层30厘米因为侵蚀过程主要发生在表土这也是USLE系列模型里的默认约定。有机碳和有机质之间还有一个高频换算点很多土种志记录的是有机质含量而EPIC公式需要有机碳两者换算系数是1.724也就是有机碳 有机质 / 1.724。不少初算K值的人栽在这个系数上直接用有机质代入公式算出的K值偏小不说空间趋势也会被扭曲。3.3 数据生产中的质量控制点一个能全国发布的K因子产品背后必须经过若干质量检查。常见方法一是看数值范围全国尺度K值国际制合理区间一般在0.005到0.08如果出现大量超过0.1的像元说明输入数据或公式使用有问题。二是看空间分布是否符合土壤地理学常识比如黑土区、高有机质地区K值应该偏低沙性强的地区K值波动大黄土母质发育的土壤因粉粒含量高而K值偏中等偏上。三是抽样验证选取若干典型土壤剖面对比计算结果和诺模图法或文献实测值误差控制在可接受范围。数据发布时还应注意把无效值统一设为固定NoData值并在元数据里写明坐标系、单位、有效值范围、算法版本。这些看起来琐碎但对使用者来说元数据就是生命线。4. 拿到数据后的实操从栅格到侵蚀量4.1 第一步把坐标系和单位盘清楚拿到这套数据后不要急着做计算先花十分钟做三件事查看栅格属性里的坐标系、统计直方图、查看无效值编码。坐标系决定你后续叠加的其它图层是否需要重投影直方图判断单位无效值编码决定你需不需要做掩膜。实操中用QGIS或者ArcGIS皆可。打开属性表后看“源”或“栅格信息”页签记录投影类型。如果是地理坐标系WGS84像元尺寸会显示为0.000833度左右约等于3弧秒如果是投影坐标系如Albers等积投影像元尺寸会显示为90米左右。两种情况都能用但我个人更建议后续统一到研究区常用的投影坐标系特别是当研究区跨度大时Albers或Lambert等积投影比经纬度网格更适合面积计算。统计直方图方法更简单加载栅格后右键查看直方图看横轴范围。国际制K值一般不超过0.1美制一般不超过0.08两者默认范围看起来差不多所以要结合参考土的典型值综合判断。我经验里最稳妥的做法还是直接查元数据直方图只能作为佐证。4.2 裁剪、重采样与无效值处理按流域边界裁剪是常规操作。用多边形矢量边界裁剪K因子栅格时注意选择“提取掩膜”而不是“裁剪到几何范围”前者按边界精确提取后者只按外接矩形切一刀边界外会留一圈多余像元。裁剪后务必做一次无效值检查尤其当流域边界和无效值区域重合时被切出来的空洞在后期计算里会传染给结果图层。重采样最常见的问题是放大到更细分辨率。有人为了和5米DEM对齐把K因子从90米重采样到5米这是典型的资源浪费——K因子的真实空间变异尺度根本达不到5米重采样产生的细节全是虚假信息。如果必须统一分辨率建议重采样到25米或30米就足够重采样方法选“最近邻”或“双线性”均可但对类别属性强烈建议用最近邻双线性会在属性边界处产生中间值导致不存在的过渡土壤类型。无效值处理也要讲究。如果后续要做流域平均K值建议先把无效值区域单独提取出来作为掩膜参与统计的像元只保留有效值区域如果做栅格乘除计算则保持无效值编码不变待全部因子计算完毕再做统一掩膜否则某一层的空洞会污染所有结果。4.3 结合R、K、LS、C、P因子的完整计算K因子只是公式中的一个乘数实际操作时要把所有因子图层对齐到同一个像元网格。推荐工作流是这样的先选定统一的网格基准比如90米或30米以此为基础重采样R、K、LS、C、P各图层保证所有栅格的像元大小、范围、坐标系完全一致然后再用栅格计算器执行乘法。表达式以栅格计算器为例A R_factor * K_factor * LS_factor * C_factor * P_factor这里要注意两个细节。一是所有因子栅格的无效值编码必须统一否则计算时A的无效值范围会扩大二是输出结果不要直接截断为整数先保留浮点数最后再按需求取有效数字。另外一个常见问题是C因子和P因子经常被做成单值栅格即整幅图一个值这在计算时会放大K因子空间格局的影响因此K因子数据的精度对最终结果的空间分布有决定性影响。计算完毕后建议做一次合理性自检输出结果A的单位是t/(ha·a)对照研究区多年平均土壤侵蚀模数的文献值看数量级是否一致。如果单位或K因子转换错误A值往往会偏离一到两个数量级一眼就能发现。5. 常见问题速查表与实用避坑心得5.1 六个高频问题问题现象可能原因解决办法K值范围出现0.1以上的像元美制单位未换算乘以0.1317转为国际制裁剪后出现大面积空洞矢量边界与无效值区域重叠重新检查无效值掩膜先填洞再裁剪重采样后K值分布变形使用了错误的插值方法改用最近邻避免双线性平滑K因子和DEM范围对不齐投影或像元范围不一致统一重投影到目标坐标系再对齐网格计算结果A值比文献大十倍K因子单位或R因子单位存在问题逐一核对各因子量纲K值最可疑区域内K值完全均匀90米格网未有效反映属性变化检查原数据是否被重采样过低分辨率考虑局部校正5.2 几个值得单说的坑第一K因子图层不要反复重投影。每投影一次重采样插值都会改变原始像元值分布尤其是从等积投影转到高斯投影再转回等积投影来回折腾后K值直方图可能面目全非。建议把所有数据统一到目标坐标系后再计算全程只做一次投影处理。第二不要盲目用“全国均值”填到小区域。K因子的空间异质性极强即使同一县域内不同成土母质和地形部位的K值差异可能超过一个数量级。用90米数据时务必先查看研究区内的K值分布再做区域统计不要偷懒取全图平均。第三K因子不是越精确越好它需要和模型其他因子的精度匹配。如果你的R因子只有月尺度数据、C因子用植被指数粗略估算那么把K因子抠到极致也没有意义整体误差还是由最粗糙的那个因子决定。合理的做法是把K因子精度做到中等偏上把更多精力放在时间变异性更强的R和C因子上。第四关于90米数据与局部实测的衔接。我通常会在研究区内布设少量采样点采集表土样品测质地和有机碳用EPIC公式算点上的K值再与栅格值比较。偏差在20%以内的像元比例如果占多数说明数据质量可靠可以放心用于建模偏差过大则需要考虑局部土壤母质的特殊性和插值平滑的影响以实测点数据为控制进行局部校正。5.3 一个实用的检查技巧分享一个我常用的快速检查方法将K因子栅格与研究区土壤类型图叠加用分区统计工具按土壤类型分组统计K值的均值、标准差、最大值、最小值。如果某个土种的K值变异系数明显高于其它土种往往提示该区域的土壤属性插值存在问题或者这个土种分布跨越了多种成土母质区。这个方法不需要额外数据几分钟就能定位数据异常区域值得一试。这套90米可蚀性因子K数据对我而言最大的价值不是“省了算K因子的时间”而是给了我和DEM、土地利用数据天然对齐的标准化输入让模型跑起来更顺结果也更可控。在实际使用中我习惯先把K因子栅格和LS因子栅格叠加看一眼侵蚀热点一般就集中在K值偏高又恰好是高坡度的地方这个预览动作能让我在正式计算前就对结果有个预期省掉不少回头查问题的功夫。希望这篇文章能帮到正在为K因子发愁的你。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进