CN101718775A - 围垦地土壤中重金属含量的空间变异分布图生成方法 - Google Patents

Abstract

一种环境监测技术领域的围垦地土壤中重金属含量的空间变异分布图生成方法，包括：获取待测围垦地的采样点土壤样本；对采样点土壤样本中重金属元素含量数值进行数据预处理，分别获得各土壤重金属元素的拟插值子样本；对拟插值子样本进行探索性空间数据分析处理，获得采样点土壤样本的各向异性参数、步长和步数；对采样点土壤样本重金属元素含量数值进行函数拟合，获得适用于重金属元素含量插值的理论半方差函数模型；空间变异预测结果分布图的绘制。本发明所得分布图能够对围垦地土壤重金属含量的空间状况进行准确评价。

Description

围垦地土壤中重金属含量的空间变异分布图生成方法

技术领域

本发明涉及一种环境监测技术领域的方法，具体是一种围垦地土壤中重金属含量的空间变异分布图生成方法。

背景技术

上海地区围垦地土壤为盐渍土土壤，是由江河搬运泥沙与海水顶托，以及海流、潮汐、波浪等因素共同作用下的盐渍淤泥所形成，呈条带状分布在沿海及河口沿江地段，在分布和形成上属于滨海盐土类型。此类型土壤碱性大，质地粘重，有机质含量不高，而盐分含量很高。由于海滩开发的特殊性，其空间分布主要受排污口、围垦、潮滩位置和土壤粘土含量的影响，不仅空间变异大，时间上的变化也快。滩涂湿地土壤受重金属污染十分严重，带来了严重的生态风险。土壤重金属污染不仅影响农业的可持续发展而且会导致环境恶化、危害人体健康。近年来土壤重金属污染已日趋严重，尤其是在城市郊区的高投入的农田，因此预测围垦地土壤的重金属含量对当地的环境监管与治理具有极其重要的意义。

围垦地土壤重金属的分布受到诸多因素的影响，作为区域变量，其变异性中的确定性部分不是常量，而是随空间分布的，因而地统计方法结合GIS技术在土壤重金属污染调查预测中有广泛的应用。地统计方法(Geostistics)是根据已知采样点土壤样本的观测数据推测其周围未采样点土壤样本的相关特征，将不连续的点状数据形成连续的面状区域，来描述整个研究区域的土壤空间变异特征的方法。借助GIS可将土壤重金属污染的空间分布及预测结果进行可视化表达。现有的对土壤重金属含量的预测方法主要采用被称为最佳线性估值过程的泛克里格法(Universal Kriging，简称UK)。但在对具体插值模型的选择方面，缺乏对采样点土壤样本观察数据空间自相关性及向异性等空间特征的详细分析，在选择插值模型时对所涉及参数的确定过程讨论都较为模糊，这样极大地影响未采样点土壤样本预测结果的精度和可信度。另一方面，对大样本的插值过程缺乏整理和筛选，特别是对离群值的分析讨论及相关处理，使得预测结果的准确性大打折扣。

经过对现有技术的检索发现，刘庆等在《安全与环境学报》(2007，7(2)：109-113)上发表“基于GIS的农田土壤重金属空间分布研究”该研究在数据预处理上仅采用对数变换，当数据符合正态分布能满足变异函数的计算就进行插值计算，没有详细分析采样点土壤样本观测数据的空间特征，这样在插值过程中特异值的影响比较显著，预测结果的精度也会受到影响。由于围垦土重金属污染受海水侵蚀、土地利用类型及人类活动的影响，其空间变异特征明显，特别是较强的向异性，直接影响重金属含量插值模型尺度及相关参数的选择，从而也与土壤重金属预测紧密相关。

发明内容

本发明的目的在于克服现有技术的不足，提出一种围垦地土壤中重金属含量的空间变异分布图生成方法，根据围垦地土壤重金属的空间变异特征，对研究对象进行详细的探索性空间数据分析，既充分考虑各向异性及空间自相关性的影响，又通过选择合适的插值步长来确定插值的尺度参数，合理有效地处理或去除统计离群值和局部离群值，以确定最佳的插值尺度等参数，从而提高围垦地土壤中重金属含量泛克里格插值预测的精度和准确性。

