基于多臺站高頻GPS的地震預(yù)警優(yōu)化算法研究

張文浩，尹玲,2，胡文博

( 1. 上海工程技術(shù)大學(xué) 電子電氣工程學(xué)院, 上海 201620；2. 中國地震局地質(zhì)研究所, 北京 100029；3. 上海大學(xué) 通信與信息工程學(xué)院, 上海 200444 )

0 引 言

隨著我國經(jīng)濟建設(shè)的發(fā)展，大型高層建筑、人口密集城市不斷增多，我國作為地震活動最強烈的國家之一，地震造成的生命財產(chǎn)損失不計其數(shù)，如何減少地震造成的人員和財產(chǎn)損失，是目前迫切需要解決的問題. 由于地震的復(fù)雜性，當前的技術(shù)水平難以實現(xiàn)對地震的精準預(yù)測，而地震預(yù)警則是目前國際上公認的能夠減輕地震災(zāi)害的有效手段之一[1]. 傳統(tǒng)地震儀方法雖然比較成熟，能在地震發(fā)生后短時間內(nèi)確定震級和震源，但對強震震級進行估算時會出現(xiàn)超出量程、積分結(jié)果帶有誤差等現(xiàn)象，導(dǎo)致震級飽和[2].2011年3月11日，日本東北部發(fā)生Mw9.0級地震，日本氣象廳預(yù)警系統(tǒng)發(fā)布警報，但該預(yù)警系統(tǒng)的地震級別為Mw8.1，預(yù)測準確度有很大差距，同時確保正確預(yù)警并防止來自系統(tǒng)的誤報也是不小的挑戰(zhàn).2016年8月1日，日本地震預(yù)警系統(tǒng)發(fā)出Mw9.1級地震警報，預(yù)估日本大面積地區(qū)將達到烈度7度，該誤報造成了小田急電鐵和首都圈大部分地鐵緊急停運，一些重大工程緊急關(guān)停. 一旦發(fā)生誤報，會引起民眾恐慌，降低政府的信譽，造成不必要的損失.

地震預(yù)警系統(tǒng)(EEWs)基本可分為區(qū)域預(yù)警和單站預(yù)警[3]. 區(qū)域預(yù)警是基于密集的臺站檢測網(wǎng)絡(luò)，收集并分析來自多個臺站的數(shù)據(jù)估算相關(guān)源參數(shù)，用于遠距離地點提供早期預(yù)警. 單站預(yù)警則是由震源區(qū)域附近的單個傳感器組成，在某些系統(tǒng)中，為估計臺站區(qū)域的震動強度，還需確定震級和震源[4]. 單站預(yù)警通常比區(qū)域預(yù)警更快，但容易發(fā)生誤報，精確度較低[5-6]. 最合適的方式應(yīng)該是既保留區(qū)域預(yù)警的準確度，又盡可能實現(xiàn)單站快速預(yù)警.

現(xiàn)有的EEWs多是基于實時強震網(wǎng)絡(luò)框架建立的，系統(tǒng)中的強震儀以觸發(fā)模式進行工作，當達到一定強度后記錄數(shù)據(jù)并通過積分獲得位移結(jié)果. 但強震儀有基線偏差、量程有限等問題，而現(xiàn)代預(yù)警技術(shù)將高頻GPS技術(shù)引入地表位移監(jiān)控中. 它既可觀測靜態(tài)位移也可記錄由地震造成的動態(tài)位移，具有不限制量程，不積累誤差，能直接獲取位移數(shù)據(jù)等優(yōu)點. 文獻[7]闡述了高頻GPS技術(shù)在地震學(xué)中的應(yīng)用方法、實際意義及目前存在的問題. 文獻[8]利用高頻GPS數(shù)據(jù)使用最速下降法和OKADA模型進行反演，探討了高頻GPS的峰值地面加速度(PGA)、峰值地面速度(PGV)、峰值地面位移(PGD)計算地震烈度的可能性.

