EEMD-多尺度排列熵的GPS高程時間序列降噪方法

魯鐵定 謝建雄

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013 2 東華理工大學江西省數字國土重點實驗室，南昌市廣蘭大道418號，330013 3 漳州市測繪設計研究院，福建省漳州市龍溪南路5-58號，363000

GPS高程時間序列中包含各類信號和噪聲，因受測站外部環境、觀測技術誤差等因素的影響，表現出明顯的周期性變化。利用傳統周期函數模型估計的周期項振幅為常數，但實際上，時間序列周期性信號的振幅是隨時間變化的[1-2]，因此采用傳統方法獲得的分析結果不能很好地反映高程時間序列的真實運動。

將隨機噪聲進行有效剔除是獲取合理可靠的周期性變化信號的關鍵。劉煥玲等[3]利用經驗模態分解(EMD)方法分析IGS站高程時間序列，但該方法通過Hilbert頻譜圖識別信息IMF時存在較大的主觀性，缺乏統一的標準；張雙成等[4]結合相關系數原理，應用EMD方法獲得的重構信號與原始時間序列較為一致，但不同測站各分量的噪聲特性存在差異[5]，并且當EMD處理的時間序列存在中斷時會產生較為明顯的模式混疊現象[6]，導致部分站點依據相關系數首個局部極小值識別噪聲IMF分界時出現區別效果不明顯的情況，從而影響降噪結果的準確性。

針對上述問題，綜合考慮整體經驗模態分解(EEMD)能有效減弱模式混疊及多尺度排列熵(MPE)可以反映時間序列復雜程度的優勢，提出一種EEMD-MPE閾值降噪方法，并通過仿真算例和部分IGS站高程時間序列對該方法進行驗證。

1 算法原理

1.1 EEMD算法原理

EEMD[7]的本質是基于高斯白噪聲在所有頻率上具有相等能量的特性，通過在待分解信號中加入高斯白噪聲后再進行多次EMD處理，從而實現減弱或消除模式混疊的方法。EEMD算法的基本步驟如下：

1)向待分解信號中加入隨機高斯白噪聲：

xi(t)=x(t)+nωi(t)

(1)

式中，x(t)為待分解信號；n為加入白噪聲的幅值系數；ωi(t)為加入的白噪聲，i=1,2,…,N。

2)采用EMD處理每個xi(t)，獲得k個IMF和1個余項ri(t)：

(2)

式中，cij(t)為第i次分解得到的第j個IMF。

3)循環步驟1)和2)，但每次均加入不同的隨機高斯白噪聲。

4)計算分解所得的IMF的總體均值，得到EEMD的最終結果：

(3)

關于EEMD的2個重要參數，采用Wu等[7]的建議，設置整體平均次數N=100，白噪聲的幅值系數為0.2。

1.2 多尺度排列熵(MPE)

多尺度排列熵是一種能夠反映信號不確定性和不規則性的非線性分析方法，與排列熵[8]相比具有更好的穩定性和更強的抗噪聲能力，其基本思想是對時間序列多尺度化后，計算對應尺度上的排列熵。多尺度排列熵的計算步驟參考文獻[9]。

為了使多尺度排列熵值位于[0,1]，通常對其進行歸一化處理：

(4)

1.3 基于EEMD-多尺度排列熵(MPE)的降噪新方法

1)利用EEMD處理GPS高程時間序列，獲得一系列IMF和1個余項；

4)GPS高程時間序列的有用信息主要集中在混合IMF的“干凈”部分和信息IMF中，因此對高頻噪聲IMF直接剔除，然后采用改進閾值函數[11]處理混合IMF；

5)重構利用閾值降噪后所得的結果與信息IMF即可獲得最終的輸出序列。

2 仿真算例分析

由于GPS高程時間序列包含的噪聲種類較多，頻率分布范圍較廣，為保證與實際情況相符，需構造復合仿真信號來模擬GPS高程時間序列的實測信號。仿真信號的采樣頻率為365 Hz，采樣點個數為4 000，信號包括3個周期項和隨機噪聲項，其中隨機噪聲項(noise)包含高斯白噪聲和有色噪聲。仿真信號如圖1所示，其表達式為：

Y=9sin(0.6πt)+6cos(2πt)+

11sin(3πt)+noise

(5)

圖1 仿真信號Fig.1 Simulation signal

圖2 各階IMF與仿真信號的相關系數Fig.2 Correlation coefficients of each IMF and simulated signal

圖3 各階IMF的值

