掌握准确的土壤类型空间分布，能为土地规划利用、农业生产生活等提供更可靠的参考依据。数字土壤制图（Digital soil mapping）通过土壤属性与地理环境的协同变化关系，利用多种数学方法和空间分析方法进行区域内土壤属性的推理预测^[1]，并凭借其高效性、低成本、高精度等特点逐渐得到广泛应用^[2]。由其推理过程可知，选取合适的环境因子、使用恰当的推理方法是数字土壤制图研究过程中的重点内容。

环境因子选取的有效性直接决定着数字土壤制图的精确度^[3]。有较多的研究在进行环境因子选取时，直接将母质作为一种预测因子用于推理制图^[4-5]。但母质作为土壤发育形成的物质基础和植物矿物元素的最初来源，会直接影响成土的过程和速度^[6-7]，当制图点的样本量较少时直接进行环境因子筛选很有可能会降低母质在所有因子中的重要性。黄魏等^[8]按照母质类型使用5个地形因子对研究区环境因子进行规则挖掘以实现推理制图，然而该研究更侧重于利用少量地形因子进行传统土壤图的更新，同时也并未深入探究母质类型对环境因子选择的具体影响。

机器学习常用于各领域的回归和分类任务^[9-10]，在数字土壤制图中，应用较多的机器学习方法如决策树，随机森林，支持向量机，人工神经网络等^[11-13]，支持向量机算法与KNN算法在土壤类型分类预测的理论基础上具有一定的相似性^[14]。这些方法均通过环境因子与土壤类型的关系根据算法自身特点进行推理预测，在很大程度上忽略了环境因子本身存在的关联关系。本文选用了随机森林方法（Random Forest，RF）、基于C5.0决策树的土壤景观推理模型方法（Soil-land Inference Model，SoLIM）、K邻近算法（K-Nearest Neighbor，KNN）作为推理预测方法，并在此基础上使用FP-Growth（Frequent Pattern Growth，FP-Growth）算法挖掘环境因子关联关系，通过本研究设计的结合流程将环境因子关联关系分别与三种制图方法结合，并验证方法的有效性。总之，本研究在利用实地采样点对研究区土种类型进行推理的框架下，不仅探究了母质类别对环境因子选择的影响，同时也提出了新的制图方法将环境因子关联关系融入到土壤类型推理中，为数字土壤制图提供了新的制图思路。

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

研究区位于湖北省麻城市乘马岗镇，总面积约为28.14 km²，研究区中心地理坐标为31°25′30″N，114°58′30″E，高程介于80~346 m之间。乘马岗镇是一个山区农业乡镇，属于亚热带湿润季风气候，全镇地势北高南低，研究区共有3种母质类型，9种土壤类型，研究区的基础信息如图 1所示。

1.2 数据来源 1.2.1 母质及样本信息

研究区由花岗岩-片麻岩坡积或残积物、白云石英片岩-花岗片麻岩坡积或残积物以及近代河流冲积物等3种母质构成，母质信息来源于湖北省第二次土壤普查数据。本研究的样本数据通过完全随机采样、主观采样、混淆区域边界采样和均匀采样（网格采样）等4种方式相结合共采集实地样点141个，采样时间为2020年10月—2021年6月，根据分层随机抽样方法将其分为100个制图点和41个验证点^[15]，研究区各土壤类型对应的母质及样点个数信息如表 1所示。

1.2.2 环境因子

在研究区内，近代河流冲积物母质上仅有潮土田一种土壤类型，且面积占比较小，故在分区建模时将此区域与白云石英片岩-花岗片麻岩坡积或残积物母质区域合并研究。地形因子和遥感因子作为预测变量进行土壤类型推理被证明是可行且精度高的^[16]。本研究使用了地形和遥感共计18个因子，地形数据来源于湖北省测绘局的等高线数据，利用ArcGIS10.8软件生成10 m分辨率的数字高程模型（DEM），继而得到高程（Elevation）、坡度（Slope）、坡向（Aspect）、平面曲率（Plan curvature）、剖面曲率（Profile curvature）、地形湿度指数（Topographic wetness index，TWI）、地形起伏度（Relief Degree of Land Surface）、地表粗糙度（Surface Roughness）等8个地形因子。遥感数据通过欧空局（ESA）官网获取2021年12月10日的Sentinel-2A遥感影像（云量较低、地表植被覆盖度适中），通过ENVI5.3软件对遥感影像处理得到归一化植被指数（Normalized differential vegetation index，NDVI）、第一主成分（First principal component，FPC）以及均值（Mean）、方差（Variance）、协同性（Homogeneity）、对比度（Contrast）、相关性（Correlation）、相异性（Dissimilarity）、信息熵（Entropy）、二阶矩（Second Moment）8个纹理信息等共计10个遥感因子，获取的18个环境因子用于后续因子筛选。

