土壤侵蚀是指在自然界多种外力的作用下，部分土壤被搬运离开原先位置的过程。土壤侵蚀尤其是水蚀作为最主要的土地退化类型之一，会引起沟壑加剧、土层变薄、土壤肥力和质量下降、耕地面积不断减小等危害。重庆市奉节县是三峡库区重庆段中“山地，农业经济为主”地区的典型代表，区域内水土流失严重，生态脆弱敏感。为防止当地陷入发展停滞与资源破坏的恶性循环^[1]，对土壤侵蚀现象必须予以高度的重视和采取有效措施加以防治。

模型代数法是传统计算土壤侵蚀速率的方法，在过去已经发展出数十种模型，其中以RUSLE模型、CSLE模型等在国内外最为常见^[2-5]。这类模型通常根据一定的先验知识，人为地给模型因子附加不同的参数。基于模型代数法的土壤侵蚀计算已有不少学者研究，但由于土壤侵蚀具有明显的地域差异性，不同条件下土壤侵蚀形成机制和各驱动因子的主导作用均有差异，土壤侵蚀模型的参数难以确定，使得基于模型的土壤侵蚀计算方法精度较低。机器学习方法的发展给土壤侵蚀计算带来新的思路。相比于传统模型代数法，机器学习不需要明确土壤侵蚀演变过程的复杂多变关系^[6]，而是依据主要驱动因子数据和对应点实测数据来模拟量化驱动因子与土壤侵蚀之间的关系。它还可以充分挖掘遥感影像信息，从中将土壤侵蚀与地表土壤、土地利用、植被覆盖等地表地物结合起来，在最少人为干预的情况下对未知数据做出准确的预测。如黄晨璐^[7]基于随机森林(random forest, RF)模型实现黄土高原地区的土壤侵蚀速率的计算，最终R²为0.73。Sahour等^[8]在伊朗北部的卡西立安分水岭使用土壤侵蚀钉测量年侵蚀率，分别采用多元线性回归模型，RF、神经网络对土壤侵蚀数据进行训练，研究发现RF表现最优，在验证集中R²达到0.92。Rakhohorib等^[9]以印度西孟加拉邦和比哈尔邦之间的索巴流域为研究区，将土壤侵蚀速率分为低，中，高和极高4个等级，通过支持向量机(support vector machine, SVM)，RF，梯度回归树(boosted regression trees, BRT)，分类与回归树(classification and regression trees, CART)4个模型对土壤侵蚀速率进行分类预测，结果表明BRT和RF均具有出色的表现。另一方面，机器学习方法在主导因子分析研究中也有相对广泛的应用，如郝珊珊等^[10]使用BP神经网络对彭阳县土壤侵蚀主导因子的显著性进行比较分析，朱青等^[11]基于RF模型分析赣江上游流域不同子流域土壤侵蚀变化的主导因子重要性。

基于机器学习模型预测土壤侵蚀速率已有不少研究成果，但大都聚焦于单一算法模型，单一学习方法只能从单个角度观测特征数据，未能充分发挥机器学习算法的优势^[12]。为使不同学习方法能够取长补短训练出更优异的模型，进一步提升土壤侵蚀速率计算精度，笔者以三峡库区奉节县为例，综合对比当前常见的RF模型、极限梯度提升(eXtremegradient boosting, XGBoost)模型和K近邻回归(K-nearest neighbor regression, KNN)模型等9个机器学习模型，以此构建优化集成多学习器的Stacking模型，进而获取整个研究区内土壤侵蚀速率的空间分布图；最后在县级尺度下比较土壤侵蚀主导因子重要性，并绘制部分依赖图实现奉节县土壤侵蚀与其主导因子之间依赖性分析，分析结果可为三峡库区奉节县土壤侵蚀治理提供参考。

1 研究区概况

