数字土壤制图（Digital soil mapping）以土壤—景观模型为理论基础，利用具备协同空间变化的土壤类型和地理环境数据，借助于多种空间分析统计手段，推测土壤的空间分布，进而实现土壤类型的快速成图，因而被公认为是一种高效、实用的土壤数据更新方法^[1-2]，该方法是基于土壤成因理论和土壤-环境知识理论，核心是土壤空间差异与环境因子的协同关系^[3-5]。然而，将具有协同变化的地形因子作为环境变量参与成图推理过程时，在丘陵区域研究效果比较显著，而在环境变异较小的平原区域，难以有效地反映土壤的空间变异，探索和验证平缓地带的辅助推理因子成为目前本领域的研究瓶颈。已有研究分析和评估了遥感数据源与常规数据集组合，对研究典型土壤分布空间差异进行了探索，证明了遥感因子也可以作为平缓地区土壤推理制图的重要知识源^[6-7]。因此，获取平缓区域土壤空间变异信息的环境协变量成为数字土壤制图的研究重点。

Xu等^[8]探索TN（Total Nitrogen）与从Landsat 8、Sentinel 2和WorldView-2图像得出的各种光谱指数之间的相互关系，基于不同遥感图像的TN多空间预测模型，评估不同遥感传感器预测总氮的能力；Zhang等^[9]综合分析训练样本和降水数据的像素日期对基于时间序列MOD09A1图像的区域SOM预测精度的影响；Duan和Zhang^[10]利用Landsat8遥感影像证明了添加纹理特征能够提高土壤类型解译的准确性；在研究方法上，利用决策树、随机森林、支持向量机等方法进行数据挖掘推理土壤属性进行制图研究较多，但在对于参数较多的集成算法环节，目前使用较多的是运用自动调参的方法-网格搜索（Grid Search），但此方法只适用于较小的数据量级^[11]。此外，影响土壤属性的环境因子较多，如何在较多的因子中选取最优的适用于平原和丘陵的因子、科学有效地进行数字土壤制图也是需要解决的问题。综上所述，虽有研究将纹理特征作为参与遥感图像地表识别的指标，但很少考虑将纹理特征与地形因子重要性程度用于土壤分类。

本研究的主要目的是综合利用地形因子和遥感影像作为推理制图的知识源，分别运用递归特征消除算法、ReliefF算法和基于Tree的特征选择算法探究环境协变量在土壤推理制图中的重要性，利用经基于TPE（Tree Parzen Estimator）的贝叶斯优化（Bayesian optimization）算法参数调优后的梯度提升树（Gradient Boosting Decision Tree）模型进行建模，同时比较不同特征筛选算法后整体区域与按照地形区域推理制图精度结果，为提高平缓区域制图精度奠定方法基础。

1 材料与方法 1.1 研究区概况

研究区为湖北省麻城市乘马岗镇，地势北高南低，属于典型的平原向丘陵过渡的地貌类型，经纬度范围为114°54′—115°03′E，31°21′—31°30′N，面积约28.14km²。成土母质有“花岗岩片麻岩坡积或残积物”、“白云石花岗岩片麻岩坡积或残积物”和“近代河流冲积物”3种类型，研究区如图 1所示。

研究区属亚热带湿润季风气候，具有冬冷夏热、四季分明、降水充沛、光照充足、无霜期长、雨热同季的特征。农业集约化程度高。该区资源比较丰富，农业耕地历史悠久，具有优越的自然条件。

1.2 数据来源 1.2.1 地形数据与遥感数据

本研究以中国土壤发生分类土壤类型图、等间距为10m的等高线数据（来源于湖北省测绘局）和遥感影像（选取植被覆盖度适中的2020年12月15日的Sentinel-2遥感影像）为知识源。传统土壤类型图（1︰50000）来源于第二次全国土壤普查成果，研究区内共有9种土壤类型，如表 1所示。

1.2.2 样本数据

