張恒璟 ，程鵬飛

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

一、引 言

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

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

二、經(jīng)驗?zāi)Ｊ椒纸饧捌錇V波性質(zhì)

經(jīng)驗?zāi)Ｊ椒纸獾幕痉椒ㄊ牵和ㄟ^不斷的“篩選”，將信號分解為一系列頻率由高到低的本征模態(tài)函數(shù)分量(intrinsic mode function，IMF)和殘差項，見式(1)，具體的過程可參考文獻[1]。

(1)

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

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

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

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

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

1. 軟閾值法消噪

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

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

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

(2)

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

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

(3)

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

(4)

式中，L是信號長度。

最后消除IMF分量的噪聲

(5)

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

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

2. 部分IMF重構(gòu)算法

(6)

此處k如果取1，則去噪信號與原始信號相同，表明原始信號中不包含噪聲?；贓MD分解的連續(xù)均方誤差CMSE準則如下

(7)

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

(8)

(9)

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

3. 最優(yōu)信號重構(gòu)算法

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

(10)

式中，ai表示第i個IMF分量的權(quán)重系數(shù)。

(11)

J1取得最小值時，令上式對ai求導(dǎo)結(jié)果等于0，計算出每一個IMF對應(yīng)的權(quán)重系數(shù)。

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

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

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

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

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

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

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

(12)

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

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

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

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

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

四、結(jié)束語

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

本文將基于EMD分解的信號去噪算法，引入到GPS高程時間序列分析領(lǐng)域，是一個新的嘗試，為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分析與應(yīng)用[M].北京：冶金工業(yè)出版社，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] 符養(yǎng).中國大陸現(xiàn)今地殼形變與GPS坐標(biāo)時間序列分析[D].上海：中國科學(xué)院上海天文臺，2002.

[19] DONG D，F(xiàn)ANG 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] 黃立人，符養(yǎng).GPS連續(xù)觀測站的噪聲分析[J].地震學(xué)報，2007，29(2)：197-202.

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

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

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

[24] 張恒璟，程鵬飛.基于經(jīng)驗?zāi)Ｊ椒纸獾腃ORS站高程時間序列分析[J].大地測量與地球動力學(xué)，2012，32(3)：129-134.