EMD自適應(yīng)三角波匹配延拓算法

王錄雁,王 強(qiáng),魯冬林,張梅軍,毛 迪

(解放軍理工大學(xué)，南京 210007)

Huang等[1]提出了經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)技術(shù)，在處理非平穩(wěn)、非線性信號(hào)方面具有獨(dú)特的自適應(yīng)性，使其成為近年國(guó)內(nèi)外學(xué)者的研究熱點(diǎn)。經(jīng)過(guò)多年的研究，EMD技術(shù)的理論體系逐漸完善，并已在眾多領(lǐng)域得到了實(shí)踐應(yīng)用，但其端點(diǎn)效應(yīng)問(wèn)題一直制約著該方法的推廣應(yīng)用。解決EMD端點(diǎn)效應(yīng)主要有兩個(gè)途徑：數(shù)據(jù)延拓法和改進(jìn)插值法。其中數(shù)據(jù)延拓法具有理論簡(jiǎn)單、效果明顯的特點(diǎn)，因而成為主要的研究方向。文獻(xiàn)[2]對(duì)近年出現(xiàn)的典型的數(shù)據(jù)延拓算法進(jìn)行了分類(lèi)綜述，分析了各種方法的特點(diǎn)及不足。其中波形匹配方法[3-8]，立足于近端點(diǎn)處波形與內(nèi)部波形的相似性實(shí)現(xiàn)數(shù)據(jù)延拓，延拓?cái)?shù)據(jù)不但考慮了近端點(diǎn)處數(shù)據(jù)的局部變化趨勢(shì)，同時(shí)保留了波形的原始特征，具有其他方法所不可比擬的優(yōu)勢(shì)，因而成為遏制EMD端點(diǎn)效應(yīng)的一個(gè)重要研究方向。

波形匹配法的關(guān)鍵步驟之一是波形匹配程度的計(jì)算。文獻(xiàn)[3，8]及文獻(xiàn)[4-7]分別采用了三角波和基于距離的波形匹配算法，在實(shí)際應(yīng)用中尚且存在一些不足，本文在此基礎(chǔ)上，提出一種自適應(yīng)的三角波形匹配延拓算法。

1 EMD方法及其端點(diǎn)效應(yīng)

1.1 基本原理

EMD方法是將非線性、非平穩(wěn)的信號(hào)分解為有限個(gè)具有單一瞬時(shí)頻率的固有模態(tài)函數(shù)(Intrinsic mode function，IMF)，步驟如下：

(1)確定信號(hào)的極值點(diǎn)；

(2)用三次樣條曲線對(duì)所有極值點(diǎn)進(jìn)行插值，生成上、下包絡(luò)線；

(3)計(jì)算上下包絡(luò)線均值m1，求出h1

h1=x(t)-m1

(4)判斷h1是否符合IMF判據(jù)，如果h1符合IMF判據(jù)[1]，則h1即為一個(gè)IMF；若h1不符合IMF判據(jù)，則重復(fù)執(zhí)行步驟(1)～(3)，直到h1k=h1(k-1)-m1k符合判據(jù)。

(5)分離IMF后重復(fù)執(zhí)行步驟(1)～(4)，直到分離出所有的IMF(ci)直到殘留項(xiàng)rn不再滿足篩分條件[1]為止，則有：

(1)

1.2 端點(diǎn)效應(yīng)問(wèn)題

在采用三次樣條曲線對(duì)極值點(diǎn)進(jìn)行插值時(shí)，由于數(shù)據(jù)兩端端點(diǎn)處數(shù)值不能同時(shí)既為極大值又為極小值，于是在擬合包絡(luò)線時(shí)，會(huì)在端點(diǎn)處失去約束，導(dǎo)致包絡(luò)線在端點(diǎn)處發(fā)生大幅擺動(dòng)。隨著分解的進(jìn)行，誤差就會(huì)由端點(diǎn)處逐漸向內(nèi)傳播，嚴(yán)重的情況下會(huì)使分解的數(shù)據(jù)喪失物理意義。如圖1所示，在沒(méi)有采取抑制措施的情況下，包絡(luò)線在端點(diǎn)處出現(xiàn)了大幅擺動(dòng)，致使分解結(jié)果發(fā)生嚴(yán)重偏離，只分解出了兩個(gè)有意義的頻率成分便終止了循環(huán)篩分；采用鏡像極值延拓法有效遏制了端點(diǎn)效應(yīng)，獲得了三個(gè)有效的頻率成分。

