熊常亮，賀小星，孫喜文，馬下平

（1. 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院，武漢 430077；2. 中國科學(xué)院大學(xué)，北京 100049;3. 華東交通大學(xué) 土木建筑學(xué)院，南昌 330013；4. 陜西鐵路工程職業(yè)技術(shù)學(xué)院，陜西 渭南 714000；5. 西安科技大學(xué) 測繪科學(xué)與技術(shù)學(xué)院，西安 710054）

0 引言

隨著全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system, GNSS）逐步完善，全世界范圍內(nèi)廣泛建立了GNSS 連續(xù)跟蹤觀測站，GNSS 觀測站對(duì)國際地球參考框架的建立與維持發(fā)揮了至關(guān)重要的作用[1-3]。累計(jì)的長周期GNSS 時(shí)間序列，為地殼形變監(jiān)測、全球板塊運(yùn)動(dòng)等地學(xué)研究提供了基礎(chǔ)數(shù)據(jù)，而剔除GNSS 時(shí)間序列中夾雜的噪聲，已成為當(dāng)前研究的熱點(diǎn)[4-7]。針對(duì)GNSS 時(shí)間序列分析過程繁瑣的問題，文獻(xiàn)[8]開發(fā)了格格馬特拉布（GGMatlab）軟件，可對(duì)GNSS 時(shí)間序列進(jìn)行站速度、階躍等進(jìn)行分析與圖形繪制，但無法對(duì)時(shí)間序列進(jìn)行降噪分析。文獻(xiàn)[9]開發(fā)的吉特薩（GITSA），可對(duì)時(shí)間序列進(jìn)行基于多種算法的頻譜分析與回歸分析，但其對(duì)矩陣實(shí)驗(yàn)室（matrix laboratory, MATLAB）外部工具箱具有依賴性。文獻(xiàn)[10]基于派森（Python）語言，開發(fā)了特薩納雷澤爾（TSAnalyzer），但不具備基準(zhǔn)站網(wǎng)共模誤差提取的功能。文獻(xiàn)[11]所開發(fā)的軟件，雖能夠快速實(shí)現(xiàn)GNSS 時(shí)間序列降噪，但在數(shù)據(jù)處理算法與數(shù)據(jù)分析功能上較為單一。

考慮到GNSS 時(shí)間序列非線性分離及殘余噪聲影響，結(jié)合當(dāng)前GNSS 時(shí)間序列分析軟件存在的不足，本文在介紹GNSS 時(shí)間序列噪聲處理理論基礎(chǔ)上，基于MATLAB 平臺(tái)開發(fā)了一款功能完善、界面友好、操作簡單的GNSS 時(shí)間序列分析與降噪軟件。該軟件擁有完整的GNSS 時(shí)間序列降噪與共模誤差分離功能，兼顧數(shù)據(jù)預(yù)處理、統(tǒng)計(jì)分析等擴(kuò)展功能，且無需外部工具箱，為GNSS 時(shí)間序列數(shù)據(jù)處理與分析提供了一個(gè)有力工具。

1 GNSS 時(shí)間序列降噪分析

1.1 EMD 降噪數(shù)學(xué)模型

經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition,EMD）方法是一種自適應(yīng)的信號(hào)處理方法，可用于GNSS 時(shí)間序列噪聲剔除[12-14]。它將原始信號(hào)序列分解為一系列從高頻至低頻的本征模態(tài)函數(shù)（intrinsic mode function, IMF）fIMF。EMD 表達(dá)式為

式中：x(t) 為原始時(shí)間序列；t為歷元數(shù)；k為IMF分量的序號(hào)；m為分解的總IMF 個(gè)數(shù)；r(t) 為趨勢項(xiàng)。

