消除EMD邊界效應的四中點估計方法

張立振, 劉坦然, 魏恩泊

(1. 中國海洋大學 數(shù)學科學學院, 山東 青島 266100; 2. 中國科學院 海洋研究所, 山東 青島 266071)

消除EMD邊界效應的四中點估計方法

張立振1, 劉坦然1, 魏恩泊2

(1. 中國海洋大學 數(shù)學科學學院, 山東 青島 266100; 2. 中國科學院 海洋研究所, 山東 青島 266071)

為了快速高效地實現(xiàn)信號的本征模態(tài)分解, 消除分解過程中邊界效應, 在 Directly-Mean EMD方法的基礎上, 通過引入左、右待估中點和左、右延拓中點, 并根據(jù)信號首尾兩端各種可能的情況給出了四中點估計公式, 建立了可完全消除邊界效應的Directly-Mean EMD四中點估計方法。該方法不僅減少了樣條插值次數(shù), 提高分解速度, 而且還可以有效避免因插值節(jié)點過于稀疏所產(chǎn)生的大幅波動, 使分解結(jié)果更加準確。將新方法應用于日長數(shù)據(jù)序列的本征模態(tài)函數(shù)分解, 得到了滿意的分解結(jié)果。

本征模態(tài)函數(shù); Directly-Mean EMD方法; 四中點估計方法; 希爾伯特－黃變換; 日長(LOD)

1996年開始, Huang等[1]對HHT(希爾伯特－黃變換)又稱 EMD/HSA(經(jīng)驗模態(tài)分解與希爾伯特譜分析)進行了系統(tǒng)研究。HHT可以有效地分析非線性、非平穩(wěn)信號。近幾年來, 其應用越來越廣泛。分解信號為本征模態(tài)函數(shù)是 HHT分析中最為關鍵的一環(huán)。Huang給出了信號分解的sifting方法, 這種方法從信號本身的特征尺度出發(fā), 具有很好的自適應性。但是, 在用樣條函數(shù)擬合信號上、下包絡時會出現(xiàn)嚴重的邊界效應, 這將嚴重影響據(jù)此所做出的 Hilbert譜。目前有很多方法來改善邊界效應,例如已經(jīng)向NASA申請專利的特征波法, 黃大吉等[2]給出了鏡像閉合延拓法和極值延拓法, 沈路[3]研究了相似極值延拓法, 劉慧婷[4]提出了基于多項式擬合的處理方法, 胡愛軍, 孫華剛, 曹沖鋒等[5-7]也各自提出了有效的邊界處理方法。

Zhong等[8]利用極值點與線性包絡的“中點”作為型值點擬合出樣條曲線, 將其直接代替 Sifting方法中上、下包絡的平均值曲線, 從而建立一種本征模態(tài)函數(shù)分解的Directly-Mean EMD方法。該方法不僅減少了插值次數(shù), 而且因為增加了型值點數(shù), 使樣條曲線的精度得到提高。但是, Zhong等并沒有涉及對邊界效應的處理, 也沒有討論當信號極值點個數(shù)少于4時如何求取中點曲線的問題。本文便是在Directly-Mean EMD方法的基礎上, 通過引入左右待估中點和左右延拓中點的概念及其估計方法, 不僅使 Directly-Mean EMD可以完全消除分解時邊界效應, 同時針對極值點個數(shù)少于4時的各種情況, 給出了中點曲線的求取方法。

1 本征模態(tài)函數(shù)與Sifting方法

所謂本征模態(tài)函數(shù)是滿足以下兩條件的數(shù)據(jù)集: (1)在整個數(shù)據(jù)集上, 其極值點的個數(shù)和跨零點的個數(shù)相同或至多相差為1, 即在相鄰的兩個上跨零點之間只有一個極大值點和一個極小值點。(2)由極大值點所定義的上包絡與由極小值點所定義的下包絡的平均值為0。 Huang創(chuàng)立的將信號分解成IMF的sifting方法[1], 基本步驟為首先求出信號的局部極大(小)值點的包絡(樣條插值)。然后求出上、下包絡的平均值,并將其從信號中減掉。重復此操作直至得到本征模態(tài)函數(shù), 將此本征模態(tài)函數(shù)從原信號分離出去。重復上面的操作, 直至剩余信號為只有一個極值點或單調(diào)函數(shù)為止。

2 Directly-Mean EMD方法

