一種適用于流數據分析的快速EMD算法

胡利萍，宋恩亮，李寶清，袁曉兵

（中國科學院 上海微系統與信息技術研究所 無線傳感網實驗室，上海 200050）

希爾伯特-黃變換（HHT）是Huang［1］提出的一種時間序列分析新方法，包括前端處理EMD和后端Hilbert譜分析。無須先驗信息，EMD方法可自適應地將信號分解為一組IMF分量。由于HHT具備較強的非線性非平穩信號分析能力，且應用方式靈活，既可單獨采用EMD進行處理，也可完整采用HHT提取Hilbert譜或邊際譜，因此HHT已在眾多工程領域得到廣泛應用。例如基于 EMD 的信號去噪［2］，電站監測［3］，圖像處理［4］，戰場偵察傳感網［5-6］等。EMD 是 HHT 的核心處理技術，但是EMD存在計算量大、分解速度慢、實時性差的缺陷。經典EMD采用分段處理法，將長序列分段然后逐段進行分解，顯然IMF在分段處是不連貫的。另外，由于端點效應的影響，有效長度因分段變短。另一種方法是等所有數據都采集完畢后，作一次性分解，但是在長序列的情況下，該做法影響實時性，也對硬件的存儲空間和計算能力提出了更高要求。在HHT應用日益普遍的形勢下，開發一種實時性好、計算量小的快速EMD算法具有重要意義。在快速算法研究方面，目前有以下三種主要思路：

（1）采用更加簡便的擬合方法，適當放寬終止準則［7］；

（2）只對有效數據進行擬合、和進行終止準則判定［8］；

（3）建立濾波器組提取各個分量［9-10］。

但是這些方法未能擺脫經典EMD的分段處理法，它們的共同特點是：

（1）分段處理，劃分處引發端點效應；

（2）若進行重疊劃分，由于各段去均值的次數不同，劃分處的IMF不連續；

（3）IMF從低階到高階依次產生，當低階IMF完全確定后才能開始高一階IMF的計算；

（4）全局終止準則，不能體現局部偏離情況，若不滿足終止準則，須重新進行全局去均值。

本文提出一種適用于流數據分析的快速EMD算法。該算法建立在“局域”和“局域終止準則”的概念基礎之上，采用二次B樣條函數直接進行均值線擬合，采用局域終止準則判斷原型模態函數的對稱性，在局域內進行分解，同時產生各階IMF。若采用該算法，無論序列長度，端點效應僅存在于流數據的起始端和結束端，相比分段處理法不僅有效數據長度增加，而且各階IMF是連續的。

1 經典EMD方法基本原理

EMD方法能把非平穩、非線性時間序列分解成一組由序列本身決定的數據序列集，即本征模態函數（IMF）。IMF須滿足2個條件：

（1）極值點和過零點數目必須相等或至多相差一點；

（2）對稱性，由局部極大點構成的包絡線和局部極小點構成的包絡線的平均值為零。

其本質是通過分辨特征時間尺度獲得本征振動模式，然后由本征振動模式來分解時間序列的一種算法。對時間序列x（t）進行分解的基本步驟如下：

（1）初始化：r0（t）=x（t），i=1；

（3）ri（t）=ri-1（t）-imfi（t）；

（4）若ri（t）滿足繼續分解的條件，則i=i+1，轉到（2）；否則分解結束，ri（t）是殘余分量。

最終分解結果為：

即原始序列可分解為n個IMF和一個殘余項之和。

2 快速EMD算法基本原理

從上述基本原理，經典算法的復雜度主要來自于不斷地進行三次樣條插值。SD準則通過判斷原型模態函數的對稱性以判定IMF。但是通常原型模態函數僅局部出現非對稱［11］（如圖1），因此每次插值和去均值操作僅部分為有效運算。

以“局域”替代經典EMD的分段處理技術，將“局域”作為當前運算范圍，一個“局域”的長度遠小于一個分段的長度，從而減少無效運算。采用B樣條函數以滑動平均方式擬合均值，在“局域”內使用SD終止準則，即局域SD準則。

圖1 局部非對稱的原型模態函數Fig.1 A proto-mode function with local asymmetry

設U是m個非遞減數的集合，u1≤u2≤u3≤…≤um，ui稱為節點，集合U稱為節點向量，半開區間［ui，ui+1）是第i個節點區間。B樣條函數的Cox-de Boor遞推定義為：

其中，Ni，p（u）是由節點ui，…，ui+p+1確定的p次 B 樣條基函數，p為冪次。顯然p次B樣條函數僅與當前p+2個節點有關，具有局域控制特性。

根據Chen［12］的論證和實驗，二次和三次B樣條函數均適合EMD的均值擬合，其中二次B樣條函數的計算量較小。采用B樣條函數擬合均值時，將序列的極值點位置作為節點向量。假設u1，…，um為序列x（t）的極值點位置，則由二次B樣條函數擬合而成的序列x（t）的均值線為：

由于B樣條函數的局域控制特性，均值線m（t）可以滑動平均方式分段求取，即去均值操作可分段進行。基于上述概念，提出如下快速算法：

