曲保安 劉希強 蔡 寅， 范曉勇 林秀娜于慶民 趙 瑞 李 鉑 周彥文

1）中國山東泰安271000山東省地震局泰安基準地震臺

2）中國山東濟南250014山東省地震局

3）中國北京100081中國地震局地球物理研究所

引言

隨著數(shù)字地震觀測資料應用的日益廣泛和深入，地震學者對震相自動識別技術日趨關注，震相自動識別在海量資料分析和實時資料處理等領域得到了較快的發(fā)展.

1976年，Stevenson發(fā)表了其在地震自動分析方面的研究成果，其方法核心即為長短時平均能量比（short term average／long term average，簡寫為STA／LTA）（Stevenson，1976）.次年Stewart（1977）提出了用于近震實時檢測和定位的方法.1978年，Allen發(fā)表了應用單通道地震記錄進行地震自動檢測和到時識別的文章（Allen，1978）.1982年，Allen又對當時震相自動識別的應用現(xiàn)狀和未來發(fā)展進行了總結(jié)和展望，同時提出了用特征函數(shù)作為輸入?yún)?shù)的STA／LTA算法（Allen，1982）.該特征函數(shù)思想經(jīng)過不斷發(fā)展目前仍在廣泛應用（李山有等，2006；高淑芳等，2008）.迄今，經(jīng)過30多年的發(fā)展，各種新的震相自動識別算法隨著數(shù)字信號處理技術、概率論以及自動控制理論的發(fā)展被廣泛提出和應用.我們將震相自動識別分為3步.第一步：預處理.是對原始信號進行各種濾波和時頻轉(zhuǎn)換，包括傅里葉變換、Gaber變換、Wigner-Ville分布、小波變換、小波包變換、Hilbert-Huang變換等預處理算法（劉希強等，1998，2000，2002；楊配新等，2004；Granet，1983）.第二步：震相識別算法.又可分為能量成分算法（如STA／LTA算法）、頻率成分算法（如瞬時頻率算法）、波形算法（如AIC算法）、波動矢量差異算法（如偏斜角）和動態(tài)特征算法（如偏振分析、線性度、功率譜和頻譜）等5部分（Roberts et al，1989）.第三步：方法綜合判定.對第二步的算法進行各種組合和閾值設定，并依據(jù)各種情形確定相應的識別程序及方法，包括后期的用人工神經(jīng)網(wǎng)絡對判定方法進行訓練（Baer，Kradolfer，1987）.經(jīng)過上述3步的震相自動識別方法在P波震相識別方面得到了廣泛的應用，并取得了良好的效果.但是在S波自動識別方面，由于P波尾波以及各種反射震相，如PmP的強干擾，STA／LTA算法用于S波震相自動識別的效果較差，Roberts等（1989）基于S波與P波質(zhì)點振動方向的差異，提出了用偏振分析方法識別S波.Cichowicz（1993）提出了采用偏斜角、偏振度和切向能量與總能量比值3種參數(shù)相乘得到的值作為判斷依據(jù)的方法，該方法總有效性達到65%—70%.Bai和Kennett（2000）采用能量分析、自回歸表征和瞬時頻率技術并與偏振分析結(jié)合進行綜合震相識別，此方法對震中距17.3°—41.6°的各種遠震震相的識別能力較高，包括遠震S波震相識別，但用于近震識別S波震相識別精度較低.Diehl等（2009）提出了用于地震層析成像的S波自動識別方法，該方法需經(jīng)過STA／LTA判定、偏振分析、坐標系變換、自回歸Akaike信息準則（auto regressive-Akaike information criterion，簡寫為AR-AIC）等一系列綜合判定，計算量較大，可用于資料的后續(xù)分析.經(jīng)過對2 500個S波震相識別結(jié)果與人工分析結(jié)果比較，該方法最大誤差大約為0.27s.Amoroso（2012）對臺網(wǎng)資料后續(xù)分析的S波自動識別采用了偏振濾波與波形一致性方法.Kuperkoch等（2012）采用波形自回歸預測方法，對地方震、近震和遠震的S波進行了自動判定，識別的平均誤差為0.5s±0.8s.

