CN102999927B - 一种土壤污染物含量空间分布的精细分区方法 - Google Patents

Abstract

本发明的目标是提供一种土壤污染物含量空间分布的精细分区算法，解决目前空间插值方法由于出现插值跨界问题而无法实现土壤污染物分布的精细分区的难题。发明以土壤污染物含量分布与污染源之间的空间关系为参照，提取约束区域，建立基于约束区域的边界约束空间插值算法，将空间插值约束在可靠的空间区域，实现土壤污染物含量分布的精细分区。本发明的优点在于：提取污染源影响范围为边界约束条件，将对未知位点的空间插值约束在可靠的空间区域，实现对土壤污染物含量分布的分区的准确性和可靠性。本发明技术适用于土壤、生态、水文等领域的污染物含量分布的数字地图制作。

Description

一种土壤污染物含量空间分布的精细分区方法

技术领域

本发明属于土壤污染物含量的空间插值技术，更具体涉及一种土壤污染物含量空间分布的精细分区方法，适用于适用生态、水文等领域相关污染指标含量分布的数字地图制作。

背景技术

受人力和物力的影响，对研究区域土壤污染物含量的取样观测总是有限的，如何通过有限的、离散的观测样点数据准确地描述土壤污染物的面源分布，一直是人们感兴趣的难点问题。目前，描述土壤污染物空间分布信息的主要方法是空间插值。但是，现有空间插值方法都存在空间插值跨界问题，即在污染物含量分布差异的边界地带的插值结果不可靠，故不能提供土壤污染物含量的精细分区结果。用以下示例来说明该问题。图1中A、B分别代表污染物含量有差异（假定污染物含量的大小排序为：A<B）的两地块，黑色实心圆点表示监测样点，空心带数字1和2两圆点表示未监测点。用现有空间插值方法估算未知点1的污染物含量，由于插值过程中有地块B的样点参与，会导致对未知点1的污染物含量的插值结果偏高；而对未知点2的污染物含量的插值结果，由于地块A的样点参与而偏低。

发明内容

本发明的目的是在于针对现有技术存在的上述问题，提供及一种土壤污染物含量空间分布的精细分区方法，以污染物含量与污染源分布的空间关系为参照，边界约束分区算法提取空间插值的约束边界，将对未知位点的插值约束在可靠的空间区域内，实现对土壤污染物含量分布的精细分区。

为了实现上述目的，本发明采用以下技术方案：

一种土壤污染物含量空间分布的精细分区方法，包括以下步骤：

步骤1、通过对遥感影像或航片进行分类，提取被研究区域的约束区域BArea；

步骤2、记N个监测样点组成的样点集S={S(x₁,y₁),S(x₂,y₂),..,S(x_N,y_N)}，对任一监测样点S_i，1≤i≤N，其对应的空间坐标为S(x_i,y_i)、观测值为S(X_j)，监测样点集S包含于约束区域，即定义最小外接矩形MinRec为能覆盖所有监测样点且面积最小的规则矩形；

步骤3、将最小外接矩形MinRec分割成等大小的栅格单元，即待插值单元，最小外接矩形MinRec将对待插值单元的空间插值约束在监测样点能所覆盖最大区域内进行，即对仅位于最小外接矩形MinRec内的待插值单元进行插值，

步骤4、将约束区域BArea细分为非目标监测区域NMArea、非污染影响区域NPArea及不同污染影像区域PArea，PArea＝PArea₁∪PArea₂∪...,∪PArea_K，其中K为不同污染影像区域PArea的数目，确定优化的约束区域OBArea，即OBArea＝(NMArea∪NPArea∪PArea)∩MinRec，并分别对监测样点集S和待插值单元进行对应的非目标监测区域、非污染影响区域及不同污染影像区域分类；

步骤5、选择监测样点对待插值单元进行插值，监测样点和待插值单元具有步骤4中相同类别且空间邻近，待插值单元P(x_i,y_i)∈P，待插值单元对应的估算值P(X_i)基于以下公式

$P\left(X_i\right)=\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{j}S\left(X_j\right)$

其中：S(X_j)表示监测样点S(x_j,y_j)的观测值，λ_j表示权重系数，L为参与插值的样点数，1≤j≤L≤N，监测样点S(x_j,y_j)是待插值单元P(x_i,y_i)的L个最邻近点之一。

遥感影像或航片分辨率不超过5米。

步骤2中最小外接矩形MinRec的长度Length为沿水平方向监测样点之间的最大距离，最小外接矩形MinRec的高度Height为沿垂直方向监测样点之间的最大距离。

待插值单元的边长Ulength=Min（Length，Height）/250，即最小外接矩形MinRec的长度Length和高度Height的较小值除以250。

