黑土区小流域土壤速效磷空间分布模拟方法

王雪珊, 沈庆松,, 高凤杰, 张兴义, 张少良, 王 力

(1.东北农业大学 资源与环境学院, 哈尔滨 150030; 2.中国科学院 东北地理与农业生态研究所,哈尔滨 150081; 3.西北农林科技大学 水土保持研究所, 陕西 杨凌 712100)

土壤速效磷(Available phosphorus,AP)是衡量土壤供磷能力和磷素流失风险的重要指标[1-2]，该指标能够反映农田土壤肥力特征，研究其空间分布规律对于区域施肥管理和农业面源污染防控等具有重要意义[3]。受生物地球化学循环过程、自然环境和人类活动等因子影响[4-5]，农田土壤AP在空间分布上存在高度空间异质性。研究证实AP预测的精度直接影响农田土壤磷素的管理水平[6]，因此利用有限信息，改进模拟方法提高AP空间分布预测模型的精度十分迫切。

空间插值方法是利用采样点数据通过数学模型对研究区内未知点进行模拟预测，并能够获得连续的属性值空间分布特征[7]。回归克里格法[8]将多元线性回归模型与克里格方法相结合，同时引入多个环境因子参与建模，并考虑土壤属性空间分布过程中的结构性因子与随机性因子，其空间插值结果更加精确[9-10]。Li[10]和Zhang等[11]提取地形和土地利用等变量利用回归克里格法预测土壤有机质的空间分布特征，研究结果表明回归克里格方法要优于单变量空间插值方法。但是空间非平稳性的存在导致回归克里格方法不能有效地捕捉土壤属性空间变异的局部特征[12-13]。地理加权回归模型利用局部加权最小二乘法进行逐点参数估计和空间预测，能够在引入环境变量的同时基于空间非平稳性假设揭示更多的局部信息，该模型已逐渐成为土壤属性空间预测的有效手段[14-15]。地理加权回归克里格将地理加权回归模型与克里格方法相结合，既能够考虑空间信息的局部特征，还能利用克里格法对代表随机性的残差进行插值，插值结果能够揭示可能被空间非平稳性所掩盖的一些局部变化，反映出更加真实的土壤属性空间变异情况[16-18]。杨顺华等[16]引入相对高程和径流强度指数预测乡镇尺度(235 km2)土壤有机质空间分布规律；马泉来等[17]研究表明海拔和水系距离可以作为辅助变量预测黑土区小流域(119.16 km2)土壤有机质的空间分布规律，研究证实相对于普通克里格和回归克里格法，地理加权回归克里格能够显着提高土壤有机质的模拟精度；Kumar等[12]引入海拔、坡度、植被指数、土壤类型、土地利用类型等变量研究宾夕法尼亚州(117 599 km2)土壤有机碳的空间分布规律；Ye等[19]研究表明海拔、坡度和地形湿度指数能够模拟区域尺度下(16 400 km2)不同采样点密度土壤有机碳的空间分布格局，以上研究证明考虑了空间非平稳性和残差空间自相关的地理加权回归克里格方法能够加强土壤有机碳的模拟精度。然而，应用地理加权回归克里格时，不同属性值在空间插值过程选用的辅助变量不同，即使同一属性值不同土壤类型或区域所选择的辅助变量也不完全相同，因此有必要在不同区域和不同土壤类型开展空间插值方法的研究。

东北黑土区是我国重要商品粮生产基地，其粮食产量关系到国家粮食安全问题，合理和精准模拟土壤养分空间分布格局是区域土壤养分调控的关键。因此，为了提高黑土区AP空间分布的模拟精度，本研究尝试将回归克里格方法和地理加权回归克里格方法应用于AP空间模拟精度研究中，综合考虑景观格局、区域尺度、采样方法等信息，选取东北黑土区两个典型黑土小流域作为研究对象，比较不同空间插值方法(反距离权重法、径向基函数法、普通克里格法、协克里格法、多元线性回归模型、回归克里格法、地理加权回归模型和地理加权回归克里格法)对黑土区AP的空间预测精度，确定黑土区AP的最优插值方法，旨在为研究区合理施肥和农业面源污染防控等提供理论依据和技术支撑。

1 研究区概况

