1 引言     

空间异常探测是空间数据挖掘领域一项重要研究内容[1-2]，旨在从海量空间数据中挖掘小部分偏离普遍模式的空间实体。这些异常实体通常蕴含着意想不到的知识，可能代表着地理现象或地理过程的特殊发展规律。近年来，空间异常探测受到学者们的广泛关注，并已应用于地质灾害监测、环境监测与保护、公共卫生、遥感图像数据处理等领域，具有重要的研究价值。

异常（Outlier，亦称离群点）的定义最初来源于Hawkins的研究工作，给出异常的本质性定义即“严重偏离其他对象的观测数据，以至于令人怀疑它是由不同机制产生的”[3]。Shekhar等人[4]将传统异常在空间上进行了有效扩展，指出空间异常是指专题属性与其空间邻近域内的其他参考实体的专题属性显著不同，而在整体数据范围内差异可能不明显的空间实体。空间异常通常根据空间属性（即位置）确定空间邻近关系，根据专题属性确定异常程度。现有的空间异常探测方法可分为：(1)基于统计的方法[3]；(2)基于图形的方法[5]；(3)基于距离的方法[3][6][7]；(4)基于密度的方法[8][9]；(5)基于聚类的方法[10][11]；(6)基于模型的方法[12-14]。基于统计的方法要求数据服从一定的分布，适用性不强。基于图形的方法主要是采用可视化的方法（如变量云和散点图等）来查找空间异常点，利用人眼进行主观性的识别使这类方法已很少使用[16]。基于距离的方法采用专题属性值与其空间邻域的专题属性的均值或中值的差值来度量实体的异常程度，继而通过统计测试的方法识别异常，这类方法采用全局策略，仅适合探测全局异常，不适合探测局部异常实体。基于密度的方法采用不同的密度估计策略度量实体的局部异常度，异常度较大的实体识别为空间异常。这类方法依赖于空间邻近域的选择，缺乏一定的准确性和稳定性。基于聚类的方法旨在发现空间簇，探测空间异常的能力有限，且探测结果依赖于聚类算法的选择。基于模型的方法采用统计模型、机器学习模型等数学工具进行空间异常探测。这类方法需满足模型的假设条件，如数据分布，这在实际运用中很难准确获得，可能使得探测结果偏离实际情况。

对现有方法进行分析归纳可以发现，当前方法存在的主要问题是：（1）现有方法多认为所有空间实体间具有等同的相关性，而实际上空间实体间具有局部的相关性，在全局更呈现异质特性。如在一个大区域内，经常出现专题属性差异较明显的几个子区域，在探测时不加以区分，一方面可能导致空间邻域内实体不满足相关性假设，导致异常度量的偏差；另一方面可能造成由于局部出现较多异常实体，使得一些在整体上不明显而局部差异明显的空间异常难以被发现。如图1所示（Chawla & Sun（2006）采用的模拟数据集），Ⅰ、Ⅱ、Ⅲ各自区域内的实体间具有较强的相关性，而不同区域内实体间的相关性较弱，甚至是独立的。根据SLOM法共探测出5个异常，如浅绿色标识。然而，在区域Ⅱ、Ⅲ中三个浅蓝色背景标识SLOM值明显偏大，很可能是空间异常。但由于Ⅰ区出现较多明显的异常点，导致Ⅱ、Ⅲ区域内的异常被忽略的现象。因此，需要加以区分。（2）现有方法多采用空间属性确定空间邻近关系，专题属性度量异常度，缺乏一种耦合两类属性的异常度度量的指标。

针对上述问题，本文引入数据场的概念，结合空间聚类技术，发展了一种基于场论的空间异常探测方法，方法流程如图2所示。

(a) 专题属性空间分布(b) 空间异常度SLOM值

图1 模拟空间数据集（Chawla & Sun, 2006）

图2 基于场论探测空间异常算法流程

2基于场论的空间异常探测方法

本文空间异常探测方法主要分为4个步骤：（1）采用空间聚类技术获取空间簇，针对各空间簇生成Delaunay三角网，进而约束获得空间邻近域；（2）采用专题属性变化梯度修复策略处理空间邻域中潜在异常的影响；（3）引入空间数据场概念，采用似高斯势函数度量空间异常度；（4）针对各空间簇分别统计识别空间异常，并进行评价与分析。下面对每个步骤进行详细阐述。

2.1 空间邻近域构建

本文采用一种基于多约束的自适应空间聚类方法ASCDT[17]进行空间聚类。该方法通过施加不同层次、不同类型的约束可以适应空间数据分布不均匀、空间簇形态各异、位置邻近等复杂情况下的聚类，且输入参数较少，自适应能力强。通过聚类，将所有空间实体划分为若干空间相关性较强的空间簇。针对各空间簇，借助Delaunay三角网来构建空间邻近关系。Delaunay三角网是一种满足最大最小角特性、外接圆特性和唯一性的三角剖分，能自然的反映空间实体间的邻接关系。但从图3(b)可以发现，原始Delaunay三角网在边界和空洞处的边长明显偏长，如实体A与B，C与D空间邻近显然是不合理的。文献[18]通过实验证明，不合理边可采用删除超过平均边长一定倍数的边的策略有效移除。本文亦采用一种稳健平均边长处理不合理边。

