优越的土壤肥力是保证粮食作物产量的基础，在现代农业生产中有着决定性作用^[1-2]。而土壤速效钾是土壤肥力中重要的影响因素，是作物生长发育必不可少的养分元素之一，对作物产量有着重要影响^[3-4]。耕地土壤速效钾空间分布主要受空间位置、地形部位、成土母质、气候环境、人为施肥等多种因素影响^[5-7]。相对而言，安庆市耕地土壤速效钾含量整体偏低^[8]，因此准确预测其含量及分布情况对指导农业施肥、促进作物增产、提升耕地质量等级尤为重要。

传统数字土壤制图的空间预测方法主要包括地统计插值法和确定性插值法^[9]。地统计插值法的代表性方法是普通克里金插值法^[10]，确定性插值法的代表性方法是反距离权重插值法^[11]。近年来，机器学习方法，尤其是随机森林模型逐步成为数字土壤制图方向的国内外研究热点。Kalambukattu等^[12]在遥感和地形数据的基础上，利用人工神经网络模型绘制了喜马拉雅流域土壤有机碳和氮元素空间分布图；高凤杰等^[13]运用贝叶斯最大熵模型预测了黑土区小流域土壤有机质空间分布，通过与协同克吕格法比较，结果表明贝叶斯最大熵模型显著优于协同克吕格法。然而人工神经网络模型需要大量参数，学习时间过长，不能观察之间的学习过程；贝叶斯模型对数据的表达形式过于敏感；决策树^[14]则有过拟合的缺点。不同学者各自通过随机森林模型对不同研究区的土壤养分含量进行预测，结果均表明随机森林模型的预测精度十分可靠，是数字土壤制图的理想方式^[15-16]。相对于其他机器学习方法，随机森林可以输入大量变量，且学习过程较快，同时当部分数据缺失时，仍然可以维持相当高的精准度。因此，本研究采用了随机森林来训练模型，探索随机森林在数字土壤制图中的应用效果，同时与传统插值模型中最为常见的普通克里金和反距离权重方法作比较，旨在准确预测安庆市耕地土壤速效钾空间分布模式。

1 材料与方法

1.1 研究区概况

安庆，又称宜城，辖二市五县三区，位于皖西南部，是皖西南区域中心城市。安庆市属于亚热带季风性湿润气候，年平均降水量为1368 mm，年平均温度为16.5℃。地势整体西北高，东南低，呈山地——丘陵岗地——平原阶梯状分布。土壤母质类型丰富，主要成土母质为酸性结晶岩类风化物、山河冲积物、泥页岩类风化物；土壤类型多样，以水稻土、黄棕壤、黄褐土、潮土为主。

1.2 采样点数据来源

本研究采样点数据来源于2019 年安庆市耕地质量等级评价取样结果，布点要求平原地区耕地每667 hm² 不少于一个点位，山地地区耕地每133 ～ 150 hm² 一个点位，并尽可能分布至全市各乡镇、土类和成土母质。采样按照五点取样法，采取0 ～ 20 cm耕地表层土壤化验，全市共获取682 个采样点位，随机选择15%（103 个点位）的样点作为验证样点，85%（579 个点位）的样点作为测试点位。具体分布情况如图1 所示。

1.3 环境变量因子来源与选取

