利用標準時頻變換方法在強噪聲環境下無偏拾取地震P波、S波到時

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

1 中國科學院精密測量科學與技術創新研究院大地測量與地球動力學國家重點實驗室，武漢 430077 2 中國科學院大學地球與行星科學學院，北京 100049

0 引言

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

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

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

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

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

1 方法原理

1.1 標準時頻變換

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

(1)

(2)

式中，w(t)表示窗函數.在標準時頻變換中，核函數的傅里葉變換式(3)滿足條件(4)，

(3)

(4)

(5)

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

(6)

構建的標準小波變換表達為：

(7)

對一個調和信號(8)式h(t)進行標準時頻變換時，可滿足(9)式的即時無偏特性.

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

(8)

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

(9)

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

(10)

(11)

所以Ψh(τ,β)=Aexp(j(φ+βτ))=h(τ)，一個調和信號的標準時頻譜無偏的描述原始信號的特征.

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

圖1 地震信號及其時頻譜分析對比(a)地震信號；(b)莫雷特連續小波變換；(c)S變換；(d)標準時頻變換.圖(a)中紅色點為地震信號振幅最大的時刻，時頻譜中黃色虛線為地震信號紅點變換后對應的時刻.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變換可確定地震信號的即時信息，結合NTFT時頻譜特征對STA/LTA方法的特征函數進行改進，通過NTFT確定的即時頻率對地震信號發生頻段進行約束，選取地震信號發生的即時頻段[ωp,ωq]進行震相拾取，即時頻段的判別過程參見(姚彥吉等，2020)，頻段約束定義如下：

(12)

(13)

(14)

因此，文中改進特征函數表示為

(15)

(16)

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

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

1.3 基于NTFT的AIC方法

(17)

(18)

或者表示為：

(19)

同理，根據標準時頻譜判斷信號發生的時間段，設置滑窗的時間范圍.

AIC函數定義為(劉勁松等，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)表示被時間段[m,n]與NTFT即時頻率段[ωp,ωq]約束的地震信號；var(G[m,k])和var(G[k+1,n])分別表示滑窗內數據集合G[m,k]與G[k+1,n]的方差.

2 算例分析

2.1 合成波形測試

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

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

(21)

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

(22)

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

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

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

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

表1 不同背景噪聲下地震震相自動拾取統計表Table 1 Automatic picking results of seismic phase under the different background noise

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

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

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

圖4 強能量單頻噪聲干擾(-6 dB)下理論地震信號到時自動拾取結果(a)疊加強能量單頻噪聲的地震信號及其標準時頻譜(b)；(c)STA/LTA方法的拾取結果；(d)NTFT-STA/LTA方法的拾取結果；(e)VAR-AIC方法的拾取結果；(f)NTFT-AIC方法的拾取結果.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.

在不同環境噪聲下我們對四種方法進行了對比分析，STA/LTA方法、VAR-AIC方法與NTFT結合使波形在分界點出現明顯的突變，提高了拾取結果的精度.在STA/LTA的計算結果中，振幅、頻率突變的響應結果中包含較多的毛刺干擾，導致拾取到時時較難把控閾值的選擇.閾值過低容易引發誤拾，閾值過高則又導致漏拾，說明STA/LTA較難拾取強噪聲背景環境下的震相到時.基于NTFT的STA/LTA方法計算結果較為光滑，信號振幅響應優勢明顯，利于到時的自動拾取.相比VAR-AIC方法，基于NTFT的AIC方法具有較好的穩定性與魯棒性，在不同噪聲干擾下均保持較高的拾取準確率.

2.2 實測數據拾取結果對比

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

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

圖6 實測地震信號的到時拾取(a)滑窗范圍選定過程；(b)滑窗直接從標準時頻譜中選取信號；(c)不同滑窗拾取的到時；(d)滑窗拾取的最佳到時與VAR-AIC方法拾取結果對比；(e)STA/LTA方法與NTFT-STA/LTA方法拾取結果對比.圖中紅色實線表示信號的P波到時，紅色虛線表示信號的S波到時，綠色的虛線表示面波結束時刻.|A|表示對信號的振幅取模.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組實測地震數據進行了P波、S波到時的拾取，實驗中的余震記錄均按照人工拾取的P波到時為30 s統一截取.文中四種方法均對原始地震記錄直接進行到時拾取，前30組余震記錄到時拾取結果對比如表2所示.根據表2可知，噪聲的干擾導致STA/LTA方法漏拾次數較多，且S波到時的漏拾次數遠遠大于P波，S波到時的漏拾率為9.52%，P波為0.95%.通過NTFT對STA/LTA特征函數的改進，S波漏拾率降為1.90%，P波則沒有漏拾.VAR-AIC方法在低信噪比地震信號到時拾取中存在一定的漏拾，P波漏拾率為0.95%，S波漏拾率為1.90%；將VAR-AIC與NTFT結合對低信噪比地震信號拾取到時，P、S波到時漏拾情況得到較好的改進.可見STA/LTA與AIC方法中引入NTFT有效抑制了地震信號頻段外的強噪聲干擾.

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

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

圖7 三種方法到時拾取誤差比較(a)P波到時拾取的誤差比較；(b)S波到時拾取的誤差比較.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.

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

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

圖9 四種方法拾取結果的誤差離散分布情況統計圖中“+”表示拾取誤差較大的異常值.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波拾取結果的誤差統計Table 3 Error statistics of P and S wave pickup results of 105 groups of seismic records

3 討論

地震信號復雜多樣，實現地震信號P波、S波到時準確、穩定、快速拾取難度較大.本文研究顯示，將NTFT引入STA/LTA與AIC方法中，還存在兩個主要的問題：

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

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

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

圖10 震相漏拾情況分析(a)實測地震記錄；(b)與(a)對應的地震信號的標準時頻譜；(c)STA/LTA方法與NTFT-STA/LTA方法拾取結果與人工拾取結果對比；(d)VAR-AIC方法、NTFT-AIC方法拾取結果與人工拾取結果對比.圖中紅色和藍色的豎直實線表示人工拾取的P波到時，紅色和藍色的豎直虛線表示人工拾取的S波到時.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 結論

隨著地震臺站的不斷增設，地震數據急劇增多，強噪聲的干擾給海量地震數據的處理帶來了挑戰，針對地震P波和S波到時自動、快速、準確拾取的問題，本文在STA/LTA、AIC方法的基礎上，引入了NTFT，結合信號時間域與頻率域特征，提出了基于NTFT的STA/LTA方法，以及基于NTFT的AIC方法.

通過理論分析，合成數據及實測數據測試，基于NTFT的地震P波與S波自動拾取方法繼承了STA/LTA、AIC方法振幅響應敏感的優點，同時還具有以下特征：

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

(2)基于NTFT的AIC方法，能夠適應多種噪聲的復雜環境，具有一定的穩定性，并且地震P波、S波到時拾取誤差均較小，NTFT的引入進一步增強了AIC方法的魯棒性.基于NTFT的AIC方法，結合了AIC方法的準確性與NTFT方法的無偏特性，提高了震相自動拾取的準確率.

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

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