徐力彬，宋余慶，劉 毅

(江蘇大學(xué)計(jì)算機(jī)科學(xué)與通信工程學(xué)院，江蘇鎮(zhèn)江212013)

1 概述

Hilbert-Huang變換(Hilbert-Huang Transform，HHT)是由Huang N E等人通過(guò)深入分析和總結(jié)后提出的一種新的非平穩(wěn)非線性信號(hào)分析方法［1］。該方法主要包含2個(gè)部分:經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)和 Hilbert譜 分 析(Hilbert Spectral Analysis，HAS)，其中 EMD 為核心部分［2］。與傳統(tǒng)方法相比，EMD無(wú)須預(yù)先設(shè)定任何基函數(shù)，而是根據(jù)信號(hào)特點(diǎn)自適應(yīng)地將原始信號(hào)分解為一系列固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，去除了快速傅里葉變換方法中引入的諧波基函數(shù)，也沒(méi)有小波分析復(fù)雜信號(hào)過(guò)程中產(chǎn)生信號(hào)模糊與失真的問(wèn)題［3-6］。然而，HHT變換也有一定的局限性，例如端點(diǎn)效應(yīng)問(wèn)題［7］和模態(tài)混疊問(wèn)題［8］等。尤其是端點(diǎn)效應(yīng)問(wèn)題會(huì)嚴(yán)重妨礙信號(hào)的特征提取，從而導(dǎo)致信號(hào)分析的誤差。

EMD中使用三次樣條函數(shù)對(duì)包絡(luò)進(jìn)行擬合時(shí)，由于難以確定端點(diǎn)處是否為極值點(diǎn)，導(dǎo)致出現(xiàn)了端點(diǎn)效應(yīng)問(wèn)題。國(guó)內(nèi)外專家針對(duì)端點(diǎn)效應(yīng)已經(jīng)提出了多種有效的抑制方法，包括鏡像延拓［9］、極值延拓［10］、基于神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的數(shù)據(jù)延拓法［11］以及多項(xiàng)式擬合延拓法［12］等。這些方法在抑制端點(diǎn)效應(yīng)方面都有著良好的效果，但同時(shí)也存在弊端，如多種延拓方法存在同一種缺陷，即在對(duì)原始信號(hào)進(jìn)行延拓處理以后，其端點(diǎn)依然是不固定的，擬合的上下包絡(luò)線兩端還是會(huì)出現(xiàn)發(fā)散現(xiàn)象，只是在EMD完成后截去延拓部分，端點(diǎn)效應(yīng)相對(duì)減輕。本文在鏡像延拓的基礎(chǔ)上，利用余弦窗函數(shù)法對(duì)延拓后的信號(hào)進(jìn)行處理，使樣條函數(shù)可以更好地?cái)M合包絡(luò)線，對(duì)端點(diǎn)效應(yīng)起到抑制作用。

2 EMD算法

2．1 算法簡(jiǎn)介

EMD算法是HHT的核心算法之一，用來(lái)將一個(gè)信號(hào)的能量按照時(shí)域各種固有尺度的波動(dòng)進(jìn)行分解，得到一系列頻率從小到大的IMF。這里IMF須滿足如下2個(gè)性質(zhì)［13］:

(1)信號(hào)的極值點(diǎn)數(shù)目和過(guò)零點(diǎn)數(shù)目相等或最多相差一個(gè);

(2)由局部極大值構(gòu)成的上包絡(luò)線和由局部極小值構(gòu)成的下包絡(luò)線的平均值為0。

EMD的分解過(guò)程如下:

Step1設(shè)原始信號(hào)曲線為X(t)，求出所有的極大值、極小值。

Step2對(duì)極值使用三次樣條函數(shù)進(jìn)行處理，描繪出原始信號(hào)的上、下包絡(luò)線，分別記為 u(t)和v(t)。

Step3將上、下包絡(luò)線的平均值記為:

Step4將均值曲線m(t)從原信號(hào)X(t)中分離出來(lái)，得到的剩余部分記為:

Step5判斷h1(t)是否滿足IMF的2個(gè)性質(zhì)。若滿足，則h1(t)為第一個(gè)IMF分量;否則，記h1(t)為X(t)，重復(fù)Step1～Step4，直到第k次迭代之后得到一個(gè)IMF，將其記為c1(t)=h1k(t)。這里，將標(biāo)準(zhǔn)差準(zhǔn)則作為第k次迭代的終止原則，使其位于0．2～0．3 之間，停止篩選［14］。

Step6記r1(t)=X(t)-c1(t)為新的待處理信號(hào)，再重復(fù)Step1～Step5，得到第2個(gè)IMF分量，將其記為c2(t)。繼續(xù)重復(fù)Step1～Step6，當(dāng)所得余項(xiàng)rn(t)是一個(gè)單調(diào)信號(hào)或者rn(t)的值小于事先給定的閾值時(shí)，分解結(jié)束［15］。

綜上所述，原始信號(hào)X(t)為所有IMF分量與余項(xiàng)rn(t)的和:

2．2 端點(diǎn)效應(yīng)

在EMD中，為得到信號(hào)的瞬時(shí)平均包絡(luò)線，需要對(duì)原信號(hào)的局部極大值和極小值分別進(jìn)行三次樣條插值算法擬合出上、下包絡(luò)線后，再計(jì)算出信號(hào)的局部平均值。由于信號(hào)兩端無(wú)法同時(shí)處于極大值和極小值點(diǎn)，因此上、下包絡(luò)線在數(shù)據(jù)信號(hào)兩端不可避免地會(huì)出現(xiàn)偏移，并且這種偏移現(xiàn)象會(huì)隨EMD的進(jìn)行逐漸傳播到信號(hào)內(nèi)部，從而嚴(yán)重影響信號(hào)分解的質(zhì)量［16］，這種現(xiàn)象就是端點(diǎn)效應(yīng)，如圖1所示，其中，實(shí)心點(diǎn)為延拓后的點(diǎn)，實(shí)線為真實(shí)的包絡(luò)線，虛線為無(wú)端點(diǎn)約束產(chǎn)生的錯(cuò)誤包絡(luò)線。

圖1 端點(diǎn)效應(yīng)示意圖

2．3 抑制端點(diǎn)效應(yīng)的常用方法

2．3．1 常用延拓方法

延長(zhǎng)數(shù)據(jù)序列或者在數(shù)據(jù)兩端增加極值點(diǎn)來(lái)抑制端點(diǎn)效應(yīng)是得到普遍認(rèn)可的一種方法［17］。對(duì)于復(fù)雜的信號(hào)，準(zhǔn)確地對(duì)數(shù)據(jù)序列進(jìn)行延拓是不現(xiàn)實(shí)的，但可以使延拓后得到的平均包絡(luò)線與真實(shí)平均包絡(luò)線相對(duì)較為接近，這是到目前為止普遍認(rèn)為的抑制端點(diǎn)效應(yīng)的有效途徑。常用的延拓方法如下:

(1)鏡像延拓法

鏡像延拓法的實(shí)質(zhì)是將原始信號(hào)對(duì)稱地延拓成一個(gè)閉合的環(huán)形信號(hào)，這種環(huán)形信號(hào)沒(méi)有端點(diǎn)，其上、下包絡(luò)線完全由信號(hào)內(nèi)部確定，從而在根本上避免了端點(diǎn)效應(yīng)的產(chǎn)生。但是，當(dāng)無(wú)法確定端點(diǎn)是否為極值點(diǎn)時(shí)，可能會(huì)使延拓部分和原始信號(hào)的均值有明顯差異，從而影響EMD產(chǎn)生的最終結(jié)果。

(2)極值延拓法

