基于EMD的GNSS時間序列降噪軟件實現(xiàn)

賀小星 ，熊常亮 ，常 苗 ，班 亞 ，李姝瑩 ，鄧心茹 ，鄧麗紅

0 引言

全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system，GNSS）已經(jīng)在各行各業(yè)中獲得了廣泛應(yīng)用，近 20余年累積的全球 GNSS連續(xù)運行參考站（continuously operating reference stations，CORS）坐標(biāo)時間序列為大地測量學(xué)及地球動力學(xué)研究提供了寶貴的基礎(chǔ)數(shù)據(jù)[1-3]。GNSS CORS不僅是區(qū)域及全球高精度時空基準(zhǔn)的重要基礎(chǔ)設(shè)施，也是導(dǎo)航與位置服務(wù)、地質(zhì)災(zāi)害監(jiān)測、大地測量學(xué)與地球動力學(xué)等工程和科學(xué)研究的重要支撐[1]。GNSS CORS坐標(biāo)時間序列不僅包含測站的線性變化，也包含非線性變化[4-7]。隨著對大地測量成果所要求的精度越來越高，CORS坐標(biāo)時間序列越來越受到關(guān)注。然而，GNSS時間序列數(shù)據(jù)處理工作非常繁瑣且耗時，且很多重復(fù)性工作。另一方面，GNSS時間序列具有非線性、噪聲成分與信號不易分離等特性[8]，文獻[9]所提出的經(jīng)驗?zāi)B(tài)分解（empirical mode decomposition，EMD）方法處理此類型的信號有一定的優(yōu)越性。文獻[10]將經(jīng)驗?zāi)B(tài)分解結(jié)合小波分解應(yīng)用于 GNSS水汽時間序列，成功分析出不同周期性水汽振蕩的成因。文獻[11]對長周期地心運動時間序列進行 EMD分解，利用重構(gòu)后的時間序列進行分析，提高了各運動方向的周期貢獻率。文獻[12]基于EMD方法精確提取了 GNSS高程時間序列的周期項。綜上所述，EMD方法對于GNSS CORS坐標(biāo)序列非線性變化精密建模具有重要意義，但在實際數(shù)據(jù)處理過程中缺少相應(yīng)的分析軟件。

基于此，為了增強軟件的實用性、交互性，結(jié)合 GNSS時間序列降噪分析的基本原理及分析策略，開發(fā)出基于MATLAB的GNSS時間序列降噪軟件：介紹了基于經(jīng)驗?zāi)B(tài)分解的 GNSS數(shù)據(jù)降噪基本原理，同時也給出了數(shù)據(jù)處理流程；軟件包含多個時間序列降噪分析模塊，交互性較好，可避免繁雜的計算工作，提高 GNSS時間序列數(shù)據(jù)處理及分析的效率。

1 基于EMD的GNSS時間序列降噪原理

1.1 EMD降噪原理

EMD方法應(yīng)用于GNSS時間序列降噪，主要包括信號分解、分界本征模態(tài)分量判定、信號重構(gòu)以及降噪效果評定4部分，降噪流程如圖1所示。EMD降噪的基本思想是將信號分解為若干頻率由高至低的本征模態(tài)函數(shù)（intrinsic mode function，IMF）和1個趨勢項，趨勢項與低頻IMF分量重構(gòu)可實現(xiàn)降噪。

圖1 GNSS坐標(biāo)序列EMD方法降噪流程

通常高頻IMF中包含噪聲信息，而低頻IMF及趨勢項包含信號的特征信息及趨勢。在篩選過程中，每個IMF必須滿足2個條件[12]：①序列中的極值點數(shù)目與過零點數(shù)相差不超過 1個；②由局部極大值構(gòu)成上包絡(luò)線，局部極小值構(gòu)成下包絡(luò)線，2者均值為0。然而，在實際信號分解過程中，第2個條件很難被嚴(yán)格滿足，因此以1個閾值對條件2作替換處理，當(dāng)IMF達到給定閾值時，就認(rèn)為其滿足條件2，閾值表達式為

式中：SD為每個IMF停止篩選的閾值，通常情況下取 0.2~0.3；k為信號分量的序數(shù)；N為序列長度； dk(t)、dk?1( t )為每個IMF篩選過程中2個相鄰的數(shù)據(jù)序列。

分解完成后的原始信號序列可由下式進行表示為

