基于BME-GWR 法的景观单元土壤有机碳密度空间预测

刘欢瑶，孟 岑，邹冬生，吴金水

（1.湖南农业大学资源环境学院，湖南 长沙 410128；2.中国科学院亚热带生态农业研究所，湖南 长沙 410125）

【研究意义】土壤是陆地生物圈中最大的碳库，土壤有机碳（SOC）控制全球碳循环，是重要的大气碳源［1］。作为评价土壤质量的重要指标，土壤有机碳库的研究对防治土壤退化和发展可持续农业具有重要意义。准确估算、分析土壤有机碳的储量和空间分布，为有效调控碳源/汇方向、全球气候变化以及为土壤质量综合评估提供了理论支持。土壤有机碳密度（SOCD）的估算精确度受样点数量的影响。随着3S技术的飞速发展，与SOCD相关的环境变量数据获取变得简单、便捷，充分利用多源环境先验信息作为“软数据”是克服传统以克里格为代表的地统计方法缺陷，提高SOCD估算精度的有效技术途径。

【前人研究进展】当前已有很多插值方法结合环境辅助变量进行土壤属性的空间预测，包括多元线性回归模型、协同克里金、回归克里金、地理加权回归等，当环境因子与土壤属性相关性较强时这些方法预测精度均高于不结合环境辅助变量的地统计方法［2］。然而，由于多数插值方法均未考虑环境辅助变量的先验分布，很难处理非高斯分布变量，因而限制了利用多源辅助数据进行空间预测的能力［3］。贝叶斯最大熵（Bayesian Maximum Entropy,BME）结合了贝叶斯方法论和最大熵原理来模拟时空变量，将精度较高的“硬数据”和多渠道源得到的“软数据”区分使用，能够有效地综合利用各种不同来源和精度的数据［4］。近几年来，BME模型尝试被引入到土壤、大气污染和环境的时空预测中［4-6］。例如，费徐峰等［5］将田间测量的土壤重金属含量数据作为软数据，将实验室测定的相应土壤重金属含量作为硬数据，利用土壤重金属在土壤母质类型中的概率密度函数作为软数据，建立了土壤重金属含量BME和普通克里格（OK）插值模型，以探讨研究区内重金属污染分布情况和影响因素。

【本研究切入点】软数据本身的来源多样化，例如将另一种检测手段得到的数据集作为软数据，按环境相关法得到预测属性的概率分布作为软数据［7］，用空间预测模型结合辅助变量获得相对不确定性的数据作为软数据［4］，但较少有文献报道软数据的获取途径对BME结果的影响。受复杂的成土因素（如植被、地形）等环境因子影响，随着研究尺度的变化，土壤属性的主控因素也会随之变化［8］。因此，选择合适范围建立环境辅助数据和预测变量之间的关系，是有效控制预测模型不确定性产生的途径。【拟解决的关键问题】本研究基于BME方法框架，选取湖南省桃源县盘塘镇王家垱村内SOCD实测数据作为“硬数据”，以地理加权回归模型（GWR）结合地形、土地利用等多源环境数据获得“软数据”，比较贝叶斯最大熵结合地理加权回归模型（BMEGWR）与传统GWR方法对土壤SOCD空间分布预测精度，进一步分析基于土地利用类型估算得到的“软数据”所计算的BME-GWR模型模拟精度的变化，为综合利用环境变量提高土壤属性的空间预测模型精度提供依据。

1 材料与方法

1.1 研究区概况

综合考虑地形地貌、土地利用、土壤等因素，在我国亚热带红壤丘陵区选取具有代表性的农业生态景观单元作为研究区域［9,11］。研究区域位于湖南省桃源县盘塘镇王家垱村，总面积为355.12 hm2，地处中亚热带北缘，为季风性湿润气候（年平均降雨量为1330 mm，年平均气温为16.8 ℃）。该区域属于低山岗地地貌，高程约为81～112 m，土壤类型以地带性的红壤为主。研究区土地利用方式主要为林地（35.7%）、稻田（39.9%）、果园（15.1%）、旱地（9.3%），其中林地以人工林（马尾松、栎树、樟树）为主，稻田以双季稻为主，旱地主要作物为棉花、油菜、玉米、苎麻等经济作物，果园以橘子、柚子、桃子为主［10］。

1.2 研究方法

