噪聲環(huán)境中的EMD改進算法*

劉迎軍，楊志景，董健衛(wèi)，李淑龍

(1．南方醫(yī)科大學生物醫(yī)學工程學院，廣東 廣州 510515；2．廣東工業(yè)大學信息工程學院，廣東 廣州 510006；3．廣東藥學院基礎學院，廣東 廣州 510006)

1998年，美國工程院院士Huang 提出了一種新的信號分解方法——經驗模式分解(Empirical Mode Decomposition, EMD)方法[1]。該方法適用于非線性和非平穩(wěn)信號，能將信號自適應地分解成若干個本征模函數(Intrinsic Mode Function, IMF)之和，已經成功地被應用到各個領域[2-5]。然而，該方法理論基礎尚不完善，還存在一些問題有待進一步研究，比如邊界延拓、停止條件、模式混疊等等。對于確定性的信號，EMD的分解結果是合理的；但是對于含有噪聲的信號，經過EMD分解之后，原來的噪聲成分也被強制分解成了多個IMF之和，此時的每個IMF分量的物理意義不明確，甚至分解得到的噪聲分量個數也無法完全確定[6]，使得分解結果變得不可靠。為了解決該問題，本文提出一種改進的EMD方法：在首輪分解時，采用光滑樣條擬合來代替三次樣條插值。這樣可以避免對噪聲成分過度分解，從而使得結果更加可信。

1 EMD方法簡介

1.1 IMF的定義

EMD將信號分解成一系列的IMF之和，其中IMF的定義為滿足下列兩個條件的信號：① 整個信號中零點數與極值點數最多相差為1；② 信號上任意一點，由局部極大值點和局部極小值點確定的上、下包絡線的均值為零。

1.2 EMD算法

任給一個信號s(t), 首先對其所有極大值點和極小值點分別用三次樣條插值，得到上下兩條包絡線。用m(t)記為上下包絡線的平均值，令h(t)=s(t)-m(t)， 將h(t)當作新的s(t)，重復以上操作，直到h(t)滿足IMF的條件為止，則得到原信號s(t)的第一個IMF，即

c1(t)=h(t)

(1)

再令

r(t)=s(t)-c1(t)

(2)

然后將r(t)視為新的s(t)，重復以上過程，可依次得第2個IMFc2(t)，第3個IMFc3(t)，……，第N個IMFcN(t)，最終得分解式：

(3)

其中r(t)為余量信號，代表信號的平均趨勢。

2 改進的EMD方法

2.1 GCV-光滑樣條

設f(t)為未知的光滑函數且εi～N(0,σ2)，即εi為服從同一分布的高斯白噪聲，而數據點{(ti,yi)|i=1,2,…,n}由以下模型得到：

yi=f(ti)+εi

(4)

其中a

光滑樣條是很重要的一種數據擬合方法，最常用的定義式為：

(5)

其中，C2[a,b]表示區(qū)間[a,b]上二階連續(xù)可微的函數的集合，f″(t)為f(t)的二階導數。λ是一個光滑參數，用來權衡擬合程度與光滑程度。

如何選取光滑參數λ是個關鍵問題，已經有一些方法，比如交叉驗證(Cross Validation, CV)、廣義交叉驗證(Generalized Cross Validation, GCV)、Akaike信息準則等方法來選取合適的光滑參數。其中GCV是普遍使用的一種方法，計算公式如下

(6)

通過最小化公式(6) 中的GCV函數就能得到最優(yōu)光滑參數λ，從而得到相應的GCV-光滑樣條。

在模型(4)中，噪聲εi假設為高斯白噪聲。實際上，GCV-光滑樣條不僅適用于高斯白噪聲的情形，而且還對于一定程度內的高斯白噪聲和非白噪聲的混合情形也是適用的，而源自實際的數據序列經常可看作這種情形[7]，因此它被廣泛應用到許多領域。

2.2 噪聲IMF

