姚彥吉，柳林濤，盛敏漢，許厚澤

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430077 2 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京 100049

0 引言

地震P波、S波到時(shí)的拾取準(zhǔn)確率直接影響地震震源位置的反演精度，因此地震P波、S波到時(shí)拾取是地震數(shù)據(jù)處理的一項(xiàng)重要的基礎(chǔ)工作(何先龍等，2016；于子葉等，2018).在地震發(fā)生后需要快速實(shí)現(xiàn)P波、S波到時(shí)的精確拾取，這對分析地震發(fā)生規(guī)律、地震應(yīng)急、以及后續(xù)防震減災(zāi)工作尤其關(guān)鍵(張小紅等，2012).早期依靠人工拾取地震P波、S波到時(shí)，這種方法精度高，但是效率低、主觀性強(qiáng)、一致性差(陳金煥等，2015；趙明等，2019).隨著地震儀在全球范圍的增設(shè)以及區(qū)域密集臺(tái)陣的增加，波形數(shù)據(jù)量急劇增長，收集的數(shù)據(jù)質(zhì)量參差不齊，海量數(shù)據(jù)與強(qiáng)噪聲環(huán)境的干擾使得震相拾取的精度和速度需求在不斷的提高，研究者從時(shí)間域、頻率域與時(shí)頻域結(jié)合等方面提出了多種波形到時(shí)自動(dòng)拾取的方法(王彩霞等，2013；趙明等，2019).

典型的時(shí)間域震相拾取方法有長短時(shí)間平均(Short-Term Average and Long-Term Average ratio，STA/LTA)方法(Allen,1978,1982;Baer and Kradolfer,1987).根據(jù)波形在振幅、頻率等方面的不同特征建立特征方程，使比值靈敏地反映地震信號(hào)的到達(dá).但是在信噪比低的場景，應(yīng)對被噪聲干擾的信號(hào)時(shí)，傳統(tǒng)方法往往適應(yīng)性差且拾取精度難以達(dá)到預(yù)期(Han et al.,2009).為此，研究者引入多種特征函數(shù)(Allen,1978,1982;Baer and kradolfer,1987;Earle and Shearer,1994;高淑芳等，2008；付繼華等，2019)增強(qiáng)信號(hào)到達(dá)時(shí)信號(hào)幅值等特征的變化，以利于震相的自動(dòng)準(zhǔn)確拾取.AIC(Akaike Information Criteria)方法(Akaike,1974)利用赤池信息準(zhǔn)則尋找AIC函數(shù)局部最小值作為震相到時(shí)，能夠抑制強(qiáng)背景噪聲的干擾.Maeda(1985)提出了VAR-AIC方法，直接根據(jù)地震記錄得到AIC函數(shù)，AIC函數(shù)的局部最小值對應(yīng)的時(shí)間為地震震相到時(shí).為了提高震相的拾取精度，STA/LTA與AIC方法被不斷改進(jìn)(劉勁松等，2013；何先龍等，2016)；楊旭等(2019)使用基于參數(shù)優(yōu)化的頻帶-帶寬拾取算法、AICD拾取算法和峰度拾取算法有效拾取地震P、S波到時(shí).另外，神經(jīng)網(wǎng)絡(luò)方法也在拾取震相中得到了有效應(yīng)用(Zhou et al.,2019;趙明等，2019；李健等，2020；Chai et al.,2020).Perol等(2018)應(yīng)用CNN(Convolutional Neural Networks)對美國俄克拉荷馬州地區(qū)的連續(xù)記錄波形進(jìn)行識(shí)別與定位；Zhu和Beroza(2019)使用U-Net識(shí)別美國南加州地區(qū)P波和S波初至，利用概率分布確定初至波的到時(shí)，取得較好的效果.趙明等(2019)基于U形卷積神經(jīng)網(wǎng)絡(luò)的震相識(shí)別與到時(shí)拾取，參照人工拾取結(jié)果計(jì)算均方根誤差Pg為0.41 s，Sg為0.54 s.神經(jīng)網(wǎng)絡(luò)方法實(shí)現(xiàn)模型有識(shí)別信號(hào)的能力，需要搭建復(fù)雜的網(wǎng)絡(luò)模型及準(zhǔn)備大量已知到時(shí)的訓(xùn)練樣本.在時(shí)間域中拾取地震波到時(shí)容易受噪聲的干擾，導(dǎo)致到時(shí)的拾取存在誤差；將地震信號(hào)從時(shí)間域轉(zhuǎn)換到頻率域在一定程度上可抑制噪聲的干擾.

