典型山西黄土高原区土壤有机质的空间异质性及空间插值预测

张小美 ，高春瑞 ，闫晓斌 ，杨 莎 ，乔星星 ，王 超 ，杨武德 ，Fahad Shafiq ，冯美臣 ，李广信

（1.山西农业大学 农学院，山西 太谷 030801；2.拉合尔政府学院大学植物学系，巴基斯坦 旁遮普省 54000）

土壤作为地球陆地生态系统中最大的碳库，是陆地生态系统的核心，是碳封存的关键，也是全球碳循环的重要通道[1]。据估计，全球1 m 深度的土壤有机碳储量约为1550 Gt，大约是大气碳库的2 倍和陆地植物碳库的3 倍[2]。土壤有机碳（Soil Organic Carbon，SOC）是土壤肥力形成的基础，同时也是表征土壤质量和功能、减缓全球气候变化的重要因素[3-5]。因此，明确区域SOC 含量的空间分布及其影响因素，对于研究SOC 的空间差异性分布及变化具有重要意义[6]。

获取SOC 空间分布信息的传统方法主要通过实地采样、预处理和化学化验等过程，存在耗时、费力、低效的缺点，难以满足大样本、大尺度、快速高效的现实需求。空间插值方法可以将SOC 含量从离散的点状信息转化为面状连续信息，用来表征SOC 空间分布特征以及预测采样点分布区域的土壤属性空间信息，是较为常用的SOC 含量估算方法[7-8]。但是应用不同的插值方法得到的土壤有机碳的空间分布也会有所不同，王翔[9]研究发现，普通克里格插值可以很好地利用采样点模拟空间聚集度；谢梦姣等[10]研究发现，径向基函数插值神经网络法可以避免普通克里格插值的平滑效应现象，提高土壤有机质空间预测精度；王海菁等[11]在宜黄县城的洪水淹没分析中表示，普通克里格插值相较于反距离权重法而言能使误差的方差更小，更好地预测洪水的淹没范围。因此，探索预测精度更高的插值方法对实现土壤有机碳的高效预测具有重要意义。

SOC 的含量因易受到土壤形成因素（如环境因子、植被因子、水文因子、气候因子和地形因子）的影响而具有很强的空间异质性。因此，当研究SOC 空间分布的区域尺度较大时，可能会存在影响因素多且不容易明确哪些因素影响比较大的问题，尤其是当气象站点与采样点距离较远，极易造成较大的误差，但研究小区域则能尽可能的减少这种误差[12-13]。同时，基于单一因素的单一空间插值方法也会忽略部分环境因素带来的影响，导致制图精度降低[14-15]。为了能够更好地捕捉环境变量在空间上对SOC 的影响，文雯等[16]、张宏帅等[17]将地形因子和气象因子等辅助变量用于构建区域或国家尺度的SOC 空间插值模型，取得了较好的效果；杨顺华等[18]研究发现，结合回归克里格和空间位置的地理加权回归克里格能更好地模拟SOC 空间分布格局；QIN 等[19]结合环境推理模型与反距离权重模型也提高了SOC 的作图精度；江叶枫等[20]研究也指出，结合多个环境因素的回归克里格模型精度要优于单一的普通克里格模型。因此，利用多模型结合对SOC 的时空变异性研究具有极大的应用潜力。

本研究以山西省忻州市忻府区农田SOC 为对象，通过引入温度、降水、坡度、高程为环境辅助变量，在利用反距离权重法、径向基函数插值、普通克里格插值、多元线性回归、回归克里格、回归反距离权重法和回归径向基函数插值多种模型插值方法的基础上，探讨小尺度区域条件下SOC 与环境变量的关系，定量评价不同插值方法的表现，寻找最优的预测方法并进行空间制图，以实现小区域尺度条件下农田SOC 的精确预测制图，旨在为未来大尺度农田土壤碳库估算提供一定的理论和技术参考。

1 材料和方法

1.1 研究区概况

忻州市忻府区位于山西省中北部（北纬38°6'～39°40'、东经110°53'～113°58'），属季风型大陆性气候，占地面积为1986.53 km2，整体地貌是被黄土广泛覆盖的山地型高原，忻府区地形西高东低，逐步倾斜，北、西、南三面环山，东部开阔平坦，是忻定盆地的主体部分。降水集中于7～9月，年均降雨量在462.5 mm 左右，无霜期平均167.1 d，年均温在5.7 ℃，主要种植玉米、谷子、糜子、向日葵等作物。

1.2 数据来源

