土壤可蚀性因子K是土壤流失方程的关键参数之一。K因子是侵蚀预报模型计算和水土流失动态监测的基础数据。K值的定义是在标准小区上单位降雨侵蚀力所引起土壤流失量的多少^[1]，反映土壤对侵蚀的敏感程度，其大小与土壤本身性质有关。

我国第1次水利普查时获得了全国范围首张K值分布图^[2]。由于实测K值使用的径流小区建造昂贵，花费时间长且缺少统一的测量规范，能收集到的K值数据十分有限，因此利用经验公式对K值进行估算^[3]。然而，使用的土种志数据距今已有近40 a，难以反映目前的土壤情况。为了获得更准确的K值结果，亟需在大范围内对计算K值的基础数据进行更新。

依靠有限样点很难满足越来越高的制图精度要求，因此还需探索一种能够实现高空间分辨率的K值点面扩展方法。以往研究在利用样点制作K值分布图时，常使用土壤图斑链接法及地统计插值^[4]。然而，图斑链接法导致K值在土壤类型边界发生突变，且不能反映同一土壤图斑内存在的K值变异性^[5]。使用地统计学方法对K值进行制图时，由于采样密度对插值精度影响较大，需空间结构因素大于随机部分才能得到较准确的结果，因此研究范围大多局限在流域尺度^[4]。随着数字土壤制图的发展，遥感数据已被应用于K值空间预测，例如高程作为协变量被引入克里格插值进行K值制图^[6]，高程及植被作为空间变量可以提高K值的预测精度^[7]。机器学习算法中的随机森林模型，是一种通过自助重采样技术对多个训练集随机选取特征变量，分别进行回归建模并平均多个预测值得到最终结果的方法^[8]。该方法已被成功运用于预测土壤有机质含量、质地与土层厚度等关键土壤属性^[9]。

为了实现全国土壤可蚀性因子更新目标，笔者提出一套基于近期汇编的中国土系志数据，利用随机森林模型进行点面扩展获得全国水蚀区K值分布图的方法，解决更新的2个关键问题：第一是如何获取全国地区最新的样点K值，第二是如何将样点的更新结果扩展至区域尺度。本次更新结果将有助于改善当前全国K值图数据陈旧及制图方法落后的问题。该成果可为土壤侵蚀预报模型提供基础数据，为年度水土流失动态监测工作提供服务。

中国土系志数据是我国自第2次土壤普查以来最新且系统性的土壤调查数据^[10]。全国土系调查项目于2008—2018年基于定量标准和统一分类原则开展系统性调查研究^[11]。笔者基于各省份土系调查资料，收集并整理4 327个采样点的土壤剖面数据，查阅其调查地点(经纬度)和理化性质。其中有机碳质量分数通过除以系数0.58(采用Van Bemmelen因子，假定土壤有机质的含碳质量分数为58%)统一转换为有机质质量分数^[12]。将数据进行预处理，主要采用加权平均的方法提取各样点0~30 cm耕层土壤的有机质和机械组成信息，当土体厚度 < 30 cm时则以实际厚度进行加权平均。

采用全国第1次水利普查时K值的计算方法，使用通用土壤流失方程中的诺模图(Nomo)计算K值(Nomo-K)^[2]。

式中：K为土壤可蚀性，t·hm²·h/(MJ·mm·hm²)；O为土壤有机质质量分数，%；S为土壤结构系数(计算土壤粒径几何平均直径后查表所得，量纲为1)；P为土壤渗透性等级(根据质地分类查表所得，量纲为1)；M为质地指数，由N₁和N₂计算所得，%；N₁为粒径在0.002~0.100 mm之间的土壤颗粒质量分数比例，%；N₂为粒径 < 0.002 mm的土壤颗粒质量分数，%。

由于全国土系调查中的质地分布仅包含砂粒(≥0.050~2.000 mm)、粉粒(≥0.002~0.050 mm)和黏粒(< 0.002 mm)分级，Nomo公式计算时需要把已有土壤颗粒分析结果转化为极细砂质量分数(≥0.050~0.100 mm)。本研究比较了3种插值方法，分别是自然对数线性插值法、三次样条函数法、三次样条函数结合自然对数。计算过程如下：首先将各级别土壤颗粒质量分数累加计算得到小于某一粒径的质量分数(< 0.002 mm、< 0.050 mm、< 2.000 mm)，并计算该粒径对应的自然对数(ln(0.002)、ln(0.050)和ln(2.000))。自然对数线性插值法可建立颗粒累加质量分数与粒级自然对数之间的线性关系；三次样条函数法可建立颗粒累加质量分数与粒径之间的三次样条函数；三次样条函数结合自然对数法则可建立颗粒累加质量分数与粒径自然对数之间的三次样条函数。基于以上方法得到了 < 0.100 mm粒径所对应的质量分数，并将 < 0.050 mm和 < 0.100 mm粒径质量分数相减得到极细砂质量分数，通过分析极细砂质量分数的合理性和比较插值方程R²确定最优的插值方法。

