田優(yōu)平 趙愛華

1) 中國(guó)長(zhǎng)沙410004湖南省地震局 2) 中國(guó)北京100081中國(guó)地震局地球物理研究所

?

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

1) 中國(guó)長(zhǎng)沙410004湖南省地震局 2) 中國(guó)北京100081中國(guó)地震局地球物理研究所

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

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

引言

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

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

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

多特征法主要是基于地震波傳播的整體特征或多個(gè)特征來(lái)識(shí)別震相， 如全波震相分析法(孫進(jìn)忠等， 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)．

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

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

1 方法

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

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

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

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

(1)

(2)

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

(3)

來(lái)計(jì)算. 式中，Y(i)為i時(shí)刻的地震記錄值．

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

(4)

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

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

1.2 峰度AIC方法

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

(5)

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

(6)

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

(7)

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

1.3 一維小波包多尺度分解

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

(8)

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

離散小波包變換的表達(dá)式可寫為

(9)

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

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

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

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

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

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

2 理論應(yīng)用

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

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

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

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

SNR/dB自動(dòng)拾取的P波初至到時(shí)/sSTA/LTA方法峰度AIC方法小波包?峰度AIC方法自動(dòng)與理論到時(shí)差Δ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方法在合成理論地震記錄中加入不同噪聲時(shí)表現(xiàn)為抗噪能力更強(qiáng)且拾取的P波到時(shí)精度更高； 當(dāng)加入的實(shí)際地震噪聲SNR≥8 dB時(shí)， ΔT≤0.25 s， 效果要明顯優(yōu)于其它兩種方法．

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

3 實(shí)際應(yīng)用

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

3.1 數(shù)據(jù)

本文對(duì)研究區(qū)(97.5°E—107°E， 20°N—30°N)2009年2月—2011年6月107次ML1.4—4.9地震的722個(gè)單臺(tái)垂直向記錄進(jìn)行研究和分析， 涉及的臺(tái)站共計(jì)52個(gè)． 選取單臺(tái)記錄的震中距為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時(shí)拾取的P波初至精度較高， 當(dāng)信噪比太低時(shí)拾取的效果不甚滿意， 故本文對(duì)SNR<8 dB的原始數(shù)據(jù)進(jìn)行濾波， 以提高拾取到時(shí)的精度．

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

鑒于有限沖激響應(yīng)(finite impulse response, 簡(jiǎn)寫為FIR)數(shù)字濾波器的穩(wěn)定性好， 易實(shí)現(xiàn)嚴(yán)格的線性相位， 信號(hào)處理后不會(huì)產(chǎn)生相位畸變， 應(yīng)用范圍甚廣(萬(wàn)永革， 2012)， 故本文用FIR數(shù)字濾波器(采用等波紋最佳逼近設(shè)計(jì)線性相位的Remez算法， 是目前最優(yōu)的FIR設(shè)計(jì)算法)對(duì)[f1，f2]頻段內(nèi)的原始數(shù)據(jù)(722個(gè)地震記錄)進(jìn)行濾波， 本文將該方法稱為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的高達(dá)96.12%． 圖5b給出了濾波前后P波震相初至識(shí)別差值的變化. 可以看出， 濾波前自動(dòng)拾取與人工拾取的P波震相初至識(shí)別差值ΔT′分布較分散， 有些差值過大， 而濾波后的差值分布較為均勻， 基本都在0值附近波動(dòng)． 由此可見， 經(jīng)FIR最佳頻帶濾波后不僅增強(qiáng)了信噪比， 而且有效提高了P波識(shí)別精度．

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

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

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

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

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

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

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

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

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

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

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

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

圖9 初至不清晰時(shí)P波震相自動(dòng)識(shí)別與人機(jī)交互識(shí)別結(jié)果對(duì)比圖中數(shù)值分別為3種方法得到的平均絕對(duì)差值. (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 初至不清晰時(shí)P波震相自動(dòng)識(shí)別與人機(jī)交互識(shí)別差值分布(a) STA/LTA方法； (b) 峰度AIC方法； (c) 小波包-峰度AIC方法

4 討論與結(jié)論

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

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

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

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

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

常旭， 劉伊克. 2002. 地震記錄的廣義分維及其應(yīng)用[J]. 地球物理學(xué)報(bào)， 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. 微震信號(hào)自動(dòng)檢測(cè)的STA/LTA算法及其改進(jìn)分析[J]. 地球物理學(xué)進(jìn)展， 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).

劉希強(qiáng)， 周蕙蘭， 鄭治真， 沈萍， 楊選輝， 馬延路. 1998. 基于小波包變換的弱震相識(shí)別方法[J]. 地震學(xué)報(bào)， 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.

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

劉希強(qiáng)， 周彥文， 曲均浩， 石玉燕， 李鉑. 2009. 應(yīng)用單臺(tái)垂向記錄進(jìn)行區(qū)域地震事件實(shí)時(shí)檢測(cè)和直達(dá)P波初動(dòng)自動(dòng)識(shí)別[J]. 地震學(xué)報(bào)， 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波自動(dòng)撿拾[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).

曲保安， 劉希強(qiáng)， 蔡寅， 范曉勇， 林秀娜， 于慶民， 趙瑞， 李鉑， 周彥文. 2014. 近震S波震相實(shí)時(shí)自動(dòng)識(shí)別方法研究[J]. 地震學(xué)報(bào)， 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).

孫進(jìn)忠， 趙鴻儒， 彭一民. 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).

萬(wàn)永革. 2012. 數(shù)字信號(hào)處理的MATLAB實(shí)現(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ì)[J]. 西北地震學(xué)報(bào)， 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).

王海軍， 劉貴忠， 范萬(wàn)春. 2004. 廣義分維在地震信號(hào)初至檢測(cè)中的應(yīng)用[J]. 核電子學(xué)與探測(cè)技術(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. 流動(dòng)地震臺(tái)陣觀測(cè)初至震相的自動(dòng)檢測(cè)[J]. 地震學(xué)報(bào)， 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).

王娟， 劉俊民， 范萬(wàn)春. 2004. 神經(jīng)網(wǎng)絡(luò)在震相識(shí)別中的應(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ù)壓縮和震相到時(shí)拾取中的應(yīng)用研究[D]. 北京: 中國(guó)地震局地球物理研究所: 87--90.

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

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

Wu D P. 2004.ResearchontheReal-TimeProcessingMethodofSeismicPhaseIdentification[D]. Harbin: Institute of Engineering Mechanics， China Earthquake Administration: 11--22 (in Chinese).

楊配新， 鄧存華， 劉希強(qiáng)， 任勇， 顏其中. 2004. 數(shù)字化地震記錄震相自動(dòng)識(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]. 中國(guó)礦業(yè)大學(xué)學(xué)報(bào)， 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).

趙大鵬， 劉希強(qiáng)， 李紅， 周彥文. 2012. 峰度和AIC方法在區(qū)域地震事件和直達(dá)P波初動(dòng)自動(dòng)識(shí)別方面的應(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).

趙鴻儒， 孫進(jìn)忠， 唐文榜. 1990. 全波震相分析的應(yīng)用[J]. 地球物理學(xué)報(bào)， 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ò)自動(dòng)拾取地震記錄初至[J]. 石油地球物理學(xué)報(bào)， 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

國(guó)家自然科學(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波震相自動(dòng)識(shí)別方法. 地震學(xué)報(bào), 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.