本文将从传统土壤图中选取训练样本，采用均匀网格与随机抽样集合的方式确定了1973个训练样点集（表 1），再利用数据挖掘算法概括样本中所蕴含的土壤学家“隐形知识”^[12-15]并对其可视化，训练样点及野外独立验证样点141个空间分布（图 2）。

为探究环境因子在按地形推理制图与整体区域制图的差异，本研究依据地形起伏度分布图和土壤志等历史资料，以海拔150m为平原和丘陵区域的临界值，得到平原和丘陵区域的训练样点集。

1.2.3 土壤-环境知识获取

利用等高线数据生成高程（Elevation），从而使用ArcGIS10.7提取坡度（Slope）、坡向（Aspect）、水平曲率（Horizontal curvature）、平面曲率（Plan curvature）、剖面曲率（Profile curvature）、地形湿度指数（Topographic wetness index，TWI）、地表粗糙度（Roughness）和地形起伏度（Relief）即8个地形因子。同时，使用ENVI5.3对预处理后的遥感影像的多个波段进行主成分分析（Principal component analysis，PCA），进而基于主成分分析提取纹理特征（Textural feature）。光谱信息的主成分（PC）波段是原始波段的线性合成，主成分分析可以隔离噪声和较少数据维数，并提供更广泛的空间信息。因此本文提取归一化植被指数（Normalized differential vegetation index，NDVI）、第一主成分（First principal component，FPC）和8个纹理特征包括均值（Mean）、方差（Variance）、相异性（Dissimilarity）、信息熵（Entropy）、对比度（Contrast）、协同性（Homogeneity）、二阶矩（Second Moment）、相关性（Correlation）为后续筛选推理制图的遥感因子。

1.3 特征筛选方法 1.3.1 递归特征消除算法

本文使用Py Charm2021中的RFE和RandomForestClassifier包，通过十折交叉验证的递归特征消除来确定所选特征的最佳数量。递归特征消除算法（Recursive feature elimination，RFE）是一种贪心算法，用于找到最优特征子集并降低数据维度以获得最高精度的模型。算法的基本思想：首先构建底层模型并对初始特征集进行训练，计算每个特征的重要性并赋予权重，然后去掉重要性低的权重最小的特征，将剩余特征组合成新的特征子集，并进行训练。该过程递归重复，直到达到所需的最终特征数目^[16]。在本研究中，RFE使用随机森林（Random forest，RF）作为底层模型，选取默认参数，并选择训练结果精度最大的特征子集作为构建模型的特征组合。RF为由多棵决策树构成的集成模型，模型的最终输出结果由森林中的每棵决策树共同决定，最终输出结果为每棵决策树输出的均值。RF算法的具体过程如下：在原始训练特征集中用Bootstrap抽样方法获得n个特征子集；对每个特征子集选择m个特征，并对每个训练特征子集构建决策树，得到n个决策树模型，建立起随机森林；计算每棵决策树的结果，将n棵决策树输出结果的均值作为最终结果。

1.3.2 ReliefF算法

Relief算法只能处理两类类别的数据。Kononenko^[17]提出的ReliefF算法是在Relief算法的基础上进行改进的，它不仅能够处理有噪声的多类及回归问题，还可以处理不完整的数据。在解决多类及回归问题时，ReliefF算法会从训练样本集中随机选取一个样本R，并从与R同类的样本集中找出R的k个近邻样本（Near hits），从与R不同类的样本集中找出k个近邻样本（Ness misses），然后更新每个特征的权重，该过程会重复m次^[18]。

ReliefF具有较高的评估效率、对数据的类型没有限制要求等优点，可以很好去除无关特征，但不能去除冗余特征，是公认的效果较好的过滤式的特征评估算法，具备过滤式算法的典型特征，即省去了对特征子集的分类器进行训练的步骤，减少了计算量，且简单高效。

1.3.3 基于Tree的特征选择

