土壤有机质(soil organic matter, SOM)属土壤重要组成部分^[1], 是植物必需营养元素的主要来源^[2], 虽然在土壤总重量中所占比重不大, 但其含量是衡量土壤肥力水平的一项重要指标^[3-5]。柑橘(Citrus reticulate Blanco.)产业是湖北省重要经济支柱之一, 种植面积达全省果园总面积的45%, 然而近年来由于柑橘园SOM含量总体偏低^[6], 严重抑制橘树的生长发育, 导致柑橘产量及品质均有降低^[7]。因此, 为科学指导柑橘园精准施肥, 需要建立快速准确的方法预测SOM空间分布与含量。

受多种成土因素和生态过程影响, 即便是小区域尺度内的SOM含量也呈现出不同程度的空间异质性和非平稳性^[8-11], 单纯通过野外土壤田间调查又难以满足精确获取SOM空间分布信息的实际需求。为解决以上问题, 近年来诸多学者将预测SOM含量的关注点放在能够对未采样区域待估点进行无偏最优估计^[12]的克里格法(Kriging)上。在普通克里格插值法(Ordinary Kriging, OK)的应用方面, Singh等^[13]提出在不同果树种植园中, 通过OK对SOM进行空间估计的不确定性比非空间方法低且能在未采样位置获得满意的预测结果; 李增兵等^[14]使用不同插值方法分别对历城区不同地貌类型建模, 得到OK法更适用于平原区的结果。OK方法虽能考虑SOM的空间自相关性和异质性^[15], 但仅依赖土壤样本建模集有机质数据插值, 并不能发挥辅助变量对估算SOM的优势。协同克里格插值法(Cokriging, COK)则可基于不同样点支撑^[16-17], 借助同一时空不同变量的协同区域化特征, 用样品多的辅助变量对样品少的目标变量SOM进行估值^[18]。苏晓燕等^[19]和Yang等^[20]提出, 辅助变量的个数大于目标变量且辅助变量与目标变量的相关性越强时, 使用COK越能提高预测精度。但COK的建模精度受辅助变量插值样点数量影响很大, 且交叉变异函数能否充分考虑目标变量和辅助变量间的相关信息以及是否可在空间分布模拟中维持单变量的结构性等问题尚需进一步补充与验证。因此, 本文结合以上2种传统克里格插值法的特点, 提出建立融合辅助变量协同相关性的改进OK模型(improved OK models with cooperative correlation of auxiliary variables, CCOK)。该模型在保持目标变量与辅助变量样点支撑一致且数目有限的情况下融入COK法中的协同区域化思想, 假设目标变量及其辅助变量受到同样的区域化现象或空间过程影响, 各变量不但存在自相关性, 而且与其他变量之间存在交互相关性, 尝试将传统OK法插值数据从单一变量发展到两个相关性极显著的变量函数, 既兼顾了目标变量及其辅助变量的空间结构相似性^[21], 又可借助辅助变量提高目标变量实测值与预测值的拟合精度。

本研究以湖南省宜都市红花套镇为研究区, 基于协同区域化理论, 先后利用建模集辅助变量、全局辅助变量对该地SOM建立COK模型(COK₁、COK₂), 探索性构造CCOK模型(CCOK₁、CCOK₂), 并将插值结果与OK模型作对比, 分别通过探讨纳入辅助变量、改变辅助变量插值数量以及结合辅助变量协同相关性对SOM含量预测的影响, 获得可显著提高SOM拟合精度的空间插值方法, 以期为监测土壤农田肥力投入情况及增加农业经济效益提供参考。

1 材料与方法 1.1 研究区概况

湖北省宜都市红花套镇处于鄂西山区向江汉平原过渡的地带(30°26′45″~30°35′59″N, 111°14′13″~ 111°26′38″E), 总面积为149 km²。境内西北地势高, 以丘陵山地为主; 东南地势低, 以冲积平原为主, 土地肥沃。该镇属亚热带湿润季风气候, 年均气温16.8 ℃, 年均降雨量1 200 mm, 全年无霜期275 d, 水文、生物、矿产资源丰富, 生态良好, 适宜水稻(Oryza sativa L.)、棉花(Gossypium spp.)和柑橘生长。其中柑橘种植面积28 km², 年产量达10万余t, 素有“万亩桔乡”的美称。本研究区内高程范围为42~112 m, 园地、林地、耕地和其他土地利用类型分别占总面积的72.46%、0.59%、0.13%和26.82%(图 1)。主要土壤类型为黄棕壤、水稻土, 也有少量潮土、紫色土和石灰土。

