基于經驗模式分解的濾波去噪方法研究進展

張恒璟 ，程鵬飛

(1. 遼寧工程技術大學 測繪學院，遼寧 阜新 123000； 2. 中國測繪科學研究院，北京 100830；3. 國家測繪產品質量檢驗測試中心，北京 100830)

一、引 言

1998年黃鄂教授提出了研究非線性非平穩信號分析技術的希爾伯特-黃變換方法[1](Hilbert-Huang transformation，HHT)，它包括經驗模式分解(empirical mode decomposition，EMD)和希爾伯特變換兩個步驟，已經被廣泛應用于爆破振動信號、機械振動信號、橋梁振動信號、陣列聲波信號、地震波信號、電力諧波信號、語音信號、金融時間序列信號等社會生產建設的諸多領域。

EMD具有類似于小波分解的濾波性質，國際上已有不少關于EMD濾波去噪算法的研究結果，但EMD方法本身還不能用完整的數學公式加以描述或定義，相比較小波完善的理論體系，EMD技術發展并完善的路還很長，且需要更多的人去關注和研究。本文以EMD濾波去噪問題為切入點，總結了近年來專家學者提出的基于EMD的濾波去噪算法，分析了各自算法的適用條件及可能的改進，最后以GPS高程時間序列噪聲提取的具體問題，分析了EMD去噪算法相對于傳統周期擬合模型獲取噪聲的優勢。

二、經驗模式分解及其濾波性質

經驗模式分解的基本方法是：通過不斷的“篩選”，將信號分解為一系列頻率由高到低的本征模態函數分量(intrinsic mode function，IMF)和殘差項，見式(1)，具體的過程可參考文獻[1]。

(1)

式中，N為分解得到IMF的個數；rN(t)為分解的殘差項；ci(t)表示第i個IMF分量。

分解產生的IMF函數需要滿足兩個條件：一是待分解信號中極值點個數與過零點個數相等或至多相差1；二是在任意點上，由局部極大值定義的上包絡線與局部極小值定義的下包絡線的均值為零。這種分解基于數據自身的局部特征，與傅里葉頻譜分析分解成恒定振幅與頻率的一系列正弦函數和小波分解需預先給定基函數相比，EMD是一種自適應的信號分解方法。

由于受信號中斷的影響，EMD會產生模式混疊現象[2]，不僅混淆了信號的時頻分布特性，更重要的是使信號IMF分量的物理意義變得不清晰。為了減弱模式混疊現象的影響，Wu Zhaohua等人在EMD方法的基礎上提出了整體經驗模式分解方法[3](enable EME，EEMD)。這是一種噪聲輔助的數據分析方法，其基本思路是在分解之前，將原始序列信號加入一定信噪比的白噪聲，生成待分解信號，經過EMD分解后生成一系列的IMF分量；由于每次隨機生成的白噪聲都是不同的，這樣反復EMD分解過程，得到多次分解的IMF分量，將對應相同分解尺度的IMF分量取平均，即得到了最終信號分解的IMF分量。

Flandrin研究了分形布朗運動的增量過程，即分形高斯噪聲的EMD分解，獲得了一系列高低頻IMF分量，試驗發現EMD具有類似小波分解的二進濾波性質[4-6]；Wu Zhaohua采用EMD分解高斯白噪聲，通過數值試驗發現EMD等效于一個二進濾波器，且分解的IMF分量都服從正態分布，在以周期對數為橫軸的坐標系中，IMF分量的傅里葉頻譜曲線幾乎相同[7]；Wu Zhaohua又研究了當EMD分解的篩選次數改變時的濾波性質[8]，以白噪聲和delta函數為例，發現當篩選次數從有限次到無限次時，濾波的比率從二進濾波的2逐漸變化為1，但實際上篩選次數不可能無限多，因此比率總是大于1。經驗模式分解類似于小波的濾波性質，為其在信號去噪方面的應用提供了有利的條件。

三、基于EMD的濾波去噪算法

1. 軟閾值法消噪

