CN114167489A

CN114167489A - 一种地震前兆特征筛选分析方法、系统及存储介质

Info

Publication number: CN114167489A
Application number: CN202111504790.3A
Authority: CN
Inventors: 刘一宾; 雍珊珊; 王新安; 谢锦汉; 张馨宝; 马一中; 鲍振宇
Original assignee: Peking University Shenzhen Graduate School
Current assignee: Peking University Shenzhen Graduate School
Priority date: 2021-12-10
Filing date: 2021-12-10
Publication date: 2022-03-11
Anticipated expiration: 2041-12-10
Also published as: CN114167489B

Abstract

本发明涉及一种地震前兆特征筛选分析方法、系统及存储介质。本发明依据实际监测精度的要求在待监测地区设置多个监测点，对监测点处不同时间的地震前兆信号的监测数据和监测时间内发生地震事件的时间进行监测，并对地震事件发生前预设时间段的监测数据进行统计提取其异常特征，并分析异常特征与地震事件的相关性，不需要对被测点附近地表和建筑产生任何损害，以为地震预测模型的构建提供高质量的数据输入，进而能够显著提高地震预测的准确性。并且能通过固定监测台站对待监测地区进行长期监测，具有较好的执行可行性。

Description

一种地震前兆特征筛选分析方法、系统及存储介质

技术领域

本发明涉及数据处理技术领域，特别是涉及一种地震前兆特征筛选分析方法、系统及存储介质。

背景技术

地震是一种对于人类社会有着巨大危害的自然灾害，对人类社会造成了巨大的社会经济损失。中国地震活动频繁，是世界上受地震灾害最为严重的国家之一。为了尝试破解地震监测预测的难题，目前研制了AETA多分量地震监测系统，希望通过低成本、大区域、高密度的设备布设建立完善的监测网络，捕捉比较一致的前兆异常信号，服务于地震预测。虽然该系统目前在现有特征数据方面已经取得了丰富的研究成果，但现有特征数据存在有效信息过少以及部分特征不够合理的问题，进而导致地震预测的准确度不高。

发明内容

本发明的目的是提供一种地震前兆特征筛选分析方法、系统及存储介质，以提高地震预测的准确性。

为实现上述目的，本发明提供了如下方案：

一种地震前兆特征筛选方法，包括：

在监测区域中布设监测点；

对待监测区域中所述监测点处的地震前兆信号进行监测，获得地震前兆原始数据序列；所述地震前兆原始数据序列包括：待监测区域中监测点处不同时间的地震电磁扰动信号；

对所述地震前兆原始数据序列进行预处理，获得以小时为粒度的特征序列；

对以小时为粒度的特征序列进行异常检测得到异常序列。

优选地，所述对所述地震前兆原始数据序列进行预处理，获得以小时为粒度的特征序列，具体包括：

对所述地震前兆原始数据序列进行滤波处理，获得滤波后的地震前兆原始数据序列；

基于滤波后的地震前兆原始数据序列确定以分钟为粒度的特征序列；以分钟为粒度的特征序列包括以分钟为粒度的ULF特征序列和功率谱密度特征序列；

基于以分钟为粒度的特征序列确定以小时为粒度的特征序列。

优选地，所述基于滤波后的地震前兆原始数据序列确定以分钟为粒度的特征序列，具体包括：

以设定时间为窗口，在滤波后的地震前兆原始数据序列上滑动，获得ULF特征序列；

使用时频分析方法对滤波后的地震前兆原始数据序列进行时域分析，获得一个以分钟为粒度的功率谱密度特征序列。

优选地，所述基于以分钟为粒度的特征序列确定以小时为粒度的特征序列，具体包括：

基于以分钟为粒度的特征序列统计基本特征；

基于以分钟为粒度的特征序列根据基本特征确定频域特征；

基于以分钟为粒度的特征序列根据基本特征确定差分特征。

优选地，所述对以小时为粒度的特征序列进行异常检测得到异常序列，具体包括：

将以小时为粒度的特征序列转换为以日为粒度的特征序列；

采用滑动四分位法对以日为粒度的特征序列进行异常检测得到异常特征序列；

