基于支持向量机的地震事件类型自动识别及应用

蔡杏辉，廖诗荣，张燕明，陈惠芳，林彬华

（福建省地震局，福州 350003）

0 引言

地震事件类型分类是地震监测的必要工作，信息的迅速可靠对地方政府和公众至关重要。目前我国非天然地震的判别与分析主要靠人工进行，响应速度慢。虽然非天然地震事件直接破坏性普遍不大，但仍有一定数量的事件对社会产生重大影响，如2015年8月天津港爆炸，分别记录到ML2.3级与ML2.9级两次非天然地震，2019年3月江苏响水爆炸，记录到ML2.2级非天然地震等，这些突发事件对地震部门快速应急响应提出了新要求，非天然地震进行准实时处理，自动判定事件类型，建立规范化、标准化、自动化处理流程，可为政府的快速应急响应提供技术支撑。此外，由于事件波形复杂性，不能完全剔除非天然地震事件，会影响地震目录产出质量，导致天然地震活动性估计不准确，对疑似事件准确定性可提高地震编目产出质量和工作效率。

目前，已有许多学者对天然及非天然地震的识别进行研究，并取得相关研究成果。如BP神经网络、卷积神经网络、支持向量机、矩阵决策算法等方法被运用于天然地震与非天然地震自动识别分类[1-6]。支持向量机（SVM）是人工智能算法的一个重要分支，是一类按监督学习方式对数据进行二元分类的广义线性分类器，在人像识别、文本分类、模式识别等许多领域中得到应用[7]。已有研究者将支持向量机技术应用到地震事件类型的分类识别研究。国内，研究者[2，3，5]利用波形小波特征，采用支持向量机对地震和爆破事件进行识别。国外，Kortstrom等（2016）[7]提出了一种基于支持向量机的事件全自动分类方法，该方法通过SVM学习天然地震和人工地震在P、P波尾波、S、S波尾波等地震记录不同部分能量分布的差异来进行事件类型判别。Rabin（2016）[8]应用机器学习算法中的扩散图方法，将地震图转换为标准化的超声波图，进行天然地震事件与核爆事件的自动识别。Saad等（2019）[9]采用小波滤波器组联合支持向量机提出了天然地震事件与采石场人工爆破事件的自动分类算法。研究结果表明，支持向量机用于事件分类具有很好的效果。

应用于事件自动识别的算法主要有神经网络和支持向量机，支持向量机在解决小样本、非线性及高维模式识别中表现出许多特有的优势，并能够推广应用到函数拟合等其他机器学习问题中[10]。BP神经网络是包含多个隐含层的网络，具备处理线性不可分问题的能力，但BP神经网络也存在主要缺陷学习速度慢及容易陷入局部极小值；算法特点上支持向量机在样本一定的情况下运算量比BP要小，判识结果产出快；CNN卷积神经是深度学习的代表算法之一，其无需人为提取特定特征信息，但运算量大：非天然地震塌陷、矿震、滑坡为小样本事件，福建台网对塌陷、矿震等特殊地震动鲜有记录到；基于以上原因，因此本研究采用支持向量机作为自动识别分类器。本研究，在前人基础上充分利用福建台网记录的非天然地震、地震数字波形资料，采用波形能量分布、P/S振幅比、小波分析等多种方法，开展基于支持向量机事件类型自动识别研究及软件模块研发，以期提高地震事件类型识别的准确率和应用区域的广泛性。研发可在地震台网实时运行的软件模块，提高人工编目的效率和质量，实现地震事件类型产出自动化、智能化处理，为政府的快速应急响应提供技术支撑。

1 识别算法及基本原理

1.1 支持向量机原理

支持向量机方法的基本思想是：定义最优线性超平面，并把寻找最优线性超平面的算法归结为求解一个凸规划问题。进而基于核函数展开定理，通过非线性映射φ，把样本空间映射到一个高维乃至于无穷维的特征空间，使在特征空间中可以应用线性学习机的方法解决样本空间中的高度非线性分类和回归等问题[11]。简单地说就是升维和线性化。常见的核函数有多项式核、径向基函数核、拉普拉斯、Sigmoid核。基于RBF核可将样本空间映射至无限维空间的考虑，采用径向基函数核作为支持向量机核函数，径向基函数核有两个参数C和gamma。经过尝试多个C和gamma配对组合，采用C为100，gamma为0.49效果较好。对于测试样本，将提取的特征向量联合支持向量机构建自动识别算法，依其输出值a判定该事件的类型，判定标准为：当a=1，判定为天然地震，当a=0，判定为人工爆破（非天然地震事件）。

本文所应用的支持向量机算法（v-SVC）[12]：