不同的地震記錄，由于噪聲的干擾和介質(zhì)的復(fù)雜性，導(dǎo)致地震信號(hào)在時(shí)間域、頻率域形態(tài)上均表現(xiàn)出復(fù)雜的特征(Johnston et al.,1979).地震信號(hào)往往是一種非平穩(wěn)信號(hào)，其頻率的局部化特性對震相、弱震相(劉希強(qiáng)等，1998)的準(zhǔn)確拾取十分重要.常用的頻率域震相拾取方法有瞬時(shí)頻率法(Bai and Kennett,2001)、小波變換的主成分分析法(Anant and Dowla,1997;劉希強(qiáng)等，2000)、小波包變換的時(shí)頻分析法(劉希強(qiáng)等，1998)，小波包方法使用小波變換系數(shù)刻畫信號(hào)的突變和漸變，用小波特征進(jìn)行到時(shí)拾取.但是，在強(qiáng)噪聲環(huán)境中不同的地震記錄特征存在較大的差異，單一應(yīng)用時(shí)間域、頻率域方法提取的固定特征較難保證所有地震記錄的P波、S波到時(shí)的精確拾取.

時(shí)頻分析方法具有較好的抗噪能力，主要包括短時(shí)傅里葉變換(Gabor,1946)、小波變換(Grossmann et al.,1985)、S變換(Stockwell et al.,1996)等.Shensa(1977)用快速傅里葉變換的能量譜密度檢測震相.劉代志等(2005)采用小波包分解檢測地震波信號(hào)初至點(diǎn).Anant和Dowla(1997)、Mousavi等(2016)通過小波變換系數(shù)中的偏振和振幅信息拾取P波和S波到時(shí).Pinnegar和Mansinha(2003)、張小紅等(2012)根據(jù)S變換后時(shí)頻圖上振幅能量的突變來判斷P波的到時(shí).S變換具有較好的時(shí)頻分辨率，相比時(shí)域方法具有較強(qiáng)的抗噪能力(王彩霞等，2013).但是，定義在L2(R)空間中的小波變換存在頻率偏移現(xiàn)象(Liu and Hsu,2012；蘇曉慶，2014)，S變換中存在高頻部分頻率分辨率不準(zhǔn)確(武國寧等，2011)的問題，會(huì)導(dǎo)致震相到時(shí)的拾取產(chǎn)生偏差.柳林濤和許厚澤(2009)提出的標(biāo)準(zhǔn)時(shí)頻變換方法(Normal Time-Frequency Transform，NTFT)，是傅里葉變換、小波變換與S變換的繼承與發(fā)展.NTFT改進(jìn)了小波變換、S變換的問題，能無偏地表示信號(hào)的即時(shí)頻率、即時(shí)振幅和即時(shí)相位，最大限度的保留原始信息，具有較好的時(shí)頻分辨率(柳林濤等，2016).

考慮震相拾取時(shí)間域方法振幅響應(yīng)靈敏特性與NTFT變換無偏特性的優(yōu)點(diǎn)，本文提出強(qiáng)噪聲環(huán)境下地震P波、S波到時(shí)拾取的時(shí)頻域綜合方法.通過理論分析、合成測試以及實(shí)際數(shù)據(jù)分析，證實(shí)了本文所提出方法的有效性.

1 方法原理

1.1 標(biāo)準(zhǔn)時(shí)頻變換

NTFT時(shí)頻分析方法(柳林濤和許厚澤，2009；柳林濤等，2016)滿足量綱守恒，信號(hào)的頻率、振幅和相位隨時(shí)間的變化而變化，消除了可能存在的偽頻率和頻率偏移，提高時(shí)頻譜的分辨率.對于一個(gè)時(shí)間信號(hào)f(t)∈C，其標(biāo)準(zhǔn)時(shí)頻變換可以表達(dá)為：

(1)

(2)

式中，w(t)表示窗函數(shù).在標(biāo)準(zhǔn)時(shí)頻變換中，核函數(shù)的傅里葉變換式(3)滿足條件(4)，

(3)

(4)

(5)

式中w(t)表示高斯窗，σ表示高斯窗窗寬參數(shù).本文構(gòu)造的核函數(shù)如式(6)所示：

(6)

構(gòu)建的標(biāo)準(zhǔn)小波變換表達(dá)為：

(7)

對一個(gè)調(diào)和信號(hào)(8)式h(t)進(jìn)行標(biāo)準(zhǔn)時(shí)頻變換時(shí)，可滿足(9)式的即時(shí)無偏特性.

h(t)=Aexp(j(βt+φ)),

(8)

式中，A表示振幅，β表示頻率，φ表示初始相位；βt+φ表示即時(shí)相位.

(9)

令t-τ=x，則t=x+τ，

(10)

(11)

所以Ψh(τ,β)=Aexp(j(φ+βτ))=h(τ)，一個(gè)調(diào)和信號(hào)的標(biāo)準(zhǔn)時(shí)頻譜無偏的描述原始信號(hào)的特征.

