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

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

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

1 EMD方法簡(jiǎn)介

1.1 IMF的定義

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

1.2 EMD算法

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

c1(t)=h(t)

(1)

再令

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

(2)

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

(3)

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

2 改進(jìn)的EMD方法

2.1 GCV-光滑樣條

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

yi=f(ti)+εi

(4)

其中a

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

(5)

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

如何選取光滑參數(shù)λ是個(gè)關(guān)鍵問(wèn)題，已經(jīng)有一些方法，比如交叉驗(yàn)證(Cross Validation, CV)、廣義交叉驗(yàn)證(Generalized Cross Validation, GCV)、Akaike信息準(zhǔn)則等方法來(lái)選取合適的光滑參數(shù)。其中GCV是普遍使用的一種方法，計(jì)算公式如下

(6)

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

在模型(4)中，噪聲εi假設(shè)為高斯白噪聲。實(shí)際上，GCV-光滑樣條不僅適用于高斯白噪聲的情形，而且還對(duì)于一定程度內(nèi)的高斯白噪聲和非白噪聲的混合情形也是適用的，而源自實(shí)際的數(shù)據(jù)序列經(jīng)?？煽醋鬟@種情形[7]，因此它被廣泛應(yīng)用到許多領(lǐng)域。

2.2 噪聲IMF

如第1節(jié)所述，一個(gè)信號(hào)要成為IMF，必須滿足兩個(gè)條件，以保證IMF具有時(shí)間軸上的局部對(duì)稱(chēng)性[8]。這與通常的周期信號(hào)、準(zhǔn)周期信號(hào)，甚至復(fù)雜周期信號(hào)的特點(diǎn)是一致的；但是對(duì)于噪聲信號(hào)而言，這樣的條件就顯得過(guò)于苛刻。因?yàn)樵肼曅盘?hào)與確定性的信號(hào)不同，描述噪聲的IMF應(yīng)該從概率分布的角度出發(fā)。于是，在考慮噪聲信號(hào)時(shí)，不適合采用樣條插值來(lái)得到包絡(luò)，而應(yīng)該考慮擬合的手段，故可定義噪聲IMF滿足如下兩個(gè)條件：

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

2)噪聲IMF的極大值和極小值序列的均值曲線之和為零，即兩者關(guān)于X軸對(duì)稱(chēng)。

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

圖1的例子表明，就純?cè)肼曅盘?hào)而言，對(duì)極值點(diǎn)序列采用擬合比插值更合適。在GCV-光滑樣條擬合的方法下，純?cè)肼曅蛄胁粫?huì)被強(qiáng)制進(jìn)行分解，而是作為一個(gè)獨(dú)立的分量而存在，即為噪聲IMF。

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

2.3 改進(jìn)的EMD算法

改進(jìn)的EMD算法具體描述如下：

1)任給信號(hào)s(t)，對(duì)其進(jìn)行滑動(dòng)平均得sa(t)(滑動(dòng)窗寬取為5)，然后計(jì)算

(7)

如果T>0.05，說(shuō)明是信號(hào)含有噪聲，轉(zhuǎn)②；否則對(duì)s(t)進(jìn)行原始EMD分解，程序終止；

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

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

(8)

就得到原信號(hào)s(t)的第一個(gè)IMF，也就是噪聲IMFc1(t)；否則，將hk(t)當(dāng)作新的s(t)，轉(zhuǎn)②；

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

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

3 仿真實(shí)驗(yàn)

3.1 仿真數(shù)據(jù)

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

3.2 實(shí)驗(yàn)結(jié)果

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

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

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

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

以上仿真實(shí)驗(yàn)結(jié)果表明，本文提出的改進(jìn)EMD算法，對(duì)含有噪聲的信號(hào)分解結(jié)果有很大的改進(jìn)。

4 兩個(gè)實(shí)際數(shù)據(jù)例子

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

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

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

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

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

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

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

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

5 結(jié) 論

EMD是一種適用于非平穩(wěn)和非線性的信號(hào)分析方法，但它容易受噪聲成分的干擾，常常表現(xiàn)出不穩(wěn)定，從而導(dǎo)致分解結(jié)果不可靠。雖然文獻(xiàn)[11]提出了一種利用IMF的能量密度與其平均周期乘積的分布來(lái)判別具體的IMF是噪聲還是信號(hào)成分，但是EMD本身無(wú)法保證將噪聲成分與信號(hào)成分很好的分離，比如圖8中IMF2就屬于這種情況，“模式混疊”是EMD分解中經(jīng)常碰到的一個(gè)問(wèn)題。如果通過(guò)其它手段抑制噪聲成分，則會(huì)使得EMD分解更為有效。正是基于這一想法，本文提出了結(jié)合GCV-光滑樣條的改進(jìn)的EMD方法。該方法能克服原EMD方法對(duì)信號(hào)中噪聲成分的過(guò)度分解的弊端，不僅使得分解結(jié)果更加簡(jiǎn)潔，而且極大地抑制了噪聲對(duì)信號(hào)分解的干擾，從而使得結(jié)果更加可信。仿真實(shí)驗(yàn)和實(shí)際氣候數(shù)據(jù)例子都證實(shí)了新提出的方法的有效性和優(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 的雅魯藏布江徑流變化周期及趨勢(shì)分析[J]. 中山大學(xué)學(xué)報(bào)：自然科學(xué)版, 2010, 49(1): 125-130.

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

[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] 鐘佑明, 秦樹(shù)人. 希爾伯特-黃變換的統(tǒng)一理論依據(jù)研究[J]. 振動(dòng)與沖擊, 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.