馮紅武，顏文華，張煒超，王建昌

(陜西省地震局，陜西 西安 710068)

0 引言

震相識別是現(xiàn)代地震學研究中的重要課題之一，也是地震定位、地球內(nèi)部結(jié)構(gòu)、震源機制等一系列研究的基礎(chǔ)。隨著數(shù)字化地震臺網(wǎng)和地震科學臺陣大規(guī)模建立和發(fā)展，僅靠人工撿拾震相已滿足不了實際需要。震相的自動撿拾不但可以提高地震觀測數(shù)據(jù)的處理效率，而且提升了地震臺網(wǎng)的快速響應(yīng)能力，也是地震預(yù)警技術(shù)的一個關(guān)鍵環(huán)節(jié)。因此，研究更加有效和準確的震相自動識別方法是我們亟需解決的一個科學問題。

一般來說，地震事件的自動檢測是震相自動識別的基礎(chǔ)。目前，長短時平均比方法(STA/LTA方法)是一種既可判斷有無地震事件，又能自動撿拾初至震相到時的方法。這種方法因算法簡單、撿拾效率高、適應(yīng)性強等優(yōu)點，已被廣泛地應(yīng)用于地震波初至到時的自動撿拾[1-7]。但該方法缺點是當?shù)卣鸩ㄐ旁氡容^低或初動不清晰時，自動撿拾震相的效果會較差。除此之外，主要的震相自動識別方法還有分形分維法[8-9]、基于小波變化的主成分分析法[10-11]、人工神經(jīng)網(wǎng)絡(luò)法[12-14]、基于幅值和頻率的P波識別方法[15-16]等，都取得了一定的效果。

STA/LTA方法和AIC準則相結(jié)合的多步驟撿拾方法是目前P波震相自動撿拾的主要方法，該方法實質(zhì)是綜合了幾種單一算法的優(yōu)點，避免單一識別方法受使用范圍和條件限制的缺點，提高了自動識別的精度，是一種綜合分析法[17]。本文在STA/LTA方法和AIC準則相結(jié)合的多步驟撿拾方法的基礎(chǔ)上，針對地震波的變化特征，將地震信號的瞬時振幅和瞬時頻率作為基本參數(shù)，重新構(gòu)建識別P波和S波的STA/LTA方法的新特征函數(shù)，提出了基于新特征函數(shù)的STA/LTA和VAR-AIC多步驟撿拾P波和S波的方法，以期提高震相撿拾的準確性和可靠性，最后對流動地震科學臺陣記錄的地震震相進行離線拾取，來檢驗該方法的識別效果。

1 STA/LTA的特征函數(shù)構(gòu)建

地震波到達后，臺站記錄到的地震波初動振幅和頻率明顯區(qū)別于臺站記錄的背景噪聲，而且后至震相的振幅或頻率較初至波也有明顯的變化。為了準確撿拾P波和S波，STA/LTA特征函數(shù)應(yīng)能充分體現(xiàn)地震波到達后的振幅增大和頻率改變(變大或變小)的變化特征，同時考慮到P波和S波的變化特征一般分別在垂直向和水平向突出。為此，本文將瞬時振幅和瞬時頻率作為基本參量，構(gòu)建了拾取P波和S波的特征函數(shù)。

由于地震波是具有時變特性的非線性、非平穩(wěn)信號，我們利用美籍華裔科學家Huang等人提出的希爾伯特-黃變換方法(HHT)[18]計算地震波的瞬時頻率。HHT方法可將復(fù)雜信號分解為有限個固有模態(tài)分量，從而賦予了瞬時頻率合理定義和物理意義。

P波特征函數(shù)：

(1)

S波特征函數(shù)：

(2)

2 震相自動撿拾方法

2.1 STA/LTA法

STA/LTA方法是一種能量方法，其原理是用信號的短時平均值STA和長時平均值LTA的比值來反映信號水平或能量的變化。當?shù)卣鸩ǖ絹頃r，STA要比LTA變化的快，STA/LTA值相應(yīng)會有一個顯著的增加，若該值大于設(shè)定的觸發(fā)閾值，即可判定地震信號的到來。本文采用遞歸STA/LTA算法，計算公式如式(3)、式(4)：

(3)

(4)

式中STAi和LTAi分別為地震信號在時刻點i的短時平均值和長時平均值，CF(i)為地震信號在i時刻點的特征函數(shù)值，Nsta和Nlta分別為短時平均值和長時平均值的時間窗內(nèi)所包含的記錄點數(shù)。本文選取新構(gòu)建的P波和S波的特征函數(shù)。

2.2 基于AIC準則的震相識別方法

由于STA/LTA方法僅能粗略撿拾震相的到時，而且識別出的震相到時往往滯后于實際到時，即該方法得到的結(jié)果只是震相到達的大致位置。

