項洪，吳琛, 杜喜朋

(1.福建工程學(xué)院 土木工程學(xué)院，福建 福州 350118; 2.福建省土木工程新技術(shù)與信息化重點實驗室，福建 福州 350118)

?

基于雙重對稱延拓的HHT端點效應(yīng)抑制方法

項洪1,2，吳琛1,2, 杜喜朋1,2

(1.福建工程學(xué)院 土木工程學(xué)院，福建 福州 350118; 2.福建省土木工程新技術(shù)與信息化重點實驗室，福建 福州 350118)

摘要：經(jīng)驗?zāi)B(tài)分解( empirical mode decomposition，EMD)存在端點飛翼的固有問題，使信號兩端出現(xiàn)扭曲失真。為抑制EMD端點效應(yīng)，文章提出雙重對稱延拓法，以端部數(shù)據(jù)對稱延拓作為首次延拓，以極值點對稱延拓作為二次延拓，該方法可同時實現(xiàn)EMD分解和Hilbert變換兩階段端點效應(yīng)的抑制。仿真信號和地震響應(yīng)的工程算例通過IMF分量對比、瞬時頻率對比，以及整體正交性、相關(guān)系數(shù)等指標(biāo)的對比，驗證了該方法具有良好的端點效應(yīng)抑制效果,同時還能改善IMF分量之間的正交性，并節(jié)約計算時間。

關(guān)鍵詞：Hilbert-Huang變換; 經(jīng)驗?zāi)B(tài)分解; 端點效應(yīng); 雙重對稱延拓; 評價指標(biāo)

Hilbert-Huang變換[1](Hilbert-Huang transform , HHT) 包括經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)和Hilbert譜分析兩個階段。HHT具有完全自適應(yīng)性的優(yōu)點，處理非線性、非平穩(wěn)的信號具有明顯優(yōu)勢；自1998年Norden E.Huang等人提出后，被廣泛應(yīng)用于機械故障診斷[2-3]、地震工程[4-5]等領(lǐng)域。但在實際工程應(yīng)用中，HHT還存在一些問題：如端點問題、包絡(luò)線擬合問題、模態(tài)混疊等，這些問題導(dǎo)致HHT分析得到的物理信息不夠理想。本文主要針對HHT過程中EMD分解和Hilbert變換所產(chǎn)生的端點效應(yīng)問題進行研究。

1　端點效應(yīng)分析

HHT在處理非平穩(wěn)信號時存在以下兩類端點效應(yīng)問題：

第一類是EMD分解過程中出現(xiàn)的端點飛翼現(xiàn)象。在每一次的EMD篩分過程中，應(yīng)根據(jù)信號的上、下極值點三次樣條插值包絡(luò)線來計算信號的局部平均值。由于信號不是無限長的，所以不能同時確定信號左右兩端的局部極值點。如果直接在端點處棄值，將導(dǎo)致三次樣條插值得到的上下包絡(luò)線在信號兩端附近扭曲失真(如圖1)，產(chǎn)生端點飛翼。

圖1　EMD分解端點飛翼示意圖　　Fig.1　The end swing result of EMDdecomposition

針對EMD分解的端點飛翼現(xiàn)象，Huang等人最早提出利用原始信號進行延拓的思路。此后，許多文獻提出用不同延拓的方法抑制端點飛翼：如極值點對稱延拓[6]、多項式擬合延拓[7]、極值平移法[8]等均采用向原始數(shù)據(jù)添加極值點的思路；又如神經(jīng)網(wǎng)絡(luò)[9]、支持向量回歸機[10-11]、最大相關(guān)波形[12]等方法則采用延長原信號數(shù)據(jù)的思路。

第二類是Hilbert變換過程中出現(xiàn)的端點發(fā)散現(xiàn)象[9]。對各個IMF分量進行Hilbert變換時，由于Hilbert變換是基于傅立葉變換實現(xiàn)的，而傅立葉變換會產(chǎn)生頻譜泄露，出現(xiàn)嚴(yán)重的端點發(fā)散現(xiàn)象(如圖2)，所得Hilbert譜無法真實地反映原始信號的特征。

