運用分頻段希爾伯特黃變換進行多分量信號的頻散分析?

蔣禮

（1.中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢430074；2.華北水利水電學院數(shù)學與信息科學學院，鄭州450011）

運用分頻段希爾伯特黃變換進行多分量信號的頻散分析?

蔣禮1，2

（1.中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢430074；2.華北水利水電學院數(shù)學與信息科學學院，鄭州450011）

將分頻段希爾伯特黃變換應用于多分量信號的頻散分析中。首先，利用帶通濾波器和經(jīng)驗模態(tài)分解相結(jié)合，成功實現(xiàn)了經(jīng)驗模態(tài)頻率分解，并準確提取了經(jīng)驗頻率模態(tài)函數(shù)；然后，使用該方法準確地獲取了多分量含噪信號的時頻能量譜和時頻相位譜；最后，基于同步相差和異步相差算法，精確繪制了原信號的相速度頻散分析曲線。數(shù)值試驗表明，該算法擁有較高的時頻分辨能力和良好的抗噪性能，對于復雜且信噪比較低的信號，能獲得比傳統(tǒng)希爾伯特黃變換更準確的頻散分析結(jié)果。

希爾伯特黃變換；經(jīng)驗模態(tài)頻率分解；經(jīng)驗頻率模態(tài)函數(shù)；瞬時頻率；頻散曲線

1 引言

基于瞬時頻率的希爾伯特黃變換（Hilbert-Huang Transform）是一種進行時頻分析的有效方法［1］，因為該方法可以準確刻畫信號的瞬時頻率和瞬時相位，所以它成為進行多分量信號頻散分析的有效工具［2-3］。但由于傳統(tǒng)希爾伯特黃變換的核心是基于時間序列的經(jīng)驗模態(tài)分解，該方法在擁有較高的時間分辨能力的同時往往頻率分辨能力較差［4］。雖然對于普通的多分量信號，使用該方法能夠取得令人滿意的結(jié)果，但是對于所含頻率分量較多的情況，使用該方法往往不能令人滿意。

本文對傳統(tǒng)的希爾伯特黃變換進行改進，運用經(jīng)驗模態(tài)頻率分解理論［5］代替?zhèn)鹘y(tǒng)的經(jīng)驗模態(tài)分解，從而提出分頻段希爾伯特黃變換。該方法的頻率分辨能力較傳統(tǒng)的希爾伯特黃變換有明顯提高。將該方法用于進行多分量信號的頻散分析，可以得到令人滿意的結(jié)果。

2 原理簡介

傳統(tǒng)希爾伯特黃變換的核心為經(jīng)驗模態(tài)分解（EMD）和固有模態(tài)函數(shù)（IMF）。由于經(jīng)驗模態(tài)分解進行的是一種時域濾波，所以時間分辨率通常是令人滿意的；但由于分解后的固有模態(tài)函數(shù)個數(shù)有限，故頻率分辨率往往較低。對于頻率成分較為豐富的信號直接進行希爾伯特黃變換通常得不到理想的結(jié)果。

2.1 經(jīng)驗模態(tài)頻率分解

分頻段希爾伯特黃變換可以有效解決上述問題。分頻段希爾伯特黃變換的核心為經(jīng)驗模態(tài)頻率分解（EMFD）［5］，它是由北京大學程乾生教授于2005

年提出的一種時頻分解理論，初始目的是為了給希爾伯特黃變換提供一組正交基。但在實際應用的過程中，經(jīng)此方法分解所得的經(jīng)驗頻率模態(tài)函數(shù)（EFMF）不僅正交，而且也具有比固有模態(tài)函數(shù)更高的頻率分辨能力。

假設某一信號s（t）存在K個分量，通過經(jīng)驗模態(tài)分解后可得

式中，每個ci（t）稱為一個固有模態(tài)函數(shù)。設每個固有模態(tài)函數(shù)的頻譜為Ci（f），由于希爾伯特黃變換的頻域分辨率較差，雖然Ci（f）和Ci+1（f）的主頻并不相同，但兩者在頻譜卻有不少重疊的部分［6-7］（這也正是傳統(tǒng)的固有模態(tài)函數(shù)不是正交基的原因）。假設兩者的主頻分別位于fi的兩側(cè)，則設定一個理想帶通濾波器：Hi（f）=1，fi≤f≤fi-1，其余頻段Hi（f）=0，將該濾波器作用于原信號，可得

濾波所得信號的時間序列di（t）稱為一個經(jīng)驗頻率模態(tài)函數(shù)。上述的整個過程稱為一個經(jīng)驗模態(tài)頻率分解過程。對該多分量信號通過完整的經(jīng)驗模態(tài)頻率分解可得

