基于多变量与RF算法的耕地土壤有机碳空间预测研究*——以福建亚热带复杂地貌区为例

袁玉琦，陈瀚阅，张黎明，任必武，邢世和，童珺玥

（福建农林大学资源与环境学院，土壤生态系统健康与调控福建省高校重点实验室，福州 350002）

耕地土壤有机碳（Soil Organic Carbon，SOC）是衡量土壤质量和表征土壤肥力的重要指标，也是农业温室气体减排潜力的主要来源。掌握土壤有机碳的含量及空间分布，对提升耕地质量、评估土壤健康和缓解全球气候变化具有重要意义。实地采样和分析是获取土壤属性空间分布格局的传统方法，但因费时费力、成本高、区域通达性差限制了其在地区、国家乃至全球尺度的应用[1]。随着 3S技术的发展，越来越多学者通过建立辅助环境变量与土壤属性之间的预测模型来实现区域土壤属性空间制图[2]。

目前，用于数字土壤制图的预测模型主要包括逐步线性回归模型[3]、递归偏最小二乘算法[4]及随机森林算法（Random Forest，RF）[5]等。其中RF算法是一种机器学习算法，它由分类和回归树（Classification and Regression Tree，CART）模型构成，因其在处理多元非线性数据方面的优势，越来越多地被用于土壤属性空间预测[6]。相继有学者将RF模型用于大区域有机碳含量及储量[7-9]的空间预测，并取得了较好的预测精度（r值在0.7～0.95范围内）。但RF算法在土壤有机碳空间制图的应用起步较晚，其在不同区域的预测精度和可行性方面仍需进行深入研究，其中需要关注的一个研究重点即是环境变量的筛选[10-11]。

在基于 RF算法的土壤有机碳空间预测中，最为常用的辅助环境因子包括数字高程模型（Digital Elevation Models，DEM）和遥感植被指数[7，11]等遥感变量。土壤有机碳含量高低与地表植被覆盖状况和地形地貌密切相关[12]，而遥感变量因在地形和植被空间信息表达上具备独特优势，且较易获取，常被用于数字土壤制图。但土壤理化性质是多方面成土因素综合作用的结果，单纯利用遥感变量预测土壤有机碳有所欠缺。部分学者[8，13]在此基础上增加了相对较易获取的气候因子（如温度、降水）和土壤类型、成土母质等土壤定性属性。然而，与土壤有机碳机理上互为相关的土壤定量属性数据因获取难度较大，以往较少被用于土壤有机碳空间预测，但随着测定技术的发展，该类变量获取难度有所降低，日渐被部分学者关注。Liu等[14]在植被要素和气候因子等较易获取的环境变量基础上增加了容重、土层厚度等土壤属性，基于RF算法模拟中国北方草地土壤有机碳密度（R2达到0.73）。Were等[9]在对比三种机器学习算法（包括RF算法）用于南非地区SOC储量预测的研究中，创新性地选择Mg、Ca、P、全氮（Total Nitrogen，TN）和pH等土壤定量属性联合遥感数据和地形及衍生因子作为模型驱动因子，结果显示土壤属性中的 TN变量在所有三种模型的变量贡献率中占绝对主导地位。

由此可见，基于RF模型预测SOC所用的环境变量随科学测定技术的发展有所创新，但国内外研究少有基于不同组合的环境变量对预测精度影响的评估，尤其是获取难度较大的土壤定量属性，作为变量加入对模型精度的贡献如何尚不十分明确。此外，针对国内基于RF算法的SOC预测多应用于地貌类型较单一的小尺度区域[13，15]，而大范围复杂地貌区因地形、气候、植被、土壤肥力状况等因素在水平和垂直方向上存在显著分异[16]，导致SOC的空间分布规律及其与环境变量间的关系更为复杂，不同变量对SOC预测精度的贡献率仍需进一步评价。RF模型正因为在模拟变量间复杂非线性关系和高阶相互作用方面的突出优势[17]，更适用于本文亚热带复杂地貌区土壤有机碳含量的空间预测。