研究区地形复杂多样，地势起伏巨大，海拔落差峰值达到1700 余m，地貌类型主要包括西北大别山区、中部丘陵岗地区、东南长江沿江平原区；成土母质类型丰富，土壤类型多样；气候为亚热带季风性湿润气候，气候温和、雨量充沛，但因地形地貌不同，研究区各区域气温和降水量差异较大，直接影响岩石风化、耕地土壤矿物质的移动和分布，进而直接或间接影响研究区耕地土壤速效钾的分布情况。本研究选取的环境变量因子主要包括空间位置变量因子、地形变量因子、土壤变量因子和气候变量因子。其中，空间位置变量因子包括经度和纬度；地形变量因子包括高程（图2）、坡度（Slope）、坡向（Aspect）、剖面曲率（Profile curvature）、平面曲率（Plane curvature）和地形起伏度，数据来源于地理国监测云平台（http：//www.dsac.cn）获取的30 m分辨率研究区数字高程模型（Digital elevation model），坡度、坡向、剖面曲率、平面曲率，地形起伏度数据运用了ArcGIS 10.2 的空间分析工具（Spatial Analyst），根据数字高程模型计算得来；土壤变量因子包括成土母质和土壤类型，成土母质（图3）与土壤类型图，数据来源于安庆市第二次全国土壤普查纸质图件，经ArcGIS 10.2 矢量化形成；气候变量因子包括年平均降水量（图4）和年平均温度（图5），数据来源于全球气候数据网（http：//www.worldclim.org），获取了分辨率为1000 m的1970 ～ 2000 年的年平均降水量及年平均温度数据。

1.4 空间预测方法

1.4.1 随机森林

随机森林（Random Forest，RF）是机器学习中常用的一种方法，最早由Breiman于20 世纪80 年代提出^[17]。模型运用bootstrap随机重复选取样得到预测集，并对每个预测集组合决策树，得到决策树数目（ntree）。在各个决策树节点，选择样本预测器的数目（mtry）作为分割变量，输入矩阵值，通过汇总分类树的结果组合成随机森林，并以此预测变量。公式如下：

式中，$\overline{RF}\left(x\right)$表示预测结果值；ntree 表示决策树数目；$X_{mtry}$表示输入矩阵；$RF\left(X_{mtry};{\theta }_i\right)$表示第i棵决策树的预测结果。

本研究采用MATLAB2015 R2015b软件编程，选取地形变量因子（高程、坡度、坡向、剖面曲率、平面曲率和地形起伏度）、土壤变量因子（成土母质、土壤类型）、气候变量因子（年平均降水量、年平均温度）以及空间位置变量因子（经度、纬度），共计12 个变量因子作为预测变量。一般来说，mtry取值为变量因子总数三分之一^[18]，ntree一般选取默认值为500^[19]，但当ntree的增加不能显著提升预测精度时，ntree需要尽可能小^[20]。本研究通过3 组12 次试验（ntree=200，300，400， 500；mtry=2，3，4），对比各试验组中验证样点的R² 值。其中，R² 越接近1，拟合精度越高。结果表明，在试验组中，当ntree=400，mtry=3 时，研究区RF预测精度最高。

1.4.2 普通克里金

克里金插值法（Kriging）起源悠久，南非矿产工程师Danie G.Krige在20 世纪50 年代初开创性地采用此方法找寻金矿^[21]，法国统计学家Georges Matheron为此方法奠定了基础，并命名为Kriging^[22]。普通克里金（Ordinary Kriging，OK）插值法作为地统计学的重要内容之一，在空间自相关基础上，利用半变异函数对区域内样点进行预测分析，OK要求原始数据或经过函数变换后的数据正态分布。其公式为：

式中，s 表示不同空间位置的坐标点位，可默认为经纬度坐标；$Z\left(s\right)$是 s 处的变化量；$\mu$表示趋势值；$ϵ\left(s\right)$表示自相关误差。

1.4.3 反距离权重

反距离权重（Inverse Distance Weighting，IDW） 插值法是基于地理学第一定律的原理：即万事万物都是有联系的，两个物体距离越近，其性质越接近；两个物体距离越远，其性质差异性越明显^[23]。IDW插值法作为运用最为广泛和普遍的空间插值方法，通过已知点和预测点之间的距离作为权重并加权平均计算，离已知点越近的预测点权重越大。 IDW插值法要求插值对象分布均匀适中，并且插值面有局部因变量。其公式为：

式中，$Z\left(s_0\right)$为 s₀ 处的预测值，N 为预测点周围样点需要使用的数目；${\lambda }_i$为使用样点的权重值，$Z\left(s_i\right)$是 s_i 处获得的已知值。

1.4.4 研究区模型精度检验

