于宏旭 文漢江 劉煥玲 董 杰 藺文奇

1 遼寧工程技術(shù)大學(xué)測(cè)繪與地理科學(xué)學(xué)院，遼寧省阜新市玉龍路88號(hào)，123000 2 中國(guó)測(cè)繪科學(xué)研究院，北京市蓮花池西路28號(hào)，100036

常用的GNSS高程方向去噪方法有Kalman濾波[1]、小波去噪[2-3]、經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)及其改進(jìn)算法[4-7]等。但這幾種方法都存在不足：Kalman濾波在處理非線性時(shí)間序列時(shí)效果不佳；小波去噪在設(shè)置閾值時(shí)只考慮了高頻噪聲的影響，對(duì)于低頻噪聲的去噪效果并不理想[8]；CEEMDAN算法舍棄高頻分量，保留低頻分量，容易丟失有用信息且忽略了低頻分量中包含的噪聲。小波包多閾值分解基于頻率高低分段選取合適的閾值，能夠高效剔除各頻段的噪聲，同時(shí)保留高頻段中的有用信息[9]。

基于以上研究，本文結(jié)合CEEMDAN和小波包多閾值去噪的優(yōu)勢(shì)，引入一種CEEMDAN+小波包多閾值方法對(duì)GNSS高程時(shí)間序列進(jìn)行去噪。

1 算法原理

1.1 CEEMDAN算法原理

CEEMDAN算法原理如下[10]：

(1)

(2)

殘差為：

(3)

(4)

殘差為：

(5)

3)按照上述2個(gè)步驟進(jìn)行分解，直到待分解信號(hào)只有2個(gè)極值點(diǎn)，即為單調(diào)函數(shù)時(shí)停止計(jì)算，此時(shí)CEEMDAN中原始待分解信號(hào)可表示為：

(6)

1.2 小波包多閾值原理

小波包多閾值變換包括小波包分解、閾值處理、小波重構(gòu)等3個(gè)步驟[11]。最優(yōu)小波基、最佳分解層數(shù)和小波包去噪閾值的選取是小波包分解的關(guān)鍵因素。本文選擇較常用的dbN和symN小波基函數(shù)。分解層數(shù)的選擇對(duì)小波包去噪效果尤為重要，任何信號(hào)都存在一個(gè)去噪效果最好的分解層數(shù)，通常分解層數(shù)越大，噪聲和信號(hào)表現(xiàn)的不同特性越明顯，越有利于二者的分離；但分解層數(shù)越大，重構(gòu)信號(hào)失真越嚴(yán)重，會(huì)在一定程度上影響去噪效果。常用信噪比(SNR)和均方根誤差(RMSE)綜合確定分解層數(shù)：

(7)

(8)

式中，x(i)為原始信號(hào)，x′(i)為降噪后的信號(hào)，m為信號(hào)長(zhǎng)度。RMSE反映去噪后信號(hào)與原始信號(hào)之間的區(qū)別，RMSE越小則信號(hào)去噪效果越好；SNR反映信號(hào)與噪聲的比例，重構(gòu)信號(hào)的信噪比越高則信號(hào)的重構(gòu)效果越好。

小波包多閾值去噪方法可根據(jù)頻率的高低選擇閾值，避免了單閾值去噪不能充分考慮信號(hào)與噪聲分布、極易去掉中高頻信號(hào)中有用信息、產(chǎn)生過度去噪的缺點(diǎn)。目前小波包分解常用的閾值準(zhǔn)則有固定閾值(sqtwolog)準(zhǔn)則、自適應(yīng)閾值(rigrsure)準(zhǔn)則與極大極小閾值(minimaxi)準(zhǔn)則。sqtwolog準(zhǔn)則是將小波包系數(shù)全部置0，能夠較強(qiáng)地去除噪聲，適合處理高頻信號(hào)；rigrsure準(zhǔn)則是基于無偏似然估計(jì)原理的自適應(yīng)閾值準(zhǔn)則，僅將部分系數(shù)置0，適合處理中低頻信號(hào)；minimaxi準(zhǔn)則是固定形式的閾值準(zhǔn)則，是一種較為保守的處理方法，適合處理中低頻信號(hào)。因此，本文在GNSS高程時(shí)間序列去噪時(shí)，采用rigrsure準(zhǔn)則處理低頻段，采用minimaxi準(zhǔn)則處理中頻段，采用sqtwolog準(zhǔn)則處理高頻段[12]。對(duì)于閾值函數(shù)類型，選取最常用的軟閾值函數(shù)。本文去噪流程見圖1。

