韓振華，郭浩雨，李 宇，張 玲，馮秀芳

(太原理工大學 信息化管理與建設中心，太原 030024)

震相到時拾取是地震學研究的重要課題之一。隨著數(shù)字地震臺網(wǎng)的大量建設，海量測震數(shù)據(jù)的產(chǎn)生，自動化快速精確識別震相對于地震學的研究有著重要的意義[1]。傳統(tǒng)的自動拾取算法如STA/LTA(長短時窗比)[2]利用短時窗平均值與長時窗平均值比值的變化來反映能量的變化情況，存在依賴閾值選取的問題；模板匹配算法[3]是根據(jù)波形相似性進行判斷，存在依賴模板、泛化性不好等問題。隨著深度學習的快速發(fā)展，卷積神經(jīng)網(wǎng)絡以及循環(huán)神經(jīng)網(wǎng)絡在震相拾取方面也得到了應用：PEROL et al[4]應用8層的卷積神經(jīng)網(wǎng)絡對美國70萬條連續(xù)地震波形數(shù)據(jù)進行了識別與定位，取得了較好的結(jié)果；趙明等[5]利用卷積神經(jīng)網(wǎng)絡對地震事件-噪聲問題進行分類，該方法可識別大量人工標記遺漏的微震事件；蔡振宇等[6]提出基于小波時頻譜的卷積神經(jīng)網(wǎng)絡模型對P波到時進行拾取，輸入地震波形數(shù)據(jù)后可直接輸出P波到時，具有運算速度快的特點；于子葉等[7]提出一種基于Inception網(wǎng)絡結(jié)構的震相拾取算法，該方法通過神經(jīng)網(wǎng)絡的自適應提取波形特征對P波和S波到時進行拾取，使噪聲數(shù)據(jù)有良好的穩(wěn)定性；ZACHARY et al[8]提出一種基于門控循環(huán)神經(jīng)網(wǎng)絡的震相連接模型，該方法可將P波和S波在起始時間至12 s內(nèi)發(fā)生的地震事件相關聯(lián)；CHIN et al[9]提出了一種基于循環(huán)神經(jīng)網(wǎng)絡的實地震波檢測算法，該算法可識別出實時發(fā)生的地震事件，同時可拾取P波和S波的到時。

RNN(循環(huán)神經(jīng)網(wǎng)絡)適合處理時間序列數(shù)據(jù)，被廣泛應用于信號處理。鑒于波形數(shù)據(jù)實際上是一個在時域維度且具有前后相關性的連續(xù)變化過程的數(shù)據(jù)，并且循環(huán)神經(jīng)網(wǎng)絡在處理時序數(shù)據(jù)上具有優(yōu)越性，因此將雙向LSTM模型作為本文震相拾取的基本方法。雙向LSTM模型可克服傳統(tǒng)RNN模型梯度消失的問題[10]，同時可以將波形數(shù)據(jù)過去與未來的特征進行整合，從而輔助決策震相到時的預測問題。

本文通過對原始地震數(shù)據(jù)提取波形特征并預測P波到時，提出一種基于改進的雙向LSTM模型的震相拾取算法。通過雙向LSTM模型對波形數(shù)據(jù)進行自適應特征提取，同時采用Spatial-Dropout機制在隨機區(qū)域?qū)W(wǎng)絡進行稀疏性約束以解決過擬合問題，最后通過引入Time-Distributed機制以充分利用波形數(shù)據(jù)前后時間步的特征信息并輸出P波預測到時。

1 方法原理

1.1 數(shù)據(jù)預處理

本文采用地震臺網(wǎng)人工標注震相到時的波形數(shù)據(jù)5 404條。從中隨機挑選4 324條作為訓練集，用于優(yōu)化模型參數(shù)；685條作為驗證集，用于判斷訓練停止時刻；395條作為測試集，用于驗證模型的泛化能力。

數(shù)據(jù)預處理對于神經(jīng)網(wǎng)絡的訓練有著重要的意義，預處理過程會直接影響模型收斂效果。本文的原始數(shù)據(jù)為時間窗口長度30 s的連續(xù)波形片段，數(shù)據(jù)的采樣頻率為100 Hz，因此數(shù)據(jù)格式為3 000×3的二維矩陣。由于波形數(shù)據(jù)中存在低頻的長期背景擾動以及高頻的噪聲，這兩者會導致模型在訓練過程中接受過多的無用信號，從而降低擬合效果，因此首先對原始數(shù)據(jù)進行帶通濾波。其次，地震數(shù)據(jù)中任意時間步的振幅差距過大，為了降低數(shù)據(jù)的方差，使損失函數(shù)的收斂過程更加平緩，需要對波形數(shù)據(jù)的每一道分量進行歸一化操作：

(1)

