基于LSTM網絡的單臺儀器地震烈度預測模型

李山有,王博睿,盧建旗*,王傲,張海峰,謝志南,陶冬旺

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

0 引言

大震會造成建筑結構的大量破壞及人員傷亡.為了減輕地震災害,科學家除了采取地震動參數區劃(Zhang and Jin,2008)、建筑抗震設計(Li et al.,2008)以及減隔振(Nepal and Saitoh,2020)等震前地震災害的預防措施以外,還提出了一種實時減災手段——地震預警系統(Allen et al.,2009).該系統首先識別破裂初期產生的P波信息(馬強等,2013; Lu et al.,2020; Wong et al.,2021),進而利用該信息對此次地地震發生的位置(金星等,2012; Satriano et al.,2008; Horiuchi,2005; Long et al.,2020)、地震的大小(彭朝勇等,2013; Peng et al.,2017; 李鴻杰等,2017; 朱景寶等,2022; Yang et al.,2021),以及地震可能造成破壞的范圍進行快速估計,為處在可能遭受地震破壞的區域內居民或重大工程發送預警信息,避免次生災害及人員傷亡.目前,地震預警系統已經在多個地震多發國家和地區開始了運行或試運行,如美國ShakeAlert預警系統(Cochran et al.,2018),日本JMA預警系統(Hoshiba et al.,2008),意大利PRESTo預警系統(Satriano et al.,2011),以及我國在建的國家地震烈度速報與預警系統等.地震預警系統主要有3種類型(Allen et al.,2009):(1)區域預警,區域預警在可能發生地震的區域內均勻布置觀測臺站.當區域內發生一次地震時,區域地震預警系統首先利用少數臺站接收到的P波信息估算本次地震的震源參數.例如:利用臺站接收到P波的到時信息進行地震定位,利用P波的峰值位移等參數進行震級估算.在此基礎上,利用地震動衰減關系估計此次地震可能造成建筑物破壞的影響范圍,并向該范圍內的用戶發布預警信息.(2)異地預警,異地預警將臺網布置在已知發震斷層附近,當地震發生時,利用布置在發震斷層附近的臺站觀測數據估計震源參數,結合地震動衰減關系估計預警目標區域可能遭受的破壞,為目標城市或區域發布預警信息.例如墨西哥地震預警系統就是典型的異地預警(Aranda et al.,1995).然而,人們并非對未來潛在震源的位置了如指掌,大部分可能發生地震的位置都是未知的.隨著地震觀測臺網建設成本的降低以及臺網密度的增加,這種預警方式逐漸被區域地震預警所取代.(3)現地預警,現地預警在重大工程廠址布設臺站,當臺站接收到的地震動超過一定閾值時,重大工程根據閾值采取合適的緊急處置措施(Nakamura,1988; Carranza et al.,2013).顯然,這種預警方法操作簡單,建設成本低,便于管理和維護.但是,閾值預警能夠為重大工程提供的預警時間非常有限.為了克服這一缺點,人們發現P波的峰值位移與地震動峰值之間存在一定的比例關系,在此基礎上提出了基于P波峰值位移的現地地震動預測方法(Wu and Kanamori,2005a,b; 劉辰等,2019; 宋晉東等,2021; Caruso et al.,2017).相比閾值預警方式,這種方法能夠為重大工程提供更長的緊急處置時間.

不論現地地震預警技術還是區域地震預警技術,均采用P波峰值位移等參數估計地震的最終大小或地震動的最終大小,存在大震低估及精度低的問題,從而可能導致漏報或誤報(Kodera et al.,2018).為了克服單一預警方式存在的缺陷,人們提出了融合多種方法的預警手段(Wu et al.,2017; Minson et al.,2017).例如,2011年日本311地震以后,Kodera等針對區域預警系統存在的漏報問題,提出了一種不依賴于震源參數估計的PLUM(The Propagation of Local Undamped Motion)預警方法.Caruso等人提出了融合區域臺網及單臺現地預測的地震預警原型系統(Caruso et al.,2017).盡管融合方法在一定程度上能夠解決漏報或誤報的問題,但不能解決預測大震低估及精度低的問題.

大量相關研究表明,斷層的破裂是一個復雜的時間過程(Lo et al.,2018; Wang et al.,2019; Liu et al.,2019; Tian et al.,2019),由于基于P波峰值及周期等參數的震級及地震動預測方法僅僅使用了P波的低維信息,沒法考慮斷層破裂過程蘊含的震源破裂過程信息.尤其P波的峰值位移、周期等參數在大震級段均出現了飽和(Trugman et al.,2019; Leyton et al.,2018),難以解決大震低估的問題.想要得到更好的估計效果,需要一種能夠考慮時間變化過程的地震動估計方法.而長短時記憶神經網絡(Long short-term memory,LSTM)則是一種能夠自動學習蘊含在時間序列中的規律的一種網絡結構,是目前解決這一問題最適合的一種手段.