Python中SelectFromModel包是一个meta-transformer（元转换器），它可以用来处理任何带有coef_或者feature_importance属性的训练之后的评估器。如果相关的coef_或者feature_importance属性值低于预先设置的阈值，这些特征将会被认为不重要并且移除掉。除了指定数值上的阈值之外，还可以通过给定字符串来使用内置的启发式方法找到一个合格的阈值，可以使用的启发式方法有mean、median以及使用浮点数据乘以阈值。

本文加权模式使用梯度提升树来提取所有特征的重要性。

1.4 土壤推理制图

本文运用梯度提升树（Gradient Boosting Decision Tree，GBDT）模型进行推理制图。GBDT集成梯度提升（Gradient Boosting）与决策树（Decision Tree）两种算法思想，其中梯度提升属于boosting算法，其每次建立的单个学习器是在之前建立的模型损失函数梯度下降方向。损失函数（loss function）越大，模型越容易出错，如果损失函数持续下降，则意味着模型在不断地自我改进，而最好的改进方式就是损失函数在其梯度（Gradient）方向下降。基于梯度提升算法的学习器叫作梯度提升机（Gradient Boosting Machine）。理论上，梯度提升机可以选择各种不同的学习算法作为基础学习器。GBDT采用决策树作为基础学习器，通过构造一颗决策树，然后在已有模型和实际样本输出的残差上构造另一颗树来迭代训练。最终将所有决策树的结果作为输出。在每轮迭代中，GBDT会产生一个弱分类器，该分类器在上一轮分类器的残差基础上进行训练。为了避免过拟合，弱分类器应该足够简单，满足低方差和高偏差^[19-20]。在训练过程中，降低偏差可以提高分类器精度，因为模型越简单越容易产生欠拟合。为了提高数据的有效性，训练集采用随机抽取，以最大限度保证分布情况的一致性。当然，选择合适的超参数和建立有效的预测模型也是至关重要的。

贝叶斯优化算法采用的是一种逼近思想，当算法的计算过程非常复杂、需要较高的迭代次数时通常能够起到很好的效果，大多被用于机器学习算法模型的超参数取值确定。实现贝叶斯优化的库有很多种，本文选用基于TPE贝叶斯对GBDT进行优化。

对GBDT算法模型进行参数优化的基本描述是：首先设置超参数取值区间组合，接着利用贝叶斯最优化算法的相关理论，通过不断训练GBDT模型并使用评价函数对每组参数组合得到的分类预测结果进行评价，以此寻求最优的参数取值组合。最后，将最优参数取值组合代入GBDT算法模型中，提高预测准确度并降低算法运行耗时。

1.5 精度验证

本研究采用混淆矩阵以及Kappa系数对推理结果的精度进行评价。混淆矩阵是用矩阵形式来评价土壤分类精度的一种方式，用矩阵中的总体精度、用户精度以及生产精度三个评价指标来比较推理结果与实际采样之间的差异。

Kappa系数用于衡量分类的效果，计算结果通常介于0~1之间，0~0.2为极低一致性，0.2~0.4为一般一致性，0.4~0.6为中等一致性，0.6~0.8为高度一致性，0.8~1为几乎完全一致。

2 结果与讨论 2.1 特征筛选结果 2.1.1 基于递归特征消除算法的特征筛选

随着递归特征消除反复构建随机森林模型，精度达到最高值，分类效果最佳。此外，RF算法提供了总和为1的特征重要性接口，用来比较各个特征的重要性。经过参数优选和多次训练，确定RF的主要参数决策树数量为1000。

RF-RFE算法使用精度（Accuracy）值作为特征选择过程中筛选的原则。首先算法将删除重要性最小的特征，如下图 3所示，删除第20个特征，然后RF-RFE重新计算剩余特征的精度值，再次删除重要性最小的特征，依次循环，不断迭代删除特征^[21]。图 3表示的是每一次迭代时采用十折交叉验证所删除特征的精度值，迭代结束后，根据精度的大小选取最优特征集合。

