梁黎明， 王茂芝， 徐文皙， 譚夢(mèng)婷， 張明月， 王尚坤

(1. 成都理工大學(xué) 數(shù)理學(xué)院，成都 610059; 2. 成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室，成都 610059)

自從Huang等[1-2]于1998年提出經(jīng)驗(yàn)?zāi)B(tài)分解方法(empirical mode decomposition, EMD)以來，該方法被廣泛應(yīng)用于各學(xué)科領(lǐng)域非線性非平穩(wěn)信號(hào)處理。EMD將信號(hào)分解為頻率由高到低的本征模態(tài)函數(shù)(intrinsic mode function, IMF)，是一種完全自適應(yīng)的信號(hào)分析方法。然而，EMD在提取IMF時(shí)，依據(jù)信號(hào)局部極值點(diǎn)構(gòu)建其上、下包絡(luò)線，而信號(hào)端點(diǎn)通常不是極值點(diǎn)，進(jìn)而需要借助端點(diǎn)延拓才能使包絡(luò)線覆蓋信號(hào)并進(jìn)行后續(xù)操作。不良的端點(diǎn)延拓會(huì)導(dǎo)致包絡(luò)線兩端的誤差放大，進(jìn)而在所提取的IMF分量中 (特別是在IMF的兩端) 出現(xiàn)虛假成分，并且這些虛假成分逐漸向內(nèi)“傳播”，最終造成信號(hào)分解結(jié)果失真，產(chǎn)生所謂的端點(diǎn)效應(yīng)。當(dāng)輸入信號(hào)較短、頻率較小時(shí)，端點(diǎn)效應(yīng)更為明顯，進(jìn)而導(dǎo)致所提取的IMF失去物理意義。因此，抑制和消除端點(diǎn)效應(yīng)是EMD應(yīng)用的重要前提，受到學(xué)者的廣泛關(guān)注。

已有研究表明，抑制EMD端點(diǎn)效應(yīng)的端點(diǎn)延拓算法設(shè)計(jì)主要從增加信號(hào)兩端極值點(diǎn)來實(shí)現(xiàn)延拓包絡(luò)線的目的。如：Rilling等[3]提出的鏡像法以信號(hào)兩端的第一個(gè)極值點(diǎn)為邊界，把信號(hào)向外映射，通過獲取原信號(hào)的鏡像得到一個(gè)覆蓋原信號(hào)的包絡(luò)曲線；王紅軍等[4]提出改進(jìn)的多項(xiàng)式擬合法，利用端點(diǎn)處的三個(gè)極值點(diǎn)進(jìn)行多項(xiàng)式擬合，將計(jì)算值作為延拓點(diǎn)的近似取值；Huang等[5]提出極值法，以端點(diǎn)的一個(gè)特征波為依據(jù)，在兩端各延拓兩個(gè)極大值點(diǎn)和極小值點(diǎn)；Datig等[6]提出斜率法；Wu等[7]在Datig等的基礎(chǔ)上提出改進(jìn)的斜率方法(improved slope based method, ISBM)，根據(jù)信號(hào)端點(diǎn)相鄰極值點(diǎn)求得斜率，利用端點(diǎn)的第一個(gè)極值點(diǎn)和斜率確定延拓點(diǎn)所在的線性方程，以此確定延拓點(diǎn)的位置；Xu等[8]提出基于三次樣條的延拓算法(cubic spline based method, CSBM)，利用延拓后的包絡(luò)線必須覆蓋原包絡(luò)線的限制，確定信號(hào)兩端局部極值點(diǎn)與延拓點(diǎn)之間的關(guān)系，得到滿足條件的延拓點(diǎn)。此外，還有利用Lyapunov指數(shù)預(yù)測(cè)模型[9]和波形特征[10]來預(yù)測(cè)延拓點(diǎn)，以及采用其他類型的樣條函數(shù)[11-12]擬合包絡(luò)線進(jìn)行延拓。

CSBM由于較好地控制了延拓包絡(luò)線和原始包絡(luò)線在共同區(qū)域內(nèi)保持一致的特性，具有良好性能。但同時(shí)，CSBM在確定延拓點(diǎn)坐標(biāo)位置時(shí)依靠文中提出的兩個(gè)準(zhǔn)則及包絡(luò)線來確定，缺乏對(duì)信號(hào)趨勢(shì)的考量。另一方面，ISBM在延拓點(diǎn)位置的確定上，利用斜率來刻畫信號(hào)局部趨勢(shì)。本文提出一種基于斜率和三次樣條相結(jié)合的端點(diǎn)延拓方法(slope and cubic spline based method, S-CBM)，該方法綜合了ISBM和CSBM的優(yōu)勢(shì)，能有效抑制端點(diǎn)效應(yīng)對(duì)EMD分解結(jié)果的影響。