因此,本文采用LSTM網絡,以地震動的能量、幅值及周期參數的變化過程作為輸入參數,以中國儀器地震烈度(GB/T 17742-2020)作為預測目標,利用日本K-NET臺網記錄的102次地震學習了能量變化過程與烈度的對應關系.結果表明,該方法能夠學到能量變化過程與最終烈度的關系,可以為地震預警中影響場的快速估算提供參考.

1 數據集的建立及參數提取

機器學習是一種數據驅動算法,數據集以及特征參數的選取對機器學習結果的好壞有著非常重要的影響.

1.1 數據集

本文選取了日本K-NET所記錄的震級在2～7.4級的191次地震,10457條三分量加速度記錄,圖1為K-NET數據臺站以及震中分布圖.選取102次地震作為訓練集,89次地震作為測試集,從這兩部分中的每次地震中隨機抽取部分記錄組成驗證集,保證各數據集中的震級、震源距分布盡可能接近.最終訓練集5103條記錄,測試集3781條記錄,驗證集1573條記錄(見圖2).

圖1 K-NET數據臺站及震源分布

圖2 震級與震源距直方圖

1.2 參數提取

(1) 預測目標

烈度是衡量地震動大小及建筑物破壞的一個綜合指標,常用于震害評估.為了能夠在一次地震發生以后快速獲得建筑物破壞情況以用于應急救援,人們通過地震動參數與宏觀烈度之間的統計關系建立了儀器地震烈度,用于地震烈度速報及應急救援.在地震預警系統中,儀器地震烈度也是地震預警信息發布的決策參數.因此,本文選擇了中國儀器地震烈度作為預測目標.

(2) 輸入特征參數

無論用那種手段進行地震動峰值參數的預測,地震動峰值參數的物理空間分布規律是所有方法首先要考慮的物理背景.人們在地震動衰減規律的研究方面已經取得了大量的研究成果,對地震動空間分布規律也取得了公認的認識.例如,忽略地震動的場地效應、近場飽和效應以及震級飽和效應等一系列影響因素之外,地震動峰值加速度(PGA)、P波峰值位移(Pd)以及烈度(I)等參數隨距離(R)的衰減滿足最簡單的衰減規律:

log10(PGA)=a0+a1M+a2log10(R)±σln(PGA),

(1)

其中,PGA為地震動峰值加速度(cm·s-2),R為臺站到震源或震中的距離(km),M為震級,a0、a1、a2為回歸系數,σln(PGA)為對數標準差.

log10(Pd)=b0+b1M+b2log10(R)±σln(Pd),

(2)

其中,Pd為P波峰值加速度(cm·s-2),b0、b1、b2為回歸系數,σln(Pd)為對數標準差.

I=c0+c1M+c2log10(R)±σI,

(3)

其中,c0、c1、c2為回歸系數,σI為標準差.

通過上述公式,我們就可以直接建立地震動峰值加速度或烈度與P波峰值位移之間的關系,進而利用P波峰值位移預測該臺站處的峰值加速度或烈度:

log10(PGA)=(a0-b0)+(a1-b1)M+(a2-b2)

×log10(R)+log10(Pd),

(4)

I=(c0-b0)+(c1-b1)M+(c2-b2)log10(R)

+log10(Pd),

(5)

由公式(4)、(5)可知,無論峰值加速度還是烈度與P波峰值位移的關系中,都無法忽略震級與距離的影響.由于不同地震動參數衰減規律的不同,衰減公式中的系數如a0與b0、a1與b1、a2與b2等并不相等,無法消除震級項與距離項,仍然需要考慮震級與距離參數的影響.因此,在機器學習的入參數選擇中,我們需要考慮哪些參數間接地蘊含了震級的信息和距離信息.在此基礎上,本文提出了以下幾個用于機器學習的輸入特征參數.

