峰度和AIC方法在區域地震事件和直達P波初動自動識別方面的應用*

趙大鵬，劉希強，李 紅，周彥文

(1.中國地震局蘭州地震研究所，甘肅蘭州730000;2.山東地震局，山東濟南250014)

峰度和AIC方法在區域地震事件和直達P波初動自動識別方面的應用*

趙大鵬1，劉希強2，李 紅2，周彥文2

(1.中國地震局蘭州地震研究所，甘肅蘭州730000;2.山東地震局，山東濟南250014)

提出了一種基于直達P波信號的峰度和Kurtosis-AIC方法進行區域地震事件實時檢測和直達P波初動精細識別的新方法，并應用于山東地震臺網記錄的地震波資料處理.結果表明:(1)應用峰度方法能夠有效識別出地震事件，可有效減低地震事件的錯誤報警率和漏報率;(2)與人工識別震相到時的結果相比，根據Kurtosis-AIC震相自動識別方法得到的震相到時的平均絕對值誤差為 (0.09±0.08)s。

峰度方法;Kurtosis-AIC方法;地震事件識別;震相識別

0 引言

日本是第一個建立地震預警系統的國家，日本預警系統在2011年3月11日日本發生9.0級特大地震時發揮了作用，東京地區民眾在地震S波最大震動到達之前收到地震報警信息，得以避險(陳會忠等，2011)。隨著經濟的快速發展和社會的不斷進步，社會公眾對地震速報提出了越來越高的要求。地震預警具有極其重要的社會效益，它不僅可以減輕人員傷亡和降低次生地震災害的發生，還可為震后緊急救援和搶修提供第一手的信息 (劉希強等，2009)。

快速確定地震參數是地震預警系統的關鍵技術環節之一，關系到預警時間的長短。現場地震預警系統利用P波進行預警可以獲得更長的預警時間，即在預警的目標區安裝地震計，利用P波比S波傳播速度快的原理，用P波到達后數秒的數據來確定地震參數，對S波還沒到達的地區通過電磁波發出預警，從而達到預警的目的 (周彥文等，2010)。目前國外已經涌現出大量關于地震事件自動實時檢測、P波初動拾取、震相類型識別方面的研究方法和成果 (Allen，1978，1982;Bear，Kradolfer，1987;Sleeman，Van，1999;Bai，Kennett，2000;Diehl et al，2009;Kuperkoch et al，2010)。國內學者提出了基于分形分維 (常旭，劉伊克，2002)、自回歸模型 (Zhang et al.，2003;王海軍等，2003;楊配新等，2004;王繼等，2006)、人工神經網絡法 (Wang，Teng，1995;張范民，李清河，1998;李慧婷，黃文輝，2000;王娟等，2004)、波形變化增長法 (朱元清等，2002)、小波變換的主成分分析法 (劉希強等，2000)、定尺度比方法 (劉希強等，2002)、地震事件自動檢測的EFGLP方法 (劉希強等，2009)、初至波震相自動精確識別的TOC-AIC方法 (劉希強等，2009)、基于幅值和頻率的P波識別 (毛燕等，2011)等來識別震相，取得了較好的效果。

Kuperkoch等 (2010)根據高階統計量 (HOS)中的峰度(kurtosis)和偏斜度 (skewness)與Akaike信息準則提出了一種地震事件和震相的自動檢測和識別算法，并將該算法與Allen(1978)、Bear和Kradolfer(1987)的檢測算法對自動獲得P波震相讀取值進行了比較研究。研究表明，高階統計量方法顯示出了更佳的識別效果，尤其峰度識別效果最佳，且分別使用偏斜度和峰度對3 000多個人工測量的P波震相的到時識別其平均偏差分別為(0.38±0.75)s和 (0.26±0.64)s。筆者使用峰度方法對山東地震臺網記錄的單臺垂向記錄進行了區域地震事件自動檢測和震相識別探索研究，同時應用AIC方法對峰度方法所確定的包含地震P波在內的時間窗內記錄進行自動分析，探索進一步提高P波震相初至識別的技術和方法。

1 地震波形資料的橢圓型帶通濾波方法

