CN107729293B - 一种基于多元自适应回归的地理空间异常探测方法 - Google Patents

Abstract

本发明公开了一种基于多元自适应回归的地理空间异常探测方法，充分考虑了地理空间数据的相关性与异质性，针对获取的探测变量及其它地理变量进行数据预处理及回归变量的确定，以构建自适应空间邻近域，并基于自适应带宽的权函数选择和空间邻域实体专题属性变化梯度修复，构建地理加权回归模型，然后计算所有空间实体的稳健空间异常度，并组成集合，最后将异常度偏离均值一定倍数的采样点判别为空间异常。本发明不仅顾及空间相关性的自适应带宽选择和空间异质性的探测变量与其它地理变量间关系的定量表达，还增强了本发明的实用性和解释性。

Description

一种基于多元自适应回归的地理空间异常探测方法

技术领域

本发明属于空间数据挖掘与空间分析领域，尤其涉及一种基于多元自适应回归的地理空间异常探测方法

背景技术

异常探测最初源自统计学中的粗差探测研究，但在实践中发现有些“粗差”并不一定是错误，而往往隐含了某种特殊的规律或性质，具有重要的应用价值。空间异常探测是传统异常探测方法在地球信息科学领域的扩展和延伸，作为空间数据挖掘的一种重要手段，旨在从海量空间数据发现偏离整体或局部普遍分布模式的小部分异常实体，这部分异常实体通常蕴含着地理现象或地理过程的特殊发展规律。Hawkins最早给出异常的本质性定义，即“严重偏离其它对象的观测数据，以至于令人怀疑它是由不同机制产生的”。Shekhar等进一步拓展“异常”的定义到空间域，给出空间异常“专题属性与空间邻近域内的其它参考实体的专题属性显著不同的空间实体”，即违背地理学第一定律越近越相似的实体。

地理现象通常是多种地理要素综合作用的结果，且构成地理现象的各种地理要素都不是独立出现的，它们之间相互作用，可能是相互抑制，也可能是相互促进的。当前顾及多元数据的空间异常探测方法大多认为探测变量(即因变量，如气温)与其它地理变量(即自变量，如经纬度、海拔)之间具有等同的影响，采用马氏距离等进行等权处理，或者采用全局回归分析的方法进行处理。实际上不同地理变量对探测变量的影响存在很大差异，且全局方法无法顾及空间异质性对回归分析的影响，因此本发明采用地理加权回归的思想，定量描述探测变量与其它地理变量间的关系。

地理加权回归模型(Geographically Weighted Regression，GWR)是一种对空间子区域上自变量和因变量之间关系随着空间变化进行建模的非参数局部空间回归分析方法，其中自变量的回归系数是随着空间位置而变化的，是在局部加权回归模型的基础上引入局部光滑的思想对空间异质性进行建模，对于空间数据具有强大的局部分析能力。地理加权回归模型的核心在于空间权重的定义，包括权重函数选择与带宽选择。现有研究表明，权重函数的选择对地理加权回归模型的影响较小，而对带宽的选择非常敏感。且在参与回归的数据点包含潜在异常时会对回归结果造成较大影响。同时，现有多元空间异常探测大多没有顾及不同影响变量贡献权重的差异，且没有考虑空间异质性对回归因子的影响。

在实现本发明过程中，发明人发现现有技术中至少存在如下问题：地理空间异常探测中存在着难以顾及多元数据的稳健异常度量、在地理加权回归模型中难以处理潜在异常影响以及自适应选择带宽。

发明内容

为了解决上述技术问题，本发明提供一种基于多元自适应回归的地理空间异常探测方法，旨在解决多种地理变量相互作用下的地理空间异常模式的提取。本发明中，探测变量即指因变量，其它地理变量即指自变量，其技术方案如下：

步骤一数据预处理及回归变量的确定

1)针对获取的探测变量及其它地理变量，检验是否存在缺失，若存在，则进行插补；

2)对探测变量进行均值、中值、范围、标准差、变异系数、峰度和偏度等描述性统计；

3)针对其它地理变量，以模型拟合优度的检验统计量F1为检验统计量，采用逐步回归法(包括向前选择法、向后剔除法和逐步选择法)确定用来构建模型的回归变量。

步骤二自适应空间邻近域构建

