王衛(wèi)東，鄭 怡，張永志，趙云峰

(1.長安大學地質(zhì)工程與測繪學院，陜西西安710054;2.中國地震局第二監(jiān)測中心，陜西西安710054)

0 引言

地震波波形是具有時變特性的非線性、非平穩(wěn)信號，使得對地震波的分析方法受到某些限制。傳統(tǒng)上，傅立葉頻譜提供了一個簡單的方法求得“能量—頻率”分布，但它要求系統(tǒng)必須是線性的，而且資料必須是周期性的或是平穩(wěn)的，否則求得的頻譜就沒有物理意義。近十幾年來，使用最多的是小波分析法(Basu，Gupta，2001;Shan et al，2009;Spanos et al，2005;曹暉等，2002)，可進一步求得“能量—頻率—時間”分布，然而本質(zhì)上，小波分析法由在傅立葉頻譜分析法的一些假設(shè)基礎(chǔ)上推論而來，仍然具有原方法的限制與缺點。為了突破傅立葉分析法的限制，Huang等(1998)年提出了一種非線性、非平穩(wěn)信號分解方法——經(jīng)驗模態(tài)分解(Empirical Mode Decomposition method，簡稱EMD)，該方法基于信號本身的時間尺度特征，把復雜的信號分解為有限個固有模態(tài)分量(Intrinsic Mode Functions，簡稱IMF)和一個余項，每一項IMF分量都適用希爾伯特變換，由此可計算瞬時頻率，因此就可同時在頻率軸及時間軸上討論一事件，稱為希爾伯特—黃變換方法(Hilbert-Huang Transform method，簡稱 HHT)。HHT是一種自適應的信號處理方法，比傳統(tǒng)的傅立葉分析和依賴于先驗基的小波分析等方法更適用于非線性、非平穩(wěn)信號分析(Huang et al，1998;Wu，Huang，2009;)。目前，此方法已廣泛應用于目標識別、結(jié)構(gòu)損傷診斷、信號去噪、地震信號(Li，Chi，2008;Lin，Wang，2006;方嘉治等，2012;錢昌松等，2010;王彬等，2005)和爆破信號分析等方面(李夕兵等，2005，2006;張義平等，2005)。

筆者利用小波變換、S變換、Wigner-ville分布和Hilbert-Huang變換(HHT)對模擬信號和陜西數(shù)字地震臺網(wǎng)記錄的寬頻帶爆破振動信號進行了時頻分析，比較了這些方法的優(yōu)缺點和適用范圍。

1 時頻分析方法

1.1 小波變換

小波變換是一種窗口面積固定但形狀可以改變，即時間窗和頻率窗都可改變的時頻局部化分析方法。對于任意的函數(shù)或信號，小波變換定義為

其中，Ψ(x)為基本小波函數(shù)或母小波函數(shù)，在實際應用中需要選擇小波基，所選小波基函數(shù)在具有緊支撐性和一定正則性的同時，還需外形與被分析的對象有一定的相似性，本文選用Morlet小波。

1.2 S變換

S變換是以Morlet小波為基本小波的連續(xù)小波變換的延伸。在S變換中基本小波由簡諧波與高斯函數(shù)的乘積構(gòu)成，基本小波中的簡諧波在時間域僅作伸縮變換，而高斯函數(shù)則進行伸縮和平移，故可認為是連續(xù)小波變換的“相位校正”。函數(shù)的S變換表示為

S變換以高斯視窗為基底函數(shù)來回移動，可以分析時頻域特征，雖然直接和傅立葉相關(guān)，但S變換的優(yōu)點是其局部化的高斯視窗可以將信號詳細轉(zhuǎn)換，提高了在頻率域的解析度(Pinnegar，2005;Stockwell et al，1996;樊劍等，2008)。

1.3 Wigner-Ville分布(WVD)

Wigner-Ville分布首先由Wigner提出，用于量子力學領(lǐng)域問題研究，后由Ville引入到信號分析領(lǐng)域。它有很多優(yōu)良的數(shù)學性質(zhì)，且表達式直觀簡單。對于信號，其Wigner-Ville分布為(科恩，1998)

1.4 Hilbert-Huang變換