橢圓型帶通濾波器屬于無限脈沖響應濾波器，它在有限頻率范圍內存在傳輸零點和極點;其通帶和阻帶都具有等波紋特性，通帶、阻帶逼近特性良好;對于同樣的性能要求，該濾波器比巴特沃斯、切比雪夫這兩種濾波器所需用的階數都低，而且它的過渡帶比較窄。

圖1a給出了階數為7的橢圓形帶通濾波器，其基本參數如下:通帶波動系數為0.1，阻帶衰減系數為30、數據采樣率為100 p/s、通帶低端截止頻率5.5 Hz、通帶高端截止頻率7.5 Hz、帶通阻帶低端截止頻率為5 Hz、帶通阻帶高端截止頻率為8 Hz。從幅頻特性可以看出，橢圓型帶通濾波器具有通帶、阻帶逼近特性良好的特點。圖1b給出的是由1.5、6.5和10.5 Hz正弦波組成的合成信號，其中6.5 Hz的信號起始點在2 s。圖1c給出了根據上述濾波器進行數據濾波的結果，可以看出，濾波后結果與仿真信號相比，在振幅和相位兩方面非常一致。

圖1 橢圓形濾波器對合成模擬仿真信號的帶通濾波結果(a)橢圓形帶通濾波器的幅頻特征;(b)由3種頻率合成的仿真信號;(c)帶通濾波信號與指定頻段仿真信號的對比Fig.1 Synthesis analog simulation signal filtered by elliptical band-pass filter(a)Amplitude-frequency characteristics of the elliptical bandpass filter;(b)Analog simulation signal synthesized by sine wave of three different frequency;(c)Comparison between band-pass filter signal and simulation signal with designated frequency

在處理區域地震事件 (震中距小于120 km)的數據時，考慮到其主頻分布范圍，對橢圓型濾波器的設計參數如下:通帶低端截止頻率2 Hz、通帶高端截止頻率15 Hz、帶通阻帶低端截止頻率為1.5 Hz、帶通阻帶高端截止頻率為16 Hz。

2 區域地震事件觸發判斷的峰度方法

隨機變量X的數學期望表示為

式中，p(x)是隨機變量X的概率密度。

利用期望函數，隨機變量X的k階統計量mk可表示為

式中，當k=2時表示X的方差，m2表示X的取值與數學期望的偏離程度。若X比較集中，則m2較小;反之，若取值比較分散，則m2較大。因此，m2是刻畫X取值分散程度的一個量，它是衡量X取值分散程度的一個尺度。

峰度是指頻數分布的集中程度，也就是分布曲線的尖峭程度。隨機變量的峰度的表達式為

峰度K是四階中心動差m4除以標準差四次方。將四階動差m4除以標準差四次方，同樣也是為了得到相對數，以便進行不同頻數分配的峰度比較。當峰度K＞3時，表示頻數分布比正態分布更集中，分布呈尖峰狀態，平均數的代表性更大;K＜3時表示頻數分布比正態分布更分散，分布呈平坦峰，平均數的代表性較小。

3 初至波震相自動精確識別的Kurtosis-AIC方法

常用的震相自動識別技術是AR-AIC方法，其基本原理是假定地震波記錄可分為若干個局部穩態的過程，每個過程滿足自回歸模型，并認為新的地震波震相到達前后的信號是兩個不同的穩態過程 (Sleeman，Van，1999)。設長度為N的地震波記錄包含有背景信號和地震波信號，通過滑動變量k的不斷改變，得到長度分別為 ［1，k－M］和 ［k－M+1，N－M－k］數據段的最大似然函數的最大值，進而得到AIC函數:

式中，k表示滑動變量 (下同)，M表示AR擬合過程的階數，和表示最大似然函數的最大值，C2是常數。AR過程的階數需要通過實際記錄并經多次試驗計算得到，AR系數根據Yule-Walker方程確定 (Haykin，1996)。當所選數據中包含有背景噪聲和地震波信號時，AIC函數的最小值對應的時間就是地震震相初至。

Maeda(1985)提出了無需計算AR系數，直接根據地震圖記錄而得到AIC函數的方法，稱為VAR-AIC方法。AIC函數的最小值對應的時間就是地震震相初至。AIC函數表示為

