田優(yōu)平 趙愛華

1) 中國長沙410004湖南省地震局 2) 中國北京100081中國地震局地球物理研究所

?

基于小波包和峰度赤池信息量準(zhǔn)則的P波震相自動識別方法

1) 中國長沙410004湖南省地震局 2) 中國北京100081中國地震局地球物理研究所

基于小波包變換和峰度赤池信息量準(zhǔn)則(AIC)， 提出了一種新的自動識別P波震相的綜合方法， 即小波包-峰度AIC方法． 首先對由加權(quán)長短時窗平均比(STA/LTA)法粗略確定的P波到時前后3 s的記錄進行小波包三尺度的分解與重構(gòu)， 分別計算每個尺度重構(gòu)信號的峰度AIC曲線并將其疊加， 疊加曲線的最小值則為P波震相到時; 然后對原始地震記錄進行有限沖激響應(yīng)自適應(yīng)濾波以提高信噪比和識別精度; 最后將小波包-峰度AIC方法應(yīng)用到合成理論地震圖及實際地震記錄的P波初至自動識別中. 結(jié)果表明： 初至清晰度對識別精度的影響比信噪比對其影響更大； 與單獨使用加權(quán)STA/LTA方法和峰度AIC法相比， 小波包-峰度AIC法具有更強的抗噪能力， 識別精度更高； 當(dāng)初至清晰時， 小波包-峰度AIC法自動識別與人工識別的P波到時平均絕對差值為(0.077±0.075) s．

P波震相 有限沖激響應(yīng)(FIR)數(shù)字濾波 震相自動識別 小波包-峰度AIC方法

引言

地震震相的檢測和識別是地震學(xué)研究中的重要課題， 是地震定位、 地震預(yù)警、 地球內(nèi)部結(jié)構(gòu)、 層析成像、 震源機制等一系列地震研究的基礎(chǔ)． 以往的震相識別主要是人工識別， 由于其在很大程度上依賴經(jīng)驗， 所以費時、 費工且又繁雜． 隨著全球?qū)掝l帶地震儀的發(fā)展， 地震記錄中所包含的信息越來越豐富， 傳統(tǒng)的人工識別已不能滿足海量地震記錄的分析， 因此， 震相自動識別方法應(yīng)運而生． 該方法能節(jié)省大量時間， 大大提高地震速報和預(yù)警的效率， 為震后應(yīng)急救援贏得寶貴時間．

近年來， 基于震相的不同特征發(fā)展起來的震相自動識別方法可歸納為單特征法、 多特征法和綜合分析方法等3類．

單特征法是根據(jù)震相的能量、 頻譜、 偏振和分形維等單一特征的變化進行震相的判別和拾取， 其中最常用的是長短時窗平均比(short term average/long term average， 簡寫為STA/LTA)方法(Stevenson， 1976； Allen， 1978， 1982； Baer， Kradolfer， 1987； 劉晗， 張建中， 2014)． 該方法快速簡單， 但當(dāng)信噪比較低或初動不明顯時， 效果較差甚至出現(xiàn)誤判， 且拾取遠震信號到時的精度較低． 近年來震相識別方法以頻譜分析方法(如短時傅里葉變換、 小波及小波包變換、 S變換和Cohen類時頻分布)中的小波和小波包變換(Yomogida， 1994； 劉希強等， 1998， 2002； 王喜珍， 2004； 楊配新等， 2004)為主． 小波變換克服了短時傅里葉變換的單分辨率分析的不足， 具有可調(diào)的時頻窗， 在低頻區(qū)頻率分辨率較高， 高頻區(qū)時間分辨率較高； 小波包變換則對小波變換進行了改進， 具有將隨尺度增長而變寬的頻譜窗進一步分割變細的優(yōu)良特性， 其應(yīng)用前景較好． 此外， 對赤池信息量準(zhǔn)則(Akaike Information Criterion， 簡寫為AIC)的研究也較多， 依據(jù)AIC尋求曲線極小點作為背景噪聲與有效信號的最佳劃分點即震相的到時點. 基于AIC發(fā)展起來的一系列震相識別方法， 如基于預(yù)測模型的自回歸AIC(autoregressive-AIC， 簡寫為AR-AIC)方法(Sleeman， van Eck， 1999； 王海軍等， 2003； Küperkochetal， 2012)、 無需計算自回歸階數(shù)的方差A(yù)IC(variance-AIC， 簡寫為VAR-AIC)方法(Maeda， 1985； 曲保安等， 2014)、 直接根據(jù)地震圖記錄計算AIC函數(shù)的三階統(tǒng)計AIC(three-order cumulative-AIC， 簡寫為TOC-AIC)方法(劉希強等， 2009)、 基于峰度的AIC(Kurtosis-AIC， 簡寫為Kur-AIC)方法(Küperkochetal， 2010， 趙大鵬等， 2012)以及基于小波變換的AIC(wavelet-AIC， 簡寫為W-AIC)方法(Zhangetal， 2003)． 除上述方法外， 單特征法還有分形分維(Boschettietal， 1996； 常旭， 劉伊克， 2002； 王海軍等， 2004)、 高階統(tǒng)計量(Saragiotisetal， 2002； Küperkochetal， 2010)、 偏振分析(Cichowiczetal， 1988; Jurkevics， 1988)等方法．