圖1 EMD端點(diǎn)效應(yīng)問(wèn)題及端點(diǎn)延拓

2 波形匹配方法

近年來(lái)，波形匹配延拓法的研究主要有正弦波匹配法[5]、基于距離的波形匹配法[4,6-7]以及相似三角波的匹配法[3,8-9]。

2.1 正弦波匹配延拓

文獻(xiàn)[5]提出一種使用正弦函數(shù)進(jìn)行端點(diǎn)延拓的方法。其思想是根據(jù)端點(diǎn)至首個(gè)極值之間的波形特征，構(gòu)造出一段幅值、頻率、相位與其近似匹配的正弦波實(shí)現(xiàn)數(shù)據(jù)的外延(圖2)。

圖2 正弦波匹配延拓

該方法將近端點(diǎn)的波形近似看作是某正弦波的一部分，通過(guò)該段波形構(gòu)造出一段近似匹配的波形加以延拓，實(shí)現(xiàn)了波形的平順延拓。但該方法僅僅依靠近端點(diǎn)處的數(shù)據(jù)特征實(shí)施延拓，割裂了近端點(diǎn)處信號(hào)與內(nèi)部信號(hào)的聯(lián)系，無(wú)法真實(shí)的反映信號(hào)的變化趨勢(shì)及內(nèi)部特征。

2.2 基于距離的波形匹配法

邵晨曦等[4]將基于距離函數(shù)的波形匹配法應(yīng)用于EMD端點(diǎn)延拓，取得了較好的效果。隨后，王婷等[6-7]，也分別運(yùn)用相似的算法實(shí)現(xiàn)了EMD端點(diǎn)延拓。這類(lèi)的延拓的思想是，在波形內(nèi)部尋找一段與近端點(diǎn)處子波(S1)相似的匹配波形，截取匹配波形(Sj)前一段子波(S’)對(duì)端點(diǎn)進(jìn)行延拓。如圖3所示，兩子波的匹配距離為：

(2)

式中，s1(i)、sj(i)分別為兩子波第i個(gè)點(diǎn)的值，N為子波數(shù)據(jù)長(zhǎng)度。取d(S1，Sopt)=min(d(S1,Sj))，截取Sopt前包含兩個(gè)極值(或更長(zhǎng))的子波S’作為匹配波形平移至端點(diǎn)處實(shí)現(xiàn)延拓。

圖3 基于距離的波形匹配過(guò)程

該方法用原始波形內(nèi)部的一段子波實(shí)現(xiàn)延拓，不但充分考慮了信號(hào)近端點(diǎn)處的數(shù)據(jù)變化趨勢(shì)信息，同時(shí)兼顧了端點(diǎn)處信號(hào)與內(nèi)部信號(hào)的聯(lián)系，在獲得波形平順延拓的同時(shí)，維持了信號(hào)的原始特征。但是，由于該方法建立在子波數(shù)據(jù)逐點(diǎn)比對(duì)的基礎(chǔ)上，因此要求子波間的數(shù)據(jù)長(zhǎng)度嚴(yán)格一致，導(dǎo)致在應(yīng)用中存在一些問(wèn)題。

(1)致使無(wú)意義的子波截取

當(dāng)波形的時(shí)間尺度變化較大或子波較長(zhǎng)時(shí)，難以保證所截取的子波在參考極值兩側(cè)具有相同的增減特性。如圖4，按照等長(zhǎng)度截取子波S1及S8，其中S8包含了兩個(gè)極大值點(diǎn)，使得匹配失去實(shí)際意義。