奉节县位于E 109° 1′ 17″~109° 45′ 58″、N 30° 29′ 19″~31° 22′ 33″之间，地处重庆市东北部，三峡库区腹地，全县境内处于三峡库区水土流失重点防治区。平均海拔约为946 m，其地形以山地为主，占总面积的88.3%。奉节气候总体上属于中亚热带湿润季风气候，四季较分明，雨量充沛，年均降水量达1 132 mm。奉节县境内地形复杂、海拔变化剧烈，当地土壤侵蚀敏感性较高，至2021年侵蚀等级为轻度及以上的面积占当地土地总面积43.28%，并主要发生在耕地、林地和建筑用地上。研究区概况如图 1所示。

2 数据与方法 2.1 数据获取与特征信息提取

本研究所使用的基本数据及来源如表 1所示, 各数据经ENVI、SNAP和Excel等软件预处理后，通过七参数转换与通用墨卡尔投影方法在ArcGIS 10.5软件中将地理坐标系统一为WGS84坐标系，投影系统一为WGS_1984_UTM_Zone_49N。

以往多数学者研究表明，土壤侵蚀速率主要受到地形、降雨、土壤性质、植被覆盖以及水土保持措施等多个方面的影响，其中地形是土壤侵蚀的决定性因素之一。陡峭的斜坡会加速土壤侵蚀，而较平缓的斜坡会抑制土壤侵蚀的加剧^[13]；降雨是土壤侵蚀最直接的因素，土壤侵蚀在空间上的分布很大程度取决于降雨侵蚀力的分布^[14]；土壤性质的影响是长期的。总体而言，土壤颗粒间的结合力足够强，土壤结构体就不易分散，土壤抗侵蚀能力就比较强；植被可以通过根部对土壤的固定能力以及冠层对降雨的缓冲作用来抑制土壤侵蚀^[15-16]；水土保持措施与人为活动密切相关，主要是通过对土地的管理来间接影响土壤侵蚀^[17-18]。这一方面体现在对土地的利用情况上，另一方面体现在人类活动运行轨迹上^[19]；遥感光学影像中包含着丰富且重要的地物光谱信息^[20]，对光谱信息进行波段运算等处理，可以获取一些遥感指数来凸显植被覆盖与土地利用情况的信息；因此，本研究将从地形、降雨、土壤性质、水土保持措施、光谱波段、植被指数6方面构建初始特征集。各特征基于ArcGIS空间插值、栅格计算器、欧氏距离以及矢量转栅格等工具获取，经裁剪后统一重采样为10 m分辨率的栅格数据。初始特征集如表 2所示。

2.2 特征优选

土壤侵蚀初始特征集中各特征间相关性大，而且存在很多冗余特征。因此笔者利用RF算法得到各个特征的相对重要性，基于特征类别进行特征优选，可过滤冗余特征，构造效率更高、消耗更低的独立的预测模型。45个初始特征的重要性排序如图 2所示。

研究根据各特征集中特征重要性排序以及特征间相关性构建机器学习优选特征集如表 3所示。

2.3 样本选取

样本数据集基于ArcGIS 10.5软件获取，整体研究区每隔1 km布设采样点，朱衣镇、永安街道、永乐镇等土壤侵蚀相对严重的乡镇区域按0.5 km为间隔布设采样点，最终得到抽样点共有6 962个，其中土壤侵蚀严重乡镇采样点为2 920个。通过多值提取至点方法将土壤侵蚀速率以及优选特征数据赋值到采样点中，按7 ∶3的比例划分训练集与验证集。

2.4 Stacking模型计算土壤侵蚀速率

Stacking算法是集成学习的新方向，它可以利用元学习器来综合集成不同学习器的结果，是1种有层次的集成学习。首先将训练集划分为n份，通过交叉验证的方式分别传入第1层的基学习器中进行训练，然后将预测结果拼接起来后作为次级特征集，结合训练集标签形成次级数据集，传入下一层的学习器中，最终经元学习器训练完成后得到最终结果。

Stacking算法中层数和数据集的折数可自由设置，但从各个领域的研究和应用来看，一般2层结构且选择5-fold进行交叉验证的Stacking算法既能强化学习效果又不至于造成模型过复杂^[21-22]，该算法原理如图 3所示。

