高速鐵路地震預警P波與S波復合自動快速識別的理論方法與應用

王子珺，趙伯明

(北京交通大學 土木建筑工程學院，北京　100044)

地震預警系指在破壞性強的地震抵達目標區域前的有限時間(幾秒至數十秒)內，利用專用監測預警系統發出預警信息，進而啟動防御災害的緊急處置措施。作為近年來得到快速發展的抗震減災的重要手段之一，日、美、意等地震多發國家已經在地震預警方面開展了大量的理論與應用研究[1-5]。尤其是對于可以控制其運行的生命線工程和高危行業，例如高鐵、核電站、燃氣等領域，地震預警系統(EEWS)已在多次地震中被證明能夠有效地減輕地震災害。

高速鐵路地震預警主要包括2個重要過程：預警報警和緊急控制。預警報警是在破壞性地震動尚未到達沿線之前，基于預警判斷過程發出報警信息；緊急控制是指在接收到報警信息后緊急啟動控車措施，使高速運行中的列車實現減速、停車等目的。這種預警—控車機制，可以防止由于地震造成的列車脫軌或傾覆等災難性事故的發生，減少次生災害和人員傷亡，保障列車在運營過程中的安全。

根據《高速鐵路自然災害及異物侵限監測系統總體技術方案(暫行)》(鐵科技[2013] 35號)，高鐵采用基于沿線路布置的原位線型地震監測自備網絡，地震監測點設置于牽引變電所、AT所、分區所內，布設間距平均為20 km左右，今后新建高鐵也將同期建設此類沿線地震監測網絡。線型網絡依存于地震波的入射方位，其觸發監測點的初始波到時具有較大的分散性，同時原位預警的方式也要求將預警信息進行快速實時處理，因此利用單臺站的P波與S波初始段數據進行地震預警是高速鐵路所必須采取的基本方式。事實上，日本是最早開展高速鐵路地震預警研究的國家，在世界上率先提出了利用單臺站記錄快速進行地震預警的模式[2]。高鐵預警系統是高鐵行車安全的相關系統，一旦發生強震漏報誤報將會導致高鐵重大災害事件。中國高鐵運營速度高，采用以橋代路的里程所占比例較大，所處地震環境復雜，地震時的危險性更高，急需開發具有實時自動、精準快速的適合中國高鐵特點的地震預警理論與方法。綜上所述，基于單臺站的初始先達記錄，對地震波進行快速自動高精度地識別，是實現高鐵單臺站預警目標的首要環節和必要條件。目前，關于P波與S波的識別方法因其應用目的不同而具有不同的特點，但是均不能很好地同時達到自動、快速和準確這3項要求，難以滿足中國高鐵所需要的單臺站預警模式。因此，針對我國高鐵地震預警的實際需求，研究提出1套適合于單臺站快速準確地震預警的三步驟P波與S波復合自動快速識別理論，并選取記錄的2008年汶川Ms8.0級地震的主震及余震的地面運動加速度(以下簡稱為加速度)，對提出的理論方法進行綜合驗證。

1　地震數據與數據預處理

中國數字強震動臺網在2008年汶川Ms8.0級地震中獲得了大量的強震觀測記錄，觀測臺網所使用的儀器均為安裝于自由地表的數字強震儀，采樣頻率為200 Hz，測量范圍為±2g。除去420個臺站獲得的主震記錄外，截至2008年9月30日12時，臺網共記錄余震383次；其中Ms4.0級及以上余震611次，Ms5.0級及以上余震56次，Ms6.0級及以上余震8次[6]。由于發震地區位于青藏高原向平原過渡的區域，具有復雜豐富的地質地貌特征，因此廣泛分布于該地區的地震記錄為開展預警方法的研究以及開發更具有適用性及穩健性的P波和S波識別方法提供了重要的數據保障。

