李宏偉，劉宇航，楊 輝

（1.阜新市水土保持工作總站，遼寧阜新123000；2.中水東北勘測(cè)設(shè)計(jì)研究有限責(zé)任公司科學(xué)研究院，吉林長(zhǎng)春130061）

Hilbert-Huang變換(Hilbert-Huang transform,HHT)[1，2]是 Norden E.Huang1998年提出的一種數(shù)據(jù)處理方法,經(jīng)驗(yàn)?zāi)Ｊ椒纸?(empirical mode decomposition,EMD)是該方法的第一步。HHT比傳統(tǒng)方法更適合分析非線形、非平穩(wěn)數(shù)據(jù)。現(xiàn)在在地震工程、結(jié)構(gòu)損傷、系統(tǒng)識(shí)別、音頻分析等領(lǐng)域都得到了廣泛的應(yīng)用。

1 EMD分解

EMD分解必須滿足2個(gè)假設(shè)：（1)在整個(gè)數(shù)據(jù)長(zhǎng)度中，極值點(diǎn)的數(shù)量(包括極大值點(diǎn)和極小值點(diǎn))和過零點(diǎn)的數(shù)量必須相等或至多相差一個(gè)。（2)在任一時(shí)間點(diǎn)上，信號(hào)局部極大值確定的上包絡(luò)線和局部極小值確定的下包絡(luò)線的均值必須為零。

分解方法用局部極大值和極小值的包絡(luò)來進(jìn)行。所有的局部極大值用三次樣條函數(shù)插值形成數(shù)據(jù)的上包絡(luò)，所有的局部極小值通過插值形成數(shù)據(jù)的下包絡(luò)，求出上包絡(luò)和下包絡(luò)的平均值，根據(jù)過零點(diǎn)和極值點(diǎn)數(shù)目、包絡(luò)均值為零的條件，檢查原始數(shù)據(jù)減去平均值得到的是否是一個(gè)固有模態(tài)函數(shù)(intrinsic mode function，IMF)，如果不是，就對(duì)得到的數(shù)據(jù)再次進(jìn)行篩選，直至得到符合上述條件的IMF。原數(shù)據(jù)減去第一個(gè)固有模態(tài)函數(shù)當(dāng)作新的數(shù)據(jù)再進(jìn)行篩選，直至殘余項(xiàng)成為一個(gè)單調(diào)函數(shù)，分解停止，最后得到一系列的固有模態(tài)函數(shù)和一個(gè)殘余項(xiàng)。

隨著HHT研究的深入和應(yīng)用的拓展，還有一些關(guān)鍵的理論和實(shí)際問題需要完善和改進(jìn)。EMD中有幾個(gè)關(guān)鍵問題對(duì)分解效果有顯著影響，分別是：IMF的篩分終止條件判斷依據(jù)；邊界的端點(diǎn)處理問題；各個(gè)IMF是否正交。

2 Rilling算法

有效的IMF必須滿足在任一時(shí)間點(diǎn)上，信號(hào)局部極大值確定的上包絡(luò)線和局部極小值確定的下包絡(luò)線的均值為零，通常的做法是利用Huang提出的標(biāo)準(zhǔn)偏差法[1,2]。然而,這個(gè)判斷條件的大小直接影響 IMF 的穩(wěn)定性[3，4]。

2003年，Rilling等[5]對(duì)原始EMD分解中的濾波停止條件進(jìn)行了改進(jìn)，該方法的核心思想是在EMD 篩選中引入 3個(gè)參數(shù) θ1，θ2，α 作為判斷依據(jù)[6]。具體過程是：確定時(shí)程曲線x(t)的所有極值點(diǎn)，所有的局部極大值用三次樣條插值函數(shù)形成數(shù)據(jù)的上包絡(luò)線fmax(t)；所有的局部極小值通過插值形成數(shù)據(jù)的下包絡(luò)線fmin(t)，上包絡(luò)和下包絡(luò)的平均值仍記作m(t)。