圖4 等長(zhǎng)度子波截取時(shí)的錯(cuò)誤

(2)端點(diǎn)位于極值附近時(shí)預(yù)測(cè)準(zhǔn)確性差

按照距離匹配算法，當(dāng)端點(diǎn)數(shù)據(jù)位于極值附近時(shí)，由于可用數(shù)據(jù)相對(duì)較少，將難以準(zhǔn)確預(yù)測(cè)數(shù)據(jù)的走勢(shì)。以常用仿真信號(hào)為例，表達(dá)式為：

t∈[0,500]

(3)

由圖5所示，端點(diǎn)距極大值只有一個(gè)點(diǎn)距。整體上看，與S1對(duì)應(yīng)匹配的是S9和S17，但是計(jì)算結(jié)果(表1，精確到10-5)顯示S13與S1匹配距離最小。由小窗口可看出，匹配錯(cuò)誤是由于匹配數(shù)據(jù)太少，子波S13在形狀上和S1更加接近。

表1 各子波與S1的匹配距離分布

圖5 端點(diǎn)位于極值附近時(shí)的錯(cuò)誤匹配

三角波匹配法，是指以端點(diǎn)與鄰近兩個(gè)極值點(diǎn)所構(gòu)成的三角形波形為匹配子波，構(gòu)造三角波形的匹配算法以實(shí)現(xiàn)波形的延拓。其突出的特點(diǎn)是不要求子波長(zhǎng)度相等(圖6)。文獻(xiàn)[3,8]分別研究了三角波相似匹配的問(wèn)題，但并未提及三角波起始點(diǎn)的確定方法。文獻(xiàn)[9]將自適應(yīng)三角波匹配應(yīng)用于局部均值分解的端點(diǎn)延拓，采用等比關(guān)系確定三角子波起始點(diǎn)t(i)的方法：

(4)

或：

(5)

其中：t(mi)是第i個(gè)極大值的時(shí)刻位置，t(ni)是第i個(gè)極小值的時(shí)刻位置。

由式(4)可知，t(i)的確定建立在A、B兩個(gè)三角形相似的基礎(chǔ)上。而在實(shí)際應(yīng)用中，難以保證t(ni-1)于t(ni)之間一定存在一點(diǎn)，使得A、B兩三角形相似。即，根據(jù)式5計(jì)算所得t(i)不一定位于t(ni-1)和t(ni)之間，同樣使得截取的三角子波無(wú)意義(如圖4)，更無(wú)法保證所截取的三角子波具有最佳的匹配度。

圖6 三角波匹配法

3 自適應(yīng)三角波匹配延拓算法

綜上所述，正弦波延拓割裂了與內(nèi)部波形的聯(lián)系；距離函數(shù)匹配法要求子波長(zhǎng)度相等，易因數(shù)據(jù)較少而無(wú)法真實(shí)反映波形宏觀走勢(shì)，而增加子波長(zhǎng)度時(shí)易導(dǎo)致無(wú)意義的匹配計(jì)算；三角波匹配法，不要求子波長(zhǎng)度相等，但固定比例法截取子波，難以保證其有效性和匹配的準(zhǔn)確性。本文在此基礎(chǔ)上提出一種自適應(yīng)的三角波匹配算法。以首個(gè)極值為極大值情況下的左端點(diǎn)延拓為例說(shuō)明算法步驟。

(1)提取極值點(diǎn)

設(shè)，信號(hào)x(t)分別有M個(gè)極大值和N個(gè)極小值；x(i)為信號(hào)第i點(diǎn)的數(shù)值，m(i)，x(m(i))，分別為第i個(gè)極大值的序列位置和數(shù)值；n(i)，x(n(i))分別為第i個(gè)極小值的序列位置和數(shù)值。

(2)截取模板子波

如圖7，截取從端點(diǎn)至第二個(gè)極值點(diǎn)間的數(shù)據(jù)構(gòu)成模板子波，設(shè)為S1。

(3)局部尋優(yōu)