圖2　Hilbert變換端點發(fā)散示意圖Fig.2　The end effects result of Hilbert transform

針對Hilbert變換過程中產(chǎn)生的端點發(fā)散現(xiàn)象，目前主要通過神經(jīng)網(wǎng)絡(luò)[9]、時變參數(shù)ARMA[13]和支持向量回歸機[10]等數(shù)據(jù)延拓法先對IMF分量進行端點延拓，再進行Hilbert變換，并通過拋棄延拓的兩端數(shù)據(jù)，使端點效應(yīng)釋放到原始數(shù)據(jù)的外端，保證原始數(shù)據(jù)范圍內(nèi)Hilbert變換的有效性。

上述HHT端點處理方法對特定信號的端點抑制效果較理想，但對于非線性、非平穩(wěn)信號，由于其復(fù)雜與多變性，上述方法在計算精度和計算速度上難以同時達(dá)到最佳。本文將延長原信號數(shù)據(jù)與添加極值點相結(jié)合，提出在EMD分解過程中實施端部信號與極值點雙重對稱延拓的方法以抑制HHT的端點效應(yīng)。與現(xiàn)有的極值點對稱延拓等方法相比，這種新方法的端點效應(yīng)抑制效果更佳，并能改善信號的整體正交性。與神經(jīng)網(wǎng)絡(luò)、時變參數(shù)ARMA和支持向量回歸機的數(shù)據(jù)延拓等方法相比，本文算法更為簡單，計算速度更快。并且，這種新方法不需要在Hilbert變換之前再對IMF分量進行端點延拓，可以在一次延拓中同時抑制EMD分解和Hilbert變換產(chǎn)生的端點發(fā)散。

2　雙重對稱延拓的端點效應(yīng)處理

雙重對稱延拓首先進行端部數(shù)據(jù)對稱延拓，按照原信號在端點附近的自然趨勢延長原信號數(shù)據(jù)，把端點效應(yīng)釋放至原始信號外側(cè)，從而使端點效應(yīng)更緩慢影響到原信號內(nèi)部。但是，如圖3所示，如果僅進行端部數(shù)據(jù)延拓，三次樣條插值所形成的包絡(luò)曲線在端部依然存在端點發(fā)散，經(jīng)過EMD分解，將影響到原始信號內(nèi)包絡(luò)曲線的準(zhǔn)確性。因此，本文提出在端部數(shù)據(jù)延拓的基礎(chǔ)上進一步進行極值點對稱延拓，進一步控制包絡(luò)曲線的走勢，確保其準(zhǔn)確性。

圖3　端部數(shù)據(jù)對稱延拓端點飛翼示意圖Fig.3　The end flying result of end signal data symmetrical extension

2.1端部數(shù)據(jù)對稱延拓

設(shè)原始信號為x(t)，對x(t)執(zhí)行如下操作：

1)求x(t)的極大值indmax(1)，indmax(2)，…，indmax(n)和x(t)的極小值indmin(1)，indmin (2)，…，indmin (m)。其中，m=n，或丨m-n丨=1。

2)左端信號數(shù)據(jù)對稱延拓：

當(dāng)與左端點x(0)相鄰的極值點為極大值時，比較x(0)與indmin(l)：

a.若x(0)>indmin(1)，以indmax(l)為對稱點，將indmax(l)～indmax(3)之間的信號數(shù)據(jù)向左對稱延拓(如圖4(a)所示)。

b.若x(0)≤indmin(1)，以x(0)為對稱點，將x(0)～indmin(2)之間的信號數(shù)據(jù)向左對稱延拓(如圖4(b)所示)。

當(dāng)與左端點x(0)相鄰的極值點為極小值時，比較x(0)與indmax(l)：

c.若x(0)

d.若x(0)≥indmax(1)，以x(0)為對稱點，將x(0)～indmax(2)之間的信號數(shù)據(jù)向左對稱延拓(如圖4(d)所示)。

圖4　左端數(shù)據(jù)對稱延拓Fig.4　Left end data symmetrical extension

