蔡杏輝，廖詩榮，張燕明，陳惠芳，林彬華

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

0 引言

地震事件類型分類是地震監(jiān)測的必要工作，信息的迅速可靠對地方政府和公眾至關(guān)重要。目前我國非天然地震的判別與分析主要靠人工進行，響應速度慢。雖然非天然地震事件直接破壞性普遍不大，但仍有一定數(shù)量的事件對社會產(chǎn)生重大影響，如2015年8月天津港爆炸，分別記錄到ML2.3級與ML2.9級兩次非天然地震，2019年3月江蘇響水爆炸，記錄到ML2.2級非天然地震等，這些突發(fā)事件對地震部門快速應急響應提出了新要求，非天然地震進行準實時處理，自動判定事件類型，建立規(guī)范化、標準化、自動化處理流程，可為政府的快速應急響應提供技術(shù)支撐。此外，由于事件波形復雜性，不能完全剔除非天然地震事件，會影響地震目錄產(chǎn)出質(zhì)量，導致天然地震活動性估計不準確，對疑似事件準確定性可提高地震編目產(chǎn)出質(zhì)量和工作效率。

目前，已有許多學者對天然及非天然地震的識別進行研究，并取得相關(guān)研究成果。如BP神經(jīng)網(wǎng)絡(luò)、卷積神經(jīng)網(wǎng)絡(luò)、支持向量機、矩陣決策算法等方法被運用于天然地震與非天然地震自動識別分類[1-6]。支持向量機（SVM）是人工智能算法的一個重要分支，是一類按監(jiān)督學習方式對數(shù)據(jù)進行二元分類的廣義線性分類器，在人像識別、文本分類、模式識別等許多領(lǐng)域中得到應用[7]。已有研究者將支持向量機技術(shù)應用到地震事件類型的分類識別研究。國內(nèi)，研究者[2，3，5]利用波形小波特征，采用支持向量機對地震和爆破事件進行識別。國外，Kortstrom等（2016）[7]提出了一種基于支持向量機的事件全自動分類方法，該方法通過SVM學習天然地震和人工地震在P、P波尾波、S、S波尾波等地震記錄不同部分能量分布的差異來進行事件類型判別。Rabin（2016）[8]應用機器學習算法中的擴散圖方法，將地震圖轉(zhuǎn)換為標準化的超聲波圖，進行天然地震事件與核爆事件的自動識別。Saad等（2019）[9]采用小波濾波器組聯(lián)合支持向量機提出了天然地震事件與采石場人工爆破事件的自動分類算法。研究結(jié)果表明，支持向量機用于事件分類具有很好的效果。

應用于事件自動識別的算法主要有神經(jīng)網(wǎng)絡(luò)和支持向量機，支持向量機在解決小樣本、非線性及高維模式識別中表現(xiàn)出許多特有的優(yōu)勢，并能夠推廣應用到函數(shù)擬合等其他機器學習問題中[10]。BP神經(jīng)網(wǎng)絡(luò)是包含多個隱含層的網(wǎng)絡(luò)，具備處理線性不可分問題的能力，但BP神經(jīng)網(wǎng)絡(luò)也存在主要缺陷學習速度慢及容易陷入局部極小值；算法特點上支持向量機在樣本一定的情況下運算量比BP要小，判識結(jié)果產(chǎn)出快；CNN卷積神經(jīng)是深度學習的代表算法之一，其無需人為提取特定特征信息，但運算量大：非天然地震塌陷、礦震、滑坡為小樣本事件，福建臺網(wǎng)對塌陷、礦震等特殊地震動鮮有記錄到；基于以上原因，因此本研究采用支持向量機作為自動識別分類器。本研究，在前人基礎(chǔ)上充分利用福建臺網(wǎng)記錄的非天然地震、地震數(shù)字波形資料，采用波形能量分布、P/S振幅比、小波分析等多種方法，開展基于支持向量機事件類型自動識別研究及軟件模塊研發(fā)，以期提高地震事件類型識別的準確率和應用區(qū)域的廣泛性。研發(fā)可在地震臺網(wǎng)實時運行的軟件模塊，提高人工編目的效率和質(zhì)量，實現(xiàn)地震事件類型產(chǎn)出自動化、智能化處理，為政府的快速應急響應提供技術(shù)支撐。

1 識別算法及基本原理

1.1 支持向量機原理