光荣小流域位于黑龙江省海伦县光荣村(47.21°—47.23°N,126.50°—126.51°E)，流域总面积1.86 km2，地形为典型的漫川漫岗(高程186—242 m)，土壤类型以典型黑土为主，当前以大豆玉米轮作为主，一年一熟制，其基本理化形状见表1。该地区以温带大陆性季风气候为主，冬季寒冷干燥，夏季温热多雨，年均降水530 mm，主要集中在6—8月(65%)，年均气温1.5℃，年均有效积温2 450℃，年均日照时数为2 600～2 800 h，冻融期长达130 d[1](图1，表1)。

海沟河小流域位于黑龙江省哈尔滨市阿城区东部，地理坐标为126°55′45″—127°10′05″E,45°34′8″—45°40′50″N，属于典型的“林地—旱田—水田”景观格局，总面积约119.76 km2。流域东部为高山丘陵区，主要植被类型为原始的针阔混交林，土壤类型为森林暗棕壤；中部为典型的漫川漫岗区，主要耕作植物为玉米，土壤类型为典型黑土；西部地区位于流域下游区域，地势相对平坦，以水稻种植为主，土壤类型为黑土发育的水稻土。海沟河小流域以温带大陆性季风气候为主，年均降水600～800 mm，主要集中在4—10月，占总降雨量的90%，年均气温为3.9°C(图1，表1)[3,17]。

2 材料与方法

2.1 样点采集与测定

本研究分别于2012年和2014年秋季在光荣和海沟河小流域采用随机和网格设计法布点，兼顾不同土地利用方式，五点法采集0—20 cm耕层土壤样品各120个(图1)，同时记录采样点基本信息，包括经纬度坐标、高程、坡位、耕作方式、前茬植被等。土样经室内风干、研磨、过筛等处理，采用0.5 mol的NaHCO3浸提—钼蓝比色法测定AP含量[20]。

图1 研究区地理位置及采样点分布

2.2 环境变量提取

环境变量能够直接或间接地反映地球生物化学循环过程，如土壤发生过程、地表径流、淋溶、植被分布等进而影响土壤养分空间分布特征[11-12]。因此，环境变量(地形因子和遥感数据)常被用做辅助变量参与土壤养分空间异质性分析以提高模拟精度。本研究共提取30个环境变量通过相关性分析和回归分析揭示其与AP的关系，环境变量提取方法见表2[21]。

表1 研究区土壤基本属性及样点参数

表2 辅助变量及其提取方法

2.3 验证方法

为比较不同空间插值方法模拟精度，通过ArcGIS 10软件子集要素模块随机自动抽取训练样本集和验证样本集，样本数量比为4∶1。运用验证集样点的实测值与预测值进行精度评价，常用评价指标包括平均误差(ME)、平均相对误差(MAE)、均方根误差(RMSE)和相对提高度(RI)，相应的计算公式如下：

(1)

(2)

(3)

RI=(RMSER-RMSEE)/RMSER×100%

(4)

式中:n为验证点数量;z(xi,yi)和z*(xi,yi)分别为第i个验证点的实测值与预测值。ME是插值无偏性的量度,其结果越接近0说明结果越无偏;MAE和RMSE是插值精度的量度,其值越小说明插值方法越精确;RI用来两两对比不同插值方法精度的提高程度。

3 结果与分析

3.1 土壤速效磷描述性分析

对随机抽取的训练集样本AP数据进行经典统计学分析，结果显示光荣小流域AP含量的范围为2.56～50.01 mg/kg，平均值为14.51 mg/kg，偏度和峰度分别为1.86，3.85，结合Kolmogorov-Smirnov test验证法进行非参数检验可知研究区AP数据不符合正态分布，需进行数据转化。变异系数为51%，表明研究区AP的空间变异程度为中等(表3)。

海沟河小流域AP含量的范围为6.44～162.56 mg/kg，平均值为49.48 mg/kg，较光荣小流域AP含量高71%，偏度和峰度1.49，3.12，未通过K-S检验，需进行数据转化。AP变异系数为57%，属中等程度空间变异(表3)。

表3 研究区AP描述性统计特征值

3.2 土壤速效磷与辅助变量的相关性分析

光荣小流域AP与高程呈显著正相关关系(r=0.247)；与坡度、地形起伏度、地表粗糙度、地表切割深度呈极显著负相关关系；与交通距离和居民点距离呈极显著负相关关系。由于相关地形因子之间存在多重共线性，为减少数据冗余，提高精度，本研究对与AP显着相关的环境变量进行主成分分析，得到第一主成分(地形因子)和第二主成分(人类活动因子)，均与AP呈极显著负相关，可以作为辅助变量参与建模(表4)。海沟河流域AP与高程呈显著负相关关系(r=-0.257)；与含铁矿物指数呈显著负相关关系(r=0.221)；与其他环境变量并未达到显著相关(表4)。