標(biāo)準(zhǔn)時(shí)頻變換是小波變換與S變換的繼承與發(fā)展.如圖1所示，將一個(gè)波形進(jìn)行不同的時(shí)頻分析對比，圖1a中紅點(diǎn)對應(yīng)時(shí)刻的信號(hào)成分分別進(jìn)行小波變換、S變換與NTFT后與圖1b,1c,1d中黃色虛線時(shí)刻對應(yīng).文中小波變換采用的小波基函數(shù)為莫雷特小波，小波變換與S變換均需要先驗(yàn)信息調(diào)試參數(shù)達(dá)到最佳分辨率，時(shí)頻譜中均存在右偏現(xiàn)象，且高頻部分時(shí)頻分辨率較低.標(biāo)準(zhǔn)時(shí)頻變換無需先驗(yàn)信息，從數(shù)據(jù)本身出發(fā)且保證較高的時(shí)頻分辨率.從圖中波形最大振幅時(shí)刻進(jìn)行判斷，在莫雷特連續(xù)小波變換中存在0.05 s延遲，S變換中存在0.1 s延遲，標(biāo)準(zhǔn)時(shí)頻變換能夠無偏地描述原始信號(hào)的即時(shí)特性.標(biāo)準(zhǔn)時(shí)頻變換的無偏性減少了NTFT方法拾取到時(shí)中不必要的偏差，為準(zhǔn)確拾取到時(shí)提供了有利的基礎(chǔ)條件.

圖1 地震信號(hào)及其時(shí)頻譜分析對比(a)地震信號(hào)；(b)莫雷特連續(xù)小波變換；(c)S變換；(d)標(biāo)準(zhǔn)時(shí)頻變換.圖(a)中紅色點(diǎn)為地震信號(hào)振幅最大的時(shí)刻，時(shí)頻譜中黃色虛線為地震信號(hào)紅點(diǎn)變換后對應(yīng)的時(shí)刻.Fig.1 Comparison of seismic signal with its time-frequency spectrum analysis(a)Seismic signal;(b)Morlet continuous wavelet transform;(c)S transform;(d)Normal time-frequency transform.In figure (a),the red dot is the moment when the seismic signal amplitude is the largest,and the yellow dotted line in the spectrum is the moment when the seismic signal is transformed and corresponds to the red dot.

1.2 基于NTFT的STA/LTA方法

NTFT變換可確定地震信號(hào)的即時(shí)信息，結(jié)合NTFT時(shí)頻譜特征對STA/LTA方法的特征函數(shù)進(jìn)行改進(jìn)，通過NTFT確定的即時(shí)頻率對地震信號(hào)發(fā)生頻段進(jìn)行約束，選取地震信號(hào)發(fā)生的即時(shí)頻段[ωp,ωq]進(jìn)行震相拾取，即時(shí)頻段的判別過程參見(姚彥吉等，2020)，頻段約束定義如下：

(12)

(13)

(14)

因此，文中改進(jìn)特征函數(shù)表示為

(15)

(16)

式中，STAi為短窗口的特征函數(shù)均值，LTAi為長窗口的特征函數(shù)均值，CF為特征函數(shù)，S為短窗口寬度，L為長窗口寬度.

公式(15)特征函數(shù)能同時(shí)響應(yīng)地震信號(hào)即時(shí)振幅和即時(shí)頻率的變化，通過NTFT即時(shí)頻率對特征函數(shù)的頻率范圍進(jìn)行約束，可增強(qiáng)地震信號(hào)即時(shí)振幅和即時(shí)頻率的變化特征.將(15)式代入(16)式中，分別在長、短窗內(nèi)計(jì)算地震信號(hào)特征函數(shù)的均值，并通過STA/LTA的值檢測P波和S波的到時(shí).

1.3 基于NTFT的AIC方法

(17)

(18)

或者表示為：

(19)

同理，根據(jù)標(biāo)準(zhǔn)時(shí)頻譜判斷信號(hào)發(fā)生的時(shí)間段，設(shè)置滑窗的時(shí)間范圍.

AIC函數(shù)定義為(劉勁松等，2013)：

AIC(k)=(k-m+1)log{var(G([m,k]))}

+(n-k-1)log{var(G([k+1,n]))},

(20)

G([m,k])(k=m，m+1，…，n)表示被時(shí)間段[m,n]與NTFT即時(shí)頻率段[ωp,ωq]約束的地震信號(hào)；var(G[m,k])和var(G[k+1,n])分別表示滑窗內(nèi)數(shù)據(jù)集合G[m,k]與G[k+1,n]的方差.

2 算例分析

2.1 合成波形測試

為了測試在強(qiáng)噪聲環(huán)境下本文方法的性能，我們利用了合成波形檢驗(yàn)方法的有效性.利用frequency-wavenumber方法(Zhu and Rivera,2002)，根據(jù)CRUST1.0模型(Laske et al.,2013)，合成了2019年2月23日四川榮縣M4.7地震在EMS臺(tái)的理論波形，信號(hào)采樣頻率為20 Hz.該地震的震源參數(shù)根據(jù)Yang等(2020)的反演結(jié)果得到.用f(t)代表一個(gè)理論波形，將f(t)與噪聲信號(hào)n(t)疊加，如式(21)，構(gòu)成疊加信號(hào)s(t)以模擬含有噪聲的實(shí)際觀測地震信號(hào).合成實(shí)驗(yàn)中強(qiáng)噪聲信號(hào)n(t)主要假設(shè)為對地震波形到時(shí)干擾較大的白噪聲、脈沖干擾以及強(qiáng)能量單頻噪聲三種；脈沖干擾會(huì)使地震記錄失真，嚴(yán)重影響數(shù)據(jù)處理結(jié)果的質(zhì)量(周寶峰，2012)；疊加的強(qiáng)能量單頻噪聲為實(shí)測地震記錄中截取獲得，強(qiáng)能量單頻噪聲在時(shí)間域與地震波形相似甚至淹沒地震信號(hào)(張晟瑞等，2018).

