基于CEEMDAN和小波包多閾值的GNSS高程時間序列去噪方法

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

1 遼寧工程技術大學測繪與地理科學學院，遼寧省阜新市玉龍路88號，123000 2 中國測繪科學研究院，北京市蓮花池西路28號，100036

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

基于以上研究，本文結合CEEMDAN和小波包多閾值去噪的優勢，引入一種CEEMDAN+小波包多閾值方法對GNSS高程時間序列進行去噪。

1 算法原理

1.1 CEEMDAN算法原理

CEEMDAN算法原理如下[10]：

(1)

(2)

殘差為：

(3)

(4)

殘差為：

(5)

3)按照上述2個步驟進行分解，直到待分解信號只有2個極值點，即為單調函數時停止計算，此時CEEMDAN中原始待分解信號可表示為：

(6)

1.2 小波包多閾值原理

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

(7)

(8)

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

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

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

2 仿真實驗

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

圖2 仿真實驗的模擬信號Fig.2 The analog signal used in the simulationexperiment

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

單獨使用CEEMDAN去噪時，利用分界點將IMF分為噪聲IMF和信號IMF。然而，在噪聲IMF中仍有部分有用信息，在信號IMF中也包含部分噪聲。為克服該缺點，本文不進行信號與噪聲分界點的判斷，而是采用小波包多閾值法對每個IMF分量進行小波包多閾值去噪。首先計算每個IMF分量中的根節點能量占IMF總能量的百分比，根據能量不同將節點分為低、中、高頻3個部分，然后分別采用rigrsure、minimaxi和sqtwolog準則進行數據處理。

在利用小波包多閾值去噪時，選擇對稱性與正則性較好的sym小波(sym4、sym5、sym6)和db小波(db3、db4、db5)作為小波基函數，分解層數定為3、4、5層。利用RMSE、SNR評價最終降噪結果。經過多次對比實驗(圖4)，發現小波基選擇sym5、分解層數為4層時得到的RMSE最小，SNR最大，說明此時小波包多閾值去噪效果最優。

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

利用小波包多閾值去噪對各IMF分量去噪，將去噪后的IMF分量與殘余項進行重構，得到去噪后的信號，結果見圖5。

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

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

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

表1 去噪效果比較

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

3 實例分析

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

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

利用CEEMDAN算法將時間序列分解成12個IMF分量和1個殘余項，將各IMF分量進行小波包多閾值去噪，得到降噪后的各IMF分量，并與殘余項重構得到去噪后的信號，結果見圖8。

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

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

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

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

4 結 語

GNSS高程時間序列中含有大量噪聲，本文針對EMD、CEEMDAN、小波、小波包多閾值等傳統去噪方法去噪不徹底或去噪過度的缺點，嘗試利用CEEMDAN +小波包多閾值的方法對信號進行去噪，并分別利用仿真信號、LHAZ站高程時間序列進行實驗。結果表明，5種去噪方法均可以達到去噪效果，但本文方法去噪后各項指標均最優，效果最好。