劉?？刀≈痉謇?媛崔仁勝1）中國(guó)天津300180中國(guó)地震局第一監(jiān)測(cè)中心2）中國(guó)北京100081中國(guó)地震局地球物理研究所3）中國(guó)北京100036中國(guó)地震局地震預(yù)測(cè)研究所

?

STFT與WVD在地震波時(shí)頻分析中的應(yīng)用

劉?？?), 2)丁志峰2)李 媛1)崔仁勝3)
1）中國(guó)天津300180中國(guó)地震局第一監(jiān)測(cè)中心2）中國(guó)北京100081中國(guó)地震局地球物理研究所3）中國(guó)北京100036中國(guó)地震局地震預(yù)測(cè)研究所

摘要地震波是一種非平穩(wěn)信號(hào)，傳統(tǒng)的傅里葉變換已不適用，時(shí)頻分析方法能在時(shí)間—頻率域給出非平穩(wěn)信號(hào)不同時(shí)刻的頻率分布特征。Wigner-Ville分布（WVD）對(duì)非平穩(wěn)信號(hào)具有較好的時(shí)頻聚集性，但因其為雙線(xiàn)性變換，會(huì)產(chǎn)生嚴(yán)重的交叉干擾項(xiàng)。短時(shí)傅里葉變換（STFT）是一種線(xiàn)性變換，不存在交叉項(xiàng)干擾，但時(shí)頻聚集性較差?？紤]到二者優(yōu)點(diǎn)，在對(duì)地震波進(jìn)行時(shí)頻分析時(shí)，采用短時(shí)傅里葉變換和Wigner-Ville分布的聯(lián)合算法，既有效抑制了Wigner-Ville分布交叉項(xiàng)干擾，又保持了較高的時(shí)頻聚集性。通過(guò)3種方法，對(duì)模擬信號(hào)和實(shí)際數(shù)據(jù)處理結(jié)果進(jìn)行對(duì)比分析，發(fā)現(xiàn)聯(lián)合算法對(duì)獲取地震波信號(hào)的時(shí)頻特征和縱波到時(shí)檢測(cè)，具有更好效果。

關(guān)鍵詞短時(shí)傅里葉變換；Wigner-Ville分布；時(shí)頻分析；地震波譜

E-mail：xikangliu@126.com，電話(huà)：15501065900

本文收到日期：2015-11-02

0 引言

根據(jù)信號(hào)的統(tǒng)計(jì)特征，將信號(hào)分為平穩(wěn)信號(hào)和非平穩(wěn)信，各統(tǒng)計(jì)量（如頻率、相關(guān)函數(shù)、功率譜等）隨時(shí)間變化的信號(hào)為非平穩(wěn)信號(hào)（張帆等，2006）。傳統(tǒng)的傅里葉變換，對(duì)信號(hào)的表征要么完全是時(shí)間域，要么完全是頻率域，不能分析信號(hào)中頻率隨時(shí)間的變化關(guān)系。因此，傅里葉變換無(wú)法表述信號(hào)的時(shí)間—頻率局域性質(zhì)，而此性質(zhì)恰是非平穩(wěn)信號(hào)的根本性質(zhì)。

時(shí)頻分析或稱(chēng)時(shí)頻分布，是描述信號(hào)頻率隨時(shí)間變化的方法，采用時(shí)間—頻率聯(lián)合表示信號(hào)，將一維時(shí)域信號(hào)映射到二維時(shí)頻平面，在時(shí)頻域內(nèi)對(duì)信號(hào)進(jìn)行分析，全面觀(guān)測(cè)信號(hào)的時(shí)間—頻率聯(lián)合特征，同時(shí)掌握信號(hào)的時(shí)域及頻域信息（趙淑紅，2010；劉?？档?，2013）。