分解完成后，需判定高頻噪聲fIMF分量與低頻信號(hào)fIMF分量的分界值K。常用的判斷方法為相關(guān)系數(shù)法則，該方法計(jì)算每個(gè)fIMF分量與原始時(shí)間序列的相關(guān)系數(shù)，得到對(duì)應(yīng)的相關(guān)系數(shù)矩陣。矩陣中的第一個(gè)極小值對(duì)應(yīng)的fIMF即為分界fIMF，分界fIMF及分界fIMF之前的分量視為高頻噪聲分量，在信號(hào)重構(gòu)時(shí)不再考慮，相關(guān)系數(shù)計(jì)算公式以及信號(hào)重構(gòu)表達(dá)式為

式中：x?(t) 為降噪后時(shí)間序列；kρ為相關(guān)系數(shù)；N為總歷元數(shù)。

獲得重構(gòu)序列后，還可依據(jù)相關(guān)系數(shù)式（2）、均方根誤差（root mean square error, RMSE）SRMSE及信噪比（signal noise ratio , SNR）σSNR對(duì)重構(gòu)序列進(jìn)行精度評(píng)估，均方根誤差sRMSE的計(jì)算公式為

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

若使用真值與重構(gòu)序列計(jì)算上述指標(biāo)，則當(dāng)相關(guān)系數(shù)越接近于 1，重構(gòu)序列與原始序列擬合越好，降噪效果越好；均方根誤差越小，重構(gòu)序列與原始序列偏差越小，降噪效果越好；信噪比越大，噪聲占信號(hào)的比重越小，降噪效果越好。

1.2 空間濾波數(shù)學(xué)模型

連續(xù)GNSS 時(shí)間序列中存在一種空間相關(guān)的誤差，稱之為共模誤差（common mode error,CME）。CME 被普遍認(rèn)為是GNSS 時(shí)間時(shí)間序列的主要來誤差源之一，CME 的存在，極大地影響GNSS 解的精度與可靠性，尤其是對(duì)測站速度估計(jì)的影響較大。文獻(xiàn)[15]于1997 年提出了一種依據(jù)全球定位系統(tǒng)（global positioning system, GPS）單日解殘差時(shí)間序列進(jìn)行CME 計(jì)算的方法，稱之為堆棧濾波法。

在堆棧濾波法中，假設(shè)基準(zhǔn)站網(wǎng)中有S個(gè)站，在剔除GNSS 時(shí)間序列的均值與趨勢項(xiàng)后，得到其殘差時(shí)間序列，計(jì)算CME 的公式為

式中：εk,t為第k個(gè)基準(zhǔn)站第t個(gè)歷元的殘差時(shí)間序列；εSF(t)為堆棧濾波算得的第t個(gè)歷元的CME 值。

原始?xì)埐顣r(shí)間序列減去CME 值，即為濾波后的殘差時(shí)間序列。通常以濾波前后時(shí)間序列的均方根（root mean square, RMS）SRMS值評(píng)估時(shí)間序列濾波后的穩(wěn)定性和精確度。SRMS的計(jì)算公式為

2 軟件模塊及框架設(shè)計(jì)

2.1 基于EMD 的GNSS 時(shí)間序列降噪模塊設(shè)計(jì)

軟件中首先實(shí)現(xiàn)了基于EMD 的GNSS 時(shí)間序列降噪模型。此外，針對(duì)EMD 方法中仍存在著模態(tài)混疊、分界fIMF值判斷不準(zhǔn)等問題，影響GNSS 時(shí)間序列的精確降噪，實(shí)現(xiàn)了三種改進(jìn)的EMD 方法：

1）削弱混疊法。該方法每次僅分離出最高頻的fIMF，通過重復(fù)利用EMD 方法不斷從剩余的高頻fIMF分量提取低頻分量，對(duì)于模態(tài)混疊問題可起到一定的抑制作用[16]。

2）周期密度法。該方法基于平均周期與能量密度的乘積為常數(shù)的特點(diǎn)，推導(dǎo)出一種能夠自適應(yīng)確定分界fIMF值的算法，對(duì)于分界fIMF值的判斷可取得精確的結(jié)果[17]。

3）復(fù)合指標(biāo)法。該方法綜合考慮時(shí)間序列的平滑度與SRMSE計(jì)算出一張復(fù)合指標(biāo)，用于分界fIMF值的判斷，結(jié)果具有可靠性[18]。

