CN103869367A - 一种基于共享密度的地震事件时空聚集模式提取方法 - Google Patents

Abstract

本发明涉及一种基于共享密度的地震事件时空聚集模式提取方法。本发明根据时空窗口k邻近关系，识别时空共享邻近地震事件，进而依据时空共享邻近关系，估计时空共享密度，最后依据时空共享邻近关系将高密度地震事件聚集成簇。优点如下：不需要用户设定地震事件时空聚集模式的数目与形态，同时可以提取不同密度的地震事件时空聚集模式，可以从空间和时间耦合的视角发现地震事件的动态演变规律，将地震事件的时空聚集模式可视化表达。

Description

一种基于共享密度的地震事件时空聚集模式提取方法

技术领域

本发明属于时空数据挖掘领域，涉及一种基于共享密度的地震事件时空聚集模式提取方法。

背景技术

地震事件是影响人类社会经济发展与生命财产安全地一种主要自然灾害。依据地震事件地观测资料对地震爆发的规律与机理进行分析与监测是降低地震灾害影响的一种重要途径。目前，随着观测手段的提升，已经积累了大量的地震事件监测数据，地震事件观测的精度与时效性大大提高。地震事件时空聚集模式提取是地震事件分析中一个重要内容，已经成为研究地震事件发展、演化规律的一个重要手段。

现有的地震事件聚集模式提取方法大致可以分为三类：（1）基于时空距离的方法；（2）时空扫描统计方法；以及（3）基于时空密度的方法。基于时空距离的方法首先定义时空实体间的时空距离表达（如临界距离、k阶最邻近距离及时空耦合距离），进而采用假设检验或传统的聚类方法挖掘时空聚集模式。时空耦合距离法常用于地震爆发聚集性探测，难以推广到其它领域。根据临界距离与k阶最邻近距离定义的时空距离常用于检验时空实体的聚集性，而难以获得准确的时空簇。时空扫描统计法首先需要假设数据的概率分布，进而采用一定的扫描窗口（即由空间距离定义半径、时间间隔定义高度的圆柱体），以每个时空实体为中心进行统计检验，确定时空聚集模式。时空扫描统计的方法在疾病爆发热点分析领域进行了成功的应用。然而，时空扫描统计的方法同样无法准确地描述时空聚集模式的范围、大小和位置信息，且扫描窗口的大小难以自动适应时空实体分布密度的变化。基于时空密度的方法是在DBSCAN算法基础上进行的扩展，其核心在于时空邻近域的定义。主要方法包括：在空间邻近的基础上，引入时间半径的概念，定义了圆柱形的时空邻近域；借助Delaunay三角网与时空统计方法定义不规则的时空邻近域。此外，若假定数据的先验分布（如均匀泊松分布），可以将时空事件数据分解为特征与噪声两部分，进而借助密度相连的思想聚集成簇。基于密度的方法可以详细描述时空簇的位置、形状及尺寸信息。然而，基于时空密度的方法不可避免地继承了传统的基于密度方法的局限，即全局性的阈值设置难以识别不同密度的簇，且对于邻近的簇很难区分。综上分析，目前尚缺乏一种能够有效识别不同密度、不同形态地震事件时空聚集模式提取的高效方法。

发明内容

本发明的目的在于提供一种基于共享密度的地震事件时空聚集模式提取方法，主要解决目前缺乏一种能够有效识别不同密度、不同形态地震事件时空聚集模式提取的高效方法。

本发明采取的技术方案是包括以下步骤：

步骤1：针对每个时空实体搜索其时空k邻居；

步骤2：通过时空k邻居，针对每个实体搜索时空共享近邻，并计算其时空共享近邻密度；

步骤3：选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇；

步骤4：将所有能够标记的时空核标记，将所有未加入任何时空簇的时空实体标记为噪声。