当有机质质量分数＞12%时，使用Nomo公式可能会导致负值出现，因此对有机质>12%的样点采用EPIC模型计算K值^[3]

式中：S_a为粒径≥0.050~2.000 mm的砂粒质量分数，%；S_i为≥0.002~0.050 mm的土壤粉粒质量分数，%；C_l为粒径 < 0.002 mm的土壤黏粒质量分数，%；C为土壤有机碳质量分数，%，通过有机质乘以转换系数58%所得；S_n1为土壤碳酸钙质量分数，%。

由于Nomo和EPIC公式的计算过程不同，采用EPIC公式得到的K值(EPIC-K)不能直接替代Nomo-K值。因此，笔者使用线性、指数、幂函数和3阶多项式函数分别建立2种公式计算K值的转换关系。比较4种拟合方程的R²，除线性方程拟合精度略低以外(0.795 5)，其他3种非线性方程的拟合精度十分接近(0.806 4~0.807 2)。从趋势来看，Nomo-K和EPIC-K值之间并非简单线性关系，随着EPIC-K值的增大，Nomo-K值出现先缓慢增加后快速增加的趋势。这与指数方程所显示的规律是一致的。总之，当土壤有机质≤12%时，使用Nomo公式计算K值；当有机质＞12%时，采用EPIC公式计算K值并利用其与Nomo-K值的指数关系式进行修正。

沿用在太湖流域片预测K值时使用的遥感方法获取的环境因子指标作为预测变量，包括气候、地表温度、植被、光波、地形和母质^[13]。由于随机森林模型具有避免过拟合的优势，且模型训练集样本量足够大，没有对环境因子进行剔除。

各项环境指标的空间分辨率从30 m到1 km不等，通过最邻近法将所有指标统一重采样为30 m分辨率后进行叠加运算。由于气候和地表温度在较小尺度下不存在大的差异，重采样后在1 km范围内气候及地表温度值不变符合实际情况。各样点对应的环境变量值通过坐标位置提取。

采用随机森林回归方法建立K值预测模型。首先对样点的环境要素及K值进行训练，然后利用多源遥感影像将建立起K值预测模型推广至全国地区，完成K值由点及面的反演计算。其中模型的2个关键参数通过调参选择：模型误差随着决策树数量(n_tree)的增加而降低，当决策树数量超过400以后，模型误差趋于稳定，因此选择使用默认参数n_tree=500；随机森林节点变量数(m_try)一般选择环境变量数的1/3(本研究对应m_try=20)，通过对节点变量数在m_try=20左右逐个比较模型精度，确定最佳的m_try值为17。

K值预测模型的建模精度通过外部验证的方法进行评估，采取常用的5折交叉验证法，重复100次取平均值作为评价结果。评价指标包括决定系数(R²)，均方根误差(E₁)以及平均绝对误差(E₂)，计算公式如下：

式中：n为训练样本个数；Y_i为观测K值；$\hat{Y}_i$为模型预测K值；Y_i为平均观测K值。

基于K值预测模型，以61项环境因子空间分布图作为自变量，使用随机森林方法进行预测，输出全国尺度K因子栅格图，K值更新图的空间分辨率为30 m×30 m。

研究利用Excel 2019计算样点Nomo-K值与EPIC-K值，并建立两者的回归关系；利用ArcGIS 10.8处理空间数据，包括投影转换、空间重采样，提取样点对应的环境变量，统计栅格极值、平均值及标准差等；利用R语言编程实现土壤粒径的批量插值转换，同时构建随机森林回归模型，进行全国范围K值的空间预测。在分析K值分布规律时，土壤类型图使用1995年编制的《1∶100万中华人民共和国土壤图》。

插值获得诺模图计算所需极细砂的质量分数，以江西省关山系样点为例，使用3种插值方法得到粒径质量分数结果如表 1所示。其中，自然对数线性插值法拟合精度较差，获得的 < 0.100 mm粒径质量分数低于 < 0.050 mm粒径，导致极细砂质量分数为负值；三次样条函数法受样本数量过少(n=3)的影响，获得的累加极细砂粒径(< 0.100 mm)质量分数高于100%，明显不合理；三次样条函数结合自然对数法克服了函数分布不合理的缺点，以高精度方程拟合了不同土壤粒径的质量分数(R²≈1)，同时避免了质量分数高于100%的情况(图 1)。