3)同理，對右端信號對稱延拓，形成信號記為x1(t)。

2.2極值點對稱延拓[6]

對經(jīng)過端部數(shù)據(jù)對稱延拓的信號x1(t)執(zhí)行極值點對稱延拓。比較信號的端點值與端部極值，當(dāng)端點值大于相鄰的第一個極大值或者是小于相鄰的第一個極小值時，以端點作為對稱點，當(dāng)端點值在相鄰極值中間時，以相鄰極值點作為對稱點，向信號外側(cè)鏡像添加兩個極值點。

2.3按雙重延拓后信號的包絡(luò)線執(zhí)行EMD篩分

1)確定信號x1(t)延拓的所有局部極值點，用三次樣條曲線連接信號x1(t)所有的局部極大值和局部極小值，分別形成上包絡(luò)線x1max(t)和下包絡(luò)線x1min(t)。

2)求出x1max(t)和x1min(t)的平均值，記為m11(t)，將原數(shù)據(jù)序列x1(t)去掉該平均值后得到新數(shù)據(jù)序列h11(t)：

3) 判斷h11(t)是否滿足極值點個數(shù)與過零點個數(shù)相等或至多相差一個、局部極大值和局部極小值定義的包絡(luò)均值為零這兩個IMF的基本要求[1]，如果不滿足則將h11(t)作為原始數(shù)據(jù)，重復(fù)上述1)～3)處理過程，直到新數(shù)據(jù)序列

滿足IMF的條件，得到第一個IMF分量c1(t)：

4) 從x1(t)中分離出c1(t)，得到剩余序列r1(t)：

5) 將r1(t)作為一個新的原序列，按照以上1)～4)步驟，依次提取第2、3直至第n個固有模態(tài)函數(shù)cn(t)。當(dāng)殘量rn(t)成為一個單調(diào)函數(shù)或小于某一預(yù)定值時，分解結(jié)束。截取原始信號x(t)時域內(nèi)IMF分量。

2.4EMD分解得到IMF分量直接進行Hilbert變換

信號x1(t)經(jīng)EMD分解得到IMF分量ci(t)，不截取原始信號時域內(nèi)的數(shù)據(jù)，保留延拓部分，直接進行Hilbert變換，將Hilbert變換的端點效應(yīng)釋放到原始數(shù)據(jù)外端：

構(gòu)造解析信號

進一步取得瞬時頻率

最后截取原始信號x(t)時域內(nèi)瞬時頻率。

3　仿真信號HHT分析

3.1EMD分解端點效應(yīng)處理

已知簡諧疊加信號x(t)=cos(2πt/50)+ 0.6cos(2πt/25) + 0.5cos(2πt/200)，t∈[1,300]。如圖5所示?，F(xiàn)按雙重對稱延拓方法進行HHT分析。

圖5　簡諧疊加信號Fig.5　Simple harmonic superimposed signal

圖6是未經(jīng)端點處理，直接以端點值作為極值點得到的IMF分量，其中虛線表示真實分量。該圖顯示，IMF1的端部信號明顯收縮；隨著頻率減小，低頻分量由于其時間尺度大，極值點之間的距離大，分量信號的內(nèi)部更容易受到端部邊緣效應(yīng)的影響，導(dǎo)致EMD分解得到的低頻分量失真。

圖6　采用端點值作為極值點得到的IMF分量Fig.6　The IMF components using endpoint value as extreme point

圖7為僅對該信號進行極值點對稱延拓所獲得的IMF分量。雖然端點效應(yīng)在高頻成份中得到了抑制，但在低頻分量中，IMF2和IMF3右端點與真實值的誤差達(dá)到30%和51%。

圖7　采用極值點對稱延拓得到的IMF分量Fig.7　The IMF components with extreme symmetrical extension

圖8為僅對該信號進行端部數(shù)據(jù)對稱延拓獲得的IMF分量，該方法由于利用了更多的信號數(shù)據(jù)參與EMD分解，因此抑制效果略優(yōu)于極值延拓的情況。但隨著頻率減小，低頻分量的端點效應(yīng)逐漸傳播到信號內(nèi)部，IMF3分量端點與真實值的最大誤差仍達(dá)到30%。