支持向量機方法的基本思想是：定義最優(yōu)線性超平面，并把尋找最優(yōu)線性超平面的算法歸結(jié)為求解一個凸規(guī)劃問題。進而基于核函數(shù)展開定理，通過非線性映射φ，把樣本空間映射到一個高維乃至于無窮維的特征空間，使在特征空間中可以應用線性學習機的方法解決樣本空間中的高度非線性分類和回歸等問題[11]。簡單地說就是升維和線性化。常見的核函數(shù)有多項式核、徑向基函數(shù)核、拉普拉斯、Sigmoid核?；赗BF核可將樣本空間映射至無限維空間的考慮，采用徑向基函數(shù)核作為支持向量機核函數(shù)，徑向基函數(shù)核有兩個參數(shù)C和gamma。經(jīng)過嘗試多個C和gamma配對組合，采用C為100，gamma為0.49效果較好。對于測試樣本，將提取的特征向量聯(lián)合支持向量機構(gòu)建自動識別算法，依其輸出值a判定該事件的類型，判定標準為：當a=1，判定為天然地震，當a=0，判定為人工爆破（非天然地震事件）。

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

（2）選取適當?shù)暮撕瘮?shù)：K（xi，yj）和參數(shù)υ，構(gòu)造并求解最優(yōu)化問題。

得最優(yōu)解a*=（a*1，…，a*m）T。

（4）構(gòu)造決策函數(shù)

1.2 分類識別器設(shè)計

自動識別系統(tǒng)關(guān)鍵算法有兩個模塊：一為特征提取模塊，從數(shù)據(jù)中提取出識別所需特征；其二為識別判定模塊，在特征空間中用模式識別方法把識別對象歸為某一類別。在對事件類型的分類識別中，如何提取有效的識別特征是識別的關(guān)鍵。P/S振幅比是研究最為深入的一個判別量，在地震與爆破識別中效果較好；小波分析算法被廣泛應用于識別研究；Kortstrom（2016）[13]方法是將地震信號通過數(shù)據(jù)處理，轉(zhuǎn)換成“數(shù)值譜圖”，在事件分類識別上取得了較好的效果。本文采用P/S振幅比、小波分析、波形能量分布特征等多種特征組合聯(lián)合支持向量機（SVM）進行事件類型判定。對多種特征組合進行測試，得出識別效果較好的分類識別器。圖1為分類識別器（SVM）設(shè)計示意圖。

圖1 分類識別器設(shè)計示意圖Fig.1 Schematic diagram of classification recognizer design

1.3 小波分析

一維離散小波變換（DWT）能量比特征：

若S為原始信號，其長度為J，信號采樣點序號為j，Si為信號S分解后的第i個小波系數(shù)，其長度為K，k為其樣點序號，則小波系數(shù)的能量比（Ewt）按（1）式定義[3]。

小波包變換（WPT）小波包對數(shù)能量熵特征：

若S為原始信號，對它進行n層小波包分解后，得到第n層的小波包系數(shù)總共為N個。Si為信號S分解后的第i個小波的系數(shù)，其長度為J，小波系數(shù)的結(jié)點序號為j，則從第i個小波的系數(shù)中提取出小波包對數(shù)能量熵（Elg）特征按下式定義[3]。

1.4 波形能量分布

采用20個窄帶通道過濾器過濾地震記錄（濾波器為零相位二階巴特沃斯濾波器，頻帶1~41 Hz），并將事件波形分為四個階段的窗口：P和P coda，S和S coda。然后，計算了每個濾波器通道和相位窗口的短期平均（STA）值，共80個判別參數(shù)[7]。

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

其中N為樣本中STA窗口的長度，yi為第i個樣本經(jīng)過過濾的時間序列y。

2 數(shù)據(jù)資料及特征組合方式

2.1 數(shù)據(jù)資料

為使研究成果能應用于臺網(wǎng)日常監(jiān)測與編目工作，需使學習樣本泛化能力強。采用福建臺網(wǎng)2016年1月至2019年6月所有觸發(fā)事件，不進行震中區(qū)域、臺站、震級的劃分，事件數(shù)6570個，事件震級范圍ML1.1～4.0，包括極低震級單臺事件，每個事件各提取出一個單臺波形記錄（Δ＜100 km），訓練集和測試集按約50%：50%的比例，隨機將6570個事件中的3253個事件作為訓練集，其余3317個事件作為測試集。訓練集907個為天然地震記錄，另2346個為人工爆破記錄；測試集2054個為天然地震記錄，另1263個為人工爆破記錄。

