一種適用于流數(shù)據(jù)分析的快速EMD算法

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

（中國科學(xué)院 上海微系統(tǒng)與信息技術(shù)研究所 無線傳感網(wǎng)實(shí)驗(yàn)室，上海 200050）

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

（1）采用更加簡(jiǎn)便的擬合方法，適當(dāng)放寬終止準(zhǔn)則［7］；

（2）只對(duì)有效數(shù)據(jù)進(jìn)行擬合、和進(jìn)行終止準(zhǔn)則判定［8］；

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

但是這些方法未能擺脫經(jīng)典EMD的分段處理法，它們的共同特點(diǎn)是：

（1）分段處理，劃分處引發(fā)端點(diǎn)效應(yīng)；

（2）若進(jìn)行重疊劃分，由于各段去均值的次數(shù)不同，劃分處的IMF不連續(xù)；

（3）IMF從低階到高階依次產(chǎn)生，當(dāng)?shù)碗AIMF完全確定后才能開始高一階IMF的計(jì)算；

（4）全局終止準(zhǔn)則，不能體現(xiàn)局部偏離情況，若不滿足終止準(zhǔn)則，須重新進(jìn)行全局去均值。

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

1 經(jīng)典EMD方法基本原理

EMD方法能把非平穩(wěn)、非線性時(shí)間序列分解成一組由序列本身決定的數(shù)據(jù)序列集，即本征模態(tài)函數(shù)（IMF）。IMF須滿足2個(gè)條件：

（1）極值點(diǎn)和過零點(diǎn)數(shù)目必須相等或至多相差一點(diǎn)；

（2）對(duì)稱性，由局部極大點(diǎn)構(gòu)成的包絡(luò)線和局部極小點(diǎn)構(gòu)成的包絡(luò)線的平均值為零。

其本質(zhì)是通過分辨特征時(shí)間尺度獲得本征振動(dòng)模式，然后由本征振動(dòng)模式來分解時(shí)間序列的一種算法。對(duì)時(shí)間序列x（t）進(jìn)行分解的基本步驟如下：

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

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

（4）若ri（t）滿足繼續(xù)分解的條件，則i=i+1，轉(zhuǎn)到（2）；否則分解結(jié)束，ri（t）是殘余分量。

最終分解結(jié)果為：

即原始序列可分解為n個(gè)IMF和一個(gè)殘余項(xiàng)之和。

2 快速EMD算法基本原理

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

以“局域”替代經(jīng)典EMD的分段處理技術(shù)，將“局域”作為當(dāng)前運(yùn)算范圍，一個(gè)“局域”的長(zhǎng)度遠(yuǎn)小于一個(gè)分段的長(zhǎng)度，從而減少無效運(yùn)算。采用B樣條函數(shù)以滑動(dòng)平均方式擬合均值，在“局域”內(nèi)使用SD終止準(zhǔn)則，即局域SD準(zhǔn)則。

圖1 局部非對(duì)稱的原型模態(tài)函數(shù)Fig.1 A proto-mode function with local asymmetry

設(shè)U是m個(gè)非遞減數(shù)的集合，u1≤u2≤u3≤…≤um，ui稱為節(jié)點(diǎn)，集合U稱為節(jié)點(diǎn)向量，半開區(qū)間［ui，ui+1）是第i個(gè)節(jié)點(diǎn)區(qū)間。B樣條函數(shù)的Cox-de Boor遞推定義為：

其中，Ni，p（u）是由節(jié)點(diǎn)ui，…，ui+p+1確定的p次 B 樣條基函數(shù)，p為冪次。顯然p次B樣條函數(shù)僅與當(dāng)前p+2個(gè)節(jié)點(diǎn)有關(guān)，具有局域控制特性。

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

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

（1）初始化：周期閾值T，SD閾值T_SD，IMF個(gè)數(shù)I=0；

（2）從頭或接上一循環(huán)遍歷x（t），直到找到k個(gè)極值點(diǎn)；

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

（4）在當(dāng)前“局域”內(nèi)計(jì)算前I個(gè)IMF：

① 初始化：i=0；

② （a）令i=i+1，若i＞I，轉(zhuǎn)到（2）；若i≤I，進(jìn)行以下操作；