從 Sifting方法的過程可以看到, 在求信號的上下包絡時用到了兩次樣條插值, 然后再求平均值。Directly-Mean EMD方法[8]只用一次樣條插值, 具體的方法如下:

步驟 1: 求取信號x(t)的極值點(極大值點和極小值點)。

步驟2: 求取中點曲線。用直線段將極大值點順次連接起來, 稱為信號的線性上包絡(圖1)。用直線段將極小值點順次連接起來, 稱為信號的線性下包絡。每個“▽點”位于“線性上包絡”上, 發(fā)生的時刻與相應的極小值點一致, 例如“下三角點 D”與相應的極小點C發(fā)生在同一時刻。類似地, 每個“△點”位于“線性下包絡線”上, 發(fā)生的時刻與相應的極大值點一致。設極大值點 A, 極大值點B, 極小值點 C三點的坐標分別為:。容易推出“下三角點 D”的坐標為:

將“下三角點D”與極小值點C連線的中點E稱為信號與極小值點C對應的中點, 其坐標為:

每個極值點都可以類似地定義相對應的中點(圖1中用o表示的點)。易知用直線段將所有中點連接起來得到的折線恰為上、下線性包絡的平均值。由于信號的真實包絡不一定正好是線性的, 所以, 將所有中點作為型值點擬合出的樣條曲線m1(t)作為真實包絡平均值的估計, 并稱其為信號的中點曲線。

圖1 Directly-Mean EMD方法幾個主要概念Fig. 1 Important concepts in Directly-Mean EMD method

于是得到x(t)的本征模態(tài)函數(shù)分解式:

上述分解過程即為Directly-Mean EMD方法。

甜高粱與苜蓿混合不同比例的青貯飼料，甜高粱采自塔里木大學動物科學學院實驗站，苜蓿采自第一師五團。甜高粱與苜蓿青貯原料的化學成分見表1。

3 消除邊界效應的四中點估計方法

在運用Directly-Mean EMD方法進行本征模態(tài)函數(shù)分解時, 求出信號的中點曲線是非常關鍵的一步。而要得到中點曲線首先要確定所有中點的坐標。在許多情況下, 由于信號的起點或終點并不適合當作極值點來處理, 所以信號兩端的中點不能如Directly-Mean EMD方法中的步驟2確定。例如: 圖2中與第一極大值點B相對應的中點F, 其縱坐標只能通過某種方法給出估計值。類似地, 還有與最后極大值點對應的中點H。以下分別將這兩個中點稱為“左待估中點F”和“右待估中點H”。利用左右待估中點以及其他中點(圖 2中用“o”表示的點)便可擬合出樣條曲線, 即中點曲線m1(t)。由于這些中點所跨橫坐標區(qū)間一般短于原信號橫坐標區(qū)間,這樣得到的中點曲線在信號的兩端常常會發(fā)生所謂的“飛翼”現(xiàn)象(圖1中m1(t)在左、右兩端分別發(fā)生了過度下行和上翹的現(xiàn)象, 均嚴重偏離了其平均變化趨勢)。鑒于此, 有必要在信號兩端再補充兩個中點(圖2中用“◇”表示), 分別稱之為“左延拓中點G”和“右延拓中點J”。以所有中點(用“o”表示的點)以及左、右待估中點, 左、右延拓中點擬合出的樣條曲線m1(t)將會完全消除“飛翼”現(xiàn)象(圖 2短虛線)。下面討論左、右待估中點及左、右延拓中點坐標的估計方法, 以便較準確地求取中點曲線m1(t)。

3.1 含有四或四個以上極值點信號的中點曲線求取方法

分四種情況加以討論: (I)第一極值為極大值, 最后極值為極大值; (II)第一極值為極大值, 最后極值為極小值; (III)第一極值為極小值, 最后極值為極大值; (IV)第一極值為極小值, 最后極值為極小值。首先討論第一種情況(I), 見圖2。

圖2 邊界處理中幾個重要概念Fig. 2 Important concepts in boundary treatment

關于右待估中點H與右延拓中點J坐標的估計方法可仿此進行。

如此, 根據(jù)所有中點“o”, 左、右待估中點及左、右延拓中點共同擬合可得到樣條曲線。一般地, 這樣得到的樣條曲線會比原信號長, 將樣條曲線介于起點A和終點K之間的部分作為中點曲線m1(t), 此時的m1(t)將會完全消除“飛翼”現(xiàn)象。對于(II)(III)(IV)三種情況, 其估計方法可仿此給出。