s(t)=f(t)+n(t).

(21)

另外，通過改變噪聲信號(hào)的功率pn，來研究三種不同背景噪聲下地震震相自動(dòng)拾取方法的實(shí)際性能.式(22)為所求疊加信號(hào)的SNR.

(22)

式中，ps為信號(hào)f(t)的功率，pn為噪聲n(t)的功率.

理論地震波形P波到時(shí)TP為32.10 s，S波到時(shí)TS為47.16 s.分別在白噪聲(0 dB)、脈沖干擾(2.5 dB)、強(qiáng)能量單頻噪聲(-6 dB)三種不同的背景噪聲下，采用STA/LTA方法、AIC方法與本文方法對疊加信號(hào)s(t)進(jìn)行地震震相自動(dòng)提取，并根據(jù)拾取結(jié)果分析不同背景噪聲下各方法的實(shí)際性能.STA/LTA的短窗寬度S設(shè)為1 s，長窗寬度L設(shè)為10 s，P波和S波到時(shí)自動(dòng)拾取的閾值設(shè)為2.5.為了更好地與STA/LTA方法比較，基于NTFT的STA/LTA方法也選擇同樣的長、短窗寬度.

在上述給定的條件下，不同類型的噪聲信號(hào)n(t)(白噪聲0 dB，脈沖2.5 dB)分別隨機(jī)產(chǎn)生20組，強(qiáng)能量單頻噪聲(-6 dB)從實(shí)測地震數(shù)據(jù)中截取20組分別與合成地震信號(hào)疊加，地震P波和S波自動(dòng)拾取結(jié)果如表1所示.表中誤拾率是指地震P波和S波到時(shí)拾取誤差大于0.1 s的事件所占比例；漏拾率指未拾取到地震P波或S波到時(shí)的事件所占比例.

不同背景噪聲下地震P波、S波拾取效果如表1所示，在0 dB的白噪聲環(huán)境下，STA/LTA方法P波出現(xiàn)了20%的誤拾率與10%的漏拾率，正確率為70%.VAR-AIC方法P波出現(xiàn)了5%的誤拾率，對于S波的拾取，VAR-AIC方法的誤差小于STA/LTA方法；基于NTFT的STA/LTA、AIC方法在白噪聲干擾的情況下能較好地自動(dòng)提取地震P波、S波到時(shí)，基于NTFT的STA/LTA方法P波拾取結(jié)果平均偏差和標(biāo)準(zhǔn)偏差均為0.02 s，S波拾取結(jié)果平均偏差和標(biāo)準(zhǔn)偏差均為0.01 s，優(yōu)于STA/LTA方法與VAR-AIC方法；由于部分合成波形P波初至較弱且被強(qiáng)噪聲干擾，導(dǎo)致P波拾取誤差偏大.在脈沖干擾(2.5 dB)的環(huán)境下，STA/LTA方法受脈沖信號(hào)干擾失效，脈沖信號(hào)所發(fā)生的時(shí)刻直接影響STA/LTA計(jì)算峰值出現(xiàn)的時(shí)刻；對于VAR-AIC方法而言，脈沖信號(hào)對到時(shí)拾取干擾嚴(yán)重，P波和S波的拾取誤差均大于0.1 s.而基于NTFT的STA/LTA的方法能有效拾取P波(平均偏差0.02 s，標(biāo)準(zhǔn)偏差0.02 s)、S波(平均偏差0.02 s，標(biāo)準(zhǔn)偏差0.03 s)，基于NTFT的AIC方法明顯改善了S波(平均偏差0.01 s，標(biāo)準(zhǔn)偏差0.01 s)拾取結(jié)果的偏差.在-6 dB的強(qiáng)能量單頻噪聲環(huán)境下，STA/LTA方法P波誤拾率為65%；VAR-AIC方法P波誤拾率為35%，S波漏拾率為10%.相比之下，在強(qiáng)噪聲干擾的復(fù)雜環(huán)境中，基于NTFT的STA/LTA方法與基于NTFT的AIC方法仍保持著較高的拾取準(zhǔn)確率.

表1 不同背景噪聲下地震震相自動(dòng)拾取統(tǒng)計(jì)表Table 1 Automatic picking results of seismic phase under the different background noise