定义1稳健平均边长：给定各空间簇C，C中所有实体构成Delaunay三角网的N条边构成边长集合E，E中所有边长按升序排列，序列中位于上四分位数和下四分位数之间所有边长的均值称为稳健平均边长，记为RAE(E)，表示为：

（1）

式中：Q1为上四分位数，Q3为下四分位数，n表示上下四分位数之间所有边的数量。

定义2不合理边：边长集合E中，与稳健平均边长相比明显偏大的边定义为不合理边，所有不合理边构成不合理边集合EIC，表示为：

 （2）

式中：为调节因子，用于调整不合理边的判断阈值，其取值范围可通过启发式策略确定[18]，本文通过大量实验发现，取值[2,3]时较为适宜。

在空间簇Delaunay三角网中打断不合理边后，可获得每个实体的空间邻近实体。如图3（c）所示，经打断操作后（=2）空洞和边界处的不合理边被有效移除，据此建立的实体间的空间邻近关系更为合理、稳定。

定义3空间邻域：对于空间簇中任一空间实体Pi，通过打断不合理边Delaunay三角网的边直接相连的空间实体构成Pi的空间邻域NN(Pi)。

(a) 空间聚类结果 (b) 各空间簇原始Delaunay三角网 (c) 约束后结果图

图3 空间邻近域构建

2.2 专题属性变化梯度修复

在度量空间异常度时，需有效消除空间邻近实体中异常值的影响，借鉴文献[10]的研究策略，本文采用一种专题属性变化梯度修复策略进行消除。

定义4专题属性变化梯度：给定空间实体P，其空间邻域为NN(P)={X1,X2,…,Xn}，f(Xi)表示实体Xi的专题属性值，专题属性变化梯度定义为空间实体P与其某一邻近域实体Xi的专题属性差值的绝对值与二者间欧氏距离的比值，记为G(P,Xi):

 （3）

针对任一空间实体P进行邻域异常值修复的步骤为：

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

（2）针对任一空间邻域Xi，计算专题属性变化梯度偏离GD(Xi)，，按升序排列获取序列GD(P)。

（3）将邻域实体按专题属性变化梯度偏离划分为大、中、小三个等级，处于最大等级的[(n+1)/3]个实体组成待修复集合R(P)。进而，采用专题属性变化梯度序列的中位数M(P)进行修复：

 （4）

式中，表示专题属性修复值。

异常值修复策略旨在消除空间邻域内潜在异常值对度量异常度的影响，保证空间邻域内实体间专题属性的局部平稳性假设。且这种修复是暂时的，仅在空间异常度量的过程中进行，并不改变实体的固有专题属性值。

2.3 空间异常度量

现有的空间异常探测方法多数是从几何的角度出发，采用各种距离指标度量异常度，缺乏实际的物理意义。且缺乏一种耦合空间属性和专题属性的度量指标。本文受物理学中场论思想的启发，从数据场的角度对空间异常探测进行解释，采用似高斯势函数度量空间异常度。

场的概念最早是由英国物理学家法拉第于1837年提出，用于描述物质粒子间的非接触相互作用。随着场论思想的发展，人们将其抽象为一个数学概念，用于描述某个物理量或数学函数在空间内的分布规律，分为矢量场和标量场。势场是一种重要的标量场，势函数是位置或距离的函数，可以叠加。数据场是指数据通过数据辐射将其数据能量从样本空间辐射到整个母体空间，接受数据能量并被数据辐射所覆盖的空间[19]。如图4描述的是借助Voronoi图和Delaunay三角网的空间数据的凝聚场。实验研究表明，短程场作用更有利于揭示数据分布的凝聚特性，因此数据场的作用范围必须在有限范围内迅速衰减[21]。高斯函数则可以满足迅速衰减的特性，在充分考虑空间数据自身的特点以及空间数据辐射的特性，本文在基于Delaunay三角网获取空间邻域的基础上，采用似高斯势函数来定义空间实体的空间异常度。

定义5空间异常度：给定空间实体P，其空间异常度SOM表达为：

，

 （5）

式中：NN(P)为实体P的空间邻域，为空间邻域数目，表示实体P与Xi间的欧氏距离，表示实体P与Xi间的专题属性距离，本文采用归一化的闵氏距离，为辐射因子。

图4 数据场—凝聚场示意图

2.4 空间异常识别

针对空间聚类获取的若干空间相关性较强的空间簇，采用统计判别法进行识别，分别探测异常。