（2）选取适当的核函数：K（xi，yj）和参数υ，构造并求解最优化问题。

得最优解a*=（a*1，…，a*m）T。

（4）构造决策函数

1.2 分类识别器设计

自动识别系统关键算法有两个模块：一为特征提取模块，从数据中提取出识别所需特征；其二为识别判定模块，在特征空间中用模式识别方法把识别对象归为某一类别。在对事件类型的分类识别中，如何提取有效的识别特征是识别的关键。P/S振幅比是研究最为深入的一个判别量，在地震与爆破识别中效果较好；小波分析算法被广泛应用于识别研究；Kortstrom（2016）[13]方法是将地震信号通过数据处理，转换成“数值谱图”，在事件分类识别上取得了较好的效果。本文采用P/S振幅比、小波分析、波形能量分布特征等多种特征组合联合支持向量机（SVM）进行事件类型判定。对多种特征组合进行测试，得出识别效果较好的分类识别器。图1为分类识别器（SVM）设计示意图。

图1 分类识别器设计示意图Fig.1 Schematic diagram of classification recognizer design

1.3 小波分析

一维离散小波变换（DWT）能量比特征：

若S为原始信号，其长度为J，信号采样点序号为j，Si为信号S分解后的第i个小波系数，其长度为K，k为其样点序号，则小波系数的能量比（Ewt）按（1）式定义[3]。

小波包变换（WPT）小波包对数能量熵特征：

若S为原始信号，对它进行n层小波包分解后，得到第n层的小波包系数总共为N个。Si为信号S分解后的第i个小波的系数，其长度为J，小波系数的结点序号为j，则从第i个小波的系数中提取出小波包对数能量熵（Elg）特征按下式定义[3]。

1.4 波形能量分布

采用20个窄带通道过滤器过滤地震记录（滤波器为零相位二阶巴特沃斯滤波器，频带1~41 Hz），并将事件波形分为四个阶段的窗口：P和P coda，S和S coda。然后，计算了每个滤波器通道和相位窗口的短期平均（STA）值，共80个判别参数[7]。

短期平均（STA）值计算公式：

其中N为样本中STA窗口的长度，yi为第i个样本经过过滤的时间序列y。

2 数据资料及特征组合方式

2.1 数据资料

为使研究成果能应用于台网日常监测与编目工作，需使学习样本泛化能力强。采用福建台网2016年1月至2019年6月所有触发事件，不进行震中区域、台站、震级的划分，事件数6570个，事件震级范围ML1.1～4.0，包括极低震级单台事件，每个事件各提取出一个单台波形记录（Δ＜100 km），训练集和测试集按约50%：50%的比例，随机将6570个事件中的3253个事件作为训练集，其余3317个事件作为测试集。训练集907个为天然地震记录，另2346个为人工爆破记录；测试集2054个为天然地震记录，另1263个为人工爆破记录。

2.2 特征提取及组合方式

对信号特征提取并截取有效波形，首先要确定P到时，对离线数据识别训练和测试，直接从Jopens震相数据表中获取台站震相P到时数据；处理波形窗长选择，选取P波到时前0.5 s至P波到时后19.5 s，共20 s的单台Z分向和三分向事件波形（需包含完整的P与S波列）进行处理。

P/S振幅比特征的提取：信号的S波和P波最大振幅比是天然地震和人工爆破波形信号中的一个重要特征，Colin等[14]研究表明6～8 Hz和8～10 Hz的滤波频带对爆破的识别效果较好，因此采用一阶Butterworth对原始波形进行滤波，共选取3个滤波频带，分别为6.0～8.0、8.0～10.0、6.0～10.0，分别量取P与S波列的最大值，获取3个P/S振幅比特征值。小波特征的提取：对信号进行4层小波及小波包分解。对小波分解得到小波系数利用公式（1）提取出小波能量比（Ewt）特征组成5维特征向量；对第4层的16个波包系数利用公式（2）提取小波包对数能量熵波形特征（Elg）分别组成16维特征向量。波形能量分布特征的提取：事件波形四个阶段的窗口，通过20个窄带通道过滤器滤波，利用公式（3）计算短期平均（STA）值，获得80个特征值。

将提取的特征进行组合，在特征组合方式选取上，采用以下3种组合方式：①P/S振幅比+小波包对数能量熵；②P/S振幅比+小波能量比+小波包对数能量熵；③P/S振幅比+小波包对数能量熵+波形能量分布（P，P code，S，S code）。分别组成19维、24维、99维特征向量，利用支持向量机检验各特征组合的分类能力。

3 识别能力检验及分析

3.1 各特征组合训练及测试结果