Hilbert-Huang變換首先用經(jīng)驗模態(tài)分解(EMD)，將信號分解為若干IMF，然后進行Hilbert變換，將實數(shù)信號配上虛部，得到解析信號，再根據(jù)解析信號的相位隨時間的改變而得到瞬時頻率。

EMD方法根據(jù)數(shù)據(jù)本身的特征進行分解，將信號分解成一組穩(wěn)態(tài)的數(shù)據(jù)序列集，即IMF分量模式。EMD方法通過特征時間尺度來獲得IMF分量篩分過程，將原始信號分解成n個imf分量和一個殘余分量，將信號從高頻到低頻的順序依次分解。即:

對于任一IMF分量Cj(t)作Hilbert變換，得到一個復IMF分量Zj(t)，aj(t)表示復IMF的振幅，θj(t)表示相位，瞬時頻率因此可得到Hilbert譜為:

可見，HHT譜精確地描述了信號的幅值隨時間和頻率的變化規(guī)律，HHT方法具有較好的自適應性和頻率瞬時性，對信號的非線性反映能力較好且對穩(wěn)態(tài)和非穩(wěn)態(tài)信號都能進行分析，適合分析具有非線性和非平穩(wěn)動態(tài)變化特征的信號。

2 仿真實驗

假設(shè)模擬信號s為4個余弦疊加信號，其表達式為

其中，f1=1 Hz，f2=10 Hz，f3=20 Hz，f4=50 Hz，數(shù)據(jù)采樣頻率為500 Hz，樣本數(shù)據(jù)為1000個，對 s(t)進行 EMD分解，得到結(jié)果如圖1所示。

圖1 信號s(t)及其EMD分解Fig.1 Simulative signal s(t)and its EMD components

由圖1可見，C0為原始信號，C1、C2、C3、C4為用EMD分解得到的4個分量。可見，EMD分解是信號本身所決定的一個自適應分解過程，能很快地提取信號特征并分解出信號的分量，分解出的4個IMF分量正是仿真信號的4個原始信號，表明了分解的可靠性、高效性，體現(xiàn)了IMF分量本身的物理含義。

圖2為信號s(t)的HHT頻譜圖、S變換頻譜圖、小波變換頻譜圖和WVD圖。由圖可見，HHT準確地提取了信號的時頻特征，通過S變換亦較好地提取了信號的特征，較好地分解出了仿真信號的4個原始信號，但其分辨率低于HHT，小波變換則存在漏能現(xiàn)象，頻率分辨率明顯低于HHT和S變換，Wigner-Ville分布具有良好的時頻聚集性，但由于交叉干擾項的影響，出現(xiàn)了虛假頻率成分。

圖2 信號s(t)轉(zhuǎn)換頻譜圖(a)HHT頻譜;(b)S變換頻譜;(c)小波變換頻譜;(d)WVD圖Fig.2 Transform spetrum of singnal s(t)(a)HHT spectrum;(b)S-transform spectrum;(c)Wavelet transform spectrum;(d)WVD

3 爆破振動信號時頻分析

實際記錄取自華陰地震臺記錄到的爆破振動信號，采用FBS-3B型寬頻帶數(shù)字地震計，頻帶范圍為0.05～20 s，采樣率為每秒50個樣點。圖3為爆破振動信號的原始曲線(C0)和EMD分解圖(C1～C6)。可見，EMD分解較好地分離了P波段和Rg波段，爆破振動信號可分為4個變化比較明顯的階段:①爆破發(fā)生前的地脈動記錄階段。此時爆破還沒有發(fā)生，在時頻圖上能量譜值接近于0。②P波階段。能量從無到有，特別對P波初動的突變性，被精確地刻劃出來，但能量值不是很大。③Rg波階段。這個階段是一個漸變過程，Rg波初動的突變性被較精確地刻劃出來，能量遠大于其它部分，是地震能量釋放的主要階段。④ 尾波持續(xù)階段。含有一定能量，相對較弱，變化不穩(wěn)，可能是介質(zhì)散射造成的。⑤ 爆破振動結(jié)束階段。爆破振動能量基本消失，又回到正常地脈動記錄階段。

圖3 爆破振動信號及其EMD分解Fig.3 Blasting vibration signal and its EMD components