1)根据待检测实体的空间位置(X、Y，若是经纬度则转换为投影坐标)构建Delaunay三角网，表达实体的粗空间邻接关系；

2)把Delaunay三角网所有边长按升序排列，构成边长序列E，序列中位于上、下四分位数之间所有边长的平均值作为稳健平均边长，利用下列公式进行计算：

其中，表示边e_i的边长，n表示Delaunay三角网的边长数目。

3)采用如下公式获取并打断Delaunay三角网中所有不合理边Unreasonable_E，表达采样点的细空间邻接关系：

其中，α为调节因子，用于调整不合理边的判断阈值，取值范围为[2，4]，默认设为3。

4)在打断不合理边的Delaunay三角网中，针对每个空间实体，采用如下公式计算其空间邻域半径R(Pi)：

其中，Mean²(P_i)表示点P_i的2阶邻域内所有边长的平均值；STD(P_j)表示点P_j1阶邻域内边长的标准差。

5)落在每个空间实体对应的空间邻域半径范围内的其它实体即为该实体的空间邻域。

步骤三地理加权回归模型构建

1)针对探测变量及其它独立地理变量，构建如下地理加权回归模型：

其中，(u_i，v_i)表示第i个采样点的地理坐标，y_i与x_i1，x_i2，...，x_ip表示自变量y和因变量x₁，x₂，...，x_p在位置(u_i，v_i)处的观测值；β_j(u_i，v_i)(j＝1，2，...，p)是第i个采样点上的第j个回归参数；ε_i～N(0，σ²)，Cov(ε_i，ε_j)＝0，(i≠j)。

2)基于自适应带宽的权函数选择：针对每个空间实体，根据距离回归点的距离远近，把实体空间邻域中的回归实体重新排序，空间邻域中距离回归点最远的空间实体的距离d_ij设为带宽b。设最近数据点(该点可以是回归点自身)的权重为1(该点可以是回归点自身)，其他数据点的权重按如下公式衰减：

3)空间邻域实体专题属性变化梯度修复

a)针对每个实体P，根据其空间邻域，采用下式计算专题变化梯度：

其中，f(P)表示采样点P的专题属性值，即探测变量值，D(P，P_i)表示实体P与邻域实体P_i间的欧氏空间距离。

b)令f(P)＝0，分别计算实体P与其空间邻域实体P_i间的专题属性变化梯度G(P，P_i)，并按升序排列获取序列G(P)，计算序列中位数M(P)；

c)针对任一邻域实体，采用下式计算专题属性变化梯度偏离GD(P_i)，并升序排列获取序列GD(P)：

GD(P_i)＝|G(P，P_i)-M(P)|

d)将邻域实体按专题属性变化梯度偏离划分为大、中、小三个等级，处于最大等级的[(n+1)/3]个实体组成待修复集合R(P)，采用如下公式进行修复：

f_R(P_i)＝M(P)*D(P，P_i)，

说明：修复是暂时的，是在每个实体地理加权回归的过程中进行，并不改变采样点的固有专题属性值。

4)采样下式进行回归参数估计：

其中，W_i为位置(u_i，v_i)处的权重矩阵，W_i＝diag(w_il，w_i2，...，w_in)，y_R为经过专题属性变化梯度修复的因变量矩阵。

步骤四稳健异常度量与异常判别

1)针对每个空间实体P，探测变量观测数据与地理加权回归预测数据间的差异即为实体P的稳健空间异常度RSOM(Robust Spatial Outlier Measure)，采用下列公式计算：

其中，f(P)表示空间实体P的探测变量观测数据，表示空间实体P的探测变量的地理加权回归预测值。

2)针对所有空间实体的稳健空间异常度，组成集合S_RsoM＝{RSOM(P₁)，RSOM(P₂)，...，RSOM(P_n)}，其中异常度偏离均值一定倍数的采样点判别为空间异常，采用下列公式计算：

S_outlier＝{P_i|RSOM(P_i)-μ＞kσ，P_i∈SDB}

其中，μ表示稳健异常度均值，σ表示稳健异常度标准差，k表示判别系数，默认设置为2。

本发明的有益效果在于：

本发明提出的一种基于多元自适应回归的地理空间异常探测方法，充分考虑了空间数据的相关性与异质性，主要技术优点如下：

1)顾及空间相关性的自适应带宽选择，增强了本发明的实用性；