Donoho等設計了基于小波的軟閾值去噪算法，分離了信號中混合的加性高斯白噪聲[9-10]。基于對EMD濾波性質的試驗結果，類似于小波分解的濾波去噪功能，Boudraa于2005年提出了基于EMD的軟閾值去噪算法[11]，基本原理是：將原始信號x(t)經過EMD分解后，計算每一個IMF分量消噪的軟閾值，分別進行去噪。

首先計算每一個IMF序列的絕對中位差(median absolute deviation，MAD)

MADj=median{|cj(t)-median(cj(t))|}

(2)

式中，median()表示取序列元素從小到大排列后的中位值，如果序列長度為奇數，則取中間的元素；如果序列長度為偶數，則取中間兩個元素的算術平均值。

然后計算IMF所含的噪聲水平(即噪聲中誤差)

(3)

再計算每一個IMF消噪的軟閾值

(4)

式中，L是信號長度。

最后消除IMF分量的噪聲

(5)

軟閾值法消噪的公式中，下標j表示第j個IMF分量，j=1，2，…，N，N是EMD分解得到的IMF分量的個數。

軟閾值法消噪的濾波器設計是基于分解信號中包含有高斯白噪聲的前提，去噪性能比較優越；如果信號中包含有色噪聲的干擾，則設計的濾波器能否有效去噪，還需要進一步試驗研究。

2. 部分IMF重構算法

(6)

此處k如果取1，則去噪信號與原始信號相同，表明原始信號中不包含噪聲。基于EMD分解的連續均方誤差CMSE準則如下

(7)

由式(7)可知，CMSE值相當于第k個IMF分量的能量密度，從數學公式上稱為該IMF分量的均方誤差(mean squared error，MSE)。當k取2～N時可依次計算出CMSE值，將CMSE值第一次顯著發生變化處的k值(記為臨界點js)作為干凈信號的起點，見式(8)，重構后的信號表示見式(9)，臨界點以前的IMF為信號中包含的噪聲，即從原始信號中扣除重構信號得到了部分IMF重構算法的噪聲。

(8)

(9)

部分IMF重構算法在判斷臨界點時，依據的是CMSE值第一次顯著發生變化的位置，這一顯著性是通過試驗觀察得到的，缺少數學規則和證明。EMD分解得到一系列頻率由高到低的IMF分量，每一個IMF的CMSE值如果比較接近，就帶來了判斷臨界點的困難。單純依靠CMSE準則的判斷就需要改進，此時可以考慮信號的實際物理意義，根據信號EMD分解后每一個IMF分量的時頻特性，并結合其他的準則(如平均周期)對部分IMF重構算法作進一步完善。

3. 最優信號重構算法

Weng在2008年提出基于EMD的最優信號重構算法[14](optimal signal reconstruction，OSR)。具體做法是：將原始不含噪聲的信號d(t)加入一定量的高斯白噪聲w(t)后，得到含噪信號x(t)；將含噪信號采用EMD分解成一系列IMF分量和殘差項。假設對IMF分量線性加權，得到去噪后的重構信號如下

(10)

式中，ai表示第i個IMF分量的權重系數。

(11)

J1取得最小值時，令上式對ai求導結果等于0，計算出每一個IMF對應的權重系數。

后半部分原始信號用來測試去噪效果，將后半部分信號進行EMD分解，將權重系數回代到式(10)中，通過計算后半部分重構信號與后半部分原始信號的MSE來衡量去噪性能。

Weng還提出了雙向最優信號重構算法(bidirectional OSR，BOSE)，對不同層次IMF分量加權的同時考慮每個IMF分量中信號之間的水平相關性(即IMF分量各個采樣值之間的相關性)。從Weng采用仿真信號的試驗結果來看，高頻IMF分量的權重系數較小，低頻權重系數較大，可知噪聲主要包含在信號的高頻部分。

OSR方法的局限性是必須知道原始干凈信號，這樣才能基于MSE最小準則估計每個IMF分量的權重系數。但一般情況下試驗得到的原始信號已經包含了噪聲，原始干凈信號無法獲得，因此基于EMD方法的OSR去噪算法在應用中就受到了限制。