从下表 2可以得出，在20个环境因子中，母质在整体区域（Whole）、平原（Plain）和丘陵区域（Hill）均是重要性第一的因子，但其重要性不同，在丘陵区域重要性最大，重要性值为0.253，平原区域重要性较低，重要性值为0.170。遥感因子中的第一主成分、均值和相关性在平原区域的重要性评分明显高于整体和丘陵区域。协同性、对比度和相异性纹理特征在整体区域、平原区域和丘陵区域重要性较低，均不参与后续建模。

由图 3可知，通过十折交叉验证，得出丘陵的最佳特征数量为7，对应精度最高为0.655，整体区域特征数量达到17时精度最高，精度为0.553，当特征数量达到14时平原区域精度为0.453，且精度最高。在每次筛选得到的精度结果中，整体区域精度和丘陵区域精度高于平原区域。表 4显示丘陵区域选取特征排名前7的因子、整体区域选取特征排名前17的因子以及平原区域选取特征排名前14的因子。

2.1.2 基于ReliefF算法的特征筛选

母质相比其他因子重要性较大，基于ReliefF算法的特征权重图中使得其他因子重要性显示不明显，因此，本文将得到去除母质的特征权重（如图 4所示），特征权重越大，表示该特征对样本的区分效果越好。本文计算所有特征权重均值作为阈值选择最优特征子集^[22]。整体区域中，特征权重的均值为0.07。高程、归一化植被指数、坡向、坡度、地形起伏度、地形湿度指数的权重超过7%，其余特征的权重低于7%。因此，权重高于均值的特征包括高程、归一化植被指数、坡向、坡度、地形起伏度、地形湿度指数共6个特征，加上母质共计7个用于构建特征子集。平原区域中，特征权重的均值为0.031。高程、坡向、归一化植被指数的权重超过8%，地形湿度指数、坡度、地形起伏度、均值的权重值在3.1%~8%之间，其余特征的权重低于3.1%。因此，权重高于均值的特征包括高程、坡向、归一化植被指数、地形湿度指数、坡度、地形起伏度、均值共7个特征，加上母质共计8个特征用于构建特征子集。丘陵区域中，特征权重的均值为0.081。高程、坡向、归一化植被指数、坡度、地形起伏度超过8%，其余特征的权重低于8.1%。因此，权重高于均值的特征包括高程、坡向、归一化植被指数、坡度、地形起伏度共5个特征，加上母质共计6个特征用于构建特征子集。

综上，根据筛选结果表明，高程、坡向、归一化植被指数、坡度、地形起伏度为公共特征因子，遥感因子中的纹理特征均值只出现在平原区域中，可知均值在平原区域中重要性程度大于整体和丘陵区域，在土壤类型推理制图中均值在平原区域占有一定的比重。

2.1.3 基于Tree的特征选择结果

利用Py Charm2021软件中GradientBoostingClassifier和SelectFromModel包，对特征因子进行重要性评分，结果如下表 3所示，通过给定内置阈值分别去除不重要因子，得出整体、平原和丘陵区域最优特征数量分别为11、15和15。

依据整体区域重要性将环境因子从高到低进行排序，结果如表 3所示，在整体区域中被筛选出的遥感因子占所有因子的5%，而在平原和丘陵区域中被筛选出的遥感因子所占比重均为25%，但其因子不同，平原区域重要性较高的遥感因子为归一化植被指数、第一主成分、相关性、均值和信息熵，而丘陵区域重要性较高的遥感因子为归一化植被指数、第一主成分、相关性、均值和方差，纹理特征中协同性、对比度、二阶矩和相异性均为不重要的因子。

2.2 特征筛选算法分析

本文基于递归特征消除算法、ReliefF算法以及基于Tree的特征选择构建整体区域、平原以及丘陵区域的特征子集。

