基于生成對抗網絡與隨機森林組合模型的地震與地脈動區(qū)分研究

劉赫奕， 宋晉東， 李山有

(1.中國地震局工程力學研究所 地震工程與工程振動重點實驗室，哈爾濱 150080;2.地震災害防治應急管理部重點實驗室，哈爾濱 150080)

隨著地震儀在全球范圍的大規(guī)模部署，人類已經進入地震大數據時代。如何處理地震網絡每天收集的大量數據是一個非常棘手的問題，特別是從質量參差不齊的連續(xù)波形記錄中分離出地震和噪聲。由于不同區(qū)域下不同類型地震事件的特征不同，地震事件的分類與篩選具有很大的地域性差異。目前這些工作仍主要靠低效的人工完成，因此在地震大數據的時代背景下，發(fā)展高效、高精度、普適性強的自動波形檢測算法顯得越來越重要。

傳統(tǒng)的地震檢測方法在推進地震數據自動檢測方面取得了很大進展，如長短時窗方法(short term averaging/long term averaging,STA/LTA)[1]、自回歸方法[2]和高階統(tǒng)計量法[3]，但由于僅采用波形的部分特征設定閾值，所以到時拾取精度比人工處理結果低，誤拾取率較高。根據波形相似性提出的波形自相關方法[4]、模板匹配方法[5]等，雖然可以檢測出大量微震，檢測靈敏度較高，但需要事先建立模板，且互相關計算量大，處理效率低，難以用于地震數據的實時處理。

隨著海量地震臺站的建設，人工智能已經引入地震事件識別的研究中。目前，人工智能在地震檢測[6-8]中已經取得了很多優(yōu)秀的成果。Kong等[9]利用人工神經網絡方法識別出Myshake系統(tǒng)中10 km內98%的地震記錄，從而高效準確地區(qū)分出地震事件與人類行為。江汶鄉(xiāng)等[10]利用觸發(fā)后3 s內的波形特征值訓練決策樹模型，建立了特征值和事件類型之間的映射關系，從而正確地區(qū)分強震記錄和干擾信號，達到了提高預警系統(tǒng)可靠性的目的?；诓煌貐^(qū)發(fā)生的不同類別的地震，Perol等[11]提出了針對地震識別和定位的卷積神經網絡，該方法成功在美國俄克拉何馬地區(qū)檢測出了超過原始編目17倍的地震事件。為了檢測不同時長的地震事件，Wu等[12]設計了一種名為DeepDetect的級聯(lián)卷積神經網絡，其識別事件的準確率可以達到63.8%，顯著超過模板匹配方法的準確率。相比于傳統(tǒng)的地震檢測手段，這些機器學習方法可以自動提取更多更接近數據本質的特征，從而得到更為準確的地震識別效果。

為了探究機器學習中無監(jiān)督學習算法對地震識別的效用，本文提出了一種利用生成對抗網絡(generative adversarial network，GAN)與隨機森林(random forest，RF)組合模型區(qū)分地震記錄與地脈動噪聲的方法。GAN是一種無監(jiān)督學習算法，與人工標注不同的是，它可以自動提取波形特征，簡化訓練過程。隨機森林是一種以決策樹為基學習器的二分類模型，它可以利用GAN提取出的波形特征進行分類，并在做出準確分類的基礎上有效解決過擬合的問題。通過該組合模型，本文將復雜的地震識別問題轉化為簡單的二分類問題，從而避免了識別過程中可能會出現的問題，旨在探索無監(jiān)督學習方法在地震監(jiān)測與地震預警中的應用可行性。

1 數 據

對于訓練集數據，本文選取了2007年10月—2017年9月共79 762條日本K-NET和KiK-net井上豎直方向的強震記錄，其中包括39 097條K-NET記錄和40 665條KiK-net井上記錄，震級范圍為3級≤M≤8級。首先，本文篩選出了信噪比大于5的記錄(信噪比定義為P波到時后2 s與前2 s的峰值速度比)，然后利用馬強等[13]的方法對記錄進行P波到時自動拾取，并對拾取后的到時進行人工校核。接著，本文篩選出了到時后1 s標準差大于前1 s標準差的加速度記錄。