多特征法主要是基于地震波傳播的整體特征或多個特征來識別震相， 如全波震相分析法(孫進忠等， 1985； 趙鴻儒等， 1990)、 相關(guān)法(王繼等， 2006； Satrianoetal， 2008)和人工神經(jīng)網(wǎng)絡(luò)法等， 其中人工神經(jīng)網(wǎng)絡(luò)中的誤差反向傳播神經(jīng)網(wǎng)絡(luò)應(yīng)用較廣(莊東海等， 1994； Dai， MacBeth， 1997； 王娟等， 2004)．

綜合分析方法一般先利用某種單特征法或其改進算法對震相到時進行粗略估計， 在此基礎(chǔ)上結(jié)合其它一種或多種單特征法得出各震相比較精確的到時． Bai和Kennett(2000)聯(lián)合STA/LTA、 自回歸分析、 瞬時頻率和偏振分析等4種方法來自動識別寬頻帶地震記錄的P波和S波等多種震相并拾取其到時， 大大提高了震相拾取的能力． 毛燕等(2011)通過STA/LTA、 幅值和頻率特征初步確定P波到時， 再用最小二乘法和AIC方法精確地確定P波到時．

在進行不同震相的識別時, 每種方法都各有其獨特的優(yōu)勢， 但也都或多或少存在一些缺陷． 單特征法在地震記錄信噪比高時能取得較好的結(jié)果， 但信噪比低或初動不清晰時拾取精度較低； 多特征法中全波震相分析法和人工神經(jīng)網(wǎng)絡(luò)法需要較大的工作量， 相關(guān)法則主要受續(xù)至波和時窗范圍的影響較大； 綜合分析方法拾取震相的精度相對較高， 但往往計算量較大． 震相自動識別是一項傳統(tǒng)而又持續(xù)發(fā)展的研究課題， 目前還沒有完全實現(xiàn)高效的地震震相自動識別． 為此， 本文在前人研究的基礎(chǔ)上， 提出了小波包-峰度AIC方法來自動識別P波震相到時， 并將該方法應(yīng)用到合成理論地震圖和云南地區(qū)近震P波到時的自動識別中， 以進一步探索提高P波自動識別精度的方法．

1 方法

本文提出小波包-峰度AIC方法來自動識別P波震相到時， 其基本思路為： 首先利用加權(quán)STA/LTA方法自動檢測出有效的地震事件并粗略地拾取P波初至到時； 然后對粗略拾取的到時前后各推3 s的時間窗內(nèi)的信號進行小波包三尺度分解重構(gòu)； 最后分別計算3個尺度重構(gòu)信號的峰度AIC曲線并進行疊加， 疊加的AIC曲線的最小值即為最終拾取到的精細P波初至到時．

1.1 加權(quán)STA/LTA方法粗略拾取P波到時

STA/LTA方法是一種常用的震相識別方法， 具有計算速度快、 簡單易行的優(yōu)點， 其基本原理是用信號短時窗平均值(STA)與長時窗平均值(LTA)之比來反映信號能量或振幅的變化． 當(dāng)信號到達時， STA比LTA變化快， 相應(yīng)的STA/LTA值明顯增加， 當(dāng)該比值大于事先設(shè)定的閾值時， 則視為地震事件發(fā)生， 該點被判定為初至點， 其對應(yīng)的時間為初至到時．

STA/LTA的計算方法分為標(biāo)準(zhǔn)STA/LTA和遞歸STA/LTA兩種， 其中后者速度快且結(jié)果平滑， 因而應(yīng)用普遍． 遞歸STA/LTA計算公式為

(1)

(2)

式中， STAi和LTAi分別表示信號在i時刻的短、 長時窗平均值，fc(i)為信號在i時刻的特征函數(shù)值，NSTA和NLTA分別表示短、 長時窗包含的記錄點數(shù)． 特征函數(shù)fc選用能同時反映振幅和頻率變化的公式

(3)

來計算. 式中，Y(i)為i時刻的地震記錄值．

當(dāng)信噪比低或初動不明顯時， STA/LTA算法的結(jié)果不太理想． 為此， 武東坡(2004)引入加權(quán)系數(shù)α對STA/LTA方法進行改進. 設(shè)STA/LTA比值為R， 改進后的計算式為

(4)

式中，d為信號的采樣間隔，t為時間間隔(t取1—2 s可達理想結(jié)果)． 當(dāng)有效信號未到時，α值接近1； 當(dāng)信號到達時,α值會突然增大， 從而對式(4)中的R值產(chǎn)生顯著的放大作用， 有利于更好地識別信噪比低或初動不明顯的信號．

本文取短時窗窗長為0.2 s， 長時窗窗長為2 s， 滑動步長為1個采樣點， 閾值為8． 運用加權(quán)STA/LTA方法可拾取P波初至到時， 并將其作為粗略到時．

1.2 峰度AIC方法

隨機變量x的k階統(tǒng)計量可表示為

(5)

隨機變量的峰度(Kurtosis)可衡量分布曲線的尖峭程度， 其表達式為

(6)

用兩個時間段數(shù)據(jù)的峰度作為特征函數(shù)fc， 可得峰度AIC的表達式為

(7)

式中，N為特征函數(shù)fc的長度．

1.3 一維小波包多尺度分解

Wickerhauser(1992)提出了離散小波包變換公式． 定義算子F0和F1為

(8)

式中， {hm}m∈Z和{gm}m∈Z分別為低通和高通濾波器．

離散小波包變換的表達式可寫為

(9)

圖1 小波包三尺度分解樹 Fig.1 Three-scale wavelet packet decomposition tree