此時不僅由于分解所得的di（t）頻譜部分沒有重疊，從而構(gòu)成一組正交基；而且因為分解過程中K個帶通濾波器的作用，所以可以準確無誤地提取每個頻率分量，故該方案具有極高的頻率分辨能力。

2.2 時頻分析過程

經(jīng)驗模態(tài)頻率分解的數(shù)學原理雖然十分簡單，但是按照上述的推導流程來實現(xiàn)整個分解過程將會存在下列難點：

（1）對于一個未知的信號，由于無法預知其頻率分量，從而對于濾波器的數(shù)目及參數(shù)難以確定；

（2）雖然可以在進行經(jīng)驗模態(tài)頻率分解前先運用經(jīng)驗模態(tài)分解來確定濾波器的數(shù)目及參數(shù)，但這樣會浪費分解時間，降低分解效率。

為了高效實現(xiàn)經(jīng)驗模態(tài)頻率分解過程，可首先在整個頻段上設置一系列的頻點，將相鄰兩個頻點組成帶通濾波器；然后將待分析的信號通過每個帶通濾波器可以得到一系列的頻段信號；最后對每個分頻段信號進行經(jīng)驗模態(tài)分解，可以得到每個頻段的固有模態(tài)函數(shù)。這時每個頻段的固有模態(tài)函數(shù)實際上就是經(jīng)驗頻率模態(tài)函數(shù)（每個經(jīng)驗頻率模態(tài)函數(shù)的頻譜不會有重疊的部分），而對于分頻段信號的經(jīng)驗模態(tài)分解過程實際上就是經(jīng)驗模態(tài)頻率分解過程。

完整的分頻段希爾伯特黃變換實現(xiàn)流程如圖1所示。與數(shù)學推導流程的不同之處在于，數(shù)學推導是先進行經(jīng)驗模態(tài)分解，而后再制定帶通濾波器濾波；而此處對于未知信號是先做帶通濾波再做經(jīng)驗模態(tài)分解，但兩者最終所得結(jié)果完全一樣（K個固有模態(tài)函數(shù)即為原信號的K個頻率分量，而且頻譜部分不會重疊完全正交）。

圖1 分頻段希爾伯特黃變換算法流程Fig.1 The algorithm flow of sub-band Hilbert-Huang transform

2.3 頻散分析過程

對原始多分量信號進行分頻段希爾伯特黃變換可以得到時頻能量譜和時頻相位譜［3］；對上述兩譜進行分析，可以獲得任意分量信號的瞬時頻率、瞬時相位、瞬時能量、持續(xù)時間等信息［3］。根據(jù)已獲知的信息運用同步相差和異步相差理論可以求取任意頻率分量的相速度值［3］。完整的頻散分析過程如圖2所示。

圖2 運用分頻段希爾伯特黃變換進行頻散分析流程Fig.2 The algorithm flow of dispersion analysis by sub-band Hilbert-Huang transform

3 仿真實例

為了驗證分頻段希爾伯特黃變換的正確性以及較傳統(tǒng)希爾伯特黃變換的優(yōu)越性，此處設定一個較為復雜的多分量信號模型。

3.1 模型建立

假設瑞利波源由10個分量信號疊加而成，初始信號持續(xù)時間為0.2 s；各分量信號的群速度和相速度隨著頻率的增加而減小，且群速度小于相速度；信號能量集中在50 Hz、60 Hz，隨著頻率的增加和減少，能量依次遞減。各分量信號的具體參數(shù)詳見表1。

表1 瑞利波源分量參數(shù)列表Table 1 Rayleigh wave source component parameter list

兩個檢波器分別距離震源12 m和18m，檢波時間為1.023 s，采樣時間間隔為1ms。為了檢驗該算法的抗噪性能，檢波信號除了添加10%的高斯白噪聲模擬檢波干擾外，再整體添加5%的高斯白噪聲作為環(huán)境噪聲。12m處獲取的信號如圖3所示，18m處獲取的信號如圖4所示。

圖3 信號1Fig.3 Signal1

圖4 信號2Fig.4 Signal2

3.2 譜分析

3.2.1 希爾伯特黃變換

運用傳統(tǒng)希爾伯特黃變換所得譜分析結(jié)果如圖5所示。圖5（a）和（b）分別為12m和18m處信號的時間-頻率-能量譜；而圖5（c）和（d）則分別為12m和18m處信號的時間-頻率-相位譜。從上述各圖可以發(fā)現(xiàn)，無論是時頻能量譜還是時頻相位譜都不能準確刻畫出10個頻率分量的完整信息。雖然由于希爾伯特黃變換良好的時間分辨能力，對比12m處和18m處從高到低的各頻率分量能量譜可以清晰地看出頻散效應；但是由于不能將10個頻率分量的時頻信息完整地反映出來，故而無法求取每個分量的準確相速度值。