式中：Ai為每一道分量中第i個采樣點的振幅。

為避免原始數(shù)據(jù)較少而產(chǎn)生過擬合問題，在清洗后的波形數(shù)據(jù)中分別加入振幅最大值5%和10%的高斯噪聲形成加噪數(shù)據(jù)來擴增訓練數(shù)據(jù)量，最終形成16 212個地震事件數(shù)據(jù)集。圖1表示對原始波形數(shù)據(jù)進行加噪處理，這里繪制出隨機選取的兩個地震數(shù)據(jù)分別加入5%，10%高斯噪聲的對比圖。

此外，本文提取出原始波形數(shù)據(jù)中P波到時，并將其前后0.3 s作為允許誤差內(nèi)的P波時間窗標簽，其余時間步均標記為噪聲。這樣可將較長的時間序列數(shù)據(jù)預測某一時刻的回歸問題轉(zhuǎn)化為事件-噪聲的二分類問題，從而提高模型的準確率。圖2表示隨機選取一個數(shù)據(jù)的標簽處理前后對比，圖2(a)為原始數(shù)據(jù)，紅色線為P波到時；圖2(b)為處理后的P波到時標簽，可以看出處理后的P波到時信息與背景噪聲對比明顯，利于模型進行特征學習。

圖1 數(shù)據(jù)加噪前后對比Fig.1 Comparison of data before and after adding noise

圖2 數(shù)據(jù)標簽處理前后對比Fig.2 Comparison of data before and after label processing

1.2 模型

1.2.1LSTM模型

長短期記憶(long-short term memory，LSTM)模型[11]是一種優(yōu)化改進后的循環(huán)神經(jīng)網(wǎng)絡模型，可有效解決傳統(tǒng)循環(huán)神經(jīng)網(wǎng)絡在訓練過程中由長期依賴導致的梯度消失等問題[12]。LSTM由若干個記憶單元構成，其中每個記憶單元由細胞狀態(tài)、遺忘門、輸入門和更新門4部分組成。通過這些門的結(jié)構可讓信息選擇性的通過，用以去除或者增加信息到下一個細胞狀態(tài)。圖3表示LSTM單個記憶單元結(jié)構。

圖3 LSTM單個記憶單元結(jié)構Fig.3 Single memory cell structure of LSTM

1) 細胞狀態(tài)(Ci)，存儲當前時間步的特征信息，在所有鏈式記憶單元中傳遞。

2) 遺忘門，用于控制遺忘上一層細胞狀態(tài)的內(nèi)容，根據(jù)上一時間步的輸出hi-1和當前時間步的輸入xi，通過sigmoid函數(shù)將其映射為[0,1]之間的向量fi，表示上一層細胞狀態(tài)信息的保留程度，0代表全部遺忘，1代表全部保留。其公式為：

fi=σ(Wf[hi-1,xi]+bf) .

(2)

式中：W和b分別是LSTM模型的權值和偏置，均為模型的訓練參數(shù)；σ是激活函數(shù)sigmoid.

3) 輸入門，通過處理當前時間步的輸入確定需要更新的信息并更新細胞狀態(tài)。首先，利用sigmoid函數(shù)生成更新程度it；其次，利用tanh函數(shù)生成含有待更新信息的候選向量Ct；最后用遺忘門的輸出ft與上一時間步的細胞狀態(tài)Ct-1相乘，丟棄需遺忘的信息；再將it與候選向量Ct相乘，得到更新后的信息，二者相加得到本時間步的細胞狀態(tài)。其公式為：

it=σ(Wi[ht-1,xt]+bi) .

(3)

Ct=tanh(Wc[ht-1,xt]+bc) .

(4)

Ct=ft?Ct-1+it?Ct.

(5)

其中，?為逐點乘積。

4) 輸出門，基于細胞狀態(tài)保存的信息選擇性的輸出細胞狀態(tài)Ct的內(nèi)容。通過sigmoid函數(shù)輸出信息被過濾的程度Ot，再利用tanh函數(shù)將當前細胞狀態(tài)中的信息映射在[-1,1]之間，二者相乘后即得到當前時間步的輸出h.其公式為：

Ot=σ(Wo[ht-1,xt]+bo) .

(6)

ht=Ot?tanh(Ct) .

(7)

1.2.2改進的雙向LSTM模型

盡管LSTM模型可解決傳統(tǒng)RNN長期依賴的問題，但是沒有利用未來的信息[13]。為同時考慮數(shù)據(jù)的過去以及未來信息，本文采用雙向長短期記憶(Bi-Directional LSTM，Bi-LSTM)模型，其網(wǎng)絡結(jié)構如圖4所示。該模型的工作原理是通過前向LSTM層與反向LSTM層疊加，將原始序列數(shù)據(jù)輸入前向LSTM層并把數(shù)據(jù)的反向副本輸入給反向LSTM層，由此得到相反的兩個序列隱藏層狀態(tài)ht和h＇t，再將其連接得到預測輸出yt.