表4 研究区AP与辅助变量相关性分析

3.3 土壤速效磷与辅助变量的回归分析

多元线性逐步回归分析能够保证与AP显着相关的环境变量进入回归模型的同时去除各环境变量之间的共线性。结果表明，坡度和交通距离是光荣小流域AP空间分布规律的最佳解释变量，引入主成分分析得到的地形因子和人类活动因子参与回归分析可提高回归方程模拟精度。高程和含铁矿物指数是海沟小流域AP空间分布规律的最佳解释变量，各变量之间不存在多重共线性问题(表5)。

表5 研究区AP多元线性回归方程及其相关参数

为了便于对不同插值方法模拟精度进行比较，在进行地理加权回归模型建模过程中同样选取上述因子作为环境变量(表6)。

表6 研究区AP地理加权回归模型参数

3.4 土壤速效磷的空间变异特征分析

通过GS+9.0软件对AP进行地统计分析揭示其空间变异特征，通常当C0/(C0+C)<25%时，变量具有较强的空间自相关性，且主要受地形等结构性因子影响；当25%75%时，变量空间自相关性较弱，主要受施肥等随机因子的影响。本结果表明光荣小流域AP的变程为274 m，可作为该区域AP空间分布特征研究的参考采样距离；C0/(C0+C)为32.5%，表明AP主要受结构因子影响，同时在一定程度上受随机因子影响，处于中等程度空间自相关水平。而海沟河流域C0/(C0+C)为7.2%，表明海沟河流域土壤AP空间自相关性较强，主要受结构性因子影响，变程为1 650 m(表7)。

表7 研究区AP地统计学模型和参数

3.5 土壤速效磷的空间分布特征

不同空间插值方法所预测的光荣小流域AP空间分布格局整体趋势基本一致(图2)。AP含量在流域上游居民点和道路附近含量较高，并沿水系方向从上游向流域出口处递减，与研究区高程变化情况较吻合。从制图效果来看，反距离权重法和径向基函数法存在明显的斑块效应，但其预测值范围更接近实测值变化范围。多元线性回归模型虽然能够反映相似的AP空间分布格局，但是其预测值范围存在明显的低估现象，插值结果并不理想，相比之下，考虑了空间非平稳性的地理加权回归模型制图效果与模拟精度均高于多元线性回归模型。普通克里格、回归克里格和地理加权回归克里格插值范围及空间预测格局相近，其中地理加权回归克里格得到的预测值范围更接近于实测值，同时考虑了环境变量的回归克里格和地理加权回归克里格的高低值界线趋于模糊化，地理加权回归克里格还能够揭示流域内AP的局部分布特征。

海沟河流域AP低值区出现在流域东部林地种植区，在中部道路两侧出现明显的条带状富集区，并沿水系方向在流域出口处达到最大值，其空间分布特征与高程变化及土地利用方式密切相关。从制图效果看(图3)，反距离权重法和径向基函数法存在明显的斑块效应，但其预测值范围更接近实测值变化范围。多元线性回归模型和地理加权回归模型出现明显的低估现象。回归克里格和地理加权回归克里格模拟的AP空间分布情况更加精细化，能够反映更多的空间信息。

3.6 模拟精度比较

以普通克里格法作为参照标准对比不同空间插值方法的评价参数，结果表明光荣小流域反距离权重法和径向基函数法的各项验证指标均低于普通克里格。多元线性回归模型和地理加权回归模型得到的均方根误差较普通克里格方法分别降低了12%，2.1%，而回归克里格和地理加权回归克里格的各项验证指标均优于普通克里格，模拟精度分别提高了4.4%，4.6%。此外，引入主成分分析结果得到的插值方法评价参数均有一定程度的提高，其中引入主成分分析的回归克里格方法模拟精度较普通克里格和回归克里格分别提高5.7%，4.5%，引入主成分分析的地理加权回归克里格方法模拟精度较普通克里格和地理加权回归克里格分别提高了6.9%，2.4%(表8)。

在海沟河流域反距离权重法、径向基函数法和协克里格的各项验证指标均低于普通克里格法；而回归克里格和地理加权回归克里格的各项验证指标均优于普通克里格法，模拟精度分别提高了0.7%，1.7%；地理加权回归模型较多元线性回归模型的模拟精度提高了3.6%(表8)。综上所述，考虑了环境变量和残差空间自相关的回归克里格和地理加权回归克里格能够提高AP的空间预测精度，地理加权回归克里格法最优。