圖5 傳統(tǒng)希爾伯特黃變換的譜分析圖Fig.5 Spectrum analysis diagram of Hilbert-Huang transform

圖6 分頻段希爾伯特黃變換的譜分析圖Fig.6 Spectrum analysis diagram of the sub-band Hilbert-Huang transform

3.2.2 分頻段希爾伯特黃變換

分頻段帶通濾波器必須預設兩個參數(shù)：有效頻段和濾波器帶寬。有效頻段反映的是波源信號的物理特性，而帶寬則反映了頻率分辨能力。理論情況下，無窮大的有效頻段和無限小的帶寬可以滿足任何多分量信號的分解；但在實際應用中，有效頻段設定越小，帶寬設定越大，則運算效率越高。此處有效頻段設為5～105 Hz，以10 Hz為帶寬進行窄帶濾波。

運用分頻段希爾伯特黃變換所得譜分析結(jié)果如圖6所示。從上述各圖可以發(fā)現(xiàn)，無論是在時頻能量譜還是時頻相位譜中，10個頻率分量都顯示得十分清晰，而且可以看出各分量的瞬時頻率非常接近于真實的物理頻率，延遲時間、持續(xù)時間和能量分布也基本正確。

從圖6（a）和（b）中可以看出12m處及18m處的各分量信號均通過0.2 s；18m處的信號各分量均通過0.25 s。為了分析方便，只測量這兩個時刻的瞬時頻率和瞬時相位，所測量的數(shù)據(jù)如表2所示（表2中的平均頻率為某一分量所有時刻瞬時頻率的加權(quán)平均［3］）。從表2中可以看出，無論是瞬時頻率還是平均頻率都非常接近接近于該分量信號的真實物理頻率。

表2 實測時間-頻率-相位Table 2 Measured values of time-frequency-phase

3.3 計算結(jié)果分析

對表2所測得的數(shù)據(jù)進行相速度計算，所得結(jié)果及誤差如表3所示。從計算結(jié)果來看，所得相速度較為理想。第三分量至第十分量基本接近于相應模型分量的真實速度；第一分量和第二分量雖然誤差較大，但仍在可接受的范圍之內(nèi)。

表3 相速度及誤差計算結(jié)果Table 3 Results of phase velocity and error

上述兩個方法出現(xiàn)的誤差主要產(chǎn)生于分母部分的相差和頻率的比值項Δφ／ω0［3］。理論上Δφ／ω0的誤差會隨著相位解纏時頻率的升高而逐漸減小。

3.4 頻散曲線

根據(jù)表3的計算結(jié)果可以繪制4條頻散曲線。若在位于24m的測點處增加一個檢波器記錄，則可以繪制12條頻散曲線。其中8條是道間距為6m的頻散曲線，4條為道間距為12m的頻散曲線。所有的相速度頻散曲線如圖7所示。

圖7 相速度頻散曲線Fig.7 Phase velocity dispersion curves

從圖7中可以看出，頻散曲線正確反應了頻散模型的相速度變化規(guī)律。隨著頻率的升高，頻散曲線的收斂性越來越好；而且道間距為12m的頻散曲線誤差要整體小于道間距為6m的頻散曲線誤差。

4 結(jié)論

經(jīng)驗模態(tài)頻率分解在保留經(jīng)驗模態(tài)分解較高的時間分辨能力同時，又提高了頻率分辨能力，以其為基礎的分頻段希爾伯特黃變換具有極為準確的時頻-相位定位能力。將該方法用于進行多分量信號的頻散分析，不僅計算結(jié)果準確，而且具有較強的抗噪能力。

使用該方法時有兩個細節(jié)需要引起注意。

（1）進行分頻段處理時選擇合適的頻窗，不同窗型及窗口的大小會影響最終的結(jié)果。因為進行窄帶濾波必然會降低希爾伯特譜的時間分辨率。此時應注意，頻窗不應取得太小，窗型不應選擇劇烈跳變類型。

（2）如果對源信號各分量的頻率成分較為熟悉，則可采用頻段交叉的頻窗進行窄帶濾波，即在不減小窗口大小的情況下使窗口頻段交叉（頻窗的主頻部分不可交叉）。這樣可以在相同有效頻段的基礎上增加帶通濾波器的數(shù)量，即能夠在不降低時間分辨率的同時進一步提高頻率分辨率。

