基于支持向量機的現(xiàn)地地震預警地震動峰值預測

余 聰，宋晉東，李山有

(1.中國地震局工程力學研究所，哈爾濱 150080;2.中國地震局地震工程與工程振動重點實驗室，哈爾濱 150080)

地震是對人類生存安全危害巨大的一種突發(fā)性自然災害。地震預警是在地震發(fā)后，利用震源附近地震臺站觀測到的地震波初期信息，快速估計地震參數(shù)并預測地震對周邊地區(qū)的影響，搶在破壞性地震波到達震中周邊地區(qū)之前，發(fā)布各地地震動強度和到達時間等預警信息，能夠給重要的工程或人員密集地區(qū)提供幾秒到數(shù)十秒的預警時間，以減輕地震災害[1]。地震預警已被證明是實時減災的有效手段之一，經過幾十年的發(fā)展，世界上多個地震活躍的國家已運行或正在發(fā)展地震預警系統(tǒng)，例如墨西哥[2]、日本[3-5]、美國[6-9]、意大利[10-12]、中國[13-14]。

地震預警可以分為異地預警與現(xiàn)地預警[15]。異地預警在潛在震源區(qū)布設高密度的地震臺網，利用近震源臺(P波或S波)記錄信息來確定地震參數(shù)，進而估計預警目標區(qū)地震動，并對地震潛在破壞區(qū)發(fā)出警報；現(xiàn)地預警原理是利用初至P波和后至破壞性地震波(S波或者面波)的到時差，Kanamori[16]認為P波攜帶地震信息，因此可從P波中預測后續(xù)的峰值地震動，如PGA與PGV，通常的方法是建立單個P波特征參數(shù)與目標地震動參數(shù)的線性關系進行后續(xù)的地震動預測，例如Wu等[17]、Zollo等[18-19]使用的Pd方法，即采用P波觸發(fā)后3 s內的垂直方向位移峰值Pd預測后續(xù)地震動峰值PGV，這是地震預警最常用的一種現(xiàn)地地震動預測方法。此外，Nakamura[20]采用破壞烈度DI方法，F(xiàn)esta等[21]、Brondi等[22]采用的速度平方積分IV2方法等預測后續(xù)地震動。但是，利用單個P波參數(shù)預測后續(xù)地震動只利用了單一的地震動初期信息，B?se等[23]利用人工神經元網絡、基于多頻帶多參數(shù)實現(xiàn)了地震預警PGA、PGV的持續(xù)預測方法，為綜合利用多種地震波初期信息進行現(xiàn)地地震動預測提供了新的思路。

本文以0.1～10 Hz帶通濾波后三分向矢量合成地震動峰值PGA與PGV為預測目標參數(shù)，利用日本K-net強震臺網P波觸發(fā)后3 s數(shù)據(jù)，基于加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里葉譜幅值AMmax、速度平方積分IV2、破壞烈度DI、累積絕對速度CAV、阿里亞斯烈度Ia這8種特征參數(shù)構建支持向量機(Support Vector Machine，SVM)預測模型，旨在利用支持向量機樣本需求量小、泛化能力強和抗干擾能力強等優(yōu)點，探索多參數(shù)預測后續(xù)地震動的新方法，為現(xiàn)地地震預警儀器地震烈度預測提供服務。

1 數(shù) 據(jù)

本研究從日本K-net強震觀測臺網中選取地震數(shù)據(jù)，并依據(jù)以下原則篩選[24]:

(1)時間為自2007年10月至2019年8月；

(2)發(fā)生于陸地的地震，因為我國的地震多發(fā)生于陸地；

(3)震源深度為10 km以內的地震，因為我國大部分地區(qū)發(fā)生的是淺源地震，同時，在震級相同的情況下，淺源地震會比深源地震造成更大的破壞性；

(4)震中距取100 km以內的強震數(shù)據(jù)，對5級以上的地震，補充100～200 km范圍的數(shù)據(jù)。

基于上述原則共篩選出2 068次地震共25 308組(75 924條)強震動數(shù)據(jù)，豎直向原始加速度記錄未濾波幅值范圍0.14～665.78 cm/s2，水平向原始加速度記錄未濾波幅值范圍0.78～1 538.16 cm/s2。隨機選取訓練數(shù)據(jù)集與測試數(shù)據(jù)集，比例為7:3。地震震中及臺站的分布，如圖1所示。震級、震中距、記錄數(shù)量分布，如圖2所示。

