土壤中不同大小直径的矿物颗粒的组合状况称为土壤质地。土壤质地是土壤最基本的物理性质之一，也是十分稳定的土壤自然属性，对土壤性状如通气透水性、保水保肥性、保温导温性、养分含量以及耕作性状都有很大影响，质地相似的土壤往往具有类似的肥力特征。因此，土壤质地是土壤肥力最重要的指标^[1-2]。由于传统土壤机械组成测定方法复杂、用时长，以往土壤调查中，区域土壤质地分布成图较少。随着现代农业对区域土壤质量的信息化、精准化和智能化的管理要求，对土壤物理性质空间分布信息的需求日渐迫切。

基于 20 世纪 60 年代地统计学原理，发展了多种数字土壤预测制图方法^[3-4]。近 20 年来，随着机器学习方法的深入研究和应用，数字土壤制图方法越来越完善^[5-6]，对不同地形区域土壤物理性状的制图研究也逐渐增多。对于山区，随机森林法（RF）使用较多^[7-8]，原因在于机器学习算法可以处理数据之间的非线性关系并且数据只需做较少预处理^[9-10]，并且 RF 对样本量不敏感，适用性较广。Wang 等^[11]采用 RF、增强回归树（BRT）和经典地统计学回归克里格法（RK）模型，对黑河流域土壤质地进行制图，RF 模型具有更好的精度性能。申哲等^[12]在预测黄土区土壤机械组成中， RF 取得比 RK 和普通克里格法（OK）更优的预测精度。对于复杂地形的土壤养分制图，也有相似的研究结果^[13-15]。支持向量机（SVM）作为机器学习方法，在小样本的分类及回归方面表现卓越，但其要求辅助变量为连续型变量，无法将土壤类型、土地利用类型等与土壤属性关系密切的分类变量充分利用，影响辅助变量与土壤属性间的协同变化关系^[16]，目前该法在数字土壤属性制图中应用相对较少；神经网络模型是机器学习方法中使用较多的模型，对于变量类型适用性强，且对数据友好，能够很好地描述土壤属性与辅助变量间的非线性关系，但它对于样本数量敏感，模型的参数较多，不易确定，并且容易过拟合^[17]。在平缓地区，地形因子等通常不能有效用于预测土壤制图，所以地统计学制图方法应用最广泛，以 OK 和 RK 最具代表性。张世熔等^[18]在河北省曲周县冲击平原地区，对 124 个耕层土壤机械组成使用 OK 插值取得较好的预测结果。张世文等^[19]在北京市平谷区采用传统统计和地统计学相结合的方法预测土壤质地空间分布，OK 插值结果较好。对于平缓地区的土壤养分制图，也得到相似结果^[20]。地理加权回归（GWR）作为空间局部回归技术，也常用于平缓地区预测值图，它将数据的空间位置嵌入线性回归模型中，可以很好探测空间关系的非平稳性。但 GWR 本质上是局部地理回归的结合，也有研究认为，反复使用数据来估计模型参数，容易产生多重比较，并且模型易受共线性影响^[21]。随着研究的深入，近年有一些新的环境变量，用于土壤空间预测模型。Zhu 等^[22]开发降水后基于遥感图像的地表动态反馈模式（LSDF）用于推导环境协变量，可以帮助预测低起伏地区的土壤质地变化。Liu 等^[23] 根据地表对太阳辐射的动态反馈构造环境协变量，以粗略的分辨率预测这些地区的土壤空间变化，在土壤质地预测制图中提供了良好的结果。Loiseau 等^[24]在法国西部地质变化较大的区域预测表层土壤机械组成，使用伽马射线光谱替代母质图，提高了预测性能指标，并且地图的噪音更小、更易解释。Fathololoumi 等^[25]在预测土壤有机碳、砂粒含量、碳酸钙含量的空间分布中，加入多时相卫星图像（动态协变量）作为新的辅助变量，结果显示新的辅助变量改善了土壤性质的预测。目前还有研究选择高光谱数据^[26-27]等作为辅助变量进行预测，取得较优结果。

