EMD自適應三角波匹配延拓算法

王錄雁,王 強,魯冬林,張梅軍,毛 迪

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

Huang等[1]提出了經驗模態分解(Empirical Mode Decomposition，EMD)技術，在處理非平穩、非線性信號方面具有獨特的自適應性，使其成為近年國內外學者的研究熱點。經過多年的研究，EMD技術的理論體系逐漸完善，并已在眾多領域得到了實踐應用，但其端點效應問題一直制約著該方法的推廣應用。解決EMD端點效應主要有兩個途徑：數據延拓法和改進插值法。其中數據延拓法具有理論簡單、效果明顯的特點，因而成為主要的研究方向。文獻[2]對近年出現的典型的數據延拓算法進行了分類綜述，分析了各種方法的特點及不足。其中波形匹配方法[3-8]，立足于近端點處波形與內部波形的相似性實現數據延拓，延拓數據不但考慮了近端點處數據的局部變化趨勢，同時保留了波形的原始特征，具有其他方法所不可比擬的優勢，因而成為遏制EMD端點效應的一個重要研究方向。

波形匹配法的關鍵步驟之一是波形匹配程度的計算。文獻[3，8]及文獻[4-7]分別采用了三角波和基于距離的波形匹配算法，在實際應用中尚且存在一些不足，本文在此基礎上，提出一種自適應的三角波形匹配延拓算法。

1 EMD方法及其端點效應

1.1 基本原理

EMD方法是將非線性、非平穩的信號分解為有限個具有單一瞬時頻率的固有模態函數(Intrinsic mode function，IMF)，步驟如下：

(1)確定信號的極值點；

(2)用三次樣條曲線對所有極值點進行插值，生成上、下包絡線；

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

h1=x(t)-m1

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

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

(1)

1.2 端點效應問題

在采用三次樣條曲線對極值點進行插值時，由于數據兩端端點處數值不能同時既為極大值又為極小值，于是在擬合包絡線時，會在端點處失去約束，導致包絡線在端點處發生大幅擺動。隨著分解的進行，誤差就會由端點處逐漸向內傳播，嚴重的情況下會使分解的數據喪失物理意義。如圖1所示，在沒有采取抑制措施的情況下，包絡線在端點處出現了大幅擺動，致使分解結果發生嚴重偏離，只分解出了兩個有意義的頻率成分便終止了循環篩分；采用鏡像極值延拓法有效遏制了端點效應，獲得了三個有效的頻率成分。

圖1 EMD端點效應問題及端點延拓

2 波形匹配方法

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

2.1 正弦波匹配延拓

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

圖2 正弦波匹配延拓

該方法將近端點的波形近似看作是某正弦波的一部分，通過該段波形構造出一段近似匹配的波形加以延拓，實現了波形的平順延拓。但該方法僅僅依靠近端點處的數據特征實施延拓，割裂了近端點處信號與內部信號的聯系，無法真實的反映信號的變化趨勢及內部特征。

2.2 基于距離的波形匹配法

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

(2)

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

圖3 基于距離的波形匹配過程

該方法用原始波形內部的一段子波實現延拓，不但充分考慮了信號近端點處的數據變化趨勢信息，同時兼顧了端點處信號與內部信號的聯系，在獲得波形平順延拓的同時，維持了信號的原始特征。但是，由于該方法建立在子波數據逐點比對的基礎上，因此要求子波間的數據長度嚴格一致，導致在應用中存在一些問題。

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

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

圖4 等長度子波截取時的錯誤

(2)端點位于極值附近時預測準確性差

按照距離匹配算法，當端點數據位于極值附近時，由于可用數據相對較少，將難以準確預測數據的走勢。以常用仿真信號為例，表達式為：

t∈[0,500]

(3)

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

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

圖5 端點位于極值附近時的錯誤匹配

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

(4)

或：

(5)

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

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

圖6 三角波匹配法

3 自適應三角波匹配延拓算法

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

(1)提取極值點

設，信號x(t)分別有M個極大值和N個極小值；x(i)為信號第i點的數值，m(i)，x(m(i))，分別為第i個極大值的序列位置和數值；n(i)，x(n(i))分別為第i個極小值的序列位置和數值。