1.2 数据来源与预处理

在已有红花套镇1:10 000比例尺地形图、第2次土壤普查相关资料的基础上, 综合地形地貌、土地利用类型、土壤类型以及植被覆盖等相关因素预先设置研究区范围。当该地区柑橘收获且未施用肥料之时, 进行实地调查并结合实际因素对样点位置做适当调整。使用差分式全球定位系统(DGPS)记录经纬度及高程, 采集样点处0~20 cm表层土。采样在平原的园地、耕地区样点内以梅花形法进行, 在丘陵的林地区样点内以蛇形法进行, 共获得329个土壤样本。

地形平缓的小区域内SOM含量分布主要受其他土壤属性影响^[20], 因此根据研究区主要土壤养分的分布情况、主要耕种作物的土壤环境营养元素需求特性以及当地农业生产等方面的实际情况, 获取表 1中土壤属性数据作为影响因子参与研究, 获取方法均参考文献[22]。

使用SPSS 19.0中3倍标准差法检测样本数据中并无需要剔除的异常值, 再通过ArcGIS 10.2中子要素集功能对样点数据集进行均匀随机划分, 最终获得建模集数据263个(80%), 验证集数据66个(20%)。

1.3 研究方法 1.3.1 空间自相关分析

空间自相关分析是用来反映空间邻近区域某变量属性值的相似程度和空间分布的聚集特征, 全局Moran’s I是其中最常用的检验统计量^[23-25]。

$ I = \frac{{\sum\nolimits_{i = 1}^n {\sum\nolimits_{j \ne i}^n {{w_{ij}}\left( {{x_i} - \bar x} \right) \times \left( {{x_j} - \bar x} \right)} } }}{{{S^2}\sum\nolimits_{i = 1}^n {\sum\nolimits_{j \ne i}^n {{w_{ij}}} } }}\;\;\;\left( {i \ne j} \right) $

(1)

式中: I代表Moran’s I指数值, n代表变量x的样本数, x_i、x_j代表位置i和j的样本实测值, S²代表其方差, $ {\bar x}$代表平均值, w_ij代表对称二项分布空间权重矩阵元素。Moran’s I的取值范围为–1~1, I > 0表示空间正相关性, 其值越大, 空间相关性越明显; I < 0表示空间负相关性, 其值越小, 表示空间差异越大; 否则, I=0, 空间呈随机性。

1.3.2 OK模型

OK是单个变量的局部线性最优无偏估计方法, 也是Kriging中最基础的方法之一。它能够根据给定区域内相邻实测点的加权平均值, 为变量在未测量位置提供预测值。该方法通过考虑实测点表示的固有随机函数来简要描述^[26]。

$ Z_{{\rm{OK}}}^*\left( {{x_0}} \right) = \sum\limits_{i = 1}^n {{\lambda _i} \times Z\left( {{x_i}} \right)} $

(2)

式中: $ Z_{{\rm{OK}}}^*\left( {{x_0}} \right)$代表使用OK插值方法得到的待估点x₀处的预测值, λ_i代表分配给每个实测样点的权重且$\sum {{\lambda _i} = 1} $满足无偏性条件, Z(x_i)代表实测值, n代表参与预测的实测样点数量。

1.3.3 建模COK模型与全局COK模型

在统计意义及空间位置上均具有某种程度相关性, 并且定义于同一空间域中的区域化变量, 即具有协同区域化特征^[16]。一般情况下, 结合了辅助变量信息的协同克里格能够比普通克里格进行更多的预测。本文提出的建模COK模型与全局COK模型, 其目标变量均为建模集SOM数据, 区别在于建模COK模型选取了建模集中辅助变量参与插值, 而全局COK模型则选取了全局数据的辅助变量参与插值, 以此比较改变辅助变量插值数量对SOM预测精度的影响。

$ Z_{{\rm{COK}}1}^*\left( {{x_0}} \right) = \sum\limits_{a = 1}^m {{\lambda _a} \times Z\left( {{x_a}} \right) + \sum\limits_{b = 1}^p {{\lambda _b} \times U\left( {{x_b}} \right)} } $

