一种顾及土壤有机碳密度空间异质性的克里格插值模型

摘要

本发明公开了一种顾及土壤有机碳密度空间异质性的克里格插值模型，结合土地利用信息降低空间异质性和空间异常值对普通克里格法的影响，且根据不同的趋势项函数，所述克里格插值模型包括虚拟变量回归克里格法、均值中心化克里格法和中位数中心化克里格法。并利用模型实现更高精度的SOCD预测和制图，最后根据模型评价指标值选择最佳预测方法及其预测结果。本发明的模型可靠性、预测合理性、预测精度和决定系数R

说明书

技术领域

本发明属于土壤质量评价技术领域，涉及一种土壤有机碳密度(SOCD)空间分布预测和制图方法，特别涉及一种顾及土壤有机碳密度空间异质性的土地利用-克里格插值模型。

背景技术

土壤有机碳密度(SOCD)是评价土壤质量的重要指标，深刻影响着土壤的物理化学性质。准确描述区域土壤有机碳密度的空间分布特征，完成SOCD空间分布的精确制图，对实现精确农业、指导农业生产和提高农业生产效率有重要意义。普通克里格法(OK)作为线性无偏最优的空间插值方法，可以根据采样点属性值的空间统计规律，得到区域化变量在全域的空间分布特征。但是，由于外界因素的综合干扰，使得区域化变量SOCD在全域内存在空间异质性，甚至出现空间异常值，违背了OK的平稳性前提假设，导致模型的可靠性和预测精度下降。

现有研究通常引入辅助环境变量提升模型对SOCD空间变异的解释能力，主要包括气候因素(张慧智等，2008年)、植被覆盖状况(邱乐丰等，2010年)、土地利用类型(顾成军等，2014年)、地形地貌(杨顺华等，2015年)、和土壤类型(Ku Wang等，2013年)等；近年来，学者们还引入了模糊坡位信息(Qin等，2012年)、高时间分辨率遥感数据(刘峰等，2011年)和近地高光谱信息(刘艳芳等，2017年)用于数字土壤制图。但是，在人类活动对土壤属性影响较大的平原农田地区，由于地形起伏度小、区域内自然植被已经基本被农作物取代、区域范围内的气候降水差异较小等原因，使用传统自然要素作为辅助变量时与土壤属性值的协同程度较低，难以提升模型预测精度。

由于土地利用类型可以由遥感影像直接解译得到、获取方便，因此结合土地利用类型提升土壤属性空间插值模型的精度成为当前平原农田地区精细数字土壤研究的重要方向。在结合土地利用信息进行土壤属性克里格插值的研究中，有学者根据土地利用类型进行空间分区，分别实现子区域空间插值后将结果合并，但这种方法多用于土地利用方式在空间上大面积连续的情况(Wu Qian等，2017年)；此外还有学者根据土地利用类型属于名义变量的特性，使用均值中心化与克里格结合的方法(顾成军等，2014年)或虚拟变量回归克里格法(Wen Wen等，2015年)进行土壤属性值空间预测，但是目前还没有一种模型可以将这两种方法统一起来，并进行两种方法预测精度的比较；此外，当校正集存在空间异常值时，这两种方法的预测精度还有待提升。

发明内容

为了解决上述技术问题，本发明提出了一种顾及空间异质性的土地利用-克里格模型，应用于受人类影响深刻的平原农田地区，目的是实现更高精度的SOCD空间分布预测及制图，指导农业生产、实现精确农业。

本发明提供的一种顾及土壤有机碳密度空间异质性的克里格插值模型，其特征在于：结合土地利用信息降低空间异质性和空间异常值对普通克里格法的影响，且根据不同的趋势项函数，所述克里格插值模型包括虚拟变量回归克里格法(DV_RK)、均值中心化克里格法(MC_OK1)和中位数中心化克里格法(MC_OK2)。

所述不同的趋势项函数是指，虚拟变量回归克里格法的趋势项是连续函数，而均值中心化克里格法和中位数中心化克里格法的趋势项是分段函数。其中，中位数中心化克里格法，是将不同类别的SOCD值减去该类别的中位数，用得到的残差进行普通克里格插值，则任一位置的SOCD值等于残差克里格插值与该点所属类别中位数之和。虚拟变量回归克里格的趋势项由SOCD值和土地利用类型利用普通最小二乘(OLS)拟合得到，而土地利用类型需要转换为虚拟变量：以水浇地地类为参照组，设其值为0，则水田地类值为1。均值中心化克里格法和中位数中心化克里格法的趋势项是仅与土地利用类型的类别有关的分段函数，其值分别等于各类别的均值/中位数。

本发明提供的一种利用顾及土壤有机碳密度空间异质性的克里格插值模型进行SOCD预测和制图的方法，其特征在于，包括以下步骤：

步骤1：原始校正集和验证集的建立；

从SOCD总样本集中随机选出一部分样点作为原始校正集，其余作为外部验证集。

步骤2：单因素方差分析；

根据单因素方差分析结果判断土地利用方式的差别是否导致SOCD值产生显著差异，如果存在显著差异，则使用土地利用-克里格插值模型可以提高预测精度。

步骤3：三种土地利用-克里格方法校正集的建立和正态性检验；

计算出校正集中不同土地利用类型SOCD的均值和中位数；以土地利用类型作为自变量(虚拟变量)，以SOCD为因变量，得到两者的回归方程；则可以得到原始校正集经过虚拟变量回归、均值中心化和中位数中心化后的残差值集合，分别作为三种方法的校正集；通过QQ图定性观察三组结合土地利用的校正集的正态分布特征、使用K-S检验定量确定三组校正集是否呈正态分布。

步骤4：SOCD残差空间预测；

