Hilbert-Huang變換與寬頻帶爆破信號(hào)的時(shí)頻分析*

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

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

0 引言

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

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

1 時(shí)頻分析方法

1.1 小波變換

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

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

1.2 S變換

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

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

1.3 Wigner-Ville分布(WVD)

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

1.4 Hilbert-Huang變換

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

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

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

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

2 仿真實(shí)驗(yàn)

假設(shè)模擬信號(hào)s為4個(gè)余弦疊加信號(hào)，其表達(dá)式為

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

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

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

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

圖2 信號(hào)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 爆破振動(dòng)信號(hào)時(shí)頻分析

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

圖3 爆破振動(dòng)信號(hào)及其EMD分解Fig.3 Blasting vibration signal and its EMD components

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

圖4 爆破振動(dòng)信號(hào)的HHT時(shí)頻圖Fig.4 HHT spectrum of the blasting vibration signal

4 討論與結(jié)論

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

曹暉，賴明，白紹良.2002.基于小波變換的地震地面運(yùn)動(dòng)仿真研究［J］.土木工程學(xué)報(bào)，35(4):41－46.

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

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

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

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

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

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

張義平，李夕兵，趙國(guó)彥.2005.基于HHT方法的爆破地震信號(hào)分析［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.時(shí)頻分析:理論與應(yīng)用［M］.白居憲，譯.西安:西安交通大學(xué)出版社.

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.