1 EMD算法

與傅里葉變換和小波變換等傳統(tǒng)方法相比，EMD是一種新型的時(shí)間序列分析方法。該方法依據(jù)數(shù)據(jù)自身時(shí)間尺度特征進(jìn)行信號(hào)分解，無須預(yù)設(shè)基函數(shù)，克服了傳統(tǒng)時(shí)頻分析方法需要確定基函數(shù)的局限性。理論上，EMD可以應(yīng)用于任何類型的信號(hào)處理，特別在處理非平穩(wěn)非線性數(shù)據(jù)上具有明顯優(yōu)勢(shì)。通過EMD篩選出來的IMF需滿足兩個(gè)準(zhǔn)則：一是上包絡(luò)線和下包絡(luò)線均值為零；二是極值點(diǎn)個(gè)數(shù)與過零點(diǎn)個(gè)數(shù)最多相差一個(gè)。

任給一個(gè)測(cè)試信號(hào)x(t)，IMF的提取過程如表1所示。

表1 EMD算法Tab.1 EMD algorithm

若篩選過程在n次迭代后停止，測(cè)試信號(hào)x(t)可以被表示為

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

由于EMD算法中篩選IMF的兩個(gè)條件過于苛刻，因此常用另一個(gè)有效判斷工具[13-14]標(biāo)準(zhǔn)偏差(standard deviation, SD)來替代。當(dāng)SD小于某個(gè)指定閾值(本文SD=0.05)時(shí)，終止篩選。其定義為

2 基于斜率和三次樣條的端點(diǎn)延拓算法

為描述方便，表2給出了有關(guān)符號(hào)定義和變量說明。

表2 符號(hào)定義和變量說明Tab.2 Symbol definition and variable description

由于延拓點(diǎn)的“好壞”將直接決定EMD分解的端點(diǎn)效應(yīng)程度，如何最大程度利用已有的信息合理確定延拓點(diǎn)位置是解決端點(diǎn)效應(yīng)的關(guān)鍵。本文的設(shè)計(jì)思路是：從信號(hào)及其包絡(luò)變化趨勢(shì)兩個(gè)角度共同來控制延拓點(diǎn)位置。其中信號(hào)變化趨勢(shì)度量是利用在信號(hào)兩端相鄰局部極值計(jì)算兩個(gè)斜率值，并以斜率變化率作為延拓點(diǎn)所在直線的斜率，進(jìn)而構(gòu)造延拓點(diǎn)的直線方程。另一方面，采用CSBM方法確定延拓點(diǎn)需滿足的另一個(gè)方程。最后通過聯(lián)立求解兩個(gè)方程來確定延拓點(diǎn)的坐標(biāo)位置。這種設(shè)計(jì)思路同時(shí)充分考慮了信號(hào)及其包絡(luò)的變化趨勢(shì)，進(jìn)而一定程度也克服了CSBM僅依靠包絡(luò)線來確定延拓點(diǎn)位置的局限性。

該算法需在信號(hào)兩端各延拓一個(gè)局部極大值點(diǎn)和一個(gè)局部極小值點(diǎn)。下面以信號(hào)左端延拓點(diǎn)為例，介紹并描述算法，算法示意圖如圖1所示。右側(cè)延拓點(diǎn)可在此基礎(chǔ)上推導(dǎo)得到。

圖1 基于S-CBM的信號(hào)左端延拓Fig.1 Extension of the signal’s left end of by S-CBM

表3 綜合斜率和CSBM的端點(diǎn)延拓方法Tab.3 Algorithm description of end condition method combined slope and CSBM

表3(續(xù))

圖1示例了按照表3算法流程得到阻尼正弦信號(hào)x(t)=sin(2πt)·e0.1t左端延拓極值點(diǎn)的計(jì)算情況(此時(shí)第一個(gè)極值點(diǎn)是極小值點(diǎn))。

3 試驗(yàn)結(jié)果與性能分析

