一种基于协同克里金插值法的土壤锰元素含量预测方法与流程

本发明属于农业信息化技术领域，涉及一种土壤金属元素含量预测方法，具体地说，涉及一种基于协同克里金插值法的土壤锰元素含量预测方法。

背景技术：

锰元素是植物生长的必要营养元素，能促进植物的光合作用，提高品质，并且能够促进植物的新陈代谢，增强磷和钙的吸收，提高植物的呼吸强度，但若锰含量过多会导致土壤重金属污染。因此，对锰元素含量的预测研究具有重要的现实意义。然而，目前研究只考虑单一因素作为辅助变量进行协同克里金插值，本发明在果园地、菜地、水浇地、灌溉水田和旱地这五种耕地类型下，利用主成分分析法将土壤锌、铁和铜三种金属元素组合成一个综合因素作为辅助变量对土壤锰元素含量进行协同克里金插值预测。

土壤重金属数据一般通过采样点实测得到，但人力、物力、财力限制使得获取的采样点在数量和空间分布上及其有限。因此，若要了解区域尺度上土壤重金属空间分布，就需要对点尺度上的土壤重金属测量值进行空间插值。土壤重金属含量有很多预测方法，但插值结果受制于采样点数据的稀疏程度和空间均匀性，具有不确定性、离散性高，很难得到精确的插值效果。现有研究中，往往不重视空间数据分析、插值方法的筛选和插值参数的优化，且对其空间相关性和其他因素的影响考虑较少，不适的插值方法和插值参数会造成结果的偏差，使数据具有的某种空间变化趋势难以辨别。以往的研究一般将单一因素作为辅助变量考虑，而实际情况很多变量之间相互影响，将这些变量综合考虑，再运用协同克里金插值更有意义。

技术实现要素：

本发明提出了一种基于协同克里金插值法的土壤锰元素含量预测方法，考虑多个因素作为辅助变量进行协同克里金插值。

其技术方案如下：

一种基于协同克里金插值法的土壤锰元素含量预测方法，包括以下步骤：

步骤1、方差分析

为了探讨耕地类型对土壤重金属含量的影响，使用方差分析法研究了果园地、菜地、水浇地、灌溉水田和旱地五种耕地类型对土壤重金属的影响。五种耕地类型为自变量，是分类型变量；重金属锌、铁、铜、锰、硼、硫的含量为因变量，是数值型变量。

步骤2、协同克里金插值法

使用协同克里金插值法cokriging预测了土壤重金属锰的含量；为了解决协同克里金插值时土壤锌、铁、锰这三个金属辅助变量的权重问题，使用主成分分析法pca将这三个存在相关性的变量通过降维技术转换为一个综合因子。然后，利用综合因子作为辅助变量对主变量锰进行协同克里金插值。

根据主成分分析原理，将土壤中重金属zn、fe、cu作为一个综合土壤重金属变量，主成分构成公式如下：

fi＝α1izn+a2ife+a3icu

式中，fi为第i主成分变量，aji(j＝1，2，3)分别表示为第i主成分对应原始变量zn、fe、cu的系数，度量了相应原始变量对fi的重要性。zn表示zn的含量，fe表示fe的含量，cu表示cu的含量。

综合因子的构成公式如下：

y＝∑λifi

式中，y为综合因子，fi为第i主成分，λi为主成分的贡献率。

协同克里金法插值的计算公式如下：

式中，z^*(l，k)0为(l，k)0处土壤重金属mn含量的估值，n和m为主变量和辅助变量的采样点数，z1(l，k)1i为各点土壤重金属mn的含量，λi为赋予各点土壤重金属mn含量的一组权重系数，且∑λi＝1；z2(l，k)2j为各点土壤综合重金属含量，λi为赋予各点综合土壤重金属含量的一组权重系数。

半方差函数公式如下：

式中，z1(l，k)为土壤中重金属mn的含量值，z2(l，k)为土壤中重金属综合金属含量值，且z1(l，k)、z2(l，k)在同一位置上，l+hl为与l距离为h的空间位置，k+hk为与k距离为h的空间位置。