2.2 特征提取及組合方式

對信號特征提取并截取有效波形，首先要確定P到時，對離線數(shù)據(jù)識別訓練和測試，直接從Jopens震相數(shù)據(jù)表中獲取臺站震相P到時數(shù)據(jù)；處理波形窗長選擇，選取P波到時前0.5 s至P波到時后19.5 s，共20 s的單臺Z分向和三分向事件波形（需包含完整的P與S波列）進行處理。

P/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振幅比特征值。小波特征的提?。簩π盘栠M行4層小波及小波包分解。對小波分解得到小波系數(shù)利用公式（1）提取出小波能量比（Ewt）特征組成5維特征向量；對第4層的16個波包系數(shù)利用公式（2）提取小波包對數(shù)能量熵波形特征（Elg）分別組成16維特征向量。波形能量分布特征的提取：事件波形四個階段的窗口，通過20個窄帶通道過濾器濾波，利用公式（3）計算短期平均（STA）值，獲得80個特征值。

將提取的特征進行組合，在特征組合方式選取上，采用以下3種組合方式：①P/S振幅比+小波包對數(shù)能量熵；②P/S振幅比+小波能量比+小波包對數(shù)能量熵；③P/S振幅比+小波包對數(shù)能量熵+波形能量分布（P，P code，S，S code）。分別組成19維、24維、99維特征向量，利用支持向量機檢驗各特征組合的分類能力。

3 識別能力檢驗及分析

3.1 各特征組合訓練及測試結(jié)果

將測試集記錄用上述方法進行特征提取并組合，用對應特征組合訓練得到的支持向量機進行地震事件類型分類，檢驗各特征組合的分類能力。測試數(shù)據(jù)輸入采用單臺Z分向和三分向數(shù)據(jù)作為輸入，首先采用單臺Z分向數(shù)據(jù)測試各個特征組合分類效果，獲得分類效果最好的特征組合，其次對該特征組合采用單臺三分向數(shù)據(jù)測試分類效果，對比不同數(shù)據(jù)的輸入方式對分類效果的影響。訓練及測試結(jié)果見表1。

表1 各特征組合對訓練集及測試集的分類效果Table 1 Classification effect of each feature combination on training set and test set

由表1結(jié)果可知，3種特征組合訓練集和測試集準確率達到80%以上，說明選用的特征指標可對事件進行有效分類，有效特征進行組合有助于識別率的提高；其中P/S振幅比+小波包對數(shù)能量熵+波形能量分布特征組合的識別率最高，自驗證準確率達到96.8%，分類準確率達到91.3%；對該特征組合采用單臺三分向數(shù)據(jù)輸入，自驗證準確率達到98.3%，分類準確率達到94.5%，分類識別率比單臺單通道輸入提高了3.2%。

3.2 本文最優(yōu)算法性能評估及事件實列

表2給出了本文最優(yōu)算法性能評估指標結(jié)果，結(jié)果表明，靈敏度和特異度比較均衡，精確率達到96.0%，說明算法比較好，對地震和爆破多數(shù)都能有效識別；靈敏度稍高于特異度，說明對地震識別效果好于爆破；該算法對地震分類的判別準確率為95.5%，對爆破分類的判別準確率為93.4%。算法測試3308個事件，共消耗的時間8270 s，平均每個事件全流程產(chǎn)出結(jié)果需耗時2.5 s；測試環(huán)境：使用英特爾處理器核心E5-26090@2.40 GHz，內(nèi)存8 GB，windows 10，64位操作系統(tǒng)。

表2 本文最優(yōu)算法性能評估指標Table 2 The performance evaluation index of optimal algorithm given in this paper

圖2、3分別是福建仙游爆破正確和錯誤識別事件仙游西苑臺波形，表3是其波形特征計算參數(shù)，每個事件有99個計算參數(shù)。采用同一類型事件同一臺站對正確和錯誤識別事件波形特征計算參數(shù)進行對比更具有可比性，對比發(fā)現(xiàn)特征序列的部分特征值存在較大的差別；對兩者波形進行對比，錯誤識別事件波形背景噪聲較大，可能事件波形背景噪聲較大對特征計算參數(shù)結(jié)果有影響，導致事件錯誤識別。

