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

魯鐵定 謝建雄

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

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

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

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

1 算法原理

1.1 EEMD算法原理

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

1)向待分解信號(hào)中加入隨機(jī)高斯白噪聲：

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

(1)

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

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

(2)

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

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

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

(3)

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

1.2 多尺度排列熵(MPE)

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

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

(4)

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

1)利用EEMD處理GPS高程時(shí)間序列，獲得一系列IMF和1個(gè)余項(xiàng)；

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

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

2 仿真算例分析

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

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

11sin(3πt)+noise

(5)

圖1 仿真信號(hào)Fig.1 Simulation signal

圖2 各階IMF與仿真信號(hào)的相關(guān)系數(shù)Fig.2 Correlation coefficients of each IMF and simulated signal

圖3 各階IMF的值

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

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

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

表1 不同方法的降噪評(píng)價(jià)指標(biāo)

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

3 實(shí)例分析

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

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

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

圖6 各階IMF與原始序列的相關(guān)系數(shù)Fig.6 Correlation coefficients of each IMFand original sequence

圖7 各階IMF的值

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

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

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

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

dnSNR=10 lg(Ps/Pg)

(6)

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

表2 各種方法的降噪效果評(píng)價(jià)指標(biāo)

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

4 結(jié) 語

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