極值延拓法不需要對(duì)信號(hào)本身延拓，可以只對(duì)信號(hào)的極值進(jìn)行延拓。該方法在鏡像延拓的基礎(chǔ)上，對(duì)延拓后的對(duì)稱點(diǎn)是否為極值點(diǎn)進(jìn)行了判斷，提高了處理效果［18］。但該方法也存在缺陷，即若將端點(diǎn)都作為極值點(diǎn)進(jìn)行處理，端點(diǎn)處的包絡(luò)線會(huì)產(chǎn)生收縮，從而導(dǎo)致包絡(luò)線失真。

(3)多項(xiàng)式擬合延拓法

通過(guò)多項(xiàng)式擬合法算出信號(hào)端點(diǎn)處對(duì)應(yīng)的函數(shù)值，將其作為端點(diǎn)處極值的近似值。用正交多項(xiàng)式擬合的方法會(huì)比一般的多項(xiàng)式擬合更加準(zhǔn)確，算法復(fù)雜度也不搞。該方法對(duì)準(zhǔn)周期信號(hào)延拓效果比較好，但對(duì)隨機(jī)信號(hào)等變化規(guī)律不明顯的信號(hào)在端點(diǎn)效應(yīng)抑制方面的效果一般。

2．3．2 窗函數(shù)法

窗函數(shù)法是一種能夠有效抑制端點(diǎn)效應(yīng)的方法，對(duì)信號(hào)進(jìn)行窗函數(shù)處理可以增強(qiáng)函數(shù)中心處的信號(hào)，同時(shí)抑制遠(yuǎn)離中心處的信號(hào)［19］。在EDM 分解過(guò)程中，端點(diǎn)效應(yīng)從端點(diǎn)逐漸向信號(hào)內(nèi)部發(fā)散，誤差不斷積累，導(dǎo)致發(fā)散程度越來(lái)越嚴(yán)重。

在對(duì)信號(hào)進(jìn)行窗函數(shù)處理后，信號(hào)端點(diǎn)處的值變?yōu)?，包絡(luò)線也會(huì)相對(duì)平滑。因此，三次樣條函數(shù)能夠更真實(shí)地?cái)M合包絡(luò)線，最終上、下包絡(luò)線收逐漸收斂于端點(diǎn)，不再發(fā)散。

在窗函數(shù)的選取問(wèn)題上，筆者不需要自己設(shè)計(jì)窗函數(shù)，國(guó)內(nèi)外數(shù)學(xué)家已經(jīng)提出了很多能適合不同需要的窗函數(shù)。在處理端點(diǎn)效應(yīng)的問(wèn)題上，應(yīng)用較多的有漢寧窗、海明窗、矩形窗和余弦窗等。

3 本文算法

當(dāng)原信號(hào)進(jìn)行窗函數(shù)處理后，原信號(hào)的包絡(luò)線將被改變，可能會(huì)對(duì)接下來(lái)的過(guò)程帶來(lái)誤差，這時(shí)可以考慮先將信號(hào)進(jìn)行延拓處理，再進(jìn)行加窗，最后在EMD過(guò)程結(jié)束后舍去延拓部分。

仿真信號(hào)表達(dá)式如下:

其中，f1=200 Hz;f2=100 Hz;f3=50 Hz。采樣頻率為1 000 Hz。信號(hào)由3個(gè)正弦信號(hào)組成。圖2為原信號(hào)x(t)分解出來(lái)的3階IMF，分解時(shí)沒(méi)有對(duì)原始信號(hào)進(jìn)行任何處理，可見(jiàn)IMF1，IMF2，IMF3都產(chǎn)生了明顯的端點(diǎn)效應(yīng)，而且隨著分解過(guò)程的進(jìn)行，端點(diǎn)的發(fā)散程度越來(lái)越嚴(yán)重，在IMF3中甚至已經(jīng)無(wú)法表達(dá)出原信號(hào)中的低頻部分。

圖2 未經(jīng)處理的EMD分解結(jié)果

本文方法的具體步驟如下:

(1)將原信號(hào)x(t)進(jìn)行鏡像延拓處理，得到延拓后的信號(hào)y(t)。

(2)將延拓后的信號(hào)y(t)用一定周期ΔT的余弦窗函數(shù)α(t)進(jìn)行處理。即將y(t)與α(t)進(jìn)行內(nèi)積運(yùn)算，得到信號(hào) σ(t)=［x(t)，α(t)］，確定處理信號(hào)σ(t)的所有極值點(diǎn)。

(3)計(jì)算由局部極大值和極小值點(diǎn)確定的上下包絡(luò)線Mmax(t)和Mmin(t)。

(4)計(jì)算包絡(luò)線均值:

m(t)=(Mmax(t)+Mmin(t))/2

(5)求出h(t)=y(t)-m(t)。

(6)如果h(t)不滿足IMF的條件，則重復(fù)以上循環(huán);如果h(t)是一個(gè)IMF，將信號(hào)y1(t)=y(t)-h(t)作為原始信號(hào)重復(fù)以上循環(huán)。當(dāng)yk(t)成為一個(gè)單調(diào)信號(hào)時(shí)，循環(huán)結(jié)束。

4 仿真實(shí)驗(yàn)及結(jié)果分析

4．1 實(shí)驗(yàn)分析

4．1．1 延拓處理

采用余弦窗函數(shù)對(duì)原始信號(hào)進(jìn)行處理后，有可能會(huì)改變信號(hào)的特征，因此，改進(jìn)算法先對(duì)信號(hào)進(jìn)行鏡像延拓，使延拓信號(hào)和原始信號(hào)的交界處變得光滑。

本文將仿真信號(hào)x(t)作為延拓對(duì)象，在從左向右的第L個(gè)極值點(diǎn)處和從右向左的第R個(gè)極值點(diǎn)處分別設(shè)置2個(gè)光滑的反射面，將信號(hào)向外反射，得到數(shù)據(jù)序列長(zhǎng)度兩倍于反射信號(hào)長(zhǎng)度的周期性信號(hào)。這時(shí)，信號(hào)作為一個(gè)環(huán)形的閉合曲線，滿足了原信號(hào)曲線在端點(diǎn)處一階、二階倒數(shù)存在的要求。然后將閉合曲線作為一個(gè)整體進(jìn)行EMD分解，最終取上半部分作為分解結(jié)果。圖3所示結(jié)果為延拓后的信號(hào)，可以看出，此時(shí)信號(hào)兩端點(diǎn)處已呈一定的周期性。

圖3 鏡像延拓后的EMD分解結(jié)果

4．1．2 加窗處理

余弦窗是信號(hào)處理中應(yīng)用較為廣泛的一種窗函數(shù)，在本文中，余弦窗函數(shù)可被定義為:

式(6)的曲線圖如圖4所示，可見(jiàn)，在窗函數(shù)中間部分的值為1，窗函數(shù)兩端的值逐漸減少至0。這樣，信號(hào)端點(diǎn)處的值變?yōu)?，包絡(luò)線變得相對(duì)平滑，擬合出的包絡(luò)線也變得更為真實(shí)。隨著EMD分解的進(jìn)行，上下包絡(luò)線逐漸收斂于端點(diǎn)，發(fā)散現(xiàn)象就不會(huì)再發(fā)生。圖5為延拓后信號(hào)進(jìn)行余弦窗函數(shù)處理后的情況，可見(jiàn)IMF1，IMF2，IMF3的端點(diǎn)效應(yīng)均得到了明顯的抑制，而且IMF4分量中的低頻部分仍能表達(dá)出來(lái)。

圖4 余弦窗函數(shù)α(t)

圖5 加窗處理后的EMD分解結(jié)果

4．2 評(píng)價(jià)指標(biāo)