本发明的技术特点还在于步骤1具首先确定地震时空事件的时空窗口k邻居，针对每个时空实体搜索其时空k邻居：ST_i(L_i,t_i)表示一个地震时空事件，其中L_i表示空间位置，t_i表示时间点，对于给定的时间窗口ΔT，任一地震时空事件ST_i(L_i,t_i)的时空k邻居为一系列满足下列两个条件的地震时空事件集合（ST_i+1(L_i+1,t_i+1),ST_i+2(L_i+2,t_i+2),…,ST_i+k(L_i+k,t_i+k)），记为STN^k(ST_i)，

①∣t_i+k-t_i∣≤ΔT;

②D_wd(ST_i,ST_i+1)≤D_wd(ST_i,ST_i+2)≤...≤D_wd(ST_i,ST_i+k)，其中D_wd表示欧氏空间距离。

步骤2中时空共享近邻密度的计算方法为时空共享密度计算的具体方法为，首先定义时空直接可达：给定任一地震时空事件ST_i，ST_j∈STN(ST_i)，若N_SNN(ST_i,ST_j)>K_T（K_T为一个大于0的整数），则称时空直接可达，记为

，时空共享近邻密度：一个地震时空事件ST_i的时空共享近邻密度定义为STN(ST_i)中所有满足与ST_i时空直接可达关系的地震时空事件数目，记为STD(ST_i)，时空共享密度计算的具体实现是在Matlab软件中编程实现。

步骤3中首先定义共享时空邻居：给定任一时空事件ST_i，所有与其满足时空直接可达的时空实体构成ST_i的共享时空邻居，记为RSTN(ST_i)，定义时空核：对于任一时空事件ST_i，若其满足STD(ST_i)≥MinPts（MinPts为一个大于0的正整数，STD(ST_i)表示时空共享近邻密度，则称ST_i为一个时空核，定义时空可达：对于一组时空实体ST₁,ST₂,…,ST_i,…,ST_n，且ST_n为一个时空核，若满足

，则称ST₁与ST_n时空可达，记为ST_n→ST₁，定义时空相连：若时空实体ST_i与ST_j对于另一时空实体ST_k，同时满足ST_i→ST_k，ST_j→ST_k，则称ST_i与ST_j时空相连，记为ST_i～ST_j，定义边界时空实体：对于一个非时空核实体ST_i，若其共享时空邻居内存在至少一个时空核，则称之为边界时空实体，最好计算时空簇：对于时空数据集STD，时空簇STC可以描述为满足如下两个条件的非空时空实体集合：

1） $\∀{\mathrm{ST}}_i\∈\mathrm{STD},{\mathrm{ST}}_j\∈\mathrm{STD}:{\mathrm{ST}}_i\→{\mathrm{ST}}_j;$

2） $\∀{\mathrm{ST}}_i\∈\mathrm{STD},{\mathrm{ST}}_j\∈\mathrm{STD}:{\mathrm{ST}}_i~{\mathrm{ST}}_j.$

步骤4中，首先给定聚类参数ΔT、k、K_T与MinPts，选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇，并给予簇内时空核相同标记，直到所有时空核都被扫描且标记后，聚类过程结束，未加入任何时空簇的时空事件标记为噪声，最后，生成的时空簇（或时空聚集模式）将在Matlab软件中可视化展示给用户，并输出最终聚类结果文件。聚类参数k值为12～30；MinPts取值为0.45k～0.55k；k_T取值为0.5k。

本发明的有益效果是能够适应更加复杂的地震事件时空聚集模式提取操作，可以探测复杂结构的时空簇，对噪声点稳健，且具有较高的运行效率。

附图说明

图1是摘要流程图；

图2是时空窗口k邻居识别示意图；

图3是时空共享邻居识别示意图；

图4是参数识别结果；

图5是模拟数据：a：STD₁；b：STD₂；

图6是模拟数据聚类结果：a：STD₁聚类结果；b：STD₂聚类结果；

图7是地震时空数据：a：STD₃；b：STD₄；

图8是地震数据聚类结果：a：STD₄聚类结果；b：STD₄聚类结果。

具体实施方式