圖2為在0 dB強(qiáng)噪聲環(huán)境下地震信號(hào)及其以上四種不同方法的提取結(jié)果，圖中紅色實(shí)線表示地震信號(hào)P波的理論到時(shí)，紅色虛線表示S波的理論到時(shí).STA/LTA方法如圖2e所示，求得的STA/LTA值在地震P波到時(shí)處并無顯著優(yōu)勢，背景噪聲干擾導(dǎo)致在P波到達(dá)時(shí)發(fā)生漏拾，在P波到達(dá)之后S波到達(dá)之前誤拾噪聲信號(hào)；基于NTFT的STA/LTA方法如圖2f所示，繼承了STA/LTA方法的靈敏性，同時(shí)噪聲得到了較好地抑制，振幅的響應(yīng)變化明顯增強(qiáng).VAR-AIC方法拾取震相如圖2g所示，該方法在強(qiáng)噪聲情況下拾取震相具有一定的魯棒性，但是在P波初至較弱且被強(qiáng)噪聲淹沒的情況下，P波到時(shí)處的AIC函數(shù)值變化較弱，導(dǎo)致P波到時(shí)拾取存在誤差；在相同的背景噪聲下，S波到時(shí)處AIC函數(shù)值變化特征明顯，相比P波而言，S波到時(shí)的拾取誤差較小.基于NTFT的AIC方法拾取震相如圖2h所示，除去了地震信號(hào)所在頻段外的噪聲干擾，震相到達(dá)時(shí)刻的AIC函數(shù)值變化特征顯著，有利于震相到時(shí)的準(zhǔn)確拾取.

圖3為脈沖干擾環(huán)境下一組地震信號(hào)及其以上四種不同方法的拾取結(jié)果.脈沖信號(hào)對震相拾取的準(zhǔn)確率影響較大，脈沖干擾產(chǎn)生的原因有記錄器故障、周圍環(huán)境的突發(fā)性變化和干擾、以及操作人員的過失等(周寶峰，2012).如圖3c為STA/LTA方法計(jì)算結(jié)果，脈沖信號(hào)發(fā)生的位置引起了較大的振幅響應(yīng)變化，很難從閾值的角度加以剔除，容易導(dǎo)致地震信號(hào)震相拾取發(fā)生錯(cuò)誤.脈沖信號(hào)也使得VAR-AIC方法有較高的誤拾率與漏拾率，如圖3e所示，AIC計(jì)算的局部最小值并非地震信號(hào)的到時(shí).脈沖信號(hào)相對于地震信號(hào)分布于高頻段，在STA/LTA方法與AIC方法中引入NTFT可選取地震信號(hào)所在頻段拾取到時(shí)，有效抑制脈沖對地震信號(hào)的干擾.

進(jìn)一步測試本文方法在不同強(qiáng)噪聲干擾下的穩(wěn)定性能，我們分別從實(shí)測數(shù)據(jù)中截取強(qiáng)能量單頻噪聲與理論地震信號(hào)疊加進(jìn)行測試.單頻噪聲產(chǎn)生的原因有儀器噪聲、高壓線或人為活動(dòng)、以及機(jī)器運(yùn)行等，在時(shí)間域中極易干擾甚至淹沒地震信號(hào)(張晟瑞等，2018).如圖4所示，在時(shí)間域中地震信號(hào)的P波初至完全被強(qiáng)噪聲掩蓋，在STA/LTA方法中如圖4c漏拾P波到時(shí)；VAR-AIC方法如圖4e在P波到時(shí)處也無顯著優(yōu)勢.在噪聲淹沒的情況下，NTFT-STA/LTA方法如圖4d所示能較好的拾取P波與S波到時(shí)，但是在面波結(jié)束之后，與地震信號(hào)處于同一頻段的噪聲容易發(fā)生誤拾取情況，NTFT-STA/LTA方法的拾取精度受閾值影響；相對于NTFT-AIC方法而言，如圖4f所示P波與S波到時(shí)特征明顯，地震信號(hào)P波與S波的到時(shí)拾取不依賴于閾值的選取.

圖4 強(qiáng)能量單頻噪聲干擾(-6 dB)下理論地震信號(hào)到時(shí)自動(dòng)拾取結(jié)果(a)疊加強(qiáng)能量單頻噪聲的地震信號(hào)及其標(biāo)準(zhǔn)時(shí)頻譜(b)；(c)STA/LTA方法的拾取結(jié)果；(d)NTFT-STA/LTA方法的拾取結(jié)果；(e)VAR-AIC方法的拾取結(jié)果；(f)NTFT-AIC方法的拾取結(jié)果.Fig.4 Automatic picking of theoretical seismic phase under high energy single frequency noise interference (-6 dB)(a)Seismic signals superimposed with high energy single frequency noise and their Normal Time-Frequency spectrum (b);(c)Picking result of STA/LTA method;(d)Picking result of NTFT based STA/LTA method;(e)Picking result of VAR-AIC method;(f)Picking result of NTFT based AIC method.