圖1 本研究訓練數(shù)據(jù)集與測試數(shù)據(jù)集的地震震中分布(空心圓直徑與震級大小成正比)K-net臺站(灰色三角形)分布

圖2 本研究訓練數(shù)據(jù)集與測試數(shù)據(jù)集的震中距及震級分布圖(三角形為訓練數(shù)據(jù)，圓圈為測試數(shù)據(jù))

對上述數(shù)據(jù)進行如下處理：

(1)采用馬強等[25]的方法自動識別P波震相及到時，即先利用長短時平均法對P波震相初步識別、再利用AIC準則精確撿拾P波到時。為確保P波到時的準確性，對自動識別的地震P波到時進行人工校正；

(2)對加速度積分，積分一次得到速度記錄、積分兩次得到位移記錄。為了消除積分后的低頻漂移現(xiàn)象，對未處理加速度記錄和積分后的記錄作0.075 Hz高通濾波。

(3)對加速度記錄與積分后的速度記錄做0.1～10 Hz帶通濾波，進行三分向矢量合成計算得到地震動峰值PGA與PGV。

2 方 法

2.1 預測目標參數(shù)

地震烈度作為標定地震引起的地震動及其影響的強弱程度，具有度量地震破壞原因和后果的雙重性，其快速的評定對于震后災害評估、應急救援有著十分重要的作用。但是限于烈度本身的定義模式，震后難以組織人員對受災區(qū)域進行快速評定。隨著地震觀測技術的持續(xù)發(fā)展及地震觀測儀器數(shù)量的不斷增加，根據(jù)地震觀測儀器獲取的地震動記錄快速定量地計算儀器地震烈度成為可能。儀器地震烈度的計算是地震烈度速報的基礎，可為災情快速判斷、地震應急救援決策和行動、工程搶險修復決策等提供科學依據(jù)[26]。

《中國地震烈度表》[27]給出了儀器地震烈度的計算過程：

獲取地震動記錄所用觀測儀器應可為三方向(條件不具備時可為正交兩水平方向或單水平方向)，地震動加速度和速度記錄的每個分向均應采用數(shù)字濾波器進行0.1～10 Hz帶通濾波，PGA、PGV分別為濾波后三分向的加速度、速度矢量合成的最大值，那么應用PGA得到的儀器地震烈度計算值IPGA可由式(1)得到，應用PGV得到的儀器地震烈度計算值IPGV可由式(2)得到，進而可由IPGA和IPGV根據(jù)式(3)計算出儀器地震烈度II。

IPGA=

(1)

IPGV=

(2)

II=

(3)

2.2 特征參數(shù)

本研究采用如下8種特征參數(shù)作為支持向量機預測模型的輸入?yún)?shù)，且計算時間窗均為P波觸發(fā)后3 s。

(1)加速度幅值Pa：P波觸發(fā)后3 s內的垂直方向加速度峰值，該參數(shù)表示地震動的加速度幅值特征。

(2)速度幅值Pv：P波觸發(fā)后3 s內的垂直方向加速度峰值，該參數(shù)表示地震動的速度幅值特征。

(3)位移幅值Pd：P波觸發(fā)后3 s內的垂直方向位移峰值，該參數(shù)表示地震動的位移幅值特征。

(4)傅里葉譜幅值AMmax：P波前3 s垂直方向記錄傅里葉譜的幅值，該參數(shù)表示地震動的頻譜特征。

(5)速度平方積分IV2(Integral of the Squared Velocity)，該參數(shù)與地震早期輻射能量緊密相關，同時它也可用于預測地震動[21]。IV2的計算公式如式(4)，v為垂直向速度記錄，tc為P波到時點，Δt為積分計算的時間窗長。

(4)

(6)破壞烈度DI(Destructive Intensity)：通過地震動豎向加速度a和速度v的內積計算，該參數(shù)可估計地震的危險程度。DI的計算方法如式(5)，由式(5)可見，DI值表示的是“功率”的概念，DI取加速度和速度的內積絕對值的對數(shù)。

DI=lg|a·v|

(5)

(7)累積絕對速度CAV(Cumulative Absolute Velocity)：由地震動加速度的絕對值對時間間隔積分得到的[28]，該參數(shù)考慮了地震動的幅值與持時特征，CAV的由式(6)計算:

(6)

(8)阿里亞斯烈度Ia(Arias Intensity)由Arias[29]提出,該參數(shù)同時包含了地震動的幅值、頻率及持時特征，計算方法如下：

