區(qū)域GNSS基準(zhǔn)網(wǎng)解算及其時間序列分析

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

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

0 引言

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

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

1 基準(zhǔn)網(wǎng)解算策略及時間序列分析方法

1.1 區(qū)域GNSS基準(zhǔn)網(wǎng)解算策略

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

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

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

1.2 GNSS時間序列分析方法

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

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

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

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

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

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

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

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

2 區(qū)域GNSS網(wǎng)綜合平差及其時間序列分析

2.1 數(shù)據(jù)收集及整理

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

表1 IGS站點(diǎn)名稱及時間跨度表

圖1 14個IGS站點(diǎn)位分布圖

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

表2 改正模型名稱及時間跨度度統(tǒng)計(jì)表

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

2.2 基線解算及綜合平差

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

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

表3 單天解NRMS統(tǒng)計(jì)表

2.3 平差獲取坐標(biāo)時間序列

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

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

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

圖2 LHAZ站網(wǎng)解時間序列（剔除粗差之后）

2.4 主成分分析（PCA）

采用主成分分析（PCA）模型時：1）去除時間序列中的階躍項(xiàng)、并對序列的趨勢項(xiàng)進(jìn)行改正，改正后的序列只剩下年周期項(xiàng)和半年周期項(xiàng)；2）對序列殘差進(jìn)行去除周期項(xiàng)分析，序列周期信息包括了由于地殼質(zhì)量負(fù)荷引起的周期項(xiàng)，地殼質(zhì)量負(fù)荷改正主要包括了如：大氣改正、積雪改正、土壤濕度改正、非潮汐改正，本文采用全球大氣資料Atmospheric-M、全球積雪資料Snow-M、全球土壤濕度資料Soil-M、非潮汐改正資料Ecco-M等模型，根據(jù)坐標(biāo)序列中的坐標(biāo)值進(jìn)行差值計(jì)算，獲取所有站點(diǎn)的坐標(biāo)時間序列值對應(yīng)的質(zhì)量負(fù)荷改正量，并實(shí)現(xiàn)對時間序列的改正，改正后獲得的所有站點(diǎn)坐標(biāo)序列將作為一個整體序列群輸?shù)街鞒煞址治瞿Ｐ椭腥ァ?/p>

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

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

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

2.5 PCA濾波結(jié)果分析

通過地表負(fù)荷改正和 PCA濾波處理剔除了CME之后，各個站點(diǎn)時間序列曲線的周年運(yùn)動的振幅有所減?。浩骄暾穹兓瘻p少最小為3.22%；最大減少了65.23%；平均減少了32.52%。

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

圖3 LHAZ站PCA濾波前后時序

在 ITRF2014框架下計(jì)算站點(diǎn)的框架坐標(biāo)和速度場，并將結(jié)果與IGS網(wǎng)站公布的ITRF2014框架的坐標(biāo)和速度進(jìn)行比較，差異情況如圖4、圖5所示。

圖4 濾波后坐標(biāo)值與IGS成果對比

圖5 濾波后速度值與IGS結(jié)果對比

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

3 結(jié)束語

本文利用 GAMIT/GLOBK軟件實(shí)現(xiàn)了區(qū)域GNSS網(wǎng)絡(luò)的綜合解算和平差，并采用 PCA方法對長期坐標(biāo)時間序列進(jìn)行了穩(wěn)定性分析。重點(diǎn)介紹了如何采用濾波的方法來探測和剔除共模誤差，從而獲得純凈的坐標(biāo)序列，提高了坐標(biāo)精度。得到結(jié)論如下：

1）在全球GNSS網(wǎng)絡(luò)平差中，需對約束站給予合適的坐標(biāo)和速度約束才能保證平差結(jié)果的精確。

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

3）實(shí)施負(fù)荷改正和主成分分析等區(qū)域?yàn)V波后，基站坐標(biāo)時間序列趨于平穩(wěn)，波動振幅有所減小，大大提高了時間序列的信噪比。

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