隨著地震檢測儀器在全球范圍的大規(guī)模部署，數(shù)據(jù)規(guī)模增大，如何完成從大量的連續(xù)波形數(shù)據(jù)中實現(xiàn)到震拾取是一個棘手的問題. 近年來人們發(fā)明了各種震相識別算法，比較典型的有：基于短時窗平均/長時窗平均算法(STA/LTA)的方法[9]. 信號到達時，STA/LTA值會發(fā)生突變，當比值大于某個閾值后，則判定事件發(fā)生，該方法應(yīng)用于信噪比(SNR)高的數(shù)據(jù)可達到較高準確率，SNR低的數(shù)據(jù)容易發(fā)生遺漏. 赤池信息量準則(AIC)方法[10]根據(jù)AIC尋求序列波形的極小點作為震相的到時點，該方法對震相和頻率的變化較敏感，常被用于震相到時的精確拾取，與STA/LTA類似，AIC的結(jié)果也強烈依賴于SNR以及檢測區(qū)間. 神經(jīng)網(wǎng)絡(luò)算法[11]可以利用神經(jīng)元捕捉輸入特征和輸出值之間的非線性關(guān)系. 文獻[12]利用神經(jīng)網(wǎng)絡(luò)對強震三分量數(shù)據(jù)進行到時提取，并加入噪聲提高模型泛化能力，展現(xiàn)出神經(jīng)網(wǎng)絡(luò)在自動化震相拾取方面的優(yōu)越性. 文獻[13]針對地震事件和噪音分類這一問題訓(xùn)練卷積神經(jīng)網(wǎng)絡(luò)，對序列片段進行分類.

本文探索用深度學(xué)習的方法對高頻GPS數(shù)據(jù)進行分析，將區(qū)域信息通過長短期記憶網(wǎng)絡(luò)(LSTM)[14]融合到單站預(yù)警中，盡可能降低誤報率同時保證快速反應(yīng). 分別使用融合區(qū)域特征的網(wǎng)絡(luò)模型、未融合區(qū)域特征的網(wǎng)絡(luò)模型、STA/LTA方法對新西蘭地區(qū)多個臺站無震高頻GPS序列進行誤報檢測，比較每種方法的誤報次數(shù)，另外為檢測算法的有效性對新西蘭凱庫拉地震進行檢出. 結(jié)果表明：融合區(qū)域特征后的模型不僅能夠?qū)Φ卣鹗录M行檢出，而且對比未融合區(qū)域特征的網(wǎng)絡(luò)模型和STA/LTA能有效地減少誤報次數(shù)，在應(yīng)用于高頻GPS時有很好地表現(xiàn)，可對當下的地震預(yù)警工作做出一定補充.

1 高頻GPS數(shù)據(jù)獲取與預(yù)處理

在選擇臺站時，盡可能選擇離震中距離均勻、記錄數(shù)據(jù)相對完整的臺站，然后進行數(shù)據(jù)解算. 文章使用2016年11月13日凱庫拉Mw7.8級地震[15]作為檢出事件，由于臨近該地震區(qū)域無震數(shù)據(jù)較少，本文選擇周邊區(qū)域記錄大量無震數(shù)據(jù)的臺站，利用6 h的無震數(shù)據(jù)去微調(diào)網(wǎng)絡(luò)，達到模型遷移的效果.

1.1 實驗數(shù)據(jù)