本发明通过以下技术方案实现，本发明包括如下步骤：

第一步、获取待测围垦地的采样点土壤样本。

第二步、对采样点土壤样本中重金属元素含量数值进行数据预处理，分别获得各土壤重金属元素的拟插值子样本。

所述的数据预处理是指：

2.1)通过重金属元素含量数值的偏峰度计算及检验，考察其是否符合标准正态分布条件H₀，当符合标准正态分布条件H₀，则进行第三步操作；否则，执行步骤2.2。

所述的符合标准正态分布条件H₀是指：重金属元素含量数值的三阶矩(g₁)即偏度和四阶矩(g₂)也即峰度满足如下式的条件时，表明数据为标准正态分布：

$g_1=\frac{n\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^3}{(n-1)(n-2){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{n-1}\right)}^{\frac{3}{2}}}\≤u_{\mathrm{\α}}\sqrt{\frac{6n(n-1)}{(n+1)(n-2)(n+3)}};$

$g_2=\frac{n(n+1)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^4}{(n-1)(n-2)(n-3){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{(n-1)}\right)}^2}-\frac{3{(n-1)}^2}{(n-2)(n-3)}\≤u_{\mathrm{\α}}\sqrt{\frac{24n(n-1)}{(n-3)(n-2)(n+3)(n+5)}}$

式中u_α为在显著性水平为α时的u值(查U检验表可得)。

2.2)计算样本均值μ、标准差σ，去除区间(μ-1.96σ，μ+1.96σ)以外的值后再次执行步骤2.1，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.3。

2.3)对采样点土壤样本的数据进行对数变换，再次执行步骤2.2，其中变换后样本的μ，σ需重新计算，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.4。

2.4)对土壤样本数据，反复执行步骤2.3，当符合标准正态分布条件H₀则进行第三步，否则采用Box-Cox变换后再执行步骤2.2，当不符合标准正态分布条件H₀，则所使用样本数据不符合本专利所述的克里格插值的基本要求，停止计算重新采集待测围垦地的采样点土壤样本，否则执行第三步。

所述的Box-Cox变换是指：

$X^{\′}\left(\mathrm{\λ}\right)=\left\{\begin{array}{cc}\frac{X^{\mathrm{\λ}}-1}{\mathrm{\λ}}& \mathrm{\λ}\≠0\\ \ln X& \mathrm{\λ}=0\end{array}\right.;$

其中：X表示原样本，λ是一个待定变换参数。针对不同的λ包括：对数变换(λ＝0)、平方根变换(λ＝0.5)或倒数变换(λ＝-1)。

第三步、对拟插值子样本进行探索性空间数据分析处理，获得采样点土壤样本的各向异性参数、步长和步数。

所述的探索性空间数据分析处理是指：

3.1)计算拟插值子样本不同步长，即土壤样本点之间的距离下的半变异函数。

所述的半变异函数具体是指：

$\mathrm{\γ}\left(h\right)=\frac{\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2}{2N\left(h\right)},$

其中：区域化变量Z(x)在点x_i和x_i+h处的值Z(x_i)与Z(x_i+h)差的方差的一半为区域化变量Z(x_i)的半变异函数，即半方差函数。

3.2)分析半方差函数计算结果并将其中局域离群值及其所对应的采样点剔除；

所述的局域离群值是指：采用给定步长l，将所有点对距离h的集合分为n组，若半方差值高于μm+2.33σm(μm，σm分别为第m步时半方差值样本的均值和标准差，m＝1，2，…，n)，则记为离群半方差，并记录相应的点对；当n组离群半方差的相应点对中的一个点与5个以上离群半方差有关且这类点的个数小于总样本数的5％，那么这些点就被认为是离群值。否则就不认为是局域离群值。

第四步、对采样点土壤样本重金属元素含量数值进行函数拟合，获得适用于重金属元素含量插值的理论半方差函数模型，具体包括以下步骤：

4.1)选择半方差函数模型中一种模型(如球状模型)，基于最小块金效应原则，分别计算分析不同步长所对应的基台值与块金值，通过绘制基台值-步长曲线和块金效应-步长曲线，将基台值平稳区间内的块金效应达到最小值时的步长作为插值的参数。