综上，由3种特征筛选算法筛选结果可以看出在平原区域纹理特征的重要性程度高于整体与丘陵区域，并且遥感因子中的归一化植被指数和均值相较于其他的遥感因子所占比重较大，原因是平原区域中采样点之前的地形因子差异比较小，平原区域地势起伏较小，难以有效地反映土壤的空间变异，所以需要遥感因子参与其中进而对平原区域土壤类型分类进行推理。

2.3 GBDT参数优化

GBDT作为一种集成算法，不同参数对算法结果的影响力不同。运用Py Charm2021中的hyperopt包，基于TPE对GBDT模型对如下参数进行优化。

以整体区域为例：

在进行n_estimators参数调优时，首先取值25~200，步长为25，得到最佳平均精度对应的值为125。再通过逐步逼近法，取n_estimators为75~150，步长为5，再次得到最佳平均精度所对应的值为80；进行多次尝试，将范围缩小，步长为1，得到最优的值为71。在进行learning_rate参数调优时，首先取值0.05~2.05，步长为0.05，得到最佳平均测试分数为0.35；再通过逐步逼近法，取learning_rate为0.05~0.65，步长为0.001，得到最优的值为0.096。在进行Subsample参数调优时，首先取值0.1~0.8，步长为0.1，得到最佳平均测试分数为0.6；再通过逐步逼近法，取Subsample为0.4~0.8，步长为0.05。进行多次尝试，将范围缩小，步长为0.0025，得到最优的值为0.782；进行max_depth调参，首先取值2~30，步长为2，得到最佳平均测试分数为28；再通过逐步逼近法，取max_depth：14~32，步长为1，得到最优的值（18）。min_impurity_decrease在进行参数调优时，此参数较不稳定，所以不再进行参数调整。

如果其中有的参数趋于稳定，将不进行继续调参，在目标函数中，直接写入固定值。通过实验，将基于递归特征消除算法、ReliefF和基于Tree的特征选择的筛选结果对整体、平原及丘陵区域的GBDT算法的参数设置为当前场景下的最佳参数取值。

2.4 模型制图及精度验证

采用整体区域以及面向不同地形区域分别进行土壤类型预测，基于3种特征筛选算法的整体区域以及按照地形的推理结果如图 5所示。对比基于3种特征筛选算法的整体区域以及分区域的推理图，可以发现整体上呈现出一致趋势，相比原始土壤图，整体区域制图以及按照地形推理图图斑更加破碎，沙泥土、潮土田、沙泥田以及林地泥沙土分布范围基本一致，沙泥土、沙泥田和林地泥沙土均由花岗岩片麻岩坡积物或残积物母质发育而来，沙泥土、沙泥田和林地泥沙土主要分布在东部地区，相比整体区域推理图，林地泥沙土按照地形推理效果更好。潮土田唯一由近代河流冲积物发育而来，较其他8种土壤类型而言推理效果最佳。硅沙泥田在平原区域分布，并与硅沙土相互夹杂，《麻城市土壤志》记载，这两种土壤类型的成土环境相似，土壤具有相似特征，因此这两种土壤往往相伴而生，虽如此，按照地形推理制图分类情况有所提升。局部细节中，相比于整体区域推理的土壤图，硅麻骨土和林地薄层硅麻骨土在按照地形的土壤图中图斑明显增多，空间细节更加丰富，由此可见土壤类型的空间分布与景观特征相吻合，这也进一步验证了土壤类型空间分布的准确性。

比较3种不同特征筛选算法的推理图，按照地形推理结果较整体区域推理结果更加破碎，利用野外实地采集的独立验证点对推理图进行精度评估，用于精度检验的野外实地样点共141个，对比基于递归特征消除算法、ReliefF算法以及基于Tree的特征选择算法的整体区域以及按照地形区域推理图，按照地形区域的推理图正确分类的样点个数分别为107个、88个和100个，总体精度分别为75.89%、62.41%和70.92%，分别高于整体区域推理图（64.54%、60.28%和60.99%）；而Kappa系数分别为0.727、0.565和0.666，分别高于整体制图的0.137、0.022和0.116（整体区域分别为0.589、0.543和0.55），能较好地反映研究区的土壤分布情况。