将通过K-S检验的校正集输入GS+9软件，可以得到校正集的理论半变异函数模型及相关参数，将结果输入ArcGIS10.2软件地统计模块中进行普通克里格插值，可以得到全域的残差预测结果，并转为栅格图像。

步骤5：SOCD空间分布预测和制图；

使用采样同时期遥感影像，通过计算机解译得到研究区土地利用类型图，且土地利用栅格图的分辨率应与克里格插值栅格图的分辨率一致。则可以使用ArcGIS10.2中rastercalculator工具，将每个栅格按照不同地类加上相应的趋势项函数值，即可得到SOCD的空间预测结果和空间分布栅格图像。

本发明提供的一种顾及土壤有机碳密度空间异质性的克里格插值模型评价方法，其特征在于：从内部舍一交叉验证和外部验证两个方面，使用平均克里格方差(AKV)评价模型可靠程度、根据SOCD预测范围判断模型预测合理性、根据均方根误差(RMSE)判断模型的预测精度、根据R²值判断模型对SOCD总方差的解释能力(具体计算见公式(1)、(2)、(3))。综合以上模型指标评价结果，从三种土地利用-克里格法中选择预测精度最好的结果，作为该研究区SOCD空间分布预测的最终结果；

其中，n表示校正集样本数量，O_i表示SOCD第i个实测值，表示SOCD实测值的均值，P_i表示SOCD第i个预测值，表示位置i处的克里格方差。

本发明具有的有益效果是：

(1)在前人研究的基础上提出中位数中心化克里格法，并根据DV_RK、MC_OK1、MC_OK2这三种方法的共性，将它们置于土地利用-克里格模型之中。其中，三种方法的共同思想是，将SOCD的实测值拆分为趋势项和残差项，利用辅助变量消除趋势项后，得到代表随机性的残差项，可以通过普通克里格法得到全区域的残差预测值，则位于x₀处的SOCD的预测值等于自身趋势项函数值加上残差插值结果；不同的是，虚拟变量回归克里格法的趋势项是连续函数，由最小二乘法拟合得到，而均值中心化和中位数中心化克里格法的趋势项是分段函数，其趋势项函数值仅和土地利用类别值有关。

(2)三种结合土地利用的克里格法相比于普通克里格法，减弱了空间异质性对模型的影响，其模型可靠性、预测合理性、预测精度和对SOCD总方差解释能力均显著提高，相比于普通克里格法实现了更高精度的SOCD空间制图，有利于指导农业生产；中位数中心化克里格法在SOCD存在空间异常值时，是一种有效的预测方法，且预测精度高于DV_RK和MC_OK1；可以根据SOCD现有的空间分布特征，可以进行成因分析等更深层次的研究。

附图说明

图1是本发明实施例的研究区区位和采样同时期研究区土地利用遥感解译图；

图2是本发明实施例的OK、DV_RK、MC_OK1和MC_OK2四组校正集样本分布QQ图；

图3是本发明实施例的OK、DV_RK、MC_OK1和MC_OK2四种方法SOCD空间分布预测结果。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本实施例的研究区位于江汉平原，毗邻洪湖，地形起伏度小，受到人类活动的深刻影响，自然植被已基本被农作物取代，土地利用类型主要是水田、水浇地、坑塘水面建设用地，且土地利用斑块破碎度较大。在这种情况下，前人提出的方法都不完全适用。下面将通过这个案例，具体说明土地利用-克里格模型在平原农田SOCD预测及制图的应用：

(1)总样本集共计242个点，其中水田114个、水浇地128个(附图1)。从中随机选出193个点作为校正集(占总数的80％)，其余49个点作为验证集(占总数的20％)。

(2)单因素方差分析结果显示，水田和水浇地的SOCD存在显著差异(p<0.01)，即土地利用方式的差异是造成SOCD存在空间异质性的原因之一。

(3)计算出水田SOCD的均值和中位数分别为6.037和6.140kg/m²，水浇地SOCD的均值和中位数分别为4.697和4.355kg/m²，虚拟变量回归方程为：Y_SOCD＝1.341X_LU+4.697，则可以得到原始校正集和三组土地利用-克里格模型的校正集QQ图(附图2)，可以发现：结合土地利用类型的三组校正集在高值区和低值区的点群比OK校正集更贴近对角线，且MC_OK2校正集在低值区贴近效果更好。

(4)K-S检验的结果表明四组校正集均符合正态分布(p<0.05)，使用GS+9.0获取四种预测方法的理论半变异函数模型及参数，将结果输入ArcGIS10.2软件地统计模块中进行克里格插值。

(5)根据采样同时期Landsat 8遥感图像解译得研究区土地利用图，结合土地利用图和四组校正集克里格插值图，可以得到最终的SOCD空间分布预测图，见附图3(为了使得它们在目视评价上具有可比性，将四图的SOCD取值范围均设置为0～12kg/m²)。发现图2a中存在空间异常值点，而在图2b、2c和2d中部分已被消除。

本实施例分内部舍一交叉验证和外部验证，统计OK、DV_RK、MC_OK1和MC_OK2的模型评价指标值，见表1，可以发现：三种土地利用-克里格插值方法的预测效果优于OK，且MC_OK2预测效果更优，故选择附图3d作为SOCD空间分布预测最终结果。

表1：模型预测精度评价

以上实验结果表明：在研究区地形因子差异不大的平原地区，土地利用类型作为一个易获取的辅助变量，可以有效减弱空间异质性和空间异常值对模型的影响，与克里格法结合，降低模型不确定性、提升预测精度，并生成更高精度的SOCD空间分布栅格图像，帮助揭示SOCD空间分异规律，指导农业生产，验证了顾及土壤有机碳密度空间异质性的土地利用-克里格插值模型的有效性。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。