1.2.1 数据来源与处理 通过GPS准确定位采样点，并记录采样点的经纬度、高程，通过ArcGIS 10.3软件经投影转化，得到以米为单位的平面坐标；再将矢量等高线和高程控制点转化为分辨率为5 m的数字高程模型（图1），并提取1:10000航拍的土地利用信息，得到分辨率为5 m土地利用类型分布图，在数字高程模型中提取坡度、地形湿度指数等数据。于2003年4—5月，以样区中心为横纵坐标轴的交叉点，划定1.0 cm×1.0 cm的网格，按采样密度耕地（稻田和旱地）5个土样/ hm2、果园2～3个土样/ hm2、林地1个土样/ hm2的原则随机布点，并采用“S”形路线多点取样法（不少于15点）采集表层土样（0～20 cm），共采集土壤样品523个［11］，其中，稻田、旱地、果园、林地样品数分别为249、66、192、16个。

土壤容重采用环刀法测定，土壤有机碳含量采用碳氮元素分析仪（Vario-MAX C/N，德国）测定。SOCD指单位面积内一定厚度土层的土壤有机碳储量（kg/m2），计算公式为：

式中，SOCDi为第i层土壤有机碳密度（kg/m2）；SOCi为第i层土壤有机碳含量（g/kg）；ρi为第i层土壤容重（g/c m3）；Di为土壤厚度（cm），本研究取值20 cm。

图1 研究区样点分布、高程、坡度、地形湿度指数、土地利用类型Fig.1 Sample plot distribution,elevation,slope,topographic wetness index（TWI）and land use types in the research area

1.2.2 地理加权回归（GWR）GWR模型是普通线性回归的拓展，可被视作局部的加权最小二乘回归模型，允许局部参数估计，将观测数据空间特性的变化参数嵌入到模型中，并假定线性回归模型中的回归参数是观测点处的任意函数，得到GWR模型，反映样本对回归方程贡献在空间上的分异［13］，其计算公式为：

式中，y（u）为位置u的因变量值；β0（u）为截距值；βk（u）为第K个协变量在位置u上的回归系数，回归系数由每个观测点周围数据估算而得；xk（u）为在位置u上的第k个协变量的值；ε（u）为在位置u上的随机误差值。GWR回归模型中的参数采用加权最小二乘法计算，本研究把SOCD测定数据作为因变量，高程、坡度、地形湿度指数、土地利用类型数据作为自变量进行GWR模型拟合。

1.2.3 贝叶斯最大熵模型（BME）BME模型结合了最大熵原理和广义贝叶斯条件公式，可以利用不同精度、质量、表示方式（如区间、分类等）的多源数据，这些数据总体上可以分为广义知识（General knowledge，G）和特定知识（Special knowledge，S）两类。广义知识用来描述自然规律、物理法则和硬数据的统计规律（如数学期望、半方差）等。特定知识包括硬数据和软数据，硬数据是特定观测点上的误差可以忽略的数据，如实测数据、时间序列数据等，本研究以样点SOCD实测数据作为硬数据；而软数据是相对模糊的、有误差的数据，包括历史数据、粗测量数据、模型模拟数据等。

BME的基本原理是利用广义知识（G）通过最大熵原理获得先验概率密度函数pdf〔fG（ymap）〕，再基于结合了硬数据和软数据的特定知识（S）、贝叶斯条件概率将研究区未测点变量分布的先验pdf转换为后验pdf〔fT（ymap）〕，并可计算待预测点xk的预测均值或最大值。

式中，I=Imn+1U Imn+2U…Im是软数据集，ydata=[y1，y2，…ym]是特定知识。为预测结果的平均值；为预测结果的最大值；yk为测量值；f（yk）为后验pdf；dyk为测量值的变化量。

本研究地理加权回归结合贝叶斯最大熵模型（GWR-BME）是以GWR模型所预测连续栅格表面为基础，将在样点空间位置上提取的被预测的SOCD属性值作为软数据；而按土地利用类型估算的贝叶斯最大熵模型（GWR-BMEL），则是先提取不同的土地利用类型，在每个土地利用类型单独以SOCD实测数据作为因变量，高程、坡度、地形湿度指数作为自变量进行GWR模型预测，提取样点位置的GWR模型预测值作为软数据。最后，结合“软数据”和“硬数据”进行BME预测，GWR-BMEL则要将各土地利用类型的预测结果合并得到整个区域的SOCD空间分布。

2 结果与分析

2.1 土壤有机碳密度与环境变量描述性特征统计

从表1 可以看出，研究区域的SOCD 分布范围为1.30～4.71 kg/m2，平均值为3.22 kg/m2，标准差为0.77 kg/m2，变异系数为24%。不同土地利用类型下SOCD差异显著，表现为稻田〔3.84（±0.44）kg/m2〕>旱地〔2.92（±0.57）kg/m2〕>果园〔（2.58（±0.44）kg/m2〕>林地〔2.31（±0.55）kg/m2〕，变异系数为11%～24%，均属于中等变异（10%～100%）。