目前，所有提出的端點(diǎn)效應(yīng)抑制方法都需要一個(gè)定性與定量相結(jié)合的評(píng)價(jià)指標(biāo)來(lái)綜合考慮方法的優(yōu)劣。本文主要采用均方根有效值評(píng)價(jià)法來(lái)定性分析端點(diǎn)效應(yīng)抑制方法的效果:可以先計(jì)算出原始信號(hào)和通過(guò)篩選得到的各個(gè)IMF分量的均方根有效值。如式(7)所示。

其中，RMS表示信號(hào)有效值;s(i)為信號(hào)序列;n為信號(hào)的采樣點(diǎn)數(shù)。

按式(7)比較各IMF分量有效值的總和與原信號(hào)有效值，得到一個(gè)評(píng)價(jià)指標(biāo)θ。

其中，RMSoriginal表示原始信號(hào)的有效值;RMSi表示第i個(gè)IMF分量的有效值;n為IMF分量的個(gè)數(shù)，其中包括EMD的殘留項(xiàng)。

根據(jù)定義，θ≥0，當(dāng)端點(diǎn)效應(yīng)對(duì)EMD分解過(guò)程沒(méi)有影響時(shí)，θ=0。θ值越大，說(shuō)明端點(diǎn)效應(yīng)對(duì)EMD分解的影響越大。θ值越小，表明IMF分量的發(fā)散程度越小，端點(diǎn)效應(yīng)的影響也就越小。

4．3 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)給出了未加任何處理的EMD算法、極值對(duì)稱延拓法、鏡像延拓法、窗函數(shù)法和改進(jìn)算法對(duì)處理端點(diǎn)效應(yīng)的效果評(píng)價(jià)，給出了EMD分解前后的效果對(duì)比參數(shù)值 θ。從表1中可以看出，未加處理的EMD算法的θ值明顯偏大，極值對(duì)稱延拓法、鏡像延拓法、窗函數(shù)法的θ值都有不同程度的減小，而本文改進(jìn)算法的θ值較其他4種算法都小，說(shuō)明對(duì)端點(diǎn)效應(yīng)的抑制效果更好。

表1 多種算法的端點(diǎn)效應(yīng)抑制評(píng)價(jià)結(jié)果

5 結(jié)束語(yǔ)

本文提出了一種將鏡像延拓和余弦窗函數(shù)結(jié)合的端點(diǎn)效應(yīng)抑制方法。首先將原信號(hào)進(jìn)行鏡像延拓處理，然后在對(duì)延拓后信號(hào)進(jìn)行IMF提取之前，在其兩端加上余弦窗函數(shù)，使得上下包絡(luò)線能夠更好地?cái)M合。該方法彌補(bǔ)了鏡像延拓可能使延拓部分信號(hào)與原始信號(hào)的均值有明顯差異的不足，同時(shí)也解決了單純窗函數(shù)法因?yàn)楦淖冊(cè)盘?hào)而可能帶來(lái)其他問(wèn)題的缺陷。另外，本文利用θ值作為評(píng)價(jià)指標(biāo)綜合評(píng)價(jià)了多種抑制算法的效果，結(jié)果顯示出改進(jìn)方法較傳統(tǒng)方法略有提升。但是，在EMD分解過(guò)程中，由于每次篩選得出的新信號(hào)發(fā)散程度不同，如果使用單一窗函數(shù)對(duì)篩選后的新信號(hào)進(jìn)行處理可能還會(huì)帶來(lái)信號(hào)失真等問(wèn)題，隨著EMD的進(jìn)行，這種信號(hào)失真現(xiàn)象可能越來(lái)越嚴(yán)重。因此，對(duì)于EMD分解中的端點(diǎn)效應(yīng)問(wèn)題還有待進(jìn)一步研究。

［1］ Huang N E，Shen Zheng，Long S R，et al．The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis［J］．Proceedings ofRoyalSociety London，1998，454:903-995．

［2］ 王黎黎．基于希爾伯特-黃變換的時(shí)頻分析算法研究［D］．西安:西安電子科技大學(xué)，2009．

