經(jīng)驗?zāi)Ｊ椒纸飧倪M(jìn)算法的比較

李宏偉，劉宇航，楊 輝

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

隨著HHT研究的深入和應(yīng)用的拓展，還有一些關(guān)鍵的理論和實際問題需要完善和改進(jì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分解中的濾波停止條件進(jìn)行了改進(jìn)，該方法的核心思想是在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)幅值：

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

篩選停止條件有2個：

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

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

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

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

分別用標(biāo)準(zhǔn)偏差法和Rilling濾波方法進(jìn)行計算比較。取 θ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ù)中的上極值點和下極值點分別進(jìn)行三次樣條插值擬合再求平均得到包絡(luò)均值，在樣條插值時常常不能確定端點處的極值，就會在樣條插值時產(chǎn)生數(shù)據(jù)的擬合誤差。不斷累積的誤差就會由端點處向內(nèi)逐漸傳播，信號的兩端會出現(xiàn)嚴(yán)重的扭曲。

3.1 多項式擬合算法

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

取出原極值點序列最左端的3個極值點(如果極值點序列的個數(shù)小于3個則取序列中所有元素),對所取的極值點求出擬合多項式,計算出多項式對應(yīng)數(shù)據(jù)序列左端點處的函數(shù)值,并把此函數(shù)值作為極值點序列在該端點處的近似取值,同理求出極值點序列在右端點處的近似取值。最后利用3次樣條函數(shù)對新極值點序列進(jìn)行插值得到上下包絡(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ù)進(jìn)行處理

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

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

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

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

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

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

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

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

各分量的能量為

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

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

4.2 IMF的正交化處理

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

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

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

這樣就可以從第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ù)進(jìn)行處理

對式(3)中的信號使用正交化EMD方法進(jìn)行分解。把各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分量之間的正交化程度要遠(yuǎn)好于原始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是嚴(yán)格正交的，使得各個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]李彬彬,袁中凡,楊春生.改進(jìn)HHT算法及在心音信號分析中的應(yīng)用[J].四川大學(xué)學(xué)報,2007,39(4):160-163.

[7]鄭天翔,楊力華.經(jīng)驗?zāi)B(tài)分解算法的探討和改進(jìn)[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.