本文所使用的1 Hz高頻GPS數(shù)據(jù)來源于新西蘭地震信息網(wǎng)(Geonet)，該網(wǎng)絡(luò)包括Strong motion sensor、檢測區(qū)域形變的GPS sensor等臺站，剔除掉部分觀測數(shù)據(jù)質(zhì)量較差和數(shù)據(jù)缺失嚴重的臺站，選擇2016年2月中數(shù)據(jù)較為完整的3個臺站數(shù)據(jù)用于訓(xùn)練網(wǎng)絡(luò)，處理時間范圍為00:00:00—23:45:00 (UTC).選擇臨近地震區(qū)域3個臺站數(shù)據(jù)用于遷移學(xué)習，處理時間為2016-11-13 T 00:00:00—06:00:00 (UTC)，為了驗證本文所提出的方法在實際地震檢出應(yīng)用中的效果和可行性，本文對受凱庫拉Mw7.8級地震影響的多個臺站進行到震檢出，該地震發(fā)生于2016-11-13 T 11:02:59 (UTC)的新西蘭南島中部地區(qū)，震中位置位于(173.064 7°E，42.724 5°S)，震源深度為15.0 km. 選取2016-11-13 T 06:00:00—09:00:00 (UTC)的多個臺站無震時間序列作為測試數(shù)據(jù). 表1列出了用于遷移學(xué)習的3個臺站信息.

表 1 高頻GPS臺站信息表

1.2 數(shù)據(jù)處理方法

文章使用GAMIT track模塊對新西蘭南部地區(qū)1 Hz GPS數(shù)據(jù)進行解算，該模塊采用差分定位模式，實際獲取到的是目標站對于參考站的相對位置，所以參考站的選取會影響實際結(jié)果[15]. 對于震時數(shù)據(jù)的解算，如果基線距離過小，地震對參考站的影響會疊加在流動站位移序列上. 為盡可能避免同震干擾，選擇距震中較遠且觀測質(zhì)量較好的觀測站作為參考站，并檢查數(shù)據(jù)質(zhì)量，準備精密星歷文件，進行解算可得到所需臺站的時間序列數(shù)據(jù).

2 基于LSTM神經(jīng)網(wǎng)絡(luò)的地震預(yù)測模型

LSTM神經(jīng)網(wǎng)絡(luò)是一種特殊的循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)主要是解決長序列訓(xùn)練過程中的梯度消失和梯度爆炸問題[16]，在長序列中有更好地表現(xiàn).

2.1 LSTM神經(jīng)網(wǎng)絡(luò)原理

LSTM與RNN相比不再是單一 t anh 層，而是多個門結(jié)構(gòu)，分別是遺忘門、輸入門、輸出門，LSTM單元結(jié)構(gòu)如圖1所示.

圖 1 LSTM單元結(jié)構(gòu)圖

2.1.1 遺忘門

遺忘門會使用當前時刻的輸入xt和上一時刻的輸出ht-1通 過σ函數(shù)計算上一時刻細胞狀態(tài)Ct-1中信息保留的比重，0表示不保留，1表示保留. 其中Wf表示遺忘門的權(quán)重，bf為偏置項.

2.1.2 輸入門

輸入門可以利用當前時刻的輸入xt和上一時刻的輸出ht-1通 過σ函數(shù)決定更新哪些信息. 再利用xt和ht-1通 過 t anh 層得到臨時細胞狀態(tài)Wi表示輸入門的權(quán)重，bi為偏置項，Wc表示細胞單元的權(quán)重，bc為偏置項.

將舊的細胞狀態(tài)Ct-1更 新為新的細胞狀態(tài)Ct.

2.1.3 輸出門

輸出門可以利用當前時刻的輸入xt和上一時刻的輸出ht-1通過σ函數(shù)判斷輸出細胞的哪些特征，并通過與 t anh 層的細胞狀態(tài)Ct相乘得到單元輸出ht.Wo表示輸出門的權(quán)重，bo為偏置項.

2.2 網(wǎng)絡(luò)結(jié)構(gòu)

一個LSTM網(wǎng)絡(luò)包含多個隱藏層，每個隱藏層由大量LSTM神經(jīng)元組成，展開后等效于在不同層之間傳遞參數(shù)的前反饋神經(jīng)網(wǎng)絡(luò). 每一層的輸出作用于下一層的輸入，所以它可以處理時間序列數(shù)據(jù)并提取特征，LSTM預(yù)測模型整體框架如圖2所示. 為適應(yīng)地震檢出任務(wù)使用單層LSTM，使得網(wǎng)絡(luò)能夠?qū)W習高層次的時域特征，輸入層為多個臺站三分量數(shù)據(jù)，輸出為單臺站下一秒預(yù)測值.

