土壤是地球生命和农业生产的重要基础资源，了解其类型及空间分布对于土壤资源的可持续利用和管理至关重要^[1]。随着计算机科学，遥感、数据挖掘算法的快速发展，数字土壤制图(Digital Soil Mapping，DSM)成为一种新兴的、高效表达土壤空间分布的技术方法。该方法运用数据挖掘算法建立土壤-景观关系模型，以栅格数据为表达方式，在计算机环境下辅助制图^[2-3]。

常用的数字土壤制图方法主要分为地统计法和机器学习算法，线性模型^[4]、克里格法^[5]和判别分析^[6]为常见的地统计法，由于土壤与环境之间的关系过于复杂，为了更好地挖掘土壤-景观关系，机器学习算法被引入到构建土壤类型或土壤属性制图模型中，如决策树(Decision Tree，DT)^[7]、随机森林(Random Forest，RF)^[8]、梯度提升树(Gradient Boosting Decision Tree，GBDT)^[9]和极端梯度提升(Extreme Gradient Boosting，XGBoost)^[10]等。已有较多研究对比了多种机器学习算法在土壤属性制图方面的效果^[11-12]，但在土壤类型制图中尚不多见。地形因素是土壤形成与发育的主要影响因素之一，在不同的地形条件下，土壤的形成过程和发育阶段会有所不同，从而导致土壤类型和性质差异^[13]。特别是在平原-丘陵混合区域，地形的多尺度变化使得土壤类型的空间异质性难以精确捕捉，进而土壤类型推理制图精度不高。此外，土壤类型制图精度还受到多种环境因素的综合影响，这增加了平原-丘陵区域推理制图的复杂性。因此，综合考虑地形条件和土壤-景观关系模型的建立过程，比较不同数字土壤制图模型在平原-丘陵区域的土壤类型推理制图效果，并选出最优的数字制图模型，显得尤为重要。

本研究以湖北省麻城市乘马岗镇北部平原-丘陵区域为研究区，将其划分为平原和丘陵两个地形单元，通过多种机器学习算法优选各地形的制图方法，依据RF模型内嵌的因子重要性排序，筛选出最优环境变量集，分别构建各地形最优土壤-景观关系模型进行土壤类型推理制图。接着，对比整体与按地形分区制图的精度，探究制图方法优选、环境变量筛选以及制图方式对平原-丘陵区域土壤类型推理制图精度的影响，以期为复杂地形区域土壤类型推理制图提供实践案例和技术支持。

研究区位于湖北省麻城市乘马岗镇北部(图 1)，中心地理坐标为31°25′30″N，114°58′30″E，总面积约为28.14 km²，山地、平畈交错，属于典型的平原-丘陵混合区域，高程范围为80~346 m，为探究制图方法和参与制图的环境变量在整体区域、平原区域和丘陵区域的差异，依据土壤志及地形起伏度分布图，以150 m高程为分界将整体研究区分为平原区域和丘陵区域。

图 1 研究区概况 Fig. 1 Overview of the study area

根据第二次全国土壤普查成果的1：50 000传统土壤类型图可知，研究区内共有9种土壤类型和3种母质类型，其关系如表 1所示。

表 1 Table 1

表 1 研究区土壤类型、母质类型、训练点和验证点信息 Table 1 Information of soil type, parent material type, training point and validation point in the study area 序号 Number 土壤类型 Soil type 母质类型 Parent material type 土壤类型代号 Soil type code 训练点 Training points 验证点 Validation points 1 沙泥土 Silt soil 花岗岩-片麻岩的坡积或残积物 Slope deposit or remnant of granite gneiss 1 566 16 2 硅沙泥土 Silica silt soil 白云石英片岩-花岗片麻岩坡积或残积物 Dolomite quartz schist granite-gneiss slope deposit or remnant 4 607 25 3 林地泥沙土 Forest silt soil 花岗岩-片麻岩的坡积或残积物 Slope deposit or remnant of granite gneiss 1-7 97 9 4 沙泥田 Silt field 花岗岩-片麻岩的坡积或残积物 Slope deposit or remnant of granite gneiss 22 192 14 5 硅沙泥田 Silica silt field 白云石英片岩-花岗片麻岩坡积或残积物 Dolomite quartz schist granite-gneiss slope deposit or remnant 24 248 19 6 潮土田 Tidal soil field 近代河流冲积物 Modem river alluvium 30 48 13 7 硅沙土 Silica sand soil 白云石英片岩-花岗片麻岩坡积或残积物 Dolomite quartz schist granite-gneiss slope deposit or remnant 4-4 152 12 8 硅麻骨土 Silica hemp bone soil 白云石英片岩-花岗片麻岩坡积或残积物 Dolomite quartz schist granite-gneiss slope deposit or remnant 4-5 116 11 9 林地薄层硅麻骨土 A thin layer of silica hemp bone soil in woodland 白云石英片岩-花岗片麻岩坡积或残积物 Dolomite quartz schist granite-gneiss slope deposit or remnant 4-6 434 22