短時(shí)傅里葉變換（short-time Fourier transform，STFT）是一種線(xiàn)性變換，對(duì)多分量的非平穩(wěn)信號(hào)不會(huì)產(chǎn)生交叉干擾項(xiàng)，但由于采用固定窗，使得時(shí)間和頻率分辨率難以同時(shí)達(dá)到最優(yōu)，時(shí)頻分辨率較差。Wigner-Ville分布（Wigner-Ville distribution，WVD）具有理論上的最高時(shí)頻分辨率和許多優(yōu)良的數(shù)學(xué)性質(zhì)（王見(jiàn)等，2013），但因其為雙線(xiàn)性變換，對(duì)多分量的非平穩(wěn)信號(hào)存在嚴(yán)重交叉項(xiàng)干擾。因此，可以通過(guò)某種算法，將STFT 與WVD關(guān)聯(lián)起來(lái)，即STFT與WVD的聯(lián)合算法，在不失去WVD良好時(shí)頻聚集性的情況下，消除交叉項(xiàng)干擾，得到同時(shí)具有較好時(shí)頻聚集性且無(wú)交叉項(xiàng)的地震波譜圖，提高地震波時(shí)頻分析的可讀性。目前，二者聯(lián)合表示方法主要是，利用短時(shí)傅里葉變換線(xiàn)性變換無(wú)交叉項(xiàng)的特性，抑制Wigner-Ville分布交叉項(xiàng)干擾。呂宙等（2009）主要應(yīng)用于航天電子對(duì)抗雷達(dá)線(xiàn)性調(diào)頻信號(hào)檢測(cè)研究；趙淑紅（2010）利用二者直接相乘方法，提取模擬線(xiàn)性調(diào)頻地震波信號(hào)的瞬時(shí)頻率；王見(jiàn)等（2013）利用二者聯(lián)合方法進(jìn)行機(jī)械轉(zhuǎn)輪故障檢測(cè)研究。本文利用二者聯(lián)合表示方法，用于地震波時(shí)頻分析和縱波到時(shí)檢測(cè)。

1 算法原理

1.1 WVD的定義

假設(shè)實(shí)際信號(hào)x(t)的解析信號(hào)為s(t)，其WVD定義為

WVD是一種雙線(xiàn)性變換，不滿(mǎn)足疊加原理。對(duì)于2個(gè)分量的信號(hào)：x(t) = x1(t)+x2(t)，其WVD結(jié)果為

其中，信號(hào)自項(xiàng)為Wx1、Wx2，而交叉干擾項(xiàng)為Wx1x2。WVD的交叉干擾項(xiàng)是無(wú)法避免的，嚴(yán)重影響區(qū)分信號(hào)項(xiàng)和交叉項(xiàng)，多分量信號(hào)交叉項(xiàng)干擾更為復(fù)雜、嚴(yán)重。

圖1(a)所示為包含4個(gè)高斯元信號(hào)的時(shí)域信號(hào)，采樣率為1 Hz，振幅均為1。在32 s處，兩個(gè)高斯元信號(hào)頻率中心分別為0.15 Hz和0.35 Hz，持續(xù)時(shí)間20 s；在96 s處，兩個(gè)高斯元信號(hào)頻率中心分別為0.15 Hz和0.35 Hz，持續(xù)時(shí)間也為20 s。圖1(b)為該信號(hào)的WVD分布，可以看出其具有明顯的交叉干擾項(xiàng)，無(wú)法分辨原信號(hào)分量和交叉干擾項(xiàng)。

因而，為了減小交叉項(xiàng)干擾，對(duì)設(shè)計(jì)不同核函數(shù)提出很多方法，例如：偽Wigner-Ville分布（PWVD）、平滑Wigner-Ville分布（SWVD）、Born-Jordan分布、Choi-Williams分布（CWD），等等（趙淑紅，2010；劉希康等，2013），統(tǒng)稱(chēng)為Cohen類(lèi)時(shí)頻分布。

1.2 STFT及其譜

信號(hào)x(t)的STFT定義為

由于信號(hào)x(t)乘以一個(gè)短窗函數(shù)w(t)等價(jià)于在分析時(shí)間點(diǎn)t附近的取一個(gè)切片，所以短時(shí)傅里葉變換可以理解為信號(hào)x(t)在時(shí)間點(diǎn)t附近的傅里葉變換，即“局部頻譜”（胡昌華等，2002；呂宙等2009）。