計(jì)算模態(tài)幅值：

并引進(jìn)新的評(píng)估函數(shù)：

篩選停止條件有2個(gè)：

（1）σ(t)≤θ1的時(shí)間(在等步長(zhǎng)離散情況下為步長(zhǎng)數(shù))與全部持續(xù)時(shí)間之比至少達(dá)到滿足1-α成立(其中α常取為0.05)，其目的是為了保證大多數(shù)數(shù)據(jù)波動(dòng)的平均值足夠小。

（2）σ(t)≤θ2，其目的是為了限制 EMD 分解時(shí)局部數(shù)據(jù)最大漂移。

與Huang的方法相比，σ(t)更能反映IMF的均值特性，且兩個(gè)相互補(bǔ)充，使信號(hào)只能在某些局部出現(xiàn)較大波動(dòng)，從而保證整體均值為零[7]。

為考查上述方法的效果，采用一個(gè)簡(jiǎn)單的解析信號(hào)：

分別用標(biāo)準(zhǔn)偏差法和Rilling濾波方法進(jìn)行計(jì)算比較。取 θ1＝0.05，θ2＝10θ1，α＝0.05，時(shí)間取 50s。

從平均頻率和平均振幅兩個(gè)方面衡量?jī)煞N方法的優(yōu)劣。Rilling算法得到的IMF和標(biāo)準(zhǔn)偏差法得到的IMF頻率與原始波頻率并沒有明顯的差異，但是從振幅來看，Rilling算法的振幅和原始波更為接近，特別在c3頻段，兩者的差別較大，所以Rilling的控制方法是優(yōu)于標(biāo)準(zhǔn)偏差法的。

3 多項(xiàng)式擬合求端點(diǎn)的極值

EMD分解中，通過對(duì)原數(shù)據(jù)中的上極值點(diǎn)和下極值點(diǎn)分別進(jìn)行三次樣條插值擬合再求平均得到包絡(luò)均值，在樣條插值時(shí)常常不能確定端點(diǎn)處的極值，就會(huì)在樣條插值時(shí)產(chǎn)生數(shù)據(jù)的擬合誤差。不斷累積的誤差就會(huì)由端點(diǎn)處向內(nèi)逐漸傳播，信號(hào)的兩端會(huì)出現(xiàn)嚴(yán)重的扭曲。

3.1 多項(xiàng)式擬合算法

端點(diǎn)常常并不是極值，如果能夠根據(jù)極值點(diǎn)序列中端點(diǎn)以內(nèi)數(shù)據(jù)的規(guī)律得出該序列在端點(diǎn)處的近似取值,則可以防止對(duì)端點(diǎn)插值出現(xiàn)較大的擺動(dòng)。

取出原極值點(diǎn)序列最左端的3個(gè)極值點(diǎn)(如果極值點(diǎn)序列的個(gè)數(shù)小于3個(gè)則取序列中所有元素),對(duì)所取的極值點(diǎn)求出擬合多項(xiàng)式,計(jì)算出多項(xiàng)式對(duì)應(yīng)數(shù)據(jù)序列左端點(diǎn)處的函數(shù)值,并把此函數(shù)值作為極值點(diǎn)序列在該端點(diǎn)處的近似取值,同理求出極值點(diǎn)序列在右端點(diǎn)處的近似取值。最后利用3次樣條函數(shù)對(duì)新極值點(diǎn)序列進(jìn)行插值得到上下包絡(luò)線。3次樣條函數(shù)在端點(diǎn)處有值可依,避免了上下包絡(luò)線的擺動(dòng)。

具體步驟如下：

（1）根據(jù)具體問題，確定擬合多項(xiàng)式的次數(shù)n；

（2）計(jì)算 Sr和 tq

式中m為終點(diǎn)的下標(biāo)；xi是離散數(shù)據(jù)點(diǎn)的橫坐標(biāo)；r=0，1，2，…，n，n＜m。

