王衛(wèi)東,鄭 怡,張永志,趙云峰
(1.長安大學地質(zhì)工程與測繪學院,陜西西安710054;2.中國地震局第二監(jiān)測中心,陜西西安710054)
地震波波形是具有時變特性的非線性、非平穩(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)缺點和適用范圍。
小波變換是一種窗口面積固定但形狀可以改變,即時間窗和頻率窗都可改變的時頻局部化分析方法。對于任意的函數(shù)或信號,小波變換定義為
其中,Ψ(x)為基本小波函數(shù)或母小波函數(shù),在實際應用中需要選擇小波基,所選小波基函數(shù)在具有緊支撐性和一定正則性的同時,還需外形與被分析的對象有一定的相似性,本文選用Morlet小波。
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)。
Wigner-Ville分布首先由Wigner提出,用于量子力學領(lǐng)域問題研究,后由Ville引入到信號分析領(lǐng)域。它有很多優(yōu)良的數(shù)學性質(zhì),且表達式直觀簡單。對于信號,其Wigner-Ville分布為(科恩,1998)
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)變化特征的信號。
假設(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
實際記錄取自華陰地震臺記錄到的爆破振動信號,采用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
利用模擬信號,對小波變換、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.