表 1 研究区土壤类型、母质类型、训练点和验证点信息 Table 1 Information of soil type, parent material type, training point and validation point in the study area

数字土壤制图的关键在于构建土壤-景观关系模型，该模型不仅受制图算法的影响，还需要选取适宜的环境变量以保证其稳健性。基于SCORPAN土壤景观模型^[14]，本文选取了28个环境因子来获取土壤-景观知识，其中包括1个母质信息、12个地形因子和15个遥感因子^[15-16](表 2)。母质信息来源于成土母质图及当地土壤专家所给意见；地形因子中的高程数据由湖北省测绘局提供的10 m等间距等高线数据插值所得，地形衍生因子利用ArcGIS 10.5和SAGA GIS 9.4.1软件计算获得；遥感因子中的植被指数通过谷歌地球引擎(Google Earth Engine，GEE)平台中的Sentinel-2影像(成像时间2020年2月9日)计算获得，同时从GEE获取同期的Sentinel-2多波段数据，包括B2、B3、B4、B5、B6、B7、B8、B8a、B11和B12。为了减少数据冗余，采用ENVI 5.3对其进行主成分分析得到第一主成分，进而提取纹理特征。为解决各环境变量的坐标系和尺度不一致问题，所有环境变量均被标准化为10 m分辨率，并统一投影坐标系为WGS_1984_ UTM_Zone_50N。

表 2 Table 2

表 2 环境变量基本信息 Table 2 The basic information about environment variables 环境变量Environment variables 空间分辨率Resolution/m 简称Abbreviation 母质信息Parent material 母质Parent material - PM 地形因子 高程Elevation 10 EL Topographic 坡向Aspect 10 Aspect factor 坡度Slope 10 Slope 谷深Valley depth 10 VD 相对坡度位置Relative slope position 10 RSP 坡度坡长因子Slope length and steepness factor 10 LS 水平曲率Horizontal curvature 10 HOC 平面曲率Plan curvature 10 PLC 剖面曲率Profile curvature 10 PRC 地形湿度指数Topographic wetness index 10 TWI 地形起伏度Topographic relief 10 TRF 地形粗糙度Topographic roughness 10 TR 遥感因子 第一主成分First principal component 10 FPC Remote sensing 均值Mean 10 Mean factor 方差Variance 10 VA 协同性Homogeneity 10 HO 对比度Contrast 10 CO 相异性Dissimilarity 10 DI 信息熵Entropy 10 EN 二阶矩Second moment 10 SM 相关性Correlation 10 CR 归一化植被指数Normalized difference vegetation index 10 NDVI 增强植被指数Enhanced vegetation index 10 EVI 差值植被指数Differential vegetation index 10 DVI 比值植被指数Ratio vegetation index 10 RVI 绿光归一化植被指数Green normalized difference vegetation index 10 GNDVI 土壤调整植被指数Soil adjusted vegetation index 10 SAVI

表 2 环境变量基本信息 Table 2 The basic information about environment variables

本研究中的样本点包括2 460个训练点和141个验证点(图 1b)，其中训练点是依据传统土壤图，使用ArcGIS 10.5软件中的“创建随机点”工具生成，设置每种土壤类型不少于300个，样点最小允许距离为50 m。为提高模型的准确性和泛化能力，本研究使用所有训练点，并充分利用数据资源，以更准确地比较不同制图方法的效果。验证点采用完全随机采样、主观采样、混淆区域边界采样和均匀(网格)采样相结合的方式进行实地采集，采样时间为2020年10月-2021年6月，覆盖整个研究区所有土壤类型且空间分布均匀。研究区各土壤类型训练点及验证点数量如表 1所示。