在正式訓練之前，本文對篩選出的加速度記錄進行了如下的預處理：① 基線校正；② 2階0.075 Hz的巴特沃斯高通濾波；③ 利用插值方法將數據統(tǒng)一為100 Hz的加速度記錄；④ 截取P波到時前1 s至后3 s的加速度記錄作為地震數據；⑤ 截取KiK-net井上豎直記錄中起始點至P波到時點間的記錄，并隨機選取4 s波形作為地脈動數據；⑥ 將所有截取后的數據歸一化(除以最大幅值的絕對值)，圖1給出了經過預處理后的地震動和地脈動數據的實例。經過上述處理，本文共得到63 052段地震記錄和81 118段地脈動噪聲記錄(時長均為4 s)，地震記錄的臺站及震中分布圖如圖2所示。根據地震記錄的數量，本文隨機選取了63 052段地脈動噪聲，并隨機選取1/6的訓練集作為驗證集，其震級震中距分布圖如圖3(a)所示。

(a) 地震動記錄

(b) 地脈動數據

圖1 文中所用4 s數據的示例(橫坐標為0的點即為P波到時點，圖中數據均已經過歸一化處理。)

Fig.1 Examples of 4 s data used in this study (The point where the abscissa is 0 is the P arrival. The data in the figure has been normalized)

對于測試集數據，本文選取了獨立于訓練集之外的數據集，包含2020年1月—6月共7 774條K-NET豎直方向的強震記錄，震級范圍為3級≤M≤8級。經過與訓練集相同的篩選和預處理方法，我們得到了5 378條地震記錄和5 378條地脈動噪聲數據，所以測試集共有10 756條數據記錄，其震中和臺站分布圖如圖2所示，震級震中距分布圖如圖3(b)所示。

(a) 本文所用訓練集的臺站分布圖

(c) 本文所用測試集的臺站分布圖

(a) 訓練集震級震中距分布圖

(b) 測試集震級震中距分布圖

2 方 法

2.1 訓練GAN

GAN是Goodfellow等[14]提出的一種無監(jiān)督機器學習算法，由生成器和判別器組成。當向模型中輸入真實波形時，生成器會利用隨機樣本生成與真實數據相似的模擬波形，判別器負責區(qū)分真實波形和生成的模擬波形。圖4顯示了GAN訓練的基本過程，圖中粗點線為真實波形，實線代表生成波形，細點線代表判別器做出的決策，下方箭頭是用于生成波形的隨機樣本。圖4(a)表示訓練剛開始時，生成波形與真實波形差異較大，判別器可以很輕松地將它們區(qū)分開；圖4(b)和圖4(c)表示訓練時，生成波形與真實數據越來越接近，判別器漸漸混淆了兩種波形；圖4(d)表示經過多個步驟的訓練，生成器與判別器的能力均已達到飽和點，判別器無法區(qū)分兩種波形，此時訓練結束。

(a)

(c)

本文中GAN的架構如圖5(a)所示。生成器的輸入為50維服從正態(tài)分布的向量，輸出為4 s的生成波形。生成器主要包含兩個全連接層(神經元數目分別為128和400)和兩個Leaky ReLU(rectified linear unit)層。為了得到更好的訓練結果，將生成器最后一層的激活函數設置為tanh函數。判別器的輸入是4 s的波形，輸出是輸入波形為真實地震動的概率。如果輸入的是真實波形，則輸出為1；如果輸入的是生成波形，則輸出為0。判別器包含一組卷積池化層、兩個全連接層(神經元數目分別為128和1)和一個Leaky ReLU函數層。為了保持訓練的平衡，本文將判別器的學習率(0.000 2)設為生成器學習率(0.000 1)的2倍，且生成器每迭代1次，判別器迭代5次。本文將批量設置為64，所以生成器共有58 128個參數，判別器共有203 073個參數。本文中的參數設置遵循了其他訓練GAN的先例[15]，且可以在驗證集上獲得相對穩(wěn)定的結果。