此外， 軟件中還整合了經(jīng)驗(yàn)?zāi)B(tài)分解（ensemble empirical mode decomposition, EEMD）方法，并將EEMD 作為核心算法之一[19]。通過上述五種基于EMD 的降噪方法搭建了軟件的降噪分析模塊。

2.2 GNSS 時(shí)間序列共模誤差分離模塊設(shè)計(jì)

基于對(duì)堆棧濾波法理論的分析易得，該方法所計(jì)算的CME 值是GNSS 基準(zhǔn)站網(wǎng)內(nèi)所有站點(diǎn)在該歷元的殘差時(shí)間序列的平均值。因此，當(dāng)基準(zhǔn)站網(wǎng)的規(guī)模過大時(shí)，空間分布顯然不平均，此時(shí)利用堆棧濾波法計(jì)算CME 值就失去了意義。

基于此，文獻(xiàn)[20]將單日解誤差引入堆棧濾波法，提出了加權(quán)堆棧濾波法；文獻(xiàn)[21]將站點(diǎn)間相關(guān)性進(jìn)行量化計(jì)算出相關(guān)系數(shù)再與堆棧濾波法相結(jié)合，提出了相關(guān)加權(quán)疊加濾波法。文獻(xiàn)[22]采用主成分分析法（principal component analysis，PCA）與卡爾胡寧-勒夫（Karhunen-Loeve）變換相結(jié)合的方法進(jìn)行CME 提取，有效抑制了部分測站的地區(qū)效應(yīng)。

綜上所述，軟件中的空間濾波分析模塊將上述四種濾波算法全部實(shí)現(xiàn)，使得該模塊擁有解算各種尺度的GNSS 基準(zhǔn)站網(wǎng)的能力。

2.3 軟件框架設(shè)計(jì)

根據(jù)前述的GNSS 模塊，結(jié)合MATLAB 圖形用戶界面進(jìn)行軟件開發(fā)。軟件設(shè)計(jì)包括數(shù)據(jù)預(yù)處理、空間濾波、降噪分析、擴(kuò)展功能等模塊，軟件結(jié)構(gòu)見圖1。

圖1 軟件總體結(jié)構(gòu)

為了提高軟件的可操作性，各模塊單獨(dú)設(shè)計(jì)界面，并將相關(guān)源程序、輸出數(shù)據(jù)等保存到對(duì)應(yīng)文件夾，既方便軟件后期的維護(hù)與擴(kuò)展同時(shí)兼顧軟件使用時(shí)的清晰簡潔。通過批處理功能提高數(shù)據(jù)處理效率。此外，計(jì)算過程中的參數(shù)、結(jié)果與圖形都對(duì)應(yīng)有獨(dú)立的函數(shù)，進(jìn)一步提高了軟件的兼容性，便于更新與升級(jí)。

2.3.1 數(shù)據(jù)預(yù)處理模塊

數(shù)據(jù)預(yù)處理模塊主要是對(duì)GNSS 時(shí)間序列數(shù)據(jù)進(jìn)行必要的處理，主要包括GNSS 時(shí)間序列的數(shù)據(jù)質(zhì)量分析（格式編輯、時(shí)間序列可視化、SRMS值計(jì)算、站間相關(guān)性分析等）。除此之外還包含粗差剔除功能，可對(duì)GNSS 時(shí)間序列進(jìn)行粗差探測與剔除。具體功能如下：

1）格式編輯模塊。支持將中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)及斯克里普斯軌道和常駐陣列中心（Scripps orbit and permanent array center, SOPAC）所發(fā)布的兩種*.neu 文件批量轉(zhuǎn)換為軟件中支持的*.mom 文件，也可將起止時(shí)間點(diǎn)設(shè)為參數(shù)、批量截取時(shí)間序列為給定年限的文件。

2）粗差剔除模塊。包含3Sigma、5Sigma、3IQR以及MAD 四種粗差剔除算法[23]，在完成粗差剔除后，軟件將自動(dòng)繪制圖形、標(biāo)定粗差點(diǎn)位置，并給出粗差剔除率報(bào)告。