式中，var(x［1，k］)和 var(x［k+1，N］)分別表示兩個時間段內數據的方差。

VAR-AIC方法屬于二階統計量方法。二階統計量 (如二階矩、相關函數、功率譜密度函數等)方法需要假設系統是最小相位系統，噪聲為高斯白噪聲，然而實際信號常常不滿足這個假設。高階統計量與二階統計量相比具有適于檢測和描述系統的非線性，在低信噪比條件下，可有效抑制加性高斯有色噪聲的影響，檢測和識別較弱的地震信號。用峰度進行精細震相識別的具體思路是用兩個時間段數據的峰度作為特征函數代替公式(5)中的方差，本文稱為Kurtosis-AIC方法。Kurtosis-AIC函數表示為

式中，L是特征函數CF的長度，k取值為從1到L。識別震相初至的具體思路是，以峰度方法自動確定的地震事件觸發判斷結果為參考，選擇一個包含有噪聲和地震信號的時間窗，Kurtosis-AIC函數的最小值對應時刻即為震相初至。

4 資料和結果分析

本文對50次地震的198個單臺垂直向記錄進行了分析與研究，涉及的臺站數量有53個 (圖2)。所選取單臺記錄的最大震中距為110.0 km、最小為6.0 km，震中距小于100 km的單臺記錄有190個，占總數的96%。50次地震中的最小震級為0.6級，最大震級為4.2級。

圖2 研究區地震和臺站分布圖Fig.2 Distribution of earthquakes and stations in study area

采用逐點滑動的方式對數據的峰度進行實時計算和分析，計算時間窗長為10 s，步長為1個采樣點。為了總結地震事件觸發閾值的分布規律，選擇大于噪聲段最大峰度值作為每次地震事件觸發的最低閾值，峰度超過閾值所對應的時間初定為P波震相初至。然后應用時間段的峰度數據進行Kurtosis-AIC方法分析，其最小值對應時刻即為精細識別的P波震相。圖3給出了梁山臺記錄的2009年6月15日22時27分范縣2.4級地震垂直向波形、事件觸發時刻及應用峰度和Kurtosis-AIC方法對初至震相到時的拾取結果。當選取峰度閾值等于3時，由峰度方法得到的地震事件觸發時刻為49.69 s，應用Kurtosis-AIC方法得到的P波震相初至結果為49.68 s，與人機交互拾取結果49.64 s相比，震相識別的誤差為0.04 s，識別結果非常理想。

圖4給出了數據濾波前后信噪比的比較及信噪比隨震中距的變化。從圖4可以看出，由于對信號進行了橢圓型帶通濾波，有效提高了信噪比。圖4a中數據濾波前信噪比的變化范圍為－3.8～47.6 dB，均值為8.7 dB，標準差為9.4 dB，濾波后信噪比的變化范圍為5.1～49.8 dB，均值為22.0 dB，標準差為9.0 dB，可見對觀測數據進行濾波后信噪比得到了明顯提高。圖4b表明對發生在山東地區2級左右小地震進行濾波處理后的信噪比隨震中距的變化并未出現明顯衰減的特點。

圖3 梁山臺記錄的2009年6月15日22時27分范縣2.4級地震垂直向波形及事件觸發參數變化(a)垂直向記錄;(b)(a)中虛框內記錄;(c)對 (b)中波形的濾波結果;(d)數據的峰度曲線;(e)數據的Kurtosis-AIC變化曲線Fig.3 Vertical component waveform of Fanxian ML2.4 earthquake at 22:27 on Jun.15，2009 recorded by Liangshan Station and variation of event-trigger parameters(a)Vertical component waveform;(b)Waveform within the dotted frame in graph(a);(c)Filter result of the waveform in graph(b);(d)Kurtosis variation curve of data;(e)Kurtosis-AIC variation curve of data

圖4 數據濾波前后信噪比的比較 (a)及信噪比隨震中距的變化 (b)Fig.4 Signal-to-noise ratio comparison of data before and after the filtering(a)and signal-tonoise ratio change with epicenter distance(b)

圖5 用于地震時間觸發判斷的峰度最低閾值 (a)和峰度最大值分布圖 (b)Fig.5 Minimum threshold of Kurtosis used to judge the seismic time trigger(a)and distribution of the maximum Kurtosis value(b)