以往的研究大多在地形属性（或其他成土因素）变化较大的地区应用数字土壤制图，使用一个或两个模型来预测特定的土壤属性^[28-30]。对于山前冲积扇地区的土壤机械组成研究较少，且缺少地统计学方法与机器学习方法的对比研究。平原地区作为粮油果蔬主产区，集约化农业的精准管理对空间上精细的定量土壤属性信息有更迫切的需求^[31]。因此，本文选择地处山前冲积扇平原的河北省灵寿、行唐、曲阳 3 个县为研究区，选用目前较常用的基于对称对数比（SLR）转换的普通克里格法 （SLR-OK）、回归克里格法（SLR-RK）、随机森林法（SLR-RF）3 种模型方法，结合地形、土壤类型、地表温度、归一化植被指数等辅助变量，进行表层土壤机械组成的空间制图，通过精度对比，探究在山前冲积扇区土壤机械组成的最优模型。研究明确河北省西部山前冲积扇地区土壤机械组成变化的主要辅助信息。

1 材料与方法

1.1 研究区概况

灵寿、行唐、曲阳 3 个县相邻，位于石家庄以西，河北省中南部，太行山东麓浅山丘陵区与华北平原的交接地带。本文选取 3 个县高程 400 m 以下区域作为研究区，北纬 38°14′—38°57′，东径 113°57′—114°53′，面积 2636.4 km²，属半干旱、半湿润季风性气候，年平均气温 12℃，≥ 10℃积温 3968.5℃，平均降水量 480.8 mm，分配为南少北多。地势自西北向东南倾斜，依次形成丘陵、平原等梯级地貌特征。研究区土壤类型最大土类均为褐土（74%），其次为粗骨土（8.56%）、石质土（5.58%）、潮土（4.71%）等，各土壤类型面积占比见表1。采样点区域母质多为冲积母质和洪冲积母质，母质类型较单一。在 3 个县交界处分布有磁河和大沙河两条大流量河流，磁河为季节性河流，大沙河为常年河。土地利用类型主要为耕地、草地、林地，种植作物主要为玉米、花生。

1.2 样品采集与分析

采用网格法布点，在此基础上综合考虑地形地貌、土壤类型等因素，在西北部等地形复杂地区加密样点，共选取具有代表性的样点 164 个 （图1），样点平均间隔 5～6 km。因研究区域为 3 个县冲积平原，样点采集集中在 3 个县东南部，高程 400 m 以下，约占 3 个县面积 2/3。各土壤类型的采样点个数见表1。土壤样品采集时间为 2020 年 10 月、2021 年 9 月，采样深度为 0～20 cm，记录样点经纬度、地形地貌、土壤类型、土地利用类型、种植作物等信息。土壤机械组成采用湿筛-吸管法测定^[32]，粒级划分采用国际制 （砂粒 2～0.02 mm，粉粒 0.02～0.002 mm，黏粒 <0.002 mm）。

1.3 辅助变量来源及处理

环境变量是反映或影响土壤属性的主要因素。土壤机械组成与成土母质、土壤类型、地面高程、归一化植被指数和土地利用方式等自然和人文环境因子有关。由于小范围内气候因素变化较均一，因此本文不考虑降水量及气温影响。山区一般选用地形因子、土壤类型及母质等作为辅助信息；对于平缓地区，地形往往不是主要影响因子。在山前冲积扇区域，成土母质多为冲积或洪冲积母质，冲积母质地形部位要低于洪冲积母质；这个特征可由不同土壤类型表现，因此土壤类型可在一定程度上表示母质数据。地表温度近年开始用于低海拔地区的土壤空间预测研究。土壤由于收支的热量和热性质不同而发生温度的升降变化，砂土因含水量少，热容量小，昼夜温度变化较大；粘质土昼夜温度的变幅比砂土小；粉粒则介于两者之间。通过土温日变化的不同，可以预测土壤机械组成。因此，本文选择土壤类型、土地利用方式、高程及地形因子、地表温度和归一化植被指数作为辅助变量。

1.3.1 辅助变量来源