目前，我国已经完成两次全国性的土壤普查、农业农村部测土配方施肥重大国家计划，也正在进行农业长期定位实验和第三次土地调查等，积累了大量的土壤样点及属性数据，这些较难获取的土壤定量属性数据可为提高模型预测精度提供极为便利的研究条件。基于此，本文以典型的亚热带复杂地貌区——闽东南地区为研究区域，采用随机森林模型，选取易获取的遥感变量和气候因子及不易获取的土壤属性作为模型输入，基于大量土壤有机碳实测样点分别训练两种不同变量组合数据集驱动的RF模型，并进行有机碳含量空间分布格局的预测和对比，以期为区域地形地貌复杂区有机碳储量的准确估算提供理论基础。

1 材料与方法

1.1 研究区概况

研究区位于福建省东南部，简称闽东南地区（23°32′～26°04′ N，116°53′～119°91′ E），由漳州市、厦门市、莆田市、泉州市（除永春县、德化县和金门县）和福州市（仅长乐区、福清市和平潭综合试验区）组成（图1），区域总面积为29 373 km2，其中耕地面积约4 384 km2。该区位于闽中大山带戴云山-博平岭段东南侧，地貌复杂、地形起伏，海拔最高1 152 m，最低0 m，从内陆向沿海递减，其中山地丘陵（海拔>200 m）面积高达44.33%。受亚热带季风气候与地形因素综合作用，年平均气温为 20.0℃，年平均降雨量为 1 526 mm，主要集中于 3—8月，占全年降雨量80%以上。从土地利用方式来看，可以用作农业用途的土地占比小（耕地仅为14.93%），且受海拔、地形等要素影响，80.15%的耕地位于福建省最主要的农业高产区——沿海平原台地（海拔≤200 m）[18]。该区土层深厚、肥力高，农业生态气候条件好，耕地利用强度大，土壤属性空间变异强烈。全区耕作田块分布较为零散，耕地面积小而破碎，土壤类型以水稻土为主，占耕地总面积的 70.97%，其次为赤红壤，占耕地总面积的19.91%，而红壤、滨海盐土、潮土、风砂土、黄壤、石灰土和紫色土的分布面积较小，合计面积仅占总面积10%左右。总之，地貌的特殊性与复杂性导致气候、植被、土壤等空间分异显著，造成不同地区、不同海拔高程的耕地自然条件差异明显[16]。

1.2 数据来源

本研究采用的土壤有机碳实测数据来源于国家农业农村部2017年末测土配方施肥调查样点数据，共计1 257个，包括调查样点地理坐标及土壤有机质（Soil Organic Matter，SOM）、有效磷（AP）、速效钾（AK）、水解性氮（HN）、交换性镁（Ex-Mg）和 pH等土壤理化属性（图1）。每个样点均按照密度控制、代表性、均匀性、优先性及适当性原则设置，采样深度为 0～20 cm。在所选地块均匀随机采集10～15个耕层土样，充分混合后采用四分法留取1.5 kg土样装袋，自然风干过筛后备用。土壤有机质含量采用重铬酸钾氧化—外加热测定，有效磷采用碳酸氢钠浸提—钼锑抗比色法测定，速效钾采用乙酸铵提取—火焰光度计法测定，水解性氮采用碱解扩散法测定，交换性镁采用乙酸铵浸提—原子吸收分光光度法测定，pH采用酸度计法测定。

选取样点中测定的土壤理化属性，联合遥感变量和气候变量作为SOC空间预测的环境变量，具体构成见表1。遥感变量中的植被指数和地形因子分别基于Landsat8 OLI影像和ASTER GDEM高程影像提取得到，数据来源于美国地质勘探局（United States Geological Survey，https://www.usgs.gov/）和地理空间数据云网站（Geospatial Data Cloud，http://www.gscloud.cn/）。影像空间分辨率均为30 m，投影坐标系为WGS_1984_UTM_Zone_50N。为使植被指数能反映采样时地表真实状况，影像选取月份与实际调查样点获取时间一致，且天气晴朗，基本无云层覆盖。气象因子基于2017年福建全省22个标准气象站点的地面气候资料日值数据集插值而来（具体站点位置见图1），数据来源于国家气象中心网站（National Meteorological Information Center，http://data.cma.cn/）。