在數(shù)據(jù)有限的情況下，訓練參數(shù)過于龐大的神經(jīng)網(wǎng)絡會產(chǎn)生過擬合，即模型在訓練數(shù)據(jù)上完美契合，但泛化能力較差[14]。針對雙向LSTM模型對波形數(shù)據(jù)進行預測存在過擬合的問題，引入Spatial-Dropout機制，在模型的隨機區(qū)域進行稀疏性約束。Spatial-Dropout機制是2015年由TOMPSON et al[14]提出的一種Dropout方法，普通的Dropout會隨機地將部分元素置零，而Spatial-Dropout會隨機地將部分區(qū)域置零，對于地震波形這種時序數(shù)據(jù)來說，將一段信號隨機置為0可以減少過擬合的風險。其具體原理是以[0,1]之間的概率(一般為0.5)對部分區(qū)域的神經(jīng)元進行隨機失活處理，失活后的神經(jīng)元會在該次迭代中無法更新權重和偏置，由此使得模型無法通過記憶固有信息特征來完成震相到時的拾取，從而將模型變得稀疏化，解決了模型的過擬合問題。其公式為：

r(l)=Bernoulli(p) .

(8)

(9)

如圖4所示，為了對地震波形在時域上實現(xiàn)“多對多”的到時預測，通過Time-Distributed機制對每個時間步應用全連接層，通過設置sigmoid為激活函數(shù)即可從時域維度進行事件-噪音的二分類，從而輸出P波到時。

圖4 基于Bi-LSTM的P波拾取網(wǎng)絡結(jié)構圖Fig.4 Network Structure of P wave picking based on Bi-LSTM

基于雙向LSTM的P波拾取算法流程如圖5所示。整個流程由2部分組成：第一部分為數(shù)據(jù)預處理部分，對數(shù)據(jù)進行濾波、加噪與歸一化；第二部分為深度神經(jīng)網(wǎng)絡部分(見圖4)，由4層構成，分別是Input層、Bi-LSTM層、Time-Distribution層和輸出層。數(shù)據(jù)經(jīng)由預處理部分處理完畢后，由輸入層輸入預處理好的3 000×3的地震信號三分量數(shù)據(jù)(波形時長為30 s，采樣率為100 Hz)；Bi-LSTM層由雙向LSTM單元組成，實現(xiàn)地震信號進行時域分析。如圖4所示，Bi-LSTM層中，內(nèi)部有2個LSTM，分別是Forward層和Backward層，表示前向和后向，在Forward層從1時刻到t時刻正向計算一遍，在Backward層從末尾開始逆向計算一遍，得到并保存每個時刻正向與逆向隱含層的輸出。最后在每個時刻結(jié)合Forward層和Backward層的相應時刻輸出的結(jié)果進行拼接得到最終的輸出。每個LSTM單元的隱藏節(jié)點個數(shù)為128，輸出128維的特征向量。為了避免過擬合，利用Spatial-Dropout機制對模型進行稀疏化約束，對雙向LSTM輸出的特征使用ADD方式結(jié)合，輸出3 000×128的特征向量；Time-Distributed層對每個時間步的128維特征進行全連接分類；每個時間步的輸出結(jié)果進行Softmax運算后，由輸出層輸出3 000×1的One-Hot標簽，實現(xiàn)P波的位置的輸出。

圖5 基于改進雙向LSTM的P波拾取算法流程Fig.5 P wave picking algorithm flow based on improved Bi-LSTM

1.3 訓練過程

在模型訓練方面，本文使用均方根梯度下降算法(root mean square prop，RMSprop)作為優(yōu)化器，初始學習率為0.001并設置為學習率衰減，使用交叉熵損失函數(shù)(cross entropy error function)來對模型進行參數(shù)尋優(yōu)。此外，通過將訓練數(shù)據(jù)劃分為訓練集和驗證集，訓練集用于計算梯度以及更新參數(shù)，用驗證集來判斷訓練停止時刻，即在近5次迭代訓練時若出現(xiàn)驗證集誤差升高則停止訓練，同時返回模型參數(shù)。

2 實驗結(jié)果與分析

2.1 評價指標

為了更好地對模型進行評估，本文使用精準率P和召回率R作為模型評價指標。具體公式如下：

(10)

(11)

式中：TP為真正例，即模型識別的事件是真實的地震事件，相反FP則為假正例；TN為真反例，即模型識別的噪聲是真實的噪聲，相反FN則為假反例。精確率反映模型的誤檢率，精確率越高說明誤檢率越低；召回率反映模型的漏檢率，召回率越高說明漏檢率越低。