现有研究表明，过多的环境因子不仅会导致数据冗余，而且会增加计算复杂度，降低模型精度^[17]。为了降低模型的计算成本，增加模型鲁棒性，本研究采用逐步法来进行环境变量筛选。首先，基于RF内嵌模型计算所有环境变量的因子重要性，通过逐一添加正向排序后的环境变量进行土壤-景观关系模型的构建，并计算每次预测制图的总体精度(Overall Accuracy，OA)。最终，选择OA最大的环境变量组合作为最优环境变量集。

本文采用决策树、随机森林、梯度提升树和极端梯度提升4种机器学习算法进行土壤-景观关系模型的构建。

决策树是一种由树模型构成的机器学习算法，易于理解和解释。它根据特征将信息划分成两个或多个子节点，逐步构成决策树，最终选择的结果由树最远的节点所决定^[18]。决策树建立流程一般包括3个步骤：特征选择、决策树生成以及剪枝。信息增益和基尼系数是两种特征选择方法，本文采用基尼系数法。

随机森林是一种以决策树为基本分类器的机器学习算法，能够处理高维数据并自动估算特征重要性。该算法通过在特征集的不同随机子集中构建多个决策树，在每个决策树的每个节点上利用一定的特征选取方法选出最优特征进行划分，取多个预测结果的平均值作为最终预测结果^[19]。

梯度提升树是一种基于决策树的迭代集成学习算法，适用于处理复杂数据关系。该算法利用梯度下降法进行迭代，即将前一个树模型拟合的残差作为下一棵树的输入，以实现损失函数最小化，从而提升模型预测性能，最终将多个决策树的预测结果相加作为最终预测^[20]。

极端梯度提升是在梯度提升树基础上改进的一种机器学习算法^[21]。该算法与梯度提升树类似，但在损失函数后加入了正则化项，既能降低模型复杂度，提升精度，又具备准确性高、训练速度快、防止过拟合和适合大规模数据等特点^[22]。

本文所有模型均在Python 3.9环境下实现，并在使用贝叶斯优化算法进行各模型参数优化时，设定随机搜索步数(init_points)为15，贝叶斯优化步数(n_iter)为500。

本研究通过验证点与推理土壤图之间的混淆矩阵，计算总体精度、Kappa系数、生产者精度(Producers’ accuracy，PA)和用户精度(Users’ accuracy，UA)4种评价指标，用于评估不同制图方法的推理精度^[23]。Kappa系数用于衡量分类一致性，取值范围为0~1，其中0~0.20为微弱一致性，0.21~0.40为弱一致性，0.41~0.60为中等一致性，0.61~0.8为显著一致性，0.81~1为几乎完全一致性。

基于所有环境变量，本文分别采用DT、RF、GBDT和XGBoost对整体、平原和丘陵区域进行推理制图。如图 2所示，4种算法在三种区域均可以实现推理制图，但不同算法在不同区域的制图精度有所不同。除DT算法外，其他三种算法在丘陵区域的制图精度均高于平原区域和整体区域。整体区域中，制图精度由高到低分别为RF、XGBoost、GBDT和DT，其中RF算法的总体精度为68.09%，Kappa系数为0.63；平原区域中，RF算法也取得了最高的制图精度，其总体精度和Kappa系数分别为68.57%和0.62，紧接其后，制图精度由高到低分别为DT、XGBoost和GBDT；丘陵区域中，RF、GBDT和XGBoost算法的制图精度均高于70%，其中XGBoost算法制图精度最高，其总体精度为73.24%，Kappa系数为0.68。最终，分别选择各区域最优制图方法用于后续研究。

图 2 各区域推理制图方法对比 Fig. 2 Comparison of inference mapping methods in each region

基于随机森林内嵌算法，分别计算整体、平原和丘陵区域所有环境变量的因子重要性，通过逐一添加正向排序的环境变量用于各区域优选制图方法并进行精度评价。由图 3可以看出，随着正向排序环境变量数目的增加，三种区域的优选制图方法精度均呈现先上升后平缓的趋势。整体区域制图精度在环境变量个数为8时最高，总体精度为73.05%；平原区域制图精度最高为70.00%，在环境变量个数为17时得到；当丘陵区域的环境变量个数为11时，制图精度达到最高，为76.06%。

