區域GNSS基準網解算及其時間序列分析

尹 彤，魯欣宇，聶建亮，蒲永峰，武慧琳

（1.自然資源部測繪標準化研究所，西安 710054；2.河南省遙感測繪院，鄭州 450003；3.陜西測繪地理信息局，西安 710054；4.自然資源部第二地形測量隊，西安 710054）

0 引言

隨著高精度全球衛星導航系統（global navigation satellite system，GNSS）網絡的不斷建設和完善，至今，在GNSS基準站上已經積累了二十多年的觀測數據，這些數據成果將形成長期GNSS坐標時間序列，其中蘊含著很多地表變化及地球框架變化的周期和非周期的運動信息。如何實現區域 GNSS網的聯合解算，從而獲取長期連續的坐標時間序列，是進行基準站穩定性分析的前提。在連續的GNSS臺站的時間序列中，包含地殼、地表運動信息，主要分為構造信息和非構造信息。對非構造運動信息需要進行分類剔除。如海潮運動、天線相位中心、地表大氣負荷、地表水遷移等，才能精確獲取純凈的時間序列，從而分析基準站的穩定性。

目前國內 GNSS網解算都停留在綜合平差階段，沒有考慮長期時間序列噪聲和誤差分析過程，故成果精度不高。本文探討使用區域GNSS數據，采用加米特（GAMIT）/格洛布克（GLOBK）軟件進行區域聯合解算，獲取基準站的坐標時間序列，然后利用主成分分析的方法提取 GNSS時序中的共模噪聲，并在時序分析中將其剔除，繼而基于純凈的 GNSS時間序列，進一步提高基準坐標和運行速度的解算精度。

1 基準網解算策略及時間序列分析方法

1.1 區域GNSS基準網解算策略

采用GAMIT/GLOBK軟件進行區域GNSS基準網解算。GAMIT/GLOBK軟件是一款精密的GNSS定位和定軌軟件。它由美國麻省理工學院（Massachusetts Institute of Technology, MIT）和加州大學圣地亞哥分校的美國斯克里普斯海洋研究所（Scripps Institution of Oceanography, SIO）共同開發和維護。

GAMIT處理雙差觀測量，利用雙差直接消除衛星及測站的時間誤差，采用最小二乘原理進行參數估計進行平差計算。軟件適合進行上千千米的超長基線解算，點位精度可以達到毫米級別，是全球公認的高精度 GNSS長基線處理軟件。GAMIT軟件可以實現超長基線的GNSS多星座數據混合解算，它由許多不同功能的模塊組成, 這些模塊可以獨立地運行。

GLOBK是一個卡爾曼濾波器，其主要目的是綜合空間大地測量和經典大地測量的初步處理結果完成數據的后處理。它的輸入數據一般是一些準觀測量，如測站坐標、地球自轉參數、衛星軌道及它們的方差和協方差。GLOBK可以融合多種基線成果格式，實現多源基線成果的聯合平差。GLOBK軟件的主要功能有：1）實現多源基線成果和多天綜合基線成果的聯合平差；2）獲取 GNSS跟蹤站的坐標分量時間序列、計算站點運動速度場。GLOBK平差計算主要分為：1）輸入文件的格式轉換；2）平差和輸出命令文件準備；3）平差命令globk或glred、成果輸出提取。通過GLOBK軟件可以獲取多時段的GNSS站點的坐標時間序列。

1.2 GNSS時間序列分析方法

GNSS時間序列分析目的就是進一步改正和過濾 GNSS數據解算過程中未考慮的噪聲，獲取更精確的基站位置和運動狀態（位置和速度）。共模誤差（common mode error，CME）是GNSS 坐標時間序列呈現周期性的主要原因之一，需對CME進行改正，從而減少臺站坐標時間序列坐標分量的上下波動，提高時間序列的信噪比，減弱坐標序列周期信號。主成分分析方法（principal components analysis，PCA）的應用，能夠識別出規模較大的區域網絡的 CME，作為第一主成分或主要成分的線性組合，并對單個站點進行“區域調整”。主成分方法可以提高CME提取的可靠性，并顧及形式坐標誤差。