近斷層的淺源大震級地震往往造成嚴重的地震災害[7]，為此本文選取了汶川地震的主震以及50個余震記錄作為研究樣本，地震記錄的篩選條件為，震級大于Ms5.0且震源距小于150 km。除去其中沒有完整得到初始先達P波的記錄，最終使用并解析了1 548條滿足要求的加速度，相關地震事件的震源位置以及觀測臺站的空間分布如圖1所示。

圖1　地震事件的震源位置以及觀測臺站的空間分布

數據預處理工作如下：首先將原始記錄減去記錄全長的平均值，完成基線調零；然后對調零后的加速度進行0.5～30 Hz的4階Butterworth數字帶通濾波。該頻帶適用于本研究區域內的地震記錄，對于其他地區，其帶通的取值范圍則需要根據數據的特性進行調整。

2　P波與S波復合自動快速識別理論

提出的P波與S波復合自動快速識別理論由三步驟構成。首先通過極化分析方法將P波與S波進行分離；然后研究利用改進后的短時窗信號平均值(STA)與長時窗信號平均值(LTA)之比的方法初步確定P波與S波的到時范圍；最后根據初步到時范圍，通過改良的基于高階統計量的AIC算法推導P波與S波的精確到時。

2.1　P波與S波的分離

P波與S波的分離主要基于Jurkevics提出的極化分析法[8]，首先通過1個移動的固定時間窗，建立預處理后的三分量加速度的協方差矩陣

(1)

式中：var(*)和cov(*, *)分別表示三分量加速度的自協方差和協方差；z表示豎直分量，n表示南北分量，e表示東西分量。

協方差矩陣C是半正定矩陣，其本質為橢球體方程的二次式系數矩陣。假設矩陣C的3個特征值分別是λ1，λ2和λ3(λ1≥λ2≥λ3)，其對應的特征向量為u=(u1，u2，u3)，表示橢球體的3個軸方向。則根據每一時間窗內的協方差矩陣可以計算得到相應的2組偏振特性參數，基于該偏振特性參數建立本文所使用的震相濾波器。

第1組偏振特性參數為偏振度，定義為

(2)

偏振度Rec的值介于0～1之間，由于純體波具有完全線性偏振的特點，即偏振度等于1；而背景噪音則不具有相關特點，即偏振度為0；則該組參數可以用于地震波與噪音的分離。

第2組偏振特性參數是表面垂直入射角，可表示為

φinc=arccos|u11|

(3)

式中：u11為特征向量u1的垂直方向余弦。

表面垂直入射角φinc的值介于0°～90°之間，對于區域性地震，由于P波沿著近乎垂直的入射線發生偏振，因此，其入射角的余弦值接近于1；而S波的偏振方向趨于水平，則其入射角的正弦值接近于1。

基于上述這2組偏振特性參數，可以建立2個極化濾波器來分離P波和S波。采用Rosenberger[9]提出的濾波方法，將每一時間步的三分量加速度分別乘以相應的濾波器方程，從而獲得極化后的地震波。

選取2008年5月12日14時43分汶川地震余震的理縣木卡臺站的原始加速度三分量記錄，如圖2所示，該余震震級為Ms6.3，震中距約為56 km；采用0.2，2.0，5.0 s不同長度時間窗進行P波和S波的濾波，如圖3所示。由圖3可知：采用較小的時間窗長度(0.2 s)不能很好地將地震信號與噪音信號分離；而隨著時間窗長度的不斷增大，濾波方程趨于平滑且對信號的變化更加靈敏，從而易于對不同類型的地震波進行判斷識別；這是因為移動時間窗的長度直接關系到偏振度和入射角的計算，進而會對濾波器的濾波效果產生影響。根據Baillard等[10]的研究成果，在進行偏振分析時，時間窗長度應大于地震波最短周期的3倍，同時要小于P波與S波之間的時間差；通過分析發現，采用1個長度為2.0 s的移動時間窗即可達到理想的效果。通過極化濾波后的南北、東西、豎直三分量加速度如圖4所示。