将测试集记录用上述方法进行特征提取并组合，用对应特征组合训练得到的支持向量机进行地震事件类型分类，检验各特征组合的分类能力。测试数据输入采用单台Z分向和三分向数据作为输入，首先采用单台Z分向数据测试各个特征组合分类效果，获得分类效果最好的特征组合，其次对该特征组合采用单台三分向数据测试分类效果，对比不同数据的输入方式对分类效果的影响。训练及测试结果见表1。

表1 各特征组合对训练集及测试集的分类效果Table 1 Classification effect of each feature combination on training set and test set

由表1结果可知，3种特征组合训练集和测试集准确率达到80%以上，说明选用的特征指标可对事件进行有效分类，有效特征进行组合有助于识别率的提高；其中P/S振幅比+小波包对数能量熵+波形能量分布特征组合的识别率最高，自验证准确率达到96.8%，分类准确率达到91.3%；对该特征组合采用单台三分向数据输入，自验证准确率达到98.3%，分类准确率达到94.5%，分类识别率比单台单通道输入提高了3.2%。

3.2 本文最优算法性能评估及事件实列

表2给出了本文最优算法性能评估指标结果，结果表明，灵敏度和特异度比较均衡，精确率达到96.0%，说明算法比较好，对地震和爆破多数都能有效识别；灵敏度稍高于特异度，说明对地震识别效果好于爆破；该算法对地震分类的判别准确率为95.5%，对爆破分类的判别准确率为93.4%。算法测试3308个事件，共消耗的时间8270 s，平均每个事件全流程产出结果需耗时2.5 s；测试环境：使用英特尔处理器核心E5-26090@2.40 GHz，内存8 GB，windows 10，64位操作系统。

表2 本文最优算法性能评估指标Table 2 The performance evaluation index of optimal algorithm given in this paper

图2、3分别是福建仙游爆破正确和错误识别事件仙游西苑台波形，表3是其波形特征计算参数，每个事件有99个计算参数。采用同一类型事件同一台站对正确和错误识别事件波形特征计算参数进行对比更具有可比性，对比发现特征序列的部分特征值存在较大的差别；对两者波形进行对比，错误识别事件波形背景噪声较大，可能事件波形背景噪声较大对特征计算参数结果有影响，导致事件错误识别。

图2 福建仙游爆破正确识别事件波形（2020-05-07）Fig.2 The waveform of the correct identification of the blasting event in Xianyou，Fujian

表3 最优算法正确和错误事件波形特征计算参数Table 3 The optimal algorithm calculates the parameters of the characteristic of the correct and wrong event waveform

图3 福建仙游爆破错误识别事件波形（2019-05-24）Fig.3 The waveform of the misidentification of the blasting event in Xianyou，Fujian

图4、图5为本文最优算法正确和错误识别事件震中分布图，正确识别的事件中共有1962个地震和1180个爆破事件；错误识别的事件中共有92个地震和83个爆破事件；爆破错识率高于地震错识率，分析错误识别的事件波形，发现个别事件波形背景噪声大和波形零漂的情况；错误识别的爆破事件分布较为分散，错误识别的地震事件则主要集中在台湾海峡区域，占错误识别的地震事件总数51%，本文研究训练集未将台湾海峡地震作为训练样本，说明学习样本覆盖，对识别效果存在影响；通过图4、图5比较，错误识别的爆破和地震均有与正确识别的爆破和地震位置重叠，反映了算法模型具有一定程度的不一致性。

图4 正确识别事件震中分布图Fig.4 Epicentral distribution of correct identification events

图5 错误识别事件震中分布图Fig.5 Epicentral distribution of misidentified events

4 自动识别模块研制及应用

选择本文最优算法研发自动识别模块，自动识别模块的研发采用Python语言编程实现，采用obspy库进行地震数据处理和分析。模块的实现还需要应用到关键的第三方库有PyMysql、redis、stomp。自动识别模块有日常人工编目应用和与自动编目系统对接的在线应用。福建台网区域自动地震编目系统已投入测试运行，初步具备多台触发地震事件的自动编目能力，自动识别模块作为自动地震编目系统子功能模块，需研发与自动编目系统的接口，实现本台网事件类型自动识别日常化应用。

4．1 自动识别模块人工编目应用

4.1.1 实现方法

事件自动导入，将前一天Jopens数据库触发波形事件（00~24 h）自动导入指定文件夹；自动识别模块定时自启动，读入指定文件夹中事件波形，从震相数据表中获取台站震相P到时，截取有效波形，经特征值提取步骤，进行事件类型判断，并与数据库人工分析结果对比，产出类型不一致结果报表文件。软件流程自动化完成，无需人为操作。图6为自动识别模块人工编目应用技术路线图。

图6 自动识别模块人工编目应用技术路线图Fig.6 Technical roadmap for manual cataloging application of automatic identification module