本文旨在研究根據(jù)實時地震波形進行地震早期預警方面的快速S波震相識別.依據(jù)馬強（2008）的地震預警動態(tài)定位的6種情形，S波到時圓包含地震預警動態(tài)定位的重要信息，對S波進行識別是地震預警動態(tài)定位的必要環(huán)節(jié).本文嘗試性地提出一種S波自動識別的方法，秉承“濾波越少，算法越好”（“the less filtering，the better the algorithm”）的思想（Douglas，1997），對原始三分向記錄波形不做任何濾波處理，由于已知由其它自動識別方法得到的P波到時，從而對S波到時可進行快速自動高精度的識別.對于P波的精確識別，推薦采用劉希強等（2009）的三階統(tǒng)計Akaike信息準則（three order cumulative-Akaike information criterion，簡寫為TOC-AIC）方法.

1 識別方法

本文方法的基本思路是根據(jù)由其它方法自動識別的P波的動態(tài)特征，求得波形的波動矢量，即P波的質(zhì)點振動方向，并選取計算窗長對后續(xù)波形進行質(zhì)點振動方向滑動計算，求得后續(xù)波形質(zhì)點振動方向與P波質(zhì)點振動方向的夾角，即偏斜角；將偏斜角和水平能量占總能量的比值的平方積作為確定識別區(qū)間的依據(jù)；然后根據(jù)方差Akaike信息準則（variance-Akaike information criterion，簡寫為VAR-AIC）判定兩個水平分向S波初動位置，并對兩個水平分向各自的識別結(jié)果進行分析判斷，確定S波震相到時.

1.1 確定窗長

鑒于震中距小于100km的地方震，其P波的周期約為0.05—0.2s，S波的周期約為0.1—0.5s（朱元清等，2002）.根據(jù)山東省地震臺網(wǎng)密度，以及Nyquist采樣定理，選取P波到時之后0.5s的波形數(shù)據(jù)，計算P波的卓越頻率（Boore，1983；Andrews，1986）

式中，fP為P波卓越頻率，m1和m2計算公式為

式中，v2（f）和D2（f）分別為速度功率譜密度和相應的位移功率譜密度.將fP／2作為S波的卓越頻率，計算窗長

式中，fS為采樣率，lw為計算波形偏振方向的窗長.

鑒于本研究所用資料為山東數(shù)字地震臺網(wǎng)記錄的數(shù)據(jù)，為速度記錄，所以在計算m2之前要將速度記錄仿真為位移記錄.本文采用的數(shù)值積分方法為Newton-Cotes求積算法

從計算精度和計算量兩方面綜合分析，本文選擇n＝4時的Newton-Cotes公式

該計算結(jié)果具有5次代數(shù)精度，其截斷誤差為

同時，對窗長進行限定，大于采樣率的1／5，小于采樣率的一半（即Nyquist頻率）.

1.2 計算偏斜角

選取P波之后lw長度的數(shù)據(jù)計算其最大特征值對應的特征向量，即P波的偏振方向（質(zhì)點振動方向），單臺三分向lw長度的數(shù)據(jù)協(xié)方差矩陣為

長度為lw的任意兩組變量x和y的協(xié)方差為

從lw長度之后開始進行逐點滑動，求得各個數(shù)據(jù)窗最大特征值對應的特征向量，并計算與P波偏振方向的夾角.

1.3 水平能量與總能量比值

水平能量與總能量比值，即窗內(nèi)兩個水平分向能量和與三分向能量和的比值，計算公式為

式中，x和y為兩個水平分向波形數(shù)據(jù)序列，z為垂直分向波形數(shù)據(jù)序列.

1.4 VAR-AIC