圖2　原始加速度三分量

圖3　不同長度移動時窗下P波與S波的濾波結果

圖4　三分量加速度的P波與S波分離結果

由圖4可知：利用建立的極化濾波器能夠很好地完成P波與S波的分離，從而有效地保證震相的初至識別。

2.2　識別算法

首先對傳統的STA/LTA方法進行優化改進研究，根據得到的極化后的P波與S波，利用該改進方法初步確定P波與S波的到時范圍，在該范圍內，利用基于高階統計量的AIC算法分別推導P波與S波的精確到時，具體過程如下。

STA/LTA算法的基本原理是利用STA與LTA比值的變化來反應信號的變化。由于短時間窗反映的主要是信號幅值的瞬時變化，而長時間窗反映的是待檢信號背景噪音的平均幅值；當震相到達時，STA會比LTA的變化頻度高，相應的比值會有一個明顯的增加，當比值超過某一觸發閾值時，則判定有震相到達。為了更加明顯地反映信號到達時頻率及幅值等特征量的變化，通常采用Allen[11]提出的特征函數CF(characteristic function)代替原始信號進行STA/LTA的計算。該特征函數為

CF(i)=x(i)2+W(x(i)-x(i-1))2

(4)

式中：x(i)(i=1，2，…，L)為ti時刻的加速度，L為隨機變量x的樣本數；W為隨著采樣率和背景噪音性質而變化的加權常數，根據經驗取值為3。

該特征函數CF考慮了實際地震的加速度以及向前一階差分的影響，能夠靈敏地反映信號到時的頻率及幅值變化特征。雖然該特征函數能夠有效降低震相識別的誤判率，但是由于觸發閾值是根據經驗設定的固定值，在低信噪比的情況下可能會造成震相的漏判或誤判。所以為了提高P波與S波的初始到時的識別率，提出導入1個調整系數δ，對傳統的STA與LTA的比值進行優化改進，該調整系數的定義為

(5)

其中，

(6)

由于標準差可以衡量某一數據序列的數值離散程度，在地震事件發生前，系數δ趨近于1；當P波或S波的信號到達時，系數δ會有一個較大幅度的增長。因此這個改進后的STA/LTA方法可以明顯地提高對信噪比低或初動信號不明顯的記錄的識別率，從而大幅度提高識別的精度和速度。

傳統STA/LTA方法與改進后STA/LTA方法的識別效果如圖5所示，可以看出，改進后的方法對震相初至識別有非常明顯的提升效果。利用該方法可以確定震相的到時范圍，其精確的到時判定將會在下一階段完成。P波的拾取使用豎直向加速度分量，S波的拾取使用水平向加速度分量，并設初始到時的時間為T初。

圖5傳統STA/LTA方法與改進后STA/LTA方法的識別效果比較

移動時間窗的長度由采樣頻率和特征信號的周期所決定。對于STA，若時間窗的長度選擇過大，則會導致特征函數過于平滑而無法檢測到信號的瞬時變化，從而產生漏觸發；若時間窗的長度選擇過小，則會對短周期的干擾信號過于敏感從而造成誤觸發。對于LTA，其窗長取值應能反映背景噪音的平均水平。基于對大量記錄的試算，本文取：短時窗tSTA=0.5 s，長時窗tLTA=3～5 s，用戶定義參數Δ=35；P波與S波拾取的觸發閾值分別設為5和10。

為了精確確定震相的到達時刻，在下一階段中采用基于改良的高階統計量的AIC函數對地震波進一步的識別。根據Maeda[12]，計算地震數據x(i)的AIC函數，其最小值對應的時間即為地震震相的精確到時，AIC函數為

AIC(k)=klog{var(x[1，k])}+(L-K-

1)log{var(x[k+1，L])}

(7)

式中：k為滑動變量，k∈[1，L-1]。