3）參數(shù)計(jì)算模塊?？捎?jì)算SRMS值與相關(guān)系數(shù)，SRMS值反映時(shí)間序列的穩(wěn)定性的重要指標(biāo)，見式（7）。相關(guān)系數(shù)則反映不同站點(diǎn)時(shí)間序列間相關(guān)性的指標(biāo)，軟件中設(shè)計(jì)有相應(yīng)模塊來計(jì)算該指標(biāo)，其計(jì)算表達(dá)式為

計(jì)算完成后，相關(guān)系數(shù)矩陣將以文本形式進(jìn)行保存。

2.3.2 降噪分析模塊

降噪分析模塊作為軟件的核心模塊之一，內(nèi)置的算法為上述的EMD、EEMD 以及最新的三種基于EMD 的改進(jìn)算法。該模塊可實(shí)現(xiàn)完整的GNSS時(shí)間序列降噪流程，包括數(shù)據(jù)導(dǎo)入與仿真、基于EMD 的時(shí)間序列降噪以及精度評(píng)估等。其中數(shù)據(jù)仿真方面支持周期性信號(hào)、時(shí)變季節(jié)性信號(hào)以及多種噪聲模型。同時(shí)可繪制一系列與降噪分析相關(guān)的圖形，如EMD 分解后的相關(guān)系數(shù)值、降噪前后的時(shí)間序列對(duì)比圖、EMD/EEMD 分解后的每個(gè)fIMF等。降噪分析模塊也將輸出計(jì)算過程中的重要參數(shù)供算法研究。除此之外，擁有基于快速傅里葉變換的頻譜分析功能，可計(jì)算GNSS 時(shí)間序列的功率譜密度，評(píng)估頻域中信號(hào)的主導(dǎo)周期。

2.3.3 空間濾波分析模塊

空間濾波分析模塊為GNSS 基準(zhǔn)站網(wǎng)時(shí)間序列空間濾波而設(shè)計(jì)，擁有完善的濾波分析功能。核心算法為堆棧濾波法、加權(quán)堆棧濾波法、相關(guān)加權(quán)疊加濾波法以及PCA 方法，因此可適用于各種規(guī)模的GNSS 基準(zhǔn)站網(wǎng)解算。該模塊擁有數(shù)據(jù)導(dǎo)入導(dǎo)出、站點(diǎn)SRMS值計(jì)算、CME 提取及圖形繪制等功能?？煽焖贉?zhǔn)確地實(shí)現(xiàn)數(shù)據(jù)可視化，輸出如濾波前后時(shí)間序列、CME、SRMS等一系列圖形與表格。

2.3.4 擴(kuò)展功能模塊

軟件的擴(kuò)展功能模塊用于增強(qiáng)GNSS 時(shí)間序列的表達(dá)與具體應(yīng)用，分為站點(diǎn)搜索、階躍分析以及統(tǒng)計(jì)分析三部分功能：

1）站點(diǎn)檢索?；谟脩糇越℅NSS 站點(diǎn)庫，根據(jù)中心站經(jīng)緯度或點(diǎn)名，搜尋附近站點(diǎn)并輸出距離。再者，GNSS 與衛(wèi)星激光測距等空間觀測技術(shù)的并置站時(shí)間序列被廣泛運(yùn)用[24-25]?；诖耍军c(diǎn)搜索模塊中設(shè)計(jì)有GNSS 站與驗(yàn)潮站并置站搜尋功能，可快速輸出站點(diǎn)數(shù)據(jù)庫中的并置站。

2）階躍分析?？蓪?duì)仿真時(shí)間序列添加階躍，也可以依據(jù)數(shù)據(jù)中心發(fā)布的參數(shù)，對(duì)含有階躍的實(shí)測時(shí)間序列進(jìn)行修正，提高GNSS 時(shí)間序列的精確度。

3）統(tǒng)計(jì)分析。運(yùn)用箱型圖與小提琴圖的算法，對(duì)GNSS 時(shí)間序列進(jìn)行可視化表達(dá)，直觀地反映數(shù)據(jù)的分布情況。