1973年日本學者Akaike提出一個基本信息量的定階準則，即AIC準則[19]。震相到時檢測的自回歸(AR)技術(shù)是假設(shè)震相到時前后的地震記錄是兩個不同的穩(wěn)態(tài)過程，基于這個假設(shè)，Sleeman和orild于1999年提出了AR-AIC方法[20]用于對地震震相進行識別。不同于AR-AIC方法，Maeda[21]建議直接根據(jù)地震波形數(shù)據(jù)而不需要計算AR系數(shù)來得到AIC函數(shù)，該方法被稱為VAR-AIC方法。對于地震記錄x(i)(i=1,2,…,L)，將AIC檢測器定義為：

AIC(k)=klog{var(x[1,k])}+

(N-k-1)log{var(x[k+1,N])}

(5)

式中，k的范圍是地震圖某窗長內(nèi)所有的采樣點，AIC函數(shù)局部最小值對應(yīng)的時刻即為震相到時。

不同于在整個地震波信號上直接應(yīng)用AIC準則，本文是利用STA/LTA方法撿拾得到的粗撿拾點，在固定窗內(nèi)對該點進行AIC精確撿拾。

2.3 基于新特征函數(shù)的STA/LTA和VAR-AIC多步驟撿拾方法

本文對P波撿拾主要分為兩步：(1)粗略撿拾：利用STA/LTA方法對P波粗撿拾，特征函數(shù)選取本文新構(gòu)建的P波特征函數(shù)，為了避免不同閾值造成的誤觸發(fā)或漏觸發(fā)現(xiàn)象，直接取STA/LTA的最大值，最大值點出現(xiàn)時刻即為P波的粗撿拾點(即P波到達的大致位置)，記為tP；(2)精確撿拾：以P波粗撿拾點tP為中心，對其前推和后推一定時間窗，前后時間窗各取Δt時間長度，對時間窗長[tP-Δt,tP+Δt]的數(shù)據(jù)應(yīng)用VAR-AIC方法對P波到時進行精確撿拾。

由于S波容易受到P波尾波以及各種反射波的干擾，與P波撿拾相比，S波撿拾過程通常較復(fù)雜而且可靠性低。為了提高對S波撿拾的可靠性，對S波的撿拾主要分為兩步：(1)粗略撿拾：首先應(yīng)用STA/LTA方法對S波進行粗略撿拾，特征函數(shù)選取本文新構(gòu)建的S波特征函數(shù)。長、短時間窗從撿拾到的P波到時點開始，同步向前逐點滑動對S波特征函數(shù)曲線進行處理，將每個點的計算結(jié)果作為該點的STA/LTA值。同時為了解決閾值選擇的問題，直接取STA/LTA值的最大值，最大值點出現(xiàn)的時刻即為S波的粗撿拾到時(即S波到達的大致位置)；(2)精確撿拾：我們將P波段和S波段認為是兩個不同的穩(wěn)態(tài)過程，應(yīng)用AIC準則對S波到時進行精確撿拾。同樣以S波粗撿拾的到時tS為中心，對其前推和后推一定時間窗，前后時間窗各取Δt時間長度，對時間窗長[tS-Δt,tS+Δt]的數(shù)據(jù)應(yīng)用VAR-AIC方法對S波到時進行精確撿拾。

3 實際應(yīng)用

為驗證本文方法的有效性和準確性，我們將以上識別方法應(yīng)用于甘東南流動地震科學臺陣的地震記錄。由于文中的震相識別方法旨在用計算機自動識別來代替人工處理地震科學臺陣大量地震數(shù)據(jù)，我們感興趣的是對近震Pg和Sg震相的自動撿拾，按照甘青區(qū)域地震波走時表[22]，選擇震中距在360km以內(nèi)臺站記錄的波形數(shù)據(jù)。以臺陣記錄的2011年11月2日甘肅岷縣ML3.6地震為例，對其進行震相識別，記錄地震較清晰臺站數(shù)為29個，P波和S波有效數(shù)分別為29個和28個。

3.1 P波自動撿拾結(jié)果及分析

在應(yīng)用STA/LTA方法對P波粗略撿拾時，選取時間窗長，主要基于以下考慮：長時平均窗LTA描述的是相對于待撿拾信號的背景噪聲水平的平均大小，其取值應(yīng)能反應(yīng)背景噪聲水平；短時平均窗STA描述的主要是信號幅值的瞬時變化，一般取長于待撿拾信號的幾個周期左右，時間窗太短則對短周期的干擾更敏感，容易產(chǎn)生誤觸發(fā)，如果時間窗太長則顯示不出待撿拾信號的瞬時特征，容易產(chǎn)生漏觸發(fā)。因此，根據(jù)已有研究結(jié)果[7,16,23]，并且多次反復(fù)試驗，取長時間窗30s，短時間窗0.5s。對于觸發(fā)閾值的選取，若閾值過大，就有可能會出現(xiàn)漏觸發(fā)，若閾值太小，就對很多干擾會誤觸發(fā)，所以在選取閾值時，要考慮臺站的背景噪聲、地震震級、地震的遠近以及儀器頻帶范圍等，經(jīng)過試驗，設(shè)定的觸發(fā)閾值為5，且基本滿足要求。