将数据集载入不同基学习器模型中，经计算后观察其在验证集中的表现，将预测结果进行精度对比与相关性度量，构建基学习器的不同组合，并验证其在不同元学习器下的表现，从而得到Stacking模型最优的基学习器和元学习器组合。进而将带有坐标信息的制图点数据载入模型文件，经计算生成结果后产生栅格影像数据，最后对土壤侵蚀主导因子进行分析。

模型基于Python编程语言的中的Scikit-learn和Numpy等第3方库实现，采用网格搜寻法来确定每个模型的最佳超参数。通过使用决定系数(R-squared，R²), 平均绝对误差(mean absolute error, M_AE)和均方根误差(root mean square error, R_MSE)，将预测的土壤侵蚀速率与验证集上的土壤侵蚀速率进行比较，评估各机器学习算法的性能。

$ {R^2} = 1 - \frac{{\sum\limits_i^n {{{\left( {{y_i} - {{\tilde y}_l}} \right)}^2}} }}{{\sum\limits_i^n {{{\left( {{y_i} - \bar y} \right)}^2}} }}; $

(1)

$ {M_{{\rm{AE}}}} = \frac{1}{n}\sum\limits_i^n {\left| {{y_i} - {{\tilde y}_l}} \right|} ; $

(2)

$ {R_{{\rm{MSE}}}} = \sqrt {\frac{1}{n}\sum\limits_i^n {{{\left( {{y_i} - {{\tilde y}_l}} \right)}^2}} } 。$

(3)

式中：y_i为测试集上记录的土壤侵蚀速率；$ {{{\tilde y}_l}}$为模型中的预测值；y为记录的土壤侵蚀速率平均值；单位均为t/(km²·a)。

2.5 土壤侵蚀主导因子

虽然初始特征的重要性排序已经证实地形与降雨是对土壤侵蚀速率影响较大的变量，但光谱特征中存在一定相关性，难以分析植被覆盖因素以及土地利用状况等因素在其中的影响程度，研究分别从降雨、地形、土壤性质、水土保持措施和植被覆盖5方面选取13个可解释性与独立性较强的特征因子：AR_IDW、DEM、Slope、Aspect、OC、SAN、SIL、CLA、DFR、DFW、DFH、LU和NDVI，建立主导因子体系并选取排序前6的主导因子进行边际依赖性分析，量化土壤侵蚀主导因子影响程度，深入探寻土壤侵蚀速率分布规律。

3 结果与分析 3.1 基学习器对比和优化 3.1.1 基学习器精度对比

笔者分别采用集成算法中的RF模型、XGBoost模型、自适应增强(adaptive boosting, AdaBoost)模型、梯度提升(gradient boosting, GDBT)模型、轻型梯度提升器(light gradient boosting machine, LGBM)以及单一算法中的线性回归(linear regression, LR)、KNN、决策树(decision tree, DT)对土壤侵蚀速率进行训练与验证。