利用按照地形推理的原则建立混淆矩阵，以生产精度和用户精度为指标验证按照地形推理的土壤图精度。整体推理图和按照地形的推理图的生产精度和用户精度统计结果如下表 4所示。生产精度是该土壤类型正确分类点个数与该土壤类型验证点个数之比，用户精度是正确分类点个数与被判定为该土壤类型的点数之比。土种级别的分类制图相较于土属级别，分类精度的提升从原理而言更加困难。从整体区域制图和按照地形推理制图结果得出，基于递归特征筛选算法的精度分别高于基于ReliefF算法和基于Tree的特征选择算法精度的4.28%和3.55%，按照地形推理中，在平原区域，由于基于ReliefF算法筛选出平原的因子较少，加入的遥感因子较少，导致总体精度低于基于递归特征消除算法以及基于Tree特征选择算法的精度。基于3种特征筛选算法推理制图中，按照地形推理制图精度高于整体制图精度。基于递归特征消除算法的制图中，9种土壤类型生产精度和用户均高于整体区域制图，在基于ReliefF算法中，潮土田、硅沙泥田、硅麻骨土和林地薄层硅麻骨土的生产精度和用户精度高于整体推理图，基于Tree的特征筛选算法中，按照地形推理的沙泥土、硅沙泥土、林地泥沙田、硅沙土和林地薄层硅麻骨土的生产精度和用户精度均高于整体制图精度。

本研究提出的基于递归特征消除算法、ReliefF算法以及基于Tree的特征筛选算法3种不同的特征挖掘算法对环境变量（包含地形因子和遥感因子）进行筛选，对比整体区域和按照地形的推理制图，分析遥感因子在推理过程中的权重，以选择稳定的遥感指标，探索进一步提高制图准确性的途径，取得了较好的效果，但仍存在不足需要进一步探讨。在小区域尺度下，能够反映土壤与环境之间关系环境因子还有很多，缺乏其他方面的因子参与筛选，比如：土地利用方式。因此，选择更多的环境因子进行科学的分析筛选并运用至推理制图中，使得因子在判别土壤类型更具准确性是今后重点研究方向之一。如今，新的机器学习方法层出不穷，选择更多更好的机器学习模型进行对比推理制图，使结论更具客观性和准确性也是今后研究者重点研究的课题之一。本文未采用SoLIM等模糊推理算法，因此不能进行推理结果的不确定性研究，但这应该是本研究今后工作的重要方向之一。

3 结论

本研究提出了一种基于特征筛选算法进行土壤类型推理的数字土壤制图获取方法，利用递归特征消除算法、ReliefF算法以及基于Tree的特征选择方法对环境协变量进行特征选择并分析，并且采用基于TPE贝叶斯优化的超参数优化后的GBDT进行建模，利用野外采样点进行独立验证。结果表明：环境因子在参与平原和丘陵区域的推理制图时具有不同的重要性，遥感因子在平原区域参与程度大于丘陵区域，其中遥感因子中的NDVI和Mean相较于其他遥感因子对土壤类型的分类重要性更大。3种筛选算法对环境因子进行重要性评分，从而判断因子对土壤类型的影响程度，能够减少大量的冗余，提高分类精度及效率，为进一步环境变量与土壤类型之间的联系建立提供了新的研究思路。基于3种筛选算法按照地形推理的制图精度均高于整体推理制图精度，递归特征算法的按照地形推理制图精度最高为75.89%，分别高于ReliefF算法和基于Tree的特征筛选算法。因此，本研究提出的不同特征筛选方法以及基于TPE贝叶斯优化对GBDT进行参数优化的方法可为土壤制图等相关研究提供理论参考。