定义6空间异常：给定一空间簇数据集SCDB，其空间异常集SCOs为：

,  （6）

式中：为异常度平均值，为标准差，k为调节因子。通常k取值2或3，文献[22]通过实验结果证明，当k=1.645时判断异常的结果较合理可靠。因此，本文中k取值1.645。所有簇的异常集SCOs的集合构成空间数据集SDB的空间异常集SOs。

空间聚类后，各空间簇异常探测的样本减少。直接采用式（6）可能降低异常探测的稳定性。因此，针对小样本数据（样本数小于30）采用中位数Median和中位数绝对偏差MAD替代和。其中，，。

（7）

2.5 算法描述和复杂度分析

基于以上相关定义，基于场论的空间异常探测方法FTSOD可描述为：

输入：包含N个实体的空间数据集SDB

输出：空间异常数据集SOs

（1）对空间数据集SDB进行空间聚类，获取若干空间相关性较强的空间簇SC；

（2）针对各空间簇分别生成Delaunay三角网，并采用稳健平均边长移除不合理边，获取空间邻域；

（3）采用专题属性变化梯度修复策略消除空间邻域内潜在异常值的影响；

（4）采用异常度量指标计算各空间实体的空间异常度SOM；

（5）针对各空间簇分别采用统计判别准则识别异常，获取空间异常集SOs。

通过对FTSOD算法进行分析，对于包含N个实体的空间数据库，非空间属性维数为d，其复杂度主要包括：ASCDT算法的复杂度为O(NlogN)；构建Delaunay三角网、打断不合理边并获取空间邻域的复杂度约为O(NlogN)+O(6N)+O(6N)；专题属性归一化复杂度为O(dN)，异常度计算的复杂度为O(NlogN)，异常判别的复杂度为O(NlogN)。于是，该算法的复杂度为：O(NlogN)+O(NlogN)+O(6N)+ O(6N)+ O(dN)+ O(NlogN)+ O(NlogN)。当时，算法复杂度近似为O(NlogN)。

3 实验与分析

设计两组实验来验证本文所提出的FTSOD算法的正确性。实验1采用模拟数据，模拟数据是由文献[23]中三组经典数据库的一部分组成，进一步添加了一维专题属性，专题属性服从正态分布，其空间分布及部分异常局部放大如图5所示。实验2采用华南某市土壤重金属Cr浓度监测数据，进行实际应用分析。2组实验结果均与经典的SLOM[8]算法进行比较。

3.1 模拟算例分析与比较

采用的模拟数据集SDB具有不同形状、不同密度的空间簇分布，且专题属性值服从正态分布。实验结果中所得到的空间异常点用“X”表示。

图5 模拟数据集SDB空间分布            图6 模拟数据空间邻近域构建

(a)  FTSOD算法空间异常度分布图(b) 异常点空间分布

图7 FTSOD算法异常探测结果

(a) SLOM算法空间异常度分布图(b) 异常点空间分布 k=7 N=25

图8 SLOM算法异常探测结果

图7、图8分别给出了本文方法FTSOD和SLOM算法异常探测结果及异常度空间分布图。比较探测结果可以发现：（1）本文提出的FTSOD算法可以正确识别空间异常点；（2）SLOM算法出现了部分误判和漏判的现象，这是因为在度量异常度时没有有效消除邻域内潜在异常的影响，使得异常度量不准确，且从全局去识别异常，没有顾忌实体的空间分异特性。

3.2 实例分析与比较

本文采用我国华南某市环保数据中土壤重金属Cr浓度监测数据验证本文算法的可行性。采样点空间分布如图9(a)所示，采样点共103个。探测结果与SLOM算法结果进行比较。依据本文算法的基本步骤，首先进行空间聚类。图9(b)为ASCDT算法聚类结果，所有空间实体均加入空间簇中，共获得11个空间簇。针对各空间簇分别生成Delaunay三角网并处理不合理边，空间邻域结果图如图9(c)所示。当空间簇个数小于等于6时，簇中所有空间实体互为空间邻域。进而采用本文提出的基于场论的方法在各个空间簇中探测空间异常，探测结果如图10(a)所示。

(a)采样点空间分布(b)空间聚类结果 (c)空间邻近域构建

图9 实际数据

(a)FTSOD算法异常探测结果(b)SLOM算法异常探测结果  k=7 N=16

图10 空间异常探测结果

SLOM算法异常探测结果如图10(b)所示。通过与SLOM探测结果进行比较分析，可以发现，SLOM算法虽能顾及局部特性，但仅能发现整体上异常度较大的实体，对于局部异常现象，如簇9中的86号采样点（表1所示），在整体上异常程度不明显，但其专题属性严重偏离空间邻域的其他实体，表现为局部异常现象。且异常识别中的小样本识别策略亦十分稳健。

表1 空间簇9中实体专题属性值

表2 空间异常采样点土地使用类型