4.2)针对半方差函数模型的圆形模型、球状模型、四球模型、五球模型、指数模型、高斯模型等六种理论模型，选取插值获得的均方根误差(RMSE)最接近于0和均方根标准预测误差(RMSSE)最接近于1的模型，由此获得最优的重金属元素含量空间插值模型。

所述的 $\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[z\left(x_i\right)-z*\left(x_i\right)]}^2},$ $\mathrm{RMSSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left\{\frac{[z\left(x_i\right)-z*\left(x_i\right)]}{\mathrm{\σ}\left(x_i\right)}\right\}}^2},$

其中z(x_i)为位置x_i处的实测值，z*(x_i)为位置x_i处的估计值。

第五步、空间变异预测结果分布图的绘制：根据重金属元素含量数值插值模型和采样点土壤样本的拟插值对采样点土壤样本重金属元素含量数值进行克里格插值，获得土壤重金属含量的空间变异分布图。同时随机选取10％该采样点土壤样本中的重金属元素进行交叉验证，即通过插值模型计算其预测值，获得随机校验均方根误差和随机校验均方根标准预测误差。

本发明通过对围垦地土壤重金属样本数据空间变异特征的深入分析，获得拟插值样本数据及插值模型，并通过插值生成土壤重金属因子的空间变异分布图。结合相应的土壤重金属的评价标准，可对围垦地土壤重金属含量的空间状况进行评价。

附图说明

图1为本发明步骤示意图。

图2实施例1中砷插值的半方差云图。

图3实施例1中砷的半变异云图在剔除离群值前后的变化；

其中：图3(a)为未剔除离群值的半变异云图，图3(b)为剔除离群值的半变异云图。

图4实施例1中不同尺度下砷的变异性；

其中：图4(a)为大尺度下半方差表面图，图4(b)为小尺度下半方差表面图。

图5实施例1中砷的空间变异参数与研究尺度的关系；

其中：图5(a)为步长与块金效应的关系图，图5(b)为步长与基台值的关系图。

图6实施例2中土壤砷的预测结果图。

图7实施例2中铜插值半方差云图。

图8实施例2中铜的半变异云图在剔除离群值前后的变化；

其中：图8(a)为实施例2中未剔除离群值的半变异云图，图8(b)为实施例2中剔除离群值的半变异云图。

图9实施例2中不同尺度下观察铜的变异性；

其中：图9(a)为实施例2大尺度下半方差表面图，图9(b)为实施例2小尺度下半方差表面图。

图10实施例2中铜的空间变异参数与研究尺度的关系；

其中：图10(a)为步长与块金效应的关系图，图10(b)为步长与基台值的关系图。

图11实施例2中土壤铜的预测结果图。

具体实施方式

下面结合附图对本发明的实施案例作详细说明：本实施案例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1：如图1所示，本实施包括以下步骤：

第一步、获取上海崇明岛地区围垦地采样点土壤样本的重金属砷的含量；

第二步、对采样点土壤样本中重金属元素含量数值进行数据预处理，分别获得各土壤重金属元素的拟插值子样本。

所述的数据预处理是指：

2.1)通过重金属元素含量数值的偏峰度计算及检验，考察其是否符合标准正态分布条件H₀，当符合标准正态分布条件H₀，则进行第三步操作；否则，执行步骤2.2。

所述的符合标准正态分布条件H₀是指：重金属元素含量数值的三阶矩(g₁)即偏度和四阶矩(g₂)也即峰度满足如下式的条件时，表明数据为标准正态分布：

$g_1=\frac{n\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^3}{(n-1)(n-2){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{n-1}\right)}^{\frac{3}{2}}}\≤u_{\mathrm{\α}}\sqrt{\frac{6n(n-1)}{(n+1)(n-2)(n+3)}};$

$g_2=\frac{n(n+1)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^4}{(n-1)(n-2)(n-3){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{(n-1)}\right)}^2}-\frac{3{(n-1)}^2}{(n-2)(n-3)}\≤u_{\mathrm{\α}}\sqrt{\frac{24n(n-1)}{(n-3)(n-2)(n+3)(n+5)}}$