以下将结合附图和具体实施例详细说明本发明的技术方案：

本发明如图1所示，主要步骤流程为识别时空共享邻居并定义每个时空事件的时空共享密度后，依据确定的聚类参数识别时空核，选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇，并给予簇内时空核相同标记。直到所有时空核都被扫描且标记后，聚类过程结束，未加入任何时空簇的时空事件标记为噪声。

本发明具体步骤如下：

步骤1：确定地震时空事件的时空窗口k邻居，针对每个时空实体搜索其时空k邻居：ST_i(L_i,t_i)表示一个地震时空事件，其中L_i表示空间位置，t_i表示时间点。对于给定的时间窗口ΔT，任一地震时空事件ST_i(L_i,t_i)的时空k邻居为一系列满足下列两个条件的地震时空事件集合（ST_i+1(L_i+1,t_i+1),ST_i+2(L_i+2,t_i+2),…,ST_i+k(L_i+k,t_i+k)），记为STN^k(ST_i)。

①∣t_i+k-t_i∣≤ΔT;

②D_wd(ST_i,ST_i+1)≤D_wd(ST_i,ST_i+2)≤...≤D_wd(ST_i,ST_i+k)，其中D_wd表示欧氏空间距离。图2表示了地震时空事件ST_i的时空窗口6邻居的识别结果。地震时空事件的时空窗口k邻居识别在Matlab软件中编程实现，其中k与ΔT为用户制定参数。搜索时空窗口k邻居，根据每个地震时空事件的时空窗口k邻居，可进一步定义时空共享邻近事件与时空共享密度。

步骤2：时空共享邻近事件识别，针对每个时空实体搜索时空共享近邻，并计算其时空共享近邻密度；基于步骤1)中每个地震时空事件的时空窗口k邻居，两个地震时空事件（ST_i与ST_j），二者的时空共享邻近事件定义二者其时空k邻居中相同的地震时空事件，记为SNN(ST_i,ST_j)，表达为：

（1）SNN(ST_i,ST_j)=STN^k(ST_i)∩STN^k(ST_j)

图3中表示了一个模拟地震时空事件数据集从水平视角观察的结果，矩形表示时空邻域，以k=8为例对时空共享邻近事件识别进行说明。图3中，除了中心时空事件外，每个邻域（即黑框）内均包含8个时空事件。以地震时空事件ST_i与ST_j为例，4个实心黑色圆表示了ST_i与ST_j的时空共享邻近事件，时空共享邻近事件的数目记为N_SNN(ST_i,ST_j)。空心圆点表示邻域外实体，点状填充圆点便是邻域中心实体，灰色圆点便是邻域内非共享近邻。时空共享邻近事件识别亦是在Matlab软件中编程实现。

时空共享密度计算的具体方法可以通过如下定义进行表达：时空直接可达：给定任一地震时空事件ST_i，ST_j∈STN(ST_i)，若N_SNN(ST_i,ST_j)>K_T（K_T为一个大于0的整数），则称时空直接可达，记为

可见时空直达关系是采用时空共享关系对时空邻近关系的进一步约束。时空共享近邻密度：一个地震时空事件ST_i的时空共享近邻密度定义为STN(ST_i)中所有满足与ST_i时空直接可达关系的地震时空事件数目，记为STD(ST_i)。时空共享密度计算的具体实现亦是在Matlab软件中编程实现，计算获得每个地震时空事件的时空共享近邻密度。基于时空直接可达与时空共享近邻密度的定义，可以进一步对时空簇进行定义。

步骤3：选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇；为了生成任意形状的时空簇，本发明定义了将空间聚类中密度连接的策略定义时空簇，首先如下几个重要的定义：

共享时空邻居：给定任一时空事件ST_i，所有与其满足时空直接可达的时空实体构成ST_i的共享时空邻居，记为RSTN(ST_i)。

时空核：对于任一时空事件ST_i，若其满足STD(ST_i)≥MinPts（MinPts为一个大于0的正整数，STD(ST_i)表示时空共享近邻密度，则称ST_i为一个时空核。