如第1節(jié)所述，一個信號要成為IMF，必須滿足兩個條件，以保證IMF具有時間軸上的局部對稱性[8]。這與通常的周期信號、準周期信號，甚至復雜周期信號的特點是一致的；但是對于噪聲信號而言，這樣的條件就顯得過于苛刻。因為噪聲信號與確定性的信號不同，描述噪聲的IMF應該從概率分布的角度出發(fā)。于是，在考慮噪聲信號時，不適合采用樣條插值來得到包絡，而應該考慮擬合的手段，故可定義噪聲IMF滿足如下兩個條件：

1)噪聲IMF中極大值序列的均值曲線大于零，而極小值序列的均值曲線小于零；

2)噪聲IMF的極大值和極小值序列的均值曲線之和為零，即兩者關于X軸對稱。

如圖1所示，給定一個高斯白噪聲序列，過其極大值和極小值點序列作三次樣條插值，就得到了上下包絡(圖1(a))。對上下包絡求平均，則得到了均值曲線(圖1(c))。顯然，對于純高斯白噪聲序列，這樣的均值曲線沒有什么物理意義，反而會影響分解的最終結果。這是傳統(tǒng)EMD方法的一個缺陷。如果使用GCV-光滑樣條分別對極大和極小值點序列進行擬合，則得到近似的水平直線(圖1(b))。再將兩者平均，得到的還是近似水平直線(圖1(d))。這與高斯白噪聲的零均值是一致的。

圖1的例子表明，就純噪聲信號而言，對極值點序列采用擬合比插值更合適。在GCV-光滑樣條擬合的方法下，純噪聲序列不會被強制進行分解，而是作為一個獨立的分量而存在，即為噪聲IMF。

圖1 圖中虛線為仿真產生的純高斯白噪聲Fig.1 The dotted curve is a pure Gaussian noise sequence(a)實線分別為過極大值和極小值點序列的三次樣條插值所求得的上下包絡； (b)實線分別為擬合極大值和極小值點序列的GCV-光滑樣條；(c)實曲線為(a)中上下包絡的均值； (d)實線為(b)上下GCV-光滑樣條的均值

2.3 改進的EMD算法

改進的EMD算法具體描述如下：

1)任給信號s(t)，對其進行滑動平均得sa(t)(滑動窗寬取為5)，然后計算

(7)

如果T>0.05，說明是信號含有噪聲，轉②；否則對s(t)進行原始EMD分解，程序終止；

2)求出s(t)的極值點序列, 分別對其極大值點和極小值點序列用GCV-光滑樣條擬合，得到上下兩條極值點序列分布的均值曲線ek1(t)，ek2(t)；

3)令mk(t):=[ek1(t)+ek2(t)]/2，hk(t):=s(t)-mk(t)，若相鄰兩次的hk(t)和hk-1(t)滿足停止條件(ρ經驗取值為0.05)

(8)

就得到原信號s(t)的第一個IMF，也就是噪聲IMFc1(t)；否則，將hk(t)當作新的s(t)，轉②；

4)令s(t):=s(t)-c1(t)，接著對s(t)進行原始EMD分解。

從上述算法不難看出，改進的EMD方法與原始EMD的區(qū)別在于對于含有噪聲的信號的第1個IMF的求法不同。

3 仿真實驗

3.1 仿真數據

基于模型(4)，每個仿真數據由高斯白噪聲和不同趨勢疊加而成，其中高斯白噪聲的均值為0，方差為1，數據長度取n=600。正弦趨勢分別?。阂环N為單頻率的正弦信號y=cos (18πt/n),t=1,2,…,n；另一種為兩個頻率的正弦信號的疊加，即y=cos (18πt/n)+sin(10πt/n),t=1,2,…,n。

3.2 實驗結果

當原信號為正弦信號y=cos (18πt/n)與高斯白噪聲疊加而成時，EMD和改進的EMD分解結果分別顯示在圖2和圖3。在圖2中，前面4個IMF成分是噪聲成分，而圖3只有IMF1才是噪聲成分。圖2中的IMF5信號與圖3中的IMF2頻率相似，它們都是對原始信號中正弦趨勢的刻畫。但不難發(fā)現，圖2中的IMF5在時間坐標區(qū)間 [100, 200] 之間已經嚴重偏離原始的正弦信號，而圖3中IMF2則基本上保留了原始的正弦信號的特征。由此可以得出，改進的EMD算法避免了對噪聲成分的過度分解，一方面使得結果變得簡潔，另一方面極大地減弱了噪聲成分對分解造成的干擾。