(7)

2.3 支持向量機

支持向量機是人工智能的機器學習方法中一種基于統(tǒng)計學習理論的監(jiān)督學習方法，能非常出色解決多變量回歸問題，其原理為將輸入數(shù)據(jù)映射到高維特征空間,將非線性問題變化為線性問題，通過變換后的線性問題的方法求解原來的非線性問題。

對于給定訓練數(shù)據(jù)D={(xi,yi)},i=1,2，…,mxi∈Rd,yi∈R，支持向量機找到一個回歸函數(shù)(式(8))，使得f(x)與y盡可能地接近。支持向量回歸模型與傳統(tǒng)回歸模型不同之處在于其容許f(x)與y之間最多有2ε誤差，即f(x)與y之間的誤差

f(x)=ωTφ(x)+b

(8)

(9)

式中：C為正則化常數(shù)，lε為ε—不敏感損失(ε-insensitive loss)函數(shù)

lε=

(10)

s.t.f(xi)-yi≤ε+ξi

(11)

(12)

(13)

上述過程需要滿足KKT(Karush-Kuhn-Tucker)條件，用核函數(shù)K(xi,xj)替代[φ(xi)·φ(xj)]，即可得到預測函數(shù)模型：

(14)

3 結 果

本文在研究中遵循以下流程：

(1)以0.1～10 Hz帶通濾波后三分向矢量合成地震動峰值PGA與PGV為預測目標參數(shù)，且計算時間窗均為P波觸發(fā)后3 s；

(2)利用訓練數(shù)據(jù)集，以加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里葉譜幅值AMmax、速度平方積分IV2、破壞烈度DI、累積絕對速度CAV、阿里亞斯烈度Ia這8種特征參數(shù)取其所有組合做為輸入?yún)?shù)，訓練后建立預測模型，并選擇其中誤差標準差最小的參數(shù)組合做為本文的支持向量機預測模型SVM-PGA與SVM-PGV；

(3)利用訓練數(shù)據(jù)集，訓練后建立Pd預測模型Pd-PGA與Pd-PGV，形式分別為如下線性形式：

lg(PGA)=k1·lg(Pd)+b1

(15)

lg(PGV)=k2·lg(Pd)+b2

(16)

式中，k1、b1、k2、b2為擬合系數(shù)。

(4)利用測試數(shù)據(jù)集，對比支持向量機預測模型與Pd預測模型的預測結果；

(5)隨機選取兩次地震事件(不在訓練數(shù)據(jù)集與測試數(shù)據(jù)集范圍內)，對比支持向量機預測模型與Pd預測模型的預測結果。

3.1 訓練數(shù)據(jù)集與測試數(shù)據(jù)集

利用本文選取的日本K-net強震觀測臺網2007年10月至2019年8月的強震數(shù)據(jù)，按照70%數(shù)據(jù)為訓練集進行模型訓練，并在剩余30%數(shù)據(jù)上進行模型性能測試。

基于python3.5環(huán)境搭建支持向量機模型，選用高斯核函數(shù)，停止訓練的誤差精度取0.001，選取正則化參數(shù)C和核函數(shù)參數(shù)gamma并采用網格搜索來選取最優(yōu)組合，其余參數(shù)選擇默認值。

由于訓練數(shù)據(jù)集是在對數(shù)坐標下得出預測模型，因此將實測值對數(shù)與預測值對數(shù)之差定義為誤差xi，如式(17)，并計算誤差標準差σ，如式(18)：

xi=lg(實測值)-lg(預測值)

(17)

(18)

表1 不同參數(shù)組合模型在訓練數(shù)據(jù)集上預測PGA的誤差標準差

圖3 不同參數(shù)組合在訓練數(shù)據(jù)集上預測PGA的誤差標準差

依據(jù)訓練數(shù)據(jù)集建立的Pd預測模型如下：

lg(PGA)=0.567lg(Pd)+2.238

(19)

lg(PGV)=0.740lg(Pd)+1.213

(20)

表2列出了基于訓練數(shù)據(jù)集的支持向量機預測模型SVM-PGA與SVM-PGV、以及Pd預測模型Pd-PGA與Pd-PGV的誤差標準差σ計算結果。支持向量機預測模型SVM-PGA與SVM-PGV的誤差標準差σ分別為0.240與0.254，Pd預測模型Pd-PGA與Pd-PGV的誤差標準差σ分別為0.320與0.309。在訓練數(shù)據(jù)集上，支持向量機預測模型SVM-PGA與SVM-PGV的誤差標準差均分別小于Pd預測模型Pd-PGA與Pd-PGV。