圖8　采用端部數(shù)據(jù)對稱延拓得到的IMF分量Fig.8　The IMF components with signal end data symmetrical extension

圖9為利用本文提出的雙重對稱延拓方法得到的IMF分量，與圖8相比，該方法不僅在高頻分量與真實分量吻合較好，低頻分量也幾乎接近真實信號。說明其更加有效地抑制了EMD分解的端點效應(yīng)。

圖9　采用雙重對稱延拓方法得到的IMF分量Fig.9　The IMF components with dual symmetrical extension

3.2EMD端點效應(yīng)抑制效果評價指標(biāo)

采用運算時間t、整體正交性指數(shù)IO[14]、信號分量相關(guān)系數(shù)ρ[15]等3個指標(biāo)來進一步綜合評價端點值作為極值點、極值點對稱延拓、端點數(shù)據(jù)對稱延拓、雙重對稱延拓等端點效應(yīng)處理方法的抑制效果。

從表1可得，本文提出的雙重對稱延拓處理法的正交性指數(shù)IO值最小，端點值作為極值點的方法得到IO值最大，說明本文提出的雙重延拓法比僅經(jīng)端部數(shù)據(jù)延拓或極值點延拓的方法更有助于改善IMF分量的正交性。通過運算時間的比較可知雖然本文方法用時最長，但運行時間僅0.312 s，具有很高的計算效率。并且，本文方法的各個分量相關(guān)系數(shù)ρ最高，說明EMD分解結(jié)果與真實分量的誤差最小。因此，經(jīng)上述綜合評價，雙重對稱延拓法是一種端點效應(yīng)抑制良好、算法簡單、計算高效、IMF分量之間正交性程度較高的改進方法。

表1不同延拓方法得到的端點效應(yīng)抑制效果評價指標(biāo)

Tab.1The evaluation indexes of end effects reduction by different extension methods

評價指標(biāo)端點值作為極值點極值點對稱延拓端部數(shù)據(jù)對稱延拓雙重對稱延拓IO值0.05760.03780.03360.0199用時/s0.28000.29700.29200.3120ρ10.93950.99510.99670.9982ρ20.95760.99160.99810.9993ρ30.88370.98100.93890.9974

3.3Hilbert變換的端點效應(yīng)處理

將上述仿真信號x(t)采用端點值作為極值點方法得到瞬時頻率如圖10所示，端點發(fā)散效應(yīng)清晰可見。

圖10　采用端點值作為極值點得到的瞬時頻率Fig.10　The instantaneous frequency with endpoint value as extreme point

利用極值點對稱延拓后得到的瞬時頻率見圖11，邊緣效應(yīng)雖有所改善，但仍然沒有完全得以抑制。

圖11　采用極值點對稱延拓得到的瞬時頻率Fig.11　The instantaneous frequency with extreme symmetrical extension

圖12是采用端部數(shù)據(jù)對稱延拓方法得到的瞬時頻率，其端點問題釋放到原始數(shù)據(jù)的外端，其端點發(fā)散有顯著的抑制，僅在瞬時頻率為0.005 Hz的曲線右端出現(xiàn)發(fā)散，發(fā)散部分與真實值的最大誤差達(dá)到45%。可見，數(shù)據(jù)延拓的思路對于抑制EMD分解和Hilbert變換兩階段的端點效應(yīng)均有效，但在低頻分量中，仍可能因延拓部分包絡(luò)失真影響到原始信號端部的準(zhǔn)確性。

圖12　采用端部數(shù)據(jù)對稱延拓得到的瞬時頻率Fig.12　The instantaneous frequency with signal ends data symmetrical extension

圖13是采用本文提出的雙重對稱延拓方法得到的瞬時頻率圖，從圖中可清楚看出3個分量頻率分別對應(yīng)仿真信號的3個頻率成份0.04、0.02、0.005 Hz。相比，瞬時頻率為0.005 Hz的曲線右端更加平穩(wěn)，端點效應(yīng)抑制得更加徹底。