研究区样点由85%预测样点（579 个）和15%验证样点（103 个）组成，模型精度检验是对研究区103 个验证样点的误差估计，通过计算验证样点的平均绝对误差（MAE）、均方根误差（RMSE） 和决定系数（R²）确定各模型预测效果，其中， MAE、RMSE越小，R² 越接近1，模型预测精度越高，预测效果越好。具体计算公式如下：

式中，$m_i$、$m_j$、$\stackrel{-}{m}$、$n$分别表示预测样点值、验证样点值、验证样点平均值和参与精度检验的样点总数。

2 结果与分析

2.1 耕地土壤速效钾描述性统计

研究区采样点随机划分为预测样点组和验证样点组。从表1 可知，预测样点速效钾值预测范围为4 ～ 691 mg·kg^-1，验证样点速效钾值范围是23 ～ 349 mg·kg^-1，平均值、中值和标准差预测样点均略大于验证样点，预测样点存在少量极端极大值。变异系数分别为74.1%和69.8%，均属于中等变异水平。预测样点偏度和峰度分别为2.58 和11.76，经过Box-Cox变换，同时经过K-S检验后符合正态分布，变换参数为0.1。

2.2 研究区模型精度评价

通过对比研究区103 个验证样点与3 种模型预测结果，MAE、RMSE和R²，确定各模型精度。其中，RF的MAE、RMSE、R² 分别为30.93 mg·kg^-1、 41.31 mg·kg^-1 和0.58，OK的MAE、RMSE、R² 分别为31.97 mg·kg^-1、44.08 mg·kg^-1 和0.49，IDW的MAE、RMSE、R² 分别为32.77 mg·kg^-1、46.21 mg·kg^-1 和0.47。从表3 可以看出，在3 种方法中，RF的MAE和RMSE均为最小，R² 最大，RF的MAE较其余两种方法分别减少了1.04 和1.84 mg·kg^-1，RMSE较其余两种方法分别减少了2.77 和4.9 mg·kg^-1，R² 较其余两种方法分别增加了0.09、 0.11；OK的MAE和RMSE、R² 均为第二，OK的MAE较IDW减少了0.8 mg·kg^-1，RMSE较IDW减少了2.13 mg·kg^-1，R² 较IDW增加了0.02；IDW的MAE和RMSE均为最大，R² 最小。可以得出结论：研究区内速效钾空间分布的3 种方法的预测精度高低顺序分别为RF>OK>IDW。

通过筛选速效钾极大值（>200 mg·kg^-1），发现OK的预测结果中的极大值最低，由表4 可得其R² 为0.13，MAE为93.92 mg·kg^-1，RMSE为105.13 mg·kg^-1，IDW虽然因“牛眼”效应（某些极值点在插值结果上所形成的次插值点为圆心的圆斑，且圆斑的值与该极值点较为接近）使得极大值较大，但预测精度在3 种方法中仍然最低，其R² 仅为0.0002，MAE为107.64 mg·kg^-1，RMSE为111.37 mg·kg^-1。RF的极大值较大，其R² 为0.30，MAE为89.90 mg·kg^-1，RMSE为101.90 mg·kg^-1，可以看出，RF对极大值的预测精度较另外2 种有较大提高，但极大值部分的预测精度较整体的预测精度低。

2.3 随机森林环境变量因子权重

为确定各个环境因子在RF模型中的权重，通过移除某个环境变量因子，核查RF模型RMSE增加的比例，从而判定各个环境因子的权重及其排序。从表5 中可以看出，纬度、年平均温度、成土母质、高程、经度、年平均降水量是影响研究区RF模型下速效钾含量预测的重要因素，平面曲率、剖面曲率、地形起伏度、坡度、坡向与土壤类型权重较低。前5 个地形变量因子与高程有相关性，土壤类型与成土母质有相关性，可以得知，此类协同变量因子在RF模型中的权重会因共线性而导致其权重降低，这与任丽等^[24]的研究结果一致。经度、纬度通过空间位置的变化对速效钾含量产生影响。成土母质是土壤形成的基础，是速效钾含量的初始影响因素。年平均温度和年平均降水量则通过风化作用，崩解、迁移成土母质中的化学物质，从而影响速效钾含量。