主成分分析是一種通過線性變換來簡化數據集的方法，它將數據變換到一個新的坐標系統中，使得數據投影的第一大方差在第一個坐標（第一主成分）上，第二大方差在第二個坐標（第二主成分）上，依次類推。主成分分析經常用減少數據集的維數，同時保持數據集的對方差貢獻最大的特征。旨在利用降維思想，把多指標轉化為少數幾個綜合指標。在GNSS基準站的長期時間序列中，運用PCA主成分分析可以把隨時間變化的序列分解成時間域的主分量和空間域的特征向量。主成分分析模型由 GNSS坐標向量殘差時間序列構造觀測矩陣R(t,r) （= 1, 2, ···,；= 1, 2, ···,）：t表示時間域中個時間元素的第個；r表示空間域里分布的個研究對象中的第個，其中坐標序列需首先經過內插實現連續化（用個站組成的 GNSS站網構造觀測矩陣）,同時分別為每個地心分量（北、東和高）創建分量矩陣。基于觀測矩陣的元素協方差的每個元素的對稱矩陣B(t,r)定義為

主成分分析方法通過對協方差矩陣(,r)進行特征值分解，分解為和（為正交矩陣，為具有個非零對角特征值{}(≥)的對角矩陣），其計算方法為

那么，采用特征值和特征向量矩陣中對應列的降階來定義各階主成分為

式中：(t)為(t,r)的次主成分；v為第個特征值對應的特征向量。

PCA是去除CME的有效算法，所得到的正交分量（主分量）按變化幅度從大到小排序。PCA方法的基礎是選擇主分量，主分量超過一半的對應歸一化響應的值大于25%。CME的計算方法為

式中：t為時間域的個時間元素；r為空間域分布的第個研究對象；(t)為式（3）中(,r)的次主成分；v為第個特征值對應的特征向量；為前個主成分。

這種方法通過正交變換重建濾波殘差，從而獲取 GNSS殘差序列的至少 6個主成分，其比重依次降低。進而獲取 GNSS坐標序列殘差中的共模誤差。PCA方法通過減去占數據方差中的主要分量來提取共模誤差，然后利用最大似然估計對其進行分析。這樣不僅數據的均方根值會減小，隨機部分也更接近白噪聲。

2 區域GNSS網綜合平差及其時間序列分析

2.1 數據收集及整理

本文采用中國境內及周邊的國際 GNSS服務（International GNSS Service, IGS）站的原始觀測數據作為解算數據，收集了14個IGS站從2013—2019年至少 6 a的數據。基線解算軟件采用GAMIT/GLOBK10.6，時間序列分析軟件采用準觀測組合分析（quasi-observation combination analysis,QOCA）軟件的姆洛德（Mload） 模塊和 PCA模塊。采用 IGS站的名稱及時間跨度如表1所示，點位分布如圖1所示。

表1 IGS站點名稱及時間跨度表

圖1 14個IGS站點位分布圖

本文收集的其他改正模型資料包括全球大氣模型（Atmospheric-M）、全球積雪模型（Snow-M）、全球土壤濕度模型（Soil-M）、非潮汐改正模型（Ecco-M）等的數據，時間跨度為2013—2019年，如表2所示。

表2 改正模型名稱及時間跨度度統計表

從表2的時間可以判斷，時間序列中如果需要加入 Mload的各種改正的話，必須選用 2 013.0~2 019.0 a的這個時段。所以本次數據處理在最終時序分析中使用的數據時段為2 013.0~2 019.0 a。

2.2 基線解算及綜合平差

本文基線解算采用6 a的連續數據，考慮到半月周期，數據間隔采用7天一個時段，共計312個時段。采用GMAIT軟件對14個IGS站進行長基線解算。

對312個時段的基線解算進行精度統計。單天解的標準化均方根（normalized root mean square,NRMS）代表了基線解算的精度情況，通常情況下，NRMS小于0.3。除了11個時段的NRMS值超限之外，其他301個時段的NRMS值統計表3所示。

表3 單天解NRMS統計表

2.3 平差獲取坐標時間序列

本文將采用 GLOBK軟件對所有符合要求的單時段基線解進行整網約束平差。平差選擇 6個IGS站作為固定站，這6個分別為YSSK、KHAJ、PIMO、KIT3、SUWM、PIMO，對固定站的平差約束量采用其坐標和速度的一倍中誤差，平差基準采用國際地球參考框架（The International Terrestrial Reference Frame, ITRF）2014，參考歷元為2010.0。

通過格爾雷德（GLRED）解算整體綜合解，獲得站點在各個時段的坐標值，并形成坐標分量的長期坐標時間序列（N、E、U）。此處僅呈現出我國境內的LHAZ（拉薩站）站的單天剔除粗差之后的時間序列，如圖2所示。