在不同環(huán)境噪聲下我們對四種方法進(jìn)行了對比分析，STA/LTA方法、VAR-AIC方法與NTFT結(jié)合使波形在分界點(diǎn)出現(xiàn)明顯的突變，提高了拾取結(jié)果的精度.在STA/LTA的計(jì)算結(jié)果中，振幅、頻率突變的響應(yīng)結(jié)果中包含較多的毛刺干擾，導(dǎo)致拾取到時(shí)時(shí)較難把控閾值的選擇.閾值過低容易引發(fā)誤拾，閾值過高則又導(dǎo)致漏拾，說明STA/LTA較難拾取強(qiáng)噪聲背景環(huán)境下的震相到時(shí).基于NTFT的STA/LTA方法計(jì)算結(jié)果較為光滑，信號(hào)振幅響應(yīng)優(yōu)勢明顯，利于到時(shí)的自動(dòng)拾取.相比VAR-AIC方法，基于NTFT的AIC方法具有較好的穩(wěn)定性與魯棒性，在不同噪聲干擾下均保持較高的拾取準(zhǔn)確率.

2.2 實(shí)測數(shù)據(jù)拾取結(jié)果對比

為了進(jìn)一步檢驗(yàn)本文方法對于實(shí)際地震P波和S波到時(shí)拾取的可靠性，將基于NTFT的STA/LTA方法與基于NTFT的AIC方法應(yīng)用于實(shí)測地震數(shù)據(jù)，我們選用了105條信噪比較低的波形(圖5)，均為汶川余震記錄的垂直分量.數(shù)據(jù)來自“余震捕捉AI大賽”提供的汶川地震余震信號(hào)(Fang et al.,2017).應(yīng)用四種方法對105條余震信號(hào)進(jìn)行到時(shí)拾取，并將四種算法拾取的結(jié)果與人工拾取的結(jié)果進(jìn)行了比較.實(shí)測數(shù)據(jù)為多個(gè)臺(tái)站的不同事件，波形記錄中P波與S波到時(shí)受到了強(qiáng)噪聲的干擾，甚至到時(shí)信息被噪聲淹沒，波形特征如圖5所示.105條波形P波的平均信噪比為2.9 dB，S波為2.5 dB.為了方便比較分析，實(shí)驗(yàn)中STA/LTA方法與NTFT-STA/LTA方法計(jì)算短窗口寬度S均設(shè)為1 s，長窗口寬度L均設(shè)為10 s，地震P波和S波到時(shí)拾取閾值均設(shè)為3.

圖5 105條余震記錄的波形特征Fig.5 The waveform characteristics of 105 aftershock records

圖6 實(shí)測地震信號(hào)的到時(shí)拾取(a)滑窗范圍選定過程；(b)滑窗直接從標(biāo)準(zhǔn)時(shí)頻譜中選取信號(hào)；(c)不同滑窗拾取的到時(shí)；(d)滑窗拾取的最佳到時(shí)與VAR-AIC方法拾取結(jié)果對比；(e)STA/LTA方法與NTFT-STA/LTA方法拾取結(jié)果對比.圖中紅色實(shí)線表示信號(hào)的P波到時(shí)，紅色虛線表示信號(hào)的S波到時(shí)，綠色的虛線表示面波結(jié)束時(shí)刻.|A|表示對信號(hào)的振幅取模.Fig.6 Picking-up arrival time of the measured seismic signal(a)The selection process of sliding window areas;(b)The sliding window selects the signal directly from the Normal Time-Frequency spectrum;(c)Arrival times of seismic phases picked up by different sliding windows;(d)The optimal pickup results of sliding window and VAR-AIC method;(e)Results picked up by STA/LTA method and NTFT-STA/LTA method.In the figure,the solid red line represents the P-wave arrival time,the dotted red line represents the S-wave arrival time,and the dotted green line represents the end of the surface wave.|A|means the modulo operation of the amplitude of the signal.

基于本文提出的方法對105組實(shí)測地震數(shù)據(jù)進(jìn)行了P波、S波到時(shí)的拾取，實(shí)驗(yàn)中的余震記錄均按照人工拾取的P波到時(shí)為30 s統(tǒng)一截取.文中四種方法均對原始地震記錄直接進(jìn)行到時(shí)拾取，前30組余震記錄到時(shí)拾取結(jié)果對比如表2所示.根據(jù)表2可知，噪聲的干擾導(dǎo)致STA/LTA方法漏拾次數(shù)較多，且S波到時(shí)的漏拾次數(shù)遠(yuǎn)遠(yuǎn)大于P波，S波到時(shí)的漏拾率為9.52%，P波為0.95%.通過NTFT對STA/LTA特征函數(shù)的改進(jìn)，S波漏拾率降為1.90%，P波則沒有漏拾.VAR-AIC方法在低信噪比地震信號(hào)到時(shí)拾取中存在一定的漏拾，P波漏拾率為0.95%，S波漏拾率為1.90%；將VAR-AIC與NTFT結(jié)合對低信噪比地震信號(hào)拾取到時(shí)，P、S波到時(shí)漏拾情況得到較好的改進(jìn).可見STA/LTA與AIC方法中引入NTFT有效抑制了地震信號(hào)頻段外的強(qiáng)噪聲干擾.

表2 前30組地震記錄P波和S波自動(dòng)拾取結(jié)果(單位：s)Table 2 Automatic picking results of P and S waves in 30 groups of seismic records (unit：s)