如圖7，須在n(i-1)至m(i)尋找一個(gè)時(shí)刻點(diǎn)t(i)=j作為三角波的起點(diǎn)。對(duì)于每一個(gè)固定的i∈[2,n]，取j∈[n(i-1)m(i)]，則以x(j)、x(m(i))、x(n(i))為頂點(diǎn)共可構(gòu)成m(i)-n(i-1)個(gè)三角波，設(shè)為Si(j)。求取一點(diǎn)t(i)=j，使得Si(j)與S1的匹配度最大。

圖7 自適應(yīng)三角波匹配算法

為了在全面反映子波整體形狀特征的同時(shí)突出延拓的平順性，文章將單邊波斜率及端點(diǎn)幅值差引入到波形匹配計(jì)算中。

令：Xi(j)為第i個(gè)極大值時(shí)，以j時(shí)刻為起點(diǎn)，左邊波的斜率匹配度：

i∈[2,M]，j∈[n(i-1),m(i)]，下同。

令：Yi(j)為左邊波的幅值匹配度：

令：Oi(j)為右邊波幅值匹配度：

令：Pi(j)為右邊波時(shí)刻匹配度：

令：Qi(j)為端點(diǎn)幅值差匹配度：

其中：a=[x(1)-x(n(1))]×[x(t(i))-x(n(i))]

則，令三角波Si(j)與S1的匹配度Mi(j)為：

Mi(j)=w(Xi(j),Yi(j),Oi(j),Pi(j),Qi(j))

(6)

其中w為各匹配度的權(quán)重向量，默認(rèn)為等權(quán)重。

取t(i)=j=find[Mi(j)=maxMi(j)]，t(i)即為第i段三角子波的起始點(diǎn)。

(4)全局尋優(yōu)

記第i個(gè)極大值對(duì)應(yīng)的由t(i)截取的最佳三角子波記為Mi。取r=find(Mi=maxMi)，則r為最佳匹配子波對(duì)應(yīng)的極值點(diǎn)位置。按照2.2節(jié)所述方法可實(shí)現(xiàn)左右兩端點(diǎn)的匹配延拓。

(5)細(xì)節(jié)處理

在實(shí)際應(yīng)用中，由于波形長(zhǎng)度限制，或難以找到匹配度較高的子波，設(shè)置匹配閾值λ∈(0,1)，當(dāng)M(i)<λ時(shí)，可采用相似極值延拓[10]或鏡像延拓法[11]處理，文中采用鏡像延拓法。一般地，信號(hào)越長(zhǎng)、規(guī)律越強(qiáng)，λ取值應(yīng)越大。

4 算例驗(yàn)證分析

為驗(yàn)證自適應(yīng)三角波形匹配延拓的有效性，分別采用正弦波延拓法、基于距離函數(shù)的波形匹配延拓法以及支持向量回歸機(jī)(Support Vector Regression，SVR)預(yù)測(cè)延拓法[12]對(duì)仿真信號(hào)與實(shí)測(cè)信號(hào)進(jìn)行對(duì)比驗(yàn)證。

4.1 仿真信號(hào)分析

采用本方法對(duì)式3中信號(hào)進(jìn)行匹配計(jì)算，由表2可看出(精確到10-5)，該方法實(shí)現(xiàn)了準(zhǔn)確的波形匹配。由圖8～9可看出，正弦波匹配法雖然一定程度上抑制了端點(diǎn)效應(yīng)，但是沒(méi)能正確反映信號(hào)變化趨勢(shì)，分解過(guò)程出現(xiàn)較大波動(dòng)；采用距離函數(shù)法匹配時(shí)，出現(xiàn)了錯(cuò)誤匹配現(xiàn)象，致使各分量在端點(diǎn)處出現(xiàn)了變形；自適應(yīng)三角波匹配延拓與支持向量機(jī)預(yù)測(cè)延拓法，均實(shí)現(xiàn)了較為準(zhǔn)確延拓。其中w取默認(rèn)等權(quán)重，λ取0.9；SVR延拓點(diǎn)數(shù)為100個(gè)，支持向量個(gè)數(shù)為200。