圖13　采用雙重對稱延拓得到的瞬時頻率Fig.13　The instantaneous frequency with dual symmetrical extension

4　工程應(yīng)用

現(xiàn)通過結(jié)構(gòu)地震響應(yīng)的分析進一步驗證本文方法的有效性。已知單自由結(jié)構(gòu)體系質(zhì)量為107kg，剛度為5×107N/m，結(jié)構(gòu)自振頻率為0.355 9 Hz，結(jié)構(gòu)阻尼比為0.05。輸入荷載為cos 2πt，采樣時間間隔為0.01 s，采樣的個數(shù)為1 000。經(jīng)彈性時程分析，該結(jié)構(gòu)絕對加速度響應(yīng)如圖14所示。

圖14　結(jié)構(gòu)絕對加速度響應(yīng)Fig.14　Structural absolute acceleration response

信號經(jīng)過雙重對稱延拓處理后的EMD分解結(jié)果如圖15所示。由圖可知絕對加速度響應(yīng)由高頻至低頻依次分解得到3階IMF分量和一個剩余量，說明EMD分解中分辨率是自適應(yīng)的。各個分量兩端未出現(xiàn)大幅度的擺動現(xiàn)象，這說明提出的雙重延拓方法能較好的抑制EMD端點飛翼。

圖15　雙重對稱延拓處理的EMD結(jié)果Fig.15　EMD result with dual symmetrical extension

將保留延拓數(shù)據(jù)的IMF分量進行Hilbert變換得到的瞬時頻率如圖16所示，從圖中可以清晰看出前三階分量的頻率分別約為1、0.36、0.2 Hz，同時頻率信號端部也無明顯的端點發(fā)散現(xiàn)象，說明雙重對稱延拓對于實際應(yīng)用來說也是有效的。EMD分解出來的IMF分量可以得到原始信號無法得到一些物理信息。圖15中的前兩個分量描述了體系動力響應(yīng)的本征振動模態(tài)：第一個分量描述的是體系動力響應(yīng)中與輸入荷載有關(guān)的穩(wěn)態(tài)振動部分；第二個分量描述的是體系動力響應(yīng)中與結(jié)構(gòu)自振特性有關(guān)的伴生自由振動部分。

圖16　前三階IMF分量對應(yīng)的瞬時頻率Fig.16　The instantaneous frequency of the first three order IMFs

表2給出了4種分解方法對于結(jié)構(gòu)絕對加速度響應(yīng)處理后的兩個性能評價結(jié)果，其中雙重對稱延拓的IO值最小，EMD分解的各個分量正交性最好。同時運行的時間只有0.3 s，運行速度快。

表2絕對加速度響應(yīng)信號不同延拓對比

Tab.2The comparison of extension among absolute acceleration response(signals)

評價指標(biāo)端點值作為極值點極值點對稱延拓端部數(shù)據(jù)對稱延拓雙重對稱延拓IO值0.06810.06330.06520.0520用時/s0.22000.28000.26000.3000

5　結(jié)論

以雙重對稱延拓的方法較端點值作為極值點、極值點對稱延拓、端部數(shù)據(jù)對稱延拓等方法能更好地抑制EMD分解帶來的端點飛翼問題，同時改善了信號的整體正交性，使Hilbert變換的端點效應(yīng)釋放到原始信號的外端，從而抑制了Hilbert變換的端點發(fā)散效應(yīng)。與神經(jīng)網(wǎng)絡(luò)、時變參數(shù)ARMA和支持向量回歸機的數(shù)據(jù)延拓等方法相比，本文提出的方法更為簡單，計算速度更快，可以在一次延拓中同時抑制EMD分解和Hilbert變換產(chǎn)生的端點發(fā)散。同時該方法在實際信號中也能有效地抑制端點效應(yīng)問題。

參考文獻：

[1] Huang N E, Shen Z, Long S R.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of Royal Society of London,1998,454: 903-995.

[2] 鄭近德，程軍圣.改進的希爾伯特-黃變換及其在滾動軸承故障診斷中的應(yīng)用[J].機械工程學(xué)報,2015(1):138-145.