表2 支持向量機預測模型與Pd預測模型在訓練數(shù)據(jù)集上誤差標準差

圖4為支持向量機預測模型SVM-PGA與SVM-PGV的誤差概率密度曲線、以及Pd預測模型Pd-PGA與Pd-PGV的誤差概率密度曲線。可以看到，這四種預測模型的誤差概率密度曲線均呈正態(tài)分布。圖4中可以得到，支持向量機預測模型SVM-PGA與SVM-PGV的誤差概率密度曲線形狀比Pd預測模型Pd-PGA與Pd-PGV的誤差概率密度曲線形狀更高、更窄，這表明支持向量機預測模型SVM-PGA與SVM-PGV比Pd預測模型Pd-PGA與Pd-PGV的誤差更小、誤差更趨近于零。

(a)PGA

基于訓練數(shù)據(jù)集得到預測模型后，使用支持向量機預測模型SVM-PGA與SVM-PGV、Pd預測模型Pd-PGA與Pd-PGV在測試數(shù)據(jù)集上進行測試。圖5表示了預測值與實測值的分布情況，1∶1的實線表示預測值等于實測值，虛線表示訓練模型得到的1倍標準差。從圖5(a)和圖5(c)可以得到，支持向量機預測模型SVM-PGA與SVM-PGV的誤差σ分別為0.240與0.257，與訓練數(shù)據(jù)集得到的結果幾乎一致。從圖5(b)和圖5(d)可以得到，Pd預測模型Pd-PGA與Pd-PGV的誤差標準差σ分別為0.313與0.304，小于訓練數(shù)據(jù)集得到的結果。支持向量機預測模型SVM-PGA與SVM-PGV的誤差標準差均分別小于Pd預測模型Pd-PGA與Pd-PGV，這一點與訓練數(shù)據(jù)集得到的結果一致。

同時，圖5的結果也顯示，Pd預測模型在整體震級范圍內小數(shù)值區(qū)域出現(xiàn)高估現(xiàn)象、在高數(shù)值區(qū)域出現(xiàn)低估現(xiàn)象，即“小值高估、大值低估”，支持向量機預測模型在6.5級以下地震時整體趨近1∶1線性比例關系、6.5級以上地震出現(xiàn)了輕微的低估現(xiàn)象。

(a)SVM-PGA

圖6顯示了支持向量機預測模型誤差以及Pd預測模型誤差隨震中距變化的分布。

圖6(a)與圖6(c)分別顯示了支持向量機預測模型SVM-PGA與SVM-PGV的預測誤差隨震中距變化的趨勢，可以看到，支持向量機預測模型SVM-PGA與SVM-PGV的預測誤差在0～200 km距離區(qū)間均勻的分布在零誤差(虛線)兩側，這表明支持向量機預測模型的預測結果不受到震中距變化的影響，且整體誤差位于(-0.6，0.8)區(qū)間，這表明依據(jù)本文提出的支持向量機預測模型計算儀器地震烈度，最大誤差為0.8烈度單位。

(a)SVM-PGA

圖6(b)與圖6(d)分別顯示了Pd預測模型Pd-PGA與Pd-PGV的預測誤差隨震中距變化的趨勢，可以發(fā)現(xiàn)，Pd預測模型Pd-PGA的預測誤差明顯隨震中距逐漸增大而減小，Pd預測模型Pd-PGV的預測誤差也有隨震中距逐漸增大而減小的趨勢，且Pd預測模型Pd-PGA與Pd-PGV的整體誤差位于(-1.0，1.0)區(qū)間，這表明依據(jù)本文提出的支持向量機預測模型計算儀器烈度，最大誤差為1.0烈度單位。

3.2 震例分析

為了解支持向量機預測模型在具體地震事件上預測PGA、PGV的性能，隨機選取兩次不在訓練數(shù)據(jù)集和測試數(shù)據(jù)集范圍內的地震事件進行震例分析，這兩次地震的震中位置及臺站分布如圖7所示，其中，圖7(a)為2002年9月16日Mj5.3級地震，共計60組強震動記錄，圖7(b)為2007年3月25日Mj6.9級地震，共計93組強震動記錄。

