传统基于划分^[7]、层次^[8]、密度^[9]和格网的聚类方法能灵活定义聚类形状，并对海量和高维数据进行处理。然而这些聚类方法主要用来发现数据中自然存在的聚类，对于异常的识别，仅以附属的形式给出。并且不对数据在空间上是随机分布还是聚集分布进行判断，不能给出聚类在空间分布上的显著性程度，也不能很好地综合空间和非空间因素进行分析。20世纪80年代，首先在流行病学领域出现了相关的空间异常聚类模式的研究与应用。一些空间聚类方法被提出，并用来估计疾病的危险程度、识别显著的危险区域^[10]。这些方法主要是基于空间扫描窗口技术进行的。

典型的扫描统计方法有Openshaw等提出的地理分析机^[11]，Kulldorff提出的空间扫描统计方法SaTScan^{[12, 13]}。SaTScan不仅能识别出聚类的准确位置，并且能计算聚类的显著性程度，成为目前普遍使用的空间聚类分析手段。其原理是：通过扫描以给定位置(如规则格网的中心)为圆心、可变半径的圆形或椭圆窗口，遍历整个数据集空间。每次落入扫描窗口内的案例点个数可视为一个贝努利(或泊松)计数过程。当扫描窗口的大小固定时，其扫描统计量是在整个区域中窗口所含案例的最大个数。通过使用局部似然函数来估计扫描统计量的对数似然度比率(log likelihood ratio, LLR)的极大值，然后利用蒙特卡罗模拟对扫描统计量进行显著性检验，最终判别扫描统计量对应的窗口内是否存在异常聚类。

在时空异常聚类挖掘研究方面，扫描窗口变成圆柱形或椭圆柱形。按扫描窗口的时间区间可将时空异常聚类挖掘方法分为回顾型^[14]和前瞻型^[15]两种。这些方法也使用了基于Poisson分布的时空扫描统计量，但此统计量在计算过程中需要提供各区域的背景控制数据(如背景人口)。由于精确的背景控制数据很难获取到，并且数据可能会随时间而不断变化，Kulldorff等^[16]提出一种新的时空重排扫描统计量，认为当圆柱内部的案例个数相比总的案例个数非常小时，可采用广义似然函数来替代原似然函数进行聚类识别和显著性检验，并利用重排算法避免使用控制数据。

以上时空异常聚类方法采用规则扫描窗口，无法准确获取异常聚类的真实形状。由于算法所需计算的不规则区域的组合太多，是典型“NP-难”问题，不规则形状的时空异常聚类挖掘方法研究尚无好的解决办法。目前，仅在空间异常聚类挖掘领域，出现了3种化简搜索空间，求解最可能的空间异常聚类区域的方法。

1) 基于降维策略的方法^{[17, 18]}。该分法是将二维空间点的空间分布按照一定规则(如最近邻居)转换成一维的点过程；针对一维点过程数据进行分析，获得一维数据的异常聚类序列；最后再映射到二维空间点的空间异常聚类区域。此类方法受一维转换顺序的影响很大，准确性不高，只能用来做初级分析。

2) 基于压缩解空间的方法^{[19, 20]}。利用某些标准(如点密度、K邻近)将解空间缩小到容易计算的范围，再针对缩小后的解空间进行扫描统计计算，得到一个准最优解。此类方法受压缩标准的影响很大，容易陷入局部最优解，而无法获得最大的目标函数值。

3) 智能最优化搜索方法。采用的是启发式的搜索策略，利用限定约束条件或限定次数的重复迭代过程以使搜索收敛到全局最优解或准最优解，并且在搜索过程中加入随机概率扰动机制以避免目标函数陷入局部最优解。Duczmal等提出基于模拟退火的空间扫描聚类方法SaScan^[21]和基于遗传算法的空间扫描聚类方法GaScan^[22]。Janeja等^[23]提出一种基于随机行走的任意形状空间扫描统计方法(free-form spatial scan statistic, FS³)。文献[24]提出了基于蚁群最优化的空间扫描聚类方法AntScan，在处理不规则形状空间聚类的局部弱连通区上取得了比SaScan和GaScan更好的效果。文献[25]在AntScan的基础上提出多区域不规则形状空间异常聚类方法ACOMCD。方法依据限定聚类大小参数的SaTScan聚类结果获取多个候选的聚类区域位置，然后再利用AntScan对各个候选聚类区域的形状、大小进行准确识别，能够很好地识别出多个共存的不规则形状空间异常聚类区域；并且通过构造聚类区域间的禁忌表解决了在对数似然度函数LLR最大化的过程中导致的异常聚类区域形状不受控制、容易分叉的问题。

本文不规则形状时空异常聚类模式挖掘方法stAntScan是将不规则形状的异常聚类模式挖掘和时空异常聚类挖掘方法相结合，开展不规则形状的时空异常聚类区域识别研究。其方法分为以下4个步骤。