（b）在當(dāng)前“局域”內(nèi)對(duì)ri-1（t）進(jìn)行均值擬合并去局域均值；

（c）微調(diào)：計(jì)算局域標(biāo)準(zhǔn)差local_SD，若local_SD＞T_SD，則繼續(xù)均值擬合和去均值，直到滿足local_SD＜T_SD，從而得到ri（t）=ri-1（t）-imfi（t）；

（d）從頭或接上一循環(huán)遍歷ri（t），直到找到k個(gè)極值點(diǎn)；

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

（f）轉(zhuǎn)到（a）；

（5）當(dāng)r（t）遍歷完成后，結(jié)束。

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

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

3 EMD快速算法誤差分析

擬采用一個(gè)調(diào)頻調(diào)幅非線性仿真信號(hào)進(jìn)行分析，其解析表達(dá)式為：

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

為集中研究該算法的效果，這里不對(duì)信號(hào)進(jìn)行延拓，完成分解后，將受端點(diǎn)效應(yīng)影響的數(shù)據(jù)置零，比較分解結(jié)果和實(shí)際信號(hào)分量的誤差。另外，采樣率定為4 000 Hz。圖3和圖4分別是由快速算法和經(jīng)典算法分解得到的各分量及其誤差，其中誤差定義為各分量和實(shí)際信號(hào)分量的差：

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

圖4 由經(jīng)典EMD算法分解產(chǎn)生的各分量及其誤差Fig.4 Components and errors by classical EMD

這里定義相對(duì)誤差以衡量分解的準(zhǔn)確性。各階IMF的相對(duì)誤差定義為：

當(dāng)原信號(hào)存在趨勢(shì)項(xiàng)時(shí)，余項(xiàng)的相對(duì)誤差定義同各階IMF的相對(duì)誤差，即：

當(dāng)原信號(hào)無趨勢(shì)項(xiàng)時(shí)，將余項(xiàng)視為誤差本身，余項(xiàng)與原信號(hào)的比值作為相對(duì)誤差，即：

總相對(duì)誤差為各階IMF相對(duì)誤差和余項(xiàng)相對(duì)誤差之和：

圖5為快速算法和經(jīng)典算法的相對(duì)誤差比較，快速算法的精確度略有下降，但與經(jīng)典EMD算法保持在同一數(shù)量級(jí)。另外，相對(duì)誤差隨局域?qū)挾仍黾佣鴾p少。

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

4 EMD快速算法復(fù)雜度分析

4.1 時(shí)間復(fù)雜度分析

顯然，兩種算法的問題規(guī)模均為L(zhǎng)，L為序列的長(zhǎng)度。對(duì)序列從頭至尾進(jìn)行一次“去均值”作為基本操作，將兩種算法的基本操作的時(shí)間復(fù)雜度記為T1和T2。

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

其中，O+（）+O×（）分別表示加法時(shí)間復(fù)雜度和乘法時(shí)間復(fù)雜度。

表1 各次B樣條函數(shù)的時(shí)間復(fù)雜度Tab.1 Time complexity of B splines of each order

對(duì)于經(jīng)典算法，擬合的關(guān)鍵是求出三次樣條函數(shù)的分段系數(shù)，通常采用三彎矩法。三彎矩法實(shí)現(xiàn)上包絡(luò)線擬合可分為四個(gè)步驟：

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

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

（3）計(jì)算三次樣條函數(shù)的分段系數(shù)；

（4）計(jì)算每個(gè)采樣點(diǎn)對(duì)應(yīng)的包絡(luò)值。

前三個(gè)步驟的綜合時(shí)間復(fù)雜度和極值點(diǎn)總數(shù)m有關(guān)，記為T*。假設(shè)由步驟（3）得到上包絡(luò)線的分段系數(shù)為ai、bi、ci和di，即上包絡(luò)線為：

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

一般認(rèn)為一次乘法操作的耗時(shí)是一次加法操作的4倍，所以T2＞T1。另外，三彎矩法的實(shí)現(xiàn)過程復(fù)雜，因此事實(shí)上T2比T1大得多。