圖4為爆破振動信號的HHT時頻圖。由圖可見，HHT時頻圖清晰而詳細地顯示了能量隨時頻的具體分布。可較清晰地識別爆破振動信號在時頻域中的幾個峰值，最大峰值頻率為2.5 Hz，發(fā)生的時間是7.5 s左右，峰值頻率為1Hz和6Hz的發(fā)生的時間在是9 s左右，峰值頻率為10 Hz信號的發(fā)生時間在是2.5 s左右，峰值頻率為15Hz信號的發(fā)生時間則在2.5、7.5和9 s左右;還可看出，大致在17 s后，低頻段的能量超過了高頻段，在25 s，高頻段的能量顯著減弱。

圖4 爆破振動信號的HHT時頻圖Fig.4 HHT spectrum of the blasting vibration signal

4 討論與結(jié)論

利用模擬信號，對小波變換、S變換、Wignerville分布(WVD)和Hilbert-Huang變換(HHT)的時頻分析進行了比較，結(jié)果表明基于經(jīng)驗模態(tài)分解的HHT具有更強的適應性，無需事先確定分解階次，不受人為因素影響，能較為準確地提取信號的時變特征。利用HHT對實際記錄的寬頻帶爆破振動信號進行了時頻分析，其結(jié)果較為準確地描述了信號的時變特征。

曹暉，賴明，白紹良.2002.基于小波變換的地震地面運動仿真研究［J］.土木工程學報，35(4):41－46.

樊劍，呂超，張輝.2008.基于S變換的地震波時頻分析及人工調(diào)整［J］.振動工程學報，21(4):381－386.

方嘉治，趙志偉，蔡宗文，等.2012.古田地震水口水電站重力壩強震觀測記錄HHT分析［J］.地震研究，35(2):236－239.

李夕兵，張義平，劉志祥，等.2005.爆破震動信號的小波分析與HHT變換［J］.爆炸與沖擊，25(6):528－534.

李夕兵，張義平，左宇軍，等.2006.巖石爆破振動信號的EMD濾波與消噪［J］.中南大學學報(自然科學版)，37(1):150－154.

錢昌松，劉代志，劉志剛，等.2010.基于遞歸高通濾波的經(jīng)驗模態(tài)分解及其在地震信號分析中的應用［J］.地球物理學報，53(5):1215－1225.

王彬，楊潤海，郭夢秋，等.2005.昆明高層建筑強震觀測記錄HHT分析［J］.地震研究，28(1):78－81.

張義平，李夕兵，趙國彥.2005.基于HHT方法的爆破地震信號分析［J］.工程爆破，11(1):l－7.

Basu B，Gupta V K.2001.Wavelet-based stochastic seismic response of a duffing oscillator［J］.Journal of Sound and Vibration，245(2):251－260.

Huang N E，Shen Z，Long S R.1998.The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceeding of the Royal Society of London，Series A，454:903－995.

科恩L.1998.時頻分析:理論與應用［M］.白居憲，譯.西安:西安交通大學出版社.

Li X，Chi H.2008.Echo analysis of objects on the seafloor［J］.The Journal of the Acoustical Society of America，123(5):4803 －4807.

Lin Z S，Wang S G.2006.EMD analysis of solar insolation［J］.Meteorol Atmos Phys，93:123 －128.

Pinnegar C R.2005.Time-frequency and time-time filtering with the S-transform and TT-transform［J］.Digital Signal Processing，15(5):604－620.

Shan H，Ma J，Yang H.2009.Comparisons of wavelets，contourlets，andcurvelets in seismic denoising［J］.Journal of Applied Geophysics，69(2):103－115.

Spanos P D，Tezcan J，Tratskas P.2005.Stochastic processes evolutionary spectrum estimation via harmonic wavelets［J］.Computer Methods and Application Mechanics Engineering，194(12):1 367 －1 383.

Stockwell R G，Mansinha L，Lowe R P.1996.Localization of the complex spectrum:The S transform［J］.IEEE Transactions on Signal Processing，44(4):998 －1001.

Wu Z，Huang N E.2009.Ensemble empirical mode decomposition:noise assisted data analysis method［J］.Advance in Adaptive Data Analysis，1:1 － 41.