步骤3、估值验证及精度评价

为了检验锰含量预测值的准确性，采用随机抽样的方式用80％的样本点作为训练集，20％的样本点数据作为测试集，用交叉验证方法对模型的精度进行检验，即锰含量的实测值与预测值之间的误差。采用(1)标准平均值误差ms表示评价结果的无偏性，其值越接近于0越好；(2)均方根误差rmse是估值方法精确性的一种度量，值越小越好；(3)标准化均方根误差rmsse其值应该越接近于1越好，小于1，说明对预测值高于真实值，大于1，则说明对预测值低于真实值；

式中，z(xi)为实测值，z(xoi)为预测值，n为检验采样点的数值。σ(xi)为在xi处的方差平方根。

步骤4、软件平台分析

r3.3.1：基本统计分析，正态性检验，果园地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.36；旱地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.23；灌溉水田地土壤重金属锰含量数据经log转换为最优正态分布；水浇地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.3；菜地土壤重金属锰含量数据经log转换为最优正态分布；方差分析、相关性分析和主成分分析。

arcgis10.2.2：地统计学分析，构建变异函数、cokriging插值。

进一步，采用2个半方差函数的模型来求解模型，即2个计算模型方程分布为：

(1)球形计算模型

式中，c为偏基台值，表示半方差函数在h大于变程时的值；a为变程，表示区域化变量在空间上相关性的范围；h表示为模型计算滞后系数。

(2)指数计算模型

式中，c为偏基台值，表示半方差函数在h大于变程时的值；a为变程，表示区域化变量在空间上相关性的范围；h表示为模型计算滞后系数。

根据无偏最优估计和∑λi＝1得：

式中，cij＝c(0)-γ(h)，(i，j＝1，2......n)，c(0)为基台值，表示半方差函数在h大于变程时的值，是块金值c0和提高的和，提高表示在取得有用数据标准上时，可观测得到的变异幅度大小；块金值c0表示在很短距离内有较大的空间变异性，不论h多小，两个随机变量之间不相关为0时的值。

从而得到一组权重系数λi值，估算出z^*(l，k)0。

本发明的有益效果为：

本发明通对过模型精度的结果验证分析，得出本方法对土壤锰元素含量预测的准确度更高。因而，基于协同克里金法土壤金属锰含量的预测具有重要的实际意义。

附图说明

图1基于协同克里金插值法的土壤锰元素含量预测方法的流程图；

图2研究区样本分布图；

图3基于cokriging法的各耕地类型土壤锰预测分布图，其中：

a)果园地；b)菜地；c)水浇地；d)灌溉水田；e)旱地。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案作进一步详细地说明。如图1所示。

1.数据与方法

1.1数据

1.1.1研究区概况

北京市房山区位于东经115°25′—116°10′与北纬39°30′—39°50′之间，属于暖温带半湿润季风大陆性气候区，境内地貌较为复杂，山区与平原之间相对高差悬殊，气候也有明显的差异。由于经常受到强烈的人为活动和自然生态过程的影响，区域内土壤有较强的变异空间。全区主要有五种地物类型，分别是菜地、果园地、旱地、水浇地和灌溉水田(图2)。

1.1.2数据来源

本研究收集了北京市房山区1497个样本点，其中果园地390个样本点，水浇地755个样本点，菜地53个样本点，灌溉水田21个样本点，旱地278个样本点，样本点数据包括样本点的经纬度，土壤锌、铁、铜、锰、硼和硫的含量。

1.2方法

1.2.1方差分析

为了探讨耕地类型对土壤重金属含量的影响，本发明使用方差分析法研究了北京市房山区果园地、菜地、水浇地、灌溉水田和旱地五种耕地类型对土壤重金属的影响。在本研究中，五种耕地类型为自变量，是分类型变量；重金属锌、铁、铜、锰、硼、硫的含量为因变量，是数值型变量。

1.2.2协同克里金插值法