2)顾及空间异质性的探测变量与其它地理变量间关系的定量表达，增强了本发明的解释性。

附图说明

图1是本发明实施例的方法流程图；

图2是本发明实施例土壤重金属采样点的空间分布及其高程空间分布图；

图3是本发明实施例土壤重金属Cr元素反距离插值空间分布图；

图4是本发明实施例土壤重金属采样点：(a)原始Delaunay三角网；(b)打断不一致边后的Delaunay三角网图；

图5是本发明实施例地理加权回归参数的空间分布：(a)常数参数；(b)高程参数；(c)距河流距离参数；(d)距工厂距离参数；

图6是本发明实施例空间异常采样点的空间分布图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

本发明适用于气象、环保等诸多领域，如极端气候事件检测、土壤重金属异常。本发明以土壤重金属采样数据为实施例进行说明，土壤重金属浓度异常分布模式是指某采样点的重金属浓度与其空间邻近采样点的重金属浓度具有明显差异。

土壤重金属浓度异常分布一方面可能预示异常采样点更接近污染源，有助于指导污染源头的排查，另一方面可能预示土壤重金属污染特殊的传播方式，有助于重金属污染传播的治理。

实施例采用我国华南某市的环保数据，共包含103个采样点，其空间分布如图2所示。土壤重金属含量主要受到自然和人为因素的影响，其中重金属铬Cr元素对人体危害较大，因变量选择Cr元素，其土壤重金属含量反距离插值空间分布如图3所示。解释变量选择采样点高程Ele、距相关工厂的最近距离Dist_Factory、距河流的最近距离Dist_River，单位均为米m。

如图1所示，基于多元自适应回归的地理空间异常探测方法，包括步骤：

步骤1、土壤重金属采样数据预处理及回归变量确定

1)针对获取的土壤重金属采样点数据集，检验是否存在缺失，若存在，则采用一定范围内采样点数据的反距离加权值进行插补；

2)对土壤重金属采样数据Cr进行均值、中值、范围、标准差、变异系数、峰度和偏度等描述性统计；

3)以模型拟合优度的检验统计量F₁为检验统计量，采用逐步回归法确定用来构建模型的回归变量，经分析，选择的参与回归变量为采样点高程Ele、距相关工厂的最近距离Dist_Factory、距河流的最近距离Dist_River。

步骤2、土壤重金属采样点自适应空间邻近域构建

1)根据土壤重金属采样点的空间位置(地理坐标需转换为投影坐标)构建Delaunay三角网，表达采样点的粗空间邻接关系，如图4(a)所示；

2)把Delaunay三角网所有边长按升序排列，构成边长序列E，序列中位于上、下四分位数之间所有边长的平均值作为稳健平均边长，利用下列公式进行计算：

其中，表示边e_i的边长，n表示Delaunay三角网的边长数目。

3)采用如下公式获取并打断Delaunay三角网中所有不合理边Unreasonable_E，表达采样点的细空间邻接关系，如图4(b)所示，α取值为3：

α为调节因子，用于调整不合理边的判断阈值，取值范围为[2，4]，默认设为3。

4)在打断不合理边的Delaunay三角网中，针对每个土壤重金属采样点，采用如下公式计算其空间邻域半径R(P_i)：

其中，Mean²(P_i)表示点P_i的2阶邻域内所有边长的平均值；STD(P_j)表示点P_j1阶邻域内边长的标准差。

5)落在每个土壤采样点对应的空间邻域半径范围内的其它采样点即为该土壤采样点的空间邻域。

步骤3、地理加权回归模型构建

1)针对各个土壤采样点重金属浓度值及其独立影响因素，构建如下地理加权回归模型：

y_i＝β₀(u_i，v_i)+β₁(u_i，v_i)x_i1+β₂(u_i，v_i)x_i2+β₃(u_i，v_i)x_i3+ε_i

其中，(u_i，v_i)表示第i个采样点的地理坐标，y_i表示(u_i，v_i)处重金属Cr元素的观测值，x_i1表示(u_i，v_i)处的高程值，x_i2表示(u_i，v_i)处距河流的距离，x_i3表示(u_i，v_i)处距工厂的距离，β_j(u_i，v_i)(j＝0，1，2，3)表示第i个采样点上的第j个回归参数，其空间分布如图5(a)-(d)所示。