圖 2 LSTM時間序列預(yù)測框架

2.3 模型計算流程

在網(wǎng)絡(luò)訓(xùn)練前需要對數(shù)據(jù)進行處理，主要包括劃分數(shù)據(jù)集、將高頻GPS數(shù)據(jù)轉(zhuǎn)換成可以輸入LSTM神經(jīng)網(wǎng)絡(luò)的格式，即樣本、時間步、特征的三維數(shù)據(jù).對輸入的多臺站數(shù)據(jù)進行歸一化處理，將高頻GPS時間序列轉(zhuǎn)換化成監(jiān)督學(xué)習數(shù)據(jù)，每個時間步前t時刻的數(shù)據(jù)作為訓(xùn)練樣本x，當前t時刻的數(shù)據(jù)作為輸出值y. 數(shù)據(jù)處理完畢后，由輸入層輸入處理好的多臺站地震信號數(shù)據(jù)，隱藏層LSTM單元實現(xiàn)對地震信號進行時域分析，通過Dropout層避免過擬合，全連接層對每個時間步的特征進行全連接，再將每個時間步的輸出結(jié)果經(jīng)過激活函數(shù)后，由輸出層及反歸一化操作輸出預(yù)測值，期間數(shù)據(jù)不斷迭代，以降低損失函數(shù)值為目標直至收斂. 最后根據(jù)預(yù)測值動態(tài)生成閾值，通過真實值與閾值相比較，判斷是否發(fā)生地震事件，算法流程如圖3所示.

圖 3 LSTM模型計算流程圖

2.4 模型訓(xùn)練

訓(xùn)練預(yù)測模型的目的是讓模型學(xué)習無震時的數(shù)據(jù)特征，并且用來預(yù)測未來時刻無震的數(shù)據(jù)值，通過預(yù)測值與觀測值相比較來檢出異常. 使用2016年2月中篩選的無震時間序列片段，最終選擇使用YALD、VWMR、VBRM 3個臺站數(shù)據(jù)構(gòu)造數(shù)據(jù)集，選取85%的數(shù)據(jù)作為訓(xùn)練集，15%的數(shù)據(jù)作為測試集，隱藏層使用單層LSTM網(wǎng)絡(luò)，并在每層之間加Dropout層，在訓(xùn)練時以抑制部分神經(jīng)元的輸出，避免在訓(xùn)練過程中過擬合. 最后一層是全連接層，使用線性激活函數(shù)，用于輸出模型下一秒的預(yù)測值. 網(wǎng)絡(luò)的優(yōu)化算法選取為自適應(yīng)運動估計算法(Adam)算法，該算法能使收斂過程更加平穩(wěn)，可以為不同的參數(shù)自適應(yīng)學(xué)習率. 使用均方根誤差(RMSE)作為衡量模型性能的指標，度量真實值與預(yù)測值之間的相對偏差，其計算如式(8)所示. 其中yi(t) 是 在t時刻的實際觀測值，(t) 是在t時刻模型輸出的預(yù)測值，n是數(shù)據(jù)樣本的總數(shù).

訓(xùn)練所用硬件為AMD Ryzen 3900x 3.8 GHz CPU、16 GB內(nèi)存、NVIDIA GeForce GTX2080ti 11 GB顯存. 在選取batch size參數(shù)的過程中，batch size過小會導(dǎo)致隨機性更大，難以達到收斂，但是能使模型學(xué)到一定噪聲，提高模型的泛化能力. 過大的 batch size會增加迭代過程中的計算量，內(nèi)存消耗大. batch size選擇的影響如表2所示.

表 2 不同batch size對結(jié)果的影響

