基于分层多元复合模型的广东省农田土壤有机碳高精度制图

任向宁，王璐，林赋英，陈淑莹，胡月明，3，4，5

（1.华南农业大学资源环境学院，广州 510642；2.广东省土地利用与整治重点实验室，广州 510642；3.广东省土地信息工程技术研究中心，广州 510642；4.自然资源部建设用地再开发重点实验室，广州 510642；5.青海-广东自然资源监测与评价重点实验室，西宁 810016）

高精度土壤有机碳（Soil organic carbon，SOC）制图是准确测度土壤固碳能力和土壤碳储量的基础，对研究SOC时空变异规律具有重要意义。目前，受样本数量与制图方法制约，高精度SOC制图研究主要集中在空间小尺度上。而对区域碳管理、碳交易等宏观需求更迫切的国家级、省级中大尺度SOC制图精度总体不高[1]。为此，中大尺度高精度SOC 制图方法的研究显得至关重要。

农田土壤作为自然土壤的衍生物，其发生发育过程必然符合地理第一定律。同时在长期农耕过程下，农田土壤发生发育过程中受到更为频繁的人工干预与扰动，导致其强空间变异度[2]。因此，农田SOC 高精度制图难度远远大于自然土壤。当前，国内外对于土壤有机碳制图的方法主要分为4 类：第一类为基于土壤属性的土壤类型连接法[3]；第二类方法为地统计学法[4]；第三类为多元线性回归方法、空间地理回归等相关关系统计法[5]；第四类方法为人工神经网络模型、分类回归树、随机森林、支持向量机等机器学习算法模型[6]。其中土壤类型连接法适用于空间尺度较大且土壤类型比较单一的区域[7]，但对土壤类型的多样性反映不足[8]。地统计学法与相关关系统计法则较为依赖空间关联性[9-10]，在情景复杂、局部变异较强的地区，其空间特征表达能力和预测精度较差[9-13]。机器学习算法模型的预测精度对空间尺度和采样密度的依赖性较强，大尺度复杂情景下测算误差较大[14]。因此，有关学者针对大尺度复杂农田土壤环境，作出了非插值算法的探索，如构建基于因子自组织映射聚类的土壤有机质与相关因子之间的非线性关系[15]，通过Cubist 模型引入了相对坡位、高程等变量的贡献率等，以此避免辅助因素增加后插值测算模型复杂化和结果偏离度的增加，降低插值噪声强度等[16]。基于非插值算法的思路，通过测定区域SOC 空间变异的变量结构，协同不同变量分层进行克里格插值，依据变量贡献力大小通过构建分层多元复合模型[17]编制SOC 空间分布图。该方法突破了一般插值预测方法对协同变量的数量限制，能够充分体现大尺度环境中多维复杂变量的影响，提升了拟合精度。同时，通过引入变量贡献力与分层协同克里格插值，降低了复杂变量相互干扰导致的插值噪声，具有更好的制图精度和空间表达效果。

本研究以广东省为研究区，协同区域综合特征分区及大密度样本数据集，采用地理探测器确定农田SOC的变量结构，构建分层多元复合模型（Hierarchical multivariate composite model，MCM），编制研究区农田SOC高精度分布图，为理清复杂情景下中大尺度农田SOC空间分布，探究提升区域农田固碳潜力提供了科学基础。

1 材料与方法

1.1 研究区域