2)采用一种自适应权函数：针对每个土壤重金属采样点，根据采样点离回归点的距离远近，把采样点空间邻域中的回归采样点重新排序，空间邻域中距离回归点最远的空间实体的距离d_ij为带宽b。设最近数据点(该点可以是回归点自身)的权重为1，其他数据点的权重按如下公式衰减：

3)土壤采样点重金属浓度变化梯度修复

a)针对每个采样点P，根据其空间邻域采样点，采用下式计算重金属浓度变化梯度：

其中，f(P)表示采样点P的重金属浓度值，D(P，P_i)表示采样点P与邻域采样点P_i间的欧氏距离。

b)令f(P)＝0，分别计算采样点P与其空间邻域采样点P_i间的重金属浓度变化梯度G(P，P_i)，并按升序排列获取序列G(P)，计算序列中位数M(P)；

c)针对任一邻域采样点，采用下式计算重金属浓度变化梯度偏离GD(P_i)，并升序排列获取序列GD(P)：

GD(P_i)＝|G(P，P_i)-M(P)|

d)将邻域采样点按重金属浓度变化梯度偏离划分为大、中、小三个等级，处于最大等级的[(n+1)/3]个实体组成待修复集合R(P)，采用如下公式进行修复：

f_R(P_i)＝M(P)*D(P，P_i)，

说明：修复是暂时的，是在每个采样点地理加权回归的过程中进行，并不改变采样点的固有重金属浓度值。

4)采样下式进行回归参数估计：

其中，W_i为位置(u_i，v_i)处的权重矩阵，W_i＝diag(w_i1，w_i2，...，w_in)，y_R为经过重金属浓度变化梯度修复的因变量矩阵。

步骤4、土壤采样点重金属浓度异常度量与判别

1)针对每个土壤重金属采样点P，重金属浓度Cr观测数据与地理加权回归预测数据间的差异即为采样点P的稳健空间异常度RSOM(Robust Spatial Outlier Measure)，采用下列公式计算：

其中，f(P)表示采样点P的重金属Cr元素浓度观测数据，表示采样点P的重金属Cr元素浓度地理加权回归预测值。

2)获取所有采样点的稳健空间异常度，组成集合S_RsOM＝{RSOM(P₁)，RSOM(P₂)，...，RSOM(P_n)}，异常度偏离均值一定倍数的采样点判别为空间异常，采用下列公式计算：

S_outlier＝{P_i|RSOM(P_i)-μ＞kσ，P_i∈SDB}

其中，μ表示稳健异常度均值，σ表示稳健异常度标准差，k表示判别系数，默认设置为2。土壤重金属Cr元素的异常探测结果如图6所示。

虽然，上文中已经用一般性说明、具体实施方式及试验，对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (7)

1.一种基于多元自适应回归的地理空间异常探测方法，其特征在于，包括如下步骤：

步骤一：针对获取的探测变量及其它地理变量进行数据预处理及回归变量的确定；

步骤二：构建自适应空间邻域；

步骤三：针对探测变量及其它独立地理变量，构建地理加权回归模型；

步骤四：计算所有空间实体的稳健空间异常度，并组成集合，其中异常度偏离均值一定倍数的空间实体判别为空间异常；

所述步骤二包括：

1)根据待检测实体的空间位置构建Delaunay三角网，表达空间实体的粗空间邻接关系；

2)把Delaunay三角网所有边长按升序排列，构成边长序列E，序列中位于上、下四分位数之间所有边长的平均值作为稳健平均边长，利用下列公式进行计算：

其中，Len_ei表示边e_i的边长，m表示Delaunay三角网的边长数目；

3)采用如下公式获取并打断Delaunay三角网中所有不合理边Unreasonable_E，表达空间实体的细空间邻接关系：

其中，α为调节因子，用于调整不合理边的判断阈值，取值范围为[2,4]；

4)在打断不合理边的Delaunay三角网中，针对每个空间实体，采用如下公式计算其空间邻域半径R(P_i)：

其中，Mean²(P_i)表示空间实体P_i的2阶邻域内所有边长的平均值；STD(P_j)表示空间实体P_j的1阶邻域内边长的标准差；n表示数据集中空间实体的数目；

5)落在每个空间实体P对应的空间邻域半径范围内的其它实体即为该实体的空间邻域，记为NN(P)。

