土壤普查是了解土壤资源空间分布及合理利用土壤资源的前提，是农业生产的基础。传统的土壤普查主要是土壤学家以大比例尺航片和地形图为底图，通过实地采样获得土壤-景观模型，以此判断不同土壤类型的空间分布并手绘成图[1-2]。此方法存在的精度低、人为主观性强及时间和人工成本高[3]等问题推动了数字土壤制图的诞生与发展。与传统土壤普查不同的是，数字土壤制图方法是以空间分析和数学方法为技术手段，用栅格形式来表达土壤空间分布信息，此种方法可更好地表达出土壤分布的渐变特征[4-6]。土壤-环境知识的获取和土壤空间推理方法的选择是当前精细数字土壤制图的关键。

目前，数字土壤推理制图过程常用到的协同因子多为由数字高程模型(digital elevation model，DEM)衍生的地形因子。遥感影像因其时效性强、数据分辨率高、易获取及含有丰富的下垫面信息等优势逐渐应用于土壤制图，其包含的地表光谱信息可用于反演地物属性、地物分类，提取的多种遥感指数可用于定量或定性地表述植被生长状况。遥感影像可以多种形式参与到土壤推理制图过程中，如利用土壤有机质的遥感响应波段以实现其含量的估算[7]，借助影像建立特定时段的地表动态反馈模型[8-10]，在时序遥感影像数据上提取植被的轮作模式[11-12]，或直接在遥感影像数据中提取地物的分类特征[13-15]等。但在以往的研究中，土壤制图使用的遥感影像多为单时相数据[16]，不同时间点遥感影像产生的光谱信息具有较大差异，制图结果也存在一定差异。且单一时间点的遥感数据易受天气、季节等因素的影响，在遥感影像时期的选择上也存在着人为主观性，这些环节都将直接影响制图结果。

土壤分类常用的数据挖掘算法主要有模糊聚类[17-18]、决策树[19-20]、支持向量机[21]、判别分析[22]、随机森林[23-24]等。其中，随机森林算法可以有效避免原始数据的缺失及噪声、异常值造成的精度低等问题，在分类和回归等多方面具有一定的优势，在国内数字土壤制图领域中该方法被广泛应用于获取土壤类型[23]或属性[25-27]的空间分布信息。本研究拟将多时相遥感影像数据和随机森林算法应用于土壤分类，为土壤分类提供新的技术与方法。

本研究使用图斑面积加权法在传统土壤图上设计采样点，在母质和地形等环境协同因子的基础上，加入多时相遥感指数，通过随机森林算法建立土壤类型与环境协同因子组合的对应关系，并推理得到研究区各土壤类型的空间分布信息，进而对传统土壤图进行更新，同时探究多时相遥感影像数据在土壤制图中的可行性。

研究区域位于华家河镇滠水河流域(图 1)，隶属于湖北省黄冈市红安县(114°23′~114°49′E，30°56′~31°35′N)，地处大别山南麓，境内地势北高南低，为低山丘陵区。土地利用类型以林地和农耕地为主。根据全国第二次土壤普查结果，研究区土壤类型共12种，母质类型共6种，其中石英片岩、基性岩分别与硅沙泥田、乌沙泥田相对应，因此这两种土壤类型不参与后续建模和推理过程。为方便后续随机森林建模和预测，数据预处理时将石灰岩、片麻岩、近代河流冲积物、泥质岩4种母质分别编号为1、2、3、4。研究区各土壤类型和其对应的母质类型如表 1所示：