由此可見， 小波包具有將隨尺度增長而變寬的頻譜窗進一步分割變細的優(yōu)良特性， 它克服了小波變換隨尺度增大而導(dǎo)致的頻譜局部性變差的缺陷， 能夠更精細地對信號進行分析， 更好地描述由于新震相產(chǎn)生而導(dǎo)致地震信號漸變或突變的局部特征．

圖2 利用小波包-峰度AIC方法拾取P波到時(a) 垂直向地震波形； (b) 尺度1； (c) 尺度2； (d) 尺度3； (e) 3個尺度疊加

1.4 小波包-峰度AIC方法精確拾取P波到時

圖2a給出了云南省元江臺記錄的2009年7月2日13時2分峨山ML2.8地震垂直向波形. 通過加權(quán) STA/LTA 方法粗略拾取該記錄的P波初至到時tP為12.74 s， 與人機交互拾取的到時12.66 s相差0.08 s．

本文首先對地震信號進行小波包三尺度分解， 利用分解的系數(shù)對原始信號進行重構(gòu)； 然后根據(jù)式(7)計算3個尺度下[tP-3 s，tP+3 s]時間窗內(nèi)重構(gòu)信號的峰度AIC值， AIC曲線最低點即對應(yīng)P波到時． 圖2b--d分別給出了3個尺度下的拾取到時， 分別為12.69 s， 12.67 s和12.67 s． 可以看出， 3個尺度的到時相差很小， 說明地震信號在多個尺度下是相關(guān)的(噪聲一般隨著尺度增加而衰減或不相關(guān))． 將3個尺度的峰度AIC曲線疊加得到的曲線最低點作為最終拾取到的精細P波到時. 如圖2e所示， 小波包-峰度AIC方法自動拾取的P波到時為12.67 s， 與人機交互拾取的到時12.66 s相差0.01 s． 可見， 在STA/LTA方法的基礎(chǔ)上， 本文方法進一步提高了P波識別的精度．

2 理論應(yīng)用

將本文提出的小波包-峰度AIC方法應(yīng)用到理論地震記錄中， 以檢驗該方法的有效性． 圖3a給出了云南省楚雄臺記錄的2008年8月31日16時31分地震的垂直向波形， 其合成理論地震圖(Wang， 1999； Wang， Wang， 2007)如圖3b所示． 運用射線追蹤法(趙愛華等， 2003； Zhaoetal， 2004)計算出合成地震記錄(圖3b)的理論走時為23.37 s． 分別利用加權(quán)STA/LTA方法、 峰度AIC方法和小波包-峰度AIC方法(本文方法)自動拾取圖3b中P波初至到時， 結(jié)果分別為23.54 s， 23.42 s， 23.36 s， 如圖4所示， 其與理論到時的絕對差值分別為0.17 s， 0.05 s， 0.01 s．

圖3 實際地震記錄(a)與合成理論地震圖(b) Fig.3 Real seismic recording (a) and theoretically synthetic seismogram (b)

向圖3b所示的合成理論地震記錄中逐漸加入不同信噪比(signal-to-noise ratio， 簡寫為SNR)的高斯白噪聲和實際地震噪聲， 以檢測各方法的抗噪能力及其拾取P波初至精度的效果． 3種方法拾取P波到時的情況如表1和表2所示． 可以看出， 在合成理論記錄中， 無論加入高斯白噪聲還是實際地震噪聲， 3種方法所拾取的P波到時精度隨著信噪比的降低均有一定影響． STA/LTA方法受信噪比影響最大， 隨著信噪比降低， 自動拾取與理論計算的P波到時差值ΔT逐漸增大， SNR≥20 dB時ΔT基本在0.3 s以內(nèi)； 峰度AIC方法隨著信噪比減小， ΔT也有一定波動， SNR≥13 dB時ΔT基本在0.3 s以內(nèi)； 小波包-峰度AIC方法受信噪比影響最小， SNR≥10 dB時ΔT基本在0.3 s以內(nèi)． 但當(dāng)信噪比降低到一定值(表1中加入地震噪聲SNR<5 dB， 表2中加入高斯白噪聲SNR<7.5 dB)時， 由于STA/LTA曲線很難達到事先設(shè)定的閾值致使檢測不到地震信號， 因而該方法失效．

圖4 合成地震記錄的P波到時自動拾取(a) 地震垂直向記錄； (b) STA/LTA方法； (c) 峰度AIC方法； (d) 小波包-峰度AIC方法

SNR/dB自動拾取的P波初至到時/sSTA/LTA方法峰度AIC方法小波包?峰度AIC方法自動與理論到時差ΔT/sSTA/LTA方法峰度AIC方法小波包?峰度AIC方法25．0023．5523．4323．380．180．060．0122．0023．5723．4423．390．200．070．0221．0023．5823．4523．390．210．080．0220．0023．5923．4623．390．220．090．0217．0023．6423．5123．420．270．140．0515．0023．6723．5523．470．300．180．1013．6123．7023．5923．510．330．220．1413．0023．7123．6123．540．340．240．1712．0023．7823．6323．560．410．260．1910．0023．9823．6923．600．610．320．239．0024．0423．7123．610．670．340．248．0024．0823．7523．620．710．380．257．0024．1923．7823．640．820．410．276．0024．2023．8123．670．830．440．305．0024．3323．8523．720．960．480．35

總之， 本文提出的小波包-峰度AIC方法在合成理論地震記錄中加入不同噪聲時表現(xiàn)為抗噪能力更強且拾取的P波到時精度更高； 當(dāng)加入的實際地震噪聲SNR≥8 dB時， ΔT≤0.25 s， 效果要明顯優(yōu)于其它兩種方法．

