基于多源辅助变量和随机森林模型的表层土壤全氮分布预测*

周 洋，赵小敏†，郭 熙

（1. 江西农业大学国土资源与环境学院，南昌 330045；2. 江西省鄱阳湖流域农业资源与生态重点实验室，南昌 330045）

作为植物生长期间所必需的养分元素，土壤氮素是决定土壤质量和肥力的主要因素和指标，但是，土壤中的氮素如果过量则会导致水体富营养化等农业面源污染问题，且氮素能形成温室气体进入大气。因此，出于精准农业和环境质量管理的需要，了解和掌握土壤氮素含量的空间分布显得尤为重要。

土壤全氮（Soil total nitrogen，STN）是指土壤中各种形式氮素的总和，由于受到人为活动、土壤成土因素和地形因子等因素的影响，STN 往往表现出较强的变异程度，准确估计STN 的空间分布及其受到的主控因子可以为精准施肥等农业管理措施的实施提供理论依据和指导。随着“3S 技术”的发展，地形、遥感、气候等环境因子变得易于获取，因此有学者提出了数字土壤制图方法来实现土壤属性在区域范围内的精准预测。在数字土壤制图中，多元线性回归、普通克里格法、回归克里格法、地理加权回归等方法被广泛用于预测STN 等土壤属性的空间分布，但土壤与环境变量往往表现为复杂的非线性关系，多元线性回归和基于地统计的克里格方法不能很好地体现这种复杂关系，因此有学者将人工神经网络、回归树和随机森林等机器学习模型应用在土壤属性预测中，其中随机森林模型（RF）可以防止过度学习和过拟合，相比人工神经网络计算更加简便，且预测效果稳定，被证明是一种有效的土壤属性预测方法，还有研究指出，将RF 模型的预测残差与预测结果相结合的随机森林残差克里格法（RFRK）能提高随机森林模型的预测精度。Guo 等对海南岛橡胶人工林的土壤有机质进行了预测，研究表明，相比逐步线性回归（SLR）和RF，RFRK 表现出了更好的预测效果。Hengl 等对比了普通克里格（OK）、多元线性回归（MLR）、RF、梯度提升树（GBT）和3 种混合模型对希腊北部地区土壤有机质的预测效果，结果表明RFRK 的预测精度要优于其他方法。Scarpone 等使用了广义线性模型（GLM）、RF 及其与地统计学相结合的混合方法对加拿大哥伦比亚省中南部的土壤厚度进行了建模预测，发现混合方法的预测精度要优于普通的机器学习算法。作为一种有效的土壤属性预测方法，RFRK 得到了越来越广泛的应用，但在国内的土壤属性预测研究中较少采纳。寻乌县位于赣南地区，地形以山地丘陵为主，属于典型的南方丘陵区，果园产业发达，而STN 对土壤肥力和土壤碳循环均有着很大的影响，因此，很有必要对寻乌县STN 的空间分布进行分析和预测，为此本研究基于地形因子、遥感因子等多源辅助变量，使用了随机森林（RF）和随机森林残差克里格（RFRK）两种预测模型对寻乌县表层STN 的空间分布进行了预测，为寻乌县的土壤肥力评估和精准施肥等农业措施提供理论指导和数据支撑。

1 材料与方法

1.1 研究区概况

寻乌县位于江西省南部，为赣州市所辖县，位于闽粤赣三省相连处，地理坐标介于24°30′40″—25°12′10″N 和115°21′22″—115°54′25″E 之间，面积为2 352 km。该地区属于亚热带季风气候，阳光充足，雨量充沛，年平均气温为18.9 ℃，在7 月达到最大值27.2 ℃，1 月最低，为9.1 ℃，而年降水量则介于1 651 mm 到1 750 mm 之间，其中4—6 月降水量最多。寻乌东北、西北与东南部地势高，且向西南部倾斜，海拔高度介于132～1 468 m 之间，地貌类型以山地丘陵为主（图1a）。县内植被覆盖率高，主要植被类型为亚热带常绿阔叶林和低山针叶林。寻乌属于南方亚热带红壤区的南部，土壤肥力较高，pH 普遍呈酸性，土壤类型主要为红壤、黄壤和水稻土；母质类型主要有酸性结晶岩类风化物和泥质岩类风化物，还有少量的石英岩类、泥质岩类和红砂岩类风化物；土地利用类型则以林地、园地和耕地为主，其中林地占土地总面积的83.8%（图1b）。

图1 研究区DEM、土壤样点分布及土地利用类型Fig.1 Digital Elevation Model，distribution of soil sampling points and land use types in the study area

1.2 样点数据获取与测定