式中u_α为在显著性水平为α时的u值(查U检验表可得)。

2.2)计算样本均值μ、标准差σ，去除区间(μ-1.96σ，μ+1.96σ)以外的值后再次执行步骤2.1，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.3。

2.3)对采样点土壤样本的数据进行对数变换，再次执行步骤2.2，其中变换后样本的μ，σ需重新计算，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.4。

2.4)对土壤样本数据，反复执行步骤2.3，当符合标准正态分布条件H₀则进行第三步，否则采用Box-Cox变换后再执行步骤2.2，当不符合标准正态分布条件H₀，则所使用样本数据不符合本专利所述的克里格插值的基本要求，停止计算重新采集待测围垦地的采样点土壤样本，否则执行第三步。

所述的Box-Cox变换是指：

$X^{\′}\left(\mathrm{\λ}\right)=\left\{\begin{array}{cc}\frac{X^{\mathrm{\λ}}-1}{\mathrm{\λ}}& \mathrm{\λ}\≠0\\ \ln X& \mathrm{\λ}=0\end{array}\right.;$

其中：X表示原样本，λ是一个待定变换参数。针对不同的λ包括：对数变换(λ＝0)、平方根变换(λ＝0.5)或倒数变换(λ＝-1)。

如图2、图3和图4所示，第三步、对拟插值子样本进行探索性空间数据分析处理，获得采样点土壤样本的各向异性参数、步长和步数。

所述的探索性空间数据分析处理是指：

3.1)计算拟插值子样本不同步长，即土壤样本点之间的距离下的半变异函数。

所述的半变异函数具体是指：

$\mathrm{\γ}\left(h\right)=\frac{\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2}{2N\left(h\right)},$

其中：区域化变量Z(x)在点x_i和x_i+h处的值Z(x_i)与Z(x_i+h)差的方差的一半为区域化变量Z(x_i)的半变异函数，即半方差函数。

3.2)分析半方差函数计算结果并将其中局域离群值及其所对应的采样点剔除；

所述的局域离群值是指：采用给定步长l，将所有点对距离h的集合分为n组，若半方差值高于μm+2.33σm(μm，σm分别为第m步时半方差值样本的均值和标准差，m＝1，2，…，n)，则记为离群半方差，并记录相应的点对；当n组离群半方差的相应点对中的一个点与5个以上离群半方差有关且这类点的个数小于总样本数的5％，那么这些点就被认为是离群值。否则就不认为是局域离群值。

如图4和图5所示，第四步、对采样点土壤样本重金属元素含量数值进行函数拟合，获得适用于重金属元素含量插值的理论半方差函数模型，具体包括以下步骤：

4.1)选择半方差函数模型中一种模型(如球状模型)，基于最小块金效应原则，分别计算分析不同步长所对应的基台值与块金值，通过绘制基台值-步长曲线和块金效应-步长曲线，将基台值平稳区间内的块金效应达到最小值时的步长作为插值的参数。

4.2)针对半方差函数模型的圆形模型、球状模型、四球模型、五球模型、指数模型、高斯模型等六种理论模型，选取插值获得的均方根误差(RMSE)最接近于0和均方根标准预测误差(RMSSE)最接近于1的模型，由此获得最优的重金属元素含量空间插值模型。

表  六种理论模型

式中：c₀表示块金方差(间距为0时的半方差)，c为结构方差，c₀+c为基台值(半方差函数随间距递增到一定程度后出现的平稳值)，a为变程(半方差达到基台值的样本间距)。如球状模型，a表示观测点之间的最大相关距离，而高斯模型的最大相关距离为(3)^1/2a，指数模型的最大相关距离为3a。

所述的 $\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[z\left(x_i\right)-z*\left(x_i\right)]}^2},$ $\mathrm{RMSSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left\{\frac{[z\left(x_i\right)-z*\left(x_i\right)]}{\mathrm{\σ}\left(x_i\right)}\right\}}^2},$

其中z(x_i)为位置x_i处的实测值，z*(x_i)为位置x_i处的估计值。