图 3 各区域优选制图方法总体精度随环境变量数目的变化 Fig. 3 Variation of overall accuracy for optimal mapping methods in each region with the number of environment variables

将各区域所有变量和筛选后变量用于各区域最优制图方法，进行土壤类型推理制图，结果如表 3所示。与使用所有变量相比，使用筛选后变量，各区域制图精度都有了一定的提升。其中，整体区域制图精度提升效果最好，总体精度和Kappa系数分别提高了4.96%和0.06；其次是丘陵区域，总体精度提升了2.82%，Kappa系数增加了0.03；平原区域制图精度的提升效果最小，总体精度和Kappa系数分别增加了1.43%和0.02。因此，环境变量筛选能够有效提升各区域土壤类型推理制图精度。本研究分别选取正向排序前8、17和11个环境变量作为整体、平原和丘陵区域的最优环境变量集。

表 3 Table 3

表 3 变量筛选前后优选制图模型的性能 Table 3 The performance of optimal mapping models before and after variable selection 优选模型Optimal model 所有变量All variables 筛选后变量Filtered variables 总体精度 卡帕系数 总体精度 卡帕系数 Overall accuracy/% Kappa coefficient Overall accuracy/% Kappa coefficient 整体-随机森林Whole-RF 68.09 0.63 73.05 0.69 平原-随机森林Plain-RF 68.57 0.62 70 0.63 丘陵-极端梯度提升Hill-XGBoost 73.24 0.68 76.06 0.72

表 3 变量筛选前后优选制图模型的性能 Table 3 The performance of optimal mapping models before and after variable selection

根据随机森林内嵌模型计算的各区域最优环境变量集中环境变量的因子重要性，进行归一化处理，结果如图 4所示。母质是整体、平原和丘陵区域中影响力相对最大的环境变量，但其重要性有所不同，在整体、平原和丘陵区域的重要性依次为0.35、0.20和0.32。谷深、高程、相对坡度位置、坡向和坡度坡长因子为三个区域均必需的地形因子，平原和丘陵区域还需要的地形因子为地形湿度指数；遥感因子中的均值为三个区域所必需的，第一主成分为整体和平原区域所必需的因子，丘陵区域还需要的遥感因子为比值植被指数、土壤调整植被指数和归一化植被指数。与整体和丘陵区域相比，平原区域需要更多的遥感因子，包括绿光归一化植被指数、增强型植被指数、相关性、差值植被指数和方差，这可能是由于平原区域地形较为平缓，导致土壤类型分布的地带性特征不明显，仅依靠地形因子和少量遥感因子无法提供足够的信息来进行推理制图。

图 4 各区域筛选后环境变量的因子重要性 Fig. 4 Feature importance of filtered environment variables in each region

基于各区域最优制图方法和最优环境变量集，构建各区域最优土壤-景观关系模型，进行各区域土壤类型推理制图。将平原区域和丘陵区域土壤类型推理制图结果合并为按地形分区推理制图结果，并与整体区域土壤类型推理制图结果进行对比，结果如图 5所示。在使用所有变量进行推理制图时，整体推理制图的总体精度为68.09%，Kappa系数为0.6305，按地形推理制图的总体精度和Kappa系数分别为70.92%和0.6643，各自提高了2.83%和0.0338；在使用筛选后变量进行推理制图时，虽然按地形推理制图和整体推理制图的总体精度一致，皆为73.05%，但Kappa系数由整体推理制图中的0.6895提升到按地形推理制图中的0.6897。综上，按地形分区的制图方式也能有效提升平原-丘陵区域土壤类型推理制图精度，但在进行环境变量筛选后，其提升效果不明显。本研究选取基于筛选后变量进行的按地形推理制图为研究区最优推理制图结果。

图 5 不同推理制图方式之间的对比 Fig. 5 Comparison of different inference mapping ways