表2 合成理論記錄中加入高斯白噪聲時3種方法在不同信噪比下拾取P波到時的效果對比

3 實際應(yīng)用

將小波包-峰度AIC方法應(yīng)用到實際地震中， 以檢驗其在實際中的應(yīng)用效果．

3.1 數(shù)據(jù)

本文對研究區(qū)(97.5°E—107°E， 20°N—30°N)2009年2月—2011年6月107次ML1.4—4.9地震的722個單臺垂直向記錄進行研究和分析， 涉及的臺站共計52個． 選取單臺記錄的震中距為0.1°—2.29°(約11.12—254.64 km)， 大部分分布在0.1°—1.6°(約11.12—177.91 km)， 均為近震； 選取事件的震源深度為2.5—11.8 km， 主要分布在5—10 km， 均為淺源地震．

3.2 濾波方法

如上所述， SNR≥8 dB時拾取的P波初至精度較高， 當(dāng)信噪比太低時拾取的效果不甚滿意， 故本文對SNR<8 dB的原始數(shù)據(jù)進行濾波， 以提高拾取到時的精度．

濾波的主要步驟為： ① 給定一個濾波器組， 頻帶分別為1.5—3.6， 3.6—8.3， 8.3—10， 10—15， 15—20 Hz； ② 設(shè)定信噪比閾值為8 dB， 計算原始信號的信噪比； ③ 若原始信噪比大于閾值， 則認為原始信噪比較高， 無需濾波直接利用原始數(shù)據(jù)進行處理(Leach Jretal， 1999)； ④ 若原始信噪比小于閾值， 則對①中各頻帶分別進行濾波， 計算濾波后各頻帶對應(yīng)的信噪比， 確定信噪比最大時所對應(yīng)的濾波頻帶[f1，f2]．

鑒于有限沖激響應(yīng)(finite impulse response, 簡寫為FIR)數(shù)字濾波器的穩(wěn)定性好， 易實現(xiàn)嚴格的線性相位， 信號處理后不會產(chǎn)生相位畸變， 應(yīng)用范圍甚廣(萬永革， 2012)， 故本文用FIR數(shù)字濾波器(采用等波紋最佳逼近設(shè)計線性相位的Remez算法， 是目前最優(yōu)的FIR設(shè)計算法)對[f1，f2]頻段內(nèi)的原始數(shù)據(jù)(722個地震記錄)進行濾波， 本文將該方法稱為FIR最佳頻帶濾波方法．

圖5a給出了濾波前后原始記錄的信噪比變化． 可以看出： 濾波前信噪比的變化范圍為-5.82—37.54 dB， 均值為(9.12±6.76) dB， SNR>8 dB的占50.28%； 濾波后信噪比的變化范圍為0.9—38.04 dB， 均值為(14.15±5.03) dB， SNR>8 dB的高達96.12%． 圖5b給出了濾波前后P波震相初至識別差值的變化. 可以看出， 濾波前自動拾取與人工拾取的P波震相初至識別差值ΔT′分布較分散， 有些差值過大， 而濾波后的差值分布較為均勻， 基本都在0值附近波動． 由此可見， 經(jīng)FIR最佳頻帶濾波后不僅增強了信噪比， 而且有效提高了P波識別精度．

圖5 濾波前、 后信噪比(a)及P波初至識別差值ΔT′(b)的對比

3.3 P波初至識別差值、 信噪比及初至清晰度的關(guān)系

利用FIR最佳頻帶濾波方法對722個原始地震記錄進行自適應(yīng)濾波(SNR≥8 dB時不濾波)， 然后用小波包-峰度AIC方法自動拾取P波初至到時， 自動拾取與人工拾取的到時差結(jié)果見圖6a． 可以看出： P波初至識別差值為-1.4—1.2 s， 平均絕對差值為(0.234±0.271) s， 其中絕對差值在0.3 s以內(nèi)的占75.07%； 隨著信噪比的增大， P波到時拾取的精度在一定程度上得到了提高， 說明P波拾取的精度與信噪比有一定關(guān)系， 但這種關(guān)系并非是信噪比越高P波初至識別差值越小的線性關(guān)系．

為了進一步探求P波初至識別差值的影響因素， 本文將722個地震事件按初動清晰度分為初至清晰(初動突出， 肉眼能明顯識別)和初至不清晰(初動不突出， 肉眼難以識別)兩種情況： 初至清晰時清晰度設(shè)為1， 共322個地震記錄； 初至不清晰時設(shè)為-1， 共400個地震記錄． 圖6b給出了信噪比與初至清晰度的關(guān)系， 可以看出： 初至清晰時信噪比變化范圍為8.1—38.04 dB， 平均信噪比為(15.73±5.24) dB； 初至不清晰時信噪比變化范圍為0.9—30.66 dB， 均值為(12.88±4.48) dB． 由此可見， 初至清晰時對應(yīng)的信噪比往往較高， 而信噪比高時初至卻不一定清晰．

圖6c給出了P波初至識別差值與初至清晰度的關(guān)系． 可以看出： 當(dāng)初至清晰時， P波初至到時差的變化范圍為-0.27—0.46 s， 平均絕對差值為0.077 s， 標(biāo)準(zhǔn)差為0.075 s， 平均絕對差值在0.2 s和0.3 s以內(nèi)的分別占94.41%和98.45%； 當(dāng)初至不清晰時， P波到時差變化范圍為-1.4—1.2 s， 平均絕對差值為0.36 s， 標(biāo)準(zhǔn)差為0.303 s， 平均絕對差值在0.2 s, 0.3 s和0.5 s以內(nèi)的分別占38%， 56.25%和74%． 由此可見， 初至清晰度對P波初至識別差值的影響較大， 且初至越清晰， P波到時拾取的精度越高．