其中 yi為離散數(shù)據(jù)點(diǎn)的縱坐標(biāo)，q=0，1，…，n。

（3）由方程組(6)得到 a0，a1，…，an最后寫出擬合多項(xiàng)式[8]：

3.2 運(yùn)用多項(xiàng)式擬合算法對(duì)數(shù)據(jù)進(jìn)行處理

將式（3）的解析信號(hào)通過多項(xiàng)式擬合算法分解，可以看到在c3和c5分量（圖略）的端點(diǎn)處，發(fā)生了較大的區(qū)別，這是因?yàn)樵瓉聿捎玫亩它c(diǎn)處理手段是忽略內(nèi)部數(shù)據(jù)信息的方法，只是強(qiáng)制解決了端點(diǎn)飛翼，但同時(shí)也引入了較大的端點(diǎn)誤差，所以在端點(diǎn)處出現(xiàn)了較大偏差。進(jìn)一步對(duì)照原始信號(hào)圖（圖略）中的對(duì)應(yīng)c3頻段波形可以證明，經(jīng)過多項(xiàng)式擬合算法得到的結(jié)果，更加符合實(shí)際。

4 正交化經(jīng)驗(yàn)?zāi)B(tài)分解方法

Huang提出的EMD分解方法不能從理論上保證嚴(yán)格的正交性，僅僅能從數(shù)值上表明各IMF是近似正交的，不正交意味著解不唯一，也意味著信號(hào)冗余度的增加和不能嚴(yán)格遵守能量的守恒定律。

4.1 固有模態(tài)函數(shù)的正交性判斷依據(jù)

Huang定義了整體的正交性指標(biāo)和兩個(gè)分量

之間的正交性指標(biāo)，離散數(shù)據(jù)的整體正交性指標(biāo)[8]：

分量的正交性指標(biāo)：

此外還可以通過一個(gè)能量指標(biāo)來進(jìn)一步表明各分量之間的正交性[9]，設(shè)原始信號(hào)的能量為Ex：

各分量的能量為

如果分解出來的各個(gè)IMF分量是正交的，那么信號(hào)分解后的總能量應(yīng)保持不變且各分量之間的泄漏能量為零，即：

當(dāng)各IMF分量之間完全正交時(shí)，IOT和IOjk的值應(yīng)等于零。

4.2 IMF的正交化處理

對(duì)已經(jīng)得到的各階IMF分量進(jìn)行正交化,可以得到完全正交的各階IMF分量，步驟如下：

（2）正交化的IMF分量記做ci(t)，對(duì)第一個(gè)IMF分量，暫不作處理，即令c1(t)=c1(t)。

（3）從原來EMD分解得到的第二個(gè)IMF分量開始進(jìn)行正交化處理：

這樣就可以從第j+1階開始，消除各IMF中含有的非正交成分，可以得到信號(hào)X(t)的第j+1階真正的正交化 IMF分量 cj+1(t)(j=1，2，…，n)。

（4）對(duì) cj(t)做線性變換：

其中

當(dāng) i=j時(shí)，βi,j=1。

通過這樣的線性變換得到的各階c*j(t)是完全正交的，而且是完備的。

4.3 應(yīng)用正交化EMD方法對(duì)數(shù)據(jù)進(jìn)行處理

對(duì)式(3)中的信號(hào)使用正交化EMD方法進(jìn)行分解。把各IMF分量之間的正交性指標(biāo)列于表中（表略），由于正交性是具有對(duì)稱性的，表中數(shù)字是采用原始EMD分解得到的各IMF分量之間的正交性指標(biāo),表中下三角的數(shù)字是用新的方法分解的IMF分量之間的正交性指標(biāo)。

按照原始EMD分解得到的整體正交化指標(biāo)是0.064122，按正交化處理得到的整體正交化指標(biāo)是0.0006443，相差近100倍，可以看出原始EMD分解得到的各IMF之間正交化程度是比較差的。基于正交化分解得到的各階IMF分量之間的正交化程度要遠(yuǎn)好于原始EMD分解，其正交性指標(biāo)都在10-16以下。從能量分析，真實(shí)信號(hào)的總能量是Ex＝3392.2，原始EMD分解的誤差為8.89%,新的正交化EMD分解的誤差僅為2.48%(見表1)。