1.3 研究方法 1.3.1 基础方法

基于随机森林的递归特征消除算法被证明较适用于土壤类型推理，而且该方法可以评估因子的重要性^[17]。本研究采用此方法按照母质类型分别筛选最合适的环境因子用于制图，并与未按母质类型筛选环境因子的制图结果进行比较。

随机森林算法、基于C5.0决策树的SoLIM方法、KNN算法均能对土壤类型推理预测。随机森林算法在土壤类型或属性预测方面均得到了广泛应用并展现出了优越性能^[18-19]，这得益于其两次的随机选择并建立多个决策树（随机选取样本点和随机选择特征子集），增加了模型的多样性^[20]，但是随机选择就会产生误差，如何既保留算法的多样性，又降低其随机性带来的误差，成为了随机森林进一步发展的瓶颈；SoLIM方法通过隶属度函数得到各环境因子的隶属度，再根据限制因子学原理确定环境因子组合相似度，继而利用模糊最大算子确定空间内各点的土壤相似度，将相似度硬化，最终得到该区域的土壤图^[6]，该方法在对像元土壤类型判定时，某一影响程度最大的环境因子能够更多地决定该像元的土壤类型；KNN算法简单地利用预测点与已知类别点之间的距离远近进行土壤类型的推理^[21]。因此，对于三种方法而言，若能在其推理预测中加入每种土壤类型对应的环境因子内部关联关系，即可以更好地弥补由于制图推理原理本身形成的不足。

1.3.2 结合方法

频繁项集即频繁出现在数据集中的集合，具体是指在数据集中出现的频率不低于用户指定阈值的集合/项集，能够表示挖掘数据中的相关性、关联关系等^[22]。在进行频繁项集挖掘前，参考关于低频关联信息挖掘的分级方式^[23]，并结合分级结果挖掘的频繁项集最终确定将研究区的各环境因子使用自然间断点分级法分为5级，每个级别对应一个编号，由于因子过多，仅展示部分分级结果，如表 2所示。将环境因子数据分级是挖掘环境因子关联信息的基础，例如挖掘结果显示某种土壤类型挖掘的频繁项集为（d2，s1）的含义为：该土壤类型在研究区常分布在高程为133.94~167.6，坡度为0~3.97的范围内，（d2，s1）即为该土壤类型的典型环境因子组合，可以用来反映环境因子的关联关系。

对于本研究在数据分级时使用的是整体区域像元点的18个环境协变量，而非样本点对应的协变量数据用于分类统计做如下解释：（1）由于预测制图的对象是研究区整体，样本点数据并不能全部覆盖各种环境因子数值上下限，且样点数量较少时，统计分类结果可能未必精确。（2）有研究在进行土壤样本代表性修正时，将样本点的环境因子朝着与研究区整体环境因子更相近的方向进行样本代表性修正^[24]，表明整体区域的环境因子更能体现土壤的空间变化。

FP-Growth算法是目前常用的关联关系挖掘算法之一，采用树的结构遍历两次完成频繁项集挖掘^[25]。在该算法中，某个频繁项集关联关系的强弱主要通过支持度和置信度来衡量，支持度表示该频繁项集在所有历史交易数据中同时出现的次数，置信度表示在某个频繁项集出现的情况下，另一个项出现的概率^[21]，即关联关系。为了进一步获取更加准确的频繁项集，本研究参考邱小倩等^[23]关于衡量耕地质量数据库低频关联关系的指标公式，提出关于土壤类型分类频繁项集高频出现的过滤公式，称为误分概率P。误分概率越大，表示该种频繁项集越不容易将这种土壤类型误分给其他类，具体公式表达如下：

$ \text{P=}\frac{频繁项集在一种土壤类型中出现的次数}{此频繁项集在所有土壤类型中出现的次数} $

(1)

本研究将FP-Growth算法挖掘到的频繁项集，经误分概率筛选和排序后与文中使用的机器学习制图方法相结合，再次推理制图。数据的预处理流程如下：获取各种制图方法的推理结果，将制图结果的所有像元转为矢量点，继而提取研究区环境因子分级结果至点，且与每个像元点的预测土壤类型一一对应，形成制图结果环境因子分级数据集，以和利用采样点挖掘的每种土壤类型的频繁项集数据进行匹配判断。