表2 各子波與S1的匹配度

圖8 各種延拓方法延拓結(jié)果對(duì)比

圖9 采用各種延拓方法處理的EMD結(jié)果

表3顯示(精確到10-6)，采用自適應(yīng)三角波匹配延拓法，使EMD處理精度得到顯著提高，對(duì)于這種規(guī)律信號(hào)，這些誤差基本上是由于包絡(luò)擬合誤差所產(chǎn)生。

4.2 實(shí)驗(yàn)信號(hào)分析

圖10為某軸承振動(dòng)信號(hào)截取的片段，從小窗口可看出，a,b兩段波形具有較高的相似性，若能實(shí)現(xiàn)準(zhǔn)確的匹配延拓，將有效抑制EMD的端點(diǎn)效應(yīng)問(wèn)題，提高分解精度。

表3 仿真信號(hào)EMD分解指標(biāo)

圖10 某軸承振動(dòng)信號(hào)片段

由圖11～12可知，采用距離函數(shù)匹配延拓法，在EMD處理過(guò)程中，在左端點(diǎn)處出現(xiàn)了較明顯的波動(dòng)；自適應(yīng)三角波匹配延拓法，較好的抑制了EMD端點(diǎn)效應(yīng)。表4顯示，自適應(yīng)三角波匹配延拓法，分解精度較高，且不需要大量樣本學(xué)習(xí)計(jì)算，運(yùn)行速度比SVR方法有大幅提升。其中w取默認(rèn)等權(quán)重，λ取0.8；SVR延拓點(diǎn)數(shù)為50個(gè)，支持向量個(gè)數(shù)為100。

圖11 各方法延拓結(jié)果對(duì)比

表4 實(shí)測(cè)信號(hào)的EMD分解指標(biāo)

圖12 兩種延拓方法EMD結(jié)果的前5項(xiàng)

5 結(jié) 論

本文在總結(jié)現(xiàn)有波形匹配延拓方法特點(diǎn)與不足的基礎(chǔ)上，提出了一種自適應(yīng)的三角波形匹配延拓方法。該方法有以下有點(diǎn)：

(1)采用三角波匹配法實(shí)現(xiàn)延拓，規(guī)避了子波等長(zhǎng)度匹配所導(dǎo)致的問(wèn)題；

(2)通過(guò)改進(jìn)三角波匹配算法，自適應(yīng)地截取最佳匹配三角子波，改善了固定比例截取子波時(shí)產(chǎn)生的問(wèn)題，保證了子波截取的有效性與合理性。

(3)通過(guò)局部尋優(yōu)與全局尋優(yōu)相結(jié)合的方法實(shí)現(xiàn)匹配子波的截取與定位，提高了波形匹配延拓的準(zhǔn)確性。

(4)匹配度算法中加入了單邊斜率匹配度和端點(diǎn)幅值差匹配度，突出了延拓的平順性，強(qiáng)化了近端點(diǎn)波形與內(nèi)部波形間的聯(lián)系。

仿真及實(shí)驗(yàn)數(shù)據(jù)分析表明，該方法提高了波形匹配的準(zhǔn)確性，可有效遏制EMD端點(diǎn)效應(yīng)問(wèn)題，顯著提高分解精度。文中僅對(duì)匹配度閾值λ的取值方法做了定性分析，尚需進(jìn)一步研究其取值及調(diào)整方法，另外權(quán)重向量w的取值及調(diào)整方法同樣需要進(jìn)一步深入研究。

參 考 文 獻(xiàn)

[1]Huang N E，Shen Z，Long S R．The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[C]．Proceedings of the Royal Society of London Series．1998，454：903-995．

[2]賀 智，王 強(qiáng)，沈 毅，等．希爾伯特—黃變換端點(diǎn)效應(yīng)抑制算法綜述[J]．軟件，2011，32(10)：1-7．