1.3 环境变量的获取、组合与筛选

为了探索多源协同变量解释 SOC空间变异的可能性及土壤属性的加入对提高模型预测精度的贡献程度，本研究将表1中16个环境变量根据获取的难易程度分成两种不同的组合：（1）仅基于易获取的遥感变量（Remote Sensing Variables，RS）和气象因子（Climate Factor，CF）训练的模型Training model based on simple-to-obtain variables（RF-S）；（2）基于遥感变量、气象因子和土壤属性（Soil Attribute，SA）所有变量训练的预测模型 Training model based on all variables（RF-A）。

表1 环境变量的构成Table 1 The composition of environmental variables

影响土壤有机碳含量的环境变量众多，模型训练前需利用 RF算法预测所产生的袋外误差的大小对部分变量进行剔除[10]，即依据逐次剔除某一变量后RF模型袋外得分（Out-of-bag Score，OOB Score）的增减判断该变量是否保留，OOB Score值增加则变量剔除，反之保留[11]。

1.4 RF模型构建和验证

以各模型筛选后的变量数据集为输入，利用RF模型进行有机碳含量的回归预测。RF模型是建立在决策树基础上的一种集成学习方法，通过多次bootstrap抽样获取多个随机样本，并通过这些样本子集分别构建相应决策树，从而构建随机森林[10]。当模型用于回归预测时，取所有决策树预测结果的均值作为最终的预测结果[19]。模型的运算过程中需设定两个关键参数：n_estimators和max_depth，其中 n_estimators为决策树的数量，即使用 bootstrap重抽样的次数，max_depth为决策树的最大深度。依据预测过程中产生的OOB Score的大小，RF-S和RF-A模型设定的（n_estimators，max_depth）分别为（1 400，9）和（1 400，12）。

为衡量两种不同变量组合下模型的表现，将所有训练样本点作为模型验证点进行验证。此外，使用不同抽样百分比（20%、30%、40%、50%、60%和70%）的验证数据集独立评估模型性能（例如，80%样点训练，余下20%样点进行验证，以此类推），并与常用的地统计插值方法——普通克里格插值模型（Ordinary Kriging，OK）结果进行对比。使用平均绝对误差（MAE）、均方根误差（RMSE）、相关系数（r）和变异系数（CV）评估模型预测的绝对误差表现，使用相对误差（Relative Error，RE）和相对均方根误差（Relative RMSE，RRMSE）定量化模型的准确性程度，数值越小，模型准确性越佳。

1.5 数据处理方法

数据处理主要为环境变量提取、RF模型构建和验证，以及耕地土壤有机碳空间分布图的生成3个方面，具体如下：

环境变量提取。土壤定量属性采用ArcGIS10.2地统计插值模块中普通克里格法分别进行插值以获取 30 m× 30 m栅格数据。遥感变量中的 NDVI和TVI植被指数通过Landsat8 OLI影像波段运算得到，计算公式分别如下：

式中，近红外（Near Infrared，NIR）和红光（Red）波段反射率为Landsat8 OLI影像利用ENVI5.3软件进行大气校正、镶嵌、裁剪等预处理后获取的研究区地表反射率。地形因子中的DEM由福建省23幅ASTER GDEM影像数据经坐标转换、镶嵌及研究区腌膜运算提取。而其他4个地形因子（坡度、坡向等）是利用ArcGIS10.2软件Spatial Analyst模块，基于 DEM 计算得到。气象因子栅格数据由全省气候资料数据集使用ArcGIS10.2地统计插值模块中反距离权重法（Inverse Distance Weighted，IDW）插值后，再通过研究区腌膜提取得到，空间分辨率与遥感变量一致，为30 m。将上述环境变量形成的栅格数据集，利用 ArcMap提取到与土壤有机碳实测样点空间匹配的训练数据集中（如图1所示），用于RF模型的构建和验证。

RF模型构建和预测的实现均通过 Python scikit-learn库中RandomForestRegressor包实现。变量相对重要性排序可直接调用工具包中 feature_importances属性实现。