基于數學分析理論，高階統計量與二階統計量相比具有適于檢測和描述系統的非線性的特點。在低信噪比的條件下，可有效抑制背景噪音的影響，易于對較弱的地震信號進行檢測和識別，因此本文導入峰度函數代替二階方差函數，并應用于AIC函數以確定震相的精確到時。峰度函數K為樣本的四階中心矩，即

(8)

在震相初始到時T初的基礎上，通過確定1個時間窗[T初-0.5，T初+0.5]，應用峰度-AIC函數對P波與S波的到時進行精確識別。選擇0.5 s既可以保證將觸發時刻控制在該時間窗內，又能夠滿足地震預警所需的實時性和時效性要求。

圖6表示應用峰度-AIC函數對震相進行二次拾取的結果。由圖6可以看出，利用第1階段改進后的STA/LTA方法初步確定P波與S波的到時后，再求取峰度-AIC函數的最小值，就能夠準確地判斷震相的實際到時。

圖6　峰度-AIC函數的最小值對應震相的精確到時

由此可見，由于改進后的SLA/LTA方法具有適應性強、拾取率高的特點，在使用該方法進行震相初始識別的基礎上，通過改良的基于高階統計量的AIC函數進一步保證了識別的穩定性和精確性。因此，本文命名上述由3個步驟構成的識別方法為P波與S波復合自動快速識別方法。

3　汶川地震數據的應用與分析

將提出的P波與S波復合自動快速識別方法應用于記錄的2008年5月12日14時43分汶川Ms8.0級地震的主震以及其后50個余震的1548條地震的加速度，識別震相的精確到時，并與人機交互處理結果和改進后的STA/LTA方法處理結果進行對比，驗證提出的復合方法的準確性和有效性。在判別過程中，P波的拾取針對豎直方向的加速度分量，S波的拾取針對南北以及東西2個水平方向的加速度分量。以2.1節使用的余震為例，選取茂新南新、理縣木卡和理縣沙壩3個臺站所獲得的地震加速度。其P波的到時識別過程如圖7所示，S波的到時識別過程如圖8和圖9所示。

圖7不同臺站的P波的到時識別過程(采用豎直方向加速度分量)

圖8不同臺站S波的到時識別過程(采用南北方向水平加速度分量)

為進一步研究P波與S波復合自動快速識別方法的準確性，以人機交互結果為標準，對使用的1 548個三分量加速度的自動識別結果進行誤差計算，其偏差分布如圖10所示。P波到時的自動識別的結果為：拾取正確率達到91%，其中89%的誤差在0.1 s以內，平均偏差為-0.021 s，標準差為0.068 s。S波到時的自動識別結果為：按照識別誤差小于0.5 s的精度要求，識別有效率高達85%，其中有92%的誤差在0.2 s以內，平均偏差為-0.025 s，標準差為0.169 s。

圖9不同臺站S波的到時識別過程(采用東西方向水平加速度分量)

圖10P波與S波的自動識別與人機交互識別結果的偏差分布

P波與S波自動識別的平均偏差基本一致，但P波識別的標準差以及識別率優于S波，因為P波屬于縱向壓縮波，其傳播過程比橫向S波相對單純，而S波的產生及傳播過程比較復雜，受傳播路徑的影響較大，對其識別相對有較大難度。另外，本文對P波的拾取只針對豎直方向的加速度，而S波的拾取則使用了水平2個方向的記錄，S波的樣本數量是P波的2倍，因此也會造成S波識別率要低于P波的結果。考慮到每個觀測點都會受到復雜

的震源、傳播路徑以及場地效應的影響，根據既有研究無法實現百分百的識別率。

此外本文也對改進后的STA/LTA方法與人機交互的結果進行了比較：對于P波到時的識別，兩者之間偏差的平均值為-0.217±0.326 s；對于S波到時的識別，兩者偏差的平均值為-0.191±0.379 s。該結果表明，通過對傳統的STA/LTA方法進行改進，大幅度提高了識別精度和速度，但由三步驟構成的P波與S波復合自動快速識別方法的效果更優。