由定義可知，短時(shí)傅里葉變換是一種線(xiàn)性變換，對(duì)多分量信號(hào)來(lái)說(shuō)，不存在交叉項(xiàng)。由不確定原理可知，短時(shí)傅里葉變換（短時(shí)傅里葉變換英文縮寫(xiě)STFT，為了描述連貫采用中文全稱(chēng)）的時(shí)間分辨率和頻率分辨率不可能同時(shí)變小，其時(shí)寬—帶寬乘積存在一個(gè)下限，短時(shí)傅里葉變換的時(shí)頻聚集性由此受到限制。由于高斯窗函數(shù)達(dá)到時(shí)寬—帶寬積的下界（科恩，2000），因此本文進(jìn)行短時(shí)傅里葉變換使用高斯窗函數(shù)。

WVD是一種二次型變換，為了與WVD對(duì)比，一般使用短時(shí)傅里葉變換譜。短時(shí)傅里葉變換譜的定義為STFT模的平方

圖1(c)為上述高斯元信號(hào)的短時(shí)傅里葉變換譜。與圖1(b)相比，短時(shí)傅里葉變換譜無(wú)交叉干擾項(xiàng)，但時(shí)間頻率分布范圍較寬，時(shí)頻聚集性比WVD差得多。

圖1　四分量高斯元信號(hào)、能量譜密度及其時(shí)頻分布（a）高斯元信號(hào)能量譜密度及其理論時(shí)頻；（b）信號(hào)（a）的WVD分布；（c）信號(hào)（a）的STFT譜分布；（d）信號(hào)（a）的STFT譜和Wigner-Ville聯(lián)合方法分布Fig.1 Signal of 4 Gauss atoms, energy density and the time-frequency diagrams

1.3 STFT與WVD的聯(lián)合算法

從STFT和WVD各自?xún)?yōu)缺點(diǎn)出發(fā)，可以得出，由于STFT是線(xiàn)性變換，不存在交叉干擾項(xiàng)，能夠正確反映信號(hào)的各分量信息，但時(shí)頻聚集性不好，分辨率較低；對(duì)于WVD，具有很好的時(shí)頻聚集性及較高時(shí)頻分辨率，但存在嚴(yán)重交叉干擾項(xiàng)。因此，可以對(duì)二者進(jìn)行聯(lián)合計(jì)算，結(jié)合二者優(yōu)點(diǎn)，可得到既不含交叉干擾項(xiàng)而時(shí)頻分辨率又高的分布。

當(dāng)?shù)玫叫盘?hào)的STFT譜和WVD后，二者聯(lián)合算法可定義（王見(jiàn)等，2013）為

式（6）為取STFT譜和WVD中數(shù)值較小者，按此處理，將WVD中有交叉項(xiàng)部分用相應(yīng)STFT譜中的數(shù)值代替，以達(dá)到消除交叉項(xiàng)的目的。式（7）中c為STFT譜的交叉項(xiàng)消除閾值。當(dāng)STFT譜的值小于該閾值時(shí)，返回0；如果大于該閾值，則返回1。WVD中有交叉項(xiàng)對(duì)應(yīng)STFT譜部分的數(shù)值肯定小于該閾值，用數(shù)字0與WVD相乘以消除交叉項(xiàng)。式（8）中“·”代表矩陣的點(diǎn)乘，α、β稱(chēng)為冪調(diào)節(jié)系數(shù)。該式的作用是增強(qiáng)兩變換數(shù)值較大部分而消弱有交叉項(xiàng)部分。一般隨著α、β的增大，STFT和WVD聯(lián)合變換的時(shí)頻聚集性會(huì)提高，在實(shí)際應(yīng)用中，需靈活調(diào)試α、β取值，以免丟失某些有用信息。式（8）較式（6）、式（7）結(jié)果更加理想，使用更為靈活。

圖1(d)為STFT譜和WVD聯(lián)合算法對(duì)圖1(a)信號(hào)的處理結(jié)果，可以看出，該方法既克服短時(shí)傅里葉變換固定時(shí)頻分辨率的缺陷，又在消除WVD交叉項(xiàng)干擾的同時(shí)保留了良好的時(shí)頻聚集性，與信號(hào)的理論頻率保持一致，具有較高的時(shí)間—頻率分辨率。