batch size大小對模型學(xué)習有很大影響，當batch size足夠大時才可以對整個數(shù)據(jù)的梯度進行穩(wěn)定的計算，但過大的batch size會在迭代的過程中增加計算消耗. 過小的batch size能學(xué)習到一定噪聲，可以避免局部最優(yōu)，但可能會由于樣本不足，網(wǎng)絡(luò)權(quán)重突然波動而導(dǎo)致收斂緩慢甚至不收斂，從而對計算造成負面影響. 由表2可以看出，在batch size小的情況下單次迭代時間長，batch size大的情況雖能加快計算，但是考慮到精度，選擇batch size為64.

網(wǎng)絡(luò)進行訓(xùn)練得到的Loss圖像如圖4所示，實線為訓(xùn)練集Loss值變化，虛線為驗證集Loss值變化.訓(xùn)練Loss曲線從0.009開始減少，初始時收斂較快，經(jīng)過約15次迭代后收斂速度降低，最終收斂于0.000 6.

圖5為模型應(yīng)用于測試數(shù)據(jù)集得出的YALD站預(yù)測圖，計算得出RMSE為0.003 429，可以看到LSTM神經(jīng)網(wǎng)絡(luò)可以對臺站時間序列進行有效預(yù)測，可以充分挖掘區(qū)域各個臺站之間的相關(guān)性.

得到滿足預(yù)測精度的模型后，用地震周邊區(qū)域臺站無震數(shù)據(jù)微調(diào)模型. 選擇臨近地震區(qū)域MRBL、WITH、WEST 3個臺站6 h數(shù)據(jù)構(gòu)建數(shù)據(jù)集用于遷移學(xué)習. 應(yīng)用于測試集的預(yù)測結(jié)果如圖6所示，實線為真實觀測數(shù)據(jù)，虛線為來自預(yù)測模型的預(yù)測結(jié)果.由圖可知預(yù)測模型可以準確預(yù)測到觀測數(shù)據(jù)的變化趨勢，最終測試集的RMSE為0.002 989.

圖 4 LSTM網(wǎng)絡(luò)模型Loss值變化曲線

圖 5 YALD臺站預(yù)測圖

圖 6 MRBL臺站預(yù)測圖

圖 7 MRBL臺站地震檢出

圖 8 WITH臺站地震檢出

圖 9 WEST臺站地震檢出

3 模型應(yīng)用

利用本文設(shè)計的地震檢測模型可以對包含地震事件的GPS時間序列進行檢出，并對多條無震長序列進行檢測統(tǒng)計誤報次數(shù).

3.1 方法有效性和準確性

為驗證本文模型可對異常事件進行檢出，對包含凱庫拉地震事件的時間序列進行檢測，模型學(xué)到的特征為無震時的區(qū)域特征，而發(fā)生地震后序列會產(chǎn)生突變，突變?nèi)暨B續(xù)超出置信區(qū)間，則此次突變可作為異常實現(xiàn)檢出. 在對多個受地震影響的臺站進行檢出，MRBL臺站序列檢出結(jié)果如圖7所示，在第203時間步時觀測值發(fā)生明顯變化且超出置信區(qū)間，判定為異常實現(xiàn)地震事件，WITH臺站序列檢出結(jié)果如圖8所示，在第141時間步時觀測值發(fā)生異常，后續(xù)序列波動明顯且多次超出置信區(qū)間，實現(xiàn)到震檢出. WEST臺站序列檢出結(jié)果如圖9所示，在第202時間步觀測值異常變化超出置信區(qū)間，此時模型生成的閾值也跟隨觀測值發(fā)生一定突變，但模型結(jié)合區(qū)域之前的地表位移特征，使閾值變化不如實際值大，當實際值超出閾值時實現(xiàn)到震檢出.

對無震時序進行誤報檢測時，正常情況下觀測值不超出置信區(qū)間，若超出置信區(qū)間則判定為一次誤報，文章對多個臺站時間序列進行誤報檢測. 使用本文方法對MRBL臺站10 800 s無震序列進行檢測，檢測所用時間為31.74 ms，最終未見誤報；對WITH臺站10 800 s無震序列進行檢測，未見誤報；對WEST臺站10 800 s無震序列進行檢測，出現(xiàn)1次誤報. 另使用本文模型嘗試對周邊其他臺站3 h無震序列進行檢測，結(jié)果HANM、KAIK、LKTA臺站均未發(fā)生誤報，檢測圖如圖10所示.