(a)

2.2 區(qū)分P波與地脈動

因為GAN的輸入只有地震動記錄，所以它只能用于區(qū)分真實的地震動和生成的地震動，因此本文只利用GAN提取P波特征，并進一步利用RF分類器區(qū)分地震動與地脈動噪聲，這樣地脈動數據也可作為RF的輸入參與組合模型的訓練。RF是一種基于決策樹的分類器[16]，它利用bootstrap sampling[17]方法多次重復抽樣構成子集[18]，并融合多個子集決策樹的結果確定原始數據集的測試結果。因為它具有良好的二分類性能，所以可以充分利用GAN提取的特征對地震動和地脈動噪聲進行區(qū)分，從而得到更為準確的分類結果。

因此，當GAN訓練結束后，本文去掉了判別器最后一個全連接層，并將其余部分(圖5(a)中虛線框中的部分)組成一個P波特征提取器，如圖5(b)所示。本文假設訓練后的判別器已具有識別地震的能力，所以P波特征提取器可以提取出地震波的關鍵特征[19]。隨后，我們將10 756條測試數據輸入到P波特征提取器中提取波形特征，其輸出的128維向量即為RF的輸入數據。經過RF的訓練，本文可以得到最終的分類結果。如果輸入樣本為地震動數據，那么輸出為1；如果輸入樣本為地脈動數據，則輸出為0。

訓練時，RF中樹的個數為100，樹深為45，其余參數設置均遵循sklearn.ensemble.RandomForestClassifier中的初始設置[20]。除此之外，當RF的參數probability設置為True時，模型可以輸出一個對分類結果“自信程度”的概率，我們可以通過它了解測試結果的可靠性。

3 結 果

經過訓練，GAN與RF組合模型的分類性能如表1所示，其中的標簽含義如下：① 真陽性，組合模型將真實的P波判定為P波；② 真陰性，組合模型將真實的地脈動判定為地脈動；③ 假陽性，組合模型將真實的地脈動判定為P波，即為地震預警中的誤報；④ 假陰性，組合模型將真實的P波判定為地脈動，即為地震預警中的漏報。

從表1可以看出，組合模型可以識別出測試集中99.80%的P波和99.11%的地脈動，證明該方法可以準確地區(qū)分地震事件與地脈動噪聲，是一種較為可靠的識別地震事件的方法。此外，因為表1中假陽性樣本的比例為0.89%，表明該方法可以有效地降低由地脈動造成的誤報數量。

表1 測試集中組合模型的分類性能

本文還啟用了RF中的輸出概率，從而了解組合模型分類結果的可靠性。輸出概率是一個介于0和1的概率值，其中1表示真樣本(即P波)，0表示假樣本(即地脈動)。輸出概率越接近于1，該樣本為P波的可能性越大；反之，輸出概率越接近于0，該樣本為地脈動噪聲的可能性越大。本文將測試集中的輸出概率全部導出，并統(tǒng)計了輸出概率在各間隔之間(假設概率間的間隔為0.1)的樣本百分比，如圖6所示。橫坐標為P波的輸出概率，縱坐標為測試數據的百分比，所以在0.9～1.0概率段的P波樣本數越多，同時在0～0.1概率段的地脈動樣本數也越多時，模型的分類性能越好。本文可以通過圖中兩端的百分比檢驗分類結果的可靠性，并確定組合模型閾值與分類結果的關系。由圖6可以看出，測試集中輸出概率為0～0.1的地脈動數據占97.38%，輸出概率為0.9～1.0的地震波數據占94.51%，這表明絕大多數地震數據和地脈動數據的分類結果都具有較高的可信度，證明該組合模型可以準確地區(qū)分地震動與地脈動噪聲，從而識別地震事件。需要注意的是，本文將分類閾值設定為0.1和0.9，但在預警系統(tǒng)的實際應用過程中也可以通過臺站噪聲的觸發(fā)概率來調整此閾值。