本发明使用协同克里金插值法(cokriging)预测了土壤重金属锰的含量。为了解决协调克里金插值时土壤重金属锌、铁、锰这三个辅助变量的权重问题，使用主成分分析法(pca)将这三个存在相关性的变量通过降维技术转换为一个综合因子。然后利用综合因子作为辅助变量对主变量锰进行协同克里金插值。

根据主成分分析原理，本发明将土壤中重金属zn、fe、cu作为一个综合土壤重金属变量，主成分构成公式如下：

fi＝α1izn+a2ife+a3icu

式中，fi为第i主成分变量，aji(j＝1，2，3)分别表示为第i主成分对应原始变量zn、fe、cu的系数，度量了相应原始变量对fi的重要性。本发明中，zn表示zn的含量，fe表示fe的含量，cu表示cu的含量。

综合因子的构成公式如下：

y＝∑λifi

式中，y为综合因子，fi为第i主成分，λi为主成分的贡献率。

协同克里金法插值的计算公式如下：

式中，z^*(l，k)0为(l，k)0处土壤重金属mn含量的估值，n和m为主变量和辅助变量的采样点数，z1(l，k)1i为各点土壤重金属mn的含量，λi为赋予各点土壤重金属mn含量的一组权重系数，且∑λi＝1；z2(l，k)2j为各点土壤综合重金属含量，λi为赋予各点综合土壤重金属含量的一组权重系数。

半方差函数公式如下：

式中，z1(l，k)为土壤中重金属mn的含量值，z2(l，k)为土壤中重金属综合金属含量值，且z1(l，k)、z2(l，k)在同一位置上，l+hl为与l距离为h的空间位置，k+hk为与k距离为h的空间位置。

本发明主要采用2个半方差函数的模型来求解模型，即2个计算模型方程分布为：

(1)球形计算模型

式中，c为偏基台值，表示半方差函数在h大于变程时的值；a为变程，表示区域化变量在空间上相关性的范围；h表示为模型计算滞后系数。

(2)指数计算模型

式中，c为偏基台值，表示半方差函数在h大于变程时的值；a为变程，表示区域化变量在空间上相关性的范围；h表示为模型计算滞后系数。

根据无偏最优估计和∑λi＝1可得：

式中，cij＝c(0)-γ(h)，(i，j＝1，2......n)，c(0)为基台值，表示半方差函数在h大于变程时的值，是块金值c0和提高的和，提高表示在取得有用数据标准上时，可观测得到的变异幅度大小；块金值c0表示在很短距离内有较大的空间变异性，不论h多小，两个随机变量之间不相关为0时的值。

从而得到一组权重系数λi值，估算出z^*(l，k)0

1.2.3估值验证及精度评价

为了检验锰含量预测值的准确性，本发明采用随机抽样的方式用80％的样本点作为训练集，20％的样本点数据作为测试集，用交叉验证方法对模型的精度进行检验，即锰含量的实测值与预测值之间的误差。交叉检验的方法有很多，包括均差、相关系数、均方差、相对误差、绝对误差、估计优度等，本发明采用的是(1)标准平均值误差(ms)表示评价结果的无偏性，其值越接近于0越好；(2)均方根误差(rmse)是估值方法精确性的一种度量，值越小越好；(3)标准化均方根(rmsse)其值应该越接近于1越好，小于1，说明对预测值高于真实值，大于1，则说明对预测值低于真实值。

式中，z(xi)为实测值，z(xoi)为预测值，n为检验采样点的数值。σ(xi)为在xi处的方差平方根。

1.2.4软件平台

r3.3.1：基本统计分析，正态性检验(果园地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.36；旱地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.23；灌溉水田地土壤重金属锰含量数据经log转换为最优正态分布；水浇地土壤重金属锰含量的最优正态分布经box-cox数据变换的拟合参数λ＝0.3；菜地土壤重金属锰含量数据经log转换为最优正态分布)；方差分析、相关性分析和主成分分析。

arcgis10.2.2：地统计学分析(构建变异函数、cokriging插值)。

2.结果与分析

2.1房山区土壤重金属统计特征