研究区位于中国大陆南部广东省，濒临中国南海，地处20°13′～25°31′N、109°39′～117°19′E 之间，北回归线横穿研究区，陆地面积约1.8×105km2。研究区地势总体为北高南低，最高海拔1 902 m，最低为-52 m[18]。地貌类型复杂多样，包括山地、丘陵、台地和平原等，其中山地和丘陵占研究区总面积的50%以上。研究区属南亚热带季风气候，常年平均气温21.9 ℃，日照时间1 755 h，年均降雨量1 790 mm。受地势变化与距海远近影响，研究区北部山区和南部临海平原区的年平均气温相差4.5 ℃（19.6～24.1 ℃），日照时间相差1 108 h（1 316～2 424 h），降水量相差2 084 mm（1 099～3 183 mm）。研究区天然水系发育旺盛，河渠纵横，水网稠密。2018年研究区地表水资源量1.885×1011m3，集雨面积在100 km2以上的干支流543 条，河流总长2.86×104km。根据第二次土壤普查结果，研究区陆地农田土壤包括水稻土、赤红壤、潮土等16 个土类，212 个土种。在复杂气候、地形、水文、植被、成土母质及人工活动耦合作用下，研究区土壤类型的空间分布具有显著的空间地带性和复合区域性，由北向南水平地带分布红壤、赤红壤、砖红壤的铁铝土土纲系列，以及平原地区、红岩盆地和石灰岩山地3 种复合区域土壤组合类型。2018 年末研究区常住人口1.14 亿，其中珠三角核心区、沿海经济带（东西两翼）和北部生态发展区分别占研究区总人口的55.53%、29.60%和14.87%。研究区主要包括北部山地农业、中部南亚热带农业、西南部热带农业3 个农业空间带，农田面积2.259 ×105hm2，其中水田占61.89%，水浇地占28.24%，旱地占9.87%。

研究区具有显著的南北气候差异，复杂的地貌、水文、土壤等地理环境背景及不同人口压力下的农田利用方式及扰动频率，是复杂情景下农田SOC高精度制图极具代表性的研究区域。

1.2 数据来源

本研究数据来源主要包括区域地理环境数据和农田实地采样数据，其中区域地理环境数据中土地利用/覆被数据、河流沟渠及行政区划界线主要根据研究区2018 年土地利用现状变更调查成果分离、解析取得。DEM 数据采用ASTER GDEM 数据（国家基础数据中心），用于研究区海拔高度、地形坡度、坡向等信息提取，空间分辨率30 m。农田样点土壤有机质含量（Soil organic matter,SOM）、土壤类型、表土质地、土壤容重、砾石含量、土壤pH等土壤理化属性及当地降水量等数据则采集于2010—2018 年广东省耕地地力评价调查与耕地质量年度更新调查数据集，共采集208 768 个样点数据。样点布设密度按照大田每6.67×105m2（1 000 亩）、蔬菜地每3.33×105m2（500亩）～6.67×105m2（1 000 亩）布设一个，并根据地貌类型、土壤类型、利用方式及质量等级等条件适当增加样点密度，采样深度约20 cm。

经分析，剔除离群无效样点265 个，根据筛选出的样点构建研究区208 503个有效农田土壤有机碳高精度制图样本集（表1和图1）。从样本集中随机抽取70%的样本（145 952 个，不含最高值和最低值）作为空间插值集，剩余的30%（62 551 个，空间不连续）作为验证集，验证制图精度。

图1 研究区农田土壤有机碳样点分布图Figure 1 Distribution map of field survey sample points of SOC in the study area

表1 研究区农田土壤有机质含量样本统计表Table 1 Statistical table of SOM samples of typical paddy fields in the study area

1.3 研究方法

1.3.1 土壤有机碳的测度

农田土壤有机碳包括土壤耕层有机碳含量（Soil organic carbon content,SOCC）和土壤耕层有机碳密度（SOCD），其中SOCC 在国际上普遍采用58%作为SOM 中碳含量的转换系数，因此本研究中采用公式（1）测算SOCC。样本数据集中样本采样深度基本为20 cm，故统一选取耕层厚度为20 cm，本研究中采用公式（2）测算SOCD。

式中：SOCD20为农田耕层（20 cm）土壤有机碳密度，kg·m-2；SOCC 为土壤有机碳含量，g·kg-1；BD 为土壤容重，g·cm-3；H为耕层土壤厚度，本研究取20 cm；g为＞2 mm石砾含量，%。