圖6 測試數據百分比與P波輸出概率的關系(需要注意的是，地震動數據集中在0.9～1.0，地脈動數據集中在0～0.1。)

另外，本文分析了分類準確率與震級之間的關系，如圖7所示。通過觀察可知，所有震級段的準確率均超過了97%，且4～5級地震的準確率達到了100%，6級及6級以上地震的準確率稍差。所以，本文著重分析了6級及6級以上地震的分類效果。測試集共有267次6級≤M<7級的地震，34次M≥7級的地震，圖8表示了這301次地震的震中距與震級的關系。由圖7可知，6級≤M<7級和M≥7級地震的分類準確率分別為98.13%和97.06%，所以組合模型將其中的6次地震誤識別為地脈動，其中包含5次6級≤M<7級地震，1次M≥7級地震。表2顯示了這6次地震的地震名稱、震級和震中距信息，圖8中的叉表示了這6次地震震中距和震級之間的關系。由表2和圖8可知，在這6次誤識別的地震中，3次地震的震中距超過了200 km，且其它3次地震的震中距也超過了150 km。因為地震預警所研究的范圍多為近場地震，所以這6次地震的誤分類均在可允許的范圍內，所以證明該方法能夠較為準確地識別地震事件，減少由地脈動造成的誤報，而且在大震的近場臺站記錄中也同樣適用。

4 與傳統(tǒng)方法的對比

為了凸顯出本方法相比于傳統(tǒng)方法的優(yōu)越性，本文選取了使用范圍最為廣泛的STA/LTA方法與本方法進行對比，從而突出本方法在識別地震事件方面的優(yōu)勢。STA/LTA方法是Allen提出的一種能量方法，主要用于地震P波、S波初至時間的自動拾取，其主要原理是根據地震波形特征函數的長短時均比值等特征拾取初至[21]，目前廣泛應用于地震信號的檢測中，特別是對于地震弱信號。本文中采用了馬強等研究的計算公式

圖7 分類準確率與震級之間的關系(直方柱的顏色對應于該范圍內測試數據的數量)

圖8 測試集M≥6地震的震中距與震級的關系(圖中點為真陽性地震，叉為假陰性地震)

表2 測試集中誤識別為地脈動(假陰性樣本)的地震信息

(1)

式中：i為當前時刻點；k1和k2為i前某一時刻點；且k2

根據式(1)，在P波到來后，短窗內的平均值STA變化快，長窗內的LTA變化稍緩，所以STA/LTA刻畫了記錄幅值的瞬時變化，使特征變得更為明顯。如果設定合適的觸發(fā)閾值，當STA/LTA值超過設定閾值，則認為P波觸發(fā)。本文采取了馬強等關于STA/LTA方法的參數設置，相關閾值如表3所示。

表3 本文中STA/LTA窗長及閾值參數選取

因為STA/LTA方法的計算過程需要較長時間(至少30 s)的地震動記錄，所以本文并未對輸入數據進行截取。換句話說，本文只對輸入數據進行了前三步的預處理，隨后將整條地震動記錄直接輸入到STA/LTA方法中。而且，STA/LTA方法無需經過訓練過程，所以只需輸入測試數據就可以得到最終的結果。需要注意的是，因為STA/LTA方法輸出的是P波到時的具體時間，并不是分類的結果，所以本文將STA/LTA方法拾取的P波到時位于實際P波到時前1 s至后3 s內的記錄視為地震識別成功(輸出為1)的數據，其余位置均視為識別失敗(輸出為0)。