為了分析不同方案的處理效果，選用均方根誤差(RMSE)和信噪比(SNR)[12]指標進行降噪質量評價。一般認為，RMSE值越小，降噪后的信號與原始參考信號越接近，降噪質量越好；而SNR則相反，其值越高，降噪效果越顯著。

3種方案降噪后的殘差結果如圖4所示，表1統計了3種方案降噪后的RMSE和SNR值。

圖4 3種方案降噪后的信號與原參考信號的殘差結果Fig.4 Residual results of the signal after thenoise reduction of the three schemes and the original reference signal

表1 不同方法的降噪評價指標

由圖4和表1可知，方案3獲得的降噪后信號與原始參考信號的殘差結果最小，即RMSE值最小，SNR值最大，表明該方案降噪性能最優；而方案1的降噪效果最差，這是由于方案1未能準確識別出噪聲層數；盡管方案2的降噪性能略優于方案1，但該方案將混合IMF當成噪聲IMF進行了剔除，在降噪的同時將分量信號中的有效成分也一并舍去，造成信號失真。

3 實例分析

為進一步檢驗本文方法的有效性，以BJFS和NVSK站高程時間序列(圖5)為例進行降噪分析，數據來源于SOPAC (Scrips Orbit and Permanent Array Center)提供的去除了均值、趨勢項、粗差、同震跳變和非地震跳變影響后的GPS單天高程時間序列。

圖5 BJFS和NVSK站高程時間序列Fig.5 Elevation time series of BJFS and NVSK stations

從圖5可以看出，兩個IGS站高程時間序列中仍存在部分粗差，因此采用“3σ”準則進行粗差識別和剔除，以提高數據的質量。為能有效并準確地提取序列中的周期項，采用3種基于EEMD的降噪方法對高程時間序列進行處理和分析，以驗證本文方法的優越性。

圖6 各階IMF與原始序列的相關系數Fig.6 Correlation coefficients of each IMFand original sequence

圖7 各階IMF的值

圖8 基于不同方法獲取的高程時間序列對比Fig.8 Comparison of elevation time series obtained by different methods

圖9 基于不同方法濾除的噪聲序列對比Fig.9 Comparison of noise sequences filtered by different methods

圖8和9表明，對于BJFS站，MPE法降噪效果優于相關系數法，而本文MPE-閾值法降噪后的序列與原始時間序列的變化趨勢更吻合，在最大限度濾除噪聲的基礎上較好地保留了時間序列的有用成分，與另外兩種方法相比，降噪效果最佳。對于NVSK站，由于相關系數法無法準確判斷出噪聲分界層數，該方法失效；而MPE法盡管能夠準確判斷出噪聲分界層數，但不能識別出混合IMF，導致降噪效果欠佳；本文MPE-閾值法能夠將混合IMF中殘留的噪聲去除，同時保留信號的有用成分，降噪的結果優于MPE法。

由于實測信號的有效信號和噪聲功率均未知，仍采用信噪比評價降噪效果是不準確的，因此本文引入降噪誤差比(dnSNR)[13]指標評價降噪質量，dnSNR值越小降噪效果越顯著。

dnSNR=10 lg(Ps/Pg)

(6)

式中，Ps為含噪信號的功率，Pg為濾除的噪聲功率。同時，選用平滑度(r)[12]作為降噪質量評價指標，r值越小信號序列越光滑，降噪效果越好。表2統計了各種方法的降噪效果評價指標。

表2 各種方法的降噪效果評價指標

由表2可知，僅從平滑度(r)來看，采用本文方法處理BJFS站數據結果的r值最小，僅為0.001 6，表明本文方法的降噪性能最好，同時對比3種方法的dnSNR值也驗證了這一觀點；但對于NVSK站，MPE法和本文方法的r值均為0.003 9，無法體現降噪質量的優劣，而本文方法的dnSNR值小于MPE法，說明本文方法的降噪性能優于MPE法，同時在評價降噪性能方面也反映了dnSNR指標比平滑度(r)指標更具有適用性。

4 結 語

本文綜合EEMD算法和MPE的優勢，提出一種基于改進閾值函數的EEMD-MPE聯合降噪方法。該方法對EEMD處理后的混合IMF進行降噪，并重構降噪后的序列與余下分量，獲得最終降噪后的序列。通過對仿真信號進行詳細的降噪分析，驗證本文方法的降噪性能優于相關系數法和MPE法。將該方法應用于IGS站GPS高程時間序列的降噪處理中，結果表明，本文方法在降噪性能上具有明顯的優勢，獲得的降噪后時間序列更符合基準站的實際運動情況，可為進一步研究地殼運動及分析形變模式提供更加穩定、可靠的數據基礎。