(2)截取模板子波

如圖7，截取從端點至第二個極值點間的數據構成模板子波，設為S1。

(3)局部尋優

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

圖7 自適應三角波匹配算法

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

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

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

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

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

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

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

其中：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為各匹配度的權重向量，默認為等權重。

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

(4)全局尋優

記第i個極大值對應的由t(i)截取的最佳三角子波記為Mi。取r=find(Mi=maxMi)，則r為最佳匹配子波對應的極值點位置。按照2.2節所述方法可實現左右兩端點的匹配延拓。

(5)細節處理

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

4 算例驗證分析

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

4.1 仿真信號分析

采用本方法對式3中信號進行匹配計算，由表2可看出(精確到10-5)，該方法實現了準確的波形匹配。由圖8～9可看出，正弦波匹配法雖然一定程度上抑制了端點效應，但是沒能正確反映信號變化趨勢，分解過程出現較大波動；采用距離函數法匹配時，出現了錯誤匹配現象，致使各分量在端點處出現了變形；自適應三角波匹配延拓與支持向量機預測延拓法，均實現了較為準確延拓。其中w取默認等權重，λ取0.9；SVR延拓點數為100個，支持向量個數為200。

表2 各子波與S1的匹配度

圖8 各種延拓方法延拓結果對比

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

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

4.2 實驗信號分析

圖10為某軸承振動信號截取的片段，從小窗口可看出，a,b兩段波形具有較高的相似性，若能實現準確的匹配延拓，將有效抑制EMD的端點效應問題，提高分解精度。

表3 仿真信號EMD分解指標

圖10 某軸承振動信號片段

由圖11～12可知，采用距離函數匹配延拓法，在EMD處理過程中，在左端點處出現了較明顯的波動；自適應三角波匹配延拓法，較好的抑制了EMD端點效應。表4顯示，自適應三角波匹配延拓法，分解精度較高，且不需要大量樣本學習計算，運行速度比SVR方法有大幅提升。其中w取默認等權重，λ取0.8；SVR延拓點數為50個，支持向量個數為100。

圖11 各方法延拓結果對比

表4 實測信號的EMD分解指標

圖12 兩種延拓方法EMD結果的前5項

5 結 論

本文在總結現有波形匹配延拓方法特點與不足的基礎上，提出了一種自適應的三角波形匹配延拓方法。該方法有以下有點：

(1)采用三角波匹配法實現延拓，規避了子波等長度匹配所導致的問題；

(2)通過改進三角波匹配算法，自適應地截取最佳匹配三角子波，改善了固定比例截取子波時產生的問題，保證了子波截取的有效性與合理性。

(3)通過局部尋優與全局尋優相結合的方法實現匹配子波的截取與定位，提高了波形匹配延拓的準確性。

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

仿真及實驗數據分析表明，該方法提高了波形匹配的準確性，可有效遏制EMD端點效應問題，顯著提高分解精度。文中僅對匹配度閾值λ的取值方法做了定性分析，尚需進一步研究其取值及調整方法，另外權重向量w的取值及調整方法同樣需要進一步深入研究。

參 考 文 獻

[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]賀 智，王 強，沈 毅，等．希爾伯特—黃變換端點效應抑制算法綜述[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]蓋 強，馬孝江，張海勇，等．一種消除局域波法中邊界效應的新方法[J]．大連理工大學學報，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]邵晨曦，王 劍，范金鋒，等．一種自適應的EMD端點延拓方法[J]．電子學報，2007，35(10)：1944-1947．

SHAO Chen-xi，WANG Jian，FAN 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]．機械工程學報，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端點效應的新方法[J]．哈爾濱理工大學學報，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]楊 慶，陳桂明，薛冬林.基于最小平方距離相關的 EMD改進算法及應用[J]．振動與沖擊，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]鐘佑明，趙 強，周建庭.一種基于本征波匹配的EMD邊界處理方法[J]．振動與沖擊，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]張 亢，程軍圣，楊 宇.基于自適應波形匹配延拓的局部均值分解端點效應處理方法[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變換中的一種端點延拓方法[J]．振動與沖擊，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，Flandrin 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．