圖1 本文去噪流程Fig.1 Denoising flow chart in this paper

2 仿真實(shí)驗(yàn)

采用MATLAB開展仿真實(shí)驗(yàn)。因高程時(shí)間序列季節(jié)性變化近似于正弦曲線，且包含周期項(xiàng)、半周期項(xiàng)和趨勢(shì)項(xiàng)等信息，因此將y=6sin(0.01t)+2cos(0.03t)-0.005t+f(t)+w(t)作為原始信號(hào)，其中，f為有色噪聲，w為白噪聲，信號(hào)長(zhǎng)度為3 650。根據(jù)高程時(shí)間序列的噪聲特點(diǎn)，將有色噪聲振幅設(shè)置為白噪聲振幅的2倍。利用CEEMDAN對(duì)信號(hào)進(jìn)行分解，在添加高斯白噪聲時(shí)，為削弱添加的白噪聲對(duì)降噪結(jié)果的影響，設(shè)置整體平均次數(shù)為100，添加的白噪聲信噪比為0.2[13]。信號(hào)、添加的白噪聲和有色噪聲分別如圖2(a)～2(c)所示，加入噪聲后的混合信號(hào)如圖2(d)所示。

圖2 仿真實(shí)驗(yàn)的模擬信號(hào)Fig.2 The analog signal used in the simulationexperiment

利用CEEMDAN分解得到10個(gè)頻率逐漸降低的IMF分量和1個(gè)殘余項(xiàng)(圖3)。從圖3看出，隨著信號(hào)分解層數(shù)的增加，噪聲逐漸減少。

單獨(dú)使用CEEMDAN去噪時(shí)，利用分界點(diǎn)將IMF分為噪聲IMF和信號(hào)IMF。然而，在噪聲IMF中仍有部分有用信息，在信號(hào)IMF中也包含部分噪聲。為克服該缺點(diǎn)，本文不進(jìn)行信號(hào)與噪聲分界點(diǎn)的判斷，而是采用小波包多閾值法對(duì)每個(gè)IMF分量進(jìn)行小波包多閾值去噪。首先計(jì)算每個(gè)IMF分量中的根節(jié)點(diǎn)能量占IMF總能量的百分比，根據(jù)能量不同將節(jié)點(diǎn)分為低、中、高頻3個(gè)部分，然后分別采用rigrsure、minimaxi和sqtwolog準(zhǔn)則進(jìn)行數(shù)據(jù)處理。

在利用小波包多閾值去噪時(shí)，選擇對(duì)稱性與正則性較好的sym小波(sym4、sym5、sym6)和db小波(db3、db4、db5)作為小波基函數(shù)，分解層數(shù)定為3、4、5層。利用RMSE、SNR評(píng)價(jià)最終降噪結(jié)果。經(jīng)過多次對(duì)比實(shí)驗(yàn)(圖4)，發(fā)現(xiàn)小波基選擇sym5、分解層數(shù)為4層時(shí)得到的RMSE最小，SNR最大，說明此時(shí)小波包多閾值去噪效果最優(yōu)。

圖4 利用RMSE和SNR判斷最佳小波基與分解層數(shù)Fig.4 Using RMSE and SNR to determine the best wavelet basis and the number of decomposition layers