為了檢驗(yàn)本文方法的有效性，我們將相同的30組余震記錄通過Vaezi和Van Der Baan(2015)的PSD(Power Spectral Density)方法進(jìn)行到時(shí)拾取，參照人工拾取結(jié)果統(tǒng)計(jì)誤差，將NTFT-STA/LTA方法、NTFT-AIC方法與PSD方法的P波、S波到時(shí)拾取誤差進(jìn)行了比較如圖7所示，其中在PSD方法中有11組數(shù)據(jù)發(fā)生P波或S波到時(shí)漏拾的情況.30組數(shù)據(jù)為不同臺(tái)站的多個(gè)事件，在背景噪聲復(fù)雜多樣的情況下，PSD方法對P波初至不敏感拾取誤差較大，相比之下，NTFT-STA/LTA方法與NTFT-AIC方法對P波、S波到時(shí)的拾取精度較高.

圖7 三種方法到時(shí)拾取誤差比較(a)P波到時(shí)拾取的誤差比較；(b)S波到時(shí)拾取的誤差比較.Fig.7 Error comparison of arrival time pickup of three methods(a)Error comparison of P-wave pickup;(b)Error comparison of S-wave pickup.

另外，我們應(yīng)用了STA/LTA方法、NTFT-STA/LTA方法、AIC方法與NTFT-AIC方法對105條余震記錄拾取震相到時(shí)，其中11條余震記錄包含P波或S波發(fā)生漏拾取的情況，其余的余震記錄拾取結(jié)果對比分析如圖8所示.將四種方法拾取的到時(shí)結(jié)果與人工拾取結(jié)果對比分析，NTFT的引入較好的改善了震相到時(shí)的拾取精度；背景噪聲復(fù)雜多變，STA/LTA方法與AIC方法拾取精度受噪聲影響較大，NTFT分別引入兩種方法后，拾取精度均得到了較好的改善，但相比于NTFT-STA/LTA方法而言NTFT-AIC方法拾取精度更為穩(wěn)定.

圖8 四種方法到時(shí)拾取結(jié)果(a)P波到時(shí)拾取的結(jié)果；(b)S波到時(shí)拾取的結(jié)果.Fig.8 Picking results of the arrival times with four methods(a)Picking results of P-wave arrival times;(b)Picking results of S-wave arrival times.

如表3所示，我們計(jì)算并統(tǒng)計(jì)了105組實(shí)測地震記錄到時(shí)拾取結(jié)果的誤差.表3中誤拾次數(shù)是指地震P波和S波到時(shí)拾取誤差大于1 s的事件次數(shù)；漏拾次數(shù)指未拾取到地震P波或S波到時(shí)的事件次數(shù).STA/LTA方法到時(shí)拾取精度依賴于閾值的選取，在強(qiáng)噪聲干擾的環(huán)境下，STA/LTA方法因強(qiáng)噪聲干擾導(dǎo)致較高的誤拾率或漏拾率.相比于STA/LTA方法，AIC方法拾取準(zhǔn)確率較高，但在強(qiáng)噪聲淹沒地震信號(hào)時(shí)，AIC計(jì)算窗口的劃分具有不確定性，NTFT的引入可直觀設(shè)置AIC計(jì)算窗口初始參數(shù).在給定相同參數(shù)下，NTFT的引入不需要通過閾值反復(fù)權(quán)衡誤拾率和漏拾率便可以獲得較為理想的地震P波或S波到時(shí).基于NTFT改進(jìn)的STA/LTA方法P波拾取結(jié)果的均方根誤差為0.36 s，S波拾取結(jié)果的均方根誤差為0.56 s；相對容易遺漏較難拾取的S波，NTFT引入后S波的漏拾率降低了1/5；NTFT-AIC方法優(yōu)于NTFT-STA/LTA方法，P波拾取結(jié)果的均方根誤差為0.25 s，S波拾取結(jié)果的均方根誤差為0.35 s.四種方法拾取P波與S波的誤差離散分布情況如圖9所示，相比S波的拾取，P波拾取誤差相對較小.NTFT的引入使STA/LTA方法與AIC方法到時(shí)拾取準(zhǔn)確率均有提高，并較好的改善了P波、S波漏拾的情況，這是本文方法具有比較大實(shí)用性的地方，充分發(fā)揮了NTFT的無偏性與STA/LTA、VAR-AIC方法的靈敏性.

圖9 四種方法拾取結(jié)果的誤差離散分布情況統(tǒng)計(jì)圖中“+”表示拾取誤差較大的異常值.Fig.9 The error discrete distribution statistics of seismic phase pickup results by four methodsThe “+”in the figure means the outlier with a large pickup error.

表3 對105組地震記錄P波和S波拾取結(jié)果的誤差統(tǒng)計(jì)Table 3 Error statistics of P and S wave pickup results of 105 groups of seismic records

3 討論

地震信號(hào)復(fù)雜多樣，實(shí)現(xiàn)地震信號(hào)P波、S波到時(shí)準(zhǔn)確、穩(wěn)定、快速拾取難度較大.本文研究顯示，將NTFT引入STA/LTA與AIC方法中，還存在兩個(gè)主要的問題：