1.2.1 SOC 数据及预处理 本研究的SOC 数据来源于山西自然科技资源共享平台（http://www.sxzrzy.ac.cn/），共获取907 个土壤样点数据。该数据是通过GPS 定点并实地考察，通过土种图与土地利用现状图叠加进行评价，根据评价单元的个数及相应面积，在样点总数的控制范围内进行合理布局，并根据图斑大小、种植制度、作物种类、产量水平等因素确定评价单元布点数量和点位，同时根据土壤类型、生态条件、种植面积等因素，选择有代表性的土壤，采取定点取样和随机取样相结合的方法进行布局及采样，采样深度为20 cm，利用外加热-重铬酸钾氧化法测定获取土壤有机质含量。为检验不同插值方法在研究区域的预测精度，运用Arc-GIS10.2 地统计工具下的Subset 将907 个实测样本按照3∶1 的比例分为预测集和验证集。

1.2.2 环境变量数据及预处理 选择气象因子（降水量和温度）和地形因子（高程和坡度）作为辅助SOC 预测的环境变量（图1），其中，气象数据来自中国气象数据网（https://data.cma.cn/），高程及坡度数据来源于山西自然科技资源共享平台（http://www.sxzrzy.ac.cn/），将气象数据进行整理和预处理，用Excel 计算出年均温和年均降水量，利用Arc-GIS 10.2 提取出各采样点的温度和降水量。

图1 环境变量空间分布Fig.1 Spatial distribution map of environmental variables

1.3 研究方法

本研究用到了非地统计、地统计、互结合三大类方法，其中，非地统计主要包括径向基函数、反距离权重插值、多元线性回归；地统计为普通克里格插值；互结合是基于非地统计和地统计的拓展方法，主要包括回归克里格法、回归反距离权重法和回归径向基函数插值[19]。本研究中用到的7 种插值方法都可以通过ArcGIS 10.2 实现。

1.3.1 反距离权重法（Inverse Distance Weighting Method，IDW） 其是对插值点与样本点之间的距离进行加权平均，基于相近相似的原理距离越短，权重越大[16]。

1.3.2 径向基函数插值（Radial Basis Function Method，RBF） 其是某种沿径向对称的标量函数，通常定义为样本到数据中心之间径向距离（通常是欧氏距离）的单调函数[21]。

1.3.3 普通克里格插值（Ordinary Kriging，OK）其是依靠所选采样点数据在空间上加以插值，其插值效果较依赖于采样点的个数、密度，以及其数据的准确性[22]。

1.3.4 多元线性回归（Multiple Linear Regression，MLR） 其是对2 个或者多个变量进行回归分析，本研究用预测集在SPSS 中进行逐多元逐步线性回归，可以定量描述环境变量对SOC 的相关关系，并且避免了多重共线性[23]。

1.3.5 回归克里格（Regression Kriging，RK） 其是MLR 与OK 的结合，首先需要通过MLR 得到残差值，再用OK 对残差进行插值，并最终将MLR 预测的趋势项与OK 的插值结果相加，得到RK结果[22]。

1.3.6 回归反距离权重法（Regression Inverse Distance Weighting Method，MIDW） 其是MLR 与IDW 的结合，通过MLR 模型得到残差值，再用IDW 对残差进行插值，最后将MLR 预测的趋势项与IDW 的插值结果相加，从而得到MIDW 结果。

1.3.7 回归径向基函数插值（Regression Radial Basis Function Method，MRBF） 其是MLR 与RBF的结合，通过MLR 模型得到残差值，再用RBF 对残差进行插值，最后将MLR 预测的趋势项与RBF的插值结果相加，从而得到MRBF 结果。

1.4 检验方法

为检验7 种插值方法的预测精度，本研究通过预测集中的样点进行插值预测得到结果图，再根据验证集的位置提取出验证集的预测值，对验证集的实测值和预测值进行相关性分析，通过计算得到预测值和实测值的Pearson 相关系数（r）、平均绝对误差（MAE）、均方根误差（RMSE）作为检验插值精度的标准。MRE 和RMSE 值越小，说明结果越精确。

2 结果与分析

2.1 SOC 含量的描述性统计特征

分别对总样本（n=907）、预测集（n=680）和验证集（n=227）的SOC 含量进行描述性统计分析，结果如图2 所示。

图2 土壤有机碳含量基本描述特征Fig.2 Basic description characteristics of soil organic carbon content