1.3.2 区域综合特征分区方法

由于人类长期频繁扰动与自然地理要素的耦合作用，形成了研究区农田土壤复杂的系统环境。综合区划通过耦合自然、经济、社会系统要素，叠加上区域划分方法[19]，有利于收敛研究区内农田SOC 空间差异，为构建稳定的研究模型奠定基础。目前，空间分区方法主要分为定性、定量和GIS 空间模型三大类。其中定性方法有特尔斐法、主导因素法等，其精度主要取决于专家知识与经验判断，指标及其重要性的判断离散度较大。定量方法有空间叠置法、矩阵判断法[20]、状态空间法[21]、聚类分析法[22]等，其对分区影响要素的空间表达细微、准确，但分区结果的空间连续性较差。GIS 空间模型法则融合了各分区要素的空间统计、分析与处理，具有优秀的空间数据处理和分析能力，逐渐成为综合分区发展的重要方向。

本研究采用定量分析与GIS 空间模型相结合的方法，通过分析比较生态、经济、社会各子系统中主要影响因素的关系及其重要性，构建综合区划的多层次结构体系[23]。继而采用变异系数法和熵值法测算综合区划主导因子影响权重，引入基于加权Ripley′sK函数的多距离空间聚类进行分区。空间聚类中空间约束选择“Contiguity_Edges_Corners”，距离算法选择“Eucludean”，分组分析使用K均值算法。

1.3.3 地理探测器

地理探测器[24]是研究某现象空间分布不一，即空间变异性的统计学理论方法。该模型假设自变量对因变量产生了影响，在空间上的分布就会相应产生一定的相关性；主要包括分异及因子探测、交互作用探测、风险区探测和生态探测4 个探测器，其中分异及因子探测器即是探测地理属性Y的空间变异性和探测因子X1、X2……对Y空间变异的解释力（即贡献力），用q值度量。其分异及因子探测的表达式为：

式中：L表示因子X的分层，即自变量分类分级；Nh和N分别为层h和全区的单元数和σ2分别为层h和全区Y值的方差。SSW 和SST 分别为层内方差之和和全区总方差。q∈[0，1]，q值越大表示自变量X对属性Y的解释力越强，反之则越弱。

地理探测器现已被广泛应用于地理空间相似性与分异性研究中，如地理空间影响因素识别[25]、区域分异与空间优化[26]、生态环境与公共健康[27]等。地理探测器擅长分析类型变量，可用于测度地理环境、土壤类型、土壤理化性质及农田利用方式等类型变量对农田SOC 空间变异的贡献力，较主成分分析、多元回归等传统线性统计分析方法更为直观和快捷。

1.3.4 区域土壤有机碳空间测算方法

（1）基于地理探测器的多元复合模型

区域农田土壤类型复杂、农田利用方式多样性和多变性强[28]，导致SOC 空间分异的变量众多，且各变量对SOC 空间分异的贡献力不同。为实现区域农田SOC 高精度制图的目标，以协同克里格[29]为基础构建分层多元复合模型。选择高程、坡度、降水量、土壤类型、农田利用方式等因子[30-32]作为协同变量，分别进行协同克里格插值。同时，采用地理探测器探测不同变量贡献力大小，并根据影响力大小确定各综合特征分区主导变量及其权重系数。采用加权求和法，综合各变量协同克里格插值结果测算出SOCC[17]。其公式表达为：

式中：Z0,ck(x0)为x0处最终的SOCC 估算值；Zi(x0)为x0处第i个协同变量（影响因子）克里格插值的SOCC估算值；λi为第i个协同变量的权重系数。

（2）传统区域插值预测方法