首先,震級是一次地震發震斷層破裂釋放的能量大小的一種量度,地震的震級越大、破裂面積越大、釋放的能量越多、能量的釋放速率也越大、斷層破裂的時間越長.因此,能量是間接體現地震震級大小的一個參數.本文選取了累積絕對速度(Cumulative absolute velocity,CAV)(Campbell and Bozorgnia,2012)作為一個臺站接收到由震源輻射的能量的衡量指標(如公式(6)所示).同時,我們選取了CAV增長率作為衡量斷層能量釋放速率的指標(如公式(7)).除此之外,本文還選擇了表征地震動能量大小的速度平方積分(Squared velocity integral,IV2) (Brondi et al.,2015)作為輸入特征參數(如公式(8)).另外,大量研究表明,地震的震級越大,產生地震動的卓越周期越長、P波的峰值位移越大,并建立了大量利用卓越周期或P波峰值位移的震級估算方法(Allen and Kanamori,2003; Olson and Allen,2005; Ziv,2014; Hildyard and Rietbrock,2010; Peng et al.,2017).因此,本文也選擇了卓越周期及峰值位移作為間接反映震級信息的輸入特征參數,其中,本文卓越周期采用了Hildyard等人提出的卓越周期實時計算方法(Hildyard and Rietbrock,2010).

其次,考慮到現在地震預警系統在快速定位方面的技術已經比較成熟,地震預警臺網的密度也在不斷增加,利用斷層破裂初期3～5個臺站接收到P波到時后就可以得到合理的定位結果(盧建旗等,2021),這為單臺地震動預測中引入距離信息提供了可能.因此,本文選取了震源距作為衡量能量空間衰減的特征參數.

最后,考慮到中國儀器地震烈度是由地震動加速度峰值和速度峰值計算得到,本文也選取了速度幅值、加速度幅值的絕對值作為輸入特征參數,計算方法見公式(9)—(10).

假設x=[x1,x2,…,xn]表示一條加速度記錄幅值的時間序列,則該記錄的CAV的時間序列CAV=[cav1,cav2,…,cavn]可以表示為

(6)

其中,Δt是加速度記錄的采樣間隔 (s).

CAV增長率的時間序列CAVr=[cavr1,cavr2,…,cavri]可以表示為

cavri=(cavi+0.3/Δt-cavi)/(0.3),

(7)

其中,cavi表示i時刻的CAV值.

速度平方積分定義為

(8)

速度幅值時間序列Pa的取值方法為

(9)

速度幅值時間序列Pv的取值方法為

(10)

v(t)為豎向加速度記錄的速度時程(cm·s-1),即對加速度時程進行一次積分后,再在進行4階巴特沃斯0.075～25 Hz帶通濾波得到速度時程,a(t)為豎向加速度記錄經4階巴特沃斯0.075～25 Hz帶通濾波得到的加速度時程;T為時窗長度.

2 方法

LSTM神經網絡屬于循環神經網絡的一個分支,是一種用于處理時間序列特征的網絡模型.這種網絡能夠學習蘊含在時間序列中的特征,可以建立時間序列特征與目標輸出之間的映射關系,在機器翻譯、情感分析以及交通流量分析等領域得到了廣泛的應用.為了解決傳統循環神經網絡在序列較長時容易出現梯度消失和梯度爆炸的問題,該模型引入了記憶單元和門控單元(Hochreiter and Schmidhuber,1997),可以實現對時間序列前段信息的記憶以及對冗余信息的遺忘,以避免梯度消失和爆炸的問題(Gers et al.,1999; Bengio et al.,1994).本文建立了一個用于儀器烈度預測的神經網絡模型(見圖3).該模型由3個LSTM層和1個全連接層組成,第一層由8個LSTM神經元組成,用于處理100個樣點的時間序列特征.為了加強神經網絡對前幾秒地震動信息的學習能力,我們采用了不等間隔的時間序列.即,0.2～10 s采用0.2 s的樣點間隔,10.5～20 s采用0.5 s的樣點間隔,21～50 s采用1 s的采樣間隔,共100個樣點,序列總長度為50 s.

圖3 用于預測烈度的LSTM神經網絡結構圖

本文LSTM神經元中采用的記憶單元及門控單元函數如公式(11)所示.采用均方誤差作為損失函數(公式(12)),利用“NAdam”優化器(Ruder,2017) 進行優化.

(11)

其中,ft,it和ot分別表示遺忘門、輸入門和輸出門的門限函數;ct表示狀態函數;ht表示隱含狀態函數;⊙表示逐元相乘算子;W表示系數矩陣,其下標第一項表示當前狀態,第二項表示所通過的門限;bf,bi,bo,和bg表示偏置向量;σ(*)=1/(1+e-*).

(12)

(13)

表1 部分時窗長度(WL)下評價指標統計表

(14)

(15)

其中,scorek表示第k個模型的得分.

圖4a及附表1展示了所有模型在測試集上的均方損失以及平均預警時間,圖4b為各個模型在測試集上的綜合評分.通過模型的綜合評分,最終確定的神經網絡模型參數為:層數L=3,輸出神經元個數U=8,批量大小B=100,優化器O為NAdam(Ruder,2017),最終確定的神經網絡結構如圖3所示.