Maeda（1985）提出了直接根據(jù)地震圖記錄而得到AIC函數(shù)的方法，稱為VAR-AIC方法.AIC函數(shù)的最小值對應的時間就是地震震相初至.AIC函數(shù)表示為

式中，var（x［1，k］）和（var（x［k＋1，N］）分別表示兩個時間段內(nèi)數(shù)據(jù)的方差.

1.5 判定標準

對于S波的初步判定，選取偏斜角和水平能量與總能量比值的平方積作為特征函數(shù)（cf），即

閾值選擇自P波開始至計算點位置的均值的5倍加方差的5倍.以此特征函數(shù)作為判定S波識別區(qū)間的判據(jù)，對S波進行初步判定.之后選取初判位置前3個lw長度和后3個lw長度的數(shù)據(jù)作為精確判定的區(qū)間，采用VAR-AIC方法分別對NS和EW分向進行S波的精確識別.對于識別的兩個分向各自的S波位置，如果時間差小于0.1s，則認為二者的差值是由于S波分裂造成，則取二者的平均值作為識別的S波位置；如果二者的差值大于0.1s，則認為某一個分向存在誤識別，分別計算兩個分向各自的識別位置前后lw長度信號的信噪比，取信噪比較大的一個分向的判定結(jié)果作為S波的識別位置.

2 資料和結(jié)果分析

結(jié)合山東數(shù)字地震臺網(wǎng)記錄的數(shù)據(jù)，將本文震相識別方法的結(jié)果與人機交互處理結(jié)果進行對比，最后確定其有效性.

按照震級不同選取29個地震事件，每個地震事件選擇震中距100km范圍內(nèi)的臺站記錄，要求記錄可以精確識別P波到時，共選取56個臺站的118個三分向記錄.通過這118個三分向記錄波形對本文方法進行了測試.本研究所涉及的地震及臺站分布見圖1.本文所選取的地震事件參數(shù)見表1.地震選取考慮以下3方面：① 信噪比隨震級的變化；② 信噪比隨震中距變化；③ 信噪比受臺站背景噪聲及地震射線傳播路徑影響.此方法統(tǒng)籌同一臺站對不同震級地震的識別能力，又兼顧全省臺站各自背景噪聲的不同，從而在全省不同位置選取地震事件，同時也在相近震中位置選取不同震級的事件（煙臺萊州震群和濮陽范縣附近4次不同震級地震）.

圖1 本研究涉及的地震震中和臺站分布圖實心三角形表示臺站，實心圓形表示地震震中位置Fig.1 Distribution of earthquakes and stations used in this studySolid triangles denote seismic stations，solid circles denote the epicenters of earthquake

表1 本文所選取的地震事件參數(shù)Table 1 Detailed information of selected earthquake events

為了體現(xiàn)該方法的實時性，在對三分向記錄進行S波識別時，采用與現(xiàn)有Jopens系統(tǒng)類似的方法，將三分向數(shù)據(jù)以數(shù)據(jù)包形式分割發(fā)送，每秒1個包，模擬從數(shù)采實時獲取數(shù)據(jù).

圖2為118個三分向記錄的自動識別結(jié)果與人機交互識別結(jié)果的對比，即以人機交互結(jié)果為標準，自動識別結(jié)果的誤差分布圖.圖中橫軸為二者的差值，縱軸為S波對P波尾波的信噪比.對118個三分向記錄的識別結(jié)果按照S波對P波尾波的信噪比不同進行了統(tǒng)計（表2）.對于誤差大于0.1s的震相識別作為無效識別.本文所采用的118個三分向記錄，其信噪比均小于12.最高信噪比為11.56，最低信噪比為－2.28.

如表2所示，本文所提出的S波實時自動識別方法對于信噪比大于5的地震事件，按照識別誤差小于0.1s的精度要求，識別有效率高達89.39%，有效識別平均誤差為0.05s.對于低信噪比大于0小于5的三分向地震記錄，按照識別誤差小于0.1s的精度要求，識別有效率達到65.79%，有效識別平均誤差為0.07s.綜上，本文識別方法在信噪比相對較低的情況下能夠?qū)崿F(xiàn)較高的有效率和較小的識別誤差.