HE Zhi，WANG Qiang，SHEN Yi，et al．Survey on end effects Mitigation of Hilbert-Huang Transform [J]．Software，2011，32(10)：1-7．

[3]蓋 強(qiáng)，馬孝江，張海勇，等．一種消除局域波法中邊界效應(yīng)的新方法[J]．大連理工大學(xué)學(xué)報(bào)，2002，42(1)：115-117．

GAI Qiang， MA Xiao-jiang， ZHANG Hai-yong， et al.A new method to eliminate the boundary effect in local wave method [J]．Journal of Aalian University of Technology， 2002，42(1)：115-117.

[4]邵晨曦，王 劍，范金鋒，等．一種自適應(yīng)的EMD端點(diǎn)延拓方法[J]．電子學(xué)報(bào)，2007，35(10)：1944-1947．

SHAO Chen-xi，WANG Jian，F(xiàn)AN Jin-feng，et al．A self-adaptive method dealing with the end issue of EMD [J]．Acta Electronica Sinca，2007，35(10)：1944-1947．

[5]黃先祥，李勝朝，謝 建．新型經(jīng)驗(yàn)?zāi)Ｊ椒纸舛它c(diǎn)效應(yīng)消除方法 [J]．機(jī)械工程學(xué)報(bào)，2008，44(9)：1-5．

HUANG Xian-xiang，LI Sheng-chao，XIE Jian．New approach to dealing with the end effect of empirical model decomposition [J]．Chinese Journal of Mechanical Engineering，2008，44(9)：1-5．

[6]王 婷，楊莘元，李冰冰.一種改善EMD端點(diǎn)效應(yīng)的新方法[J]．哈爾濱理工大學(xué)學(xué)報(bào)，2009，14(5)：24-26．

WANG Ting，YANG Xin-yuan，LI Bing-bing．A new method for end effect of EMD [J]．Journal of Harbn University of Science and Technology，2009，14(5)：24-26．

[7]楊 慶，陳桂明，薛冬林.基于最小平方距離相關(guān)的 EMD改進(jìn)算法及應(yīng)用[J]．振動(dòng)與沖擊，2011，30(6)：62-66．

YANG Qing，CHEN Gui-ming，XUE Dong-lin．Improved method for empirical mode decomposition based on SSDA and its application[J]．Journal of Vibration and Shock,2011，30(6)：62-66．

[8]鐘佑明，趙 強(qiáng)，周建庭.一種基于本征波匹配的EMD邊界處理方法[J]．振動(dòng)與沖擊，2012，31(1)：16-19．

ZHONG You-ming，ZHAO Qiang，ZHOU Jian-ting．A method for EMD boundary processing based on intrinsic wave matching [J]．Journal of Vibration and Shock，2012，31(1)：16-19．

[9]張 亢，程軍圣，楊 宇.基于自適應(yīng)波形匹配延拓的局部均值分解端點(diǎn)效應(yīng)處理方法[J]．中國(guó)機(jī)械工程，2010，21(4)：457-462．

ZHANG Kang，CHENG Jun-sheng，YANG Yu．Processing method for end effects of local mean decomposition based on self-adaptive waveform matching extending [J]．China Mechanical Engineering，2010，21(4)：457-462．

[10]沈 路，周曉軍，張志剛，等．Hilbert-Huang變換中的一種端點(diǎn)延拓方法[J]．振動(dòng)與沖擊，2009，28(8)：168-171．

SHEN Lu,ZHOU Xiao-jun，ZHANG Zhi-gang．A endpoint extension method for Hilbert-Huang Transform [J]．Journal of Vibration and Shock，2009，28(8)：168-171．

[11]Rilling G，F(xiàn)landrin P， Gon?alvès P．On empirical mode decomposition and its algorithms[J]．IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03，2003，3：8-11.

[12]Cheng J S，Yu D J，Yang Y．Application of support vector regression machines to the processing of end effects of Hilbert-Huang [J]．Mechanical Systems and Signal Processing，2007，21：1197-1211．