在结合方法中，每个像元最终土壤类型确定的核心即判断制图点筛选后的每种土壤类型的频繁项集分级结果（a）是否是全局每个像元点的环境因子分级结果（b）的子集，具体过程如下：（1）设某个像元的土壤类型为A，读取该像元对应的b，检索土壤类型A挖掘的a，判断土壤类型A的所有a是否是该像元对应b的子集，若是，则判定此像元的最终土壤类型为A；若不是，则依次判断其他土壤类型的a是否是该像元对应b的子集，若是，则将a对应的土壤类型赋予该像元，完成该像元最终的土壤类型判断。（2）若遍历了所有土壤类型的a，均未找到某个a是该像元对应b的子集，则将此像元的土壤类型判定为原类型A，完成该像元的土壤类型判断。依次对全局的所有像元根据此流程进行土壤类型的再次判定，即完成频繁项集与某种制图方法相结合的过程。

1.3.3 精度验证

本研究基于41个验证点使用混淆矩阵以及Kappa系数对制图结果进行精度评价，各土壤类型精度通过混淆矩阵的用户精度和生产精度评价，整体制图效果使用总体精度以及Kappa系数衡量。Kappa系数通常在0~1之间，当Kappa系数在0~0.2，0.2~0.4，0.4~0.6，0.6~0.8，0.8~1时分别表示预测结果与实际结果一致性程度较差、一般、中等、较强和很强^[26]。

2 结果与讨论 2.1 环境因子筛选

高程、坡度、平面曲率、剖面曲率、地形湿度指数等5个指标在小区域内可以较好地反映土壤的发生和形成^[27-28]，本研究将这5个因子作为确定性环境因子直接参与制图，不再进行筛选。环境因子经递归特征消除算法计算后的模型精度如图 3所示。

由图 2可知全局、花岗岩-片麻岩坡积或残积物母质区域（代号1）、白云石英片岩-花岗片麻岩坡积或残积物和近代河流冲积物母质区域（代号4，代号10）的环境因子组合精度最高时对应的环境因子个数分别为6个、4个、2个，结合直接使用的5个环境因子，最终分别确定三种方式筛选和使用的具体环境因子结果见表 3。由图 2可知，母质花岗岩-片麻岩坡积或残积物区域的模型验证精度普遍高于母质为白云石英片岩-花岗片麻岩坡积或残积物和近代河流冲积物母质区域以及全局筛选时的模型精度，白云石英片岩-花岗片麻岩坡积或残积物和近代河流冲积物母质区域的模型精度整体上基本超过全局因子筛选的精度，表明将母质分开制图可能获得更好的效果。

2.2 频繁项集挖掘

本研究使用FP-Growth算法挖掘100个制图点的环境因子分级数据，确定FP-Growth算法挖掘的最小支持度为3，置信度为1，即表示在所有样点中至少有3个及以上的样点同时拥有这种频繁项集组合。地理学第三定律提出：两点的地理环境越相似，土壤属性越相似^[29]。于本研究而言，频繁项集出现的频次越高，越能代表某土壤类型的典型环境条件，即越可以用来推理土壤类型。在这一原则下，确定了频繁项集的置信度为1，即此环境因子组合一定会一起出现，为了尽可能使每种土壤类型均被挖掘到可以使用的频繁项集，确定了支持度为3。

由于挖掘的频繁项集数据较多，此处以某个频繁项集及其对应的误分概率为例进行阐述，如林地泥沙土挖掘的某个频繁项集为（m2，r3，p2），且误分概率为1，意味着此土壤类型在高程范围值为133.94~167.6，剖面曲率为-0.39~0.37，第一主成分为–208.02~267.74的环境因子组合下不会被误分为其他土壤类型。本研究获取所有土壤类型的频繁项集共计141个，经误分概率公式计算后统计，统计结果如图 3所示。

由图 3可知，部分土壤类型挖掘到了相对较多的频繁项集（如林地薄层硅麻骨土），除硅麻骨土外，其余几种土壤类型均挖掘到误分概率 > 0.75的频繁项集。因此本研究在不断试验并结合图 3统计结果的基础上最终确定将误分概率 > 0.75的频繁项集作为制图分类的依据（硅麻骨土使用误分概率最高的频繁项集），以获得相对准确的推理结果。

2.3 调参过程