本发明与现有技术相比，具有以下优点和效果：

1、本发明充分考虑了污染物含量分布变异与污染源分布的关联性，将土壤污染实地调研资料（先验知识）与土壤污染物含量分布的空间分区过程进行了有机融合，由先验知识驱动空间分区过程。

2、本发明以污染物含量分布与污染源的空间关系为参照，提取约束区域图，将空间插值操作约束在有效、可靠的空间区域内，保障了对土壤污染物含量分布的空间分区结果的客观性和准确性。

3、本发明适用范围更广，不仅对简单地理环境下污染物含量分布的空间分区有效，而且可对复杂地理环境下土壤污染物含量分布实现精细的空间分区。

附图说明

图1为现有技术插值跨界问题的分析示意图；

图2为本发明方法提取约束区域的示意图（A显示了研究区域及监测样点分布，B显示了经分类细化后的约束区域）；

图3为本发明方法的最小外接矩阵栅格化及标识待插值单元（A显示了栅格化的最小外接矩形与约束区域进行空间交互操作后的结果，B显示了栅格化的最小外接矩形与细化后的约束区域进行空间交互操作后的结果）；

图4为两种经典空间插值方法的插值结果示意图（A为反距离加权插值方法，B为克里格插值方法）；

图5为本发明方法的插值结果示意图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

实施例1：

一种土壤污染物含量空间分布的精细分区方法，包括以下步骤：

步骤1、通过对遥感影像或航片进行分类，提取被研究区域的约束区域BArea；

约束区域BArea为空间插值操作所能覆盖的空间范围。约束区域的提取包括提取非监测目标区域NMArea和监测目标区域MArea两个部分。非目标监测区域指土壤污染物监测过程中通常不涉及的监测区域，如水域、建筑区、林地等。监测目标区域指土壤污染取样监测中实际的监测区域，通常是指农业用地。约束区域、非目标监测区域及监测目标区域的三者服从如下空间关系：

$\mathrm{BArea}=\mathrm{MArea}\&cap;\stackrel{\&OverBar;}{\mathrm{NMArea}}$

这里，为NMArea的补集。

由于空间分布上非目标监测区域与监测目标区域是交叉并存的，非目标监测区域的存在会影响土壤污染分区结果，故需要剔除非目标监测区域保障土壤污染分区结果的准确性。可利用遥感影像分类方法初步提取研究区域的约束区域。具体方法是：利用研究区域内空间分辨率不超过5米的遥感影像或者彩色航片，通过对影像/航片进行分类，分类结果图，即为所需的约束区域图。

步骤2、记N个监测样点组成的样点集S={S(x₁,y₁),S(x₂,y₂),..,S(x_N,y_N)}，对任一监测样点S_i，1≤i≤N，其对应的空间坐标为S(x_i,y_i)、观测值为S(X_j)，监测样点集S包含于约束区域，即定义最小外接矩形MinRec为能覆盖所有监测样点且面积最小的规则矩形；

在空间关系上，监测样点集S包含于约束区域，即定义最小外接矩形（MinRec）为能覆盖所有监测样点、且面积最小的规则矩形。这表明：最小外接矩形的长度Length为沿水平方向监测样点之间的最大距离，Length=Max(|x_i-x_j|)（i≠j,1≤i,j≤N）；最小外接矩形的高度（Height）为沿垂直方向监测样点之间的最大距离，Height=Max(|y_i-y_j|)（i≠j,1≤i,j≤N）。

步骤3、将最小外接矩形MinRec分割成等大小的栅格单元，即待插值单元，最小外接矩形MinRec将对待插值单元的空间插值约束在监测样点能所覆盖最大区域内进行，即对仅位于最小外接矩形MinRec内的待插值单元进行插值，

栅格化是将最小外接矩形MinRec分割成多个等大小的栅格单元待插值单元，Unit。待插值单元是空间插值的最小单元，在空间插值过程中，需要对所有待插值单元进行插值。

最小外接矩形将对待插值单元的空间插值约束在监测样点能所覆盖最大区域内进行，即对仅位于最小外接矩形内的待插值单元进行插值。待插值单元中土壤污染物含量被认为是均一；且在空间插值过程，用待插值单元的中心点坐标P(x_i,y_i)来表征待插值单元。用待插值单元中心点集P表征待插值单元集，则P={P(x₁,y₁),P(x₂,y₂),..,P(x_T,y_T)}（1≤i≤T）。与监测样点集S一样，中心点集P包含于约束区域，即