房山区五种耕地类型的土壤重金属锰含量的统计特征值如表1所示，旱地锰元素含量的平均值最大为4.94mg/kg，水浇地锰元素含量的平均值最小为3.30mg/kg。从均值看，均属于低含量。灌溉水田的变异系数0.619，属于强变异系数，果园地与水浇地变异系数相差不大均属于中强度变异系数，菜地与旱地变异系数相差不大，属于稍强度变异系数。

表1.耕地类型土壤锰含量统计特征

由表2五种耕地类型土壤锰元素含量与其他元素含量的相关系数知在果园地中mn与fe、cu和综合因子具有很强的正相关性，且与综合因子(y)的相关性最高；在水浇地中mn与fe、cu和y有显著的正相关性；在旱地中，mn与cu和y成显著正相关性；在灌溉水田中，mn与zn、fe和y成显著负相关性，负相关性y>zn>fe；在菜地中，mn与fe成极显著正相关性，与y显著正相关性。由协同克里金插值法可知，辅助变量与主变量之间存在相关性才能进行协同克里金插值。因而本发明中的锰元素含量预测可以通过综合因子作为辅助变量进行插值。

表2.耕地类型土壤锰元素含量与其他元素含量的相关系数

注：*表示p＜0.05，**表示p＜0.01的极显著水平。

2.2地物类型对土壤元素含量的影响

用方差分析分析五种耕地类型对土壤重金属元素的影响(表3)，结果表明：地物类型对土壤重金属锰和铁含量有显著影响(p<0.01)；而对土壤重金属硼、硫没有显著影响(p>0.01)，地物类型对土壤锰元素影响最大，然后依次为铁>锌>铜>硼>硫。因而，本发明研究北京市房山区五种地物类型锰元素的含量。

表3.耕地类型对土壤重金属的显著性影响

2.3主成分分析的综合因子

用pca方法分别提取房山区五种耕地类型土壤重金属zn、fe、cu变量的主成分，本发明应用主成分主要目的为使用综合总金属变量，能够很好的解决在协同克里金法插值时多因子辅助变量权重的“黑匣子效应”。

得各耕地类型的主成分构成如下：

y果园地＝0.432zn+0.613fe+0.631cu

y水浇地＝0.654zn+0.173fe+0.647cu

y菜地＝0.689zn+0.535fe+0.471cu

y灌溉水田＝0.788zn+0.422fe+0.366cu

y旱地＝-0.635zn+0.341cu+0.712cu

2.4协同克里金预测结果与分析

2.4.1土壤重金属锰含量的空间变异特征

以重金属含量的综合因子为辅助变量和以单变量因子为辅助变量的交互变量最优半方差函数参数如表4所示，果园地、菜地、旱地土壤重金属mn的最优模型为指数模型，灌溉水田的最优模型为球形模型。以综合因子作为辅助变量时，除灌溉水田的块金系数为0.24，其余耕地类型的块金系数都在0.63-0.77之间，均为高强度的空间相关性，表明土壤重金属mn的空间变异性，受结构性因素和随机因素共同作用。且综合因子作为辅助变量比单一辅助变量的模型拟合精度要高点。

表4.土壤重金属mn半方差函数参数及其模型拟合精度

注：e表示指数模型；sph表示球形，()内数值表示单变量辅助插值结果。

2.4.2土壤重金属锰协同克里金插值分析

本发明采用主成分分析方法将铁、锌、铜三种与锰含量有相关关系的土壤重金属综合为一个变量，果园地锰含量与铜含量相关系数为0.24(p＜0.01)和综合因子含量相关系数0.283(p＜0.01)相关性都很强，与铁含量相关系数0.233(p＜0.05)相关性较强；水浇地锰含量与铜和铁含量相关系数分别为0.1(p＜0.05)、0.166(p＜0.01)；旱地锰含量与铜含量相关性很强0.5(p＜0.01)，与锌含量有较强的负相关性0.23，(p＜0.01)；灌溉水田的锰含量与锌、铁含量有极显著负相关性，相关系数分别为-0.546、-0.517(p＜0.01)；菜地锰含量也铁含量也极显著正相关0.553(p＜0.01)，综合后的变量都与锰含量相关性显著。通过主成分分析整合的综合重金属变量对土壤重金属锰进行协同克里金插值，能够很好的解决协同克里金插值时的多辅助变量权重的“黑匣子效应”。