本文設(shè)計(jì)了四個(gè)模擬信號(hào)和兩個(gè)實(shí)際工程應(yīng)用信號(hào)對(duì)算法性能進(jìn)行測(cè)試。四個(gè)模擬信號(hào)簡(jiǎn)述如下：第一個(gè)測(cè)試信號(hào)采樣頻率為200 Hz；第二個(gè)是由兩個(gè)調(diào)幅信號(hào)組成的采樣頻率為300 Hz的信號(hào)；第三個(gè)由高、中、低頻信號(hào)構(gòu)成，采樣頻率為400 Hz；第四個(gè)例子為一個(gè)非平穩(wěn)信號(hào)，采樣頻率為200 Hz。四個(gè)信號(hào)的數(shù)學(xué)表達(dá)式定義如下，信號(hào)波形如圖2所示。兩個(gè)實(shí)際工程應(yīng)用信號(hào)來自美國(guó)凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心的一段正常電機(jī)驅(qū)動(dòng)端加速度以及風(fēng)扇端加速度的數(shù)據(jù)[16]，分別記為x5(t)和x6(t)，如圖3所示。

圖2 模擬信號(hào)Fig.2 Analog signals

圖3 實(shí)際工程信號(hào)Fig.3 Practical engineering signals

3.1 算法性能度量指標(biāo)

本文考慮三個(gè)度量指標(biāo)對(duì)延拓效果進(jìn)行評(píng)估。

第一個(gè)度量指標(biāo)是由Huang等提出的正交指數(shù)(index of orthogonality, IO)，其定義為

式中：x(t)為原始信號(hào)；IMFi(t),IMFj(t)分別為原始信號(hào)x(t)的第i、第j個(gè)IMF分量，IO評(píng)估IMFs的正交性。理論上，IO越小，信號(hào)分解的效果越好。

第二個(gè)度量指標(biāo)是均方誤差(mean squared error, MSE)，該指標(biāo)對(duì)具有確定分量的輸入信號(hào)進(jìn)行評(píng)估，MSE定義為

第三個(gè)指標(biāo)定義為能量誤差θ。在EMD分解過程中，所有IMF的能量總和應(yīng)等于原信號(hào)的能量，但由于能量泄漏等原因，IMF的能量總和往往無法與原信號(hào)保持一致。因此可以通過比較EMD分解前后的能量，評(píng)估端點(diǎn)效應(yīng)的影響程度。能量誤差θ可描述為

式中：si(t)為原始信號(hào)或IMF分量；RMS為信號(hào)的有效值；m,n分別為IMF的總個(gè)數(shù)和樣本數(shù)；θ值越小，信號(hào)分解效果越好。

3.2 試驗(yàn)結(jié)果對(duì)比

本節(jié)僅以信號(hào)x1(t)為例，對(duì)比ISBM，CSBM和S-CBM三種不同延拓方法得到的極值點(diǎn)和包絡(luò)線，如圖4所示。

圖4 基于三種不同方法的信號(hào)x1(t)的延拓點(diǎn)、上下包絡(luò)線和平均包絡(luò)線Fig.4 Extended points, upper envelops, lower envelopes and mean envelopes of signal x1(t) based on three different methods

從圖4中可以看出，在該時(shí)間序列中，基于ISBM延拓的左端極大值點(diǎn)和右端極小值點(diǎn)沒有很好地反映信號(hào)變化趨勢(shì)。這是因?yàn)樵谇蠼庋油攸c(diǎn)所在直線的斜率時(shí)，ISBM沒有考慮斜率變化過大帶來的影響，導(dǎo)致延拓點(diǎn)出現(xiàn)較大誤差。CSBM在該信號(hào)中延拓效果較好，但從信號(hào)左端放大的圖像可以發(fā)現(xiàn)，CSBM延拓的局部極大值點(diǎn)位于信號(hào)第一個(gè)點(diǎn)(t1,x1)的下方，這雖然較好地反映了信號(hào)左端變化趨勢(shì)，但卻導(dǎo)致拓展的上包絡(luò)線在局部區(qū)域位于信號(hào)下方，也就是不能很好地規(guī)避包絡(luò)線過沖和欠沖[17]。而S-CBM由于綜合了ISBM和CSBM兩種方法的優(yōu)勢(shì)，延拓點(diǎn)則更好地反映了信號(hào)的變化趨勢(shì)。

3.3 性能對(duì)比及分析

本節(jié)展示了六組測(cè)試信號(hào)在三種端點(diǎn)延拓方法下的EMD分解效果以及評(píng)價(jià)指標(biāo)結(jié)果。為了更好地分析三種方法的性能，將三種端點(diǎn)延拓方法提取的IMF展現(xiàn)在同一張圖中，其中最后一個(gè)IMF分量為剩余分量，代表信號(hào)的平均趨勢(shì)。需要注意的是，由于EMD延拓算法的不同，其篩選得到的IMF數(shù)量可能不一樣。表4、表5和表6分別示例了三種方法的正交指數(shù)、均方誤差和能量誤差對(duì)比情況。