基于 RF模型和 OK模型生成耕地土壤有机碳空间分布图，用于评价不同模型在SOC空间异质性表达上的优劣。针对RF模型，将空间分辨率为30 m的遥感变量、气候因子栅格和土壤属性栅格依据对应的变量组合分别输入 RF-S与 RF-A模型，得到SOC空间分布格局。OK模型则是基于所有SOC样本点，使用普通克里格方法插值后重采样为30 m获取。利用ArcGIS10.2制图模块完成SOC空间分布专题制图。

2 结 果

2.1 同变量组合下RF模型预测精度对比

如表2所示，基于全部样本点，RF-A和RF-S模型精度较好，均呈现高度相关（r>0.8）和中等变异水平（CV>23%），相对误差RE小于10%。但与RF-S模型相比，加入了土壤属性变量的RF-A模型预测误差有显著下降，RMSE和 MAE分别下降45.13%和42.68%，表明加入土壤属性变量（N、P）有利于提升模型拟合度及预测精度。

表2 两种不同变量组合下RF模型的预测精度Table 2 Prediction accuracies of RF models under two different combinations of variables

如表3所示，两种RF模型预测值的均值与SOC实测值的均值非常接近，约等于 14 g·kg–1，但预测结果范围明显被压缩，变异系数和标准差值均变小。由图2也可发现，两种模型在SOC低值区（累积百分比40%以下）略高估实测值，而在SOC高值区（累积百分比75%以上）明显低估实测值，但高值样点数（SOC>26 g·kg–1）所占比重较小，所以两种 RF模型的预测结果能够解释SOC大部分空间变异。总体而言，加入了土壤属性变量的 RF-A模型预测值的累积分布图更接近实际结果，能更好地表征区域SOC的动态变化范围。

表3 两种不同变量组合下RF模型的预测结果与SOC实测值对比Table 3 Comparison SOC measured value with prediction results of RF models under two different combinations of variables

2.2 RF模型环境变量重要性

各模型最终筛选出的用于土壤有机碳预测的变量如表4所示，各类别均有环境变量被保留参与模型构建。由表4可知，在RF-S模型中，重要性最高的环境变量为DEM（23.22%），气候因子分列重要性排序第二到五位，虽然单个变量重要性稍弱，但累积贡献率达58.20%，依旧占据主导地位。在加入土壤属性变量的RF-A模型中，N的重要性超过DEM，位列第一，说明这两个变量是影响闽东南地区耕地SOC空间变异的主要协同因子，且N和SOC的关系更为密切。值得注意的是，遥感变量中的植被指数NDVI和TVI，全部通过变量筛选被保留参与模型预测，尽管单因素贡献较弱，对空间结果的预测也是不可或缺的。各类环境变量中贡献率最高的因子分别为地形因子DEM、植被指数NDVI、气象因子Mint和土壤属性N。

表4 RF模型特征变量重要性排序Table 4 Ranking of relative importance of environmental variables of RF model

2.3 基于不同抽样百分比的精度检验

在不同梯度抽样百分比下，对两种 RF模型以及OK模型的训练和预测精度进行检验，结果如表5所示。在训练数据集中，不同抽样百分比下，RF-S与RF-A模型预测的SOC值与实测值的相关系数r分别在0.75～0.85和0.89～0.92之间，而 OK模型r处于0.60左右，个别情况下出现低值0.55，整体精度明显低于RF模型。由误差MAE、RMSE、RE和变异系数CV值亦可发现，RF-A模型在各梯度下计算误差最小，变异程度最大，OK模型表现最弱。显然，RF-A模型拟合精度最高，RF-S次之，而OK模型最低。三种模型的拟合能力随训练样本量呈现不同的变化趋势：RF-A模型受样本数量影响小，整体表现稳定；RF-S模型小幅波动，规律性不强，在训练样本为70%时精度最低，其余情况较稳定；而OK模型波动较大，在训练样本为 80%时表现最佳（r=0.65），30%时表现最差（r=0.55）。由于多种环境变量的协同作用，RF模型鲁棒性较好，而OK模型对采样点要求较高，且插值结果随取样空间尺度增大会产生明显的平滑效应，这与陈飞香等[20]使用克里格法对土壤原始样点在不同采样密度下的插值结果相一致。所以，当采样点较少时，应选择 RF模型，并优先选含有土壤属性变量的 RF-A模型。在验证数据集中，相较于其他两个模型，RF-S整体精度不高，RF-A模型计算的r除在20%和40%验证数据集较OK模型低0.1，其他情况下均略高于OK模型。总体而言，RF-A模型预测精度普遍优于OK模型，表现最好。