4. GPS高程時間序列的濾波去噪

GPS高程時間序列噪聲分析時，首先需要建立合適的序列擬合模型獲取殘差序列，將擬合殘差作為GPS高程序列包含的噪聲。統一的GPS高程時間序列周期擬合模型如下

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

Tkj)ti/τj)H(ti-Tkj)+vi

(12)

式中，ti(i=1，2，…，N)為以年為單位的時間，在此以每日解構成坐標時間序列；待求系數a為序列的平均值；b為線性速率；c、d和e、f分別為年周期和半年周期項的系數；gj則表示由于各種原因引起的階躍式坐標突變(如遠場大地震引起的同震位移、由于儀器或天線變更引起的位移、甚至由于某種有清楚原因引起的點位變化等)；hi和kj表示更加復雜的階躍變化改正項；最后一項vi為GPS高程時間序列周期函數擬合模型的殘差，即下一步分析的序列觀測噪聲。

國內外學者采用式(12)的截斷模型或完整模型，獲得了GPS高程時間序列的噪聲，并以此為基礎，采用極大似然估計方法定量分析噪聲類型和含量[15-22]。筆者也曾分別采用抗差功率譜估計和EMD分解的方法，初步分析了GPS高程時間序列的周期特性[23-24]，功率譜估計結果表明序列中并不存在精確的半年、年或兩年周期項；EMD分解后各個IMF分量的周期也不是一個恒定的數值，而是一個明顯的時變序列，說明了GPS高程時間序列信號的非平穩性特征，這與傳統的研究結果一致。

筆者在文獻[24]中研究了EMD分解方法在GPS高程時間序列分析中的初步應用。GPS高程時間序列信號經過EMD分解為一系列頻率由高到低的本征函數模態分量和殘差項，通過殘差項與鄰近的IMF分量合并，并結合各個IMF方差貢獻率指標，識別GPS高程時間序列信號的非線性趨勢項并分離。中低頻的IMF分量是序列明顯的季節性變化項，包括明顯的近似月、雙月、三個月、半年、年、雙周年信號，但并不存在嚴格精確的周期。因此采用周期擬合模型獲得的序列殘差中仍包含著周期成分，以此作為噪聲的分析結果就存在著偏差。

基于EMD的軟閾值去噪算法，是假設信號中包含加性高斯白噪聲條件下提出的，GPS高程序列信號具有非平穩性特征，序列中既包含白噪聲也包含有色噪聲的影響，這已經被大量研究所證實，因此軟閾值去噪算法應用于GPS高程序列信號去噪，理論上具有局限性，去噪效果需要進一步的試驗驗證。

基于EMD的部分IMF重構算法，通過構建一個低通濾波器，基于CMSE準則，識別出信號中包含的噪聲。通過筆者在文獻[24]中的試驗可知，GPS高程序列信號的低頻部分為序列的非線性趨勢項，中低頻的IMF分量是序列明顯的季節性周期信號，高頻的IMF分量既包含有序列噪聲，還有在目前研究的序列長度區間上暫時無法識別的短期周期項，也可以一并作為噪聲處理。這樣構建的低通濾波器就為濾波GPS序列噪聲提供了方法上的可行性。

四、結束語

采用周期函數模型擬合GPS高程時間序列時采用固定的半年和年周期的做法，與序列實際的周期運動情形不符，傳統功率譜技術的譜分辨率低，小波分解的基底固定的也違背了序列本身的特征，因此傳統的法擬合序列周期項進而獲得序列噪聲的做法存在著缺陷，似乎沒有一個最佳有效的解決方案。

本文將基于EMD分解的信號去噪算法，引入到GPS高程時間序列分析領域，是一個新的嘗試，為GPS高程時間序列信號的周期特性和噪聲研究開辟了一個新的途徑。

參考文獻：

[1] HUANG N E，WU M L，LONG S R，et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J].Proceedings of the Royal Society，1998，454(1971)：903-995.