(3)

$ Z_{{\rm{COK2}}}^*\left( {{x_0}} \right) = \sum\limits_{c = 1}^m {{\lambda _c} \times Z\left( {{x_c}} \right) + \sum\limits_{d = 1}^q {{\lambda _d} \times U\left( {{x_d}} \right)} } $

(4)

式中:$Z_{{\rm{COK1}}}^*\left( {{x_0}} \right) $、$ Z_{{\rm{COK2}}}^*\left( {{x_0}} \right)$分别代表使用建模集、全局辅助变量协同建模集主变量进行COK方法插值得到的待估点x₀处的预测值, λ_a、λ_c分别代表分配给建模COK模型、全局COK模型主变量Z的实测样点权重且$ \sum {{\lambda _a} = \sum {{\lambda _c} = } 1} $, λ_b、λ_d分别代表分配给建模COK模型、全局COK模型辅助变量U的实测样点权重且$\sum {{\lambda _b} = \sum {{\lambda _d} = } 0} $, Z(x_a)、Z(x_c)均代表建模集中主变量Z的实测值, U(x_b)、U(x_d)分别代表建模集和全局数据集中辅助变量U的实测值, m代表参与预测的主变量实测样点数量, p、q分别代表参与预测的建模集和全局辅助变量实测样点数量。

1.3.4 CCOK模型

CCOK模型是传统OK插值方法与COK插值原理相结合的改进模型, 它考虑了协同相关性引入辅助变量, 是对代数量化后的主变量与辅助变量进行OK插值的方法。为充分说明与主变量存在极显著相关的辅助变量对提高插值精度具有重要意义, 同时验证对不同单位的变量不设置数据标准化预处理也不会影响最终预测效果, 特列以下公式完成CCOK建模。

$ Z_{{\rm{CCOK1}}}^*\left( {{x_0}} \right) = U_{{\rm{OK}}}^*\left( {{x_0}} \right) \times \left( {\frac{Z}{U}} \right)_{_{{\rm{CCOK1}}}}^*\left( {{x_0}} \right) $

(5)

$ Z_{{\rm{CCOK2}}}^*\left( {{x_0}} \right) = U_{{\rm{OK}}}^*\left( {{x_0}} \right) \times \left( {\frac{U}{Z}} \right)_{_{{\rm{CCOK2}}}}^*\left( {{x_0}} \right) $

(6)

式中: $Z_{{\rm{CCOK1}}}^*\left( {{x_0}} \right) $、$ Z_{{\rm{CCOK2}}}^*\left( {{x_0}} \right)$分别代表经插值计算得到的待估点x₀处主变量Z的预测值, $U_{{\rm{OK}}}^*\left( {{x_0}} \right) $代表使用OK插值方法得到的待估点x₀处辅助变量U的预测值, $ \left( {\frac{Z}{U}} \right)_{_{{\rm{CCOK1}}}}^*\left( {{x_0}} \right)$、$ \left( {\frac{U}{Z}} \right)_{_{{\rm{CCOK2}}}}^*\left( {{x_0}} \right)$则分别代表使用CCOK₁、CCOK₂插值方法得到的待估点x₀处$ {\frac{Z}{U}}$、$ {\frac{U}{Z}}$的预测值。其中:

$ \left( {\frac{Z}{U}} \right)_{_{{\rm{CCOK1}}}}^*\left( {{x_0}} \right) = \sum\limits_{j = 1}^n {{\lambda _j} \times } \frac{Z}{U}\left( {{x_j}} \right) $

(7)

$ \left( {\frac{U}{Z}} \right)_{_{{\rm{CCOK2}}}}^*\left( {{x_0}} \right) = \sum\limits_{k = 1}^n {{\lambda _k} \times } \frac{U}{Z}\left( {{x_k}} \right) $

(8)