本研究使用的原始数据主要包括母质类型图、等高线、遥感影像和传统土壤图。根据红安县土壤图和地质图及当地土壤普查专家的意见，研究区域可以划分成6种母质类型。地形因子数据集是在湖北省测绘局的10 m等高距的等高线数据的基础上衍生而来，在ArcGIS中内插得到研究区的数字高程模型，提取常用的高程、坡度、坡向等地形因子。遥感影像数据来源于哨兵二号卫星(https://scihub.copernicus.eu/)，分辨率为10 m，选择一个完整的作物生长周期，成像时间由2016年11月至2017年10月，每月选择一景遥感影像，借助ENVI软件平台计算归一化植被指数，具体成像时间如表 2所示。

研究区原始土壤类型图的比例尺为1:50 000，由第二次全国土壤普查结果经数字化后得到。传统土壤图是当地土壤普查专家进行野外调查的制图成果，其往往蕴含着土壤学家的“隐性知识”[3]，即土壤-环境知识，利用决策树、随机森林等数据挖掘算法可有效提取隐性知识并加以利用，更新传统土壤图。

采样点的设计方式对数字土壤制图的结果具有重要影响。本研究中采样点数据依据传统土壤图图斑面积，采用面积加权法进行布设采样点，并保证每个图斑中最少采样点数量不低于10个，最后共生成采样点数量6 686个。

随机森林(random forest，RF)分类算法是Breiman[28]于2001年提出的一种机器学习方法，是一种以决策树为基分类器的Bagging集成算法，广泛应用于回归和分类。随机森林算法基于自助采样法(bootstrap sampling)[29]，即有放回的采样或重复采样，从初始数据集中有放回的反复抽取数据作为训练样本集，在采样过程中约有36.8% 的数据始终没有被抽中，这部分数据为袋外(out-of-bag，OOB)数据，用以对模型性能进行估算，称为袋外估计。OOB误差是随机森林模型泛化误差的无偏估计[28]。随机森林模型中有两个重要参数：节点分裂次数mtry和决策树数量ntree，这两个参数值的选择直接决定了随机森林模型的性能。

本研究利用横截面采样、主观采样和均匀采样3种方式得到的264个验证点对推理制图结果的准确性进行评价。验证点集包含研究区存在的所有土壤类型，且在整个研究区均有分布，如图 2所示。建立实地验证点与推理土壤图的混淆矩阵，可得到总体分类精度、生产精度、用户精度及Kappa系数等评价指标。其中生产精度是预测土壤图某土壤类型正确分类的点数与野外采样点中该土壤类型的总点数的比值，而用户精度是预测土壤图某土壤类型正确分类的点数与预测土壤图中该土壤类型的总点数之比[30]。

本研究提取的原始环境因子如表 3所示，包括母质类型、地形因子和遥感因子。基于采样点数据分别建立土壤类型与地形因子、遥感因子之间的关系(图 3和图 4)。从图 3可以看出，除高程、坡度和坡向这3个地形因子外，其余因子在不同土壤类型之间差异较小，因此本研究中仅选取高程、坡度和坡向参与到后续随机森林建模和推理过程中。从图 4可以看到，不同土壤类型的全年植被指数变化规律相似，夏季7月和8月NDVI达到峰值，冬季则最低；林地细沙泥土(21)、林地沙泥土(31)和林地棕色石灰土(51)这3种土壤类型的归一化植被指数全年均处于较高水平，这可能与其土地利用类型有关；此外，4—6月及9—10月，各土壤类型对应的NDVI值差异较大，可能是因为不同类型的植被生长及秋季叶片凋谢的速度不一致，从而导致植被指数出现较大差异。

本研究通过R语言中random forest包来建立随机森林模型，环境因子为输入变量，土壤类型为输出变量，采用遍历的方式，以OOB误差的大小来确定模型的最优参数组合。由此得到mtry和ntree分别取6和600时，随机森林模型性能最佳，进而利用该模型进行预测。

对比原始土壤图(图 5A)和随机森林预测推理土壤图(图 5B)，整体上两者的各类土壤空间分布具有一定相似性，林地沙泥土、沙泥土、细沙泥土、沙泥田、棕色石灰土等在两个土壤图中的空间分布基本一致，但棕色石灰土、浅潮沙田、沙泥土、沙泥田的面积有所减少。推理土壤图的图斑数量显著增多，展现出更精细的土壤空间细节信息，且不同土壤类型之间不再存在明显的“分界线”，过渡平滑，呈现出自然界中土壤空间分布的渐变特征。

为进一步评价土壤分类的准确性，我们采用野外独立样本点对推理后土壤图进行评价。评价结果以混淆矩阵的形式展示，其中包括总体分类精度、各土壤类型的生产精度和用户精度及Kappa系数等精度评价指标，如表 4所示。总体来说，推理结果图的总体分类精度达到了86%，分类结果较好；Kappa系数为83%，表现出高度一致性，这说明验证样点与土壤图之间的契合程度高，推理土壤图所表达的土壤类型空间分布信息是可信的。大部分土壤类型的分类精度较高：除浅潮沙田外，其余土壤类型均有良好的用户精度(≥70%)；除棕色石灰土和细沙泥田外，其余土壤类型的生产精度良好(≥70%)。因此，基于多时相植被指数的土壤推理方法具有较好的制图结果。

从混淆矩阵中可以看到，首先，林地沙泥土、林地棕色石灰土的生产精度和用户精度均较高。这是由于这两种土壤类型主要分布于研究区北部高程较大的丘陵处，土地利用类型为林地，与同母质发育而来的其他土壤类型相比，其归一化植被指数值较大，容易区分；沙泥田和潮沙田与研究区内河流分布关系密切，土壤养分较丰富，肥力水平高，宜耕性强，作物长势稳健，亦具有较高的植被指数，故分类精度较高。其次，浅潮沙田的生产精度和棕色石灰土的用户精度均较低，混淆矩阵中可以看出，这两种土壤类型容易被错分为同种母质下的潮沙田和林地棕色石灰土，其原因可能是同一母质发育而来的两种土壤类型分布区域间地形起伏程度差异较小，且土地利用类型相近，作物生长规律相似，在时序遥感影像上表现为数值相近、无明显差异性的归一化植被指数，因此利用地形因子和多时相的遥感植被指数均不能有效区分，从而产生分类误差。此外，细沙泥田的生产精度较低，分类时易与细沙泥土混淆，从而影响到整体分类精度。

综上所述，经264个野外独立样点验证结果表明，基于多时序遥感影像推理得到的土壤图分类精度达到了86%，土壤图斑更为破碎，数量变多，空间详细程度更高，Kappa系数为83%，验证样点与推理之间具有显著的一致性水平，因此利用这种方法来推理得到研究区的土壤类型空间分布信息是可行的。

本研究提出了一种基于多时相遥感影像和随机森林算法的土壤推理制图方法。利用时序哨兵二号遥感影像数据提取归一化植被指数，作为遥感因子与常用的高程等地形因子、母质类型组成完整的环境因子集，并通过随机森林算法提取土壤-环境知识后推理制图，经野外独立验证点验证其结果分类精度高达86%，更新后土壤图表达出的研究区内各土壤类型的分布信息可信度较高。

已有研究大多采用单时相遥感影像数据，通过提取植被指数、纹理信息、主成分等遥感光谱指数用于土壤分类及土壤属性制图，并取得了较好的结果，但同时我们也应看到不同季节影像的制图结果具有较大的差异，且在影像时间的选择上多存在主观性。相比于单时相遥感影像数据易受天气、季节等因素的影响，本研究采用多时序遥感影像数据参与土壤制图，避免了依赖单时相遥感影像数据的局限性。周紫燕等[23]在同一区域的研究中，同样采用了随机森林算法进行数据挖掘，但不同的是其使用的环境因子中除了常用的地形因子外，只加入了单时相的遥感因子，制图精度为76%。与其相比，本研究加入的是多时相遥感因子，推理土壤图的准确性有了明显提升，精度提高了10%。本研究结果具有较高的预测精度，各土壤类型空间分布的详细程度有较大提升，可以满足未来对高精度数字土壤制图的要求，为土壤普查工作提供参考依据。

本研究虽取得了较好的分类精度，但仍有一些需要改进的地方。在遥感影像的数据基础上，仅提取了时序归一化植被指数，未对其他可从影像提取的光谱信息加以利用，而归一化植被指数实际上由植被和土壤两部分组成，土壤背景对归一化植被指数干扰强烈，影响较大，且在植被稀疏的地方，归一化植被指数的表征能力较差。未来可引入更多能间接反映土壤发生发育的遥感光谱指数，同时，为避免环境因子较多带来的数据冗余现象，可采用主成分分析对因子进行降维处理，以减少不必要的信息冗余；此外，可对研究区进行地貌分区，针对不同区域采用不同的推理协同因子组合和权重，分区制图。