将原始土壤图与研究区最优推理土壤图进行对比分析，根据图 6可知，相比较于原始土壤图，推理土壤图的空间分布更加破碎，展现出更多的细节信息。整体上，两者空间分布大体一致，其中潮土田在两者中的空间分布完全一致，这可能是因为该土壤是唯一由近代河流冲积物发育而来。然而，部分土壤类型在推理土壤图中的空间分布范围与原始土壤图相比，发生了较大变化，如硅沙土和硅麻骨土，该两种土壤类型在推理土壤图中的分布面积明显减少，并且很大一部分被推理为硅沙泥土，这可能是由于硅沙土、硅麻骨土和硅沙泥土均由白云石英片岩花岗片麻岩坡积或残积物发育而来，成土环境类似，且在研究区中往往相伴而生，具有较强的相似性。

图 6 原始土壤图与最优推理土壤图 Fig. 6 Original soil map and optimal inference soil map

为了进一步评估最优推理土壤图的分类准确性，表 4给出了各土壤类型在推理土壤图中的生产精度和用户精度。生产精度表示验证点土壤类型在推理土壤制图中被正确分类的概率，用户精度则表示在推理土壤制图中，落在该土壤类型上的验证点被判定为该土壤类型的概率。根据推理制图的原理，土种级别的推理制图相比较于土属级别更为复杂，因此难以取得较高的精度。通过混淆矩阵可以看出，所有土壤类型的用户精度均较为理想(均超过60%)；除林地泥沙土、硅沙土和硅麻骨土外，其他土壤类型的生产精度均超过70%。尽管林地泥沙土和硅麻骨土的生产精度仅为33%和27%，但它们的用户精度均为100%，这表明在推理土壤图中，其他土壤类型不会被错误分类为这两种类型。

表 4 Table 4

表 4 实地验证点与最优推理土壤图间的混淆矩阵 Table 4 Confusion matrix between validation points and optimal inference soil map 实地验证点 Validation points 最优推理土壤图Optimal inference soil map 沙泥土 Silt soil 硅沙泥土Silica silt soil 林地泥沙土 Forest silt soil 沙泥田 Silt field 硅沙泥田 Silica silt field 潮土田 Tidal soil field 硅沙土 Silica sand soil 硅麻骨土 Silica hemp bone soil 林地薄层硅麻骨土 A thin layer of silica hemp bone soil in woodland 合计 Total 生产精度Producers’ accuracy /% 沙泥土 Silt soil 13 0 0 3 0 0 0 0 0 16 81 硅沙泥土 Silica silt soil 0 22 0 0 2 0 1 0 0 25 88 林地泥沙土 Forest silt soil 2 0 3 4 0 0 0 0 0 9 33 沙泥田 Silt field 4 0 0 10 0 0 0 0 0 14 71 硅沙泥田 Silica silt field 0 1 0 0 17 0 0 0 1 19 89 潮土田 Tidal soil field 0 0 0 0 0 13 0 0 0 13 100 硅沙土 Silica sand soil 0 5 0 0 0 0 5 0 2 12 42 硅麻骨土 Silica hemp bone soil 0 3 0 0 5 0 0 3 0 11 27 林地薄层硅麻骨土 A thin layer of silica hemp bone soil in woodland 0 1 0 0 4 0 0 0 17 22 77 合计Total 19 32 3 17 28 13 6 3 20 141 用户精度Users’ accuracy /% 68 69 100 59 61 100 83 100 85

表 4 实地验证点与最优推理土壤图间的混淆矩阵 Table 4 Confusion matrix between validation points and optimal inference soil map

本研究首先对比了DT、RF、GBDT和XGBoost算法在整体、平原和丘陵区域的推理制图效果，然后基于因子重要性进行了环境变量筛选，并比较了整体和按地形推理制图结果。最后，结合优选的制图方法、筛选后的环境变量和按地形制图方式，进行了研究区土壤类型推理制图。研究结果为平原-丘陵区域推理制图中的环境变量筛选、制图算法选择和制图方式构建提供了一定的技术参考，有助于提升推理制图精度。此外，利用该方法生成的土壤类型图可以作为参考，与第二次全国土壤普查成果的土壤类型图进行对比，以识别潜在的错分区域，并为实地调查提供建议。