时空可达：对于一组时空实体ST₁,ST₂,…,ST_i,…,ST_n，且ST_n为一个时空核。若满足

，则称ST₁与ST_n时空可达，记为ST_n→ST₁。时空相连：若时空实体ST_i与ST_j对于另一时空实体ST_k，同时满足ST_i→ST_k，ST_j→ST_k，则称ST_i与ST_j时空相连，记为ST_i～ST_j。

边界时空实体：对于一个非时空核实体ST_i，若其共享时空邻居内存在至少一个时空核，则称之为边界时空实体。

时空簇：对于时空数据集STD，时空簇STC可以描述为满足如下两个条件的非空时空实体集合：

（1） $\∀{\mathrm{ST}}_i\∈\mathrm{STD},{\mathrm{ST}}_j\∈\mathrm{STD}:{\mathrm{ST}}_i\→{\mathrm{ST}}_j;$

（2） $\∀{\mathrm{ST}}_i\∈\mathrm{STD},{\mathrm{ST}}_j\∈\mathrm{STD}:{\mathrm{ST}}_i~{\mathrm{ST}}_j.$

步骤4：时空聚集模式生成与可视化。给定聚类参数ΔT、k、K_T与MinPts，选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇，并给予簇内时空核相同标记。直到所有时空核都被扫描且标记后，聚类过程结束，未加入任何时空簇的时空事件标记为噪声。最后，生成的时空簇（或时空聚集模式）将在Matlab软件中可视化展示给用户，并输出最终聚类结果文件。

本发明中聚类参数设置：本发明需要四个输入参数，即ΔT、k、K_T与MinPts。首先需要参数选择。通过分析时空聚集模式的识别过程可以发现，时空密度、时空可达与时空相连关系都是建立在时空直接可达的基础上，而时空直接可达又是对时空k邻近关采用K_T的一种约束，进而可以发现k、K_T与MinPts三个参数间具有密切联系，即k_T与MinPts受k的直接约束，必须满足MinPts<k,k_T<k的条件。

为此，本发明具体实施时将k_T与MinPts分别设为k的一个比值。

为了检验k、MinPts与k_T三个参数的变化规律，这里设计了一组模拟数据，其中包含一个由328个时空实体组成的立方体时空簇和489个噪声（最大时间范围为10天）。为能够识别共享近邻，k的取值不应过小，取值为12～30。首先将k_T设为0.5k，在不同ΔT取值下，分别设置一系列的MinPts取值进行实验分析，错误识别的时空实体数如图3所示。

分析图4所示的结果可以发现：

①当时间窗口设置过小（ΔT=0.32，0.43）或过大（ΔT=2.56，3.14）时聚类结果的准确性较差，但是亦存在一个较长的区间（ΔT=0.64～1.71），误判实体均低于数据总数的5%，因而聚类结果较好，并且比较稳定；

②k值取值过大时（大于30），聚类结果的准确性较差，例如ΔT分别取值0.32、0.43时，误判实体总数约占数据总数的25%-30%；

③MinPts的取值超过0.55k时，大部分结果的准确性迅速降低，MinPts取值介于0.45k～0.55k时，大部分聚类结果具有较高的准确性。基于上述实验结论，本发明对参数设置做出如下约定：

①k值一般取为12～30之间，对于部分小数据集可将k设为5～10之间；

②k_T与MinPts一般设为0.5k左右，本发明设置为0.5k；

③时间窗口ΔT的取值选取多次运行比较稳定的结果作为最终的聚类结果。

本发明的优点：本发明给借助时空共享邻居的概念，给出了一种新的时空密度估计方法。与现有的基于k邻近距离的密度估计及基于固定窗口的密度估计相比，本发明方法其相当于一种归一化的概率密度，即针对不同密度的簇，其内部的共享近邻密度是类似的；在各个簇的边缘及噪声处，共享邻近密度较低。因此，实际中只需要设置一个全局的密度阈值，就可以对不同密度的时空簇以及噪声进行区分，从而有效避免混合了密度分解的复杂性。因此，本发明与传统的经典方法比较，本发明能够适应更加复杂的地震事件时空聚集模式提取操作，可以探测复杂结构的时空簇，对噪声点稳健，且具有较高的运行效率。

