蔣禮

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

運(yùn)用分頻段希爾伯特黃變換進(jìn)行多分量信號(hào)的頻散分析?

蔣禮1，2

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

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

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

1 引言

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

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

2 原理簡介

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

2.1 經(jīng)驗(yàn)?zāi)B(tài)頻率分解

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

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

假設(shè)某一信號(hào)s（t）存在K個(gè)分量，通過經(jīng)驗(yàn)?zāi)B(tài)分解后可得

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

濾波所得信號(hào)的時(shí)間序列di（t）稱為一個(gè)經(jīng)驗(yàn)頻率模態(tài)函數(shù)。上述的整個(gè)過程稱為一個(gè)經(jīng)驗(yàn)?zāi)B(tài)頻率分解過程。對(duì)該多分量信號(hào)通過完整的經(jīng)驗(yàn)?zāi)B(tài)頻率分解可得

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

2.2 時(shí)頻分析過程

經(jīng)驗(yàn)?zāi)B(tài)頻率分解的數(shù)學(xué)原理雖然十分簡單，但是按照上述的推導(dǎo)流程來實(shí)現(xiàn)整個(gè)分解過程將會(huì)存在下列難點(diǎn)：

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

（2）雖然可以在進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)頻率分解前先運(yùn)用經(jīng)驗(yàn)?zāi)B(tài)分解來確定濾波器的數(shù)目及參數(shù)，但這樣會(huì)浪費(fèi)分解時(shí)間，降低分解效率。

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

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

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

2.3 頻散分析過程

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

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

3 仿真實(shí)例

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

3.1 模型建立

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

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

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

圖3 信號(hào)1Fig.3 Signal1

圖4 信號(hào)2Fig.4 Signal2

3.2 譜分析

3.2.1 希爾伯特黃變換

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

圖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 分頻段希爾伯特黃變換

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

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

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

表2 實(shí)測(cè)時(shí)間-頻率-相位Table 2 Measured values of time-frequency-phase

3.3 計(jì)算結(jié)果分析

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

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

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

3.4 頻散曲線

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

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

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

4 結(jié)論

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

使用該方法時(shí)有兩個(gè)細(xì)節(jié)需要引起注意。

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

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

運(yùn)用譜分析方法計(jì)算多分量信號(hào)的頻散特性曲線關(guān)鍵在于兩點(diǎn)：準(zhǔn)確的時(shí)頻分析和精確的速度計(jì)算。從本文和文獻(xiàn)［3］的對(duì)比可以看出，S-HHT與傳統(tǒng)的HHT相比雖然增加了計(jì)算量，但提高了頻域分辨能力及抗噪性能。當(dāng)待分析信號(hào)較為復(fù)雜或信噪比較差時(shí)，使用S-HHT能夠獲得更準(zhǔn)確的時(shí)頻譜。而同步相差和異步相差算法，在進(jìn)行相速度計(jì)算時(shí)精度相似，都是進(jìn)行頻散分析的優(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］蔣禮.運(yùn)用希爾伯特黃變換進(jìn)行多分量信號(hào)的頻散分析［J］.電訊技術(shù)，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］程乾生，武連文.時(shí)間序列的經(jīng)驗(yàn)?zāi)B(tài)頻率分解EMFD［J］.數(shù)學(xué)的實(shí)踐與認(rèn)識(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

國家自然科學(xué)基金資助項(xiàng)目（40974079）

TN911

A

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

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

Email：jlncwu＠126.com

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