張恒璟 ，程鵬飛

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

一、引 言

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

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

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

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

(1)

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

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

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

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

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

1. 軟閾值法消噪

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

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

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

(2)

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

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

(3)

再計(jì)算每一個(gè)IMF消噪的軟閾值

(4)

式中，L是信號(hào)長度。

最后消除IMF分量的噪聲

(5)

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

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

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

(6)

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

(7)

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

(8)

(9)

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

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

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

(10)

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

(11)

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

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

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

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

4. GPS高程時(shí)間序列的濾波去噪

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

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

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

(12)

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

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

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

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

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

四、結(jié)束語

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

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

參考文獻(xiàn)：

[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] 張義平.爆破振動(dòng)信號(hào)的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)時(shí)間序列分析[D].上海：中國科學(xué)院上海天文臺(tái)，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é)報(bào)，2007，29(2)：197-202.

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

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

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

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