（1）初始化：周期閾值T，SD閾值T_SD，IMF個數I=0；

（2）從頭或接上一循環遍歷x（t），直到找到k個極值點；

（3）若是第一次找到k個極值點，將其確定的周期與T比較，若周期小于T，則I=I+1，即IMF個數增加1，并將由步驟（2）遍歷的x（t）賦值給r0（t），否則終止分解；

（4）在當前“局域”內計算前I個IMF：

① 初始化：i=0；

② （a）令i=i+1，若i＞I，轉到（2）；若i≤I，進行以下操作；

（b）在當前“局域”內對ri-1（t）進行均值擬合并去局域均值；

（c）微調：計算局域標準差local_SD，若local_SD＞T_SD，則繼續均值擬合和去均值，直到滿足local_SD＜T_SD，從而得到ri（t）=ri-1（t）-imfi（t）；

（d）從頭或接上一循環遍歷ri（t），直到找到k個極值點；

（e）若這是第一次找到k個極值點，將其確定的周期與T比較，若周期小于T，則I=I+1；

（f）轉到（a）；

（5）當r（t）遍歷完成后，結束。

局域標準差local_SD的計算方法與經典EMD相同，但僅涉及“局域”范圍內序列，因此local_SD集中體現了原型模態函數的局部對稱性。若不滿足對稱性要求，可在“局域”內對序列進行微調。由于k值一般較小（4～20），當檢測到新到達的序列有k個極值時，“局域”便向前移動，即快速算法具有邊采集數據變作分解的特點，無須像經典算法累積至一定長度。另外，三次樣條函數兩端呈開放式，而由式（3）和式（4）B樣條函數和均值線兩端均收斂于零值，從而相鄰“局域”產生的IMF是連續的。因此，快速算法不僅實時性好，而且IMF不會因分段產生不連續，適合流數據分析。

圖2 梯形“局域”示意圖Fig.2 The trapeziform local area

3 EMD快速算法誤差分析

擬采用一個調頻調幅非線性仿真信號進行分析，其解析表達式為：

該信號由一基頻為30 Hz、調制頻率為15 Hz的調頻調幅非線性信號和一頻率為120 Hz正弦信號疊加而成。在理想的分解狀態下，一階和二階IMF應該無限接近x1（t）和x2（t）。

為集中研究該算法的效果，這里不對信號進行延拓，完成分解后，將受端點效應影響的數據置零，比較分解結果和實際信號分量的誤差。另外，采樣率定為4 000 Hz。圖3和圖4分別是由快速算法和經典算法分解得到的各分量及其誤差，其中誤差定義為各分量和實際信號分量的差：

圖3 由快速EMD算法分解產生的各分量及其誤差Fig.3 Components and errors by fast EMD

圖4 由經典EMD算法分解產生的各分量及其誤差Fig.4 Components and errors by classical EMD

這里定義相對誤差以衡量分解的準確性。各階IMF的相對誤差定義為：

當原信號存在趨勢項時，余項的相對誤差定義同各階IMF的相對誤差，即：

當原信號無趨勢項時，將余項視為誤差本身，余項與原信號的比值作為相對誤差，即：

總相對誤差為各階IMF相對誤差和余項相對誤差之和：

圖5為快速算法和經典算法的相對誤差比較，快速算法的精確度略有下降，但與經典EMD算法保持在同一數量級。另外，相對誤差隨局域寬度增加而減少。

圖5 快速EMD和經典EMD的分解誤差比較Fig.5 Errors by fast EMD and classical EMD

4 EMD快速算法復雜度分析

4.1 時間復雜度分析

顯然，兩種算法的問題規模均為L，L為序列的長度。對序列從頭至尾進行一次“去均值”作為基本操作，將兩種算法的基本操作的時間復雜度記為T1和T2。

根據Cox-de Boor遞推定義，若要求得二次B樣條函數，必須從零次B樣條函數遞推得到。對于長度為L的序列，各次B樣條函數的計算量如表1所示。其中零次B樣條函數默認為1，無須計算。根據式（4），將所有二次B樣條函數合成為均值曲線需要2L次加法操作和3L次乘法操作，去均值需要L次加法操作，所以快速算法的時間復雜度為：

其中，O+（）+O×（）分別表示加法時間復雜度和乘法時間復雜度。

表1 各次B樣條函數的時間復雜度Tab.1 Time complexity of B splines of each order

對于經典算法，擬合的關鍵是求出三次樣條函數的分段系數，通常采用三彎矩法。三彎矩法實現上包絡線擬合可分為四個步驟：

（1）建立三彎矩方程組；

（2）采用追趕法求解上述方程組；

（3）計算三次樣條函數的分段系數；

（4）計算每個采樣點對應的包絡值。

前三個步驟的綜合時間復雜度和極值點總數m有關，記為T*。假設由步驟（3）得到上包絡線的分段系數為ai、bi、ci和di，即上包絡線為：

其中m1+m2=m，m1是極大值點數，m2是極小值點數。那么擬合上包絡線的計算量是6L次乘法操作和3L次加法操作。插值下包絡線的計算量與上包絡線相等，然后擬合均值線并去均值，因此，一次基本操作的復雜度為：