2 資料處理結(jié)果

2.1 地震波譜時(shí)頻分析

以2013年9月4日莆田ML4.8地震為例，選取福州日溪地震臺(tái)（FZRX）南北向分量記錄進(jìn)行地震波譜時(shí)頻分析。臺(tái)站距震中約102 km，采用CMG-3ESP型寬頻帶地震儀，采樣頻率100 Hz。

由于聯(lián)合算法對(duì)計(jì)算機(jī)內(nèi)存需求較高，為降低計(jì)算量，截取2013年9月4日06時(shí)23 分39秒至2013年9月4日06時(shí)24分29秒，時(shí)長(zhǎng)為50 s的地震波形進(jìn)行時(shí)頻分析（假設(shè)截取波形開(kāi)始時(shí)刻為0秒）[圖2(a)]，能量譜密度見(jiàn)圖2(b)，可以發(fā)現(xiàn)，其能量主要集中在小于5 Hz的頻率范圍內(nèi)，少部分在5—10 Hz，說(shuō)明該記錄為低于5 Hz頻率成分為主的地震波。圖2（c）為莆田ML4.8地震波形的Wigner-Ville分布，能夠較好反映該地震波的時(shí)頻分布特征，頻率和能量譜密度分布一致，15 s前后橫波到達(dá)，能量達(dá)到最強(qiáng)，由于縱波能量相對(duì)較弱，時(shí)頻域內(nèi)縱波未達(dá)到最大幅值的5%以上，圖中沒(méi)有很好地顯示出來(lái)。另外，由于地震波的非平穩(wěn)特性，受噪聲和自項(xiàng)相互干擾等影響，不難發(fā)現(xiàn)，時(shí)頻圖像存在較多干擾成分，與有用信息混合在一起。

圖2　福州日溪地震臺(tái)（FZRX）南北向記錄及其能量譜密度和Wigner-Ville分布（a）南北向記錄波形；（b）能量譜密度；（c）Wigner-Ville分布Fig.2 Energy density and Wigner-Ville distribution of FZRX Seismic Station NS component waveform record

2.2 縱波到時(shí)檢測(cè)

為清楚驗(yàn)證STFT和WVD聯(lián)合方法對(duì)提取地震波波譜時(shí)頻特征和縱波到時(shí)檢測(cè)的有效性，截取上述地震臺(tái)站記錄2013年9月4日06時(shí)23分39秒至2013年9月4日06時(shí)23分54秒，包含P波前后共15 s的三分量波形進(jìn)行分析（假設(shè)波形開(kāi)始時(shí)刻為0秒），圖3（a）從左向右依次為N、E、U三方向，虛豎線(xiàn)表示P波初動(dòng)時(shí)刻，在該時(shí)段內(nèi)，初動(dòng)時(shí)間讀取為3.85 s。

圖3　P波三分量記錄波形及其不同時(shí)頻分布（a）三分量P波記錄波形；（b）STFT譜和Wigner-Ville聯(lián)合方法時(shí)頻分布；（c）Wigner-Ville分布；（d）STFT譜Fig.3 Three component of P waveform and various time-frequency distribution diagrams

圖3（b）為利用二者聯(lián)合算法，由式（8）所得結(jié)果，冪調(diào)節(jié)系數(shù)取α = 1.5，β = 1，可以發(fā)現(xiàn)，聯(lián)合算法的時(shí)頻聚集性較好，同時(shí)消除了交叉項(xiàng)干擾，地震記錄頻率范圍主要集中在1—5 Hz，少量分布在5—10 Hz，在P波初動(dòng)時(shí)刻前后能量較強(qiáng)，且P波頻率成分和分布形態(tài)在三方向上具有較好的一致性，具有從低頻向高頻變化的明顯趨勢(shì)。

圖3（c）為三方向的STFT譜，選取高斯窗函數(shù)長(zhǎng)度為129，結(jié)果顯示，因沒(méi)有交叉項(xiàng)干擾，STFT譜具有較清晰的時(shí)頻分布圖，但是其時(shí)頻聚集性明顯較差，各時(shí)刻頻率分布散亂，分辨率較差，不能有效判讀P波初動(dòng)時(shí)刻，頻率成分主要在0—10 Hz。