2.4 耕地土壤有机碳含量空间分布

如图3a、图3b、图3c所示，RF-S、RF-A和OK模型预测结果的空间分布均呈现为东部沿海较低、西部内陆较高。虽然三种模型的总体趋势比较相似，但在研究区南部和中部偏北地区，OK模型与RF-S模型预测的SOC高值区明显小于RF-A模型。在上述预测差异较大区域随机选取子区域（图3a～图3c中a、b框）放大显示（图3a～图3c右边的子图a、b），可发现RF-A模型的SOC含量分级区间数明显多于其他两个模型，且空间变异更强。总体来看，RF-A模型无论在模型精度或空间异质性表达上均为最优模型，以下仅对最优模型预测的SOC空间分布格局进行分析。

RF-A模型反演得到闽东南区SOC均值为14.70±2.95 g·kg–1，范围为 3.63～25.51 g·kg–1，其中 13～19 g·kg–1区间的面积占比最高，超过研究区耕地总面积的65%，主要分布在西部内陆闽中大山带戴云山-博平岭段东南侧；小于 10 g·kg–1和大于 19 g·kg–1的面积占比较低，不足10%，分别分布在闽东南地区三大平原（漳州平原、泉州平原、莆仙平原）和西部海拔最高地；10～13 g·kg–1区间所占面积在19%左右，位于高低值过渡带。

通过数值分析发现，SOC空间特征与各类环境变量中贡献率最高的四个因子（DEM、NDVI、Mint、N）呈现明显的相关性（如图4所示）。针对 DEM和Mint两个因子，一般海拔越高，相应温度越低，而SOC含量正是与海拔高度呈正相关（图4a），而与年最低温度Mint呈负相关（图4c）。这与杨顺华等[21]的观点一致，认为平原丘陵过渡带土壤有机碳与高程等稳定因素呈极显著正相关。NDVI值越大，植被覆盖度越高，SOC含量越高（图4b），这是由于生物量也是土壤有机碳最重要的来源。值得注意的是，SOC含量在一定范围内随着N含量的升高而增加，但当 SOC含量均值>17 g·kg–1时，对应N含量呈下降趋势（图4d），可能原因是土壤碳氮比（C/N）是一个相对固定的数值，土壤N素含量极大地影响SOC含量，然而在高海拔地区，虽氮素平均含量有所减少，但鉴于高海拔地区低温降水少的条件下有机质分解速率较低，SOC含量并无下降。

3 讨 论

3.1 闽东南地区土壤有机碳空间预测及主要环境变量影响

基于 RF-A模型预测得到的闽东南地区耕地SOC含量东部沿海较低、西部内陆较高，这与刘素真[22]利用全国第二次土壤普查数据和近期野外采样数据开展的福建省有机碳含量模拟结果基本一致。西部地区因海拔高度较高伴随温度较低，用作农业生产的可能性降低，植被逐渐转化为枯落物较多的自然植被，同时土壤微生物分解有机质的速度减慢，矿化作用减弱，从而导致有机质的积累量逐渐增加[10]，而东部低海拔区域，温度较高，农业生产较便利，频繁的耕作加剧土壤扰动，促进有机碳分解，使土壤有机碳周转速率加快，积累量减少，导致SOC含量较低，这与上述模型计算结果相一致。由此可见，RF-A模型预测的研究区SOC含量空间分布格局是合理的，从这方面证明了该制图方法的可行性。

特征变量贡献率分析显示，水解性氮（N）是预测土壤有机碳含量最重要的环境变量，这与预期一致。大量研究已证实陆地生态系统碳氮循环存在紧密耦合[23]，碱解氮与土壤有机碳存在显著正相关关系[24]，但因数据获取不易，以往较少用于土壤碳空间预测。Were等[9]和谢恩泽等[25]在RF模型研究中也证明全氮TN是解释SOC空间变异最重要的变量，这与本文研究结论基本一致。