1) 研究区的时空格网划分。将研究区的时空范围按格网划分，如空间维按公里格网划分，时间维按天划分，形成一个三维的时空格网(图 1(a))。每个格网单元格就是一个最小的计算单元，首先统计各单元内的案例点个数和背景控制数据，然后扫描统计方法对当前格网及其时空邻近的格网内、外区域进行局部似然度的计算和异常聚类的显著性检验，获取时空异常聚类区域。

图  1  stAntScan方法示意图

Figure 1.  Sketch Map for stAntScan Method

2) 创建三维时空邻接单元格。时空邻接单元格的设定关系到时空异常聚类区域的扩张方式，并影响最终的聚类大小和形状。在这里，定义了时空立方体内格网单元格的26方位时空邻近单元格(图 1(a))。依据此26方位时空邻接单元格，可构建研究区时空立方体的邻接矩阵，作为时空异常聚类区域扫描的路径参考。

3) 改进的蚁群最优化时空扫描统计方法。在构建了研究区的时空立方体邻接矩阵后，利用多区域不规则形状空间异常聚类方法ACOMCD来识别多个共存的不规则形状时空异常聚类模式。与原ACOMCD方法不同的是，蚂蚁在进行路径选择时，除了空间邻近的单元格外，还可能选择在三维时空内的邻近单元格。此外，由于增加了时间维而导致的计算量大幅增加，本文对原ACOMCD方法进行了如下的改进策略。

(1) 参考SaTScan方法，增加聚类规模限制参数(ε=50%)。此参数用于在聚类区域的单元格扩张过程中，限制单元格的总长度。

(2) 设置蚁群默认的蚂蚁个数，并缩小蚂蚁初始位置的随机选择范围。蚁群默认的蚂蚁个数设置为案例点数大于平均值的单元格总数。并且蚂蚁初始位置不再是从全部单元格中随机选取，而是从以上大于平均值的单元格中随机选择。此项改进的依据是假设时空异常聚类区域必然包含一个大于平均值的单元格，并假设蚂蚁从该单元格出发去寻找最佳路径。这一假设符合单元格聚类的最基本要求，并且根据聚类区域的连通性，通过此单元格可以获取到完整的时空异常聚类区域。

(3) 为避免蚁群最优化过程提前收敛，增加了一个计数参数nochange。若蚁群重复获得大于nochange次相同的最大LLR，则保留当前的最优蚁群路径作为候选最佳路径，对原信息素矩阵进行一次初始化后再进行蚁群迭代，直到达到总的迭代次数为止。

以上改进使得stAntScan具有了更好的适应性和鲁棒性，但是算法的时间复杂度并没有因此而增加，和基本的蚁群算法一样^[26]。改进后的stAntScan运行需要设定以下5个参数：迭代次数iterationTime、信息素重要程度Alpha、启发式因子重要程度Beta、信息素蒸发系数Rho和信息素增强系数Q。它们的取值直接关系到蚁群最优化的搜索效果和收敛速度。这里不对如何设定这些参数进行详细的论证说明，因为其分析过程与传统的蚁群算法中参数设定是类似的，相关研究可参考文献[26]。经反复实验验证，5个参数的设定均可设定为默认值。其中iterationTime=单元格总个数/候选聚类区域个数，Alpha=1，Beta=1，Rho=0.9，Q=1。候选聚类区域个数由限定参数的SaTScan方法获得。从后面的实验效果看，stAntScan的收敛结果很稳定。

4) 时空异常聚类的可视化展示。时空异常聚类的结果需在三维空间对时空异常聚类的区域位置和形状展示。如图 1(b)所示, 采用三维坐标系，并以时间维分层的形式展示不同时间段内异常聚类的区域范围。此外，还可采用二维坐标系，并以分时间段的二维图层展示异常聚类的区域范围。此种方式对异常聚类的时空显示效果上不如前一种方式，但可以更好地显示每个时间段内异常聚类区域的形状。

本文将使用模拟的时空格网数据和真实微博签到数据对stAntScan方法进行实验验证，并将实验结果与经典时空扫描统计方法SaTScan进行比较分析。

实验模拟了一个持续4 d的时空聚类事件，来验证stAntScan对不规则时空异常聚类模式的识别能力。同时，实验选取了基于回顾型时空重排模型的SaTScan方法做为比较。如图 1(b)所示，模拟数据的空间范围包含20×20共400个二维格网区域，时间区间为0~3。数据给定聚类区域外格网的案例数与背景数据分别为10和40(共1 540个单元格)；给定聚类区域内格网的案例数与背景数据为20和40(共60个单元格)。聚类的时空范围是以三维格网中间4×4个单元格为中心(用红色星形标识)，随时间推移分别向北、东、南、西4个方向扩张，并且区域面积逐渐增大。为检验stAntScan对不同时间长度的聚类识别效果，实验分两组，分别对前3天1 200个、全部4天1 600个格网单元格进行扫描统计，获取LLR最大的格网单元格区域，并利用蒙特卡罗模拟数据进行聚类的显著性P值计算(P=rank/(n+1), rank是实际聚类LLR值在n次蒙特卡罗模拟数据迭代的LLR值序列中的逆序排列编号，P < 0.01则认为实际聚类的显著性很高)，实验参数如前所述。

