基于標準時頻變換的準勒夫波信號識別

程 威 柳林濤

1 中國科學院精密測量科學與技術創新研究院，武漢市徐東大街340號，430077 2 中國科學院大學，北京市玉泉路19號甲，100049

剪切波分裂測量是研究地震各向異性最成熟的方法之一[1]，但由于剪切波分裂的結果是其傳播路徑上不同深度范圍內所有各向異性信息的疊加，因此無法獲知各向異性的深度分布信息[2]。在利用散射面波研究傳播介質的各向異性時，為提高識別度，使其他高頻干擾信號(如交通信號等)被最大程度地壓制，通常使用低頻(>70 s)準勒夫(Quasi-Love)波信號[3-4]。目前Quasi-Love波已被應用于日本千島群島、阿留申島中部、湯加群島等環太平洋俯沖帶地區的上地幔各向異性梯度特性研究[5-7]，同時根據中國西藏地區臺站記錄的Quasi-Love波，證實了中國西藏地區存在地幔和地殼垂直耦合的邊界狀態[8]。然而，同剪切波方法相比，Quasi-Love波的應用還不夠廣泛。

由于小波變換更適合處理如地震信號類的非平穩信號，本文提出使用基于標準時頻變換[9-10](normal time-frequency transform，NTFT)的方法探測Quasi-Love波。NTFT是一種新的小波變換，與傳統的小波變換相比具有“無為原理”[11]屬性和強抗噪聲能力，目前已被成功應用于海潮信號分析和預測[12]、長周期地球自轉信號預測[13]及GPS鐘差信號分析[14]等地學領域。本文基于NTFT譜系數定義了評價2個信號相關程度的相似度概念，并使用該方法估計Quasi-Love波和Love波的相關性及二者的延時，首先在NTFT的基礎上引出相似度的概念，通過仿真實驗評價其抗噪聲能力，并將其與廣泛使用的互相關法進行比較，再將NTFT方法應用到Quasi-Love波的識別中。

1 方 法

本文首先介紹NTFT的基本理論，相似度的概念是基于NTFT時頻譜定義的，所以在NTFT的基礎上引出相似度的定義，并通過仿真信號測試相似度法在高背景噪聲下的性能。

1.1 標準時頻變換

對于一個復時間函數f(t)∈C，其標準時頻變換為：

(1)

(2)

(3)

(4)

1.2 相似度

估計信號相關性最常用的方法是互相關(cross correlation，CC)法，本文基于NTFT時頻譜定義類似的算法，使用信號的頻率和時域特征求2個信號的相關性。由于NTFT譜系數的模等于信號在時域中的實際振幅，所以這種基于NTFT時頻譜的算法直接反映了信號的“相似程度”，即為相似度(similarity coefficient，SC)法。對于2個時間函數f1(t)和f2(t)，其相似度函數ρ(s)定義為：

(5)

式中，ψf1和ψf2為函數f1(t)和f2(t)的標準時頻變換，Reψ為NTFT譜系數ψ(t,?)的實部，s為平滑時間因子，S為感興趣的區域。在實際應用中，首先將2個信號進行NTFT處理，得到NTFT譜系數。根據式(5)在二者的NTFT時頻譜中分別選定感興趣的區域計算相似度函數ρ(s)的值，在時間軸上移動選定的區域，當ρ(s)取得極值時，該最大值即為2個時間函數的相似度，對應的s即為時間延遲。

1.3 仿真測試

為了測試相似度法的抗噪聲能力，仿真2個信號的時間序列，在無噪聲和加入噪聲的情況下分別使用本文提出的SC法和經典CC法計算兩者的相關性，根據計算結果評價SC法的抗噪性能。測試使用的仿真信號表達式為：

(6)

式中，ε1和ε2為高斯白噪聲，n為時間延遲。設置初始頻率ω=1，采樣周期為0.1 s，時間延遲n=800 s，并加入噪聲水平為-2.5 dB的高斯白噪聲。2個仿真信號具有相同的表達式，在時域中有相同的波形，仿真信號的初始頻率為1 Hz，隨著時間的增加頻率逐漸增大。圖1(a)為仿真信號的FFT振幅譜，可以看出，仿真信號的主頻為1～2 Hz；圖1(b)為仿真信號的NTFT時頻譜，可以看出，時頻譜隨時間逐漸增大，2個信號的延時為800 s。在原始信號中加入較高水平的高斯白噪聲，圖2(a)為加入-2.5 dB高斯噪聲后仿真信號的FFT振幅譜，可以看到，加入噪聲后影響較大，信號的主頻已經不明顯；由圖2(b)可以發現，加入噪聲后NTFT時頻譜也受到一定程度的影響，但基本還能分辨出有2個周期隨時間衰減的信號。分別利用SC法和CC法對2個信號的相關性進行估計，圖3(a)和3(b)分別為無噪聲和較高噪聲情況下的結果，由圖可知，無噪聲時2種方法計算的結果幾乎重合，在延時為800 s處均取得極值，最大值約為1，與理論值一致，說明2種方法的性能一致；而在較高噪聲水平下，SC法估計2個仿真信號的相關系數約為0.8，最大值對應的時間延遲仍為800 s，而CC法已無法判斷2個信號是否為同一個信號。實驗結果表明，SC法具有較強的抗噪聲能力，在極低的信噪比情況下依然能得出正確的結論，相較于CC法，使用SC法有利于提高對相似信號識別的準確度。