3 軟件實(shí)現(xiàn)與算例分析

3.1 軟件界面

在MATLAB 的圖形用戶界面面板，分別給每個(gè)功能模塊建立獨(dú)立界面，再通過接口與主界面相連，各模塊界面如圖2 所示。

3.2 算例分析

3.2.1 時(shí)間序列降噪分析

選用一組仿真時(shí)間序列，在本軟件中進(jìn)行時(shí)間序列的降噪分析。降噪方法選用前述的五種方法，以綜合比對(duì)不同方法的降噪差異。因仿真數(shù)據(jù)真值已知，可基于三個(gè)傳統(tǒng)評(píng)價(jià)指標(biāo)定量反映各方法的降噪效果以及降噪分析模塊的穩(wěn)定性與可靠性。仿真數(shù)據(jù)表達(dá)式為

式中：1f、f2、f3為三個(gè)周期分量；δ為6 dB 的白噪聲。采樣點(diǎn)設(shè)置為2 000 個(gè)，各個(gè)分量及噪聲的圖形見圖3。

圖3 仿真數(shù)據(jù)各分量及噪聲

采用軟件中的統(tǒng)計(jì)分析功能對(duì)數(shù)據(jù)進(jìn)行表達(dá)。圖4 為基于箱型圖與小提琴圖的統(tǒng)計(jì)結(jié)果，圖中直觀地表達(dá)了時(shí)間序列的典型代表值、異常值及數(shù)據(jù)整體分布。

圖4 仿真信號(hào)的統(tǒng)計(jì)分析結(jié)果

除統(tǒng)計(jì)分析方法外，頻譜分析也是GNSS 時(shí)間序列常用的分析方法之一。軟件的頻譜分析功能可利用快速傅里葉變換計(jì)算仿真信號(hào)的功率譜密度，如圖5 所示。

圖5 仿真信號(hào)的功率譜密度

在降噪分析模塊中，使用五種降噪方法，對(duì)信號(hào)作降噪分析。其中，EEMD 的整合數(shù)設(shè)為300，噪聲與信號(hào)的標(biāo)準(zhǔn)差之比設(shè)定為0.2[19]。在全部計(jì)算完成后，從精度評(píng)估報(bào)告中讀取其相關(guān)系數(shù)、均方根誤差及信噪比三個(gè)評(píng)價(jià)指標(biāo)值，如表1 所示。

表1 降噪后信號(hào)與真實(shí)信號(hào)的各評(píng)價(jià)指標(biāo)

因仿真信號(hào)真值已知，故表1 中各指標(biāo)值是由降噪后數(shù)據(jù)與真實(shí)信號(hào)計(jì)算而得。相關(guān)系數(shù)越大，則降噪后序列與原始序列擬合越好；sRMSE值越小，則信號(hào)偏差越小；σSNR值越大，則降噪效果越好?；诖耍瑢?duì)表1 分析可得，五種方法對(duì)于序列降噪都取得了優(yōu)異的效果。其中，效果最好的是周期密度法與復(fù)合指標(biāo)法，二者精確地判定了分界fIMF值，精度指標(biāo)值最高，其次是EEMD 法、削弱混疊法及EMD 法。

除使用各評(píng)價(jià)指標(biāo)進(jìn)行定量反映降噪效果外，軟件繪制相關(guān)圖形，直觀反映各方法的降噪效果，見圖6。其中，該實(shí)驗(yàn)中周期密度法與復(fù)合指標(biāo)法的結(jié)果一致，故僅列周期密度法的時(shí)間序列對(duì)比圖。上述結(jié)果驗(yàn)證了降噪分析模塊對(duì)于GNSS 時(shí)間序列降噪分析的準(zhǔn)確性。

圖6 不同方法降噪前后信號(hào)序列對(duì)比

3.2.2 基準(zhǔn)站共模誤差分離

以美國南加州的一個(gè)小區(qū)域GNSS 基準(zhǔn)站網(wǎng)，時(shí)間序列跨度為2005—2011 年的數(shù)據(jù)進(jìn)行基準(zhǔn)站共模誤差分離，以測試軟件空間濾波模塊的可靠性，站點(diǎn)名稱及具體經(jīng)緯度如表2 所示。