式中: Z(x_j)、Z(x_k)代表样点处主变量Z的实测值且$Z\left( {{x_j}} \right) = Z\left( {{x_k}} \right) $, U(x_j)、U(x_k)代表样点处与主变量极显著相关的辅助变量实测值且$U\left( {{x_j}} \right) = U\left( {{x_k}} \right) $, n代表参与预测的实测样点数量。为满足插值数据具有无偏性, 模型需符合预测值与实测值间偏差的数学期望为0的条件, 即$ E\left[ {{{\frac{Z}{U}}^*}\left( {{x_0}} \right)} \right] = E\left[ {\frac{Z}{U}\left( {{x_0}} \right)} \right]$和$E\left[ {{{\frac{U}{Z}}^*}\left( {{x_0}} \right)} \right] = E\left[ {\frac{U}{Z}\left( {{x_0}} \right)} \right] $, 因此要求$\mathop \sum \nolimits {\lambda _j} = \mathop \sum \nolimits {\lambda _k} = 1 $。同时, $\frac{Z}{U} $、$\frac{U}{Z} $的整体平均值还要在局部预测变程范围内维持基本不变, 才能满足插值模型的二阶平稳性。

1.4 预测模型精度验证

$ {\rm{MAE}} = \frac{1}{n}\mathop \sum \limits_{i = 1}^n \left\{ {\left| {Z\left( {{x_i}} \right) - {Z^*}\left( {{x_i}} \right)} \right|} \right\} $

(9)

$ {\rm{RMSE}} = \sqrt {\frac{1}{n}\mathop \sum \limits_{i = 1}^n {{\left\{ {Z\left( {{x_i}} \right) - {Z^*}\left( {{x_i}} \right)} \right\}}^2}} $

(10)

$ {\rm{RI}} = \frac{{{\rm{RMS}}{{\rm{E}}_{x1}} - {\rm{RMS}}{{\rm{E}}_{x2}}}}{{{\rm{RMS}}{{\rm{E}}_{x1}}}} \times 100{\rm{\% }} $

(11)

式中: n代表样本数, Z(x_i)代表样本的实测值, Z^*(x_i)代表样本的预测值, RMSE_x1代表使用OK插值得到的RMSE值, RMSE_x2代表使用COK或CCOK插值得到的RMSE值。由于离差被绝对值化, 不会出现正负相抵消的情况, 因而平均绝对误差MAE能更好地反映预测值误差的实际情况; 均方根误差RMSE是插值精度的量度, 其值越小则说明插值方法越精确; 相对精度改进值RI为正值代表COK或CCOK方法较OK方法的预测精度更高, 值越大说明提高越多, RI为负值则表示COK或CCOK方法预测精度低于OK方法。

1.5 软件平台

利用SPSS 19.0完成原始数据非空间分析(剔除异常值、基本特征统计和最相关辅助变量选取), 通过GeoDa证明变量间存在的区域协同化特征(全局空间相关性分析), 采用ArcGIS 10.2进行空间分析(栅格计算)、地统计学分析(随机划分数据集、构建最优半方差模型、完成OK、COK与CCOK插值)。

2 结果与分析 2.1 研究区样点SOM含量数据的特征统计

表 2中全部数据集的SOM含量均值为16.97 g·kg^-1, 范围在6.07~29.18 g·kg^-1, 最大值和最小值间差异显著且变异系数为24.87%, 均属中等变异, 说明SOM在该研究区范围内具有一定变异性, 适合空间局部估计。建模集中的有机质数据(偏度=0.07, 峰度=-0.08)通过K-S(Kolmogorov-Smirnov Z)检验, 符合正态分布, 可直接用于插值建模。

2.2 研究区样点SOM最相关辅助变量的选取

表 3显示了该研究区中与SOM相关性由高到低依次排序的前7种影响因子及其Pearson相关系数。有机质与碱解氮、有效铁、全氮和全磷表现为正相关, 与交换性镁、阳离子交换量和速效钾表现为负相关。

其中, SOM与碱解氮间的相关系数为0.815, 呈极显著正相关, 相关系数高于其他因子。这是因为SOM主要由C、H、O、N等元素构成, 碱解氮中则包含无机的矿物态氮和部分有机质中易分解的、比较简单的有机态氮^[27], 氮素是二者重要组成成分中的共同元素。在受同一区域下水热条件、成土特征及微生物反复作用下, SOM与碱解氮之间的耦合关系逐渐趋于动态平衡^[28], 理论上有机质含量丰富, 熟化程度高, 碱解氮的含量就高, 反之则含量降低。但实际上, 碱解氮作为一种速效养分, 其含量及空间分布还与人为投入氮素化肥数量等外界因素有一定关联^[29], 可以说碱解氮中的氮素主要来源于有机质, 但其含量、组成、形态和有效性又不完全等同于有机质。此外, 相关研究表明协同变量与自变量间的相关系数大于0.5时更能在预测结果中体现协同区域化优势^[30], SOM与碱解氮间的相关性符合该条件并远高于其他因子。因此, 研究选择碱解氮作为预测SOM的最相关辅助变量。