第五步、空间变异预测结果分布图的绘制：根据重金属元素含量数值插值模型和采样点土壤样本的拟插值对采样点土壤样本重金属元素含量数值进行克里格插值，获得土壤重金属含量的空间变异分布图。同时随机选取10％该采样点土壤样本中的重金属元素进行交叉验证，即通过插值模型计算其预测值，获得随机校验均方根误差和随机校验均方根标准预测误差。

如图6所示，为本实施例最终获得的重金属各指标的空间变异分布图。

表1  砷交叉验证插值结果误差分析表

注：表1中Source_ID为采样点的代码，As的实测值为该点上砷的实际测量值，As的预测值为经过插值以后的所得的预测值，As的误差为实测值与预测值之间的误差。

实施例2：如图1所示，本实施例包括以下步骤：

第一步、获取上海崇明岛地区围垦地采样点土壤样本的重金属铜的含量；

第二步、数据预处理：对采样点土壤样本重金属铜的数据进行偏峰度检验。检验结果为α＝0.01水平，表示拒绝标准正态分布假设，然后对采样点土壤样本的数据进行对数变换。变换后再进行检验，结果仍不服从正态分布，再进行二次对数变换，变换后得到拟插值的数据。

第三步、探索性空间数据分析：通过计算两点之间的半变异函数，生成半变异函数云图，如图7所示。对重金属元素含量数值半变异函数云图进行仔细观察，检验数据集中空间自相关的局部特征并查找离群值。结合第一步的统计分析特征，确定并剔除离群值，如图8所示。通过对变异函数云图的观察分析获得数据的向异性、周期性等空间分布状况及相关参数，如图10所示。

第四步、重金属元素含量数值各向异性建模：利用最小块金效应原则：在一定的步数条件下，分别绘制基台值-步长、块金效应-步长曲线，获得重金属元素含量数值插值模型。

如图11所示。观察在基台值平稳区间内的块金效应达到最小值时的步长，以此作为插值的参数。在20米的步数下，选择的步长为1750米。其他参数分别为偏基台值为0.01067；块金值为0.030159；块金效应为0.378542；变异角度为东北偏东10度；基台值为0.082118。

第五步、空间变异预测结果分布图的绘制：根据均方根误差和均方根标准预测误差，选取10％的采样点土壤样本对其进行交叉验证，重金属元素铜插值所选用的模型为高斯模型。得出克里格插值的平均误差为9，未经探索性空间数据分析的克里格插值的平均误差为9.994564。所以在对数据进行插值时，本实施例的准确性要优于未经过探索性空间分析的方法。然后，根据所选模型的插值结果，借助GIS获得围垦地土壤重金属各指标的空间变异分布图。

如图9和图11所示，为本实施例最终获得的重金属各指标的空间变异分布图。

表2  铜交叉验证插值结果误差分析表

注：表2中Source_ID为采样点的代码，Cu的实测值为该点上砷的实际测量值，Cu的预测值为经过插值以后的所得的预测值，Cu的误差为实测值与预测值之间的误差。

Claims (8)

1.一种围垦地土壤中重金属含量的空间变异分布图生成方法，其特征在于，包括如下步骤：

第一步、获取待测围垦地的采样点土壤样本；

第二步、对采样点土壤样本中重金属元素含量数值进行数据预处理，分别获得各土壤重金属元素的拟插值子样本；

第三步、对拟插值子样本进行探索性空间数据分析处理，获得采样点土壤样本的各向异性参数、步长和步数；

第四步、对采样点土壤样本重金属元素含量数值进行函数拟合，获得适用于重金属元素含量插值的理论半方差函数模型；

第五步、空间变异预测结果分布图的绘制。

2.根据权利要求1所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的数据预处理是指：

2.1)通过重金属元素含量数值的偏峰度计算及检验，考察其是否符合标准正态分布条件H₀，当符合标准正态分布条件H₀，则进行第三步操作；否则，执行步骤2.2；

2.2)计算样本均值μ、标准差σ，去除区间(μ-1.96σ，μ+1.96σ)以外的值后再次执行步骤2.1，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.3；

2.3)对采样点土壤样本的数据进行对数变换，再次执行步骤2.2，其中变换后样本的μ，σ需重新计算，当符合标准正态分布条件H₀则进行第三步，否则执行步骤2.4；