圖1 無噪聲信號Fig.1 Signal without noise

圖2 加入噪聲后信號Fig.2 Signal with gaussian noise

圖3 SC法和CC法估計信號的相關性Fig.3 Correlation of signals estimated by SC method and CC method

2 數據和結果

本文利用2004-12-26蘇門答臘MW9.0(Global CMT)地震激發的面波數據研究Quasi-Love波，選用意大利半島的ELBR臺和CSTR臺記錄的地震數據，數據來源于IRIS數據管理中心，采樣率為50 Hz，2個臺站的震中距分別為84°和83°。為比較Quasi-Love波與Love波能量的大小，本文將水平分量數據旋轉至徑向(R)和切向(T)，同時為減小計算量，將原始波形數據降采樣至1 Hz。 為在NTFT時頻譜中觀測Quasi-Love波，分別將2個臺站記錄的Z向和T向數據作NTFT處理，周期范圍為50～200 s，核函數設計成標準Morlet小波[15]，并選擇Love波到時附近的數據段進行觀察，結果見圖4。可以看到，2個臺站T向的NTFT時頻譜十分相似，均出現了相同信號，該信號即為Love波；而Z向的NTFT時頻譜差異較大，ELBR臺在1 300～1 400 s(Love波附近)出現了一個較強的信號，CSTR臺則缺少該信號，推斷該信號為Quasi-Love波。為驗證該推斷，利用常規高階FIR濾波器對2個臺站的面波數據進行帶通濾波處理，濾波器階數為110，濾波范圍為10～20 mHz(50～100 s)，并將結果與NTFT時頻譜進行對比，圖5為2個臺站Z向和T向的地震波數據經低頻濾波后的結果。從圖5可以清晰地觀測到Love波和Rayleigh波，并且在ELBR臺Z向數據中發現有微弱的異常波形出現在Rayleigh波之前，根據對Quasi-Love波的觀測經驗可知，該異常波形即為Quasi-Love波，CSTR臺則無類似異常波形信號。

圖4 ELBR臺和CSTR臺Z向和T向的NTFT時頻譜Fig.4 NTFT time-frequency spectra of the Z componentand T component from ELBR station and CSTR station

圖5 ELBR和CSTR臺Z向和T向數據濾波后的波形Fig.5 The filtered waveforms of the Z and T components of ELBR station and CSTR station

NTFT時頻譜中觀測到的Quasi-Love波和由傳統方法得到的結果一致，但NTFT時頻譜中的更加明顯，有利于信號的識別。另外，NTFT時頻譜顯示出Quasi-Love波的時頻特征，利用傳統帶通濾波器對信號進行濾波后容易引起信號的相位偏移，而NTFT時頻譜中的Quasi-Love波不會出現相位偏移，這對估計Quasi-Love波的延時十分重要。使用SC法估算Quasi-Love波的延時，并利用式(5)估計ELBR臺T向數據中Love波和Z向數據中Quasi-Love波的相關性，結果見圖6。由圖可知，2個信號的相關性峰值為0.6，對應的時間延遲為2 s，說明Quasi-Love波和Love波幾乎同時到達臺站，且Quasi-Love波的散射體應在臺站附近，說明臺站附近存在各向異性梯度介質；而長周期面波通常對上地幔敏感，說明ELBR臺附近的上地幔存在各向異性梯度。

圖6 SC法估計ELBR臺Quasi-Love波和Love波的相關性Fig.6 SC method estimates the correlation between Quasi-Love wave and Love wave of ELBR station

3 結 語

本文利用NTFT方法識別Quasi-Love波信號，基于NTFT時頻譜系數定義了估計2個信號相關性的SC法，并使用SC法估計了Quasi-Love波的延時。與傳統方法相比，NTFT方法在未對信號作任何濾波處理的情況下成功識別出了Quasi-Love波信號，并能在NTFT時頻譜中觀察到信號的頻率隨時間的變化。另外，本文通過NTFT方法研究了意大利半島2個臺站記錄的2004年蘇門答臘大地震的面波數據，在ELBR臺識別出Quasi-Love波，并根據SC法計算了Quasi-Love波的延時，但未對其產生機理作詳細研究。Quasi-Love波的出現同上地幔各向異性梯度存在緊密聯系，而地震各向異性是地球動力學研究的重要內容。目前各向異性探測工具以剪切波分裂技術為主，但隨著高質量地震波數據的輕松獲取，Quasi-Love波在地球動力學研究方面的應用具有較好的前景。

致謝：本文使用的地震波數據來源于地震學研究聯合會(IRIS)數據管理中心(DMC)，部分圖件由GMT軟件繪制，在此表示感謝。