普通克里格（Ordinary Kriging,OK）、地理加权回归模型插值（Geographically weighted regression model-Kriging,GWRK）和径向基函数神经网络（Radial basis function neural network,RBFNN）是目前应用频率最高和精度较高的3 种传统区域插值预测方法。其中，OK利用区域化变量的原始数据和变异函数的结构特点，能够实现空间预测值的线性无偏最优估计[33]。本研究采用ArcGIS 10.4 空间分析模块直接进行OK 空间插值预测制图。GWRK 是对地理加权回归模型[34]的延伸与扩展，本研究通过将SOCC 不同结构变量代入地理加权模型，以邻近样本的值和距离预测SOCC，重构插值集后进行空间插值预测制图。RBFNN 同样是一种较为有效且精确的空间插值方法[35]，与BP 神经网络相比，结构简单、训练简洁，具有较快的学习速度和逼近能力。RBFNN不仅能处理数值型变量，也能有效处理类型变量，同时能够避免过拟合，正逐渐被应用到土壤空间信息等研究中。本研究采用SPSS 22.0中的RBFNN 模型分析，输入SOCC 及其结构变量，选择自组织学习方法作为径向基函数（RBF）中心的确定方法，选择一般径向基函数（即高斯函数模型）作为隐含层激活函数，模型输出预测SOCC 值后绘制预测结果图。

2 结果与分析

2.1 农田土壤有机碳高精度制图的区域综合特征分区

面向农田SOC 高精度制图的综合特征分区主要由主题层、目标层、指标层三级构成。其中主题层包括自然地理特征和社会经济特征；目标层中自然地理特征主要包括气候状况、地形、土壤状况和水文条件等4 个因素；社会经济特征主要包括人口聚集度、土地资源状况、区位状况、交通状况、社会发展水平和经济发展水平等6 个影响因素。为避免因素间，尤其是社会经济要素之间的多重线性干扰，对指标因素进行降维收敛处理。指标层主要确定为海拔高度、土壤类型、水网密度、人口密度和耕地占比5个主导分区因子。根据变异系数法和熵值法测算分区因子贡献度，继而采用算数平均法测算出最终的分区因子综合权重，海拔高度、土壤类型、水网密度、人口密度、耕地占比的综合权重分别为0.19、0.07、0.20、0.38、0.16（表2）。

表2 研究区综合特征分区因子综合权重（贡献力）测算表Table 2 Calculation table of index comprehensive weight of comprehensive characteristic zoning in the study area

以研究区内1 624个行政乡镇（街道）为综合特征分区的基本单元，采用基于加权Ripley′sK函数的多距离空间聚类进行分区。通过对农田SOC 高精度制图样本集比较，分析不同分区策略（7～15 个区）的SOC 样本集统计特征，结果显示：当研究区划分为13个区时，样本标准偏差均值为5.17，小于其他分区策略。样本方差均值27.91，小于大部分分区策略，Moran′sI指数0.29，大于其他分区策略（表3）。综合比较结果，确定13 个分区策略为最终综合特征分区结果，其SOC 标准偏差均值、方差均值较未分区前分别下降0.55、3.53，Moran′sI指数上升0.08。13 个综合特征分区中，各分区面积在0.27×104～3.36×104km2之间，占研究区总面积的1.50%～18.66%（图2）。

图2 研究区农田综合特征分区结果Figure 2 The result of comprehensive characteristic zoning of farmland in the study area

表3 不同聚类分区策略中土壤有机质含量样本统计特征表Table 3 Statistical characteristics of SOM samples in different clustering zoning strategies

2.2 基于地理探测器的农田土壤有机碳结构变量影响力探测

农田SOC 空间变异受自然和人为因素的耦合影响，其结构变量（影响因子）复杂多样，采用文献法遴选中大空间尺度SOC 空间分异影响因子。其中气候因素中主要选取区域降雨量因子，地形条件选取海拔高度、地形坡度、坡向3 个影响因子，土壤属性选取土壤类型、表土质地、土壤pH 3 个影响因子，人类活动则主要选取土地利用方式因子。地理探测器擅长分析类别变量，将以上8个变量进行分类分级（表4）。