[3] 熊炘.基于希爾伯特-黃變換的故障轉(zhuǎn)子振動模式分析方法研究[D].杭州：浙江大學(xué),2012.

[4] 吳祚菊，張建經(jīng)，王志佳，等.改進HHT算法在地震動特性分析中的應(yīng)用研究[J].地震工程與工程振動，2015(2):86-93.

[5] 吳琛，周瑞忠.Hilbert-Huang變換在提取地震信號動力特性中的應(yīng)用[J].地震工程與工程振動,2006(5):41-46.

[6] Rilling G, Flandrin P, Goncalves P.On empirical mode decomposition and its algorithm[C]//Proceeding of IEEE-EURASIP Workshop on Nonlinear Signal and Image (NSIPL2003).Trieste: Institute of Electrical and Electronics Engineers,2003:8-11.

[7] 劉慧婷，張旻，程家興.基于多項式擬合算法的EMD端點問題的處理[J].計算機工程與應(yīng)用,2004,40(16):84-86,100.

[8] 王學(xué)敏，黃方林.EMD端點效應(yīng)抑制的一種實用方法[J].振動、測試與診斷，2012(3):493-497.

[9] 鄧擁軍，王偉，錢成春，等.EMD方法及Hilbert變換中邊界問題的處理[J].科學(xué)通報，2001(3):257-263.

[10] 程軍圣，于德介，楊宇.基于支持矢量回歸機的Hilbert-Huang變換端點效應(yīng)問題的處理方法[J].機械工程學(xué)報，2006(4):23-31.

[11] 曠歡，王如龍，張錦，等.基于SVM的EMD端點效應(yīng)抑制方法研究[J].計算機工程與應(yīng)用，2015(11):196-200.

[12] 高強，段晨東，趙艷青，等.基于最大相關(guān)波形延拓的經(jīng)驗?zāi)Ｊ椒纸舛它c效應(yīng)抑制方法[J].振動與沖擊，2013(2):62-66.

[13] 程軍圣，于德介，楊宇.Hilbert-Huang變換端點效應(yīng)問題的探討[J].振動與沖擊，2005(6):40-42.

[14] 肖韻.Hilbert-Huang變換端點問題處理方法的分析與研究[D].成都：西南交通大學(xué),2012.

[15] 王秋生，段丹輝.經(jīng)驗?zāi)B(tài)分解的邊界效應(yīng)處理技術(shù)[J].計算機測量與控制，2006(12):1673-1675.

(責(zé)任編輯： 陳雯)

A dual symmetrical extension method for end effects reduction of Hilbert-Huang transform

Xiang Hong1,2, Wu Chen1,2, Du Xipeng1,2

(1.College of Civil Engineering, Fujian University of Technology, Fuzhou 350118, China;2.Fujian Provincial Key Laboratory of Advanced Technology and Informationization in Civil Engineering, Fuzhou 350118, China)

Abstract：End effects reduce the precision of empirical mode decomposition (EMD) greatly．A dual symmetrical extension method was proposed to decrease the end effects of EMD.Signal end data symmetrical extension was regarded as the first extension, whereas extreme symmetrical extension was regarded as the second extension.The proposed method can also decrease the end effects of Hilbert transform.The harmonic superimposed signal and absolute acceleration response were applied to test the performance of the proposed method.The results show that the proposed method can reduce the end effects of HHT and can improve its precision significantly.

Key words：Hilbert-Huang transform; empirical mode decomposition; end effect; dual symmetrical extension method; evaluation index

doi:10.3969/j.issn.1672-4348.2016.04.004

收稿日期:2016-05-31

基金項目：國家自然科學(xué)基金項目(51108083); 福建省自然科學(xué)基金項目(2015J01631)

第一作者簡介：項洪(1991-)，男，福建龍巖人，碩士研究生，主要從事結(jié)構(gòu)抗震理論與工程應(yīng)用研究.

中圖分類號：TN911.72

文獻標(biāo)志碼：A

文章編號：1672-4348(2016)04-0324-08