采用公里网格布设土壤采样点，于2018 年7 月至9 月完成采集，采样深度为0～20 cm，共采集土壤样品572 个（图1a）。土壤样品经风干，研磨和过筛后测定全氮（TN）、硫（S）、pH 等元素和指标含量。其中TN 采用半微量凯氏定氮法测定；S 采用燃烧碘量法测定；pH 采用电极法测定。

1.3 环境协变量及其数据来源

初步选择了22 种环境协变量对研究区STN 进行预测，包括地理坐标（、坐标）、地形因子、距离因子、气候因子、遥感因子和土壤理化因子。利用ArcGIS10.2 软件平台，将所有的协变量栅格数据重采样为100 m 分辨率。地形因子是土壤属性预测中使用最广泛的环境因子之一。本研究使用数字高程模型（Digital elevation model，DEM）来得到各种地形因子。从地理空间数据云平台下载研究区30 m 分辨率DEM 数据，使用ArcGIS 10.2 软件平台处理后得到高程、坡度、坡向、曲率、地形起伏度5 个地形因子，使用SAGA GIS 2.3.2 软件平台处理后得到地表粗糙度TRI、地形湿度指数TWI 和多尺度谷底平坦度MrVBF 这3 个地形因子。距离因子包括距河流距离（DTR）和距公路距离（DTH），利用ArcGIS 10.2 软件平台的距离分析工具计算得到。气候因子使用从中国科学院资源环境科学数据中心下载的1 km×1 km 分辨率的年平均气温（MAT）和年平均降水量（MAP）栅格数据集（2006 年至2015年均值）。遥感因子从地理空间数据云平台上下载的2017 年11 月1 日的寻乌县Landsat 8 OLI 遥感影像中提取，利用ENVI 5.3 软件平台进行辐射定标和大气校正后，提取归一化植被指数（NDVI）、绿带反射率（B3）、裸土指数（BI）、绿色叶绿素指数（CIg）、绿度指数（GI）和对波段进行主成分分析后的第一主成分（PCA1）。其中NDVI 和绿度指数是反应植被生长和覆盖情况的常用指标，裸土指数体现了地表的裸露情况，而绿色叶绿素指数则反映了植物冠层氮的情况。土壤理化因子选择了pH 和硫，采用普通克里格法进行空间插值后，提取100 m 分辨率的栅格数据。为了避免预测变量之间的多重共线性影响模型预测效果，在将协变量代入模型之前，计算了各环境变量的方差膨胀因子VIF，将VIF ＞ 10的变量剔除后再代入模型进行计算。

1.4 建模方法

随机森林（RF）是一种基于多个决策树的数据挖掘算法，通过在训练数据集中应用bootstrap 方法随机抽取样本构建单独的树模型，生成大量的树后组成随机森林，与CART 决策树模型相比，RF 模型对过度拟合不敏感。决策树的数量ntree 和每个分隔节点的预测变量数mtry 是RF 模型中的两个关键参数。本研究在随机划分100 组训练集和验证集数据后，使用R3.6.1 中的“caret”包对每个RF 模型分别进行五折交叉验证调参，之后使用“randomForest”包，代入不同训练集对应的最优参数来构建随机森林模型。随机森林残差克里格（RFRK）是RF 模型的一种拓展，利用普通克里格法，对RF 模型的预测残差进行空间插值后将其与RF 模型的预测结果相加，从而提高模型预测精度。RFRK 模型的公式如下：

式中，（）为RFRK 在x位置的STN 预测值，（）为RF 模型拟合的趋势项，则为普通克里格法拟合的误差项。

1.5 模型精度与预测不确定性评估

为了获取稳定的模型预测值，对RF 模型重复迭代了100 次，每次均随机选择70%的土壤样点（400 个）作为训练集参与建模，其余30%的样点（172个）作为验证集，最后取100 次迭代结果的平均值作为 RF 模型的预测结果。选取了平均相对误差（MAE）、均方根误差（RMSE）、决定系数（）和一致性相关系数（CCC）4 个模型精度指标评价模型预测结果，计算公式如下：

2 结 果

2.1 土壤全氮描述性统计特征

寻乌县STN 的描述性统计结果表明，土壤TN含量范围为0.42～2.52 g·kg，均值为1.21 g·kg，低于全国的平均水平；标准差为0.35 g·kg，偏度系数为 1.10，峰度系数为 3.39，变异系数为26.88%，属于中等程度的空间变异。

不同土壤类型、土地利用类型和成土母质下样点STN 含量的分布情况如图2 所示。不同土壤类型下STN 的分布差异不大，黄壤的整体STN 含量分布要高于红壤和水稻土，不同土地利用类型STN 含量中位数大小依次为旱地＞水田＞林地＞园地，5 种成土母质的STN 分布有较明显的差异，紫色岩类风化物母质的STN 含量要低于其他4 种母质。

图2 不同土壤类型、土地利用类型和成土母质下土壤全氮含量小提琴图Fig. 2 Violin plot of soil TN content in different soil types，land use types and parent materials