圖2 福建仙游爆破正確識別事件波形（2020-05-07）Fig.2 The waveform of the correct identification of the blasting event in Xianyou，F(xiàn)ujian

表3 最優(yōu)算法正確和錯誤事件波形特征計算參數(shù)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，F(xiàn)ujian

圖4、圖5為本文最優(yōu)算法正確和錯誤識別事件震中分布圖，正確識別的事件中共有1962個地震和1180個爆破事件；錯誤識別的事件中共有92個地震和83個爆破事件；爆破錯識率高于地震錯識率，分析錯誤識別的事件波形，發(fā)現(xiàn)個別事件波形背景噪聲大和波形零漂的情況；錯誤識別的爆破事件分布較為分散，錯誤識別的地震事件則主要集中在臺灣海峽區(qū)域，占錯誤識別的地震事件總數(shù)51%，本文研究訓練集未將臺灣海峽地震作為訓練樣本，說明學習樣本覆蓋，對識別效果存在影響；通過圖4、圖5比較，錯誤識別的爆破和地震均有與正確識別的爆破和地震位置重疊，反映了算法模型具有一定程度的不一致性。

圖4 正確識別事件震中分布圖Fig.4 Epicentral distribution of correct identification events

圖5 錯誤識別事件震中分布圖Fig.5 Epicentral distribution of misidentified events

4 自動識別模塊研制及應用

選擇本文最優(yōu)算法研發(fā)自動識別模塊，自動識別模塊的研發(fā)采用Python語言編程實現(xiàn)，采用obspy庫進行地震數(shù)據(jù)處理和分析。模塊的實現(xiàn)還需要應用到關(guān)鍵的第三方庫有PyMysql、redis、stomp。自動識別模塊有日常人工編目應用和與自動編目系統(tǒng)對接的在線應用。福建臺網(wǎng)區(qū)域自動地震編目系統(tǒng)已投入測試運行，初步具備多臺觸發(fā)地震事件的自動編目能力，自動識別模塊作為自動地震編目系統(tǒng)子功能模塊，需研發(fā)與自動編目系統(tǒng)的接口，實現(xiàn)本臺網(wǎng)事件類型自動識別日?；瘧?。

4．1 自動識別模塊人工編目應用

4.1.1 實現(xiàn)方法

事件自動導入，將前一天Jopens數(shù)據(jù)庫觸發(fā)波形事件（00~24 h）自動導入指定文件夾；自動識別模塊定時自啟動，讀入指定文件夾中事件波形，從震相數(shù)據(jù)表中獲取臺站震相P到時，截取有效波形，經(jīng)特征值提取步驟，進行事件類型判斷，并與數(shù)據(jù)庫人工分析結(jié)果對比，產(chǎn)出類型不一致結(jié)果報表文件。軟件流程自動化完成，無需人為操作。圖6為自動識別模塊人工編目應用技術(shù)路線圖。

圖6 自動識別模塊人工編目應用技術(shù)路線圖Fig.6 Technical roadmap for manual cataloging application of automatic identification module

4.1.2 人工編目應用

事件類型自動識別模塊已應用于人工編目，軟件自動對每天產(chǎn)出觸發(fā)事件進行檢測，產(chǎn)出報表文件供編目分析人員參考，極大減小人工識別事件類型的誤判；軟件運行以來，檢測出多個人工編目誤判事件及對可疑事件性質(zhì)進行判識。2019年8月1日至2020年2月29日，福建臺網(wǎng)共實時觸發(fā)1531個事件，事件震級范圍ML1.0~3.3，其中天然地震911個、爆破事件620個。從檢測結(jié)果來看，事件類型識別準確率為93.9%，實際應用識別準確率接近測試集的測試準確率，共檢測出8個事件類型人工識別錯誤，錯誤提交數(shù)據(jù)庫；其中天然地震的識別準確率為95.7%，錯誤識別的地震事件中有20個為臺灣海峽地震，占地震錯識總數(shù)51.3%，有3個為人工難以判識別事件，事件波形與爆破極為相似；2019年9月26日至27日共發(fā)生仙游小震群共152個地震事件，對該仙游小震群識別準確率達到100%；爆破的識別準確率為92.1%，錯誤識別的爆破事件中有2個信噪比低；地震識別效果優(yōu)于爆破，與測試集的測試結(jié)果一致。表4為福建臺網(wǎng)2019年8月至2020年2月實時觸發(fā)事件識別情況統(tǒng)計表。