3.2 含有三極值點信號的中點曲線求取方法

分兩種情況加以討論: (I)第一極值為極小值; (II)第一極值為極大值。

對于第一種情況(I), 見圖 3。首先按(7)(8)和(9)(10)給出左、右待估中點G和H的坐標估計:

然后按(11)(12)和(13)(14)分別給出左、右延拓中點J和K的坐標估計:

圖3 含有三個極值點時求取中點曲線Fig. 3 Midpoint curve calculation in case three extreme points existed

這樣利用左、右延拓中點J、K, 左、右待估中點G、H以及中點F五點共同擬合得到樣條曲線, 并將樣條曲線介于原信號起點和終點的部分作為中點曲線m1(t)。關于第二種情況(II), 可仿此給出。

3.3 含有兩極值點信號的中點曲線求取方法

當信號恰有兩個極值點時, 也分兩種情況加以討論: (I)第一極值為極小值; (II)第一極值為極大值。對于第一種情況(I), 見圖 4。先按(15)(16)和(17)(18)分別給出左、右待估中點E和F的坐標估計:

然后, 將過E、F兩點且介于起點A和終點D之間的直線段作為中點曲線m1(t)。對于第二種情況(II)可類似規(guī)定。

3.4 對于只含有一個極值點或沒有極值點的信號, 這時不必求取中點曲線, 可直接將其劃歸剩余信號

為敘述方便, 將根據(jù)上述各種求取中點曲線并按Directly-Mean EMD所述步驟進行本征模態(tài)函數(shù)分解的過程稱為: Directly-Mean EMD四中點估計法。

無論用 Sifting方法, 還是 Directly-Mean EMD方法, 之所以出現(xiàn)“飛翼”現(xiàn)象——邊界效應, 其原因就是擬合樣條曲線所用型值點兩端橫坐標區(qū)間是原信號橫坐標區(qū)間的真子集造成的。而根據(jù)上述求取中點曲線的過程可以看出, 左右延拓中點所跨橫坐標區(qū)間是原信號橫坐標區(qū)間的真子集的情況永遠不會發(fā)生。因此, Directly-Mean EMD四中點估計法可以完全消除邊界效應。

圖4 含有兩個極值點時求取中點曲線Fig. 4 Midpoint curve calculation in case two extreme points existed

圖5 利用Directly-Mean EMD四中點估計法分解日長信號Fig. 5 LOD signal decomposition by using Directly-Mean EMD four-midPoint estimation method

4 用Directly-Mean EMD四中點估計法對日長數(shù)據(jù)序列進行分解

日長(LOD, length of day)數(shù)據(jù)序列是綜合利用激光衛(wèi)星(SLR, Satellite Laser Ranging),激光測月(LLR, Lunar Laser Ranging), 長基線干涉術(VLBI, Very-Long-Baseline Interferometry)和全球定位系統(tǒng)(Global Positioning System)等多項先進技術測得的結(jié)果。覆蓋時間從1962年1月20日到2001年1月6日, 采樣間隔為1d, 共包含14 232 d記錄(圖5(1))。該數(shù)據(jù)序列的縱軸值代表測得的日長與標準日長(86 400 s)的差(單位ms)。從圖上容易看出, 這些年來的日長都比所謂的標準日長要長幾個毫秒。將Directly-Mean EMD四中點估計法應用于此數(shù)據(jù)序列的分解, 共得到 11個本征模態(tài)函數(shù)(圖 5(2)～(12))和一個剩余項(圖 5(13))。其中, 第一本征模態(tài)函數(shù)(圖5(2))的平均周期約為13.6 d, 為半個交點月周期。第二本征模態(tài)函數(shù)(圖5(3))的平均周期約為27.55 d,為一個近點月周期。這兩模態(tài)的物理意義十分明顯。第六本征模態(tài)函數(shù)(圖 5(11))較好地反映了其低頻變化趨勢。其他模態(tài)項的物理意義有待進一步研究。從分解結(jié)果來看, Directly-Mean EMD四中點估計法不僅從原數(shù)據(jù)序列中分離出兩個具有明顯物理意義的本征模態(tài)函數(shù), 而且邊界效應也處理得很好。

5 結(jié)論與討論