圖2 EMD分解結果Fig.2 The results of EMD最上方為原信號，由正弦信號y=cos(18πt/n)與高斯白噪聲疊加而成。下面依次是分解產生的各個IMF，最后一個是余量

圖3 改進的EMD分解結果Fig.3 The results of the proposed improved EMD最上方為原信號，由正弦信號y=cos(18πt/n)與高斯白噪聲疊加而成。下面依次是分解產生的各個IMF，最后一個是余量

圖4和圖5則給出了原信號為正弦信號y=cos (18πt/n)+sin (10πt/n)與高斯白噪聲疊加而成時，EMD和改進的EMD分解結果。圖4中的IMF5和IMF6分別與圖5中的IMF2和IMF3對應，它們都是對原始信號中高低頻率的兩個正弦趨勢的刻畫。仔細觀察，容易發(fā)現，圖4中的IMF5在末端處有衰減，而IMF6在起始處也有衰減，其原因可能是受到了噪聲成分的干擾；而圖5中IMF2和IMF3與原始的兩個正弦信號保持很強的一致性。該仿真實驗進一步證實了改進的EMD算法比原EMD算法的更為可靠。

以上仿真實驗結果表明，本文提出的改進EMD算法，對含有噪聲的信號分解結果有很大的改進。

4 兩個實際數據例子

EMD方法很早就被應用于氣候時間序列的分析(如文獻[1,9])。然而，由于影響氣候的因素過于復雜，在數據分析處理時需要十分謹慎。為了驗證本文方法對EMD的改進效果，我們兩種方法同時應用到兩個實際的氣候數據例子上。一個數據例子是1979-2012年緯向平均的每月的500 mb的溫度異常值序列，另一個則是相同時間范圍內的200 mb赤道緯向風(165°W-110°W)的原始數據。它們都源自網站http://www.cpc.ncep.noaa.gov/data/indices/，具體分解結果如圖6至圖9所示。

對比圖6和圖7，可知改進的EMD方法結果更簡潔，所有噪聲成分都集中在噪聲IMF中(圖7中IMF1)，而EMD方法結果中前兩個IMF都是噪聲成分。通過進一步觀察，不難發(fā)現圖7中IMF3-6與圖6中的IMF4-7大致相對應，但是周期性規(guī)律更為顯著。圖6中的IMF5不同周期的長度變化劇烈，尤其在1985年到1995年之間的單個周期跨度太大，明顯不合理，很可能是受到了噪聲的嚴重干擾；而圖7中的IMF4表現出合理的周期節(jié)律，而且幅度值更大，結果更為可信。

圖4 EMD分解結果Fig.4 The results of EMD最上方為原信號，由正弦信號y=cos (18πt/n)+sin (10πt/n)與高斯白噪聲疊加而成；下面依次是分解產生的各個IMF，最后一個是余量

圖5 改進的EMD分解結果Fig.5 Results of the proposed improved EMD最上方為原信號，由正弦信號y=cos (18πt/n)+sin (10πt/n)與高斯白噪聲疊加而成；下面依次是分解產生的各個IMF，最后一個是余量

圖7 緯向平均的500 mb的溫度異常值序列的改進的EMD分解結果圖Fig.7 The decomposition of the zonally average 500 mb temperature anomalies signal by the proposed improved EMD最上面是原始信號序列，下面依次是分解得到的一系列IMF分量，最后一個是余量

圖8和圖9分別顯示了不同方法對200 mb赤道緯向風的原始數據分解的結果。從圖8最上方的原始數據來看，有明顯的1年的周期性規(guī)律，EMD的分解結果中IMF2與該規(guī)律最為相近，但差異較大，主要是其中夾雜了更快的小幅度震蕩波。顯而易見，這個現象是典型的“模式混疊”[10]；相比之下，圖9中改進的EMD方法中的IMF2忠實地反映了周期為1年的基本規(guī)律。