2.2 环境协变量的选择

为避免变量共线性对模型预测精度的影响，计算了所有环境协变量的方差膨胀因子VIF，并去除了VIF＞10 的变量（坐标、地形起伏度、年平均气温、归一化植被指数、绿带反射率和裸土指数），之后计算了STN 与剩余变量之间的Pearson 相关系数。结果表明，STN 与土壤S（=0.72）、高程（=0.32）、坡度（=0.19）、地形粗糙度（=0.17）绿色叶绿素指数（=0.12）等变量表现出正相关关系，与坐标（= -0.17）、年平均降水（= -0.12）等因子则有着负相关关系。

2.3 不同模型预测精度及不确定性评价

RFRK 需要对RF 模型的预测残差进行普通克里格插值，因此，对STN 的RF 模型100 次预测结果的残差分别进行半方差分析，得到理论最优半方差函数模型及参数后，再使用普通克里格法对其进行空间插值，得到RF 模型100 次预测结果残差的空间分布。

在每次随机抽取30%的样点数据作为验证集，并进行100 次模型迭代后，对RF 模型和RFRK 模型的预测精度进行了评估，结果表明RF 模型的预测精度指标的均值（MAE=0.1570 g·kg，RMSE=0.210 8 g·kg，CCC=0.761 3，=0.629 1）均优于RFRK 模 型（ MAE=0.168 2 g·kg， RMSE=0.226 7 g·kg，CCC=0.688 1，=0.571 9），这表明将预测残差作为误差项加入RF 模型并没有改善其预测效果，反而使RF 模型的预测精度下降。

使用RF 模型迭代100 次的预测值标准差来评估模型预测的不确定性。RF 模型预测值标准差的范围为0.01～0.22 g·kg，均值为0.04 g·kg，总体上模型的预测不确定性较低，仅在研究区的西北部、西南部、东部以及南部等海拔较大的区域表现出较大的预测不确定性，这可能是因为环境协变量在高海拔区域有着较强的空间变异性。

2.4 环境协变量重要性评估

对RF 模型中的环境协变量计算了其对STN 预测结果的相对重要性（Relative importance，RI），并将其转化为百分比（图3）。结果表明，土壤硫，高程，年平均降水，绿色叶绿素指数，坐标，多尺度谷底平坦度、坡度、地形粗糙度这些因子是影响研究区STN 含量的主要协变量，占RI 的87.44%；此外，在所有类别的协变量中，土壤理化因子的解释能力最高（占RI 的48.52%），其次为地形因子（24.06%）、遥感因子（9.76%）、气候因子（6.88%）和地理坐标（5.45%）。

图3 RF 模型各变量相对重要性Fig. 3 Relative importance of variables in the random forest（RF）model

2.5 土壤全氮对各环境协变量的响应变化

为探究STN 与环境协变量之间的响应变化关系，随机选择了100 次迭代模型其中的1 个RF 模型，绘制了各环境协变量对STN 的响应变化曲线（图4），该曲线代表了STN 与每个环境变量之间的变化关系，可以看出STN 与各变量的响应关系差异较大，且大致趋势与相关性分析结果一致，其中，STN 含量随土壤硫、高程、绿色叶绿素指数、坡度、地形粗糙度、距河流距离、pH 和波段第一主成分等因子值的增大而表现出非线性增加，与年平均降水则表现为相反趋势的非线性变化，而与y 坐标、多尺度谷底平坦度、距公路距离、绿度指数等因子表现为波动变化趋势。在影响STN 含量主要协变量的响应变化曲线中，高程、坡度、地形粗糙度、多尺度谷底平坦度和pH 增大到一定值后，STN 含量不再响应其变化。

图4 土壤全氮含量与各环境协变量的响应曲线Fig. 4 Response of soil total nitrogen content to the variation of each environmental predictive factor

2.6 土壤全氮空间分布特征

研究区 STN 的空间分布表现出较强的空间异质性（图5），RF 模型和RFRK 模型两种方法预测的STN 空间分布的总体变化趋势基本一致，均表现为西北部、东北部和南部高，中部和中北部低的空间格局，且与高程的变化特征大致相符；RF 模型和RFRK 模型迭代100 次的STN 预测值均值的变化范围分别为0.70～2.05 g·kg和0.69～2.10 g·kg，均值分别为1.23 g·kg和1.24 g·kg，从结果中可以看出，RF 模型与RFRK 模型的STN 预测值范围基本一致，与土壤样点数据的统计范围较为接近。

图5 不同模型下研究区土壤全氮空间分布图Fig. 5 Spatial distribution of soil total nitrogen（STN）by different model

3 讨 论

3.1 影响土壤全氮分布的主要变量