使用孤立森林算法提取异常特征序列中的超低频电磁扰动特征得到所述异常序列。

基于上述筛选得到的异常序列，本发明对应提供了一种地震发生相关性确定方法，该方法包括：

基于待监测区域的异常序列确定地震发生相关性；所述异常序列由上述提供的地震前兆特征筛选方法确定得到。

对应于上述提供的地震前兆特征筛选方法，本发明还提供了一种地震前兆特征筛选系统，该系统包括：

监测点设置模块，用于在监测区域中布设监测点；

监测模块，用于对待监测区域中所述监测点处的地震前兆信号进行监测，获得地震前兆原始数据序列；所述地震前兆原始数据序列包括：待监测区域中监测点处不同时间的地震电磁扰动信号；

预处理模块，用于对所述地震前兆原始数据序列进行预处理，获得以小时为粒度的特征序列；

异常检测模块，用于对以小时为粒度的特征序列进行异常检测得到异常序列。

对应于上述提供的地震发生相关性确定方法，本发明还提供了一种地震发生相关性确定系统，该系统包括：

相关性确定模块，用于基于待监测区域的异常序列确定地震发生相关性；所述异常序列由上述提供的地震前兆特征筛选方法确定得到。

此外，本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质用于存储计算机软件程序；所述计算机软件程序用于执行上述提供的地震前兆特征筛选方法，或用于执行上述提供的地震发生相关性确定方法。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明依据实际监测精度的要求在待监测地区设置多个监测点，对监测点处不同时间的地震前兆信号的监测数据和监测时间内发生地震事件的时间进行监测，并对地震事件发生前预设时间段的监测数据进行统计提取其异常特征，并分析异常特征与地震事件的相关性，不需要对被测点附近地表和建筑产生任何损害，以为地震预测模型的构建提供高质量的数据输入，进而能够显著提高地震预测的准确性。并且能通过固定监测台站对待监测地区进行长期监测，具有较好的执行可行性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的地震前兆特征筛选方法的流程图；

图2为本发明实施例提供的对待监测区域设定监测点的示意图；

图3为本发明实施例提供的滑动四分位异常检测结果图；

图4为本发明实施例提供的符合条件的台站的位置示意图；

图5为本发明实施例提供的孤立森林算法异常值检测结果图；

图6为本发明实施例提供的珙县气象局地震相关性分析结果图；

图7为本发明实施例提供的马边气象局地震相关性分析结果图；

图8为本发明实施例提供的茂县测点地震相关性分析结果图；

图9为本发明实施例提供的峨眉山防震减灾局地震相关性分析结果图；

图10为本发明实施例提供的昭通地震台地震相关性分析结果图；

图11为本发明实施例提供的青神防震减灾局地震相关性分析结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1所示，本发明提供的地震前兆特征筛选分析方法的具体实施过程包括：

步骤S1：在监测区域中布设监测点。

首先选择待监测区域，在该区域中设置多个监测点。如图2所示，曲线202圈定的范围为待监测区域，标线203和标线204分别为定位精度P(单位米)的长和宽。依据监测精度要求在待监测区域设定多个监测点；多个监测点包括固定监测点201和移动监测点205，监测点的布设密度取决于对于当地地震的监测精度，精度要求越高，监测点布设的密度就越高。本实施例中，固定监测点201可实现长期稳定监测，基于丰富事件提供丰富的原始数据，以便进行地震前兆特征的筛查。一实施例中，数据相似性判断包括但不限于余弦距离法和/或欧氏距离法。

步骤S2：对待监测区域中监测点处的地震前兆信号进行监测，获得地震前兆原始数据序列；所述地震前兆原始数据序列包括：待监测区域中监测点处不同时间的地震电磁扰动信号。

对待监测区域中监测点处地震前兆信号进行监测，获得待监测区域中监测点处不同时间的地震电磁扰动信号和监测时间内发生地震事件。地震事件为监测点附近100km范围内发生的所有地震事件或台站周边200km以内的Ms3.5+地震事件。监测时间根据实际需求进行设置，优选时间长度至少大于1个月（例如，可以自地震前兆监测点安装之日起起算）。

步骤S3：对所述地震前兆原始数据序列进行预处理，获得以小时为粒度的特征序列；该步骤具体包括：