本研究在环境因子筛选时使用了递归特征消除算法，制图时分别使用了随机森林，基于C5.0决策树的SoLIM和KNN方法，其中基于C5.0决策树的SoLIM方法是通过Clementine软件挖掘规则集，并结合SoLIM软件实现土壤类型的推理，具体过程参考了韩浩武等^[30]的研究。其余机器学习模型的调参过程均参考张晓婷等^[17]的研究，各种算法调用及参数调节均依托于PyCharm2019软件编写python代码实现。

2.4 制图结果 2.4.1 基于所有母质直接制图结果

本研究在利用全部母质直接制图时，经因子筛选、规则集建立等对整体研究区域制图。在将母质作为因子之一直接用于制图时使用12种环境因子分别通过随机森林算法、SoLIM方法、KNN算法进行研究区土壤类型的推理预测，推理预测的结果如图 4所示。

由图 4可知，三种制图方法对于某些土壤类型的预测在空间分布上是相似的，但是仍有部分土壤类型预测差异明显，其中KNN算法预测结果各土壤类型交错分布现象最为严重。随机森林和SoLIM方法预测有明显差异的为硅麻骨土，据该区域原始土壤图记载硅麻骨土在研究区内的分布面积相对较小（湖北省第二次土壤普查数据），即SoLIM方法预测此土壤类型有较高的错误率。三种制图方法的结果具有一个共同点，即在某一母质区域内均出现了不属于该母质发育而来的土壤类型，这也表明将母质分开制图是有必要的。

2.4.2 基于母质类别制图结果

按照母质类别分别筛选环境因子并制图的推理结果如图 5所示。三种制图方法中SoLIM方法预测结果更加集中连片，这也与该方法的推理原理有一定关系，但是很少预测出硅麻骨土，对于硅砂土的预测也较为零散。随机森林对各土壤类型的推理相对均衡。KNN算法预测结果中，硅沙泥田并未被集中连片的预测，硅砂泥土被预测出了较大区域，硅麻骨土和林地泥沙土均被预测出较少面积。三种制图方法在母质花岗岩-片麻岩坡积或残积物区域的推理精度明显要好于其他母质区域的推理结果，这也与图 2因子筛选时得到的模型精度相符；在母质白云石英片岩-花岗片麻岩坡积或残积物和近代河流冲积物区域，随机森林和SoLIM方法的预测结果更为相似，在研究区的中部均预测出了硅沙泥田和硅砂泥土，三种方法在东部和西部分别推理出了林地薄层硅麻骨土。总之，三种制图方法在某些土壤类型推理上得到较为相似的分布预测，种类较未按母质类型制图时多，且各方法的推理结果更加集中连片，更能反映土壤的空间分布特征。同时，按照母质类型制图很大程度地避免了在某种母质区域内推理出不属于该母质发育而来的土壤类型的问题。所以，在使用实地采样点进行多母质区域的数字土壤制图时，应考虑母质类别进而选择最合适的环境因子进行推理预测。

根据母质类别分别筛选环境因子进行制图是本研究的探索方向之一，但该方法也有一定的弊端。这是因为在相邻土壤的边界处往往呈现一定程度的过渡带，过渡带兼具两种土壤类型的属性^[6]，因此，在母质交界处，土壤类型的过渡带可能不能被很好的推理。但在利用100个采样点进行土壤类型推理时，将母质作为一种预测因子直接参与推理会降低母质在所有数据集中的重要性，这也可以解释图 4中的某一母质区域内预测出不属于该母质发育而来的土壤类型的现象（且并不处于两种母质交界处），可能当采样点数量足够多时，数据集会更好地反映母质在其中的重要性，进而减少这种现象的出现。

2.4.3 基于结合型方法制图结果

本研究通过设计结合频繁项集的流程对研究区土壤类型再次推理，制图结果如图 6所示。随机森林对各土壤类型均预测出一定面积并显示更多细节。SoLIM方法预测硅沙土面积较小，硅砂泥土面积较大，结合频繁项集后，硅麻骨土预测面积增加，硅砂泥土预测面积减少。KNN算法预测结果中硅麻骨土和硅砂泥土的变化最为明显，结合频繁项集后，推理图预测出更多林地泥沙土和硅沙泥田，其空间分布更接近其他两种算法的预测结果。结合后的三种方法与未结合频繁项集的制图结果相比在整体分布上并未有明显差距，具体变化体现在三种方法均将硅麻骨土预测出了一定面积，且精度最差的KNN方法的部分土壤类型空间分布在朝着与另外两种方法制图结果更相似的方向推理。总之，三种制图方法结合频繁项集后对整体区域推理均有一定的改善。由三种制图模型的推理原理可知，其在确定土壤类型时对环境因子内部关联关系考虑较少，本研究在推理制图时融入关联关系，从原理上降低了预测误差。