灵寿、行唐、曲阳 3 个县高程及地形因子来自 1∶50000 等高线矢量图；对于土壤类型，灵寿、行唐县采用 1∶50000 土壤图，曲阳县采用 1∶500000 土壤图；使用处于植物生长季（5—9 月）、且云量覆盖度少（小于 5%）的卫星影像计算归一化植被指数和地表温度。2019 年 7 月 23 日和 9 月 25 日的 Landsat8 OLI_TIRS 四景卫星影像、2015 年 10 月 27 日和 2016 年 4 月 28 日的 MYD11A1 1 km 地表温度、发射率产品和 MOD11A1 1 km 地表温度、发射率产品均来源于地理空间数据云网站（http：//www.gscloud.cn）； 2020 年 30 m 分辨率遥感解译土地利用图来源于国家基础地理信息中心全球地表覆盖数据产品服务网站（DOI：10.11769）。

1.3.2 辅助变量的提取及处理

1.3.2.1 地形因子提取

用 1∶50000 等高线矢量图在 ArcGIS 10.6 中生成 10 m 分辨率的数字高程模型 （DEM）栅格影像，使用 SAGA GIS 软件基于 DEM 提取 14 个地形因子：高程（Ele）、坡度（Slo）、坡向（Asp）、综合曲率（GC）、平面曲率（PLC）、剖面曲率（PRC）、切向曲率（TC）、纵向曲率 （LC）、横截面曲率（CSC）、流线曲率（FLC）^[33]、坡度坡长因子（LSF）^[34]、流量累积（FA）^[35]、地形湿度指数（TWI）^[36]、风力作用指数（WEI）^[37-38]。

1.3.2.2 地表温度及归一化植被指数提取

MODIS 传感器安装在极地轨道卫星 Aqua 和 Terra 上，Terra 为上午星，从北向南于地方时 10：30 左右通过赤道，Aqua 为下午星，从南向北于地方时 13：30 左右通过赤道^[39]。这种高时间分辨率使传感器能够在 1 d 内捕捉到地表发生的动态和过程，这为建立地表温度日变化变量提供了可能。使用在上午 1：30 至上午 10：30 至下午 1：30 地表温度观测在顶点 10：30 形成的弧度角作为一个协变量β 10：30，土壤质地的不同影响地表吸收太阳能量的速率，变量β4281030 表示不同土壤质地下地表温度上升的速率。然后使用协变量β4281030 乘以斜率，即协变量 lp4281030，它表示从上午 1： 30 到下午 1：30 的总变暖增加指数^[23]。在 MRT4.1 中对 MODIS LST 数据重投影并转换格式，提取地表温度（LST）数据，在 ArcGIS 10.6 中计算生成β10271030、lp10271030（2015 年 10 月 27 日数据）和β4281030、lp4281030（2016 年 4 月 28 日数据）。在 ENVI5.3 中对遥感影像进行镶嵌、计算，提取归一化植被指数（NDVI）。NDVI 选取日期为 2019 年 7 月 23 日和 9 月 25 日，分别命名为 NDVI723 和 NDVI925。

1.3.2.3 土壤类型和土地利用类型提取

土壤类型（ST） 和土地利用方式（LU）为类别变量，无法直接参与研究区土壤机械组成的回归分析。本文采用“算术平均值变换”^[40]，以不同类别变量下定量因变量的算术平均值代替该类别变量。研究区采样点共涉及 9 种土壤类型 （表1），3 种土地利用方式（耕地、林地、草地），使用算术平均值变换为 9 个和 3 个数值变量。由于有砂粒、粉粒、黏粒 3 个因变量，故土壤类型和土地利用方式分别转化为 3 组定量变量，用于构建砂粒、粉粒、黏粒的预测模型^[38]。算术平均值变换用动态思维（类别自变量与定量因变量的关系）建立起自变量各水平与定量结果变量之间数量关系，比传统的哑变量变换更具合理性^[41]。

1.4 预测方法

1.4.1 数据转换方法

土壤机械组成属于成分数据，具有闭合效应，即所有成分（元素含量）的总量等于 1（100%）^[42]。为了使数据满足正态分布，且使预测结果满足非负、定和为 1 的要求，本文首先对土壤机械组成进行对称对数比（SLR）转换^[43]，转换公式为：