2.4 耕地土壤速效钾空间分布预测结果

对比3 种方法下研究区耕地土壤速效钾空间分布预测结果（图6），可以看到整体分布趋势接近，研究区耕地土壤速效钾含量在东南沿江处较高，同时研究区中南部出现小块高值区，其原因在于两部分区域地貌类型为平原，高程较低，成土母质以长江冲积物为主，易于速效钾累积。OK和IDW的低值区较为接近，主要分布在研究区的北部及中部区域，RF的低值区则主要分布在研究区中部区域。研究区中部区域在3 种空间预测方法中均为低值区，其主要原因是受成土母质影响，成土母质主要是紫色岩类风化物和山河冲积物。这部分区域母质以风化作用为主，风化物松散，胶结性低，易遭侵蚀，风化淋溶较弱，发育的土壤主要属于幼年土壤，剖面发育差，剖面构型以耕作层（A）-底土层（C）为主，土体浅薄，仅10 cm左右，不利于速效钾积累。研究区北部区域主要受成土母质和地貌影响。成土母质基本为酸性结晶岩类风化物，这类母质发育的土壤质地多为壤土、土体厚度相对较厚，保水保肥能力较研究区中部区域强。地貌类型为山地，山区坡度较大，地形变化复杂，易发生水土流失。综合成土母质与地貌影响，在RF空间分布预测图中，研究区北部区域耕地土壤速效钾含量不高，但并不处于低值区。

3 讨论

通过筛选速效钾极大值（>200 mg·kg^-1），发现其极大值部分的R² 仍低于整体值的R²，极大值部分的MAE和RMSE则大于整体值的MAE和RMSE，说明极大值部分的预测精度较整体的预测精度低。推断可能在于局部区域受人为因素影响导致局部速效钾含量增加，诸如人为施肥、农田管理等。而本研究只选取了环境变量因素，并未考虑选取反映局部变量的人为因素作为辅助变量因子，后期可以通过增加一些能够反映局部信息的非环境变量因子来提升RF的拟合效果。

通过3 种空间预测方法的极大值和整体精度比较，可以得知，RF的预测结果能够精确地反映研究区耕地土壤速效钾分布情况，其精度最高。其本质在于RF属于多变量因素分析，OK和IDW属于单变量因素分析，而耕地土壤速效钾含量是受多种环境变量因素综合影响^[25]。因此，OK和IDW仅仅从空间位置关系上预测速效钾含量的算法限制了其预测精度。RF的整体分布趋势与OK、IDW接近，但RF通过大量的决策树和样本预测器构造模型，在空间位置关系的基础上，充分运用环境变量因子，很好地刻画了研究区耕地土壤速效钾含量与环境变量因子的非线性关系，所以RF的预测精度在3 种空间预测方法中最高。

4 结论

对比RF、OK、IDW 3 种方法对安庆市耕地土壤速效钾含量空间分布预测结果，3 种方法整体趋势一致。

从预测精度上看，RF的MAE和RMSE均小于OK和IDW，RF的R² 大于OK和IDW；OK的MAE和RMSE、R² 均为第二；IDW的MAE和RMSE均为最大，R² 最小。说明RF在研究区范围内的预测精度最高，OK次之，IDW预测精度最低，OK和IDW只能在一定趋势上反映出速效钾空间分布的特点。其原因在于OK和IDW仅仅以空间位置作为空间预测的变量因子，RF则充分考虑了环境变量因子，故精度最高。

通过移除某个环境变量因子，验证RF的RMSE变化比例，可以得知12 个变量因子中，纬度、年平均温度、成土母质、高程、经度、年平均降水量是影响RF精度的主要因素，剩余变量因子在RF模型中的权重会因协同重复而导致其权重降低。

参考文献