圖7 震例分析的地震震中(五角星)及K-net臺站(三角形)分布

圖8為支持向量機預測模型和Pd預測模型在2002年9月16日Mj5.3級地震的預測值與實測值的分布，1∶1的實線表示預測值等于實測值，虛線表示基于訓練數(shù)據(jù)集得到的1倍誤差標準差。圖8(a)與圖8(c)分別為SVM-PGA與SVM-PGA預測模型的預測結果，可以看到預測值與實測值趨近1∶1線性比例關系，絕大部分數(shù)據(jù)分布在1倍誤差標準差范圍內。圖8(b)和圖8(d)分別為Pd-PGA、Pd-PGV預測模型的預測結果，可以看到Pd預測模型在小數(shù)值區(qū)域的預測值出現(xiàn)了明顯高估，在大數(shù)值區(qū)域的預測值出現(xiàn)輕微低估。相比Pd預測模型，支持向量機預測模型在2002年9月16日Mj5.3級地震上沒有明顯的“小值高估、大值低估”現(xiàn)象，預測值更趨近于實測值。

(a)SVM-PGA

圖9為支持向量機預測模型和Pd預測模型在2007年3月25日Mj6.9級地震的預測值與實測值的分布，1∶1的黑色實線表示預測值等于實測值，藍色虛線表示基于訓練數(shù)據(jù)集得到的1倍誤差標準差。圖9(a)和圖9(c)分別為SVM-PGA與SVM-PGA預測模型的預測結果，可以看到預測值與實測值趨近1∶1線性比例關系，雖然大部分數(shù)據(jù)在1倍誤差標準差范圍內，但整體出現(xiàn)了輕微的低估現(xiàn)象。圖9(b)和圖9(d)分別為Pd-PGA、Pd-PGV預測模型的預測結 果，可以看到Pd預測模型預測PGA的結果在小數(shù)值區(qū)域出現(xiàn)了明顯高估、大數(shù)值區(qū)域出現(xiàn)了明顯低估，且PGV的預測結果整體出現(xiàn)了低估。雖然支持向量機預測模型也出現(xiàn)了整體輕微低估的現(xiàn)象，但是這是由于3 s的預測時間窗與大地震事件較長的震源破裂過程持續(xù)時間不匹配所造成的[31]，相比Pd預測模型，支持向量機預測模型在2007年3月25日Mj6.9級地震沒有明顯的“小值低估”現(xiàn)象，且能夠明顯改善“大值低估”現(xiàn)象。

(a)SVM-PGA

4 結 論

本文以0.1～10 Hz帶通濾波后三分向矢量合成地震動峰值PGA與PGV為預測目標參數(shù)，利用日本K-net強震臺網2007年10月-2019年8月期間P波觸發(fā)后3 s數(shù)據(jù)，基于加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里葉譜幅值AMmax、速度平方積分IV2、破壞烈度DI、累積絕對速度CAV、阿里亞斯烈度Ia這8種特征參數(shù)構建支持向量機預測模型。得到的結論如下：

(1)無論是訓練數(shù)據(jù)集、測試數(shù)據(jù)集還是兩次隨機震例數(shù)據(jù)集，支持向量機預測模型SVM-PGA與SVM-PGV的預測值的整體誤差標準差均小于常規(guī)的Pd預測模型Pd-PGA與Pd-PGV；

(2)對比Pd預測模型Pd-PGA與Pd-PGV，支持向量機預測模型SVM-PGA與SVM-PGV的預測值與實測值更趨近于1∶1的線性比例關系，且預測值的誤差不隨震中距的變化而變化；

(3)相比較Pd預測模型Pd-PGA與Pd-PGV整體預測值明顯的“小值高估、大值低估”現(xiàn)象，支持向量機預測模型SVM-PGA與SVM-PGV的預測值在6.5級以下地震沒有明顯的“小值高估”現(xiàn)象，在6.5級以上地震中由于預測時間窗與震源破裂時間不匹配的原因而出現(xiàn)輕微低估。

綜上所述，本文構建的支持向量機預測模型可用于現(xiàn)地地震預警地震動峰值、即儀器地震烈度的預測。

致謝

日本防災科學技術研究所(NIED)為本研究提供了數(shù)據(jù)支持，所有數(shù)據(jù)均從K-net網站下載(網址:http://www.kyosh-in.bosai.go.jp/),文中圖件使用通用制圖工具GMT(Genetic Mapping Tools)繪制。