本发明将地震事件的时空耦合特征引入时空聚集模式提取，发现不同密度、不同形态的地震事件时空聚集模式所存在的技术问题。提供了一种需要较少的人为干预的时空事件聚集模式提取方法，依据地震事件本身的时空聚集表现，形成一种基于共享密度的地震事件时空聚集模式提取方法，进一步提高地震事件时空聚集模式提取方法的实用性与可靠性。

下面列举具体实施例对本发明进行说明：

实施例1：本发明实验验证中，

1)具体实施步骤如下表1所示：

表1

2)模拟数据实施实例结果：

图5中展示了两组模拟数据，STD₁和STD₂，分别如图5a和b所示。STD₁中共包含951个时空实体，预设了3个互相邻近、密度不同且近似高斯分布的三个时空簇，噪声占实体总数的15%。STD₂是一组更为复杂的时空数据，其中包含了形状复杂、密度各异、大小不同且互相邻近的5个时空簇，噪声占总数据量的12%。应用本发明所述方法提取的时空聚集模式如图6所示，可以看到各个预设的时空簇均被准确识别，由此证明本发明方法可以识别不同密度、不同密度、不同大小、形状复杂的时空簇。

3)实际应用实例：

我国是受地震灾害影响最为严重的国家之一，唐山地震、汶川地震等强烈地震对我国人民生命财产安全及国民经济发展造成了难以估量的损失。挖掘地震数据中的时空聚集模式是识别强烈地震的前震和余震的有力手段，进一步将有助于理解强震的产生与演变机理，服务于地震预报工作。进一步应用本发明方法对两个地震时空数据库（STD₃与STD₄）进行实际应用分析，如图7a与b所示。如图7a展示了1975年1月至1976年8月间，我国陆地及周边区域（70°-125°E,16°-48°N）3级以上的412次地震。如图7b显示了2008年1月至2009年12月间我国陆地及周边区域（73°-120°E,22°-48°N）4级以上的573次地震。

图8中展示了地震时空聚集模式的提取结果。针对STD₃，共发现三个时空簇（如图8a）。STC₁发生在唐山地震之后，其中包含的124次地震可以认为是唐山大地震的余震。STC₂中包含了龙陵地震爆发的两次强震，其中44次地震包含了龙陵地震的前震和余震。STC₃中包含了海城地震，是我国成功预报的一次强烈地震。海城地震的主震爆发在两次较低烈度的地震（M=3.1和4.6）之后，其余的23次地震可以视为海城地震的余震。

针对STD₄，共发现5个较大的时空簇（如图8b）。STC₅是最大的一个时空簇，主震为汶川地震，其中包含的273次地震可以视为汶川地震的余震。分析STC₅的形态可以发现，汶川地震的余震爆发的时空密度是不均匀的，而且时间持续性也较长。如果采用空间聚类方法，这些时间上的渐变特征则无法发现。其他几个时空簇，同样对应了该时间段内我国及周边爆发的强震，可以用来识别各次强震的余震序列。进而，本发明挖掘的地震时空聚集模式可以为地震产生及预测分析研究提供一定的参考。

Claims (6)

1.一种基于共享密度的地震事件时空聚集模式提取方法，其特征在于，包括：

步骤1：针对每个时空实体搜索其时空k邻居；

步骤2：通过时空k邻居，针对每个实体搜索时空共享近邻，并计算其时空共享近邻密度；

步骤3：选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇；

步骤4：将所有能够标记的时空核标记，将所有未加入任何时空簇的时空实体标记为噪声。