2.根据权利要求1所述的基于多元自适应回归的地理空间异常探测方法，其特征在于，所述步骤一，

1)针对获取的探测变量及其它地理变量，检验是否存在缺失，若存在，则进行插补；

2)对探测变量进行均值、中值、范围、标准差、变异系数、峰度或偏度统计；

3)针对其它地理变量，以模型拟合优度的检验统计量F₁为检验统计量，采用逐步回归法确定用来构建模型的回归变量；所述逐步回归法包括向前选择法、向后剔除法或逐步选择法。

3.根据权利要求1所述的基于多元自适应回归的地理空间异常探测方法，其特征在于，所述空间位置为经纬度时，则将其转换为投影坐标。

4.根据权利要求1所述的基于多元自适应回归的地理空间异常探测方法，其特征在于，所述步骤三，

1)针对探测变量及其它独立地理变量，构建如下地理加权回归模型：

其中，(u_i,v_i)表示第i个空间实体的地理坐标，y_i与x_i1,x_i2,…,x_ip表示自变量y和因变量x₁,x₂,…,x_p在位置(u_i,v_i)处的观测值；β₀(u_i,v_i)表示第i个空间实体的回归常数；β_j(u_i,v_i)(j＝1,2,…,p)表示第i个空间实体的第j个回归参数；ε表示随机误差，ε_i～N(0,σ²)，σ为标准差，且协方差Cov(ε_i,ε_j)＝0,i≠j；

2)基于自适应带宽的权函数选择；

3)空间邻域实体专题属性变化梯度修复：

a)针对每个空间实体P，根据其空间邻域，采用下式计算专题属性变化梯度：

其中，f(P)表示空间实体P的探测变量值，D(P,P_i)表示空间实体P与其空间邻域实体P_i间的欧氏空间距离，NN(P)表示空间实体P的空间邻域；

b)令f(P)＝0，分别计算空间实体P与其空间邻域实体P_i间的专题属性变化梯度G(P,P_i)，并按升序排列获取序列G(P,P_i)，计算序列中位数M(P)；

c)针对任一空间邻域实体，采用下式计算专题属性变化梯度偏离GD(P_i)，并升序排列获取序列GD(P)：

GD(P_i)＝|G(P,P_i)-M(P)|

d)将空间邻域实体按专题属性变化梯度偏离划分为大、中、小三个等级，处于最大等级的[(n+1)/3]个空间实体组成待修复集合R(P)，采用如下公式进行修复：

4)采用下式进行回归参数估计：

其中，W_i为位置(u_i,v_i)处的权重矩阵，W_i＝diag(w_i1,w_i2,…,w_in)；X为自变量矩阵；y_R为经过专题属性变化梯度修复的因变量矩阵。

5.根据权利要求4所述的基于多元自适应回归的地理空间异常探测方法，所述基于自适应带宽的权函数选择具体为：

针对每个空间实体，即回归点，根据距离回归点的距离远近，把空间实体的空间邻域中参与回归的空间实体重新排序，空间邻域中距离回归点最远的空间实体的距离d_ij设为带宽b，设距离最近空间实体的权重为1，其它空间实体的权重按如下公式衰减：

6.根据权利要求5所述的基于多元自适应回归的地理空间异常探测方法，所述距离最近空间实体是回归点本身。

7.根据权利要求1所述的基于多元自适应回归的地理空间异常探测方法，其特征在于，所述步骤四，

1)针对每个空间实体P，探测变量观测值与地理加权回归预测值间的差异即为实体P的稳健空间异常度RSOM，采用下列公式计算：

其中，f(P)表示空间实体P的探测变量观测值，表示空间实体P的探测变量的地理加权回归预测值；

2)针对所有空间实体的稳健空间异常度，组成集合S_RSOM＝{RSOM(P₁),RSOM(P₂),…,RSOM(P_n)}，其中异常度偏离均值一定倍数的空间实体判别为空间异常，采用下列公式计算：

S_outlier＝{P_i|RSOM(P_i)-μ＞kσ,P_i∈SDB}

其中，P_i指满足判别条件的空间实体，即空间异常；SDB(Spatial Database)表示空间数据集，包含n个空间实体；μ表示稳健异常度均值，σ表示稳健异常度标准差，k表示判别系数。