2.5 精度评价

由图 4可知，在使用实地采样点推理制图时，三种制图结果均出现了在某种母质区域内推理出并非该母质发育而来的土壤类型的现象，并且KNN算法和SoLIM方法的预测结果并不理想，因此不再分析其具体土壤类型的预测精度，仅对按母质类别制图以及结合方法制图的土壤类型进行精度分析，各土壤类型用户精度和生产精度如图 7所示。

由图 7可知，硅麻骨土的推理结果并不理想，三种制图方法均未较好地预测该土壤类型，这可能是因为该土壤类型在研究区内不仅面积小，而且据《麻城市土壤志》（内部资料1983）记载该土壤类型的地理环境描述与其他土壤类型也较为相似。频繁项集结合前后推理精度最高的依然为潮土田（100%），这主要是因为研究区内近代河流冲积物母质区域仅有此种土壤类型，不易误分。与频繁项集结合后的其他土壤类型各种推理方法的用户精度与生产精度累积长度均表现为不变或增加，这表明了将环境因子关联关系融入到制图当中的可靠性，验证了本研究制图方法的有效性。各种方法制图总效果通过总精度和Kappa系数衡量，具体结果见表 4。

表 4显示，三种方法按母质类型的制图精度均较未按母质类型制图的精度和Kappa系数高，结合了频繁项集的制图方法均较未结合的总精度和Kappa系数高。是否按母质类型制图对SoLIM方法的影响最大，这可能是因为C5.0决策树在进行规则挖掘时因子较多，样本点较少，挖掘的规则并不特别适合SoLIM方法，且SoLIM方法在进行硬化时^[31]，影响程度最大的某个因子很大程度上决定了土壤类型的判定，结合图 4b可知，硅麻骨土的因子规则很大程度影响了其他土壤类型的因子规则。然而，将母质分开制图时，使用的因子减少，且避免了某些因子之间的相互影响。因此，在进行少量样点（统计意义）的土壤制图时需考虑母质类型对环境因子的影响。

在所有方法中，结合频繁项集的随机森林制图精度最高，为70.73%。最低的为母质一起制图时KNN算法的制图结果，为43.9%。由结合频繁项集的精度可知，结合前后验证点的制图精度均为：随机森林 > SoLIM > KNN，结合频繁项集后，KNN算法总精度和Kappa系数均提升最多，分别提升了9.76%和11.70%。这主要是因为该算法的原始制图精度较低，故其提升空间相对较大。随机森林提升的效果最低，一方面是因为其本身精度较高，另一方面是因为其在随机选取特征子集时建立了多个决策树，在一定程度上也可以弥补这种误差。同时，本研究设计的结合方法是在原始制图结果的基础上进一步判断的，所以在一定程度上依赖于原制图方法，但从总精度和Kappa系数来看，本研究提出的结合制图方法是有效且实用的。

对于本研究提出的频繁项集与三种制图模型结合的方法而言，结合后三种方法的制图精度均有提升，这也是符合方法推理逻辑的，但是是否适合数字土壤制图中的其他方法仍需不断验证。并且，本研究某些土壤类型并未挖掘到频繁项集误分概率为100%的组合，可能增加相应土壤类型的采样点获取更多挖掘信息，即可以获取分类精度更高的频繁项集，从而更好地推理土壤类型的空间分布。

3 结论

本研究通过递归特征消除算法按母质类别筛选环境因子，探究了不同制图方法下母质类别对环境因子筛选的影响。继而将频繁项集分别与不同制图结果结合，验证了本研究提出方法的有效性。随机森林与频繁项集结合推理的精度最高（70.73%），推理结果能较好地描述研究区土壤类型的大致空间分布，可以为农业生产、土地保护等提供参考。

按母质类型制图时，不同方法的制图效果和预测精度均较将母质一起直接制图时好，制图结果更加集中连片，土壤类型空间分布更为合理。精度评价结果显示，不同结合型制图方法的精度均高于未结合时方法的精度，其中提升最多的是KNN方法，最少的为随机森林方法，验证了本研究提出方法的有效性以及使用FP-Growth算法挖掘典型环境因子组合的可靠性。本研究不仅为基于实地采样点进行多母质区域土种类型的环境因子筛选提供了参考，同时也为数字土壤制图提供了新的制图思路和方法，有利于更好地服务农业数字化和信息化。