綜上，快速算法的時(shí)間復(fù)雜度為n（x+1）T1，其中n是IMF個(gè)數(shù)，x表示由微調(diào)產(chǎn)生的計(jì)算量，大小與local_SD準(zhǔn)則有關(guān)，一般是一個(gè)小于1的正數(shù)。經(jīng)典算法的時(shí)間復(fù)雜度為n（y+1）T2，y的大小與SD準(zhǔn)則有關(guān)，一般y大于1。所以，快速算法的時(shí)間復(fù)雜度大為下降。取一段由戰(zhàn)場(chǎng)偵察探測(cè)器記錄的履帶車音頻信號(hào)（長(zhǎng)度L可變），信號(hào)帶寬為20 Hz～300 Hz，采樣率為2 000 Hz，分別采用經(jīng)典算法和快速算法作分解。其中采用經(jīng)典算法時(shí)，未作分段處理。圖6所示為兩種算法的分解時(shí)間t和信號(hào)長(zhǎng)度L的關(guān)系。隨著長(zhǎng)度增加，經(jīng)典算法和快速算法的分解時(shí)間均呈上升趨勢(shì)，但是前者的上升速度遠(yuǎn)大于后者。因此，在數(shù)據(jù)量較大的情況下，快速算法的優(yōu)勢(shì)較為明顯。另外，快速算法的局域?qū)挾认禂?shù)k對(duì)分解速度略有影響，當(dāng)局域?qū)挾容^小時(shí)速度較快。

圖6 快速算法和經(jīng)典算法的分解時(shí)間比較Fig.6 Time consumption of fast EMD and classical EMD

圖7 由k個(gè)極值點(diǎn)確定的一階B樣條函數(shù)Fig.7 B splines of one order decided by k extrema

4.2 空間復(fù)雜度分析

對(duì)于長(zhǎng)度為L(zhǎng)的序列，假設(shè)最終分解為n個(gè)IMF和一個(gè)殘余項(xiàng)之和，那么兩種算法最終分解結(jié)果所占的存儲(chǔ)空間是相同的。不同算法的空間復(fù)雜度區(qū)別在于分解過程中占用的工作空間不同。以下分析兩種算法的工作空間復(fù)雜度S1和S2。

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

經(jīng)典算法須記錄整個(gè)序列的上下包絡(luò)、均值和原型模態(tài)函數(shù)，記錄上下包絡(luò)需要2L個(gè)存儲(chǔ)單元，均值和原型模態(tài)函數(shù)的空間與上下包絡(luò)復(fù)用，那么：

5 結(jié)論

提出一種基于二次B樣條函數(shù)和局域SD準(zhǔn)則的快速EMD算法，并闡述了該算法的原理，分析了該算法的分解精度和復(fù)雜度。相比經(jīng)典算法，快速算法具備如下特征和優(yōu)點(diǎn)：

（1）無須對(duì)序列進(jìn)行分段，避免引發(fā)新的端點(diǎn)效應(yīng)；

（2）IMF自始至終連續(xù)；

（3）邊采集邊分解，實(shí)時(shí)性好，適用流數(shù)據(jù)處理；

（4）相比經(jīng)典算法，分解精度基本維持不變，復(fù)雜度顯著降低。

另外，快速算法的精度和復(fù)雜度依賴于局域?qū)挾认禂?shù)k。當(dāng)k較大時(shí)，快速算法退化為經(jīng)典算法，但仍然得到連續(xù)的IMF。復(fù)雜度隨著k的下降而降低，但會(huì)損失分解精度。在實(shí)際應(yīng)用中，應(yīng)視具體情況決定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］呂艷新，孫書學(xué)，顧曉輝.基于EMD和能量比的戰(zhàn)場(chǎng)聲目標(biāo)分類與識(shí)別［J］.振動(dòng)與沖擊，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］胡勁松，楊世錫.基于有效數(shù)據(jù)的經(jīng)驗(yàn)?zāi)B(tài)分解快速算法研究［J］.振動(dòng)、測(cè)試與診斷，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］張立振.分解信號(hào)為正交本征模態(tài)函數(shù)的方法［J］.振動(dòng)與沖擊，2007，26（5）：27-32.

［11］ Rilling G，F(xiàn)landrin 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.