将每个模型应用于优化特征后的训练数据集，训练之后将模型在验证集上预测的土壤侵蚀速率与真值进行比较。模型的性能评估如表 4所示。结果表明，AdaBoost在所有集成学习器中表现最不好，结果存在一定的过拟合现象。而LGBM模型在验证集上产生了最有利的结果，在所有模型中R²最高，M_AE、R_MSE最低(R²: 0.843 7, M_AE: 265.87 t/(km²·a), R_MSE: 560.43 t/(km²·a)，RF、XGBoost模型紧随其后。另外，单一学习器在其中的表现极差，仅有DT表现尚可，LR、KNN和SVR模型R²均小于0.4。这表明其预测结果甚至不如平均值，因此它们不适用直接用于土壤侵蚀速率计算，提供的预测结果也不适用于作为次级数据集的新特征。

3.1.2 基学习器预测结果相关性度量

对于Stacking集成算法来说，不同基学习器的预测结果差异程度越大，Stacking集成后可以优化的程度就越深。因此在选出性能优异的基学习器后，还需考察各个基学习器预测结果的相关性，尽可能选择差异大的学习器。图 4是预测结果指标表现相对较好的6个学习器之间的相关性热力图。

以上几种集成学习器都是以DT作为弱学习器，数据观测方式存在较强相似性，所以相关性均高于0.88。其中RF使用了用于减少方差的并行集成方式，除RF外其他几类集成学习器使用了用于降低偏差的串行集成方式。同时从表 4中可以看到，LGBM与XGBoost模型预测结果高于其他几种串行集成学习器，故学习器中选择RF、XGBoost和LGBM进行组合作为基学习器。分别构建组合1：LGBM+RF+XGBoost；组合2：LGBM+RF；组合3：LGBM+XGBoost。在不同元学习器下进行验证对比。

3.2 元学习器对比和选取

将相同的训练集输入每个基模型组合中，根据图 3原理对样本数据进行预测，并将预测结果组合成元模型所使用的次级数据集的新特征信息。将输出的3组新特征分别送给9个元学习器中，表 5为各组基学习器在不同元学习器下的指标表现。

结果表明，当元学习器为LR和GDBT时精度有较明显的提高，在3种组合中，LGBM与RF提供的新特征信息有一定差异，因此提升幅度相对较大。综合来看，Stacking集成算法能够充分融合不同模型，最大程度发挥机器学习方法的性能。在本研究中，以LR为元学习器，RF与LGBM为基学习器的Stacking集成算法实现了最高的预测准确率，精度表现分别为R²=0.868 7；M_AE=252.48 t/(km²·a); R_MSE=537.78 t/(km²·a)。将预测结果与单体算法中表现良好的RF、XGBoost、LGBM进行对比，其中R²分别提升4.7%、5.4%和2.5%；M_AE分别降低21.6、52.7和13.4 t/(km²·a)；R_MSE分别降低59.42、105.21和22.45 t/(km²·a)。图 5为本研究所构建的优选集成Stcaking模型与RF、XGBoost和LGBM模型的验证值与预测值对比图。

3.3 土壤侵蚀空间分布结果与精度

基于上文构建的Stacking模型进行预测整个奉节县土壤侵蚀速率，并形成空间分布数据，为直观对比模型结果与真实数据之间差异，参照SL 190—2007《土壤侵蚀分级分类标准》将土壤侵蚀速率分为: 微度(0~500)、轻度(＞500~2 500)、中度(＞2 500~5 000)、强烈(＞5 000~8 000)、极强烈(＞8 000~15 000)、剧烈(>15 000)6个强度等级，其中单位为t/(km²·a)。结果对比如图 6和表 6所示。结果表明，基于本研究构建Stacking模型的土壤侵蚀等级预测结果能实现较高程度的还原，经计算，Kappa系数为0.713，符合基本一致标准，总体精度为82.13%，主要在“极强烈”与“剧烈”等级的预测精度相对较低。分析认为，首先是由于较高土壤侵蚀等级的分布过少，使得原始数据集较高速率的土壤侵蚀样本过少，在训练过程中得不到充分训练；其次在土地利用情况分别为“阔叶林”“高覆盖度草地”“乔灌果园”的区域误差较高。这3者在影像中极为接近，而侵蚀速率差异较大；最后在奉节县部分土壤侵蚀严重区域存在水土保持林、水源涵养林等人工实生的水土保持工程，这些使训练产生一定的偏差。

3.4 主导因子

主导因子重要性排序如图 7所示。排序前6的因子分别为DEM、AR_IDW、NDVI、Slope、DFR和DFW。因子重要性均高于9 %，这证明土壤侵蚀是多方面因素共同作用的结果。与其他研究不同的是，由于奉节县88%的地形为山地，其坡度变化频繁，受当地多年人为活动如农业种植、生态工程等影响，坡度对土壤侵蚀的影响远低于高程带来的影响。土地利用类型中林地占研究区总面积77%，且缺少对灌木林、阔叶林等做进一步细分，使该因素重要性排序较低，但通过距城区道路等可达性因子的排序可看出，人为活动依然是土壤侵蚀中较为重要的因素。

部分依赖图(partial dependence plots)适用于可视化反映土壤侵蚀速率对某个主导因子的依赖程度以及之间是否存在线性、单调或更复杂的关系^[9]，当指定因子在其边际分布上发生变化时，部分依赖图会显示平均预测值的变化。借助部分依赖图，可以量化土壤侵蚀主导因子的影响。

部分依赖关系如图 8所示，侵蚀速率与主要主导因子之间均为非线性关系，随着高程从200~520 m，土壤侵蚀速率从3 186降低至2 946 t/(km²·a)，而后降低趋势开始加剧，在800 m后降低趋势趋于平缓，在DEM大于1 232 m的区域土壤侵蚀速率开始增高。其主要原因在于该区域内石灰岩与紫色土分布广泛，容易产生土壤侵蚀^[23]，而在高程大于1 535 m时，由于人类活动对土地扰动小，植被覆盖度高，使得其侵蚀现象又有所减少。降雨量对土壤侵蚀影响为正向的，随降雨量增多，土壤侵蚀速率不断加快，在降雨量高于1 250 mm时达到最大水平。土壤侵蚀速率与植被存在负向关系，且在0.24~0.27之间最高，整体上随着NDVI的增高，侵蚀速率不断增长。由于奉节县果园产业尤其是脐橙果园，在影像中难以与林地区分开来，另外对于部分土壤侵蚀严重区域，开展种植水土保持林、人工灌木林等水土保持工程使部分依赖图中出现一些反复。坡度与侵蚀速率之间整体存在一定正相关，随着坡度从15°到32°，土壤侵蚀速率提升至4 210 t/(km²·a)，并在26°~35°之间维持较高水平。距道路距离在220 m以内时侵蚀速率较高，并随距离增加而减小，在1 396 m处侵蚀速率最低值，并在此后趋于稳定。距水源距离在小于63 m处土壤侵蚀速率极低。此处主要受三峡大坝蓄水影响导致水面涨落影响，在距离106~387 m之间维持较高的水平，而后从距水源387~625 m处整体呈降低趋势，之后出现一定上升趋势。该上升区域主要为果园种植区，其与土壤侵蚀区域存在一定程度上的重叠。

4 结论

1) 在土壤侵蚀速率计算中引入Stacking算法，能够有效提升预测精度。本研究中，以LR为元学习器、RF与LGBM为基学习器的Stacking集成算法具有更好的性能表现，相较于RF和LGBM算法，R²分别提升4.7%和2.5%；M_AE分别降低8.6%和5.3%；R_MSE分别降低11.0%和4.2%。该算法顾及各子算法的多样性和差异性，并保留其优势，使预测精度更高、更稳定。