圖3（d）為三方向的Wigner-Ville分布，可以看出，其能量主要集中在5 Hz以下，并且有從低頻向高頻變化的趨勢(shì)，能夠判讀P波初動(dòng)時(shí)間，且與該時(shí)段內(nèi)波形記錄文件讀取到的P波到時(shí)一致，為3.85 s，這也為精確檢測(cè)P波初動(dòng)時(shí)刻提供了一種方法，但整體看由于噪聲干擾和交叉項(xiàng)干擾，有效信號(hào)被淹沒(méi)，存在虛假高頻信息，時(shí)頻分辨率并不理想。

通過(guò)以上結(jié)果對(duì)比可以發(fā)現(xiàn)，STFT和WVD聯(lián)合算法既保留了STFT的線(xiàn)性變換沒(méi)有交叉干擾項(xiàng)的特性，又保留了WVD的時(shí)頻高分辨特性，具有較好的時(shí)頻分布特性，且可以通過(guò)靈活調(diào)整冪調(diào)節(jié)系數(shù)，根據(jù)實(shí)際需要得到待分析波形的不同時(shí)頻分布形態(tài)。根據(jù)波形讀取的P波初動(dòng)和利用聯(lián)合方法得到的P波初動(dòng)時(shí)刻一致，說(shuō)明該方法對(duì)P波初動(dòng)精確檢測(cè)具有良好效果。

3 結(jié)論

本文通過(guò)對(duì)模擬信號(hào)和實(shí)際地震記錄處理分析，對(duì)比STFT譜、Wigner-Ville分布及聯(lián)合方法的處理結(jié)果，得出以下結(jié)論。

（1）實(shí)際應(yīng)用發(fā)現(xiàn)，短時(shí)傅里葉變換受窗口大小影響，時(shí)窗太大會(huì)導(dǎo)致時(shí)間分辨率降低，時(shí)窗太小，頻率分辨率較低（趙淑紅，2010），時(shí)間分辨率和頻率分辨率不能同時(shí)達(dá)到最優(yōu)，窗口大小及類(lèi)型的選取需結(jié)合實(shí)際資料，考慮處理目的和要求。

（2）Wigner-Ville分布不包含任何窗函數(shù)，可以達(dá)到最高的時(shí)頻分辨率，但存在嚴(yán)重的交叉干擾項(xiàng)，交叉項(xiàng)干擾程度與資料自身所包含的頻率成分有關(guān)，對(duì)實(shí)際資料來(lái)說(shuō)，由于頻率成分復(fù)雜，自項(xiàng)、噪聲、干擾項(xiàng)相互疊加，時(shí)頻分辨率更差。

（3）利用STFT譜和Wigner-Ville分布聯(lián)合表示方法，對(duì)模擬信號(hào)和實(shí)際數(shù)據(jù)進(jìn)行處理可知，此方法結(jié)合STFT和WVD的優(yōu)點(diǎn)，既抑制了WVD交叉項(xiàng)干擾，又保留了時(shí)頻高分辨率特性，能夠精確表示地震波的波譜特征，并檢測(cè)P波初動(dòng)時(shí)刻。另外，可以根據(jù)實(shí)際數(shù)據(jù)處理需要，通過(guò)選擇合適的窗函數(shù)、時(shí)窗長(zhǎng)度及冪調(diào)節(jié)系數(shù)，在一定程度上靈活顯示待分析信號(hào)的時(shí)頻分布特征。

（4）較高、較好分辨率的時(shí)頻分析方法在爆破識(shí)別（沈萍等，1999）、縱波到時(shí)檢測(cè)（鄒文等，2006）、工程抗震參數(shù)設(shè)計(jì)（魏來(lái)等，2007）等方面均具有重要的研究意義和應(yīng)用前景。

參考文獻(xiàn)

胡昌華，周濤，夏啟兵，等. 基于Matlab的系統(tǒng)分析與設(shè)計(jì)——時(shí)頻分析[M].西安電子科技大學(xué)出版社，2002：2-24.