表4 2019-08-01—2020-02-29福建臺網(wǎng)實時觸發(fā)事件識別情況統(tǒng)計表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 實現(xiàn)方法

AMQ消息中間件是Apache出品，最流行的，能力強勁的開源消息總線。AMQ的消息傳遞模式有二種，一種為點對點模式，另一種為發(fā)布/訂閱模式[10]；自動識別模塊與自動編目接口開發(fā)應用第二種模式。redis是完全開源免費的，遵守BSD協(xié)議，是一個高性能的key-value數(shù)據(jù)庫[15]。采用AMQ消息中間件為信息中介（需下載安裝AMQ并啟動AMQ服務），redis波形共享內(nèi)存中獲取事件波形，研發(fā)了自動識別模塊與自動編目系統(tǒng)接口，將自動識別模塊嵌入自動編目系統(tǒng)。具體流程：從消息中間件接收事件參數(shù)信息，獲取目標地震記錄的開始與結(jié)束信息，從redis波形共享內(nèi)存中獲取事件波形，經(jīng)過數(shù)據(jù)預處理、特征值提取等步驟，進行事件類型判斷；產(chǎn)出的事件類型信息，通過消息中間與其他軟件模塊共享。圖7為自動識別模塊與自動編目接口技術(shù)路線圖。

圖7 自動識別模塊與自動編目接口技術(shù)路線圖Fig.7 Technical roadmap of interface between automatic identification module and automatic cataloging

4.2.2 在線應用

目前自動識別模塊已具備在線運行的能力，初步實現(xiàn)事件類型的實時全自動判別。2019年11月至2020年3月已在線試驗運行5個月，其間檢測出多個人工編目漏分析地震事件，提高目錄完整性。對自動識別在線系統(tǒng)7×24運行穩(wěn)定性，進行不間斷測試，系統(tǒng)未出現(xiàn)死機、進程中斷現(xiàn)象，測試運行性能穩(wěn)定良好，功能和性能達到預期設(shè)計，驗證了自動識別模塊實用性和穩(wěn)定性。

5 討論和結(jié)論

（1）本文通過對多種特征組合聯(lián)合支持向量機進行大批量事件類型測試分析，結(jié)果顯示特征指標可對事件進行有效分類；據(jù)此得出識別效果較好的事件類型識別算法，事件類型識別準確率達到94.5%，靈敏度和特異度比較均衡，靈敏度和特異度均達到93%以上，對多數(shù)地震和爆破都能有效識別，極低震級單臺事件也能較好的分類，判識結(jié)果產(chǎn)出快，單個事件全流程識別僅需耗時2.5 s。

（2）算法數(shù)據(jù)輸入采用單臺Z分向和三分向數(shù)據(jù)作為輸入；采用同樣特征組合，三分向數(shù)據(jù)識別效果優(yōu)于單分向；可能的解釋是更多輸入數(shù)據(jù)對于支持向量機的輸入?yún)?shù)約束更好。此外，事件類型自動識別研究，訓練集中含有識別錯誤的事件也可能影響支持向量機分類性能；通常情況下，區(qū)域臺網(wǎng)地震目錄很少是100%正確；在對訓練集自驗證，最初訓練支持向量機已能準確識別出人工編目錯誤識別事件，其中地震訓練集中有3個錯誤識別事件，人工爆破中有8個錯誤識別事件；在最終測試前，對此數(shù)據(jù)進行剔除。

（3）采用本文研究最優(yōu)算法聯(lián)合支持向量機研制事件類型自動識別模塊，日常人工編目實際應用效果顯示，識別準確率達到93.9%，可檢測出人工識別錯誤事件，減少日常編目工作量和錯誤率。另外，實現(xiàn)與自動編目系統(tǒng)對接的在線自動識別系統(tǒng)接口，在線自動識別系統(tǒng)測試運行性能穩(wěn)定良好，檢測出多個人工漏分析地震事件，自動識別模塊具有實用性和穩(wěn)定性。

在后期工作中，增加其他區(qū)域事件進行完善和優(yōu)化學習樣本，并改進算法提高系統(tǒng)的泛化能力，自動識別模塊將會更具有普適性，實現(xiàn)事件類型識別的多樣性，塌陷、礦震、滑破等，可期應用于其他區(qū)域事件類型自動判別。