协同克里金引入辅助变量可以提高插值的估值精度。本发明利用协同克里金法进行插值预测北京市房山区五种耕地类型的重金属锰含量。从插值结果图来看，预测结果较为真实地反映了实际的土壤重金属锰的含量。对于无值区域，通过协同克里金插值方法，获取其中的信息，预测其含量。从模型的精确验证结果来看，对未知区域的预测具有一定的可靠性。从而，预测分布图能够有效的展示土壤重金属锰的含量分布情况。果园地锰含量较高位于房山区中西部和房山区的东南方向，含量介于9.4mg/kg—11.8mg/kg之间，锰含量较少的位于房山区的西部，含量在0.116mg/kg—0.758mg/kg之间(图3a)；菜地锰含量较高位于房山区中东部，含量介于9.71—11.9mg/kg之间，南部菜地锰含量较低，含量介于0.633—1.25mg/kg之间(图3b)；房山区的西南部水浇地锰含量较高，含量介于9.42—11.9mg/kg之间，中部水浇地锰含量较低，含量介于0.138—0.855mg/kg之间(图3c)；灌溉水田的锰含量较高位于房山区西南部，含量介于7.23—8.11mg/kg之间(图3d)；房山区中部的旱地锰含量较高，含量介于9.66—12mg/kg之间，边缘地区的旱地锰含量较低(图3e)。因而可以根据不同植物生长对锰需求量的多少来种植植物，也可以采取一些措施来增加锰的含量和减少锰的含量使得缺少锰含量的土壤和锰含量过度的土壤达到中和。

房山区的果园地、水浇地、旱地、灌溉水田、菜地的土壤锰含量与综合后的辅助变量y的相关系数分别为0.283、0.111、0.186、-0.632、0.293，其插值精度2.129、2.623、2.275、2.639、2.158。除了灌溉水田的负相关性外，基本符合相关系数越大，插值精度越高，本发明可以表明相关性越强的辅助变量，越有助于提高协同克里金插值结果的准确性。

2.4.3协同克里金插值精度验证

克里金插值精度验证由表4显示，由多变量通过主成分分析法整合的综合因子作辅助变量的协同克里金插值精度比单变量作为辅助变量的协同克里金插值精度稍高。通过表4显示，果园地、水浇地、旱地三种耕地类型的协同克里金插值法效果都较好，灌溉水田和菜地稍微次点，与采样点较少有关系。果园、菜地、水浇地、旱地的ms分别为0.0046、-0.012、0.0003和0.0087，水浇地的平均预测结果稍低于真实观测值，果园、菜地和旱地的平均预测结果稍高于真实观测值。rmse计算结果园为2.289，值最小，表明预测精度最高，灌溉水田为3.079，值最大，表明预测精确度最低。而从rmsse来看，较接近1的是果园、水浇地、旱地，且果园为1.072大于1，说明高估了预测值的不确定性，而水浇地和旱地分别为0.978、0.979都小于1，说明低估了预测值的不确定性。菜地和灌溉水田的rmsse值与1的差值分别为0.126、0.165，都大于0.1，说明菜地和灌溉水田的预测精度稍差点。

3.结论

①果园地、菜地、水浇地、灌溉水田和旱地这五种耕地类型对土壤金属锰有显著性影响。②房山区的土壤金属锰含量表现出中等强度的空间相关性，土壤重金属锌、铁、铜均与土壤锰元素存在相关关系。③在果园地、菜地、水浇地、灌溉水田和旱地这五种耕地下，利用主成分分析将锌、铁和铜多变量组合为一个综合因子作为辅助变量与锰进行协同克里金插值。④通过模型精度结果验证分析，本发明的方法对土壤锰元素含量的预测的准确更高。因而，基于协同克里金法土壤金属锰含量的预测，是有一定意义的。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。