基于雙軌跡構(gòu)圖可知本題有兩解，解法1漏解.構(gòu)造△ADB的外接圓圓心E，點E可在AB下方如圖4；也可在AB上方，如圖5；因此有必要思考圓心E還可在AB邊上，如圖6.

圖 10 引入?yún)^(qū)域信息模型無震序列誤報檢測

3.2 對比及分析

用STA/LTA方法和未引入?yún)^(qū)域特征的單站方法對多段無震時間序列進行檢測. STA/LTA方法對振幅變化敏感，但易受噪聲干擾，STA記錄的主要是振幅的瞬時變化，LTA的作用是衡量周圍噪聲的強度.當?shù)卣鹗录l(fā)生時，STA/LTA值明顯發(fā)生變化，當STA/LTA值超過閾值，則判斷此時刻為地震到時. 如果閾值設(shè)定過大則會遺漏地震事件，若閾值設(shè)定過小，則會因為噪聲或其他原因發(fā)生誤報，閾值不易設(shè)定.

由于高頻GPS SNR較高，STA/LTA方法用于高頻GPS時間序列數(shù)據(jù)效果并不理想，對多個臺站時長10 800 s的無震時序進行檢測，檢測一次所用時間約為0.73 ms. 圖11為STA/LTA方法檢測圖，虛線為觸發(fā)檢測器標記，由圖11可知STA/LTA方法在多個臺站上均誤報較多，具體誤報次數(shù)如表3所示.

未融入?yún)^(qū)域特征的單站模型數(shù)據(jù)復(fù)雜度較低，過深的LSTM網(wǎng)絡(luò)在訓(xùn)練時不易收斂. 在搭建網(wǎng)絡(luò)時減少每層網(wǎng)絡(luò)的神經(jīng)元數(shù)量，防止過擬合. 網(wǎng)絡(luò)架構(gòu)及相關(guān)參數(shù)如表4所示. 選取最優(yōu)模型對長10 800 s的多個臺站無震時間序列進行檢測，檢測一次所用時間約為23.59 ms，無震時間序列的檢測結(jié)果如圖12所示. 各個臺站的檢測結(jié)果均出現(xiàn)少量誤報，較STA/LTA方法有所減少，具體檢測結(jié)果如表3所示.

圖 11 STA/LTA檢出結(jié)果

圖 12 未引入?yún)^(qū)域信息模型無震序列誤報檢測

由表3可以看出，兩種神經(jīng)網(wǎng)絡(luò)方法在多個臺站的無震序列上均減少了誤報數(shù)量，而本文方法在引入?yún)^(qū)域特征后誤報數(shù)量再次減少，雖然相比于STA/LTA方法，神經(jīng)網(wǎng)絡(luò)模型計算所需時間變長，考慮到神經(jīng)網(wǎng)絡(luò)的參數(shù)量，計算所需時間必然會增加.

表 3 不同方法檢測所用時間及誤報數(shù)

表 4 網(wǎng)絡(luò)架構(gòu)及參數(shù)表

4 結(jié)束語

近些年來越來越多的連續(xù)運行的GPS接收器將記錄速率提升到1 Hz，高頻GPS在地震領(lǐng)域的應(yīng)用不斷發(fā)展，GPS數(shù)據(jù)不論是應(yīng)用于長期位移觀測還是短期地殼瞬時形變都有不錯的表現(xiàn)，在震級估算方面更為精確. 本文基于高頻GPS數(shù)據(jù)提出了一種包含區(qū)域信息的神經(jīng)網(wǎng)絡(luò)模型，在應(yīng)用于較高噪聲的GPS數(shù)據(jù)時能得到比傳統(tǒng)方法更好的效果.