圖2 118個三分向記錄的自動識別結(jié)果與人工交互識別結(jié)果對比Fig.2 Comparison between automatic identification results and human-computer interaction recognition results for the 118three-component records

表2 識別結(jié)果統(tǒng)計Table 2 Statistics on identification results

3 討論與結(jié)論

本文提出了一種用于地震早期預警的S波震相實時自動識別方法.該方法不對原始信號進行任何濾波處理，直接根據(jù)偏斜角和水平能量與總能量比值確定識別區(qū)間；采用VAR-AIC方法對兩個水平分向分別進行精確識別；對兩個識別結(jié)果進行判斷，確定S波初動時刻.

經(jīng)過對118個三分向記錄的實際應用驗證，本文方法具有識別有效率高、識別誤差小的優(yōu)點.由于P波尾波以及各種反射震相的強干擾，S波的信噪比整體偏低，不利于S波有效識別及識別精度的提高.本文方法對于信噪比大于5的地震記錄可以實現(xiàn)S波震相高精度的實時有效識別.雖然與目前已有的P波0.02s的識別精度還有一定差距，但0.1s的識別誤差小于Amoroso等（2012）和Kuperkoch等（2012）的識別誤差，而且本方法識別有效率高達89.39%.

本文方法是在基于P波已經(jīng)由其它方法精確識別的基礎上提出的.其P波卓越頻率、窗長及P波偏振方向的計算均受P波的識別精度影響，但是鑒于目前P波識別精度已達到0.02s，而且本文選擇窗長為0.5s，故本方法能夠達到較高的識別精度.

高淑芳，李山有，武東坡，馬強.2008.一種改進的STA／LTA震相自動識別方法［J］.世界地震工程，24（2）：37－－41.

Gao S F，Li S Y，Wu D P，Ma Q.2008.An automatic identification method of seismic phase based on modified STA／LTA algorithm［J］.World Earthquake Engineering，24（2）：37－－41（in Chinese）.

李山有，朱海燕，武東坡，宋晉東.2006.基于振幅和瞬時頻率的震相自動識別方法［J］.世界地震工程，22（4）：1－－4.

Li S Y，Zhu H Y，Wu D P，Song J D.2006.Automatic recognition of seismic phase based on amplitude and instantaneous frequency［J］.World Earthquake Engineering，22（4）：1－－4（in Chinese）.

劉希強，周蕙蘭，鄭治真，沈萍，楊選輝，馬延路.1998.基于小波包變換的弱震相識別方法［J］.地震學報，20（4）：373－－380.

Liu X Q，Zhou H L，Zheng Z Z，Shen P，Yang X H，Ma Y L.1998.Identification method of weak seismic phases on the basis of wavelet packet transform［J］.Acta Seismologica Sinica，20（4）：373－－380（in Chinese）.

劉希強，周蕙蘭，沈萍，楊選輝，馬延路，李紅.2000.用于三分向記錄震相識別的小波變換方法［J］.地震學報，22（2）：125－－131.

Liu X Q，Zhou H L，Shen P，Yang X H，Ma Y L，Li H.2000.Identification method of seismic phase in three-component seismograms on the basis of wavelet transform［J］.Acta Seismologica Sinica，22（2）：125－－131（in Chinese）.

劉希強，周蕙蘭，曹文海，李紅，李永紅，季愛東.2002.高斯線調(diào)頻小波變換及其在地震震相識別中的應用［J］.地震學報，24（6）：607－－616.

Liu X Q，Zhou H L，Cao W H，Li H，Li Y H，Ji A D.2002.Gauss linear frequency modulation wavelet transform and its application to seismic phase identification［J］.Acta Seismologica Sinica，24（6）：607－－616（in Chinese）.

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