通过数字高程图提取高程、坡度、地形湿度指数作为地形特征，与SOCD 进行相关分析。从图2 可以看出，研究区域内提取的高程、坡度、地形湿度指数对SOCD 有明显的影响，可以作为地形因子进行空间分布拟合。其中，SOCD 与高程和坡度均呈极显著负相关、相关系数分别为-0.68 和-0.62，而与TWI 呈极显著正相关、相关系数为0.34。

表1 研究区土壤有机碳统计Table 1 Statistical values of soil organic carbon density（SOCD）in the research area

图2 研究区土壤有机碳密度、高程、坡度、地形湿度指数相关分析Fig.2 Correlation analysis of soil organic carbon density（SOCD）,elevation,slope and topographic wetness index（TWI）in the research area

2.2 半方差函数比较

根据最大决定系数（R2）原则，GWR、GWR-BME、GWR-BMEL3种模型最优拟合模型均为指数模型（图3）。半方差函数模型中所得的主要参数为块金值（C0）、基台值（Sill）、变程（Range）以及块基比值（C0/Sill），可以反映一定空间范围内的空间变异性和变量的自相关性。其中，块金值反映了随机部分形成变量的变异性和测量误差，GWR、GWR-BME、GWR-BMEL3种模型块金值分别为0.32、0.24、0.16 kg/m2，结合BME方法可以减少由随机部分形成的变异；基台值是变异函数数值随着采样点间距的增大从块金值增大到最大的常数值，反映变量的最大变异性，体现了非随机部分形成的变异，GWR、GWR-BME、GWR-BMEL的基台值分别为0.61、0.66、0.52 kg/m2；块金值与基台值的比值分别为52%、36%、31%，为中等程度的空间依赖性（25%～75%），尤其是GWR-BMEL表现出的空间自相关性最强；变程从大到小依次为GWR（269 m）、GWR-BME（207 m）、GWRBMEL（138 m），GWR-BME、GWR-BMEL最优拟合模型变程均小于GWR，进一步验证了块金值与基台值比值的结果，反映出较强的空间自相关性。

2.3 预测精度检验

图3 半方差模型对比Fig.3 Comparison of semivariance model

其中，MAE从大到小依次为GWR（0.25）>GWR-BME（0.21）>GWR-BMEL（0.19），RMSE从大到小依次为GWR（0.40）>GWR-BME（0.35）>GWR-BMEL（0.33）。GWR-BMEL模型的MAE和RMSE比GWR-BME模型分别减少9.53%、5.71%，比GWR模型分别减少24%、17.50%。R2以GWR模型最小（0.72），GWRBME（0.79）比GWR提高9.72%，GWR-BMEL的R2最大、为0.81。表明在软数据的辅助下，GWRBME和 GWR-BMEL的模型精度均高于GWR模型，而GWR-BMEL模型在3种方法中插值精度最高，GWR-BME次之，GWR的模型精度最低（表2）。说明结合了土地利用类型划分估算单元所得到的软数据可以较好地体现土地利用类型内部的空间差异性，克服整体估算结果的平滑效应，提高了估算结果的精度。

表2 各土壤有机碳密度模型精度比较Table 2 Comparison of accuracy among different soil organic carbon density（SOCD）models

2.4 土壤有机碳密度空间分布

图4 土壤有机碳密度（SOCD）空间分布Fig.4 Spatial distribution of soil organic carbon density（SOCD）

利用GWR-BMEL方法预测研究区域的土壤SOCD的空间分布（图4），得到分辨率为5 m的空间预测结果。从图4可以看出，GWR-BMEL法拟合的SOCD空间分布特征的整体趋势特征明显，反映细节的能力强，对实测点邻域范围的空间结构体现较合理。农业生态景观单元内SOCD的预测结果具有明显的区域分布特征，SOCD的空间分布格局与土地利用、地形因子的相关性明显：在研究区域中部高程和坡度较小的坡度、坡下的稻田以及稻田周围的旱地的SOCD较高；而研究区域东部高程和坡度较大的坡腰、坡顶区域，主要分布林地和果园的区域是SOCD的低值区。

3 讨论

我国全国尺度土壤有机碳密度的平均值为3.35 kg/m2［14］，东部地区表层土壤的SOCD为6.8～21.4 kg/m2［15-16］，其中浙江和江西两省典型森林类型土壤平均SOCD为7.12～15.69 kg/m2［17］。而本研究区域的土壤有机碳密度比其他研究结果较低，这可能是由于农业生态景观单元地表的作物成熟后多被收割，直接归还到土壤部分的凋落物相对减少导致土壤有机碳含量较低［18］，且研究区域的林地以人工经济林为主，林分结构单一，植被净生产力小于天然林［16］，同时人工林受人为扰动较大，土壤微生物的分解作用加强，限制了当地土壤有机质的积累［19］。另外，研究区以低山岗地地貌为主，耕地多位于低山的坡中和坡下位置，水土流失较严重，导致土壤相对贫瘠［20］。