步骤4、将约束区域BArea细分为非目标监测区域NMArea、非污染影响区域NPArea及不同污染影像区域PArea，PArea＝PArea₁∪PArea₂∪...,∪PArea_K，其中K为不同污染影像区域PArea的数目，确定优化的约束区域OBArea，即OBArea＝(NMArea∪NPArea∪PArea)∩MinRec，并分别对监测样点集S和待插值单元进行对应的非目标监测区域、非污染影响区域及不同污染影像区域分类；

由于影响土壤污染物含量分布的污染源可能不一致，土壤污染物含量的空间分布在地理空间上是不连续的。对未知位点的土壤污染物含量的插值时，必须考虑参与插值的样点与未知点是否空间相似，即未知点与样点是否处于相同污染源影响范围内。因此，确认区域内污染源及其影响范围是关键。区域内的调查资料及监测样点数据可用于判断土壤污染物的污染源的类型。在确认了污染源后，可根据污染源周边的监测样点的污染物含量分布来确定污染源的潜在影响范围。

污染源类型及其潜在影响范围将作为先验知识，用于细化约束区域。细化约束区域就是将约束区域矢量图进一步细化成非目标监测区域NMArea、非污染影响区域NPArea  及若干不同污染影响区域PArea（PArea＝PArea₁∪PArea₂∪...,∪PArea_k）。考虑约束区域与最小外接矩形并不是完全空间匹配的，因此求解细化后的约束区域与最小外接矩形的空间交集，得到优化的约束区域OBArea    ，即O标识样点是基于OBArea对监测样点及待插值单元进行分类。可依据监测样点集S和待插值单元中心点集P与约束区域OBArea之间的空间关系，分别对监测样点和中心点进行分类标识，即细化成非目标监测区域、非污染影响区域及污染影响区域（包括若干污染影响区域子类）三个大类。

步骤5、选择监测样点对待插值单元进行插值，监测样点和待插值单元具有步骤4中相同类别且空间邻近，待插值单元P(x_i,y_i)∈P，待插值单元对应的估算值P(X_i)基于以下公式

$P\left(X_i\right)=\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{j}S\left(X_j\right)$

其中：S(X_j)表示监测样点S(x_j,y_j)的观测值，λ_j表示权重系数，L为参与插值的样点数，1≤j≤L≤N，监测样点S(x_j,y_j)是待插值单元P(x_i,y_i)的L个最邻近点之一。

步骤1中的遥感影像或航片分辨率不超过5米。

步骤2中最小外接矩形MinRec的长度Length为沿水平方向监测样点之间的最大距离，最小外接矩形MinRec的高度Height为沿垂直方向监测样点之间的最大距离。

待插值单元的边长Ulength=Min（Length，Height）/250，即最小外接矩形MinRec的长度Length和高度Height的较小值除以250。

步骤5中的监测样点与待插值单元均属于非目标监测区域时，权重系数λ_j=0；监测样点与待插值单元均不属于非目标监测区域时，权重系数λ_j＝M^-1×Dist'，其中，半方差矩阵 $M=\left(\begin{array}{ccccc}S\left(X_{11}\right)& S\left(X_{1j}\right)& ...& S\left(X_{1L}\right)& 1\\ S\left(X_{j1}\right)& S\left(X_{\mathrm{jj}}\right)& ...& S\left(X_{\mathrm{jL}}\right)& 1\\ S\left(X_{L1}\right)& S\left(X_{\mathrm{Lj}}\right)& ...& S\left(X_{\mathrm{LL}}\right)& 1\\ 1& 1& ...& 1& 0\end{array}\right),$

S(X_mz)用于说明半方差矩阵内参数定义，S(X_mz)为两邻近待插值单元P(x_i,y_i)的监测样点S(x_m,y_m)与S(x_z,y_z)的半方差，即为S(X_mz)=(S(X_m)-S(X_z))²/2，1≤m≤L，1≤z≤L，

向量 $\mathrm{Dist}=(\mathrm{Dist}\left(X_i^1\right),\mathrm{Dist}\left(X_i^2\right),...,\mathrm{Dist}\left(X_i^L\right),1);$

用于说明向量Dist内参数定义，为监测样点S(x_p,y_p)与待插值单元P(x_i,y_i)的欧式距离，即为1≤p≤L。

实施例2：

以下实例用于说明本发明，但不用来限制本发明的应用范围。本实例来源上海市一个有机农产品产地，产地东西长、南北宽各1000米和600米，总面积为0.6k平方米。该产地主要有两种土地利用类型：大田和设施蔬菜地。2010年按规则网格取样方法取样193份，监测该产地的土壤重金属Cu含量状况，样点分布图见（图2-A）。经调查发现该产地设施蔬菜土壤重金属Cu含量明显高于大田土壤的，这源于设施蔬菜的有机肥投入量大导致重金属Cu累积。因此，设施蔬菜地被认作该产地的重金属Cu的潜在污染区域。采用本发明制作该产地土壤重金属Cu含量的空间分区图，具体方法如下：