(1)地震信號(hào)到時(shí)的拾取仍存在誤差.P波與S波到時(shí)較相近時(shí)，S波受P波尾波與噪聲的干擾導(dǎo)致到時(shí)拾取誤差偏大，S波相對于P波較難拾取.

(2)當(dāng)噪聲與地震信號(hào)分布于同一頻段、同一時(shí)間段時(shí)，極易混淆模糊震相，導(dǎo)致震相到時(shí)誤拾或漏拾的情況.

如圖10所示，圖中豎實(shí)線表示人工拾取的地震P波的到時(shí)，豎虛線表示人工拾取的地震S波的到時(shí).在P波與S波到時(shí)相隔較近且受強(qiáng)噪聲干擾的情況下，STA/LTA方法與AIC方法均較難準(zhǔn)確拾取S波.在STA/LTA方法中，引入NTFT后，如圖10c中所示，S波雖能被拾取，但是與人工拾取的結(jié)果存在較大的誤差.如圖10d中所示，VAR-AIC方法與NTFT-AIC方法對S波到時(shí)拾取失效.對這些信噪比低，P波震相與S波震相相隔較近的地震事件，自動(dòng)精確拾取到時(shí)還存在一定的困難.我們將NTFT與CNN結(jié)合(姚彥吉等，2020)，應(yīng)用NTFT+CNN算法識(shí)別地震信號(hào)，并提高震相到時(shí)拾取精度；還可以聯(lián)合多臺(tái)站技術(shù)(Zhang and Wen,2015)與本文所提出的方法來識(shí)別并提取震相到時(shí).

圖10 震相漏拾情況分析(a)實(shí)測地震記錄；(b)與(a)對應(yīng)的地震信號(hào)的標(biāo)準(zhǔn)時(shí)頻譜；(c)STA/LTA方法與NTFT-STA/LTA方法拾取結(jié)果與人工拾取結(jié)果對比；(d)VAR-AIC方法、NTFT-AIC方法拾取結(jié)果與人工拾取結(jié)果對比.圖中紅色和藍(lán)色的豎直實(shí)線表示人工拾取的P波到時(shí)，紅色和藍(lán)色的豎直虛線表示人工拾取的S波到時(shí).Fig.10 Analysis of seismic phase leakage pickup(a)The measured seismic record;(b)The Normal Time Frequency spectrum of the seismic signal corresponding to (a);(c)Comparing the results picked up by the STA/LTA method and NTFT-STA/LTA method with the manual results;(d)Comparing the results picked up by the VAR-AIC method and NTFT-AIC method with the manual results.In the figure,the red and blue vertical solid lines represent the P-wave arrival time picked by manual,while the red and blue vertical dotted lines represent the S-wave arrival time picked by manual.

4 結(jié)論

隨著地震臺(tái)站的不斷增設(shè)，地震數(shù)據(jù)急劇增多，強(qiáng)噪聲的干擾給海量地震數(shù)據(jù)的處理帶來了挑戰(zhàn)，針對地震P波和S波到時(shí)自動(dòng)、快速、準(zhǔn)確拾取的問題，本文在STA/LTA、AIC方法的基礎(chǔ)上，引入了NTFT，結(jié)合信號(hào)時(shí)間域與頻率域特征，提出了基于NTFT的STA/LTA方法，以及基于NTFT的AIC方法.

通過理論分析，合成數(shù)據(jù)及實(shí)測數(shù)據(jù)測試，基于NTFT的地震P波與S波自動(dòng)拾取方法繼承了STA/LTA、AIC方法振幅響應(yīng)敏感的優(yōu)點(diǎn)，同時(shí)還具有以下特征：

(1)基于NTFT的STA/LTA方法，建立了NTFT與特征函數(shù)的內(nèi)在聯(lián)系，該方法能較好的壓制地震信號(hào)頻段外噪聲對震相拾取的干擾.STA/LTA方法在強(qiáng)噪聲環(huán)境下，其拾取結(jié)果的誤拾率與漏拾率顯著增大，基于NTFT的STA/LTA方法計(jì)算確定地震信號(hào)所在的即時(shí)頻段，通過提取出地震信號(hào)所在的頻段有效拾取震相到時(shí).

(2)基于NTFT的AIC方法，能夠適應(yīng)多種噪聲的復(fù)雜環(huán)境，具有一定的穩(wěn)定性，并且地震P波、S波到時(shí)拾取誤差均較小，NTFT的引入進(jìn)一步增強(qiáng)了AIC方法的魯棒性.基于NTFT的AIC方法，結(jié)合了AIC方法的準(zhǔn)確性與NTFT方法的無偏特性，提高了震相自動(dòng)拾取的準(zhǔn)確率.

(3)NTFT對STA/LTA與AIC方法的改進(jìn)，均增強(qiáng)了振幅響應(yīng)變化特征，響應(yīng)處振幅較大且兩邊較光滑，毛刺較少，減少了對到時(shí)拾取的干擾.

致謝感謝中國地震局地球物理研究所趙明副研究員提供整理的波形震相數(shù)據(jù)集.