2) 奉节县土壤侵蚀是多方面因素综合作用的结果，其中DEM、AR_IDW、NDVI、Slope、DFR和DFW等6项主导因素发挥关键作用。总体而言，奉节县土壤侵蚀速率与高程、植被覆盖程度之间呈正相关关系，与降雨量、坡度之间呈负相关关系。速率较高的土壤侵蚀倾向于发生在降雨充沛、植被覆盖度低、距道路及水源较近的低海拔陡峭区域。

5 讨论

本研究以三峡库区奉节县为研究区，选择9个机器学习模型对土壤侵蚀速率进行预测，对比分析各模型的预测精度与相关性，并以此作为基学习器构建Stacking模型。结果表明，Stacking模型能充分利用不同算法从数据进行观测，使不同算法能够取长补短，获得最优预测结果。其中基学习器的选取应当遵循优秀且异质原则，避免输入冗余或无效的新特征；而元学习器的选择则是以简单为主，从而避免模型产生过拟合。构建包含降雨、地形、土壤性质、植被覆盖和水土保持措施5方面的主导因子体系，对研究区(奉节县)土壤侵蚀速率分布规律进行分析，可以为相关部门进行土壤侵蚀治理提供参考。但研究还有进一步提升的空间：土地利用类型数据中未对针叶林、阔叶林、果园和旱地等做进一步细分，因此未来研究中可以尝试添加其他合适的特征变量，在土壤侵蚀速率较高的小尺度区域内进行实验，进一步提升预测性能。