（1）提取约束区域

对获取的该产地的0.6空间分辨率彩色航片进行分类，取得该产地的大田和设施蔬菜两种土地利用类型分布图，该图用作该产地空间插值的约束区域（见图2-B）。

（2）创建最小外接矩形及其栅格化

分别沿水平和垂直方向计算所有监测样点的最大距离，得沿水平方向为550米，垂直方向为900米。按长900米、高550米，创建监测样点的最小外接矩形，以2.5米间距对该外接矩形进行栅格化，并与前期获取的约束区域图进行空间相交，结果见图3-A。

（3）细化约束区域及标识样点

由于该产地土地利用方式简单，故前期得到的约束区域不需要细化，可以直接用来标识样点，对监测样点标识结果见图2-B，对插值单元的标识结果见图3-B。

（4）空间插值并输出结果

为了体现本发明的技术优势，这里将本发明与两种现有的两种典型的空间插值方法（反距离加权插值和克里格插值）进行比较。在插值过程中，参与插值计算的最近邻监测样点数，即参数k统一为5。反距离加权插值方法的实验结果见图4-A，克里格插值方法的实验结果见图4-B，本发明的插值结果见图5。分析图4可以发现，反距离加权插值和克里格插值的结果图中，大田和设施蔬菜地的相邻近的边界区域出现了明显的插值跨界现象。邻近设施蔬菜地的大田土壤重金属Cu含量呈现偏高趋势，而邻近大田的设施蔬菜地土壤重金属Cu呈现偏低趋势；并且相对克里格方法的插值结果，这种插值跨界现象在反距离加权方法的插值图中更显著。然而用本发明提供的方法进行空间插值，插值结果图中不再出现插值跨界现象，这表明本发明可有效地解决现有空间插值方法的插值跨界问题，实现对土壤污染物含量分布的精细分区。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种土壤污染物含量空间分布的精细分区方法，其特征在于： 

步骤1、通过对遥感影像或航片进行分类，提取被研究区域的约束区域BArea，约束区域BArea为空间插值操作所能覆盖的空间范围； 

步骤2、记N个监测样点组成的样点集S＝{S(x₁,y₁),S(x₂,y₂),...,S(x_N,y_N)}，对任一监测样点S_i，1≤i≤N，其对应的空间坐标为S(x_i,y_i)、观测值为S(X_i)，监测样点集S包含于约束区域，即定义最小外接矩形MinRec为能覆盖所有监测样点且面积最小的规则矩形； 

步骤3、将最小外接矩形MinRec分割成等大小的栅格单元，即待插值单元，最小外接矩形MinRec将对待插值单元的空间插值约束在监测样点所能覆盖最大区域内进行，即对仅位于最小外接矩形MinRec内的待插值单元进行插值， 

步骤4、将约束区域BArea细分为非目标监测区域NMArea、非污染影响区域NPArea及不同污染影响区域PArea，PArea＝PArea₁∪PArea₂∪...,∪PArea_K，其中K为不同污染影响区域PArea的数目，确定优化的约束区域OBArea，即OBArea＝(NMArea∪NPArea∪PArea)∩MinRec，并分别对监测样点集S和待插值单元进行对应的非目标监测区域、非污染影响区域及不同污染影响区域分类； 

步骤5、选择监测样点对待插值单元进行插值，监测样点和待插值单元具有步骤4中相同类别且空间邻近，待插值单元P(x_i,y_i)∈P，P为待插值单元中心点集，待插值单元对应的估算值P(X_i)基于以下公式 

其中：S(X_j)表示监测样点S(x_j,y_j)的观测值，λ_j表示权重系数，L为参与插值的样点数，1≤j≤L≤N，监测样点S(x_j,y_j)是待插值单元P(x_i,y_i)的L个最邻近点之一。 

2.根据权利要求1所述的一种土壤污染物含量空间分布的精细分区方法，其特征在于：所述的步骤1中的遥感影像或航片分辨率不超过5米。 

3.根据权利要求1所述的一种土壤污染物含量空间分布的精细分区方法，其特征在于：所述的步骤2中最小外接矩形MinRec的长度Length为沿水平方向监测样点之间的最大距离，最小外接矩形MinRec的高度Height为沿垂直方向监测样点之间的最大距离。 

4.根据权利要求3所述的一种土壤污染物含量空间分布的精细分区方法，其特征在于：所述的待插值单元的边长Ulength＝Min(Length，Height)/250，即最小外接矩形MinRec的长度Length和高度Height的较小值除以250。