利用小波包多閾值去噪對(duì)各IMF分量去噪，將去噪后的IMF分量與殘余項(xiàng)進(jìn)行重構(gòu)，得到去噪后的信號(hào)，結(jié)果見圖5。

圖5 CEEMDAN+小波包多閾值去噪結(jié)果和噪聲Fig.5 CEEMDAN and wavelet packet multi-thresholddenoising results and noise

為驗(yàn)證本文方法的去噪效果，分別采用本文方法(CEEMDAN+小波包多閾值)、EMD、CEEMDAN、小波去噪和小波包多閾值去噪分析仿真信號(hào)，結(jié)果見圖6。選擇SNR、RMSE評(píng)價(jià)5種方法的降噪效果，結(jié)果見表1。

圖6 5種去噪方法對(duì)比Fig.6 Comparison of five denoising methods

表1 去噪效果比較

從圖6可以看出，5種方法均能對(duì)模擬信號(hào)去噪，但EMD去噪結(jié)果過于平緩，丟失了部分信息；與小波包多閾值去噪、小波去噪、CEEMDAN去噪相比，CEEMDAN+小波包多閾值方法去噪后RMSE最小、SNR最大，說明該方法去噪效果最優(yōu)。

3 實(shí)例分析

選取JPL發(fā)布的拉薩站(LHAZ)2000～2020年單天解高程時(shí)間序列進(jìn)行降噪分析，原始時(shí)間序列中已剔除序列初值、地震和非地震跳變的影響，高程時(shí)序缺值或者階躍不影響高程序列降噪效果。利用3倍中誤差準(zhǔn)則去除時(shí)間序列中剩余的粗差，提高數(shù)據(jù)質(zhì)量。LHAZ站原始高程時(shí)間序列見圖7。

圖7 LHAZ站預(yù)處理后的高程時(shí)間序列Fig.7 Height time series after preprocessing at LHAZ station

利用CEEMDAN算法將時(shí)間序列分解成12個(gè)IMF分量和1個(gè)殘余項(xiàng)，將各IMF分量進(jìn)行小波包多閾值去噪，得到降噪后的各IMF分量，并與殘余項(xiàng)重構(gòu)得到去噪后的信號(hào)，結(jié)果見圖8。

圖8 LHAZ站CEEMDAN+小波包多閾值去噪信號(hào)和所剔除的噪聲Fig.8 CEEMDAN and wavelet packet multi-threshold denoising signal and eliminated noise at LHAZ station

使用上節(jié)中5種方法對(duì)LHAZ站高程時(shí)間序列進(jìn)行去噪處理，結(jié)果見圖9，去噪效果見表2。結(jié)合圖9和表2可以看出，5種方法均能很好地去除信號(hào)中的噪聲，其中，CEEMDAN+小波包多閾值去噪法的SNR最大、RMSE最小，表明其去噪效果最好。小波去噪過程中固定了小波閾值，只對(duì)低頻信號(hào)進(jìn)行分解，沒有對(duì)高頻信號(hào)進(jìn)行分解，導(dǎo)致無法保留高頻信號(hào)中的有用信息。EMD去噪時(shí)在峰值附近會(huì)產(chǎn)生過度去噪的現(xiàn)象。

圖9 5種去噪方法對(duì)比Fig.9 Comparison of five denoising methods

表2 5種去噪方法效果比較

4 結(jié) 語

GNSS高程時(shí)間序列中含有大量噪聲，本文針對(duì)EMD、CEEMDAN、小波、小波包多閾值等傳統(tǒng)去噪方法去噪不徹底或去噪過度的缺點(diǎn)，嘗試?yán)肅EEMDAN +小波包多閾值的方法對(duì)信號(hào)進(jìn)行去噪，并分別利用仿真信號(hào)、LHAZ站高程時(shí)間序列進(jìn)行實(shí)驗(yàn)。結(jié)果表明，5種去噪方法均可以達(dá)到去噪效果，但本文方法去噪后各項(xiàng)指標(biāo)均最優(yōu)，效果最好。