魯鐵定 謝建雄

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

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

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

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

1 算法原理

1.1 EEMD算法原理

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

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

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

(1)

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

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

(2)

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

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

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

(3)

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

1.2 多尺度排列熵(MPE)

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

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

(4)

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

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

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

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

2 仿真算例分析

由于GPS高程時間序列包含的噪聲種類較多，頻率分布范圍較廣，為保證與實際情況相符，需構(gòu)造復(fù)合仿真信號來模擬GPS高程時間序列的實測信號。仿真信號的采樣頻率為365 Hz，采樣點個數(shù)為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與仿真信號的相關(guān)系數(shù)Fig.2 Correlation coefficients of each IMF and simulated signal

圖3 各階IMF的值

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

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

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

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

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

3 實例分析

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

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

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

圖6 各階IMF與原始序列的相關(guān)系數(shù)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法降噪效果優(yōu)于相關(guān)系數(shù)法，而本文MPE-閾值法降噪后的序列與原始時間序列的變化趨勢更吻合，在最大限度濾除噪聲的基礎(chǔ)上較好地保留了時間序列的有用成分，與另外兩種方法相比，降噪效果最佳。對于NVSK站，由于相關(guān)系數(shù)法無法準確判斷出噪聲分界層數(shù)，該方法失效；而MPE法盡管能夠準確判斷出噪聲分界層數(shù)，但不能識別出混合IMF，導(dǎo)致降噪效果欠佳；本文MPE-閾值法能夠?qū)⒒旌螴MF中殘留的噪聲去除，同時保留信號的有用成分，降噪的結(jié)果優(yōu)于MPE法。

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

dnSNR=10 lg(Ps/Pg)

(6)

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

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

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

4 結(jié) 語

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