表4 三種不同延拓方法的IO結(jié)果對(duì)比Tab.4 IO results of the four signals by using EMD for three methods

表5 三種不同延拓方法下的MSE結(jié)果對(duì)比Tab.5 MSE results of the four signals by using EMD for three methods

表6 三種不同延拓方法下的θ結(jié)果對(duì)比Tab.6 θ results of the four signals by using EMD for three methods

從定性的角度，對(duì)比圖5～圖10中三種延拓算法的EMD分解結(jié)果，發(fā)現(xiàn)：各信號(hào)所提取的IMF分量在兩端端點(diǎn)處的波動(dòng)幅度越大，則端點(diǎn)效應(yīng)越明顯，對(duì)應(yīng)的端點(diǎn)延拓算法性能就越不足。下面以信號(hào)x1(t)為例，結(jié)合該信號(hào)的分解結(jié)果(見圖5)來展開說明。對(duì)比圖5中的IMF1，發(fā)現(xiàn)ISBM提取的IMF1在左側(cè)端點(diǎn)處波動(dòng)幅度明顯大于CSBM和S-CBM所提取的結(jié)果，這一現(xiàn)象在后續(xù)的IMF2，IMF3及IMF4中表現(xiàn)更為突出。也就是說，對(duì)于x1(t)，ISBM所提取的各IMF分量在兩端的波動(dòng)幅度明顯大于CSBM和S-CBM所提取結(jié)果的波動(dòng)幅度。相對(duì)而言，CSBM和S-CBM兩種方法提取的IMF波動(dòng)幅度較小。但由于S-CBM集成了ISBM和CSBM的優(yōu)點(diǎn)，相比較而言，S-CBM的波動(dòng)幅度要小于CSBM。其他信號(hào)提取結(jié)果可以類似分析。顯然，所提取的IMF波動(dòng)幅度越大，則端點(diǎn)效應(yīng)越明顯，所以，定性分析的結(jié)果與前面定量分析的結(jié)果是一致和吻合的。

圖5 基于三種方法的x1(t) EMD分解結(jié)果Fig.5 EMD results of x1(t) based on the three methods

圖6 基于三種方法的x2(t) EMD分解結(jié)果Fig.6 EMD results of x2(t) based on the three methods

圖7 基于三種方法的x3(t) EMD分解結(jié)果Fig.7 EMD results of x3(t) based on the three methods

圖8 基于三種方法的x4(t) EMD分解結(jié)果Fig.8 EMD results of x4(t) based on the three methods

圖9 基于三種方法的x5(t) EMD分解結(jié)果Fig.9 EMD results of x5(t) based on the three methods

圖10 基于三種方法的x6(t) EMD分解結(jié)果Fig.10 EMD results of x6(t) based on the three methods

總體而言，試驗(yàn)對(duì)比結(jié)果表明：本文所提出的方法在性能上整體優(yōu)于CSBM和ISBM。

4 結(jié) 論

本文針對(duì)CSBM延拓算法存在的缺陷，提出了基于斜率和三次樣條相結(jié)合的端點(diǎn)延拓方法。該方法根據(jù)信號(hào)端點(diǎn)的變化趨勢(shì)，并利用延拓前后包絡(luò)線在公共區(qū)域完全重合的特性唯一確定延拓點(diǎn)位置，克服了CSBM在設(shè)置延拓點(diǎn)坐標(biāo)時(shí)沒有充分考慮信號(hào)趨勢(shì)的不足，進(jìn)而使延拓點(diǎn)坐標(biāo)確定更合理。模擬和實(shí)際工程應(yīng)用測(cè)試信號(hào)以及IO，MSE和θ三個(gè)度量指標(biāo)的試驗(yàn)對(duì)比結(jié)果表明：本文提出的方法在抑制端點(diǎn)效應(yīng)上總體性能更好。這也是我們?cè)谇捌诠ぷ骰A(chǔ)上，對(duì)CSBM中延拓點(diǎn)坐標(biāo)選取的一種完善和改進(jìn)。需要說明的是，本文所提方法涉及的有關(guān)理論我們?cè)赬u等的研究中已展開了詳細(xì)論述。

針對(duì)本文工作，下一步的研究主題主要包括如下兩個(gè)方面：一是EMD提取IMF的迭代收斂理論證明，進(jìn)而為端點(diǎn)延拓算法設(shè)計(jì)提供理論指導(dǎo)；二是端點(diǎn)延拓算法設(shè)計(jì)的數(shù)學(xué)實(shí)質(zhì)揭示。