從圖2可以看出，在LHAZ站上，、及三個方向上均有顯著的周期波動，尤其在方向上波動較大：方向波動振幅約6 mm、方向振幅約7 mm、方向振幅約20 mm。

圖2 LHAZ站網解時間序列（剔除粗差之后）

2.4 主成分分析（PCA）

采用主成分分析（PCA）模型時：1）去除時間序列中的階躍項、并對序列的趨勢項進行改正，改正后的序列只剩下年周期項和半年周期項；2）對序列殘差進行去除周期項分析，序列周期信息包括了由于地殼質量負荷引起的周期項，地殼質量負荷改正主要包括了如：大氣改正、積雪改正、土壤濕度改正、非潮汐改正，本文采用全球大氣資料Atmospheric-M、全球積雪資料Snow-M、全球土壤濕度資料Soil-M、非潮汐改正資料Ecco-M等模型，根據坐標序列中的坐標值進行差值計算，獲取所有站點的坐標時間序列值對應的質量負荷改正量，并實現對時間序列的改正，改正后獲得的所有站點坐標序列將作為一個整體序列群輸到主成分分析模型中去。

對改正后的相對純凈的時間序列數據進行PCA主成分分析。從PCA分析來看，CME的存在是造成殘差序列中存在明顯的周期變化的主要原因，尤其在豎直方向，周期性表現尤為明顯。為了對 CME進行分離，在濾波過程中，首先由 PCA法分別計算得到14個IGS站集群、在同一時間域的前12個主成分對3個坐標分量（、、）的空間響應，其中：主成分 1對坐標分量（、、）的貢獻率分別為 70.83%、69.32%、68.23%，主成分2對坐標分量（、、）的貢獻率分別為12.89%、14.43%、10.03%，主成分 3對坐標分量（、、）的貢獻率為 5.43%、9.56%、3.32%。

14個IGS站集群在同一時間域的前3個主成分，累計貢獻率累計為 89.166%、85.453%、83.95%，所有第一、第二、第三主成分的累計貢獻率綜合了大部分的空間響應（>87%）。

以上統計表明，采用PCA方法可得到前3個主成分作為該數據群的CME，然后在各個站點的坐標序列中，需要針對CME進行濾波處理。這樣濾波之后得到純凈的GNSS站點坐標時間序列。最后，基于純凈的時間序列計算精確的站點坐標和速度。

2.5 PCA濾波結果分析

通過地表負荷改正和 PCA濾波處理剔除了CME之后，各個站點時間序列曲線的周年運動的振幅有所減小：平均年振幅變化減少最小為3.22%；最大減少了65.23%；平均減少了32.52%。

圖3為LHAZ站的結果。從圖3可以看出，曲線振幅減小了，收斂程度更高了。

圖3 LHAZ站PCA濾波前后時序

在 ITRF2014框架下計算站點的框架坐標和速度場，并將結果與IGS網站公布的ITRF2014框架的坐標和速度進行比較，差異情況如圖4、圖5所示。

圖4 濾波后坐標值與IGS成果對比

圖5 濾波后速度值與IGS結果對比

從圖4、圖5可以看出，14個 IGS站點經過濾波處理后的解算的坐標與 IGS成果差異小于9 mm，其速度值與IGS成果差異小于3 mm/a。

3 結束語

本文利用 GAMIT/GLOBK軟件實現了區域GNSS網絡的綜合解算和平差，并采用 PCA方法對長期坐標時間序列進行了穩定性分析。重點介紹了如何采用濾波的方法來探測和剔除共模誤差，從而獲得純凈的坐標序列，提高了坐標精度。得到結論如下：

1）在全球GNSS網絡平差中，需對約束站給予合適的坐標和速度約束才能保證平差結果的精確。

2）共模誤差是GNSS坐標時間序列程序周期性信號的主要因素，采用主成分濾波可以有效去除共模誤差的影響。

3）實施負荷改正和主成分分析等區域濾波后，基站坐標時間序列趨于平穩，波動振幅有所減小，大大提高了時間序列的信噪比。

4）時間序列分析目的就是分類各種信息，并加以改正、剔除、過濾，從而獲取純凈的時間序列。基于純凈的時間序列數據平差得到更高精度的站點位坐標和速度場。