表4 研究区农田土壤有机碳结构变量（影响因子）分级分类表Table 4 Classification table of SOC structure variables（influencing factors）in the study area

以农田SOCC 样本为因变量（Y），采用地理探测器进行探测。结果显示，8 个影响因子之间的相关性在0.001～0.095 之间，各因子间都呈增强型非线性弱质边缘相关，排除共线可能性。13 个综合特征分区内8 个因子的SOC 影响力（q值）在0.000～0.186 之间。各分区内同一因子的影响力差异较大，剔除各分区内影响力较小的干扰因子，每个分区保留5～7个具有显示度的影响因子作为该区内农田SOC 的结构变量（表5）。

表5 不同分区内农田土壤有机碳影响因子影响力汇总表Table 5 Influences of soil organic carbon structure variables in different zones

2.3 基于分层多元复合模型的农田土壤有机碳高精度制图

根据地理探测器探测结果，综合遴选出的SOC空间变异结构变量构建不同综合特征分区内农田SOC的多元复合模型（MCM）（表6）。在研究区内划分230 m×230 m 栅格单元3.33×106个，根据13 个分区的MCM 及208 503 个SOC 样本数据，绘制农田SOCD 高精度空间分布图（图3）。结果显示研究区农田SOCD最小值为1.56 kg·m-2，最大值为6.19 kg·m-2，平均值为3.64 kg·m-2。粤北清远市、梅州市的山地丘陵区农田SOCD 较高，均值在4.20 kg·m-2以上；南部滨海平原台地区农田SOCD 较低，均值在3.00 kg·m-2以下（表7）。采用自然断点法将研究区农田SOCD 划分为高（≥4.08 kg·m-2）、中（3.33 kg·m-2≤SOCD＜4.08 kg·m-2）、低（＜3.33 kg·m-2）三个等级。其中低等级农田面积1.155 5×106hm2，占研究区农田总面积的33.86%，具有较大的农田SOC储存提升潜力。

表7 基于高精度空间分布图的农田土壤有机碳密度统计分析表Table 7 Statistical analysis of SOCD based on high-precision spatial distribution map

图3 广东省农田土壤有机碳密度高精度分布图Figure 3 High precision distribution map of SOCD in Guangdong Province

表6 分区农田土壤有机碳样本及分层多元复合模型Table 5 Sample status of SOC and multivarible composite model in different zones

2.4 与传统方法制图结果的对比分析

采用MCM 和传统方法分别进行区域SOC预测制图，基于地理探测器构建的MCM 与OK 的区域农田SOCC 预测图的预测细节较丰富；GWRK 得到的空间分布模拟结果更为平滑，而RBFNN 预测图的细节表达不够清晰（图4）。同时根据样本验证集的验证结果，以OK 插值结果为基准进行比较（表8）。其中GWRK、RBFNN、OK的综合误差较MCM增大了6.07%、9.04%和4.77%，Pearson相关系数（r）降低了6.78%、11.86%和10.17%，MCM的综合精度相对GWRK、RBFNN 和OK 分别提升了6.45%、10.45%和7.50%。在中大空间尺度上MCM 的区域农田SOC 预测结果的综合相对精度较GWRK、RBFNN 和OK 都有较大幅度的提升。

图4 研究区农田土壤有机碳不同方法的预测结果Figure 4 Prediction results of different methods of farmland SOC in the study area

表8 分层多元复合模型与传统方法对比验证Table 8 Comparison and verification between multivarible composite model and traditional methods

同时，观察验证集样本实测值与预测值的普通Q-Q 图，MCM 与OK、GWRK 的实测值与预测值较为接近直线，表达农田SOC 空间分异特征的稳定性更强。RBFNN 则与直线偏离较大，说明其实测值与预测值的空间线性相关性低于其他三种方法，空间表达的随机性较强。MCM 尾部更接近直线，其预测误差相对更为准确（图5）。OK 与GWRK 均存在“翘尾”现象，尾部预测值误差较大。RBFNN 预测结果误差也相对较大，且预测高值和低值数据时偏差较大，预测值偏向于实测值的平均值和中位值。