表2 站點(diǎn)經(jīng)緯度

對(duì)GNSS 時(shí)間序列進(jìn)行去均值、趨勢項(xiàng)及粗差處理得到殘差時(shí)間序列，之后對(duì)缺失歷元進(jìn)行插值補(bǔ)全。最后采用前文所述中的四種誤差分離方法計(jì)算出其中的共模誤差，用原始?xì)埐钚蛄袦p去共模誤差得到濾波后的殘余時(shí)間序列。在軟件中選擇堆棧濾波、加權(quán)堆棧濾波、相關(guān)加權(quán)疊加濾波以及PCA 四種方法對(duì)原始序列殘差時(shí)間序列進(jìn)行CME 提取。不同方法進(jìn)行共模誤差提取前后的殘差時(shí)間序列SRMS值如圖7 至圖9 所示。

圖7 東（E）方向CME 提取前后的 SRMS值

圖8 北（N）方向CME 提取前后的 SRMS值

圖9 天（U）方向CME 提取前后的 SRMS值

在CME 提取后，8 個(gè)站點(diǎn)在E、N、U方向的SRMS值都有明顯降低，同在一個(gè)數(shù)量級(jí)，表明時(shí)間序列的穩(wěn)定性得到了提高。特別是經(jīng)PCA 方法處理后的時(shí)間序列，取得了精度更高的結(jié)果，而其他方法的效果整體相當(dāng)。限于篇幅，僅在圖10 至圖13 列出KTBW 站經(jīng)四種不同濾波方法處理后的效果圖，實(shí)線為原始?xì)埐顣r(shí)間序列，點(diǎn)線為濾波后殘差時(shí)間序列，虛線為對(duì)應(yīng)的CME 值。

圖10 KTBW 站堆棧濾波前后時(shí)間序列及CME 值

圖11 KTBW 站加權(quán)堆棧濾波前后時(shí)間序列及CME 值

圖12 KTBW 站相關(guān)加權(quán)疊加濾波前后時(shí)間序列及CME 值

圖13 KTBW 站PCA 濾波前后時(shí)間序列及CME 值

4 結(jié)束語

根據(jù) GNSS 時(shí)間序列數(shù)據(jù)降噪處理相關(guān)原理，開發(fā)了基于MATLAB 的GNSS 時(shí)間序列分析與降噪軟件。該軟件考慮GNSS 時(shí)間序列分析與降噪的功能需求，充分利用MATLAB 在科學(xué)計(jì)算與數(shù)據(jù)可視化上的優(yōu)勢，實(shí)現(xiàn)了 EMD 與EEMD 算法模型，并實(shí)現(xiàn)了三種基于EMD 的改進(jìn)算法，聯(lián)合時(shí)間序列仿真、輸入輸出、精度評(píng)估等功能一起組成降噪分析模塊。另一方面，將堆棧濾波、加權(quán)堆棧濾波、相關(guān)加權(quán)疊加濾波以及PCA 濾波法與數(shù)據(jù)批處理功能結(jié)合，組成空間濾波分析模塊。最后，添加各類數(shù)據(jù)預(yù)處理及統(tǒng)計(jì)分析功能，共同組成功能完善、結(jié)構(gòu)完整的GNSS 時(shí)間序列分析與降噪軟件。根據(jù)仿真數(shù)據(jù)與實(shí)測數(shù)據(jù)測試結(jié)果表明：1）軟件操作簡單，具有良好的穩(wěn)定性；2）所實(shí)現(xiàn)的三種改進(jìn)算法在傳統(tǒng)EMD 方法的基礎(chǔ)上得出了更精確的GNSS時(shí)間序列；3）空間濾波分析模塊快速實(shí)現(xiàn)GNSS時(shí)間序列的CME 提取，以可視化的形式給出數(shù)據(jù)處理結(jié)果，極大地提高了GNSS 時(shí)間序列分析的效率。