2.3 研究区样点SOM、碱解氮的空间自相关与交互相关分析

SOM、碱解氮的Moran’s I指数分别为0.11、0.08, 具有一定正空间自相关性, 当对两者建立简单函数关系后指数值分别提高到0.19和0.15(表 4)。同时, 综合比较代表标准差的Z得分与代表概率的P值发现, 属性的空间聚类置信度也从95%升为99%, 说明SOM与碱解氮在空间分布上表现出良好的交互相关性, 促进了其函数关系呈现出极显著的空间自相关性(P < 0.01), 这是将碱解氮作为辅助变量建立预测SOM的CCOK模型的关键依据。

2.4 不同克里格插值模型对研究区SOM半方差模型的选择与分析

表 5为使用不同克里格模型对SOM预测时选择的最优半方差模型及其显示参数。半方差函数预测误差的作用是为Kriging插值判断拟合的最优模型。遵循平均值越趋近于零越无偏、均方根误差RMSE越小越精准的原则, OK、COK₁、COK₂、CCOK₁和CCOK₂分别选择高斯模型、球状模型、指数模型、指数模型和指数模型作为插值的半方差模型。模型参数中, OK和COK₂的块金值较大, 说明存在由于采样或分析引起的各种正基底效应^[31]; CCOK₁的变程为3 627.41 m, 是由于为符合OK插值数据需正态分布的条件, 对有机质/碱解氮数据进行了对数转换, 转换后的数值为-2.15~-0.71, 极差仅为1.44, 因此在研究区的大面积连续范围内均表现出空间分布的相近性。而其他方法的变程较小, 表示SOM含量在短距离内具有空间连续性^[32]。就基台效应而言, COK₁、COK₂具有强烈的空间相关性, OK、CCOK₁、CCOK₂表现出中等空间相关性, 说明SOM的空间异质性受结构性因素影响的比重更大。图 2是对应于表 5的SOM半方差函数图。分析图像中绘制曲线通过表示半变异函数平滑变化的平均化数据情况, 可更加直观地反映各插值模型的拟合效果。其中, CCOK₁的半方差函数拟合效果相较于其他方法表现微弱甚至难以体现, 但因其变程范围远大于步长距离, 变程内的区域化变量仍存在空间相关状态。所以, 需将SOM半方差模型参数及函数图综合分析, 才能完成使用不同方法对SOM预测的合理预判。

2.5 不同克里格插值模型对研究区SOM空间分布预测的差异比较

对比图 3研究区的SOM空间分布预测图可知, 使用5种模型取得的预测结果分布情况整体相似, SOM含量范围为7.38~29.03 g·kg^-1。按照 < 5 g·kg^-1为极低、5~10 g·kg^-1为低、10~15 g·kg^-1为偏低、15~30 g·kg^-1为适宜、> 30 g·kg^-1为丰富的标准^[33]划分, 研究区大部有机质含量适宜, 小面积区域含量低, 无含量极低或丰富地区, 呈现出“西北和东南高, 中部高低错落”的分布趋势, 这与该研究区的土地利用类型分布格局有关。西北部为与小面积林地衔接的柑橘种植带, 受人为干扰少, 凋落物富集, SOM分解速度慢便于长期积累。东南部则为柑橘的集中种植区, 为促进果树生长和发育、提高果树产量及品质, 定期的人为施肥和排水增加了SOM含量^[34]。而研究区中部的柑橘园地与零星耕地、其他土地利用类型交叉分布, 并以土壤为媒介发生长期物质交换, 产生相互作用, 导致了有机质含量分布不均且高低差异显著, 因此是柑橘园中调控施肥用量的重点位置。