表 1(Tab. 1)

表 1 土壤质地插值转换结果对比(江西关山系) Tab. 1 Comparison of soil texture interpolation conversion results(Guanshan series in Jiangxi) 土壤粒径 Soil particle size/mm 自然对数法线性插值 Linear interpolation with natural logarithm/% 三次样条函数法 Cubic spline function/% 三次样条函数结合自然对数 Cubic spline function combined with natural logarithm/% < 2.000 100.0 100.0 100.0 < 0.100 78.2 119.1 85.9 < 0.050 80.1 80.1 80.1 < 0.002 40.8 40.8 40.8 注：< 0.100 mm为预测的土壤粒径质量分数，其余为已知土壤粒径质量分数。Notes: < 0.100 mm is the predicted soil particle sizemass fraction, and the rest is the known soil particle size mass fraction.

表 1 土壤质地插值转换结果对比(江西关山系) Tab. 1 Comparison of soil texture interpolation conversion results(Guanshan series in Jiangxi)

图 1(Fig. 1)

圆点为已知粒径质量分数，三角形为预测粒径质量分数。下同。 The dot represents the known particle size mass fraction, and the triangle represents the predicted particle size mass fraction. The same below. 图 1 土壤质地插值转换方法对比 Fig. 1 Comparison of soil texture interpolation conversion methods

在少数样点中，土壤 < 0.050 mm和 < 2.000 mm的粒径质量分数十分接近，使用三次样条函数结合自然对数法计算所得 < 0.100 mm的粒径结果可能略高于100%。针对这种情况，以山东省小岛河系样点为例(图 2)，利用线性函数在 < 0.050 mm和 < 2.000 mm分段之间插值得到 < 0.100 mm的粒径质量分数，线性插值使 < 0.100 mm粒径质量分数既小于 < 2.000 mm粒径质量分数，也大于 < 0.050 mm粒径质量分数，能确保结果的合理性。

图 2(Fig. 2)

实线为三次样条函数结合自然对数法拟合，虚线为自然对数线性插值法。 The solid line is cubic spline function combined with natural logarithm, and the dotted line is linear interpolation with natural logarithm. 图 2 土壤粒径插值结果过高处理 Fig. 2 Processing of excessive soil particle size interpolation results

采用Nomo和EPIC公式分别计算4 327个土系样点K值，统计结果如表 2所示。EPIC计算K因子的最小值比Nomo-K值大，最大值比Nomo-K值小，且标准差更小，EPIC计算K值的平均值及中位数比Nomo-K值大16%。因此，不能直接使用EPIC-K值替换高土壤有机质条件下的Nomo-K值。

表 2(Tab. 2)

表 2 全国土壤可蚀性K值计算结果统计 Tab. 2 Statistics of calculated results of soil erodibility K value nationwide 计算公式 Formula 最小值 Minimum 最大值 Maximum 平均值 Mean 中位数 Median 标准差 Standard deviation Nomo -0.192 6 0.098 7 0.032 5 0.032 4 0.017 0 EPIC 0.005 4 0.065 1 0.038 7 0.038 9 0.010 1

表 2 全国土壤可蚀性K值计算结果统计 Tab. 2 Statistics of calculated results of soil erodibility K value nationwide

建立2种公式的转换关系时，R²均>0.79(图 3)，但散点图显示随着EPIC-K值的增大，Nomo-K值出现先缓慢增加后快速增加的趋势，这与指数方程的规律相符。最终选择该方程对全国有机质＞12%的55个样点进行K值修正。这些样点主要分布在东北和青藏高原的湿润或者滞水条件下，采样部位在上、中、下坡和台地均有分布，表明并非沉积导致K出现负值。根据定义，K因子为在标准小区上由单位降雨侵蚀力所引起土壤流失量的多少^[1]。土壤流失量不可能为负数，所以K值不应存在负值的情况。校正前55个样点的EPIC-K均值为0.034 4，校正后均值为0.026 6，弥补了EPIC-K值大于Nomo-K值的差异。

图 3(Fig. 3)

EPIC-K is soil erodibility calculated by EPIC model. Nomo-K is soil erodibility calculated by Nomo formula. 图 3 土壤可蚀性K值计算结果拟合方程对比 Fig. 3 Comparison of fitting equations for the calculated results of soil erodibility K