Directly-Mean EMD四中點估計法與 Sifting兩種分解方法都是從原信號中減掉平均曲線, 逐漸得到上下包絡關于零點對稱的本征模態(tài)函數(shù)。只是求取平均曲線的方法不同。Sifting方法是對極大值點和極小值點分別進行兩次樣條插值, 并對所得兩樣條曲線取平均的方法得到平均曲線; 而Directly-Mean EMD四中點估計法方法則通過先求取信號的中點, 再用這些中點直接擬合出中點曲線,即平均曲線。通過實例比較, 兩種方法在不考慮邊界效應的情況下, 所得平均曲線基本一致。除此之外, Directly-Mean EMD四中點估計法還具有如下優(yōu)點:

(1)Directly-Mean EMD方法由于只需要一次樣條插值計算, 因而減少了計算量, 可以提高分解速度。(2)在進行樣條插值時, Directly-Mean EMD四中點估計法參與插值的中點數(shù)幾乎是Sifting方法參與插值的極大值點數(shù)或參與插值的極小值點數(shù)的兩倍。因而, Directly-Mean EMD四中點估計法可以有效避免因插值節(jié)點過于稀疏所產(chǎn)生的不必要大幅度波動。從而使所得中點曲線更加精確。

Directly-Mean EMD 四 中 點估計法 是 對Directly-Mean EMD的進一步完善。與Sifting方法相比不僅具有上述(1)(2)兩個優(yōu)點, 還通過引入左、右待估中點和左右延拓中點及其估計, 完全消除了分解過程的邊界效應。同時, 對信號中含有較少極值點的情況進行了討論, 為Directly- Mean EMD的具體實施制定了統(tǒng)一的處理規(guī)則。

[1] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition method and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proc R Soc, 1998, 454: 903-995.

[2] 黃大吉, 趙進平, 蘇紀蘭. 希爾伯特-黃變換的端點延拓[J]. 海洋學報, 2003, 25(1): 1-11.

[3] 沈路, 周曉軍, 張志剛, 等. Hilbert-Huang變換中的一種端點延拓方法[J]. 振動與沖擊, 2009, 28(8): 168-171.

[4] 劉慧婷, 旻張 , 程家興. 基于多項式擬合算法的 EMD端點問題的處理[J]. 計算機工程與應用, 2004, 16: 84-86.

[5] 胡愛軍, 安連鎖, 唐貴基, 等. Hilbert-Huang變換端點效應處理新方法[J]. 機械工程學報, 2008, 44(4): 154-158.

[6] 孫華剛, 馮廣斌, 曹登慶, 等. 經(jīng)驗模態(tài)分解端點波形延拓改進方法研究[J]. 電子測量與儀器學報, 2010, 24(4): 319-326.

[7] 曹沖鋒, 楊世錫, 楊將新, 等. 一種抑制EMD端點效應新方法及其在信號特征提取中的應用[J]. 振動工程學報, 2008, 21(6): 588-593.

[8] Zhong Youming, Qin Shuren, Tang Baoping. Research on theoretic evidence and realization of Directly-Mean EMD method[J]. Chinese Journal of mechanical engineering, 2004, 17(3): 399-404.

(本文編輯: 劉珊珊)

Four-midpoint estimation method removing boundary effect of EMD

ZHANG Li-zhen1, LIU Tan-ran1, WEI En-bo2
(1.College of Mathematical Science, Ocean University of China, Qingdao 266100, China; 2. Institute of Oceanology, the Chinese Academy of Sciences, Qingdao 266071, China)

Jan., 13, 2012

intrinsic mode function (IMF); directly-mean EMD method; four-midPoint estimation method; Hilbert-Huang Transform; LOD (length of day)

On the basis of directly-mean EMD method, four-midpoint estimation formula according to all kinds of possible cases at both ends of signal was given by introducing the concepts of left and right undetermined middle points, left and right extended middle points. Four-midpoint estimation method for EMD, which can eliminate boundary effect completely, was established. This new method decompose signal into intrinsic mode functions and eliminate boundary effect more quickly and efficiently. The improvement on decomposition speed was brought by the reduction of the spline interpolation number. And the large fluctuation generated due to too sparse interpolation nodes was avoid, which made decomposition results more accurate. The new method was applied to the mode decomposition of length of day (LOD) data and a satisfactory decomposition result is obtained.

O29; P183.3+1

A

1000-3096(2013)04-0103-08

2012-01-13;

2013-02-04

國家自然科學基金資助項目(40876094)

張立振(1962-), 男, 山東高唐人, 副教授, 博士, 研究方向為數(shù)學, 物理海洋學, E-mail: goldfield@ouc.edu.cn