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

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í)頻聚集性，但因其為雙線性變換，會(huì)產(chǎn)生嚴(yán)重的交叉干擾項(xiàng)。短時(shí)傅里葉變換（STFT）是一種線性變換，不存在交叉項(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è)，具有更好效果。【期刊名稱】地震地磁觀測(cè)與研究【年（卷），期】2016（037）001【總頁(yè)數(shù)】7【關(guān)鍵詞】關(guān)鍵詞短時(shí)傅里葉變換；Wigner-Ville分布；時(shí)頻分析；地震波譜E-mail:xikangliu@126.com，電話文收到日期：2015-11-020引言根據(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í)頻分析或稱時(shí)頻分布，是描述信號(hào)頻率隨時(shí)間變化的方法，采用時(shí)間一頻率聯(lián)合表示信號(hào)，將一維時(shí)域信號(hào)映射到二維時(shí)頻平面，在時(shí)頻域內(nèi)對(duì)信號(hào)進(jìn)行分析，全面觀測(cè)信號(hào)的時(shí)間一頻率聯(lián)合特征，同時(shí)掌握信號(hào)的時(shí)域及頻域信息（趙淑紅，2010；劉?？档?，2013）。短時(shí)傅里葉變換（short-timeFouriertransform，STFT）是一種線性變換，對(duì)多分量的非平穩(wěn)信號(hào)不會(huì)產(chǎn)生交叉干擾項(xiàng)，但由于采用固定窗，使得時(shí)間和頻率分辨率難以同時(shí)達(dá)到最優(yōu)，時(shí)頻分辨率較差。Wigner-Ville分布（Wigner-Villedistribution，WVD）具有理論上的最高時(shí)頻分辨率和許多優(yōu)良的數(shù)學(xué)性質(zhì)（王見等，2013），但因其為雙線性變換，對(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í)傅里葉變換線性變換無(wú)交叉項(xiàng)的特性，抑制Wigner-Ville分布交叉項(xiàng)干擾。呂宙等（2009）主要應(yīng)用于航天電子對(duì)抗雷達(dá)線性調(diào)頻信號(hào)檢測(cè)研究；趙淑紅（2010）利用二者直接相乘方法，提取模擬線性調(diào)頻地震波信號(hào)的瞬時(shí)頻率；王見等（2013）利用二者聯(lián)合方法進(jìn)行機(jī)械轉(zhuǎn)輪故障檢測(cè)研究。本文利用二者聯(lián)合表示方法，用于地震波時(shí)頻分析和縱波到時(shí)檢測(cè)。1算法原理WVD的定義假設(shè)實(shí)際信號(hào)x(t)的解析信號(hào)為s(t)，其WVD定義為式中，是s(t)的復(fù)共扼，將一維時(shí)間信號(hào)x(t)變換到二維時(shí)頻平面Ws(t,f)。使用解析信號(hào)的優(yōu)點(diǎn)是，可以消除信號(hào)傅里葉變換的負(fù)頻值，使其負(fù)頻率部分的頻譜值為0，從而得到正頻率頻譜的平均頻率和頻率擴(kuò)展，還可以消除時(shí)頻分布中負(fù)頻率區(qū)成分和正負(fù)頻率區(qū)之間的交叉干擾項(xiàng)(張帆等，2006)。WVD是一種雙線性變換，不滿足疊加原理。對(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)，采樣率為1Hz，振幅均為1。在32s處，兩個(gè)高斯元信號(hào)頻率中心分別為0.15Hz和0.35Hz，持續(xù)時(shí)間20s；在96s處，兩個(gè)高斯元信號(hào)頻率中心分別為0.15Hz和0.35Hz，持續(xù)時(shí)間也為20s。圖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；劉?？档龋?013)，統(tǒng)稱為Cohen類時(shí)頻分布。其中，核函數(shù)，Wx(t,f)為信號(hào)x(t)的WVD。式(3)說(shuō)明，核函數(shù)相當(dāng)于Wigner-Ville分布的一種加窗，雖然平滑了交叉項(xiàng)影響，但同時(shí)降低了WVD的時(shí)頻聚集性。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í)傅里葉變換是一種線性變換，對(duì)多分量信號(hào)來(lái)說(shuō)，不存在交叉項(xiàng)。由不確定原理可知，短時(shí)傅里葉變換(短時(shí)傅里葉變換英文縮寫STFT，為了描述連貫采用中文全稱)的時(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差得多。STFT與WVD的聯(lián)合算法從STFT和WVD各自優(yōu)缺點(diǎn)出發(fā)，可以得出，由于STFT是線性變換，不存在交叉干擾項(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)合算法可定義（王見等，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)乘，a、B稱為幕調(diào)節(jié)系數(shù)。該式的作用是增強(qiáng)兩變換數(shù)值較大部分而消弱有交叉項(xiàng)部分?！汶S著a、B的增大，STFT和WVD聯(lián)合變換的時(shí)頻聚集性會(huì)提高，在實(shí)際應(yīng)用中，需靈活調(diào)試a、B取值，以免丟失某些有用信息。式（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)站距震中約102km，采用CMG-3ESP型寬頻帶地震儀，采樣頻率100Hz。由于聯(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)為50s的地震波形進(jìn)行時(shí)頻分析（假設(shè)截取波形開始時(shí)刻為0秒）［圖2（a）］,能量譜密度見圖2（b）,可以發(fā)現(xiàn)，其能量主要集中在小于5Hz的頻率范圍內(nèi)，少部分在5—10Hz，說(shuō)明該記錄為低于5Hz頻率成分為主的地震波。圖2（c）為莆田ML4.8地震波形的Wigner-Ville分布，能夠較好反映該地震波的時(shí)頻分布特征，頻率和能量譜密度分布一致，15s前后橫波到達(dá)，能量達(dá)到最強(qiáng)，由于縱波能量相對(duì)較弱，時(shí)頻域內(nèi)縱波未達(dá)到最大幅值的5%以上，圖中沒(méi)有很好地顯示出來(lái)。另外，由于地震波的非平穩(wěn)特性，受噪聲和自項(xiàng)相互干擾等影響，不難發(fā)現(xiàn)，時(shí)頻圖像存在較多干擾成分，與有用信息混合在一起。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波前后共15s的三分量波形進(jìn)行分析（假設(shè)波形開始時(shí)刻為0秒），圖3（a）從左向右依次為N、E、U三方向，虛豎線表示P波初動(dòng)時(shí)刻，在該時(shí)段內(nèi)，初動(dòng)時(shí)間讀取為3.85s。圖3（b）為利用二者聯(lián)合算法，由式（8）所得結(jié)果，幕調(diào)節(jié)系數(shù)取a=1.5,P=1,可以發(fā)現(xiàn)，聯(lián)合算法的時(shí)頻聚集性較好，同時(shí)消除了交叉項(xiàng)干擾，地震記錄頻率范圍主要集中在1—5Hz，少量分布在5—10Hz,在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—10Hz。圖3(d)為三方向的Wigner-Ville分布，可以看出，其能量主要集中在5Hz以下，并且有從低頻向高頻變化的趨勢(shì)，能夠判讀P波初動(dòng)時(shí)間，且與該時(shí)段內(nèi)波形記錄文件讀取到的P波到時(shí)一致，為3.85s，這也為精確檢測(cè)P波初動(dòng)時(shí)刻提供了一種方法，但整體看由于噪聲干擾和交叉項(xiàng)干擾，有效信號(hào)被淹沒(méi)，存在虛假高頻信息，時(shí)頻分辨率并不理想。通過(guò)以上結(jié)果對(duì)比可以發(fā)現(xiàn)，STFT和WVD聯(lián)合算法既保留了STFT的線性變換沒(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)，窗口大小及類型的選取需結(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í)別(沈萍等，1999X縱波到時(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.王見，李金同，盧華玲，等.采用STFT-Wigner變換抑制Wigner-Ville分布交叉項(xiàng)[J].重慶大學(xué)學(xué)報(bào)，2013,36(8):15-25.張帆，鐘羽云，朱新云，等.時(shí)頻分析方法及地震波譜研究中的應(yīng)用[J].地震地磁觀測(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í)頻