2.根据权利要求1所述的一种基于共享密度地地震事件时空聚集模式提取方法，其特征是：所述步骤1具首先确定地震时空事件的时空窗口k邻居，针对每个时空实体搜索其时空k邻居：ST_i(L_i,t_i)表示一个地震时空事件，其中L_i表示空间位置，t_i表示时间点，对于给定的时间窗口ΔT，任一地震时空事件ST_i(L_i,t_i)的时空k邻居为一系列满足下列两个条件的地震时空事件集合（ST_i+1(L_i+1,t_i+1),ST_i+2(L_i+2,t_i+2),…,ST_i+k(L_i+k,t_i+k)），记为STN^k(ST_i)，

①∣t_i+k-t_i∣≤ΔT;

②D_wd(ST_i,ST_i+1)≤D_wd(ST_i,ST_i+2)≤...≤D_wd(ST_i,ST_i+k)，其中D_wd表示欧氏空间距离。

3.根据权利要求1所述的一种基于共享密度地地震事件时空聚集模式提取方法，其特征是：所述步骤2中时空共享近邻密度的计算方法为时空共享密度计算的具体方法为，首先定义时空直接可达：给定任一地震时空事件ST_i，ST_j∈STN(ST_i)，若N_SNN(ST_i,ST_j)>K_T（K_T为一个大于0的整数），则称时空直接可达，记为

，时空共享近邻密度：一个地震时空事件ST_i的时空共享近邻密度定义为STN(ST_i)中所有满足与ST_i时空直接可达关系的地震时空事件数目，记为STD(ST_i)，时空共享密度计算的具体实现是在Matlab软件中编程实现。

4.根据权利要求1所述的一种基于共享密度地地震事件时空聚集模式提取方法，其特征是：所述步骤3中首先定义共享时空邻居：给定任一时空事件ST_i，所有与其满足时空直接可达的时空实体构成ST_i的共享时空邻居，记为RSTN(ST_i)，定义时空核：对于任一时空事件ST_i，若其满足STD(ST_i)≥MinPts（MinPts为一个大于0的正整数，STD(ST_i)表示时空共享近邻密度，则称ST_i为一个时空核，定义时空可达：对于一组时空实体ST₁,ST₂,…,ST_i,…,ST_n，且ST_n为一个时空核，若满足

，则称ST₁与ST_n时空可达，记为ST_n→ST₁，定义时空相连：若时空实体ST_i与ST_j对于另一时空实体ST_k，同时满足ST_i→ST_k，ST_j→ST_k，则称ST_i与ST_j时空相连，记为ST_i～ST_j，定义边界时空实体：对于一个非时空核实体ST_i，若其共享时空邻居内存在至少一个时空核，则称之为边界时空实体，最好计算时空簇：对于时空数据集STD，时空簇STC可以描述为满足如下两个条件的非空时空实体集合：

1） $\&ForAll;{\mathrm{ST}}_i\&Element;\mathrm{STD},{\mathrm{ST}}_j\&Element;\mathrm{STD}:{\mathrm{ST}}_i\&RightArrow;{\mathrm{ST}}_j;$

2） $\&ForAll;{\mathrm{ST}}_i\&Element;\mathrm{STD},{\mathrm{ST}}_j\&Element;\mathrm{STD}:{\mathrm{ST}}_i~{\mathrm{ST}}_j.$

5.根据权利要求1所述的一种基于共享密度地地震事件时空聚集模式提取方法，其特征是：所述步骤4中，首先给定聚类参数ΔT、k、K_T与MinPts，选取一个未标记的时空核，搜索所有时空可达与时空相连实体标记为同一个时空簇，并给予簇内时空核相同标记，直到所有时空核都被扫描且标记后，聚类过程结束，未加入任何时空簇的时空事件标记为噪声，最后，生成的时空簇（或时空聚集模式）将在Matlab软件中可视化展示给用户，并输出最终聚类结果文件。

6.根据权利要求5所述的一种基于共享密度地地震事件时空聚集模式提取方法，其特征是：所述聚类参数k值为12～30；MinPts取值为0.45k～0.55k；k_T取值为0.5k。

Images