圖6 P波初至識別差值ΔT′、 信噪比及初至清晰度的關(guān)系(a) 初至差與信噪比的關(guān)系； (b) 信噪比與初至清晰度的關(guān)系； (c) 初至差與初至清晰度的關(guān)系

通過上述分析可知， P波到時自動拾取的精度與信噪比有一定關(guān)系， 但并非呈線性關(guān)系， 相對信噪比對其的影響， 初至清晰度對其的影響更大， 初至清晰時到時拾取的精度明顯比初至不清晰時高； 同時， 初至清晰時信噪比也較高， 但信噪比高時初至卻并不一定清晰． 因此， 可將原始記錄分為初至清晰和初至不清晰兩種情況， 以檢驗本文小波包-峰度AIC方法在不同實際情況下自動拾取P波初至到時的效果．

3.4 3種自動識別方法在實際地震中的應(yīng)用對比

為檢驗本文小波包-峰度AIC方法在近震P波震相自動識別中的應(yīng)用效果及識別精度， 將其與STA/LTA方法、 峰度AIC方法這兩種常用的震相自動識別方法進行對比分析．

圖7和圖8分別給出了初至清晰情況下3種方法自動識別與人機交互識別的P波到時結(jié)果及其差值分布. 可以看出： 本文小波包-峰度AIC方法所得震相識別結(jié)果與人機交互識別結(jié)果的平均絕對差值最小， 為(0.077±0.075) s， 到時差為-0.27—0.46 s， 322個單臺記錄的處理結(jié)果中， 絕對差值在0.1 s， 0.2 s和0.3 s以內(nèi)的分別占73.6%， 94.41%和 98.45%； 峰度AIC方法所得震相識別結(jié)果與人機交互所得到時的差值為-0.86—0.7 s， 平均絕對差值為(0.113±0.135) s， 絕對差值在0.1 s， 0.2 s和0.3 s以內(nèi)的分別占65.84%， 86.96%和93.17%； STA/LTA方法所得震相識別結(jié)果與人機交互拾取到時的差值為-0.27—0.98 s， 平均絕對差值為(0.16±0.184) s， 絕對差值在0.1 s， 0.2 s和0.3 s以內(nèi)的分別占55.28%， 77.33%和85.09%． 這表明， 在初至清晰的情況下， 3種方法自動識別近震P波震相的效果均較好， 其中本文提出的小波包-峰度AIC方法拾取的P波精度最高， 峰度AIC方法次之， STA/LTA方法精度相對低一些．

圖7 初至清晰時P波震相自動識別與人機交互識別結(jié)果對比圖中數(shù)值分別為3種方法得到的平均絕對差值. (a) STA/LTA方法；(b) 峰度AIC方法； (c) 小波包-峰度AIC方法

圖8 初至清晰時P波震相自動識別與人機交互識別差值分布 (a) STA/LTA方法； (b) 峰度AIC方法； (c) 小波包-峰度AIC方法 Fig.8 Error distribution of the P-wave onset time between automatic recognition and man-machine interaction recognition when first break is clear(a) STA/LTA method; (b) Kurtosis-AIC method; (c) Wavelet packet and Kurtosis-AIC method

圖9和圖10分別為初至不清晰的情況下3種方法自動識別與人機交互識別的P波到時結(jié)果及差值分布. 可以看出： 本文方法所得震相識別結(jié)果與人機交互識別結(jié)果的平均絕對差值最小， 為(0.36±0.303) s， 到時差值為-1.4—1.2 s， 400個單臺記錄的處理結(jié)果中， 絕對差值在0.2 s， 0.3 s和0.5 s以內(nèi)的分別占38%， 56.25%和74%； 峰度AIC方法所得震相識別結(jié)果與人機交互拾取的到時差值為-1.34—1.37 s， 平均絕對差值為(0.373±0.314) s， 絕對差值在0.2 s， 0.3 s和0.5 s以內(nèi)的分別占37.75%， 53.75%和71.75%； STA/LTA方法所得震相識別結(jié)果與人機交互拾取的到時差值為-1.25—1.49 s， 平均絕對差值為(0.388±0.341) s， 絕對差值在0.2 s， 0.3 s和0.5 s以內(nèi)的分別占37.25%， 53.75%和71.5%． 這表明，在初至不清晰的情況下， 3種方法自動識別近震P波震相的效果均不如初至清晰時好， 其中本文小波包-峰度AIC方法拾取的P波精度相對高一些， 峰度AIC方法和STA/LTA方法的精度相對較低， 反映了在人工震相識別中同樣面臨的當(dāng)P波初至不清晰時震相識別精度不高這一難題．

圖9 初至不清晰時P波震相自動識別與人機交互識別結(jié)果對比圖中數(shù)值分別為3種方法得到的平均絕對差值. (a) STA/LTA方法；(b) 峰度AIC方法； (c) 小波包-峰度AIC方法

Fig.9 Comparison of the P-wave onset time result by automatic recognition with that by man-machine interaction recognition on the condition that the first break is unclear The average absolute errors of the three methods are also given. (a) STA/LTA method;(b) Kurtosis-AIC method; (c) Wavelet packet and Kurtosis-AIC method

圖10 初至不清晰時P波震相自動識別與人機交互識別差值分布(a) STA/LTA方法； (b) 峰度AIC方法； (c) 小波包-峰度AIC方法