[2] HUANG N E，SHEN Z，LONG S R. A New View of Nonlinear Water Waves：The Hilbert Spectrum[J]. Annual Review of Fluid Mechanics，1999，31：417-457.

[3] WU Z，HUANG N E. Ensemble Empirical Mode Decomposition:A Noise-assisted Data Analysis Method[J]. Advances in Adaptive Data Analysis，2009，1(1)：1-41.

[4] FLANDRIN P，RILLING G，GONALVES P. Empirical Mode Decomposition as a Filterbank[J]. IEEE Signal Processing Letters，2004，11(2)：112-114.

[5] FLANDRIN P，GONCALVES P，RILLING G. Detrending and Denoising with Empirical Mode Decomposition[C]∥European Signal Processing Conference(EUSIPCO2004).[S.l.]:Cite SeerX.2004.

[6] FLANDRIN P. Empirical Mode Decompositions as Data-driven Wavelet-like Expansions[J].International Journal of Wavelets: Multiresolution and Information Processing，2004，2(4)：1-20.

[7] WU Z，HUANG N E. A Study of the Characteristics of White Noise Using the Empirical Mode Decomposition Method[J]. Proceedings of the Royal Society，2004，460(2046)：1597-1611.

[8] WU Z，HUANG N E. On the Filtering Properties of the Empirical Mode Decomposition[J]. Advances in Adaptive Data Analysis，2010，2(4)：397-414.

[9] DONOHO D L，JOHNSON I M. Ideal Spatial Adaptation by Wavelet Shrinkage[J].Biometrika，1994，81(3)：425-455.

[10] DONOHO D L. De-noising by Soft-thresholding[J].IEEE Transactions on Information Theory,1995，l41(3)：613-627.

[11] BOUDRAA A O，CEXUS J C，SAIDI Z. EMD-based Signal Noise Reduction[Z].World Academy of Science，Engineering and Technology，2005，2:33-36.

[12] 張義平.爆破振動信號的HHT分析與應用[M].北京：冶金工業出版社，2008：83-84.

[13] BOUDRAA A O，CEXUS J C.EMD-based Signal Filtering[J].IEEE Transactions on Instrumentation and Measurement，2007，56(6)：2196-2202.

[14] WENG B，BARNER K E. Optimal Signal Reconstruction Using the Empirical Mode Decomposition[J].EURASIP Journal on Advances in Signal Processing，2008：1-12.

[15] LANGBEIN J，JOHNSON H. Correlated Errors in Geodetic Time Series: Implications for Time-dependent Deformation[J].Journal of Geophysical Research，1997，102(B1)：591-603.

[16] ZHANG J，BOCK Y，JOHNSON H，et al. Southern California Permanent GPS Geodetic Array: Error Analysis of Daily Position Estimates and Site Velocities[J].Journal of Geophysical Research，1997，102(B8)：18035-18055.

[17] NIKOLAIDIS R. Observation of Geodetic and Seismic Deformation with the Global Positioning System[D].Califormia：University of Califormia，2002.

[18] 符養.中國大陸現今地殼形變與GPS坐標時間序列分析[D].上海：中國科學院上海天文臺，2002.

[19] DONG D，FANG P，BOCK Y，et al. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve Expansion Approaches for Regional GPS Network Analysis[J].Journal of Geophysical Research，2006，111(B3)：1-16.

[20] 黃立人，符養.GPS連續觀測站的噪聲分析[J].地震學報，2007，29(2)：197-202.

[21] 袁林果，丁曉利，陳武，等.香港GPS基準站坐標序列特征分析[J].地球物理學報，2008，51(5)：1372-1384.

[22] 蔣志浩，張鵬，秘金鐘，等.顧及有色噪聲影響的CGCS2000下我國CORS站速度估計[J].測繪學報，2010，39(4)：355-363.

[23] 張恒璟，程鵬飛，郭英.我國IGS基準站高程時間序列抗差功率譜估計[J].測繪科學，2012，37(3)：52-53.

[24] 張恒璟，程鵬飛.基于經驗模式分解的CORS站高程時間序列分析[J].大地測量與地球動力學，2012，32(3)：129-134.