以土系调查项目所建立的更新K值作为因变量，以包含气候、地形、植被、母质类型、地表温度、可见光及近红外波段在内的61项环境因子指标作为自变量，在m_try=17，n_tree=500的模型参数下，建立的全国K因子随机森林预测模型，运用5折交叉验证方法重复100次的评估结果为：R²=0.381，E₁=0.012 47，E₂=0.009 63。这与太湖流域片建立K值模型的精度类似^[13]，但低于使用随机森林预测土壤有机质和全氮的空间分布^[14]。这与本研究采样密度稀疏、空间建模尺度大以及缺乏与K值强相关性的预测变量等原因有关。

全国K因子分布范围在0.005 1~0.074 5之间，这与前人在样点尺度建立的我国主要土壤K因子分布范围0.000 8~0.070 5相似^[15]。分布趋势呈现出华北平原和黄土高原地区高，天山地区及东北平原次之，南方山地丘陵、青藏高原及内蒙古高原较低的规律。这与其他学者总结的我国土壤可蚀性分布规律相似^[16]，但该研究发现青藏高原的K值最高。也有研究发现青藏高原K值比黄土高原、川渝地区和滇东北地区的K值略低^[17]，这可能是由于青藏高原K值存在较大差异所造成的。

进一步利用水土保持一级区划分析全国更新K值的分布规律，发现更新K值在西北黄土高原区和北方土石山区较大，在青藏高原区最小(表 3)，这与自然及人为活动等因素的综合作用有关。其中，西北黄土高原区土壤有机质含量普遍很低因此K值偏高^[18]；北方土石山区长期以来受人类耕作活动的影响导致K值较高；青藏高原区由于地形落差大，K值具有垂直地带性，在少量低海拔河谷地带K值较高，而在大面积的高海拔地区由于成土作用弱K值明显降低^[17]。此外，东北黑土区在风砂土分布地区K值最大，在黑土集中的漫川漫岗地带由于农耕活动导致K值较大，暗棕壤分布的大兴安岭及长白山等地K值较小。

表 3(Tab. 3)

表 3 土壤可蚀性更新K值在水土保持一级区划的统计 Tab. 3 Statistics of soil erodibility update K value in the first level zoning of soil and water conservation 水土保持一级区划 Primary division of water and soil conservation 最小值 Minimum 最大值 Maximum 平均值 Mean 标准差 Standard deviation 东北黑土区Northeast black soil region 0.011 51 0.049 22 0.029 05 0.005 36 北方风沙区Northern sandstorm region 0.005 08 0.064 84 0.028 04 0.007 14 北方土石山区Northern earth-rock mountain region 0.011 94 0.074 53 0.038 79 0.008 60 西北黄土高原区Northwest Loess Plateau region 0.013 05 0.059 73 0.039 17 0.007 67 南方红壤区Southern red soil region 0.011 55 0.070 63 0.030 25 0.006 60 西南紫色土区Southwest purple soil region 0.017 44 0.057 01 0.034 65 0.004 65 西南岩溶区Southwest karst region 0.015 74 0.044 06 0.030 17 0.002 32 青藏高原区Qinghai-Tibet Plateau region 0.009 89 0.053 37 0.025 54 0.004 57 注：水土保持区划使用2012年水利部印发的《全国水土保持区划(试行)》一级区范围。Notes：The water and soil conservation zoning uses the scope of the first level of the National Water and Soil Conservation Zoning (Trial) issued by the Ministry of Water Resources in 2012.

表 3 土壤可蚀性更新K值在水土保持一级区划的统计 Tab. 3 Statistics of soil erodibility update K value in the first level zoning of soil and water conservation

K值的空间分异与我国主要土壤类型的分布情况有关，前人曾报道黄土高原的黄绵土和东北的黑土K值较大，南方红壤K值较小^{[15, 19]}。本研究从主要土壤类型统计来看，黄绵土和潮土的更新K值较大，平均值均超过0.04，而高山土的平均K值最小，其次为栗钙土、暗棕壤和红壤，均不超过0.03(表 4)。分析原因认为黄绵土和潮土的粉粒含量较高，土壤更容易遭受侵蚀。此外，潮土K值的标准差最大。前人研究发现，同一种土壤下不同母质和植被覆盖可以通过影响土壤的水稳性从而改变土壤的抗蚀能力^[20]，因此相同土壤类型下K值的差别也可能很大。

表 4(Tab. 4)