转回公式：

式中，Zij 代表第 i 个样点第 j 种粒级的相对含量 （%）； $Z_{ij}^{\mathrm{\text{'}}}$代表第 i 个样点第 j 种粒级含量的转换值；D 表示成分数据的维度，本文中 D=3；δj 为常数，取研究区第 j 种颗粒除 0 外最小含量的一半。

1.4.2 普通克里格法

普通克里格法（OK）通过变异函数对变化不剧烈的数据建立空间模型进行插值，是通过采样点数据对未知样点使用距离加权平均法进行估计。实际上，OK 拟合的是 n 个已知样点属性间的线性关系^[44]。变异函数作为地统计学的主要研究工具，有块金值、基台值和变程 3 个主要参数。通过这 3 个参数反映了变量在空间随机性、结构性以及空间相关性和依赖性关系。空间变异性主要由随机性变异和结构性变异组成。如果块金值较大，表明在小于步长变化域的尺度下仍存在空间变异^[45]。

1.4.3 回归克里格法

回归克里格法（RK）是应用回归模型分离出土壤属性与辅助变量间的趋势项，然后对剔除趋势后的残差应用克里格法插值，最终将趋势项分析和克里格插值结果结合起来的一种空间插值方法^[44]。该法的优点是考虑到了辅助变量对土壤属性的影响。常用的拟合趋势项方法有逐步线性回归模型、回归树模型等。残差项的拟合常采用 OK 或简单克里格法（简单克里格法与 OK 类似，但其假定数学期望是已知的）。其中，剔除趋势模型的选择很重要^[18]，同一地区选择不同的趋势模型往往会得出不同的结果。本文选用逐步线性回归模型来剔除趋势，首先用皮尔森相关分析对辅助变量进行初步筛选，为避免辅助变量间的共线性问题，基于筛选出的变量，采用逐步线性回归的方法对经 SLR 转换后的土壤机械组成进行预测（F 值的概率 P<0.05 为变量进入方程的标准，P>0.10 为剔除变量的标准），并对预测残差进行 OK 插值，最后将回归预测值与残差相加并转回得到预测结果。

1.4.4 随机森林法

随机森林法（RF）属于集成学习的范畴，集成学习下有两个重要的策略 Bagging 和 Boosting。 Bagging 算法：每个分类器都随机从原样本中做有放回的采样，分别在这些采样后的样本上训练分类器，再把这些分类器组合起来，其代表算法是 RF^[46]。RF 通过提供变量重要性的措施，提供了一些解释模型的能力^[47-48]。在 RF 中生成的树本身并不是可解释的，因为它们是与其他树一起作为一个集合体来产生预测的。相反，在训练完 RF 模型后，会计算出协变量的使用频率汇总，并作为各自变量重要性的衡量标准进行报告。本研究使用 K 折交叉验证法（CV）来评估地图的性能^[49-50]。CV 解释了使用不同的点进行训练时，结果会有多大变化，它通过测试结果对选择哪些点进行训练的依赖性来说明模型的稳健性。

RF 不要求原始数据符合正态分布，但为了更好地同其他两种方法做对比，且使预测结果满足定和为 1 的要求，本研究使用经 SLR 转换后的数据进行建模。由于 RF 不是简单的线性回归，皮尔森相关分析筛选变量的结果不一定适合 RF 模型，故用 CV 对 RF 辅助变量进行筛选：对变量进行重要性排序，初步剔除不重要变量^[51-52]；对剩下的变量，基于 100 次 10 折交叉验证辅助评估最优数量的变量。交叉验证曲线展示了模型误差与用于拟合的变量数量之间的关系。根据交叉验证曲线，保留重要的变量即可以获得理想的回归结果，因为此时的误差达到最小（变量筛选过程中，RF 参数均为默认参数）。基于筛选的辅助变量对 RF 模型进行参数优化：首先以均方根预测误差（RMSE）最小化为目标，通过逐次计算确定决策树节点分裂时所用变量个数（mtry）；通过观察决策树数量（ntree）与误差关系图，在保证误差稳定的前提下，选择较小的 ntree^[53]。

1.5 精度验证

从 164 个样点中随机选取 70% 的样点（114 个）作为训练集用于空间预测，30% 的样点（50 个）作为验证集用于精度验证，预测样点和验证样点分布见图1。此比例是使用 RF 算法中的交叉验证所得。本研究采用验证样点的平均绝对误差 （MAE）、RMSE、决定系数（R²）进行模型精度评价。MAE、RMSE 描述了模型的准确性，意味着预测值与实际值的接近程度。MAE 是一个线性分数，意味着所有的残差都是同等权重。RMSE 在计算其平均值之前对误差进行了平方，所以大的误差被赋予相对较高的权重。R² 指标表示模型所能解释的目标变量的方差比例，可以用这个指标来比较不同目标变量之间的模型性能。R² 值越接近于 1，对因变量的解释能力越强。MAE、RMSE 越小，模型精度越高。

1.6 建模与制图工具

本文中 RF 在 R3.6.2 中完成，相关性分析和逐步回归使用 SPSS 26 软件，OK 插值、回归残差的 OK 插值和空间制图在 ArcGIS 10.6 中完成，遥感影像镶嵌和 NDVI 的提取在 ENVI 5.1 软件中完成， SLR 转换后砂粒和粉粒数据的 BOX-COX 转换在Python 3.6 中完成。

2 结果与分析

2.1 土壤机械组成的描述性统计特征

表2 为研究区采样点土壤机械组成的描述性统计结果，可以看出研究区表层土壤砂粒含量最高（训练集和验证集砂粒平均含量分别为 72.00%、71.38%），粉粒含量最低（训练集和验证集的平均粉粒含量分别为 14.09%、13.39%），黏粒含量居中（训练集和验证集的平均黏粒含量分别为 14.39%、14.12%）。训练集中标准差比验证集中的大，相对应地，训练集中砂粒、粉粒、黏粒含量的变异系数（14.04%、41.53%、39.00%）大于验证集（13.93%、38.40%、38.75%）。方差分析结果表明训练集与验证集土壤 3 个粒级含量差异均不显著，表明训练数据和验证数据都具有较好的代表性。训练集中 3 个粒级含量偏度均接近于 0，但经 Kolmogorov-Smirnov 检验，只有黏粒含量 P 值大于 0.05，满足正态分布，粉粒和砂粒含量均不符合正态分布。经 SLR 转换后，粉粒和砂粒依然不符合正态分布，对粉粒和砂粒进行 BOX-COX 转换后，符合正态分布，可用于空间插值和线性回归。

2.2 OK 模型、RK 模型及 RF 模型的构建

对训练集数据进行 Kolmogorov-Smirnov 检验，经 SLR 转换后的黏粒数据 P 值大于 0.05，满足正态分布，SLR 转换后的砂粒和粉粒数据 P 值小于 0.01，均不满足正态分布，对其进行 BOX-COX 转换后 （slrsand 幂参数：-0.3838384；slrsilt 幂参数：-0.2222222）， P 值大于 0.05，满足正态分布。

由表3 可知，OK 插值中，砂粒和粉粒拟合模型的块金系数具有强烈的空间相关性，黏粒处于中等程度的空间变异，受结构性因素和随机性因素共同影响。可用于空间插值分析。变程是指变量空间自相关的距离，从表3 中可以看出，3 个粒级中最小值为 11709 m，大于本研究区的采样平均间距 5000 m，因此满足该区对土壤机械组成空间预测的要求。

由表4 可以看出，高程、土壤类型和风力作用指数 3 个环境变量进入线性回归模型，从调整决定系数来看，砂粒、粉粒和黏粒含量的线性回归模型方差解释率依次为 0.374、0.353、0.330。其中，砂粒拟合度最高，黏粒拟合度最低。

注：slr-sand、slr-silt、slr-clay 分别表示砂粒、粉粒、黏粒；Ele、ST、WEI 分别表示高程、土壤类型、风力作用指数。下同。

表5 为 RF 模型的拟合结果，对于砂粒，进入模型预测的辅助变量有土壤类型、高程、NDVI、地表温度、流量累积和风力作用指数；对于粉粒，进入模型预测的辅助变量有高程、土壤类型、地表温度、风力作用指数；对于黏粒，进入模型预测的辅助变量有风力作用指数、土壤类型、地表温度。因为 RF 可以拟合土壤机械组成和辅助变量之间的非线性关系，所以进入 RF 模型的变量相较于回归克里格有所增加。其中，对砂粒和粉粒含量的空间预测来说，高程和土壤类型是最重要的辅助变量，地表温度、风力作用指数的重要性相对较低，对于砂粒，流量累积变量也参与空间预测。对于黏粒，本研究使用随机森林的十折交叉验证来解释变量数量与误差值的关系（图2），结果显示，黏粒变量数量与误差值关系图为先升后降再上升的趋势，并且在变量数量最少时误差最低。由于随机森林预测需要变量的支持，所以选择关系图中先上升后下降的下降最低点，即相对误差较低处的变量数量为 4 个，作为黏粒预测中变量的数量选择。

注：砂粒、粉粒、黏粒为经对称对数比转换之后的数值。FA 表示流量累积。RF 建模变量按重要性由高到低列出。

2.3 土壤机械组成的空间分布预测结果

利用拟合的 SLR-OK、SLR-RK 和 SLR-RF 模型分别对研究区土壤砂粒、粉粒、黏粒含量的空间分布进行预测，得到其空间分布图（图3～5）。可以看出，3 种方法预测的土壤各粒级含量空间分布的趋势大体一致。研究区砂粒含量呈现为东南部高、西北部低的趋势，大部分区域含量为 60%～70%，其中在行唐县和曲阳县交界的大沙河流域含量较高，达到 85% 左右，在灵寿县和行唐县交界的磁河流域的横山岭水库出口和出县界处砂粒含量较高，因为磁河是季节性河流，并且河道较窄，在出县界处逐渐变宽。灵寿县西南部临近滹沱河流域及黄壁庄水库，曲阳县东南部临近唐河流域，这两处砂粒含量较高。粉粒和黏粒含量的空间分布趋势与砂粒相反，呈现东南部低，西北部高的趋势，在河流流域处，砂粒和黏粒含量在 10% 左右，其他区域为 15%～20%。3 种方法预测的土壤机械组成空间分布总体趋势比较相似，但在局部区域上不尽相同。SLR-RK 和 SLR-RF 的预测结果对细节刻画比较清晰，SLR-OK 预测结果只能反映研究区土壤机械组成的宏观分布趋势，难以体现局部变异信息。

2.4 不同预测方法的精度对比

表6 列出了 3 种方法在验证集的预测精度，对于砂粒，从 MAE 和 RMSE 来看，SLR-RF<SLR-RK<SLR-OK，从 R² 来看，SLR-RF>SLR-RK>SLR-OK，表明 SLR-RF 的解释能力和预测精度最高；对于粉粒，从 MAE 和 RMSE 来看，SLR-RF<SLR-OK<SLR-RK，从 R² 来看，SLR-RF>SLR-OK>SLR-RK，表明 SLR-RF 的解释能力和预测精度最高；对于黏粒，从 MAE 和 RMSE 来看，SLR-OK<SLR-RF<SLR-RK，从 R² 来看，SLR-OK>SLR-RF>SLR-RK，表明 SLR-OK 的解释能力和预测精度最高。从以上 3 个参数综合来看，SLR-RF 对山前冲积扇区砂粒和粉粒的预测精度最高；SLR-OK 对黏粒的预测精度最高。

3 讨论

本研究利用 3 种方法对河北中南部山前冲积扇区县域中尺度的土壤机械组成空间分布进行预测。在砂粒和粉粒中，SLR-OK 法预测精度较低，SLRRK 法预测精度居中，SLR-RF 法预测精度最高。本研究采样点集中在山前冲积扇平原区，区内有一定的地形起伏，河流众多，山前上部位置主要为质地相对轻的洪积物，下部平原地势趋缓，冲积物中细颗粒增多，表现为砂粒和粉粒含量对地形变化较敏感；且受水域冲击影响较大。SLR-RK、SLR-RF 法引入辅助变量参与预测，提高了砂粒的空间预测精度。对于模型中辅助变量的选择，均由模型的相应指标确定，RK 中采用逐步线性回归选取变量， RF 中通过交叉验证确定变量（图2）。并且可以看出，RK 通过对回归残差进行普通克里格插值进一步提高了多元线性回归的预测精度，其中地表温度、NDVI 和流量累积变量未进入 RK 的多元线性回归模型，表明这 3 个环境变量与土壤属性的线性相关性不高；而 RF 可以挖掘土壤属性同辅助变量之间的非线性关系，并且地表温度、NDVI 和流量累积变量参与了 RF 模型预测，表明这 3 个环境变量与土壤属性存在非线性相关性，非线性环境变量的加入，提高了模型预测精度，并且取得了比 RK 更好的预测效果。在黏粒中，SLR-OK 法预测精度最高，黏粒处于中等程度的空间变异，受结构性因素和随机性因素共同影响。RK 和 RF 的预测精度低于 OK 的原因可能是黏粒属性和研究选用的辅助因子之间的相关性不高。由 RF 的辅助变量数量与误差值的关系图（图2）可以看出，黏粒预测中，辅助变量个数为 0 时预测误差最小，即不使用所选环境变量进行预测时，误差达到最低，可能原因是本文所选环境变量与研究区土壤黏粒属性关联度不高，选用环境变量后影响黏粒属性的空间自相关分布，而 OK 在地形变化平缓的地区具有较好的预测能力，并且可以最大限度表达土壤属性空间自相关性，提高预测精度。从地理位置分析看，研究区处于山前冲积扇平原区，虽然地形变化平缓，但因地势相对较高，与下游平原上有一定高差，且无大范围的洼地，因此以较快速度随河流冲出山区的冲积物中，较粗颗粒先沉积，较细的颗粒随水流移动距离更远。土壤黏粒含量受研究区冲积扇地形变化影响较小。黏粒预测趋势为模型检验结果，可能还需要实地调研。

在山前冲积扇地区，对于粉粒和砂粒，即使在地形变化平缓地区，地形属性也是预测土壤性质非常有用的辅助信息，这与 Liang 等^[54]使用基于案例的数字土壤测绘协变量选择方法的结果一致； Mosleh 等^[31]在数字土壤测绘预测低地势地区土壤特性的有效性中，也表明地形属性是不同辅助信息中的主要预测因子。并且在加入地表温度变量后，预测精度有一定提升。对于山前冲积扇地形变化平缓地区，未来需要进行更多的调查来确定新的辅助变量，同时引入发现土壤数据和土壤变异性之间复杂关系的方法和手段，以提高机器学习算法在低海拔地区的表现。

4 结论

本研究在河北太行山麓山前冲积扇地区，对 114 个采样点，采用 OK、RK 和 RF 结合地形因子、土壤类型、NDVI、地表温度变量等辅助变量，对该区 0～20 cm 土层机械组成空间分布进行预测，并以 50 个验证点进行验证。主要结论如下：

（1）从空间预测图来看，砂粒呈现出西北低，东南高的空间分布趋势；粉粒和黏粒与砂粒相反。与 SLR-OK 相比，SLR-RK 和 SLR-RF 能够更好地反映局部变异并减小平滑效应。

（2）对比 3 个预测模型的精度，对于砂粒和粉粒，SLR-RF 法对验证集含量预测的 MAE 和 RMSE 均低于其他两种方法，且从 R² 来看，均高于其他两种方法，表明 SLR-RF 的预测精度最高；对于黏粒，SLR-OK 法对验证集含量预测的 MAE 和 RMSE 均低于其他两种方法，且 R² 最高，表明 SLR-OK 的预测精度最高。

（3）从土壤机械组成的辅助变量来看，线性回归预测模型的辅助变量包括高程、土壤类型和风力作用指数；RF 模型的辅助变量包括高程、土壤类型、NDVI、地表温度、风力作用指数和流量累积，对于砂粒和粉粒，土壤类型和高程是重要的辅助变量， NDVI、地表温度、风力作用指数和流量累积重要性相对较低；对于黏粒，通过空间自相关预测模型 SLR-OK 结果为最优。

参考文献