步骤S31：对所述地震前兆原始数据序列进行滤波处理，获得滤波后的地震前兆原始数据序列。

本实施例中，利用AETA系统监测终端获得地震前兆原始数据序列，地震前兆原始数据序列中的原始数据的频率范围覆盖为0.1Hz~10kHz，本发明主要关注原始数据在低频范围0.1Hz~200Hz中的低于30Hz的超低频分量的前兆特征，故需设计一数字滤波器将高于30Hz的成分滤除。

考虑该滤波步骤希望尽量不改变各频段信号的对滤波结果的相对影响，且考虑到滤波器阶数的滤波效果与AETA终端设备运算力的平衡，最终选用八阶巴特沃斯滤波器。

基于上述处理获取的滤波后的地震前兆原始数据序列实质是一种原始高频数据。

步骤S32：基于滤波后的地震前兆原始数据序列确定以分钟为粒度的特征序列；以分钟为粒度的特征序列包括以分钟为粒度的ULF特征序列和功率谱密度特征序列；该步骤的实施过程可以为：

步骤S321：以设定时间为窗口，在滤波后的地震前兆原始数据序列上滑动，获得ULF特征序列，具体包括：

以设定时间为窗口，在滤波后的地震前兆原始数据序列上截取滤波后的地震前兆原始数据，并基于截取的滤波后的地震前兆原始数据求取平均值，获得ULF特征。

在滤波后的地震前兆原始数据序列上滑动窗口，获得一个以分钟为粒度的ULF特征序列；所述ULF特征序列包括多个ULF特征。

本实施例中，将一分钟设置为设定时间，一秒采样500个点，故一分钟的时间窗口内有30000个点，对这三万个点求出来一个平均值，就是这里所说的平均值”ULF特征”，这个特征仅代表该窗口内的数据。在进行处理时，按照分钟将数据截取成互不重合的段，因此最后得到的序列是分钟粒度的。

步骤S322：使用时频分析方法对滤波后的地震前兆原始数据序列进行时域分析，获得一个以分钟为粒度的功率谱密度特征序列。

本实施例选用短时傅里叶变换进行处理，具体公式为：

其中，W(t-τ)为窗口函数，τ为希望的窗口尺度，X(τ,ω)为短时傅里叶变换。

本实施例将窗口函数视为一阶跃函数，即针对每分钟的数据段将其单独切分出来，并加以快速傅里叶变换，得到该分钟内数据的功率谱密度，并将30Hz以下的部分按照1Hz的频率分辨率提出，作为一中间结果。

到此为止，针对原始数据已得到多条以分钟为粒度的特征值序列，包括ULF特征序列和功率谱密度特征序列。

步骤S33：基于以分钟为粒度的特征序列确定以小时为粒度的特征序列，以达到抑制短时噪声影响的目的。该步骤的实施过程可以包括：

步骤S331：基于以分钟为粒度的特征序列统计基本特征；

具体的，针对ULF特征序列统计其均值MEAN、标准差STD和最大值MAX。针对功率谱密度特征序列提取各频段功率谱密度的均值PSD_MEAN。提取频段可以依据实际需求进行确定，例如，以5Hz为间隔进行提取。

步骤S332：基于以分钟为粒度的特征序列根据基本特征确定频域特征；

本实施例选用小时峰值频率MF和功率谱熵HF进行处理，采用的公式如下：

其中，PSD_f表示1Hz-30Hz各个频率区间的功率谱密度值。

步骤S333：基于以分钟为粒度的特征序列根据基本特征确定差分特征；

本实施例选用小时周期差分和日周期差分F进行处理，在数据分析过程中，发现部分台站会出现日周期的变化和年周期的变化，分别设计小时周期差分特征TD和日周期差分特征HD。小时周期差分特征指该小时某特征与上一小时该特征的差值，日周期差分特征指该小时某特征与上一日该小时该特征的差值。采用的公式如下：

到此为止，得到的是一系列以小时为粒度的特征序列{V1,V2,V3…}。其中，V{MEAN}：{MEAN1.MEAN2,MEAN3,……}，记作V1。V{STD}:{STD1, STD2,STD3,……}记作V2，以此类推。各个序列的特征如下表1所示。