4 討論與結(jié)論

本文提出了一種基于小波包和峰度AIC的P波震相自動識別方法. 在加權(quán)STA/LTA方法粗略拾取P波到時的基礎(chǔ)上， 利用該方法進一步拾取精細的P波到時并將其應(yīng)用于模擬事件的理論地震記錄， 模擬事件參照云南地區(qū)實際震例設(shè)計． 對理論地震記錄加入不同信噪比的高斯白噪聲和實際地震噪聲， 以由射線追蹤技術(shù)得到的到時為標(biāo)準(zhǔn)， 對比了STA/LTA方法、 峰度AIC方法和本文的小波包-峰度AIC方法自動識別P波的效果(表1和表2)． 結(jié)果表明， 本文方法具有更強的抗噪能力， P波識別的精度更高， 當(dāng)加入實際地震噪聲的SNR≥8 dB時， 其自動拾取與理論計算的P波差值在0.25 s以內(nèi)．

將本文方法應(yīng)用到云南地區(qū)722個單臺垂直向?qū)嶋H地震記錄中， 濾波前自動拾取與人機交互拾取的P波差值部分過大(圖5b)， 這是由于原始記錄信噪比低或初至不清晰而造成的； 通過本文FIR最佳頻帶方法自適應(yīng)濾波后， 大大提高了原始記錄的信噪比和P波拾取的精度． P波初至識別差值、 信噪比和初至清晰度關(guān)系(圖6)的分析結(jié)果表明： 初至清晰時信噪比較高， 但信噪比高時初至卻并不一定清晰； P波自動拾取的精度與信噪比有一定關(guān)系， 但并非是信噪比越高P波初至識別差值越小的線性關(guān)系， 相對信噪比對其的影響， 初至清晰度對其的影響更大， 且初至清晰時到時拾取的精度明顯比初至不清晰時高．

將地震記錄分為初至清晰和初至不清晰兩種情況， 以人工拾取的震相到時為標(biāo)準(zhǔn)， 分別對比STA/LTA方法、 峰度AIC方法和本文方法的效果， 結(jié)果表明本文的小波包-峰度AIC方法拾取P波初至的精度更高. 當(dāng)初至清晰時， 經(jīng)FIR最佳頻帶方法濾波后信噪比均大于8 dB， 本文方法自動識別的P波到時與人工識別結(jié)果的平均絕對差值為(0.077±0.075) s， 絕對差值在0.1 s和0.2 s以內(nèi)的分別占73.6%和94.41%， 精度非常高．

鑒于本文研究所選用的地震均為近震事件， 下一步將把本文方法運用到遠震事件中進一步檢驗該方法的效果， 當(dāng)初至不清晰時P波識別的精度有待于進一步提高．

感謝德國地學(xué)研究中心汪榮江博士提供合成理論地震圖的計算程序， 感謝Ludger Küperkoch教授提供幫助和有益探討．

常旭， 劉伊克. 2002. 地震記錄的廣義分維及其應(yīng)用[J]. 地球物理學(xué)報， 45(6): 839-846.

Chang X， Liu Y K. 2002. The generalized fractal dimension of seismic records and its applications[J].ChineseJournalofGeophysics， 45(6): 839--846 (in Chinese).

劉晗， 張建中. 2014. 微震信號自動檢測的STA/LTA算法及其改進分析[J]. 地球物理學(xué)進展， 29(4): 1708--1714.

Liu H， Zhang J Z. 2014. STA/LTA algorithm analysis and improvement of microseismic signal automatic detection[J].ProgressinGeophysics， 29(4): 1708--1714 (in Chinese).

劉希強， 周蕙蘭， 鄭治真， 沈萍， 楊選輝， 馬延路. 1998. 基于小波包變換的弱震相識別方法[J]. 地震學(xué)報， 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].ActaSeismologicaSinica， 11(4): 431--439.

劉希強， 周蕙蘭， 曹文海， 李紅， 李永紅， 季愛東. 2002. 高斯線調(diào)頻小波變換及其在地震震相識別中的應(yīng)用[J]. 地震學(xué)報， 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 phases identification[J].ActaSeismologicaSinica， 24(6): 607--616 (in Chinese).

劉希強， 周彥文， 曲均浩， 石玉燕， 李鉑. 2009. 應(yīng)用單臺垂向記錄進行區(qū)域地震事件實時檢測和直達P波初動自動識別[J]. 地震學(xué)報， 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].ActaSeismologicaSinica， 31(3): 260--271 (in Chinese).

毛燕， 崔建文， 鄭定昌， 李正光， 盧吉高. 2011. 地震記錄的P波自動撿拾[J]. 地震研究， 34(1): 47--51.

Mao Y， Cui J W， Zheng D C， Li Z G， Lu J G. 2011. Automatic P-wave detection of the earthquake recordings[J].JournalofSeismologicalResearch， 34(1): 47--51 (in Chinese).

曲保安， 劉希強， 蔡寅， 范曉勇， 林秀娜， 于慶民， 趙瑞， 李鉑， 周彥文. 2014. 近震S波震相實時自動識別方法研究[J]. 地震學(xué)報， 36(2): 184--192.

Qu B A， Liu X Q， Cai Y， Fan X Y， Lin X N， Yu Q M， Zhao R， Li B， Zhou Y W. 2014. Method for real-time automatic identification of S-phase: Application to local seismicity[J].ActaSeismologicaSinica， 36(2): 184--192 (in Chinese).