圖4 模型超參數優化圖

損失函數是展現模型訓練成功與否的重要信息,圖5展示了本文模型在訓練過程中在測試集的損失函數與驗證集的損失函數隨訓練次數(Epoch)的變化曲線圖.圖中訓練集與驗證集的損失隨著訓練次數的增加而減小,而且訓練集與驗證集的殘差水平一致,說明神經網絡學到了數據中的普適性規律.

圖5 MSE損失隨訓練次數變化曲線

3 結果

將通過訓練集得到的神經網絡模型應用于測試集數據,其結果如圖6所示.考慮到烈度本身具有較大的不確定性,以預測誤差為1度作為分類依據將測試集預測結果同觀測烈度的比較結果分為三個區域,即準確區(綠色),高估(Overpredicted,OP)區(藍色)以及低估(Underpredicted,UP)區(橙色).圖中左上角為當前時間窗以及高估比例,右下角為預測標準差(STD)和低估比例.由圖6可知,當序列長度T為3 s時,VI度以下烈度的估計較為準確,大于VI度的烈度出現了部分低估,烈度被高估的比例很低.隨著序列長度的增加,高烈度低估的現象逐漸消失.由于序列的長度越長,序列中含有的信息越豐富,預測精度隨著時間序列長度的增加而提高,說明神經網絡學到了時間序列中蘊含的規律.

考慮到在烈度達到VI度后會造成建筑物的破壞,地震預警系統在預測烈度達到VI度及以上都會發送預警信息.為了評價本文模型提供預警信息的準確性,本文將預測結果定義為4類:

真陽性(TP):Iobs≥VI且Ipred≥VI,表示對大于VI度成功預測;

假陽性(FP):Iobs

真陰性(TN):Iobs

假陰性(FN):Iobs≥VI 且Ipred

由于地震記錄中大于VI度的樣本數量和小于VI度的數量嚴重不平衡,小于VI度的記錄占比較大,而大于VI度的記錄占比較小,簡單使用大于VI度或小于VI度被成功識別的比例難以客觀評價模型的精度.因此,我們用準確率(Accuracy)、漏報率(Miss alarm)和誤報率(False alarm)直觀有效地表現模型的預測能力.上述評價指標的定義為

(16)

準確率表示在所有樣本中,大于等于VI度和小于VI度被成功預測的樣本所占的比例.即,根據預測烈度發布的預警次數中有效預警所占的比例.漏報率表示在觀測烈度大于等于VI度的樣本中,預測烈度小于VI度的樣本所占的比例.誤報率表示在觀測烈度小于VI度的樣本中,預測烈度大于等于VI度的樣本所占比例.

為客觀評價本文模型的性能,我們沒有采用訓練集和驗證集的數據來計算三個評價指標,而是采用測試集數據計算了三個評價指標隨序列長度的變化關系曲線(如圖7所示).從圖中可見,隨著序列長度的增加,準確率值迅速增加,最終趨近于1,漏報率和誤報率最終趨近于0,說明模型的準確性較好.

圖7 LSTM-I模型評估曲線

表1列出了部分時窗長度下LSTM-I模型的準確率、誤報率以及漏報率.可見3 s時間窗時預測準確率為94.21%,誤報率為1.25%,漏報率為46.78%,至10 s時間窗準確率上升至97.84%,誤報率為1.14%,漏報率下降至17.6%.為了評價模型的泛化能力或是否出現了過擬合的情況,本文分別繪制了訓練集、驗證集以及測試集的烈度預測標準差隨序列長度的變化曲線(如圖8所示).由圖8可知,訓練集、驗證集以及測試集的預測烈度標準差非常接近,說明模型沒有出現過擬合,具有很好的泛化能力.

圖8 標準差隨時窗長度變化曲線

為了評價模型的時效性,本文計算了觀測烈度到達時間與預測烈度到達時間差(Lead time),到時差越大說明能夠為預警用戶提供的預警時間越長(如圖9所示).由圖9可見,到時差隨震源距的增加而增加.和理論P-S波到時差相比,二者較為接近.說明對于絕大部分記錄,模型在接收到P波信息后就對大于VI度的記錄提供了正確的預測.然而,本文的預測模型也出現了一些時差小于0的樣本,說明預測時間滯后.

圖9 LSTM-I訓練模型的時效性

4 對比

