徐力彬，宋余慶，劉 毅

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

1 概述

Hilbert-Huang變換(Hilbert-Huang Transform，HHT)是由Huang N E等人通過深入分析和總結(jié)后提出的一種新的非平穩(wěn)非線性信號分析方法［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無須預(yù)先設(shè)定任何基函數(shù)，而是根據(jù)信號特點(diǎn)自適應(yīng)地將原始信號分解為一系列固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，去除了快速傅里葉變換方法中引入的諧波基函數(shù)，也沒有小波分析復(fù)雜信號過程中產(chǎn)生信號模糊與失真的問題［3-6］。然而，HHT變換也有一定的局限性，例如端點(diǎn)效應(yīng)問題［7］和模態(tài)混疊問題［8］等。尤其是端點(diǎn)效應(yīng)問題會嚴(yán)重妨礙信號的特征提取，從而導(dǎo)致信號分析的誤差。

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

2 EMD算法

2．1 算法簡介

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

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

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

EMD的分解過程如下:

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

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

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

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

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

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

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

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

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

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

2．3．1 常用延拓方法

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

(1)鏡像延拓法

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

(2)極值延拓法

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

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

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

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

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

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

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

3 本文算法

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

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

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

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

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

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

(2)將延拓后的信號y(t)用一定周期ΔT的余弦窗函數(shù)α(t)進(jìn)行處理。即將y(t)與α(t)進(jìn)行內(nèi)積運(yùn)算，得到信號 σ(t)=［x(t)，α(t)］，確定處理信號σ(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，將信號y1(t)=y(t)-h(t)作為原始信號重復(fù)以上循環(huán)。當(dāng)yk(t)成為一個(gè)單調(diào)信號時(shí)，循環(huán)結(jié)束。

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

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

4．1．1 延拓處理

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

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

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

4．1．2 加窗處理

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語

本文提出了一種將鏡像延拓和余弦窗函數(shù)結(jié)合的端點(diǎn)效應(yīng)抑制方法。首先將原信號進(jìn)行鏡像延拓處理，然后在對延拓后信號進(jìn)行IMF提取之前，在其兩端加上余弦窗函數(shù)，使得上下包絡(luò)線能夠更好地?cái)M合。該方法彌補(bǔ)了鏡像延拓可能使延拓部分信號與原始信號的均值有明顯差異的不足，同時(shí)也解決了單純窗函數(shù)法因?yàn)楦淖冊盘柖赡軒砥渌麊栴}的缺陷。另外，本文利用θ值作為評價(jià)指標(biāo)綜合評價(jià)了多種抑制算法的效果，結(jié)果顯示出改進(jìn)方法較傳統(tǒng)方法略有提升。但是，在EMD分解過程中，由于每次篩選得出的新信號發(fā)散程度不同，如果使用單一窗函數(shù)對篩選后的新信號進(jìn)行處理可能還會帶來信號失真等問題，隨著EMD的進(jìn)行，這種信號失真現(xiàn)象可能越來越嚴(yán)重。因此，對于EMD分解中的端點(diǎn)效應(yīng)問題還有待進(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變換在信號瞬態(tài)特征提取中的應(yīng)用［J］．沖擊與振動，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)混疊對比研究［J］．振動與沖擊，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)問題中的應(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］．西安:西安理工大學(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)的評價(jià)指標(biāo)及抑制端點(diǎn)效應(yīng)的窗函數(shù)法［J］．制造業(yè)自動化，2007，29(1):21-24．