圖5給出了地震事件觸發判斷的峰度最低閾值和最大峰度值分布圖。從圖5a中可以看出，在對數據樣本處理時，保證地震事件峰度觸發閾值變化范圍為2.5～8，均值為3.35，標準差為0.70。與之相對應，圖5b中的峰度值的變化范圍為4.2～691.9，其中位于4.2～8.0區間有3次，均值為137.0，標準差為118.5。為了保證地震事件既不漏報又不誤報，需要在積累的大量樣本中選擇1個經驗閾值。如果本文選定峰度閾值為8，研究樣本中有3個單臺記錄出現漏報，這3次漏報地震的震級都為1.3級，漏報率為1.5%，誤報率為0。上述出現漏報資料的信噪比都不高，分別是3.9 dB、1.08 dB和－1.7 dB。如果設置峰度閾值為8，并且地方震的單臺垂向記錄數據滿足信噪比大于4 dB，則較大地震不會出現漏報和誤報。

圖6給出了Kurtosis-AIC震相自動識別結果與人機交互識別結果對比的誤差分布及誤差分布與觀測數據信噪比的關系。由圖6a可以看出，震相識別誤差范圍為－0.19～0.5 s，平均絕對誤差為0.09 s，標準差為0.08 s。198個單臺垂直向記錄的處理結果中有127個誤差在0.1 s內，占總數的64.14%;有174個誤差在0.17 s內，占88%;有184個誤差在0.2 s內，占93%。從圖6a和圖6b對比可以看出，震相識別誤差較大的樣本相對應的垂向信噪比多數比較低，反應了在人工分析中同樣面臨的當P波信號較弱時震相初至較難識別這一難題。但總體上Kurtosis-AIC方法對較弱的P波地震信號，也具有較好的分辨能力，達到了比較高的精度。

圖6 Kurtosis-AIC震相自動識別結果與人機交互識別結果對比的誤差分布 (a)及誤差分布與觀測數據信噪比的關系 (b)Fig.6 Error distribution of comparision between the result of Kurtosis-AIC automatic recognition of seismic phase and that of man-machine interaction recognition(a)and the relationship between error distribution and signal-to-noise ratio of data(b)

5 結論

本文利用山東地震臺網記錄的單臺垂向記錄和峰度方法進行了區域地震事件自動檢測和震相識別方法的探索研究，為了減小噪聲干擾又不失地方震記錄信號，本文設計了保持信號的振幅和相位畸變最小的地震波形資料的橢圓型帶通濾波方法。通過研究取得了以下結果:

(1)橢圓型濾波器提高了含噪聲P波信號的信噪比。

(2)提出了一種基于直達P波信號進行區域地震事件實時檢測的峰度方法，并用該方法研究得到了山東地震臺網記錄判斷地震事件觸發的經驗閾值。表明選擇合適的峰度經驗閾值可同時有效降低地震事件的漏報率和誤報率。

(3)提出了峰度與AIC方法相結合的精細震相識別新方法對198個單臺垂向記錄的分析研究表明，與人工識別震相到時的結果相比，根據Kurtosis-AIC震相自動識別方法得到的震相到時的平均絕對值誤差為0.09±0.08 s，精度非常高。

(4)本文研究樣本仍偏少，經驗閾值和識別精度有待進一步檢驗。

常旭，劉伊克.2002.地震記錄的廣義分維及其應用［J］.地球物理學報，45(6):839－846.

陳會忠，侯燕燕，何加勇，等.2011.日本地震預警系統日趨完善［J］.國際地震動態，(4):10－15.

李慧婷，黃文輝.2000.應用人工神經網絡方法識別近震與遠震［J］.華南地震，20(4):71－75.

劉希強，周蕙蘭，曹文海，等.2002.高斯線調頻小波變換極其在地震震相識別中的應用［J］.地震學報，24(6):607－616.

劉希強，周蕙蘭，沈萍，等.2000.用于三分向記錄震相識別的小波變換方法［J］.地震學報，22(2):125－131.