以OK预测结果为对照, 添加了碱解氮作辅助变量的COK₁、COK₂模型插值获得的有机质分布更加破碎, 空间镶嵌结构更为复杂^[35], 一部分原因是OK模型本身具有平滑效应使得空间预测表现连续, 另一部分原因是COK₁、COK₂模型通过目标信息及其辅助信息间建立的交叉变异函数促使生成更加详尽的空间信息表达, 且COK₂较COK₁的辅助变量样点分布更密, 因此能够观察到有机质含量的明显变化, 并体现出局部变异细节。

与建立前3种模型得到的SOM预测图呈点状放射相比, CCOK₁、CCOK₂模型的插值结果为连续片状放射。这是由于OK、COK₁和COK₂模型以建模集有机质为惟一或主要数据源, 在数字制图中继续维持SOM的空间结构特征。而CCOK₁、CCOK₂方法的建模数据源选择为有机质与碱解氮之间的代数量化数据, 虽然经公式转回为对SOM的预测, 但同时使用样点支撑一致的变量参与插值, 使得预测制图能够充分体现出辅助变量对目标变量的协同相关性, 有机质含量高者更高、低者更低, 既突出了差异信息, 也更符合研究区土地利用类型分布的实际情况。

2.6 不同克里格插值模型对研究区SOM含量预测的精度验证

使用验证集SOM数据进行精度验证, 经5种模型插值得到的实测值与预测值间CCOK₁的相关系数r最高为0.70, CCOK₂次之, 与相关系数r仅为0.10的OK模型相比预测效果提升显著, CCOK₁和CCOK₂的平均绝对误差MAE最低且其预测值趋势线与1:1参考线的拟合程度最优。以OK插值结果为参照, 除COK₁无显著变化, COK₂、CCOK₁和CCOK₂的预测精度改进值RI分别提高了23.20%、15.40%和14.78%(图 4)。

其中, OK和COK₁两者间各精度指标相近, 充分证明当作为目标变量的SOM与作为辅助变量的土壤碱解氮样点支撑完全相同时, COK方差即为OK方差^[16], 是否添加辅助变量对预测目标变量无较大意义。COK₂则是在已有COK₁模型的基础上, 将建模集辅助变量改为全局数据集辅助变量参与插值, 通过增大样点数量提升了辅助预测效果。CCOK₁和CCOK₂的建模原理相同, 都在兼顾目标变量及其辅助变量空间结构相似性的同时, 最大程度地发挥了辅助变量的协同相关作用, 由二者并无较大差异的插值结果也能间接证明CCOK模型对未经标准化处理的建模数据仍然适用。综合r、MAE、趋势线与1:1参考线拟合程度、RMSE、RI这5项指标分析, SOM空间分布预测精度由高到低依次为CCOK₁≈CCOK₂ > COK₂ > COK₁≈OK, CCOK模型在很大程度上提升了SOM含量的拟合效果。本文提出的建立融合辅助变量协同相关性的改进OK模型进行SOM空间分布预测可缩小实测值与预测值间差距, 为研究区柑橘种植提供精准施肥建议。

3 结论

1) 本柑橘种植区中, 使用几种克里格插值模型预测得到的SOM含量均为7.38~29.03 g·kg^-1, 整体呈现“西北和东南高, 中部高低错落”的分布趋势。

2) 以OK模型作对照, 与辅助变量间建立交叉变异函数的SOM的空间变异性更容易被半方差模型识别, 使得COK₁、COK₂因模型参数存在强烈空间自相关在SOM空间分布上表现破碎, 并能维持单变量空间结构性。然而, 通过改变辅助信息样点数量提高SOM预测精度, 既会增加研究成本又未充分考虑每个样点上目标变量与辅助变量间的交互关系, 可能引进更多变异, 因此使用该模型存在一定局限。

3) CCOK₁和CCOK₂模型则能在目标变量与辅助变量样点支撑一致且样点数目有限的情况下, 最大程度发挥辅助变量的协同相关性, 其优势不仅体现在SOM的空间分布模拟结果最符合该研究区土地利用类型的实际分布情况, 更表现在CCOK₁和CCOK₂的模型预测精度分别提高了15.40%、14.78%, 且相关系数从0.10显著升到0.70、0.69。

4) 此研究已证明, 使用融合辅助变量协同相关性的改进OK模型可实现对研究区SOM的精准预测, 能够为区域柑橘种植区施肥提供合理建议, 而其适用范围及条件将是下一步讨论的重点。