表4顯示了組合方法和STA/LTA方法對地震的識別效果對比。在5 378條測試地震動數據中，STA/LTA方法可以準確地識別出4 100條記錄，也就是說，STA/LTA方法可以識別出76.24%的地震動。而組合方法可以識別出99.80%的地震動，與傳統(tǒng)方法相比，本文的方法對地震動的識別率提高了23.56個百分點，說明該方法是一種可靠的地震動識別方法，其準確率要優(yōu)于傳統(tǒng)的自動撿拾算法。由于數據量巨大，本文比較了兩種方法在計算時間上的優(yōu)劣，如表5所示。由表5可見，組合模型在訓練集中所有記錄的訓練時間總和較長，但因為訓練過程可以預先完成，所以訓練時間并不直接影響測試過程中組合模型對單條記錄的計算處理時間。換句話說，地震發(fā)生后，組合模型在測試集中單條記錄的計算處理時間僅為0.000 93 s(可忽略不計 )，而STA/LTA算法的單條記錄計算處理時間為0.15 s。所以，從識別效果和單條記錄的計算處理時間來看，組合模型相比于傳統(tǒng)方法可以有效地提升識別地震動的準確性和時效性，是一種較為可靠的地震動識別方法。

表4 測試集中STA/LTA方法與本文方法的識別準確率對比

表5 測試集中STA/LTA方法與本文方法的計算時間對比

5 結 論

為了準確區(qū)分地震動與地脈動噪聲，本文利用日本K-NET和KiK-net井上地震訓練GAN，并將判別器改進為P波特征提取器提取波形特征，結合RF識別地震事件。測試結果表明，本文的方法可以區(qū)分出99.80%的地震波和99.11%的地脈動噪聲，相比于傳統(tǒng)的長短時窗方法提升了23.56個百分點，證明該方法是一種較為可靠的識別地震事件的方法，且在大震的近場臺站中也同樣適用。與手動標注波形特征相比，該組合模型可以自動提取出表征地震動的向量，從而簡化訓練步驟，節(jié)約訓練成本，提升訓練效率。

本文采用GAN對比真實地震動和生成的模擬波形，提取出能夠代表P波特征的概率分布，從而獲得真實數據中的某種隱含特性。GAN的判別器可以度量這種相似性，并通過與生成器的博弈縮短兩種波形的差距。但是本文并沒有利用判別器直接區(qū)分地震動與地脈動噪聲，是因為GAN的輸入只有真實地震動，而不包括地脈動數據，所以它并沒有真正接觸過真實的地脈動噪聲，區(qū)分結果不具有說服力。但與RF結合后，組合模型的輸入包含了輸入到RF中的地脈動數據。經過大量訓練之后，RF會逐漸學習并區(qū)分GAN提取出的兩種特征，最終區(qū)分出兩種波形，從而正確地識別地震事件。

雖然本文的組合模型可以比較準確地區(qū)分地震波與地脈動噪聲，但是在訓練過程中仍存在一些不足之處。一方面，組合模型需要在預先提取波形信號特征用以訓練模型，所以與傳統(tǒng)方法相比需要較長的訓練時間,但是這是在地震發(fā)生前預先訓練完成的，在地震發(fā)生后，測試過程單條記錄的計算處理時間可忽略不計，隨著計算機性能的提高，訓練時間的影響可能會減小，訓練效率上可以有所提高；另一方面，生成器生成的模擬波形作為副產品，暫時還未得到具體的應用。未來，本文將會朝著優(yōu)化現有模型或數據增強的方向發(fā)展，從而發(fā)掘本方法在地震學和地震預警中更大的應用潛力。

致謝

日本防災科學技術研究所(NIED)為本研究提供了數據支持，所有數據均從日本防災科學技術研究所官網下載(網址：http://www.kyoshin.bosai.go.jp/(最后登陸日期：2019年5月20日))，文中圖件使用通用制圖工具GMT(Genetic Mapping Tools)和Python中的matplotlib模塊繪制。