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

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

0 引言

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

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

1 基于EMD的GNSS時(shí)間序列降噪原理

1.1 EMD降噪原理

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

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

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

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

分解完成后的原始信號(hào)序列可由下式進(jìn)行表示為

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

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

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

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

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

2）EMD分解，獲取mi個(gè)IMF及趨勢(shì)項(xiàng)ri(t)；

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

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

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

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

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

信噪比的計(jì)算公式為

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

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

2.1 軟件總體設(shè)計(jì)

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

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

2.2 軟件的實(shí)現(xiàn)

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

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

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

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

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

2.2.2 降噪分析模塊

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

2.2.3 結(jié)果輸出模塊

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

3 測(cè)試分析

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

3.1 EMD降噪模塊測(cè)試分析

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

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

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

圖6 選取文件保存路徑

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

圖8 降噪前后信號(hào)序列

表1 各評(píng)價(jià)指標(biāo)值

3.2 改進(jìn)的EMD降噪模塊測(cè)試分析

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

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

圖10 多次EMD分解前后信號(hào)序列

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

表2 各評(píng)價(jià)指標(biāo)值

4 結(jié)束語(yǔ)

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