李宏偉，劉宇航，楊 輝

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

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

1 EMD分解

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

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

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

2 Rilling算法

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

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

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

并引進新的評估函數(shù)：

篩選停止條件有2個：

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

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

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

為考查上述方法的效果，采用一個簡單的解析信號：

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

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

3 多項式擬合求端點的極值

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

3.1 多項式擬合算法

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

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

具體步驟如下：

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

（2）計算 Sr和 tq

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

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

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

3.2 運用多項式擬合算法對數(shù)據(jù)進行處理

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

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

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

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

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

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

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

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

各分量的能量為

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

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

4.2 IMF的正交化處理

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

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

（3）從原來EMD分解得到的第二個IMF分量開始進行正交化處理：

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

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

其中

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

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

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

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

按照原始EMD分解得到的整體正交化指標(biāo)是0.064122，按正交化處理得到的整體正交化指標(biāo)是0.0006443，相差近100倍，可以看出原始EMD分解得到的各IMF之間正交化程度是比較差的?；谡换纸獾玫降母麟AIMF分量之間的正交化程度要遠好于原始EMD分解，其正交性指標(biāo)都在10-16以下。從能量分析，真實信號的總能量是Ex＝3392.2，原始EMD分解的誤差為8.89%,新的正交化EMD分解的誤差僅為2.48%(見表1)。

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

5 結(jié)語

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

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

表1 信號總能量、各分量能量及各分量能量之和

（3）正交化處理能夠保證各個IMF是嚴格正交的，使得各個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]沈國際,陶利民,陳仲生.多頻信號經(jīng)驗?zāi)B(tài)分解的理論研究及應(yīng)用[J].振動工程學(xué)報,2005,18(1):91-94.

[4]丁康.平穩(wěn)和非平穩(wěn)振動信號的若干處理方法及發(fā)展[J].振動工程學(xué)報,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]李彬彬,袁中凡,楊春生.改進HHT算法及在心音信號分析中的應(yīng)用[J].四川大學(xué)學(xué)報,2007,39(4):160-163.

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

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

[9]樓夢麟,黃天立.正交化經(jīng)驗?zāi)Ｊ椒纸夥椒╗J].同濟大學(xué)學(xué)報,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.