孫進忠， 趙鴻儒， 彭一民. 1985. 全波模型的震相分析[J]. 石油地球物理勘探， 20(4): 352--362.

Sun J Z， Zhao H R， Peng Y M. 1985. Full wave model phase analysis (FPA)[J].OilGeophysicalProspecting， 20(4): 352--362 (in Chinese).

萬永革. 2012. 數(shù)字信號處理的MATLAB實現(xiàn)[M]. 第2版. 北京: 科學(xué)出版社: 298--322.

Wan Y G. 2012.DigitalSignalProcessingBasedonMATLAB[M]. 2nd ed. Beijing: Science Press: 298--322 (in Chinese).

王海軍， 靳平， 劉貴忠， 王曉明. 2003. 區(qū)域震相初至估計[J]. 西北地震學(xué)報， 25(4): 298--303.

Wang H J， Jin P， Liu G Z， Wang X M. 2003. Accurate estimation for arrival time of seismic wave[J].NorthwesternSeismologicalJournal， 25(4): 298--303 (in Chinese).

王海軍， 劉貴忠， 范萬春. 2004. 廣義分維在地震信號初至檢測中的應(yīng)用[J]. 核電子學(xué)與探測技術(shù)， 24(6): 634--637.

Wang H J， Liu G Z， Fan W C. 2004. The application of generalized fractal in arrival time detection of seismograms[J].NuclearElectronicsandDetectionTechnology， 24(6): 634--637 (in Chinese).

王繼， 陳九輝， 劉啟元， 李順成， 郭飚. 2006. 流動地震臺陣觀測初至震相的自動檢測[J]. 地震學(xué)報， 28(1): 42--51.

Wang J， Chen J H， Liu Q Y， Li S C， Guo B. 2006. Automatic onset phase picking for portable seismic array observation[J].ActaSeismologicaSinica， 28(1): 42--51 (in Chinese).

王娟， 劉俊民， 范萬春. 2004. 神經(jīng)網(wǎng)絡(luò)在震相識別中的應(yīng)用[J]. 現(xiàn)代電子技術(shù)， 27(8): 35--37.

Wang J， Liu J M， Fan W C. 2004. Application of neural network in the discrimination of seismical signal[J].ModernElectronicsTechnique， 27(8): 35--37 (in Chinese).

王喜珍. 2004. 小波變換在地震數(shù)據(jù)壓縮和震相到時拾取中的應(yīng)用研究[D]. 北京: 中國地震局地球物理研究所: 87--90.

Wang X Z. 2004.StudyonApplicationofWaveletTransforminCompressingSeismicDataandPickingtheOnsetTimeofSeismicPhase[D]. Beijing: Institute of Geophysics， China Earthquake Administration: 87--90 (in Chinese).

武東坡. 2004. 震相識別的實時方法研究[D]. 哈爾濱: 中國地震局工程力學(xué)研究所: 11--22.