可以發(fā)現(xiàn)，新的正交化EMD分解的缺點(diǎn)是：使每個(gè)IMF產(chǎn)生毛刺，而毛刺是由于扣除非正交分量的高頻信號(hào)產(chǎn)生的結(jié)果，因?yàn)槊總€(gè)后面的IMF都要減去和前面IMF正交的部分，這就使得本來比較光滑的IMF變得不那么光滑。

5 結(jié)語

（1）Rilling算法可以控制信號(hào)在全局都不能出現(xiàn)超過標(biāo)準(zhǔn)的波動(dòng)(σ（t）≤θz)，而即使在某些局部出現(xiàn)較大的波動(dòng)，這些大的波動(dòng)也只能是全部時(shí)間的一小部分，優(yōu)于標(biāo)準(zhǔn)偏差法的全局籠統(tǒng)控制。

（2）端點(diǎn)處理采用多項(xiàng)式擬合方法，能夠根據(jù)數(shù)據(jù)內(nèi)部變化的走勢(shì)來判斷數(shù)據(jù)端點(diǎn)應(yīng)該采用的包絡(luò)線位置，可以避免延拓方法帶來的主觀因素和誤差。

表1 信號(hào)總能量、各分量能量及各分量能量之和

（3）正交化處理能夠保證各個(gè)IMF是嚴(yán)格正交的，使得各個(gè)IMF分量的能量和與原始波形的能量誤差很小，符合工程需要。

[1]Norden E.Huang,Zheng Shen,Steven R.Long,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proc.R.Soc.Lond.A,1998,454:903:995.

[2]HUANGNE，SHENZ，LONGSR.Anewviewofnonlinear waterwaves:the Hilbert spectrum[J].Ann Rev Fluid Mech,1999,31:417-457.

[3]沈國(guó)際,陶利民,陳仲生.多頻信號(hào)經(jīng)驗(yàn)?zāi)B(tài)分解的理論研究及應(yīng)用[J].振動(dòng)工程學(xué)報(bào),2005,18(1):91-94.

[4]丁康.平穩(wěn)和非平穩(wěn)振動(dòng)信號(hào)的若干處理方法及發(fā)展[J].振動(dòng)工程學(xué)報(bào),2003,16(1):1-10.

[5]RILLING G,FLANDRIN P,GONCALVES P.On Empirical ModeDecompositionanditsalgorithms[C].IEEEEURASIP Workshop on Nonlinear Signal and Image Processing,Grado(I),June9-11,2003.

[6]李彬彬,袁中凡,楊春生.改進(jìn)HHT算法及在心音信號(hào)分析中的應(yīng)用[J].四川大學(xué)學(xué)報(bào),2007,39(4):160-163.

[7]鄭天翔,楊力華.經(jīng)驗(yàn)?zāi)B(tài)分解算法的探討和改進(jìn)[J].中山大學(xué)學(xué)報(bào),2007,46(1):1-6.

[8]劉慧婷,倪志偉,李建洋.經(jīng)驗(yàn)?zāi)B(tài)分解方法及其實(shí)現(xiàn)[J].計(jì)算機(jī)工程與應(yīng)用2006(32):44-47.

[9]樓夢(mèng)麟,黃天立.正交化經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒╗J].同濟(jì)大學(xué)學(xué)報(bào),2007,35(3):293-298.

[10]Qiuhui Chen,N.E.Huang,et al.,A B-spline approach for empiricalmodedecompositions,AdvancesinConmutational Mathematics,2006,24:171-195.

[11]ZHAO Jin-ping,HUANG Da-ji.Mirror extending and circularspline function for empirical mode decomposition method[J].Journalof Zhejiang University(Science),2001,2(3):247-252.