从图2 可以看出，总样本SOC 的平均值为15.05 g/kg，最大值为43.4 g/kg，最小值为1.4 g/kg；预测集的平均值为14.66 g/kg，最大值为43.4 g/kg，最小值为1.4 g/kg；验证集的平均值为14.66 g/kg，最大值为36.2 g/kg，最小值为3.9 g/kg。总样本变异系数为30%，即SOC 变化幅度较大，有较强的空间异质性。预测集的全距为42，验证集的全距为32.3，预测集的全距涵盖了验证集的全距，从统计学的角度可知，该数据划分合理。由图2 可知，该数据呈近似正偏态分布特征，表明数据符合进一步分析的需要。

2.2 SOC 与环境变量的相关关系

利用所提取的坡度（S）、高程（E）、温度（T）和降水量（P）与SOC 含量进行相关性分析，结果如图3 所示，在地形因子中，高程（-0.255**）和坡度（-0.214**）与SOC 含量呈极显著负相关（P＜0.01）；在气候因子中，降水（-0.085*）和温度（-0.246**）与SOC 呈显著（P＜0.05）或极显著负相关（P＜0.01）。结果表明，温度与SOC 相关性最高，降水与SOC 相关性最弱。

图3 环境变量与土壤有机碳的相关关系Fig.3 Correlation between environmental variables and SOC

2.3 7 种插值模型的构建及预测

将环境因素与SOC 数据进行多元线性回归分析，得到多元线性回归模型系数的基本统计特征，结果如表1 所示，MLR 模型的非标准化系数为高程、年均降水量、年均温和坡度，4 个环境变量的方差膨胀因子均小于7.5，则说明这4 个环境因子相互不存在共线性。此外，通过t值和P值可以看出，环境变量与SOC 之间存在显著相关关系（P＜0.05）；从标准系数看，温度对SOC 影响最大。MLR 模型的R2=0.484，即环境变量可以解释SOC 48.4%的变化。以此为基础，构建RK、MIDW 和MRBF 模型，利用MLR 模型得到残差值，分别将残差值进行OK、IDW、RBF 插值，利用栅格计算器加上MLR图层的有机碳含量值，最终得到RK、MIDW、MRBF的插值模型。

表1 多元线性回归模型系数的基本统计特征Tab.1 Basic statistical characteristics of the coefficients of the multiple linear regression models

图4 分别为MLR、RK、IDW、RBF、MIDW、MRBF 得到的SOC 空间插值预测情况，从整体来看，7 种插值模型所反映出SOC 含量的空间分布没有明显差异，但在采样点所在区域附近还是有所不同，其中，IDW、RBF、MIDW、MRBF 相对来说所得到的预测图不够平滑，呈现“牛眼”现象，能够很明显地看出采样点的位置；对比这4 种模型的SOC 空间预测图可以看出，结合MLR 模型的MIDW、MRBF 模型相对于另外2 种更加平滑。

图4 基于7 种模型的土壤有机碳空间分布预测Fig.4 Soil organic carbon spatial distribution prediction map based on 7 prediction models

2.4 插值结果的验证与分析

从表2 可以看出，除OK 空间插值方法外，其余6 种方法的预测值与验证值之间具有极显著相关性（P＜0.01），相关系数大小为r（MRBF）＞r（RBF）＞r（IDW）＞r（RK）＞r（MIDW）＞r（MLR）；再对比7 种模型的MAE 和RSME 值可知，其中，MLR 的RMSE 最小，小于其他6 种模型，MRBF 的MAE 最小，因此，MRBF 模型预测精度最高，在对土壤有机碳进行空间插值时可以将MLR 模型结合地统计模型来预测环境变量下土壤有机碳的变化。

表2 7 种模型的精度评价指标比较Tab.2 Comparison of precision evaluation of seven models

3 讨论

3.1 环境变量与SOC 之间的关系

本研究选用气候、地形作为环境变量应用于模型，其中，气候变化影响是SOC 浓度的重要因素[24]；温度会影响植物的光合作用速率、凋落物分解速率及土壤有机质矿化速率，对土壤碳循环有着显著影响[25]。本研究通过线性分析发现，温度与SOC 呈显著负相关，这可能是因为较低的温度更有助于SOC的积累。罗梅等[13]研究基于环境变量的SOC 空间分布特征时通过相关性分析发现，温度与SOC 呈负相关，与本研究结果一致。周莉等[26]研究指出，气候条件会通过制约植被类型、影响植被生产力进而影响输入土壤有机碳的输入量；同时还能通过影响土壤水分和温度等条件，进而影响微生物对有机碳的分解和转化。降水通过影响陆地土壤碳密度和土壤的通气性影响土壤持有的有机碳量[27]，当土壤水分充足时，则透气性差，原有机碳不易矿化，外源有机残体在水分作用下更易腐烂降解为有机质，从而提高SOC 含量；土壤水分缺失时，孔隙大，则会加速有机碳分解，不利于SOC 的积累。山西省属于黄土高原，地质干旱，降雨量小，因此，降雨量与土壤有机碳的含量呈负相关[28]，与本研究结果一致。