2.2 結(jié)果對比

經(jīng)過35個epoch的訓練，本模型最終在訓練集上精確率達到87%，召回率達到82%；在測試集上精確率達到90%，召回率達到92%，其訓練過程的精確率與召回率如圖6所示。為從不同角度驗證本文算法在震相到時預測問題上的適用性和優(yōu)越性，在測試集上與CAI et al[15]和ZACHARY et al[8]中的算法展開綜合對比，其對比結(jié)果如表1所示。

由表1可知，本模型在測試集上的精確率比LSTM高6%，比GRU高5%；在召回率上比LSTM高38%，比GRU高39%.表明本模型在不同的評估標準上均取得了更優(yōu)的結(jié)果。

此外，為了便于評價模型的泛化能力,將在測試集上預測的P波到時誤差分布做成直方圖，如圖7所示。圖7顯示395個測試集中的波形數(shù)據(jù)通過雙向LSTM模型計算得到的P波拾取結(jié)果與人工拾取到時之間的誤差分布，其中303條預測結(jié)果誤差小于0.1 s，占比76.7%；375條預測結(jié)果誤差小于0.2 s，占比94.2%；393條預測結(jié)果誤差小于0.5 s，占比99.5%；2條預測結(jié)果誤差大于1.0 s，占比0.5%.對其中到時拾取誤差小于0.1 s和大于1.0 s的數(shù)據(jù)分別繪制出波形圖進行分析，結(jié)果如圖8所示。其中紅色線為人工拾取的P波到時，藍色線為模型拾取的P波到時。

圖6 訓練過程的精確率和召回率Fig.6 Precision and recall rate of the training process

表1 不同算法在測試集上的比較結(jié)果Table 1 Comparison results of different algorithms on the test

由圖8分析可知，圖8(a)中模型拾取的P波到時比人工拾取滯后了0.02 s，由于人工拾取到時的標注精度是0.1 s，故認定該拾取結(jié)果準確；圖8(b)中模型拾取的P波到時比人工拾取提前了1.47 s，在人工拾取到時的標注精度之外，認定該拾取結(jié)果誤差較大。其中，圖8(a)的數(shù)據(jù)噪聲信號較少，即信噪比較高；而圖8(b)的數(shù)據(jù)參雜過多的噪聲干擾，即信噪比較低。因此分析認為本模型對信噪比較低的數(shù)據(jù)的拾取有一定的偏差，即在噪聲干擾過強的情況下可能會影響模型的拾取效果。

圖7 模型拾取到時與人工拾取到時的誤差分布Fig.7 Error distribution between model picking and manual picking

圖8 模型對高信噪比和低信噪比數(shù)據(jù)的到時拾取對比Fig.8 Comparison of the model's time picking of high SNR and low SNR data

為了驗證數(shù)據(jù)加噪對模型泛化能力的提升，本文做了噪聲訓練對比如圖9所示。從圖9中可以分析得出，用未加噪數(shù)據(jù)訓練的模型對加噪數(shù)據(jù)的P波拾取效果較差，而用加噪數(shù)據(jù)訓練后的模型對加噪數(shù)據(jù)的P波拾取效果更為準確。分析認為P波信號振幅較小，在噪聲干擾下容易被過多地覆蓋，因此在測試中未能有效識別其完整的波形特征，導致拾取效果不佳，而通過用加噪數(shù)據(jù)進行訓練可提高模型對P波周圍噪聲干擾的忍耐度，提高了模型的泛化能力。

3 結(jié)束語

本文提出一種基于改進雙向LSTM的震相到時檢測算法，首先對原始波形數(shù)據(jù)進行預處理，其次通過雙向LSTM模型對地震信號數(shù)據(jù)進行前向和反向的特征信息提取。利用Spatial-Dropout機制解決模型的過擬合問題，最后結(jié)合Time-Distributed機制實現(xiàn)P波到時預測。實驗結(jié)果表明，本文方法比LSTM、GRU等傳統(tǒng)RNN方法在震相到時預測問題上更具優(yōu)越性。同時，與CNN等方法不同的是，本文算法從時間維度將過去和未來的波形特征進行整合并且將振幅變化和時間的相關性聯(lián)系起來，因此大大減少了震相到時預測問題的誤差。此外，本文利用不同比例的加噪數(shù)據(jù)進行訓練測試，驗證了模型的泛化能力。最后，改進的雙向LSTM作為深度學習的方法之一，與STA/LTA、模板匹配等傳統(tǒng)波形檢測方法進行對比，本文模型不需要人工設定閾值，且對低信噪比數(shù)據(jù)有較高的準確率和召回率，其靈活性和普適性更強。

圖9 噪聲訓練對比Fig.9 Comparison of noise training