表 4 主要土壤类型下土壤可蚀性更新K值统计 Tab. 4 Statistics of updated K values of soil erodibility under main soil types 土壤类型 Soil type 最小值 Minimum 最大值 Maximum 平均值 Mean 标准差 Standard deviation 红壤Red earths 0.015 24 0.045 45 0.028 72 0.003 82 黄壤Yellow earths 0.017 61 0.044 19 0.030 48 0.003 52 紫色土Purplish soils 0.014 89 0.055 61 0.031 22 0.003 06 棕壤Brown earths 0.014 77 0.059 01 0.033 01 0.005 12 暗棕壤Dark-brown earths 0.011 51 0.047 05 0.028 59 0.003 52 黑土Black soils 0.020 25 0.048 90 0.035 79 0.004 56 黑钙土Chernozems 0.012 32 0.050 23 0.030 03 0.005 88 栗钙土Castanozems 0.014 54 0.054 69 0.027 80 0.007 33 褐土Cinnamon soils 0.018 91 0.063 53 0.039 54 0.007 28 潮土Fluvo-aquic soils 0.006 59 0.074 53 0.041 61 0.009 23 黄绵土Cultivated loessial soils 0.022 55 0.056 49 0.044 53 0.005 48 草甸土Meadow soils 0.005 69 0.063 19 0.031 70 0.006 20 高山土Alpine soils 0.009 89 0.051 47 0.025 86 0.004 74 盐碱土Saline soils 0.005 26 0.068 93 0.031 27 0.008 48 灌淤土Cumulated irrigated soils 0.021 40 0.055 85 0.037 21 0.004 13 水稻土Paddy soils 0.012 16 0.070 48 0.032 85 0.006 38 注：土壤类型指代“土壤发生分类”系统中的土类，其中高山土包括草毡土、黑毡土、寒钙土、冷钙土、棕冷钙土、寒漠土、冷漠土、寒冻土，盐碱土包括盐土、漠境盐土、滨海盐土、酸性硫酸盐土、寒原盐土和碱土。Notes：1 The soil type refers to the soil category in the “soil genesis classification” system, in which the alpine soils includes felty soils, dark felty soils, frigid calcic soils, cold calcic soils, cold brown calcic soils, frigid desert soils, cold desert soils, and frigid frozen soils, and the saline soils includes solonchaks, desert solonchaks, coastal solonchaks, acid sulphate soils, frigid plateau solonchaks and solonetzs.

表 4 主要土壤类型下土壤可蚀性更新K值统计 Tab. 4 Statistics of updated K values of soil erodibility under main soil types

值得注意的是，由于Nomo和EPIC模型是基于美国土壤实验所建立的，在中国地区存在适用性问题。张科利等^[19]利用径流小区观测资料在我国东部地区建立了K值的校正公式，分别为K=-0.033 36+0.744 88K_Nomo(R=0.721)，K=-0.013 83+0.515 75K_EPIC(R=0.613)；更进一步说，Nomo模型在黄土高原和黑土区的应用效果(R²=0.55~0.60)要优于红壤和紫色土(R²=0.11~0.37)，EPIC模型除黑土区以外(R²=0.05)在其他地区的应用效果较好(R²>0.58)^[21]。直接应用本研究结果可能会对K值产生不同程度的高估，因此需要在全国范围内进一步收集实测资料将K值校正以后用于土壤侵蚀模数计算。

1) 本次更新进一步完善了全国土壤可蚀性因子计算方法。本研究在美国制砂粒、粉粒、黏粒机械组成实测值的基础上插值得到极细砂质量分数，基础数据更为准确；通过计算土壤粒径几何平均直径来获取结构系数，更具有客观性；建立于全国统一的土系基础数据库，消除了不同省份之间的系统误差；以30 m分辨率栅格为制图单位，更为详细地表达K值的空间变异情况。

2) 全国土壤可蚀性更新K值的空间分布呈现出以下规律：黄土高原区和北方土石山区居高，西南紫色土区、南方红壤区、西南岩溶区、东北黑土区和北方风沙区次之，青藏高原区最小。这主要是由于不同土壤类型在各地区的分布所造成的。统计可知黄绵土和潮土的K值较大，红壤、暗棕壤、栗钙土的K值较小，高山土的K值最小。

3) 本研究采用的Nomo和EPIC公式都是基于美国土壤建立的K值计算经验公式，可能与实际监测的K真实值有差距。因此需要积累我国不同地区标准径流小区的长期观测资料，使用年平均侵蚀量和降雨侵蚀力因子进一步率定K值，以期建立起实测值与经验公式计算值之间的关系式或直接建立适用于我国的土壤可蚀性计算公式。