4.1.2 人工编目应用

事件类型自动识别模块已应用于人工编目，软件自动对每天产出触发事件进行检测，产出报表文件供编目分析人员参考，极大减小人工识别事件类型的误判；软件运行以来，检测出多个人工编目误判事件及对可疑事件性质进行判识。2019年8月1日至2020年2月29日，福建台网共实时触发1531个事件，事件震级范围ML1.0~3.3，其中天然地震911个、爆破事件620个。从检测结果来看，事件类型识别准确率为93.9%，实际应用识别准确率接近测试集的测试准确率，共检测出8个事件类型人工识别错误，错误提交数据库；其中天然地震的识别准确率为95.7%，错误识别的地震事件中有20个为台湾海峡地震，占地震错识总数51.3%，有3个为人工难以判识别事件，事件波形与爆破极为相似；2019年9月26日至27日共发生仙游小震群共152个地震事件，对该仙游小震群识别准确率达到100%；爆破的识别准确率为92.1%，错误识别的爆破事件中有2个信噪比低；地震识别效果优于爆破，与测试集的测试结果一致。表4为福建台网2019年8月至2020年2月实时触发事件识别情况统计表。

表4 2019-08-01—2020-02-29福建台网实时触发事件识别情况统计表Table 4 Statistical table of real-time triggered events recognition in Fujian Network from August 1st，2019 to February 29th，2020

4.2 自动识别模块自动编目在线应用

4.2.1 实现方法

AMQ消息中间件是Apache出品，最流行的，能力强劲的开源消息总线。AMQ的消息传递模式有二种，一种为点对点模式，另一种为发布/订阅模式[10]；自动识别模块与自动编目接口开发应用第二种模式。redis是完全开源免费的，遵守BSD协议，是一个高性能的key-value数据库[15]。采用AMQ消息中间件为信息中介（需下载安装AMQ并启动AMQ服务），redis波形共享内存中获取事件波形，研发了自动识别模块与自动编目系统接口，将自动识别模块嵌入自动编目系统。具体流程：从消息中间件接收事件参数信息，获取目标地震记录的开始与结束信息，从redis波形共享内存中获取事件波形，经过数据预处理、特征值提取等步骤，进行事件类型判断；产出的事件类型信息，通过消息中间与其他软件模块共享。图7为自动识别模块与自动编目接口技术路线图。

图7 自动识别模块与自动编目接口技术路线图Fig.7 Technical roadmap of interface between automatic identification module and automatic cataloging

4.2.2 在线应用

目前自动识别模块已具备在线运行的能力，初步实现事件类型的实时全自动判别。2019年11月至2020年3月已在线试验运行5个月，其间检测出多个人工编目漏分析地震事件，提高目录完整性。对自动识别在线系统7×24运行稳定性，进行不间断测试，系统未出现死机、进程中断现象，测试运行性能稳定良好，功能和性能达到预期设计，验证了自动识别模块实用性和稳定性。

5 讨论和结论

（1）本文通过对多种特征组合联合支持向量机进行大批量事件类型测试分析，结果显示特征指标可对事件进行有效分类；据此得出识别效果较好的事件类型识别算法，事件类型识别准确率达到94.5%，灵敏度和特异度比较均衡，灵敏度和特异度均达到93%以上，对多数地震和爆破都能有效识别，极低震级单台事件也能较好的分类，判识结果产出快，单个事件全流程识别仅需耗时2.5 s。

（2）算法数据输入采用单台Z分向和三分向数据作为输入；采用同样特征组合，三分向数据识别效果优于单分向；可能的解释是更多输入数据对于支持向量机的输入参数约束更好。此外，事件类型自动识别研究，训练集中含有识别错误的事件也可能影响支持向量机分类性能；通常情况下，区域台网地震目录很少是100%正确；在对训练集自验证，最初训练支持向量机已能准确识别出人工编目错误识别事件，其中地震训练集中有3个错误识别事件，人工爆破中有8个错误识别事件；在最终测试前，对此数据进行剔除。

（3）采用本文研究最优算法联合支持向量机研制事件类型自动识别模块，日常人工编目实际应用效果显示，识别准确率达到93.9%，可检测出人工识别错误事件，减少日常编目工作量和错误率。另外，实现与自动编目系统对接的在线自动识别系统接口，在线自动识别系统测试运行性能稳定良好，检测出多个人工漏分析地震事件，自动识别模块具有实用性和稳定性。

在后期工作中，增加其他区域事件进行完善和优化学习样本，并改进算法提高系统的泛化能力，自动识别模块将会更具有普适性，实现事件类型识别的多样性，塌陷、矿震、滑破等，可期应用于其他区域事件类型自动判别。