表1

步骤S4：对以小时为粒度的特征序列进行异常检测得到异常序列。该步骤具体包括：

步骤S41：将以小时为粒度的特征序列转换为以日为粒度的特征序列；

具体的，针对以小时为粒度的特征序列{V1,V2,V3…}，确定每一日24小时内的均值类数据的均值，最值（最大值、最小值）取每日的对应最值，将特征序列的粒度转化为日粒度，得到以日为粒度的特征序列{S1,S2,S3…}。

步骤S42：采用滑动四分位法对以日为粒度的特征序列{S1,S2,S3…}进行异常检测得到异常特征序列。

具体的，滑动四分位法基于四分位法对连续序列进行了改进。四分位法计算统计变量的分散情况，将数据从小到大排列，位于25%位置上的数为第一四分位数Q1,位于75%位置上的数为第三四分位数Q3，定义四分位距IQR为Q3-Q1。

在应用滑动四分位法时，指定一个窗口长度L，以MEAN特征为例，对长为L的序列计算其IQR值，并针对末尾位置的MEAN按照以下公式计算异常值Anomaly，作为该滑动区间序列中末尾位置所对应日的异常分数。

其中，Q2为第二分位数，即中位数，up和low为经计算所得的上下阈值，Anomaly则描述了被检测值偏离该上下阈值所标定的范围的程度。

将窗口不断向后滑动，则可以得到以日为粒度的滑动四分位异常值序列{S’1,S’2,S’3…}。

对陇川县防震减灾局2019年4月至2019年7月的数据进行ULF特征滑动四分位异常检测，其检测结果如图3所示。

步骤S43：使用孤立森林算法提取异常特征序列中的超低频电磁扰动特征得到异常序列。

孤立森林算法在2008年由刘飞、周志华等提出，是较为成熟的异常提取方案。在本发明的场景中，孤立森林算法会处理输入的多个特征值，给出针对所有特征值的整体异常情况的异常值。

针对以小时为粒度的特征序列{V1,V2,V3…}，最终得到以小时为粒度的异常序列V’’。以天为单位取滑窗，每个滑窗内是每天的24小时粒度的特征。对每个小时内的多个特征V1,V2,V3使用孤立森林算法进行异常检测，可以得到该小时是否存在异常。用每天的异常小时个数除24，得到当日的异常值。该异常值即为超低频电磁扰动异常D,其粒度是日。本实施例中，对青神防灾减震局应用了孤立森林的算法基于前步骤所述特征{V1,V2,V3…}进行了多特征异常检测，效果如下述Es指数公式所示。

步骤S5：基于异常序列确定地震发生相关性；

本步骤对上述得到的以日为粒度的异常序列确定地震事件评估序列与地震的相关性。被评估的序列可以随意选择，但需要满足日粒度的条件。本实施例中，本发明采用孤立森林算法评估所得以日为粒度的异常序列与地震的相关性。此外，在本实施例中，事先依据监测点周围地震事件的密度以及监测点本身数据的完整性对使用的监测点进行初步筛选，具体准则要求该台站在2017年9月1日至2019年11月10日之间地震日至少 10 个，且在该时段内完整数据的天数应高于700。在该实施例中我们最终选取12个台站，具体信息如下表2所示，其位置如图4所示。

表2

基于以日为粒度的异常序列和已有地震事件进行相关性评估。地震的震级和震中距是表征地震对某区域影响的两个重要的参数，另外，很多震级较高的地震会伴有密集的余震序列，综合以上因素，使用Es参数作为地震日的定义标准，对地震事件进行量化。

其中，M为里氏震级，r为观测点与震中的距离，Es’可看作单个地震对观测点的影响，Es为发生在同一天的地震的叠加影响值。Es的阈值这里取10⁸，Es大于该值的日期被认为是该观测点的地震日。

到此为止，得到了以日为粒度的该观测点地震日的序列，基于以日为粒度的该观测点地震日的序列确定其与观测点地震日之间的相关性。

步骤如下：

(1)获取某台站的异常序列AI(t)（日粒度）;