遥感变量中的DEM的重要性仅次于N。但与N不同，DEM是通过地形地貌差异间接影响SOC的空间分布，属于外部因素。与一般区域耕地分布于相对平缓且交通便利地带不同，本研究区位于亚热带复杂地貌区，耕地分布的海拔差异大，不同高度主要分布的土壤类型也有所不同，这就强调了DEM对SOC预测结果的影响。齐雁冰等[10]在利用RF模型反演陕西省土壤有机质的研究中得出了相似的结果，即在地形地貌复杂且耕地面积较小的区域内，DEM对SOC空间预测的贡献率相对较高。此外，高程和地形可通过影响降水情况间接影响 SOC含量的空间分布。当来自太平洋的暖湿气流进入福建沿海并向西北运行过程中首先遇到闽中大山带，由于地形对气流的抬升作用，在东南坡产生较多的降水，这与SOC的空间分布高度一致。由此可见，高程和降水、温度对SOC空间分布的影响存在重叠，并且高程的影响更大，所以降低了气候因子在模型中的重要性。

在基于遥感影像获取土壤有机碳的研究中，植被指数是最常用的变量之一，而本研究中所选取的两个植被指数在 RF模型中的重要性并不突出，但二者均通过 OOB Score的筛选被保留参与模型预测。从数值分析上看，SOC含量大致与NDVI呈正相关关系（图4b），这与以往研究相一致[26]。

3.2 RF-A模型精度

研究结果可得，RF-A模型无论在模型精度或空间异质性表达均为闽东南区 SOC空间预测的最佳选择。基于全部1 257个样点训练的RF-A模型的r=0.95，RRMSE值与RE值分别为13.89%和–5.91%（表2）。Hengl等[27]提出，RRMSE≤40%以内为可以接受的模型准确性，当 RE≤±10%时，模型模拟结果处于可接受范围。可见，RF-A模型精确度高，且模型的准确性处于可接受范围。相较于 RF-S模型，加入土壤属性因子的RF-A模型RMSE和MAE分别下降45.13%和42.68%，可见，RF模型的预测精度受目标变量与辅助环境变量之间的相关性强弱控制，当加入与SOC紧密耦合的内部影响因子（如N），RF模型的精度显著提升。谢恩泽等[25]在对比保留或移除辅助因子 TN对苏南农田土壤有机质空间分布预测精度的影响时，得到了一致的结论，进一步证实了土壤因子用于亚热带复杂地貌区有机碳空间预测的可行性。

但RF-A模型在验证数据集中的r值（表5）与以往文献[10-11]相比，没有明显的优势。可能原因在于研究区为复杂地貌区，耕地面积小而破碎，SOC预测难度较一般地势平坦、地形地貌单一的小区域大，且本文采用的大样本数据也会显著增加验证结果的误差。总体而言，RF-A模型在全部样点和不同抽样百分比的预测精度均处于准确度范围内（不同梯度下RRMSE值均小于34%，见表5），且预测结果与实际相符，可用于研究区SOC空间分布格局的获取。

4 结 论

本文采用随机森林模型，选择易获取的遥感变量和气候因子及不易获取的土壤属性作为模型输入，基于大量土壤有机碳实测样点分别训练两种不同变量组合数据集驱动的 RF模型，对比验证加入与未加入土壤属性变量构建的 RF模型精度，并与普通克里格插值模型进行对比。结果显示，联合遥感变量、土壤属性和气候因子共同构建 RF-A模型精度最高，可作为预测该研究区SOC含量的高效方法，在模型驱动因子中加入N、P等土壤属性变量，能显著改进模型拟合度与精度，更好地提升区域SOC空间异质性的预测能力。相较于普通克里格插值模型，RF-A综合模型在不同抽样百分比的预测精度普遍更优，且稳定性更好，因此在采样点较少的情况下，应优先选择 RF-A模型。变量贡献率结果表明，影响研究区耕地SOC预测结果空间分布最重要的环境变量是与SOC直接相关的土壤属性——水解性氮，其次是DEM。