注：A—L分别是反距离权重法、径向基函数法、普通克里格、协同克里格、多元线性回归模型、回归克里格、地理加权回归模型、地理加权回归克里格、多元线性回归模型(PCA)、回归克里格(PCA)、地理加权回归模型(PCA)、地理加权回归克里格(PCA)。

注：A—H分别是反距离权重法、径向基函数法、普通克里格、协同克里格、多元线性回归模型、回归克里格、地理加权回归模型、地理加权回归克里格。

表8 不同空间插值方法精度评价参数

4 讨 论

通过对比AP空间分布图和模拟精度评价指标，我们发现反距离权重法和径向基函数法制图效果较差，存在明显的斑块效应，同时，两种方法的各项验证指标均处于最低水平，研究表明，反距离权重法仅可在样点密度大且地形起伏较小的情况下适用，最好用于较小区域的栅格化，不宜在较大范围内进行插值[22]，径向基函数多适用于样本数目较多且变化平缓的数据[23]，Shen等[21]比较不同空间插值方法对黑土区土壤全磷模拟精度的影响，研究表明当样点空间分布不规则或均匀度较低时径向基函数会产生精确的插值结果。因此，反距离权重法和径向基函数法需在特殊情况下才可能应用并提高黑土区小流域土壤AP的空间模拟精度。

普通克里格方法利用采样点数据和半方差函数的结构性，对未采样点的区域化变量进行最优无偏估值，已被广泛应用于土壤养分空间异质性分析研究[1]。本研究中，所选择研究区样点均匀度或样点密度均较高，满足普通克里格插值方法需要，其模拟精度仅次于回归克里格和地理加权回归克里格方法。

地理加权回归模型是一种局部回归模型，可用于处理目标变量与解释变量之间的空间非平稳回归系数[14-15]。本研究中，地理加权回归模型模拟精度低于普通克里格方法，这可能是由于经逐步线性回归剔除后可用于建模的环境变量较少，回归模型拟合精度较低。Wang等[13]比较了地理加权回归模型对土壤有机质和全氮模拟精度的影响，结果表明地理加权回归模型能够有效提高土壤有机质的模拟精度，但是由于环境变量数据不详尽且未能解释所有相关变量之间的自相关和互相关而导致该方法对全氮模拟精度较差[24]。本研究已提取相关环境变量达30个，但并未进一步提高回归模型精度，这与Shen等[21]在黑土区土壤全磷空间模拟精度研究结果一致。此外，比较不同研究区回归分析结果，并未发现共同的环境变量用于AP建模。因此，有关如何进一步提高黑土区土壤磷素与环境变量之间的回归模型拟合精度的问题需进一步探讨。研究表明，样点密度会影响空间插值方法对土壤属性法的模拟精度[19]，但本研究发现，当样点均匀度相似且高于45%时，样点密度和采样距离对AP模拟精度无显著影响。

回归克里格和地理加权回归克里格已经逐渐成为土壤属性空间预测的优秀空间统计方法，大量研究表明辅助变量的引入和残差空间自相关的结合能够有效提高土壤属性的空间预测精度[9-10,16-19]，本研究结果表明，在不同尺度和景观格局下，回归克里格和地理加权回归克里格各项验证指标均优于普通克里格插值方法，插值结果也能够反映更多的空间信息，其中考虑了空间非平稳性的地理加权回归克里格方法可以作为黑土区小流域土壤AP空间预测的最优插值方法。此外，我们发现引入主成分分析结果能够有效提高回归克里格方法和地理加权回归克里格方法的模拟精度，因此，当辅助变量充足且变量间存在多重共线性时，可以引入主成分分析以进一步提高空间插值方法的模拟精度。

5 结 论

通过引入地形等环境因子建立的回归克里格、地理加权回归克里格方法较传统空间插值方法能在一定程度上提高黑土小流域土壤AP的空间模拟精度，将主成分分析结果与上述方法相结合能够进一步提高模拟效果。其中，地理加权回归克里格方法模拟效果最优，可作为黑土小流域土壤AP空间分布模拟的主要方法。本研究选取的辅助变量以地形和植被覆盖等环境因子为主，未能充分考虑与AP密切相关的施肥、耕作措施等农田管理因子，可作为提高黑土小流域土壤AP空间模拟效果的重要方向。研究结果可为黑土小流域尺度土壤养分精准管控提供技术支撑。