(2)获取该台站该时间段内符合条件的地震日列表，设有N个地震日;

(3)针对待评估的异常类别，取出其前后15天的异常序列，得到长度为31的异常序列{EVi}；

(4)从异常序列中随机抽取N天，将这N天的异常值加起来得到一个背景值Val；

重复步骤（1）-（4）10000次得到序列{Vali}

计算序列{Vali}的均值μ和标准差σ。

其中，异常分数s采用孤立森林方法计算，计算公式为：

图5以及图6-图11展示了该实施例中几个时序叠加分析的结果。该结果应用于孤立森林算法所得的超低频电磁扰动特征（即以日为粒度的异常序列）。图7中竖直方向的虚线代表地震发生日，水平方向从上到下的三条直线依次是μ+2σ，μ+σ，μ，μ-σ。若AETA的异常值与地震的确存在关联的话，在地震前后叠加异常指数应该会明显高于平时。当叠加异常值高于μ+2σ时，可以认为两者具有较强的相关性。

从结果可以看出珙县气象局、马边气象局、茂县测点、峨眉山防震减灾局、昭通地震台和青神防震减灾局都存在异常指数超过μ+2σ的时期，其中珙县气象局、马边气象局和峨眉山防震减灾局位于地震之前，茂县测点、昭通地震台和青神防震减灾局位于地震之后；这六个台站中除茂县测点外叠加异常指数基本都一直位于随机平均值μ之上。综上所述，可以有较大把握认为所提取的 AETA 超低频电磁扰动异常与地震事件有关。对其他种类特征，可应用该工作流程进行同样的相关性判断。

基于以上实施例，依据实际监测精度的要求在待监测地区设置多个监测点，对监测点处不同时间的地震前兆信号的监测数据和监测时间内发生地震事件的时间进行监测，并对地震事件发生前预设时间段的监测数据进行统计，按照实施例方法提取其异常特征，并使用时序叠加法分析异常与地震事件的相关性。本申请创新的提出基于AETA监测数据的异常提取和相关性分析的方法，不需要对被测点附近地表和建筑产生任何损害，且能通过固定监测台站对待监测地区进行长期监测，具有较好的执行可行性。

监测点设置模块，用于在监测区域中布设监测点；

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims

1.一种地震前兆特征筛选方法，其特征在于，包括：

在监测区域中布设监测点；

对以小时为粒度的特征序列进行异常检测得到异常序列。

2.根据权利要求1所述的地震前兆特征筛选方法，其特征在于，所述对所述地震前兆原始数据序列进行预处理，获得以小时为粒度的特征序列，具体包括：

3.根据权利要求2所述的地震前兆特征筛选方法，其特征在于，所述基于滤波后的地震前兆原始数据序列确定以分钟为粒度的特征序列，具体包括：

4.根据权利要求2所述的地震前兆特征筛选方法，其特征在于，所述基于以分钟为粒度的特征序列确定以小时为粒度的特征序列，具体包括：

基于以分钟为粒度的特征序列统计基本特征；

基于以分钟为粒度的特征序列根据基本特征确定频域特征；

基于以分钟为粒度的特征序列根据基本特征确定差分特征。

5.根据权利要求1所述的地震前兆特征筛选方法，其特征在于，所述对以小时为粒度的特征序列进行异常检测得到异常序列，具体包括：

将以小时为粒度的特征序列转换为以日为粒度的特征序列；

6.一种地震前兆特征筛选系统，其特征在于，包括：

监测点设置模块，用于在监测区域中布设监测点；

7.一种地震发生相关性确定方法，其特征在于，包括：

基于待监测区域的异常序列确定地震发生相关性；所述异常序列由权利要求1-5任意一项所述的地震前兆特征筛选方法确定得到。

8.一种地震发生相关性确定系统，其特征在于，包括：

相关性确定模块，用于基于待监测区域的异常序列确定地震发生相关性；所述异常序列由权利要求1-5任意一项所述的地震前兆特征筛选方法确定得到。

9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质用于存储计算机软件程序；所述计算机软件程序用于执行如权利要求1-5任意一项所述的地震前兆特征筛选方法，或用于执行如权利要求7所述的地震发生相关性确定方法。