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

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

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

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

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

1 方法原理

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

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

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

(1)

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

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

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

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

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

1.2 模型

1.2.1LSTM模型

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

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

1) 細(xì)胞狀態(tài)(Ci)，存儲當(dāng)前時間步的特征信息，在所有鏈?zhǔn)接洃泦卧袀鬟f。

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

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

(2)

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

3) 輸入門，通過處理當(dāng)前時間步的輸入確定需要更新的信息并更新細(xì)胞狀態(tài)。首先，利用sigmoid函數(shù)生成更新程度it；其次，利用tanh函數(shù)生成含有待更新信息的候選向量Ct；最后用遺忘門的輸出ft與上一時間步的細(xì)胞狀態(tài)Ct-1相乘，丟棄需遺忘的信息；再將it與候選向量Ct相乘，得到更新后的信息，二者相加得到本時間步的細(xì)胞狀態(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) 輸出門，基于細(xì)胞狀態(tài)保存的信息選擇性的輸出細(xì)胞狀態(tài)Ct的內(nèi)容。通過sigmoid函數(shù)輸出信息被過濾的程度Ot，再利用tanh函數(shù)將當(dāng)前細(xì)胞狀態(tài)中的信息映射在[-1,1]之間，二者相乘后即得到當(dāng)前時間步的輸出h.其公式為：

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

(6)

ht=Ot?tanh(Ct) .

(7)

1.2.2改進(jìn)的雙向LSTM模型

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

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

r(l)=Bernoulli(p) .

(8)

(9)

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

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

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

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

1.3 訓(xùn)練過程

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

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

2.1 評價指標(biāo)

為了更好地對模型進(jìn)行評估，本文使用精準(zhǔn)率P和召回率R作為模型評價指標(biāo)。具體公式如下：

(10)

(11)

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

2.2 結(jié)果對比

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

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

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

圖6 訓(xùn)練過程的精確率和召回率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，由于人工拾取到時的標(biāo)注精度是0.1 s，故認(rèn)定該拾取結(jié)果準(zhǔn)確；圖8(b)中模型拾取的P波到時比人工拾取提前了1.47 s，在人工拾取到時的標(biāo)注精度之外，認(rèn)定該拾取結(jié)果誤差較大。其中，圖8(a)的數(shù)據(jù)噪聲信號較少，即信噪比較高；而圖8(b)的數(shù)據(jù)參雜過多的噪聲干擾，即信噪比較低。因此分析認(rèn)為本模型對信噪比較低的數(shù)據(jù)的拾取有一定的偏差，即在噪聲干擾過強(qiáng)的情況下可能會影響模型的拾取效果。

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

3 結(jié)束語

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

圖9 噪聲訓(xùn)練對比Fig.9 Comparison of noise training