仔細對比圖8和圖9中其他IMF的表現，不難發(fā)現圖9中IMF周期性明顯更清晰：以連續(xù)兩個極大或極小值之間的長度近似為一個周期，可計算從IMF3到IMF7平均周期分別近似為2.5年，4.8 年，11年和23年。

圖8 200 mb赤道緯向風的原始數據的EMD分解結果圖Fig.8 The EMD decomposition result of the 200 mb Zonal Winds Equator (165°W-110°W)最上面是原始信號序列，下面依次是EMD分解一系列IMF分量，最后一個是余量

5 結 論

EMD是一種適用于非平穩(wěn)和非線性的信號分析方法，但它容易受噪聲成分的干擾，常常表現出不穩(wěn)定，從而導致分解結果不可靠。雖然文獻[11]提出了一種利用IMF的能量密度與其平均周期乘積的分布來判別具體的IMF是噪聲還是信號成分，但是EMD本身無法保證將噪聲成分與信號成分很好的分離，比如圖8中IMF2就屬于這種情況，“模式混疊”是EMD分解中經常碰到的一個問題。如果通過其它手段抑制噪聲成分，則會使得EMD分解更為有效。正是基于這一想法，本文提出了結合GCV-光滑樣條的改進的EMD方法。該方法能克服原EMD方法對信號中噪聲成分的過度分解的弊端，不僅使得分解結果更加簡潔，而且極大地抑制了噪聲對信號分解的干擾，從而使得結果更加可信。仿真實驗和實際氣候數據例子都證實了新提出的方法的有效性和優(yōu)越性。

[1] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J]. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1998, 454: 903-995.

[2] PINES D, SALVINO L. Structural health monitoring using empirical mode decomposition and the Hilbert phase [J]. Journal of sound and vibration, 2006, 294(1): 97-124.

[3] YUAN J S, ZHAO Z B, GAO Q, et al. Multimodal image registration based on empirical mode decomposition and mutual information [J]. Chinese Journal of Scientific Instrument, 2009, 10: 012.

[4] 邵駿, 袁鵬, 顏志衡, 等. 基于 HHT 的雅魯藏布江徑流變化周期及趨勢分析[J]. 中山大學學報：自然科學版, 2010, 49(1): 125-130.

圖9 200 mb赤道緯向風的原始數據改進的EMD方法分解結果Fig.9 The decomposition of the 200 mb Zonal Winds Equator (165°W-110°W) by the proposed improved EMD最上面是原始信號序列，下面依次是IMF分量，最后一個是余量

[5] WEI Y, CHEN M C. Forecasting the short-term metro passenger flow with empirical mode decomposition and neural networks [J].Transportation Research Part C: Emerging Technologies, 2012, 21(1): 148-162.

[6] FLANDRIN P, GON?ALVES P, RILLING G. Detrending and denoising with empirical mode decomposition [C]//Proceedings of the 12th European Signal Processing Conference (EUSIPCO’04), Vienna, Austria, September 2004.

[7] WAHBA G. Spline models for observational data [M]. Philadelphia: SIAM, 1990.

[8] 鐘佑明, 秦樹人. 希爾伯特-黃變換的統(tǒng)一理論依據研究[J]. 振動與沖擊, 2006, 25(3): 40-43.

[9] COUGHLIN K T, TUNG K K. 11-year solar cycle in the stratosphere extracted by the empirical mode decomposition method [J]. Advances in space research, 2004, 34(2): 323-329.

[10] HU X Y, PENG S L, HWANG W L. EMD revisited: A new understanding of the envelope and resolving the mode-mixing problem in AM-FM signals [J]. IEEE Transactions on Signal Processing, 2012, 60(3): 1075-1086.

[11] WU Z H, HUANG N E. A study of the characteristics of white noise using the empirical mode decomposition method [J]. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2004, 460: 1597-1611.