GNSS 時間序列分析與降噪軟件的實現

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

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

0 引言

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

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

1 GNSS 時間序列降噪分析

1.1 EMD 降噪數學模型

經驗模態分解（empirical mode decomposition,EMD）方法是一種自適應的信號處理方法，可用于GNSS 時間序列噪聲剔除[12-14]。它將原始信號序列分解為一系列從高頻至低頻的本征模態函數（intrinsic mode function, IMF）fIMF。EMD 表達式為

式中：x(t) 為原始時間序列；t為歷元數；k為IMF分量的序號；m為分解的總IMF 個數；r(t) 為趨勢項。

分解完成后，需判定高頻噪聲fIMF分量與低頻信號fIMF分量的分界值K。常用的判斷方法為相關系數法則，該方法計算每個fIMF分量與原始時間序列的相關系數，得到對應的相關系數矩陣。矩陣中的第一個極小值對應的fIMF即為分界fIMF，分界fIMF及分界fIMF之前的分量視為高頻噪聲分量，在信號重構時不再考慮，相關系數計算公式以及信號重構表達式為

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

獲得重構序列后，還可依據相關系數式（2）、均方根誤差（root mean square error, RMSE）SRMSE及信噪比（signal noise ratio , SNR）σSNR對重構序列進行精度評估，均方根誤差sRMSE的計算公式為

信噪比的計算公式為

若使用真值與重構序列計算上述指標，則當相關系數越接近于 1，重構序列與原始序列擬合越好，降噪效果越好；均方根誤差越小，重構序列與原始序列偏差越小，降噪效果越好；信噪比越大，噪聲占信號的比重越小，降噪效果越好。

1.2 空間濾波數學模型

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

在堆棧濾波法中，假設基準站網中有S個站，在剔除GNSS 時間序列的均值與趨勢項后，得到其殘差時間序列，計算CME 的公式為

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

原始殘差時間序列減去CME 值，即為濾波后的殘差時間序列。通常以濾波前后時間序列的均方根（root mean square, RMS）SRMS值評估時間序列濾波后的穩定性和精確度。SRMS的計算公式為

2 軟件模塊及框架設計

2.1 基于EMD 的GNSS 時間序列降噪模塊設計

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

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

2）周期密度法。該方法基于平均周期與能量密度的乘積為常數的特點，推導出一種能夠自適應確定分界fIMF值的算法，對于分界fIMF值的判斷可取得精確的結果[17]。

3）復合指標法。該方法綜合考慮時間序列的平滑度與SRMSE計算出一張復合指標，用于分界fIMF值的判斷，結果具有可靠性[18]。

此外， 軟件中還整合了經驗模態分解（ensemble empirical mode decomposition, EEMD）方法，并將EEMD 作為核心算法之一[19]。通過上述五種基于EMD 的降噪方法搭建了軟件的降噪分析模塊。

2.2 GNSS 時間序列共模誤差分離模塊設計

基于對堆棧濾波法理論的分析易得，該方法所計算的CME 值是GNSS 基準站網內所有站點在該歷元的殘差時間序列的平均值。因此，當基準站網的規模過大時，空間分布顯然不平均，此時利用堆棧濾波法計算CME 值就失去了意義。

基于此，文獻[20]將單日解誤差引入堆棧濾波法，提出了加權堆棧濾波法；文獻[21]將站點間相關性進行量化計算出相關系數再與堆棧濾波法相結合，提出了相關加權疊加濾波法。文獻[22]采用主成分分析法（principal component analysis，PCA）與卡爾胡寧-勒夫（Karhunen-Loeve）變換相結合的方法進行CME 提取，有效抑制了部分測站的地區效應。

綜上所述，軟件中的空間濾波分析模塊將上述四種濾波算法全部實現，使得該模塊擁有解算各種尺度的GNSS 基準站網的能力。

2.3 軟件框架設計

根據前述的GNSS 模塊，結合MATLAB 圖形用戶界面進行軟件開發。軟件設計包括數據預處理、空間濾波、降噪分析、擴展功能等模塊，軟件結構見圖1。

圖1 軟件總體結構

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

2.3.1 數據預處理模塊

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

1）格式編輯模塊。支持將中國大陸構造環境監測網絡及斯克里普斯軌道和常駐陣列中心（Scripps orbit and permanent array center, SOPAC）所發布的兩種*.neu 文件批量轉換為軟件中支持的*.mom 文件，也可將起止時間點設為參數、批量截取時間序列為給定年限的文件。

2）粗差剔除模塊。包含3Sigma、5Sigma、3IQR以及MAD 四種粗差剔除算法[23]，在完成粗差剔除后，軟件將自動繪制圖形、標定粗差點位置，并給出粗差剔除率報告。

3）參數計算模塊。可計算SRMS值與相關系數，SRMS值反映時間序列的穩定性的重要指標，見式（7）。相關系數則反映不同站點時間序列間相關性的指標，軟件中設計有相應模塊來計算該指標，其計算表達式為

計算完成后，相關系數矩陣將以文本形式進行保存。

2.3.2 降噪分析模塊

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