本研究区SOCD的GWR方法最优拟合模型为指数模型，这与很多相关的土壤有机碳研究结果一致［21-22］。研究区土壤有机碳密度的块金值和基台值分别为0.32、0.61 kg/m2，块金值与基台值的比值为52%。块金值与基台值的比值受结构和随机因素的影响，较高的比值说明空间变化是由施肥、耕作、其他人类活动等随机因素引起，而较低的比值说明由采样和实验分析本身的系统误差、土壤性质、地形、其他自然因素等结构因素在空间变化中起重要作用。一般用<0.25、0.25～0.75、>0.75分别表示土壤属性的强、中等、弱空间自相关性，本研究区的土壤有机碳密度为中等程度空间变异［23］。土壤有机碳密度半方差模型的变程为269 m，是土壤有机碳密度的最大相关距离，表明本研究SOCD采样的网格间距布设合理，推导出的空间结构关系与地面的空间结构关系一致［24］。与此相应，SOCD的GWR-BME、GWR-BMEL两种拟合模型的块金值、块金值与基台值的比值、变程均小于GWR，表明BME法能更强地解释SOCD的空间异质性。同时，GWRBMEL的基台值最小，体现了同一类型的景观约束下SOCD的变异性会降低［25］。

结合BME方法用于模拟SOCD空间分布的精度高于GWR，这与多数研究结果相同。杨勇等［7］以DEM生成的相关地形因子作为环境数据，分别用普通克里格（OK）和BME模型对湖北省京山县土壤有机质含量的密集样本和稀疏样本进行预测，结果表明BME的预测精度高于OK，尤其对稀疏样本预测的精度提高幅度更大。Xiao等［4］采用GWR-BME模型，选取气溶胶光学厚度、地形数据、气象数据和污染排放等辅助信息对我国大区域尺度的PM2.5空间分布进行拟合，GWR-BME预测结果精度较GWR模型更高。与传统方法相比，GWR-BME、GWR-BMEL方法结合了“硬数据”和“软数据”［6］，降低了预测结果的不确定性，尤其是通过划分土地利用类型获取软数据的GWR-BMEL模型，通过在空间范围上对每个土地利用类型单独模拟，获得其空间变化的局部特征以及空间对象本身相关性和异质性参数，进一步有效利用“软数据”的空间自相关性，可以较好地体现研究区SOCD的空间分布。

本研究利用GWR-BMEL对研究区SOCD模拟分布，结果表明在中间高程和坡度较小的稻田以及稻田周围的旱地SOCD较高。这是由于农田（稻田和旱地）的复种指数高，化肥投入量高，导致表层土壤的养分较高，导致SOCD高于林地和果园，尤其是稻田由于淹水环境会降低土壤有机碳的分解速率，进一步降低土壤有机质的分解速率，有利于土壤有机碳积累［25］。同时，研究区域的水稻和旱地作物经常种植在坡度、海拔较低的区域［26］，因而形成当地SOCD的高值区域集中分布在坡度、海拔较低的稻田和旱地的格局。高程和坡度较大的林地和果园是主要的SOCD低值区域，其中研究区域海拔和坡度较大的林地种植年限较短，以种植马尾松和油茶等人工林为主，其净生产力远低于天然林，有机物质投入量少［26］；当地的果园多种植桔科植物，因而地面凋落物较少，从而减少了土壤有机碳来源。而且，受耕作措施的影响，果园的深埋施肥对表层土壤的有机碳贡献也不明显，果园土壤连年深翻加速了土壤有机碳矿化过程［27-28］。另外，与农田地形较平坦相比，果园和林地常分布山地，具有较高的水力传导性，土壤有机质容易随着降雨而发生流失；这些原因也可能会导致SOCD的降低［29］。

4 结论

本研究中，由于BME-GWR和BME-GWRL模型结合了“软数据”作为辅助变量，对研究区SOCD的空间异质性有更强的解释能力，交叉验证结果也反映BME-GWR和BME-GWRL模型较GWR模型的拟合精度更高，其中BME-GWRL模型基于土地利用类型的空间范围对软数据进行拟合，考虑到了“软数据”估算单元的不确定性，对比不划分土地利用类型直接模拟的BME-GWR预测结果，从整体上更加接近于观测数据。本研究结果表明，考虑软数据的获取途径，是在实测数据和辅助变量有限的条件下提高空间预测精度的有效途径，为合理利用多源辅助数据提供新思路。