目前已有许多研究基于机器学习算法开展数字土壤制图，但仍存在一些不足，例如，杨雨菲等^[24]利用多种机器学习算法进行推理制图时未考虑环境变量的重要性。张晓婷等^[25]虽然考虑了环境变量在推理制图中的重要性，但未探讨不同机器学习算法在不同区域的表现差异。本文通过对比分析发现，RF算法在整体和平原区域的推理制图效果最好，丘陵区域制图精度最高的为XGBoost算法。尽管四种机器学习算法在不同地形中的推理制图能力不同，相对而言，RF与XGBoost算法的推理制图能力均较好(图 2)，与先前所得结论一致^{[10，16，22]}。此外，多种机器学习算法在丘陵区域的推理制图精度明显高于整体和平原区域，这是因为地形因素是土壤形成发育的主要影响因素之一，丘陵区域地形起伏度大，基于地形差异进一步引起其他环境变量的空间分布差异，从而导致土壤类型和性质的差异^[26]。因此，在进行土壤类型推理制图时，需要根据不同的地形条件，选择合适的环境变量来进行推理制图，以提高制图精度。

除母质外，本文选取了12个地形因子和15个遥感因子作为土壤类型推理制图的环境变量集。为了降低数据冗余，减少数据计算量，本文对环境变量集信息进行了充分挖掘，优选出各区域环境变量。在各区域筛选后的环境变量中，地形因子在整体和丘陵区域的贡献率分别为52%和48%，而在平原区域的贡献率仅为38%；遥感因子在整体、丘陵和平原区域的个数分别为2、4和10，贡献率分别为13%、20%和42%，再次证明了平原区域地形差异较小，土壤类型分布的地带性特征明显减弱或不明显，需要更多的遥感因子参与推理制图以提升精度^[27]。同时，与使用所有环境变量进行推理制图相比，使用筛选后环境变量进行推理制图时，土壤类型制图精度均得到了有效的提升，虽然提升效果在不同区域有所不同(表 3)。例如，整体区域推理制图的总体精度提升了4.96%，而平原区域推理制图精度仅提升了1.43%。此外，本研究还发现，按地形推理制图精度高于整体推理制图精度，与王苏放^[28]研究结果一致。因此，本研究表明考虑制图方法优选、环境变量筛选以及地形分区制图方式均能有效提升平原-丘陵区域土壤类型推理制图精度。

但本研究尚有一些不足之处，已有研究表明不同的样点布设方式会产生不同的推理制图效果^[29]，本文仅采用随机采样的方式布设训练点，未考虑训练点与图斑边界距离，导致图斑边缘训练点存在错分问题。因此若能够布设空间代表性高且准确分类的样点，有望进一步提高土壤类型推理制图精度。其次，随着计算机科学的发展，当前深度学习算法在推理制图中的应用逐渐受到关注，在今后的研究中可以做进一步的探索。同时，已有研究表明人类活动因子有利于提高土壤类型推理制图精度^[30]，但本文中未考虑人类活动因子，后续的研究中可以引入该类因子以进一步提升推理制图精度。最后，本文构建的各区域最优土壤-景观关系模型仅适用于本研究区，其是否具有较好的可迁移性，还有待进一步探讨。

本研究以平原-丘陵区域为整体研究区，将其划分为平原区域和丘陵区域两个地形单元，探讨了综合考虑制图方法优选、环境变量筛选以及制图方式对平原-丘陵区域土壤类型推理制图精度的影响。利用DT、RF、GBDT和XGBoost四种机器学习算法进行各区域土壤类型推理制图，各制图方法均能实现土壤类型分类，但它们在不同地形下的制图效果不同，应根据具体地形来选择适宜的制图方法。基于RF模型内嵌的因子重要性进行的最优环境变量筛选，证明了环境变量筛选不仅能有效减少数据冗余，还能提高推理制图精度，但不同地形下的提升效果存在差异，整体区域的精度提升最为显著，总体精度和Kappa系数分别提高了4.96%和0.06。遥感因子在平原区域的贡献度远高于整体和丘陵区域。此外，采用地形分区制图方式能有效缓解多尺度地形复杂性，进而提高平原-丘陵区域的土壤类型推理制图精度。依据综合考虑制图方法优选、环境变量筛选以及制图方式构建的最优土壤-景观关系模型绘制的研究区土壤类型图，发现其总体精度和Kappa系数分别为73.05%和0.69。因此，本研究为构建平原-丘陵区域适宜性土壤-景观关系模型，从而提高土壤类型推理制图精度提供了一定的理论参考和技术支撑。