将实验中每次蚁群迭代的似然度与实际聚类的似然度值进行比较，可获得一组逐渐逼近于1的数值，该组数值反映了stAntScan收敛的过程。如图 2展示了stAntScan在第二组数据蚁群迭代过程中LLR和最佳路径长度的变化情况。可以看到，随着迭代次数的增加，LLR最终收敛到1，而路径长度最终收敛到实际聚类大小。

图  2  stAntScan模拟数据的收敛过程

Figure 2.  Simulated Data Experiments for stAntScan

对每组实验各重复100次，实验结果表明stAntScan总能识别出完整的聚类区域。3 d和4 d模拟数据的平均收敛次数分别为103次和135次，由此证明了stAntScan收敛效果稳定。

与stAntScan相比，SaTScan方法也能识别出模拟数据的聚类所在区域，但识别的准确率不高。如图 3所示，灰色填充区域是SaTScan识别的聚类区域，红色星形标注的区域是实际的聚类区域。

图  3  4 d模拟数据的SaTScan实验结果

Figure 3.  Simulated Data Experiments for SaTScan

利用假正率、假负率和总准确率三项指标对stAntScan和SaTScan进行的准确性计算结果如表 1所示。可以发现，随着时空数据的增加和聚类形状越来越不规则，SaTScan的准确率在逐渐降低。而stAntScan方法准确性很好，可以不受时空数据增加和聚类形状大小的影响，收敛到实际聚类。

真实数据的实验采集了北京市区五环内2014-06-28~2014-06-30共3 d的新浪微博签到数据。在空间尺度上，研究区划分为593个公里格网；在时间尺度上，按天划分。共获得1 779个格网单元格。然后，以格网内签到数据的密度(签到数/格网面积)为案例点/控制点的比率值，识别签到密度显著高的时空异常聚类区域。原始签到数据逐日的密度分布如图 4所示。可以看到，中关村、工体、朝阳门、金台夕照和国贸五处在这期间均是签到密度最高的一类区域。五道口、什刹海、西单、前门和王府井等5处在06-28和06-29是签到密度最高的一类。

图  4  北京市区新浪微博签到数量时空密度分布图(06-28~06-30)

Figure 4.  Spatio-Temporal Density Distribution of Weibo Check in Data in Beijing (06-28~06-30)

采用时空重排模型的SaTScan方法识别出9个显著性P值 < 0.01的时空异常聚类区域，包含75个格网单元，并且全部9个聚类的时间段全都是在06-30，如图 5(a)所示。对比图 4的签到密度分布图，可发现SaTScan的聚类结果并未能包含密度最高的朝阳门和工体两处，对北四环及北五环附近的两处聚类的区域范围识别不够好，并且西南部和东南部有两处聚类的实际签到数并不多。06-30中的9个显著聚类区域的大小及似然比率值LLR如表 2所示。

图  5  SaTScan和stAntScan的时空聚类效果图

Figure 5.  Spatio-Temporal Cluster Results for SaTScan and stAntScan

stAntScan方法仅识别出1个时空连通的不规则形状聚类区域，包含06-28~06-30共261个格网单元格，如图 5(b)、5(c)、5(d)所示。并且，聚类区域的对数似然比率值达到了24 410，聚类显著性远高于SaTScan的聚类结果。将聚类结果与图 3的签到时空密度分布图进行对比，可发现stAntScan对每天的热点聚类区域都能很好地识别，聚类区域的大小、形状与实际签到密度分布图中“很高”和”高”的两类区域匹配程度高，连通性强，较好地反映签到密度的时空分布及变化趋势。不足之处是，stAntScan方法受到前一步骤SaTScan候选聚类区域的影响，对于西南部和东南部两处的候选聚类区域尚不能做到很好的甄别。

时空异常聚类模式是一种重要的空间点模式，其挖掘结果可以及时、准确地反映研究对象的异常分布和时空演化模式。传统的扫描统计方法虽然能有效地获取空间和时空异常聚类模式，但由于扫描窗口通常采用圆柱、椭圆柱等固定形状，其所能识别的聚类形状和范围也受到限制，准确性降低。本文提出一种新的时空扫描统计方法stAntScan，该方法利用蚁群最优化方法实现了对不规则形状的时空异常聚类区域的有效识别，准确性较经典的SaTScan方法高。

与传统的空间扫描统计方法所不同的是，时空扫描统计处理的数据量随时间推移而成几何级数的增长。这对时空扫描统计方法的执行效率构成了很大挑战。并且，在不规则聚类形状的扫描过程中，时空邻近单元格的数量在时空维度上增长很快，增大了聚类收敛的随机性。更深入的研究应该从提高效率和减少收敛随机性两个方面展开。