運用譜分析方法計算多分量信號的頻散特性曲線關鍵在于兩點：準確的時頻分析和精確的速度計算。從本文和文獻［3］的對比可以看出，S-HHT與傳統(tǒng)的HHT相比雖然增加了計算量，但提高了頻域分辨能力及抗噪性能。當待分析信號較為復雜或信噪比較差時，使用S-HHT能夠獲得更準確的時頻譜。而同步相差和異步相差算法，在進行相速度計算時精度相似，都是進行頻散分析的優(yōu)秀算法。

［1］Huang N E，Shen Zheng，Long SR.The empiricalmode decomposition and the Hilbertspectrum for nonlinear and nonstationary time series analysis［J］.Proceedings of the Royal Society of London，Series A-Mathematical Physical and Engineering Sciences，1998，454（1971）：903-995.

［2］Chen Chau-Huei，Li Cheng-Pling，Teng Ta-Liang.Surface-Wave Dispersion Measurements Using Hilbert-Huang Transform［J］.Terrestrial Atmospheric and Oceanic Sciences，2002，13（2）：171-184.

［3］蔣禮.運用希爾伯特黃變換進行多分量信號的頻散分析［J］.電訊技術，2011，51（7）：60-66. JIANG Li.Frequency Dispersion Analysis of Multi-component Signalusing Hilbert-Huang Transform［J］.Telecommunication Engineering，2011，51（7）：60-66.（in Chinese）

［4］Patrick F，GabrielR，Paulo G.EmpiricalMode Decomposition as a Filter Bank［J］.IEEE Signal Processing Letters，2004，11（2）：112-114.

［5］程乾生，武連文.時間序列的經(jīng)驗模態(tài)頻率分解EMFD［J］.數(shù)學的實踐與認識，2005，36（5）：151-153. CHEN Qian-sheng，WU Lian-wen.The EmpiricalMode Frequency Decomposition for Time Series Analysis［J］.Mathematics in Practice and Theory，2005，36（5）：151-153.（in Chinese）

［6］Huang N E，Shen SSP.Hilbert-Huang Transform and Its Application，Interdisciplinary Mathematical Sciences vol 5［M］.London：World Scientific Public Co Information，2005：57-74.

［7］WU Zhaohua，Huang N E.A Study of The Characteristics of White Noise Using The Empirical Mode Decomposition Method［J］.Proceedings of the Royal Society London，Series A-Mathematical Physical and Engineering Sciences，2004，460（2046）：1597-1611.

Analysis of Frequency Dispersion of M ulti-com ponent Signal using Sub-band Hilbert-Huang Transform

JIANG Li1，2
（1.Institute of Geophysics and Geomatics，China University of Geosciences（Wuhan），Wuhan 430074，China；2.College of Mathematics and Information Science，North China University ofWater Resources and Electric Power，Zhengzhou 450011，China）

In this paper，the sub-band Hilbert-Huang transform（S-HHT）is used to analyse the frequency dispersion ofmulti-component signal.Firstly，the band-pass filter and the empiricalmode decomposition are combined to achieve the empiricalmode frequency decomposition（EMFD）successfully，and accurately obtain the empirical frequencymode functions（EFMF）.Secondly，the time frequency energy spectrums and time frequency phase spectrums from multi-componentsignalwith noise are precisely drawn by themethod.Finally，by using the synchronous phase difference and asynchronous phase difference arithmetic the phase velocity dispersion curves are drawn accurately.The numerical results show that the S-HHTarithmetic hashigh time frequency resolution and good resistance to noise.For complex and low Signal-to-Noise Ratio（SNR）signal，the SHHT can get themore accurate results of dispersion analysis than the HHT.

Hilbert-Huang transform；empirical mode frequency decomposition；empirical frequency mode function；instant frequency；dispersion curves

The National Natural Science Foundation of China（No.40974079）

the B.S.degree and the M.S.degree in 2001 and 2004，respectively.He is now a lecturer and currently working toward the Ph.D.degree.His research concerns geoscience-signal analysis and processing

1001-893X（2012）04-0472-06

2011-12-21；

2012-03-09

國家自然科學基金資助項目（40974079）

TN911

A

10.3969／j.issn.1001-893x.2012.04.010

蔣禮（1979—），男，江蘇江陰人，2001年獲學士學位，2004年獲碩士學位，現(xiàn)為博士研究生，講師，主要從事地學信號的分析及處理。

Email：jlncwu＠126.com

JIANG Li was born in Jiangyin，Jiangsu Province，in 1979. He