劉希強，周彥文，曲均浩，等.2009.應用單臺垂向記錄進行區域地震事件實時檢測和直達P波初動自動識別［J］.地震學報，31(3):260－271.

毛燕，崔建文，鄭定昌，等.2011.地震記錄的P波自動撿拾［J］.地震研究，34(1):47－51.

王海軍，靳平，劉貴忠，等.2003.區域震相初至估計［J］.西北地震學報，25(4):298－303.

王繼，陳九輝，劉啟元，等.2006.流動地震臺陣觀測初至震相的自動檢測［J］.地震學報，28(1):42 －51.

王娟，劉俊民，范萬春.2004.神經網絡在震相識別中的應用［J］.現代電子技術，27(8):35－37.

楊配新，鄧存華，劉希強，等.2004.數字化地震記錄震相自動識別的方法研究［J］.地震研究，27(4):308－313.

張范民，李清河.1998.利用人工神經網絡理論對地震信號及地震震相進行識別［J］.西北地震學報，20(4):43－49.

周彥文，劉希強，李鉑，等.2010.基于單臺P波記錄的快速自動地震定位方法研究［J］.地震研究，33(2):183－188.

朱元清，佟玉霞，于海英，等.2002.數字化臺網的近震震相自動識別［J］.西北地震學報，24(1):5 －12.

Allen R V.1978.Automatic earthquake recognition and timing from single traces［J］.BSSA，68(5):1521 － 1532.

Allen R V.1982.Automatic phase pickers:their present use and future prospects［J］.BSSA，72(68):225 －242.

Bai C Y，Kennett B L N.2000.Automatic phase-detection and identification by full use of single three-component broadband seismogram［J］.BSSA，90(1):187 －198.

Bear M，Kradolfer U.1987.An automatic phase picker for local and teleseismic events［J］.BSSA，77(4):1437 －1445.

Diehl T，Kissling E，Husen S，et al.2009.Consistent phase picking for regional tomography models:application to the Greater Alpine region［J］.Geophys J Int，176(2):542 － 554.

Haykin S.1996.Adaptive Filter Theory［M］.New Jersey:Upper Saddle River:1－989.

Kuperkoch L，Meier T，Lee J，et al.2010.Automated determination of P-phase arrival times at regional and local distances using higher order statistics［J］.Geophys J Int，181(2):1159 － 1170.

Maeda N.1985.A method for reading and checking phase times in autoprocessing system of seismic wave data［J］.Zisin，38(3):385－379.

Sleeman R，Van E T.1999.Robust automatic P-phase picking:An on －line implementation in the analysis of broadband seismogram recordings［J］.Phys Earth Planet Interi，113(1 － 4):265 － 275.

Wang J，Teng T L.1995.Artificial neural network-based seismic detector［J］.BSSA，85(1):308 －319.

Zhang H J，Thurber C，Rowe C.2003.Automatic P-wave arrival detection and picking with muhiscale wavelet analysis for single-component recordings［J］.BSSA，93(5):1904 －1912.

Detection of Regional Seismic Events by Kurtosis Method and Automatic Identification of Direct P-wave First Motion by Kurtosis-AIC Method

ZHAO Da-peng1，LIU Xi-qiang2，LI Hong2，ZHOU Yan-wen2
(1.Lanzhou Institute of Seismology，CEA，Lanzhou 73000，Gansu，China)
(2.Earthquake Administration of Shandong Province，Jinan 250014，Shandong，China)

Basing on Kurtosis-AIC and kurtosis methods of P-wave signal，we put forward a new method for real-time detection of regional earthquake event and automatic identification of direct P-wave first motion and apply it to process seismic data recorded by Shandong Seismic Network.The results show as follows:(1)The kurtosis method effectively detect earthquake events，and may effectively reduce false alarm and missing report rates;(2)Compared with phase arrival time results in manual identification，average absolute error of phase arrival time in automatic identification based on Kurtosis-AIC method is(0.09±0.08)s.

kurtosis method;Kurtosis-AIC method;earthquake identification;phase identification

P315.6

A

1000－0666(2012)02－0220－06

2011－02－14.

國家科技支撐計劃課題 (2012BAK19B04)和中國地震局地震科技星火計劃攻關項目 (XH12029)聯合資助.