图5 验证集中样本实测值与预测值的普通Q-Q图Figure 5 General Q-Q diagram of measured and predicted values of samples in the verification set

3 讨论

综合可操作性、预测精度及空间表达力，MCM 有效体现了农田SOC 空间变异的多因性及不同变量的影响力差异，可以作为大尺度复杂地理情景下区域SOC预测的优先模型。

农田SOC 空间变异的影响变量众多[17]，如海拔高度和地形坡度影响不同地貌部位SOC 的累积和流失速度[36]，不同利用方式驱动下土壤改良、施肥策略、耕耙等人为扰动，及由其导致的土壤理化性质变化对农田SOC 时空变化的影响等[37]。目前研究多是通过加大采样密度[38]、复杂化预测模型与强化模型训练[11，39]等途径，实现农田SOC 制图的精度要求。但是，受地理环境与人为利用方式影响，大尺度复杂情景下农田SOC空间变异更加剧烈，单一线性模型或小尺度机器学习训练模型难以满足精度要求。通过区域综合特征分区进一步聚合区内相似性与区间差异性，降低了农田SOC空间预测的难度。同时，通过地理探测器对各影响变量的有效测度，在不同特征分区的测算模型中建立了变量贡献力序列，避免了因素间的相互干扰，抑制了插值噪声[19]。经对比，RBFNN 对环境特征和训练样本依赖性较强，空间移植性较差[39]。在大空间尺度下，OK插值预测结果较GWRK和RBFNN具有更好的内在空间关联性和精确性[40]。MCM 则清晰地将分层单变量协同克里格插值结果与基于变量影响力序列融入统一的模型，相比OK、GWRK 和RBFNN传统方法，其测算结果综合误差最小，精度提升最显著。MCM 能够在农田SOC 空间预测过程中协同更多易获取的辅助变量，同时体现不同变量对空间变异的贡献力，兼顾了测算过程中的秩序与协同，降低了预测结果的不确定性。

对于大尺度复杂情景下的农田SOC空间预测，通过集合区域综合特征分区、地理探测器、分层多元复合模型等技术手段，辅助大数据样本集[41]，可以取得较好的农田SOC 制图精度。但在构建面向农田土壤有机碳高精度制图的综合特征分区指标体系、识别SOC 结构变量及其分类分级中仍有进一步改善的空间。广东省农田SOC高精度制图，为编制更大空间尺度的农田SOC分布图探索了有效路径，为提升区域土壤固碳潜力，促进碳交易、碳中和奠定了基础。

4 结论

（1）通过构建综合区划的多层次结构体系，引入基于加权Ripley′sK函数的多距离空间聚类进行中大空间尺度综合特征分区，能够更加突显分区之间SOC 空间差异性，收敛分区内部SOC 样本的离散程度。分区后SOC样本的标准偏差均值、方差均值较未分区前分别下降0.55、3.53，Moran′sI指数上升0.08。

（2）地理探测器探测结果表明导致农田SOC空间变异的变量众多，且不同综合特征分区内结构差异较大。年均降水量、海拔高度、坡度等变量对山地丘陵区SOC 空间分异作用显著，对平原岗地区作用不显著，土地利用方式及土壤理化性质等变量对不同特征分区均具有较大的影响力。

（3）在大密度样本集支持下，协同区域综合特征分区与基于地理探测器构建的分层多元复合模型编制的高精度农田SOC空间分布图，较好地解决了中大尺度和复杂情景下SOC 空间分异规律与空间突变的同步表达矛盾，抑制了多变量插值噪声增加，其综合精度较GWRK、RBFNN和OK分别提升6.45%、10.45%和7.50%，预测结果更为准确，空间细节表达更为清晰。