一般認為一次乘法操作的耗時是一次加法操作的4倍，所以T2＞T1。另外，三彎矩法的實現過程復雜，因此事實上T2比T1大得多。

綜上，快速算法的時間復雜度為n（x+1）T1，其中n是IMF個數，x表示由微調產生的計算量，大小與local_SD準則有關，一般是一個小于1的正數。經典算法的時間復雜度為n（y+1）T2，y的大小與SD準則有關，一般y大于1。所以，快速算法的時間復雜度大為下降。取一段由戰場偵察探測器記錄的履帶車音頻信號（長度L可變），信號帶寬為20 Hz～300 Hz，采樣率為2 000 Hz，分別采用經典算法和快速算法作分解。其中采用經典算法時，未作分段處理。圖6所示為兩種算法的分解時間t和信號長度L的關系。隨著長度增加，經典算法和快速算法的分解時間均呈上升趨勢，但是前者的上升速度遠大于后者。因此，在數據量較大的情況下，快速算法的優勢較為明顯。另外，快速算法的局域寬度系數k對分解速度略有影響，當局域寬度較小時速度較快。

圖6 快速算法和經典算法的分解時間比較Fig.6 Time consumption of fast EMD and classical EMD

圖7 由k個極值點確定的一階B樣條函數Fig.7 B splines of one order decided by k extrema

4.2 空間復雜度分析

對于長度為L的序列，假設最終分解為n個IMF和一個殘余項之和，那么兩種算法最終分解結果所占的存儲空間是相同的。不同算法的空間復雜度區別在于分解過程中占用的工作空間不同。以下分析兩種算法的工作空間復雜度S1和S2。

快速算法只須記錄“局域”內的一階和二階B樣條函數、均值和原型模態函數。如圖7，有k個極值點的局域內，假設其寬度為l，記錄一階B樣條函數需要2l個存儲單元空間。同理，記錄二階B樣條函數需要3l個存儲單元空間，均值和原型模態函數的空間與B樣條函數復用，則：

經典算法須記錄整個序列的上下包絡、均值和原型模態函數，記錄上下包絡需要2L個存儲單元，均值和原型模態函數的空間與上下包絡復用，那么：

5 結論

提出一種基于二次B樣條函數和局域SD準則的快速EMD算法，并闡述了該算法的原理，分析了該算法的分解精度和復雜度。相比經典算法，快速算法具備如下特征和優點：

（1）無須對序列進行分段，避免引發新的端點效應；

（2）IMF自始至終連續；

（3）邊采集邊分解，實時性好，適用流數據處理；

（4）相比經典算法，分解精度基本維持不變，復雜度顯著降低。

另外，快速算法的精度和復雜度依賴于局域寬度系數k。當k較大時，快速算法退化為經典算法，但仍然得到連續的IMF。復雜度隨著k的下降而降低，但會損失分解精度。在實際應用中，應視具體情況決定k的取值。

［1］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition method and the hilbert spectrum for non-stationary time series analysis［J］.Proc Roy Soc London，1998，454：903-995.

［2］Yannis K，Stephen M.Development of EMD-based denoising method inspired by wavelet thresholding［J］.IEEE Transactions on Signal Processing，2009，57（4）：1351-1362.

［3］ Wu S，Huang W，Kong F，et al.Extracting power transformer vibration features by a time-scale-frequency analysis method［J］.J Electromagnetic Analysis ＆ Applications，2010，2：31-38.

［4］Abhyankar A，Schuckers S.Modular decomposition of fingerprint time series captures for the liveness check［J］.International Journal of Computer and Electrical Engineering，2010，2（3）：1793-8163.

［5］Cai S C，Zhang G W.FeatureAbstracting and identification of acoustic target in the battle field based on EMD［J］.Journal of Shanghai Jiaotong University（Science），2007，E-12（4）：525-529.

［6］呂艷新，孫書學，顧曉輝.基于EMD和能量比的戰場聲目標分類與識別［J］.振動與沖擊，2008，27（11）：51-55.

［7］Damerval C，Meignen S，Perrier V.A fast algorithm for bidimensional EMD［J］.IEEE Signal Processing Letters，2005，12（10）：701-704.

［8］胡勁松，楊世錫.基于有效數據的經驗模態分解快速算法研究［J］.振動、測試與診斷，2006，26（2）：119-121.

［9］Qin S R，Qin Y，Mao Y F.Fast Implementation of orthogonal empirical mode decomposition and its application into singular signal detection［C］.IEEE International Conference on Signal Processing and Communications.Dubai United Arab E-mirates，2007：1215-1218.

［10］張立振.分解信號為正交本征模態函數的方法［J］.振動與沖擊，2007，26（5）：27-32.

［11］ Rilling G，Flandrin P，Goncalvés P.On empirical mode decomposition and its algorithms［C］.IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03，Grado，Italy，2003：8-11.

［12］ Chen Q，Huang N，Riemenschneider S，et al.A b-spline approach for empirical mode decompositions［J］.Advances in Computational Mathematics，2006，24（1）：171-195.