［3］ 秦 毅，秦樹人，毛永芳．改進(jìn)的Hilbert-Huang變換在信號(hào)瞬態(tài)特征提取中的應(yīng)用［J］．沖擊與振動(dòng)，2008，27(11):129-130．

［4］ Guo D，Peng Z P．Vibration Analysis of a Cracked Rotor Using Hilbert-Huang Transform［J］．Mechanical Systems and Signal Processing，2007，21(8):3030-3041．

［5］ Li X L，Bassiunya A M．Flute Breakage Detection During End Milling Using Hilbert-Huang Transform and Smoothed Nonlinear Energy Operator［J］．International Journal ofMachine Toolsand Manufacture，2007，47(6):1011-1020．

［6］ Feldman M．Considering High Harmonics for Identification of Non-linear System by Hilbert Transform［J］．MechanicalSystems and SignalProcessing，2007，21(2):943-958．

［7］ Yan Ruqiang，Gao R X．A Tour of the Hilbert-Huang Transform:An Empirical Tool for Signal Analysis［J］．IEEE Instrumentation ＆ Measurement Magazine，2007，10(5):40-45．

［8］ 張 超，陳建軍．EEMD方法和EMD方法抗模態(tài)混疊對(duì)比研究［J］．振動(dòng)與沖擊，2010，29(8):87-90．

［9］ 蘇玉香，劉志剛，李科亮，等．一種改善EMD端點(diǎn)效應(yīng)的新方法及其在諧波分析中的應(yīng)用［J］．電工電能新技術(shù)，2008，27(2):33-36．

［10］ Zhao Jinping，Huang Daji．Mirror Extending and Circular Spline Function for Empirical Mode Decomposition Method［J］．Journal of Zhejiang University，2001，2(3):247-252．

［11］ 顧小麗，李培良，譚海濤，等．基于 RBF神經(jīng)網(wǎng)絡(luò)的EMD方法在海平面分析中的應(yīng)用［J］．海洋與湖沼，2009，40(5):532-534．

［12］ 朱金龍，邱曉暉．正交多項(xiàng)式擬合在EMD算法端點(diǎn)問(wèn)題中的應(yīng)用［J］．計(jì)算機(jī)工程與應(yīng)用，2006，42(23):72-74．

［13］ 甘錫林，黃韋艮，楊勁松，等．基于希爾伯特-黃變換的合成孔徑雷達(dá)內(nèi)波參數(shù)提取新方法［J］．遙感學(xué)報(bào)，2007，11(1):39-40．

［14］ 張 博．基于Hilbert-Huang變換的滾動(dòng)軸承故障診斷方法研究［D］．西安:西安理工大學(xué)，2009．

［15］ 戴桂平．基于EMD的時(shí)頻分析方法研究［D］．秦皇島:燕山大學(xué)，2005．

［16］ Zeng Kan，HeMingxia．A SimpleBoundaryProcess Technique for Empirical Decomposition［C］//Proceedings of International Geoscience and Remote Sensing Symposium，September 20-24，2004，Anchorage，USA．Washington D．C．，USA:IEEE Press，2004:4258-4261．

［17］ 許寶杰，張建民，徐小力，等．抑制EMD端點(diǎn)效應(yīng)方法的研究［J］．北京理工大學(xué)學(xué)報(bào)，2006，26(3):196-200．

［18］ 黃大吉，趙進(jìn)平，蘇紀(jì)蘭．希爾波特-黃變換的端點(diǎn)延拓［J］．海洋學(xué)報(bào)，2003，25(1):1-9．

［19］ 任達(dá)千，吳昭同，嚴(yán)拱標(biāo)．EMD端點(diǎn)效應(yīng)的評(píng)價(jià)指標(biāo)及抑制端點(diǎn)效應(yīng)的窗函數(shù)法［J］．制造業(yè)自動(dòng)化，2007，29(1):21-24．