2.3.3 空間濾波分析模塊

空間濾波分析模塊為GNSS 基準站網時間序列空間濾波而設計，擁有完善的濾波分析功能。核心算法為堆棧濾波法、加權堆棧濾波法、相關加權疊加濾波法以及PCA 方法，因此可適用于各種規模的GNSS 基準站網解算。該模塊擁有數據導入導出、站點SRMS值計算、CME 提取及圖形繪制等功能。可快速準確地實現數據可視化，輸出如濾波前后時間序列、CME、SRMS等一系列圖形與表格。

2.3.4 擴展功能模塊

軟件的擴展功能模塊用于增強GNSS 時間序列的表達與具體應用，分為站點搜索、階躍分析以及統計分析三部分功能：

1）站點檢索。基于用戶自建GNSS 站點庫，根據中心站經緯度或點名，搜尋附近站點并輸出距離。再者，GNSS 與衛星激光測距等空間觀測技術的并置站時間序列被廣泛運用[24-25]。基于此，站點搜索模塊中設計有GNSS 站與驗潮站并置站搜尋功能，可快速輸出站點數據庫中的并置站。

2）階躍分析。可對仿真時間序列添加階躍，也可以依據數據中心發布的參數，對含有階躍的實測時間序列進行修正，提高GNSS 時間序列的精確度。

3）統計分析。運用箱型圖與小提琴圖的算法，對GNSS 時間序列進行可視化表達，直觀地反映數據的分布情況。

3 軟件實現與算例分析

3.1 軟件界面

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

3.2 算例分析

3.2.1 時間序列降噪分析

選用一組仿真時間序列，在本軟件中進行時間序列的降噪分析。降噪方法選用前述的五種方法，以綜合比對不同方法的降噪差異。因仿真數據真值已知，可基于三個傳統評價指標定量反映各方法的降噪效果以及降噪分析模塊的穩定性與可靠性。仿真數據表達式為

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

圖3 仿真數據各分量及噪聲

采用軟件中的統計分析功能對數據進行表達。圖4 為基于箱型圖與小提琴圖的統計結果，圖中直觀地表達了時間序列的典型代表值、異常值及數據整體分布。

圖4 仿真信號的統計分析結果

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

圖5 仿真信號的功率譜密度

在降噪分析模塊中，使用五種降噪方法，對信號作降噪分析。其中，EEMD 的整合數設為300，噪聲與信號的標準差之比設定為0.2[19]。在全部計算完成后，從精度評估報告中讀取其相關系數、均方根誤差及信噪比三個評價指標值，如表1 所示。

表1 降噪后信號與真實信號的各評價指標

因仿真信號真值已知，故表1 中各指標值是由降噪后數據與真實信號計算而得。相關系數越大，則降噪后序列與原始序列擬合越好；sRMSE值越小，則信號偏差越小；σSNR值越大，則降噪效果越好。基于此，對表1 分析可得，五種方法對于序列降噪都取得了優異的效果。其中，效果最好的是周期密度法與復合指標法，二者精確地判定了分界fIMF值，精度指標值最高，其次是EEMD 法、削弱混疊法及EMD 法。

除使用各評價指標進行定量反映降噪效果外，軟件繪制相關圖形，直觀反映各方法的降噪效果，見圖6。其中，該實驗中周期密度法與復合指標法的結果一致，故僅列周期密度法的時間序列對比圖。上述結果驗證了降噪分析模塊對于GNSS 時間序列降噪分析的準確性。

圖6 不同方法降噪前后信號序列對比

3.2.2 基準站共模誤差分離

以美國南加州的一個小區域GNSS 基準站網，時間序列跨度為2005—2011 年的數據進行基準站共模誤差分離，以測試軟件空間濾波模塊的可靠性，站點名稱及具體經緯度如表2 所示。

表2 站點經緯度

對GNSS 時間序列進行去均值、趨勢項及粗差處理得到殘差時間序列，之后對缺失歷元進行插值補全。最后采用前文所述中的四種誤差分離方法計算出其中的共模誤差，用原始殘差序列減去共模誤差得到濾波后的殘余時間序列。在軟件中選擇堆棧濾波、加權堆棧濾波、相關加權疊加濾波以及PCA 四種方法對原始序列殘差時間序列進行CME 提取。不同方法進行共模誤差提取前后的殘差時間序列SRMS值如圖7 至圖9 所示。

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

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

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

在CME 提取后，8 個站點在E、N、U方向的SRMS值都有明顯降低，同在一個數量級，表明時間序列的穩定性得到了提高。特別是經PCA 方法處理后的時間序列，取得了精度更高的結果，而其他方法的效果整體相當。限于篇幅，僅在圖10 至圖13 列出KTBW 站經四種不同濾波方法處理后的效果圖，實線為原始殘差時間序列，點線為濾波后殘差時間序列，虛線為對應的CME 值。

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

圖11 KTBW 站加權堆棧濾波前后時間序列及CME 值

圖12 KTBW 站相關加權疊加濾波前后時間序列及CME 值

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

4 結束語

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