式中： x (t)為重構(gòu)的原始信號序列；m為所有IMF分量的個數(shù)；r (t)為趨勢項。依據(jù)相關(guān)系數(shù)準(zhǔn)則判定分界本征模態(tài)函數(shù)。當(dāng)相關(guān)系數(shù)第1次取極小值時，則其對應(yīng)的IMF為分界IMF，并將其納入至高頻IMF分量進行重構(gòu)。相關(guān)系數(shù)公式及重構(gòu)公式為

式中：ρk為相關(guān)系數(shù)；x?(t)為EMD降噪后的“干凈”的坐標(biāo)序列；K為分界IMF分量的序號； r(t)為趨勢項。

1.2 改進的EMD坐標(biāo)時間序列降噪方法

針對EMD方法存在的模態(tài)混疊問題[13-15]，進一步對EMD方法進行改進。假定原始數(shù)據(jù)序列設(shè)為χ1，運用EMD方法對其進行分解，改進的EMD降噪方法過程如下：

1）初始化：原始數(shù)據(jù)xi，i=1；

2）EMD分解，獲取mi個IMF及趨勢項ri(t)；

3）計算相關(guān)系數(shù)，分界IMF為IMFiK；

7）i=i+1，xi=a；代入第2個步驟。

當(dāng) ki=2時循環(huán)分解停止，此時，將所有低頻IMF分量及趨勢項進行重構(gòu)，得到降噪后的數(shù)據(jù)序列，可表示為

式中：x?1為降噪后的數(shù)據(jù)序列；Ki表示第i個原始數(shù)據(jù)序列的分界IMF下標(biāo)； Mi表示第i個原始數(shù)據(jù)序列分解后 IMF分量的個數(shù)；IMFimi表示第 i個原始數(shù)據(jù)序列的第 mi個IMF分量； rj(t)表示第j次EMD分解的趨勢項。