劉?？? 精密可控震源信號(hào)提取新方法研究[D].北京：中國(guó)地震局地震預(yù)測(cè)研究所，2013：4-6.

呂宙. 基于STFT的Wigner-Ville分布交叉項(xiàng)抑制[J]. 航天電子對(duì)抗，2009，26（3）：27-29.

沈萍，鄭治真. 瞬態(tài)譜在地震與核爆識(shí)別中的應(yīng)用[J]. 地球物理學(xué)報(bào)，1999，42（2）：233-240.

魏來(lái)，劉馬群，李奎. 信號(hào)分析方法及其在地震波處理中的應(yīng)用[J]. 洛陽(yáng)大學(xué)學(xué)報(bào)，2007，22（4）：69-73.

王見(jiàn)，李金同，盧華玲，等. 采用STFT-Wigner變換抑制Wigner-Ville分布交叉項(xiàng)[J]. 重慶大學(xué)學(xué)報(bào)，2013，36（8）：15-25.

張帆，鐘羽云，朱新云，等. 時(shí)頻分析方法及地震波譜研究中的應(yīng)用[J]. 地震地磁觀(guān)測(cè)與研究，2006，27（4）：17-22.

張?chǎng)?，趙擁軍. 基于短時(shí)傅里葉變換和Wigner-Ville分布的聯(lián)合變換[J]. 電子對(duì)抗，2008，3：39-42.

趙淑紅. 短時(shí)傅里葉變換與Wigner-Ville分布聯(lián)合確定地震信號(hào)瞬時(shí)頻率[J]. 西安科技大學(xué)學(xué)報(bào)，2010，30（4）：447-450.

鄒文，陳愛(ài)萍，顧漢明，等. S變換初至拾取技術(shù)及其在復(fù)雜山區(qū)中的應(yīng)用[J]. 石油天然氣學(xué)報(bào)，2006，28（2）：63-65.

[美]科恩著. 白居憲譯. 時(shí)—頻分析：理論與應(yīng)用[M]. 西安交通大學(xué)出版社，2000：35-40.

Application of short -time Fourier transform and Wigner -Ville distribution in time -frequency analysis of seismic wave

Liu Xikang1), 2），Ding Zhifeng2），Li Yuan1）and Cui Rensheng3）
1) First Crust Monitoring and Application Center， China Earthquake Administration， Tianjin 300180，China 2) Institute of Geophysics， China Earthquake Administration， Beijing 100081， China 3) Institute of Earthquake Science， China Earthquake Administration， Beijing 100036， China

Abstract

Seismic signal is a kind of non-stable signal， which can’t be analyzed suitably by classical Fourier analysis. Time-frequency analysis can analyze non-stationary signal in time-frequency domain and show out the distribution characteristics. WVD has perfect time-frequency concentration， but because it is a bilinear transformation， the existence of cross-terms is inevitable. STFT is linear transformation， has no cross-terms， but affected by the window size，it has bad time-frequency concentration. Based on the advantage of these two methods， we propose to obtain the time-frequency distribution by the combined method. This method can suppress the cross-term interference in WVD effectively， while keeping the excellent timefrequency concentration. By calculating simulation signal and actual data， the comparison result of STFT， WVD and the combined method indicates that the combination method is most suitable for computing the time-frequency distribution characteristics of earthquake waveform signal and more effective to detect the arrive time of P-waves.

Key words：short-time Fourier transformSTFT），Wigner-Ville distributionWVD），timefrequency analysis，earthquake spectrum

doi:10. 3969/j. issn. 1003-3246. 2016. 01. 003

基金項(xiàng)目：中國(guó)地震局“三結(jié)合”（No.163303），科技部公益性行業(yè)科研專(zhuān)項(xiàng)“中國(guó)地震科學(xué)臺(tái)陣探測(cè)——南北地震帶北段”（No.201308011），中國(guó)地震局第一監(jiān)測(cè)中心主任基金項(xiàng)目（FMC2015005）

作者簡(jiǎn)介：劉?？担?988—），男，助理工程師，在讀博士研究生，主要從事地震觀(guān)測(cè)技術(shù)等研究工作。