精確撿拾時，考慮到地震數(shù)據(jù)原始采樣率為40Hz，以及多次反復(fù)試驗，以P波粗撿拾點tP為中心，對其前后各取2s，對時間窗[tP-2s,tP+2s]的數(shù)據(jù)應(yīng)用VAR-AIC方法對P波到時進行精確識別。對該地震事件P波進行識別，P波自動撿拾結(jié)果與人工識別結(jié)果見表1和圖1，P波震相檢測曲線見圖2?？梢钥闯?，該方法對P波的識別效果較好，人工識別和自動識別誤差的平均值為0.034s，誤差很小。

3.2 S波自動撿拾結(jié)果及分析

在應(yīng)用STA/LTA方法對S波粗略撿拾時，同樣基于以上考慮以及相關(guān)研究成果[23]，我們?nèi)￠L時間窗4s，短時間窗取0.5s，特征函數(shù)的觸發(fā)閾值設(shè)定為5。精確撿拾時，考慮到地震數(shù)據(jù)原始采樣率為40Hz，以及多次反復(fù)試驗，以S波粗撿拾點tS為中心，對其前后各取3s，對時間窗[tS-3s,tS+3s]的數(shù)據(jù)應(yīng)用VAR-AIC方法對S波到時進行精確識別。對該地震事件的S波進行自動識別，S波自動撿拾結(jié)果與人工識別結(jié)果見表2和圖3，S波震相檢測曲線見圖4。可以看出，除308臺站識別誤差較大，該方法對S波的識別效果較好，人工識別和自動識別誤差的平均值為0.105s，誤差很小。按照S波識別誤差小于0.2s的精度要求，識別有效率高達96.43%，有效識別平均誤差0.091s。

表1 P波到時識別結(jié)果比較

圖1 P波自動撿拾與人工撿拾的結(jié)果比較誤差圖Fig.1 The error graph of comparison of results between manual and automatic P-wave pickups

圖2 112臺站的P波撿拾曲線Fig.2 P-wave picking curve of Seismic station 112

臺站代碼人工識別/s自動識別/s誤差臺站代碼人工識別/s自動識別/s誤差101 33.150 33.1500.000210 16.825 17.000-0.175102 30.850 30.8250.025211 14.200 14.325-0.125 106 30.075 30.0750.000212 14.275 14.175 -0.100107 27.200 27.325-0.125213 15.325 15.500-0.175109 22.550 22.625-0.075214 19.025 19.125 -0.100 111 22.150 22.250-0.100215 22.050 22.150 -0.100 112 21.600 21.675-0.075217 29.700 29.850 -0.150113 22.675 22.700-0.025308 38.325 37.825 0.500 114 30.450 30.400 0.050312 27.700 27.825-0.125115 32.325 32.425 -0.100313 28.675 28.825 -0.150201 47.275 47.225 0.050314 29.875 30.075 -0.200 205 30.875 30.950-0.075701 14.375 14.475 -0.100 208 21.075 21.150-0.075702 14.675 14.800-0.125 209 16.750 16.7250.025703 18.425 18.450-0.025

圖3 S波自動撿拾與人工撿拾結(jié)果比較誤差圖Fig.3 The error graph of comparison of results between manual and automatic S-wave pickups

圖4 112 臺站的S波撿拾曲線Fig.4 S-wave picking curve of seismic station 112

4 結(jié)論

本文基于瞬時振幅和瞬時頻率參數(shù)，構(gòu)建了識別P波和S波的STA/LTA方法的新特征函數(shù)，其能夠充分體現(xiàn)地震初至波和后至震相到達后地震波幅值增大和頻率變化(變大或變小)的特征。主要結(jié)論如下：

(1)本文提出的基于新特征函數(shù)的STA/LTA和VAR-AIC多步驟識別P波和S波方法，將該方法應(yīng)用于甘東南地震科學臺陣數(shù)據(jù)，結(jié)果表明該方法對P波具有較高的識別精度，也對S波自動撿拾有較好效果。

(2)本文方法對P波和S波進行自動撿拾，識別效果雖然具有較高的精度，但考慮到震相撿拾的精度直接影響地震學一系列科學研究，其識別精度還需進一步提高，尤其是S波的撿拾精度。另外，本文的震相識別方法是多步驟撿拾，雖然克服了單一方法的缺陷，提高了震相識別的精度，但是運算速度慢，也比較費時。