為定量反映方法的降噪效果，利用相關(guān)系數(shù)、均方根誤差、信噪比等傳統(tǒng)指標(biāo)做出評價。相關(guān)系數(shù)越接近 1表明降噪信號與原始數(shù)據(jù)序列的相似程度越高，均方根誤差值越小說明降噪信號與原始信號之間的偏差程度越小，信噪比值越大說明去噪效果越明顯。均方根誤差(root mean square error，RMSE）的計算公式為

信噪比的計算公式為

式中：RMSE為均方根誤差；Rsn為原始序列與降噪后序列的信噪比。

2 基于EMD的GNSS時間序列降噪軟件設(shè)計及實現(xiàn)

2.1 軟件總體設(shè)計

根據(jù) EMD理論設(shè)計基于 MATLAB的 EMD的GNSS時間序列降噪軟件，軟件總體結(jié)構(gòu)如圖2所示。軟件主要包括數(shù)據(jù)導(dǎo)入、降噪分析與結(jié)果輸出3大模塊，每個模塊下又分為幾個小模塊獨立進行實現(xiàn)。各模塊下設(shè)置若干子菜單，每個子菜單發(fā)揮不同的作用，共同完成單個模塊的功能。不同模塊間相互配合，通過融合與集成的方式，組成1個功能完整、操作簡單的GNSS時間序列降噪軟件。

圖2 軟件結(jié)構(gòu)

2.2 軟件的實現(xiàn)

該軟件充分利用MATLAB在矩陣運算、信號分析、圖形繪制等方面的優(yōu)越性，結(jié)合其圖形用戶界面良好的交互性進行開發(fā)，軟件主界面如圖3所示。

圖3 基于EMD的GNSS時間序列降噪軟件主界面

2.2.1 數(shù)據(jù)導(dǎo)入模塊

在數(shù)據(jù)導(dǎo)入模塊菜單下設(shè)有模擬數(shù)據(jù)與真實數(shù)據(jù)2種選擇，其中真實數(shù)據(jù)可外部進行導(dǎo)入，支持多種數(shù)據(jù)格式，而模擬數(shù)據(jù)可設(shè)置為包括白噪聲和多種有色噪聲模型的信號序列，數(shù)據(jù)導(dǎo)入模塊如圖4所示。真實數(shù)據(jù)模塊支持包括Excel與文本等多種格式的文件，無需其他操作即可完成數(shù)據(jù)的導(dǎo)入，真實數(shù)據(jù)導(dǎo)入界面見圖5。

圖4 數(shù)據(jù)導(dǎo)入模塊

2.2.2 降噪分析模塊

降噪分析模塊是本軟件的核心模塊之一，該模塊包括模擬數(shù)據(jù)分析、真實數(shù)據(jù)分析、結(jié)果輸出3大功能。用戶只需要點擊軟件菜單下對應(yīng)的降噪分析按鈕即可完成對載入數(shù)據(jù)的降噪處理，并自動繪制相關(guān)圖形。鑒于目前針對于 GNSS時間序列降噪的方法多樣性，本軟件將上文介紹的EMD方法與改進的EMD方法各作為1個模塊在軟件中進行實現(xiàn)。 “數(shù)據(jù)分析”依據(jù)前述的EMD、改進的EMD方法對GNSS數(shù)據(jù)進行處理，并根據(jù)相關(guān)系數(shù)、均方根誤差以及信噪比3個指標(biāo)，依據(jù)式（3）、式（6）、式（7）對降噪效果進行評價。

2.2.3 結(jié)果輸出模塊

通過結(jié)果輸出模塊可輸出包括各類圖表在內(nèi)的所有成果，同時為了提高用戶操作的自由度，數(shù)據(jù)存儲的路徑可設(shè)置為計算機上任意路徑，如圖6所示。

3 測試分析

為驗證軟件的有效性與可靠性，對于同一含噪信號，利用軟件中 2種不同的降噪分析方法對其作降噪處理。設(shè)置模擬信號噪聲參數(shù)如下：序列長度 1 200，采樣頻率 1 Hz，高斯白噪聲頻率 4 dB模擬信號可表示為

3.1 EMD降噪模塊測試分析

數(shù)據(jù)導(dǎo)入之后，選擇EMD分析菜單下的模擬分析，軟件即自動對原始信號進行EMD降噪處理，同時繪制相關(guān)系數(shù)曲線如圖 7所示，降噪前后的信號序列如圖8所示。

降噪分析之后，軟件將降噪后序列導(dǎo)出，同時依據(jù)前文所述的 3個評價指標(biāo)對降噪分析的效果進行評估，如表1所示。因原始信號的真值已知，故降噪后信號與含噪聲信號及真實信號分別計算各評價指標(biāo)以定量的反映降噪效果。由表1可知，降噪后的信號序列均方根誤差減小，信噪比與相關(guān)系數(shù)均增大，表明軟件中進行的數(shù)據(jù)處理是行之有效的。

圖5 真實數(shù)據(jù)導(dǎo)入

圖6 選取文件保存路徑

圖7 相關(guān)系數(shù)曲線

圖8 降噪前后信號序列

表1 各評價指標(biāo)值

3.2 改進的EMD降噪模塊測試分析

將同樣的數(shù)據(jù)利用軟件所提供的改進的EMD方法進行降噪分析。默認(rèn)情況下，EMD迭代分解次數(shù)設(shè)為3次，每次EMD分解所得的相關(guān)系數(shù)以及信號序列都由軟件進行記錄，并繪制圖形，迭代分解相關(guān)系數(shù)曲線見圖9，多次EMD分解前后信號序列見圖10。

圖9 迭代分解相關(guān)系數(shù)曲線

圖10 多次EMD分解前后信號序列

每次 EMD分解后重構(gòu)的信號序列都由軟件進行保存，并計算各評價指標(biāo)，以定量反映EMD分解次數(shù)對于降噪的影響，有助于選取降噪效果最好的信號序列。多次EMD分解后的評價指標(biāo)值如表2所示。分析表2可知，進行3次EMD分解所獲得的降噪信號較傳統(tǒng)EMD分解所得降噪信號RMSE的值降低了0.016 7，相關(guān)系數(shù)提高了0.002，信噪比提高了 0.538 dB，說明該方法能夠改善降噪效果，本軟件也實現(xiàn)了數(shù)據(jù)處理與降噪分析的預(yù)期目標(biāo)。

表2 各評價指標(biāo)值

4 結(jié)束語

根據(jù)GNSS時間序列數(shù)據(jù)降噪處理相關(guān)原理，結(jié)合MATLAB開發(fā)了基于EMD的GNSS時間序列降噪軟件，該軟件采用模塊化設(shè)計，功能較完善；在軟件上基于仿真數(shù)據(jù)分別對傳統(tǒng)EMD方法、改進的EMD方法進行實驗分析，并利用均方根誤差、相關(guān)系數(shù)及信噪比等傳統(tǒng)評價指標(biāo)判斷降噪效果。軟件測試實驗結(jié)果表明：改進的EMD方法降噪后的信號序列比之傳統(tǒng)EMD方法降噪的信號序列，RMSE的值降低了0.016 7，相關(guān)系數(shù)提高了0.002，信噪比提高了0.538 dB。該軟件能夠提高 GNSS時間序列數(shù)據(jù)處理及分析的效率，為進一步探索 GNSS基準(zhǔn)站坐標(biāo)序列非線性變化精密建模提供參考。