Liu X Q，Zhou Y W，Qu J H，Shi Y Y，Li B.2009.Real-time detection of regional events and automatic P-phase identification from the vertical component of a single station record［J］.Acta Seismologica Sinica，31（3）：260－－271（in Chinese）.

馬強.2008.地震預警技術研究及應用［D］.哈爾濱：中國地震局工程力學研究所：40－－46.

Ma Q.2008.Study and Application on Earthquake Early Warning［D］.Harbin：Institute of Engineering Mechanics，China Earthquake Administration：40－－46（in Chinese）.

楊配新，鄧存華，劉希強，任勇，顏其中.2004.數(shù)字化地震記錄震相自動識別的方法研究［J］.地震研究，27（4）：308－－313.

Yang P X，Deng C H，Liu X Q，Ren Y，Yan Q Z.2004.Studies on auto-distinguishing phase of digital seismic recording［J］.Journal of Seismological Research，27（4）：308－－313（in Chinese）.

朱元清，佟玉霞，于海英，宋俊高.2002.數(shù)字化臺網(wǎng)的近震震相自動識別［J］.西北地震學報，24（1）：5－－12.

Zhu Y Q，Tong Y X，Yu H Y，Song J G.2002.Automatic recognition of seismic phases of regional earthquake［J］.Northwestern Seismological Journal，24（1）：5－－12（in Chinese）.

Allen R V.1978.Automatic earthquake recognition and timing from single trace［J］.Bull Seismol Soc Am，68（5）：1521－－1532.Allen R V.1982.Automatic phase pickers：Their present use and future prospects［J］.Bull Seismol Soc Am，72（6B）：225－－242.

Amoroso O，Maercklin N，Zollo A.2012.S-wave identification by polarization filtering and waveform coherence analyses［J］.Bull Seismol Soc Am，102（2）：854－－861.

Andrews D J.1986.Objective determination of source parameters and similarity of earthquakes of different size［M］∥Das S，Boatwright J，Scholz C eds.Earthquake Source Mechanics.Washington D C：American Geophysical Union：20－－23.

Bai C Y，Kennett B L N.2000.Automatic phase detection and identification by full use of single three-component broadband seismogram［J］.Bull Seismol Soc Am，90（1）：187－－198.

Baer M，Kradolfer U.1987.An automatic phase picker for local and teleseismic events［J］.Bull Seismol Soc Am，77（4）：1437－－1445.

Boore D M.1983.Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra［J］.Bull Seismol Soc Am，73（6A）：1865－－1894.

Cichowiez A.1993.An automatic S-phase picker［J］.Bull Seismol Soc Am，83（1）：180－－189.

Diehl T，Deichmann N，Kissling E，Husen S.2009.Automatic S-wave picker for local earthquake tomography［J］.Bull Seismol Soc Am，99（3）：1906－－1920.

Douglas A.1997.Bandpass filtering to reduce noise on seismograms：Is there a better way？［J］.Bull Seismol Soc Am，87（3）：770－－777.

Granet M.1983.An automatic seismic signal detection based on linear prediction filter theory［J］.Ann Geophys，1：109－－114.Kuperkoch L，Meier T，Brustle A，Lee J，F(xiàn)riederich W，EGELADOS Working Group.2012.Automated determination of S-phase arrival times using autoregressive prediction：Application to local and regional distances［J］.Geophys J Inter，188（2）：687－－702.

Maeda N.1985.A method for reading and checking phase times in autoprocessing system of seismic wave data［J］.Zisin，38（6）：365－－379.

Roberts R G，Christoffersson A，Cassidy F.1989.Real-time event detection，phase identification and source location estimation using single station three-component seismic data［J］.Geophys J Int，97（3）：471－－480.

Stevenson P R.1976.Microearthquakes at Flathead Lake，Montana：A study using automatic earthquake processing［J］.Bull Seismol Soc Am，66（1）：61－－80.

Stewart S W.1977.Real-time detection and location of local seismic events in central California［J］.Bull Seismol Soc Am，67（2）：433－－452.