在所有的环境协变量中，土壤硫对研究区STN含量分布的影响最大，这是因为土壤S 元素和N 元素的循环有着相互耦合性，从而直接影响STN 的分布。从响应曲线（图5）中可以看出，土壤硫与STN 表现出了正向耦合性。在施氮肥时，含硫酸根离子的氮肥会使土壤中的S 含量增多，因此土壤硫和STN 的空间分布呈现出相互影响的变化趋势。高程、坡度等地形因子对STN 的影响仅次于土壤理化因子，两者与STN 含量的响应曲线变化的大致趋势也较为相似。有研究表明，地形因子在土壤属性预测中有着重要的作用，其对STN 的影响主要体现在影响地表水热交换和影响植被分布这两个方面。不同地形条件下，表层土壤受到的风化淋溶作用存在差异，从而影响STN 的含量。而植被覆盖度较高处，表层土壤有机质含量较高，STN 等土壤养分元素易积累，且地形因子中的高程变化对区域小气候特征也会造成一定的影响。随着海拔的上升，气温会呈现出下降的趋势，而STN 在不同的土壤温度下，分解速度也会存在差异，导致不同海拔条件下STN 含量存在差异。遥感因子对研究区STN 也有较大影响，这与其他研究一致。有研究表明，在植被覆盖度较高的区域，植被与表层土壤中的氮含量较为相关。因此，在植被茂密地区，波段反射率和遥感指数可以作为预测STN 的有效指标。本研究中，绿色叶绿素指数和绿度指数等遥感因子对STN 的影响较大，但与其他研究相比，影响程度较低，这可能与本研究选取的遥感指数有关。由于共线性问题，未将归一化植被指数和绿带反射率等遥感指数纳入到预测模型中；气候因子中，年平均降水对STN 的影响主要体现为，降水量的不同会导致地区的水文环境和土壤水分产生差异，影响土壤氮素的矿化速率，从而影响STN的分布。地理坐标对STN 也有一定的影响，表明研究区 STN 的空间分布有着较强的空间异质性。Pearson 相关性分析结果表明，STN 与y 坐标有着负相关关系，响应曲线中，STN 随着y 坐标的增大表现出先增后减的变化趋势，这与研究区特点有关。寻乌县属山地丘陵地形，三面环山，东北、西北和东南均有高山分布，区域的高程差近1 000 m，因此研究区地理坐标的改变会使成土母质和植被产生较大的垂直性地带变化，进而对STN 产生较大的影响。

3.2 不同模型预测效果的对比

本研究中，RF 模型100 次重复迭代的模型精度指标值（MAE、RMSE、CCC、）要优于RFRK模型，这与其他研究结果有所不同，原因在于，虽然残差代表了RF 模型在预测过程中无法解释的部分，但因为本研究中选取的环境协变量足够多，RF 模型的预测效果已较好，预测残差的空间自相关性较弱，使用普通克里格法对残差进行空间插值后并不能代表RF 模型预测误差的空间分布趋势，因此将残差的空间信息加入到RF 模型的预测结果中并没有提高其预测精度，反而使RF 模型的预测精度有所下降。

3.3 预测精度的影响因素

在对STN 进行建模预测时，影响预测精度的因素主要有两点。第一，部分协变量数据的精细化程度不高，在使用的环境协变量中，年平均降水因子使用的是1 km×1 km 分辨率的栅格数据，而研究区为县域尺度，该分辨率可能忽略了气温和降水的局部空间自相关性，从而降低了STN 的预测精度；土壤理化因子中，土壤硫和pH 的栅格数据是利用点位数据，采用普通克里格法进行空间插值得到，而普通克里格法插值中的误差和不确定性也会对STN的预测效果产生影响，且为了统一环境协变量的尺度而进行的栅格重采样也会导致协变量数据精度降低和信息丢失的问题。第二，代表人为活动因素的协变量考虑的不够多，所选取的环境协变量中，除了距公路距离这个人为因素之外，其他的多为自然因素，虽然遥感植被因子能够在一定程度上代表人为活动，但如果将更多的人为活动因素，比如耕作制度、施肥方式等因素作为协变量纳入预测模型，可能会在一定程度上提高模型预测精度和效果。

4 结 论

RF 和RFRK 模型两种方法对寻乌县表层土壤的STN 预测结果表明：从空间分布来看，两种方法预测的STN 空间分布大致相同，总体上表现为西北部、东北部和南部高，中部和中北部低的空间分布格局；从预测精度来看，RF 模型的决定系数均值和一致性相关系数均值均高于RFRK 模型，而平均相对误差均值和均方根误差均值均小于RFRK 模型，表明RF 模型的预测结果要优于RFRK 模型；在参与预测的环境协变量中，土壤硫，高程，年平均降水，绿色叶绿素指数，坐标，多尺度谷底平坦度、坡度、地形粗糙度等因子是影响研究区STN 分布的主要变量。