2.4)对土壤样本数据，反复执行步骤2.3，当符合标准正态分布条件H₀则进行第三步，否则采用Box-Cox变换后再执行步骤2.2，当不符合标准正态分布条件H₀，则所使用样本数据不符合本专利所述的克里格插值的基本要求，停止计算重新采集待测围垦地的采样点土壤样本，否则执行第三步。

3.根据权利要求2所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的符合标准正态分布条件H₀是指：重金属元素含量数值的三阶矩g₁，即偏度和四阶矩g₂满足如下式的条件时，表明数据为标准正态分布：

$g_1=\frac{n\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^3}{(n-1)(n-2){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{n-1}\right)}^{\frac{3}{2}}}\≤u_a\sqrt{\frac{6n(n-1)}{(n+1)(n-2)(n+3)}};$

$g_2=\frac{n(n+1)\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^4}{(n-1)(n-2)(n-3){\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{{(x_i-\stackrel{\‾}{x})}^2}{(n-1)}\right)}^2}-\frac{3{(n-1)}^2}{(n-2)(n-3)}\≤u_a\sqrt{\frac{24n(n-1)}{(n-3)(n-2)(n+3)(n+5)}}$

式中u_α为在显著性水平为α时的u值。

4.根据权利要求2所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的Box-Cox变换是指：

$X^{\′}\left(\mathrm{\λ}\right)=\left\{\begin{array}{cc}\frac{X^{\mathrm{\λ}}-1}{\mathrm{\λ}}& \mathrm{\λ}\≠0\\ \ln X& \mathrm{\λ}=0\end{array}\right.;$

其中：X表示原样本，λ是一个待定变换参数，针对不同的λ包括：

当λ＝0时为对数变换；

当λ＝0.5时为平方根变换；

当λ＝-1时为倒数变换。

5.根据权利要求1所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的探索性空间数据分析处理是指：

3.1)计算拟插值子样本不同步长，即土壤样本点之间的距离下的半变异函数：

$\mathrm{\γ}\left(h\right)=\frac{\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\Σ}}}{[Z\left(x_i\right)-Z(x_i+h)]}^2}{2N\left(h\right)},$

其中：区域化变量Z(x)在点x_i和x_i+h处的值Z(x_i)与Z(x_i+h)差的方差的一半为区域化变量Z(x_i)的半变异函数，即半方差函数；

3.2)分析半方差函数计算结果并将其中局域离群值及其所对应的采样点剔除。

6.根据权利要求1所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的局域离群值是指：采用给定步长l，将所有点对距离h的集合分为n组，当半方差值高于μm+2.33σm，其中μm和σm分别为第m步时半方差值样本的均值和标准差，m＝1，2，…，n，则记为离群半方差，并记录相应的点对；当n组离群半方差的相应点对中的一个点与5个以上离群半方差有关且这类点的个数小于总样本数的5％，那么这些点就被认为是离群值，否则就不认为是局域离群值。

7.根据权利要求1所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的第四步具体包括以下步骤：

4.1)选择半方差函数模型中的一种，基于最小块金效应原则分别计算分析不同步长所对应的基台值与块金值，通过绘制基台值-步长曲线和块金效应-步长曲线，将基台值平稳区间内的块金效应达到最小值时的步长作为插值的参数；

4.2)针对半方差函数模型，选取插值获得的均方根误差最接近于0和均方根标准预测误差最接近于1的模型，由此获得最优的重金属元素含量空间插值模型，所述均方根误差为：

$\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[z\left(x_i\right)-z*\left(x_i\right)]}^2},$

所述均方根标准预测误差为：

其中：z(x_i)为位置x_i处的实测值，z*(x_i)为位置x_i处的估计值。

8.根据权利要求1所述的围垦地土壤中重金属含量的空间变异分布图生成方法，其特征是，所述的第五步具体包括以下步骤：根据重金属元素含量数值插值模型和采样点土壤样本的拟插值对采样点土壤样本重金属元素含量数值进行克里格插值，获得土壤重金属含量的空间变异分布图，同时随机选取10％该采样点土壤样本中的重金属元素进行交叉验证，即通过插值模型计算其预测值，获得随机校验均方根误差和随机校验均方根标准预测误差。

Images