因此，綜合以上結果，本文提出的結合改進后的STA/LTA方法以及改良的基于高階統計量的AIC函數得到的復合方法可以顯著提高P波和S波的識別效果。

4　結　論

(1)高速鐵路地震預警的研究就是保障預警過程的實時性、準確性和高效性。針對我國高鐵的特點，以及高鐵地震預警對時間和精度的要求，基于沿線路排列的原位地震監測網絡所必需的單臺站預警模式，本文提出了1套適合于高鐵單臺站快速準確地震預警的三步驟P波與S波復合自動快速識別理論。

(2)為了校驗提案方法的泛用性，整理解析了具有代表性的地震——2008年汶川Ms8.0級地震的主震及余震記錄1 548條，通過將自動識別方法與人機交互處理的結果進行對比，驗證了該方法可以準確、有效地應用于具有復雜的震源、傳播路徑以及場地效應的地震進行快速預警。比之其他既有方法具有實時性、精度高、速度快的特點，能夠滿足高鐵單臺站P波或S波地震預警的嚴格要求。由于強震觀測點條件的多樣性，對P波與S波的識別率造成了一定影響，如果使用針對高鐵地震預警建設的臺站所獲得的地震記錄，那么成功率會進一步提高。

[1]NAKAMURA Y，SAITAJ，SATO T. On an Earthquake Early Warning System (EEW) and Its Applications [J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：127-136.

[2]NAKAMURA Y. On the Urgent Earthquake Detection and Alarm System (UrEDAS)[C]//Proceedings of 9th World Conference on Earthquake Engineering. Tokyo-Kyoto，Japan：International Conference for Earthquake Engineering (IAEE)& the 9th WCEE Organizing Committee，1988，Ⅶ：673-678.

[3]BROWN HM，ALLEN RM，HELLWEG M，et al. Development of the ElarmS Methodology for Earthquake Early Warning: Realtime Application in California and Offline Testing in Japan[J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：188-200.

[4]SATRIANO C，ELIA L，MARTINO C，et al. PRESTo, the Earthquake Early Warning System for Southern Italy: Concepts, Capabilities and Future Perspectives [J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：137-153.

[5]IANNACCONE G，ZOLLO A，ELIA，et al. A Prototype System for Earthquake Early-Warning and Alert Management in Southern Italy [J]. Bulletin of Earthquake Engineering，2010，8(5)：1105-1129.

[6]國家強震動臺網中心.2008年汶川地震數據集[DB/OL].[2015-11-01].ftp://www.csmnc.net.

[7]WU Y M，KANAMORI H. Experiment on an Onsite Early Warning Method for the Taiwan Early Warning System[J]. Bulletin of the Seismological Society of America，2005，95(1)：347-353.

[8]JURKEVICS A. Polarisation Analysis of Three-Component Array Data [J]. Bulletin of the Seismological Society of America，1988，78(5)：1725-1743.

[9]ROSENBERGER A. Real-Time Ground Motion Analysis: Distinguishing P and S Arrivals in a Noisy Environment[J]. Bulletin of the Seismological Society of America，2010，100(3)：1252-1262.

[10]BAILLARD C，CRAWFORD WC，BALLU V，et al. An Automatic Kurtosis-Based P- and S-Phase Picker Designed for Local Seismic Networks [J]. Bulletin of the Seismological Society of America，2014，104(1)：394-409.

[11]ALLEN R V. Automatic Phase Pickers: Their Present Use and Future Prospects [J]. Bulletin of the Seismological Society of America，1982，72(6)：S225-S242.

[12]MAEDA N. A Method for Reading and Checking Phase Times in Auto Processing System of Seismic Wave Data [J]. Zisin=Jishin，1985，38：365-379.