Wu D P. 2004.ResearchontheReal-TimeProcessingMethodofSeismicPhaseIdentification[D]. Harbin: Institute of Engineering Mechanics， China Earthquake Administration: 11--22 (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 recor-dings[J].JournalofSeismologicalResearch， 27(4): 308--313 (in Chinese).

趙愛華， 張中杰， 彭蘇萍. 2003. 復(fù)雜地質(zhì)模型轉(zhuǎn)換波快速射線追蹤方法[J]. 中國礦業(yè)大學(xué)學(xué)報， 32(5): 513--516.

Zhao A H， Zhang Z J， Peng S P. 2003. Fast ray tracing method for converted waves in complex media[J].JournalofChinaUniversityofMining&Technology， 32(5): 513--516 (in Chinese).

趙大鵬， 劉希強， 李紅， 周彥文. 2012. 峰度和AIC方法在區(qū)域地震事件和直達P波初動自動識別方面的應(yīng)用[J]. 地震研究， 35(2): 220--226.

Zhao D P， Liu X Q， Li H， Zhou Y W. 2012. Detection of regional seismic events by Kurtosis method and automatic identification of direct P-wave first motion by Kurtosis-AIC method[J].JournalofSeismologicalResearch， 35(2): 220--226 (in Chinese).

趙鴻儒， 孫進忠， 唐文榜. 1990. 全波震相分析的應(yīng)用[J]. 地球物理學(xué)報， 33(1): 54--63.

Zhao H R， Sun J Z， Tang W B. 1990. The application of full wave phase analysis[J].ChineseJournalofGeophysics， 33(1): 54--63 (in Chinese).

莊東海， 肖春燕， 顏永寧. 1994. 利用人工神經(jīng)網(wǎng)絡(luò)自動拾取地震記錄初至[J]. 石油地球物理學(xué)報， 29(5): 659--664.

Zhuang D H， Xiao C Y， Yan Y N. 1994. Seismic first arrival pickup using artificial neural network[J].OilGeophysicalProspecting， 29(5): 659--664 (in Chinese).

Allen R. 1978. Automatic earthquake recognition and timing from single trace[J].BullSeismolSocAm， 68(5): 1521--1532.

Allen R. 1982. Automatic phase pickers： Their present use and future prospects[J].BullSeismolSocAm， 72(6B): S225--S242.

Baer M， Kradolfer U. 1987. An automatic phase picker for local and teleseismic events[J].BullSeismolSocAm， 77(4): 1437--1445.

Bai C， Kennett B L N. 2000. Automatic phase-detection and identification by full use of a single three-component broadband seismogram[J].BullSeismolSocAm， 90(1): 187--198.

Boschetti F， Dentith M D， List R D. 1996. A fractal-based algorithm for detecting first arrivals on seismic traces[J].Geophysics， 61(4): 1095--1102.

Cichowicz A， Green R W E， van Zyl Brink A. 1988. Coda polarization properties of high-frequency microseismic events[J].BullSeismolSocAm， 78(3): 1297--1318.

Dai H C， MacBeth C. 1997. The application of back-propagation neural network to automatic picking seismic arrivals from single-component recordings[J].JGeophysRes， 102(B7): 15105--15113.

Jurkevics A. 1988. Polarization analysis of three-component array data[J].BullSeismolSocAm， 78(5): 1725--1743.

Küperkoch L， Meier T， Lee J， Friederich W， EGELADOS Working Group. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics[J].GeophysJInt， 181(2): 1159--1170.

Küperkoch L， Meier T， Brüstle A， Lee J， Friederich W， EGELADOS Working Group. 2012. Automated determination of S-phase arrival times using autoregressive prediction: Application to local and regional distances[J].GeophysJInt， 188(2): 687--702.

Leach Jr R R， Dowla F U， Schultz C A. 1999. Optimal filter parameters for low SNR seismograms as a function of station and event location[J].PhysEarthPlanetInt， 113(1/2/3/4): 213--226.

Maeda N. 1985. A method for reading and checking phase times in auto-processing system of seismic wave data[J].Zisin， 38(3): 365--379.

Saragiotis C D， Hadjileontiadis L J， Panas S M. 2002. PAI-S/K: A robust automatic seismic P phase arrival identification scheme[J].IEEETransGeosciRemoteSensing， 40(6): 1395--1404.

Satriano C， Zollo A， Rowe C. 2008. Iterative tomographic analysis based on automatic refined picking[J].GeophysProsp， 56(4): 467--475.

Sleeman R， van Eck T. 1999. Robust automatic P-phase picking: An on-line implementation in the analysis of broadband seismogram recordings[J].PhysEarthPlanetInt， 113(1/2/3/4): 265--275.

Stevenson P R. 1976. Microearthquakes at Flathead Lake， Montana: A study using automatic earthquake processing[J].BullSeismolSocAm， 66(1): 61--80.

Wang R J. 1999. A simple orthonormalization method for stable and efficient computation of Green’s functions[J].BullSeismolSocAm， 89(3): 733--741.

Wang R J， Wang H S. 2007. A fast converging and anti-aliasing algorithm for Green’s functions in terms of spherical or cylindrical harmonics[J].GeophysJInt， 170(1): 239--248.

Wickerhauser M V. 1992. Acoustic signal compression with wavelet packets[G]∥Wavelets:ATutorialinTheoryandApplication. San Diego， California: Academic Press: 679--700.

Yomogida K. 1994. Detection of anomalous seismic phases by the wavelet transform[J].GeophysJInt， 116(1): 119--130.

Zhang H J， Thurber C， Rowe C. 2003. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings[J].BullSeismolSocAm， 93(5): 1904--1912.

Zhao A H， Zhang Z J， Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing: Improvement in efficiency[J].JGeophysEng， 1(4): 245--251.

Automatic identification of P-phase based on wavelet packet and Kurtosis-AIC method

1)EarthquakeAdministrationofHunanProvince，Changsha410004，China2)InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China

Automatic identification of P-phase is of significance to the study on earthquake location， earthquake warning and structure of deep earth. Combining wavelet packet transform with Kurtosis-AIC (Akaike information criterion) technology， this paper puts forward a new synthetic method named wavelet packet and Kurtosis-AIC method for automatic recognition of first P-phase. Three scales of discrete wavelet packet transforms are applied to decompose and reconstructure the original recordings three seconds before and after the rough P-wave arrival time， which is picked up by weighted STA/LTA (short term average/long term average) method， then the Kurtosis-AIC values of the three-scale reconstruction signal are calculated respectively and superposed together， finally the minimum value of the superposed AIC curve is taken as the first P-wave arrival time. In order to test the new method， it is applied to theoretically synthetic seismograms and real seismic recording for automatic P-phase arrival time detection. Adding white Gaussian noise and real seismic noise to synthetic seismograms with different SNR， the optimal frequency band of adaptive FIR (finite impulse response) digital filtering is used to improve the SNR and P-wave recognition accuracy of the original signals. The results show that， with respect to the impact of SNR， the accuracy of P-wave identification is more affected by the clarity of first break; our method has greater noise immunity and higher P-wave recognition accuracy as compared to the weighted STA/LTA algorithm and Kurtosis-AIC method. When the first break of P-wave is clear， average absolute error of P-phase arrival time between automatic identification based on our method and manual identification is (0.077±0.075) seconds.

P-phase; FIR digital filtering; automatic identification of seismic phase; wavelet packet and Kurtosis-AIC method

國家自然科學(xué)基金(40974050， 41374098)資助.

2015-05-15收到初稿， 2015-08-17決定采用修改稿.

e-mail: ahzhao123@yahoo.com

10.11939/jass.2016.01.007

P315.3+1

A

田優(yōu)平， 趙愛華. 2016. 基于小波包和峰度赤池信息量準(zhǔn)則的P波震相自動識別方法. 地震學(xué)報, 38(1): 71--85. doi:10.11939/jass.2016.01.007.

Tian Y P, Zhao A H. 2016. Automatic identification of P-phase based on wavelet packet and Kurtosis-AIC method.ActaSeismologicaSinica, 38(1): 71--85. doi:10.11939/jass.2016.01.007.