此外，地形因子能调节热量和物质的再分配，对SOC 分布起重要作用。其中，坡度是影响土壤侵蚀的重要因素，受到土壤侵蚀的影响，不同坡度下SOC 含量存在较大差异[27]。本研究中，坡度与土壤有机碳呈显著负相关。土壤属性的空间变异受多种因子的控制，而地形作为五大成土因素中的重要因素，能调节热量和物质的再分配，对SOC 分布起重要作用[28]。李龙等[27]研究指出，坡度是影响农地固碳潜力的主要因素，对于不同的坡度单元，坡度大于25°的耕地退耕后的固碳效应最高，而坡度小于15°的效应最低。本研究中，坡度最大值为50°，大于25°的约占样本数的1/30，因此，就坡度而言该地区的固碳潜力有待大幅度提升。本研究区域的高程与SOC 含量呈显著负相关，高程自西向东呈现降低趋势。高海拔地区通过影响温度和降雨量制约土壤微生物的活性，进而影响土壤有机碳的分解速率[14]。高程条件通过影响区域的水热条件来影响SOC 含量的累积[29]。但在本研究中，SOC含量自西向东逐渐升高，主要原因是因为该地区属于黄土高原地区，较为干旱，温度较低。在小尺度上气候、植被以及地形因素在时间和空间上基本一致，没有太大变化，样点虽然很多，但是人为因素方面的影响仍然是不可准确估量的。因此，在大尺度上SOC 的空间分布仍然存在不确定性。

3.2 模型预测精度差异性对比

本研究使用MLR、RK、MRBF 和MIDW 模型结合地形和气候因子预测SOC 含量的空间分布，对比没有考虑地形和气候影响的OK、RBF 和IDW模型可知，MRBF 模型精度最高。RK、MRBF 和MIDW 模型包含权重函数的MLR 变体，考虑了SOC 与环境变量之间的关系，MLR 是通过对预测集进行统计分析建模，预测精度主要受到数据量和数据精度的影响，但是OK、RBF、IDW 直接利用预测集进行插值，只考虑了样本点自身的结构特征，忽略了环境变量的影响，因此，为了得到更精确的小尺度土壤有机碳的空间预测，应充分考虑土壤形成因素带来的各种影响，MLR 可以得出多个环境变量对SOC 的影响，将MLR 与OK、RBF、IDW 结合得到的RK、MRBF、MIDW 模型即能更综合地考虑到SOC 与环境变量的关系[17,30]，并且RK、MRBF、MIDW 模型精度相较于OK、RBF、IDW 也有所提升，江叶枫等[19]研究也指出，考虑辅助变量的模型会提高模型精度，更加证实了互结合模型在土壤有机碳空间预测方面的优势，但是RK、MRBF、MIDW插值相对OK、RBF、IDW 等插值方法，步骤繁多，参数设置要求较高，耗时耗力，在研究目的许可的情况下，当SOC 数据样点密度足够大且不需要结合环境变量进行预测时，可以直接考虑简单灵活的OK、RBF 或者IDW 插值，预测精度相对较高[31-32]。本研究表明，互结合模型在土壤有机碳空间预测方面有很好的应用潜力；随着信息技术的不断发展，这种互相结合的模型也将成为一种新的趋势，本研究采用的是MLR 与其他模型的相互结合，未来可以考虑更多算法模型的相互结合，更好地达到预测效果；此外，本研究只考虑了气候和地形因子对土壤有机碳的影响，考虑到土壤形成过程受到多种因素的影响，未来可以考虑加入更多的影响因子，更好地完善土壤有机碳的空间预测；如果是需要构建生态模型或是图形叠加等空间分析，则需要考虑更多因素在空间上的互相影响，选择插值精度更高、在空间上表现最好的空间插值方法进行预测[20]。

4 结论

本研究以高程、坡度、降水量、温度为辅助变量，比较OK、RBF、IDW、MLR、RK、MRBF、MIDW等7 种模型的建模结果，研究结果表明，SOC 与高程、坡度、降水量以及温度呈极显著负相关关系；其中，温度和高程与SOC 相关性最高，降水与SOC 相关性最弱；利用插值模型进行SOC 空间分布预测制图，发现MRBF 模型具有较好的拟合精度，且预测结果优于其余6 种模型。该结果为促进土壤固碳提供了科学依据，并对理解小尺度SOC 的空间异质性及精确制图提供了一定的理论和实践参考。