為考察本文建立的LSTM-I模型的先進性,利用本文選取的訓練集進行回歸得到P波觸發后3 s時間窗Pd-PGV與Pd-PGA模型,繼而通過《中國地震烈度表》(GB/T-17742-2020)給出的中國儀器地震烈度與PGA、PGV的換算關系得到Pd與烈度的對應關系.P波觸發后3 s的位移峰值Pd的取值方法如公式(17)所示:

(17)

其中,d(t)為豎向加速度記錄的位移時程(cm),即對加速度時程進行二次積分后,再進行4階巴特沃斯0.075～25 Hz的帶通濾波得到位移時程,Pd為P波前3 s的峰值位移(cm).

基于訓練集5147條地震動記錄,用最小二乘法回歸得到的Pd與PGV和PGA的關系(如公式(18)、(19),圖10所示),圖10中黑色實線為最小二乘法回歸結果,黑色虛線為±1倍標準差.

圖10 PGV與PGA關于Pd的回歸結果

log10(PGV)=0.823×log10(Pd)+1.757±0.427,

(18)

log10(PGA)=0.627×log10(Pd)+2.476±0.416,

(19)

圖11為利用測試集數據基于Pd方法計算得到的P波觸發后3 s時間窗的預測烈度和LSTM-I模型預測結果與實測烈度對比結果.圖11a為LSTM-I模型在P波觸發后3 s時間窗的預測結果與實際烈度的比較,圖11c為Pd預測結果與實測烈度的比較.對比a與c兩圖可見,基于Pd方法計算得到的預測結果高估率為20.93%,低估率為7.13%.而LSTM-I模型預測結果的高估率為3.52%,低估率為7.45%.說明LSTM-I模型的 “小值高估,大值低估”現象有明顯改善.

圖11 P波觸發3 s時間窗LSTM-I模型(a,b)與基于Pd模型計算預測(c,d)在測試集上的預測結果與實測烈度的比較

圖11b與圖11d分別展示了P波觸發后3 s時間窗時LSTM-I模型與Pd方法預測的誤差隨震源距的變化,以及誤差在震源距分布上的頻數直方圖和對應的概率密度分布曲線.左側散點圖內黑色實線表示誤差為0,虛線和點劃線分別代表±1倍標準差和±3倍標準差,右側灰色曲線表示概率密度分布.從圖11b可看出,LSTM-I模型預測烈度在0～500 km震源距區間內的誤差均勻分布于均值-0.1兩側,整體誤差位于±3倍標準差(-2.64,2.64)內.圖11d可看出,基于Pd方法得到的預測結果在0～500 km震源距區間內的誤差均勻分布于均值0.42兩側,整體誤差位于±3倍標準差(-3.57,3.57)內.繼而考察預測誤差在震源距分布上的頻數直方圖以及對應的概率密度分布曲線,兩方法對應的誤差均呈正態分布,且均值均在0附近,但LSTM-I模型的標準差小于Pd模型的標準差.

5 結論與討論

本文提出了一種基于時間序列特征的單臺儀器地震烈度持續預測LSTM神經網絡模型.該模型以累計絕對速度、累計絕對速度增長率、地震波卓越周期、加速度幅值、速度幅值、速度平方積分以及震源距為輸入參數,以單臺最大觀測儀器烈度為預測目標,利用日本K-NET強震臺網記錄的102次地震數據訓練了LSTM-I模型,用89次地震的數據作為測試集檢驗了模型的泛化能力.模型在訓練數據集與測試數據集上的預測誤差標準差變化趨于一致,這表明本文構建LSTM-I模型具備較好的泛化性.

從準確度分析,研究結果表明,P波觸發后3 s,預測烈度出現部分高烈度低估,很少出現低烈度高估現象.隨著時窗長度的增加,預測烈度的精度不斷提高.相較于基于常規Pd-PGV和Pd-PGA模型的計算方法準確性有所提升,可用于現地預警儀器地震烈度預測.

從時效性分析,本文計算了每條記錄VI度到達的時間與預測VI度到達的時間差,在不考慮系統延時的情況下可以看做該模型為用戶能夠提供的有效預警時間.結果表明,本文模型能夠提供的預警時間與P-S波到時差接近,說明模型具有較高的時效性.

由于烈度本身存在一定離散性,LSTM-I模型在P波觸發后3 s時窗的預測結果仍存在部分“大值低估,小值高估”的現象.雖隨著時間窗長度增加,低估和高估的現象大幅度減少,可得到較為準確的預測結果.

致謝使用開源平臺 Anaconda、Python、TensorFlow和Keras來設計神經網絡.日本防災科學技術研究所為本研究提供數據支持.

附表1 各個模型的訓練損失與平均預警時間