基于鏡像延拓和窗函數(shù)的端點效應(yīng)抑制方法

徐力彬，宋余慶，劉 毅

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

1 概述

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

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

2 EMD算法

2．1 算法簡介

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

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

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

EMD的分解過程如下:

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

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

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

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

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

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

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

2．2 端點效應(yīng)

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

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

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

2．3．1 常用延拓方法

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

(1)鏡像延拓法

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

(2)極值延拓法

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

(3)多項式擬合延拓法

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

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

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

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

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

3 本文算法

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

仿真信號表達式如下:

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

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

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

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

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

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

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

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

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

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

4 仿真實驗及結(jié)果分析

4．1 實驗分析

4．1．1 延拓處理

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

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

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

4．1．2 加窗處理

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

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

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

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

4．2 評價指標

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

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

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

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

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

4．3 實驗結(jié)果

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

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

5 結(jié)束語

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

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

［3］ 秦 毅，秦樹人，毛永芳．改進的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端點效應(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］ 朱金龍，邱曉暉．正交多項式擬合在EMD算法端點問題中的應(yīng)用［J］．計算機工程與應(yīng)用，2006，42(23):72-74．

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

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

［15］ 戴桂平．基于EMD的時頻分析方法研究［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端點效應(yīng)方法的研究［J］．北京理工大學(xué)學(xué)報，2006，26(3):196-200．

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

［19］ 任達千，吳昭同，嚴拱標．EMD端點效應(yīng)的評價指標及抑制端點效應(yīng)的窗函數(shù)法［J］．制造業(yè)自動化，2007，29(1):21-24．