顧及共模誤差影響的澳大利亞GPS站長時間序列分析

凌聰聰，廖超明，2，康傳利，周聰林，楊翼飛，張振宇

(1.桂林理工大學 測繪地理信息學院，廣西 桂林 541006；2.南寧師范大學 自然資源與測繪學院，南寧 530001；3.廣西壯族自治區(qū)自然資源信息中心，南寧 530028)

GPS基準站是提供全球導航定位、 衛(wèi)星精密定軌、 動態(tài)參考框架維護的重要基礎設施, 同時也為大陸構造形變監(jiān)測、 地震預測、 地球動力學研究積累寶貴的觀測數(shù)據(jù)[1-4]。由于GPS站受儀器、 環(huán)境等因素影響, 導致觀測結果中含有多種誤差, 其中區(qū)域GPS站網中存在著一種與時空相關的誤差, 被稱為共模誤差(common mode errors, CME), 因此開展GPS時間序列共模誤差分析, 對獲取區(qū)域地殼現(xiàn)今構造運動特征具有重要意義。

對區(qū)域GPS站網共模誤差研究已取得許多的成果: 胡守超等[5]對南加州區(qū)域GPS站的共模誤差進行提取, 發(fā)現(xiàn)堆棧法(staking)、 主成分分析(PCA)和Karhunen-Loeve展開(KLE)均能有效提取GPS站網的共模誤差; 伍吉倉等[6]、 賀小星等[7]通過結合PCA和KLE方法增強空間濾波穩(wěn)健性, 有效剔除GPS站網的共模誤差, 進而提高GPS站坐標時間序列精度; 馬超等[8]使用staking、 PCA、 KLE法對南極半島地區(qū)的GPS站進行了空間濾波, 結果表明PCA濾波效果優(yōu)于其他兩種方法, 濾波后殘差時間序列振幅、 RMS都有所減小; 此外, 王健等[9]使用PCA估計并剔除歐洲大區(qū)域GPS站網的共模誤差, 發(fā)現(xiàn)空間濾波后部分站點最優(yōu)噪聲模型發(fā)生改變; 劉宗強等[10]對四川23個CORS站坐標時間序列進行研究, 通過堆棧法剔除CORS站的共模誤差, 發(fā)現(xiàn)空間濾波后的部分站點出現(xiàn)隨機游走噪聲。由此可見, 對于區(qū)域GPS站網共模誤差的提取方法已比較成熟, 為基于GPS觀測數(shù)據(jù)開展區(qū)域地殼現(xiàn)今構造運動監(jiān)測研究提供了技術支撐。

目前, 姜衛(wèi)平等[11]、 汪浩等[12]對澳大利亞區(qū)域進行了相關研究, 前者研究時間跨度對GPS站水平方向噪聲模型的影響, 后者利用GPS和GRACE數(shù)據(jù)研究區(qū)域地殼垂向季節(jié)變化。 本文選擇澳大利亞區(qū)域19個GPS站長跨度坐標時間序列作為研究對象, 采用PCA方法提取共模誤差, 利用譜指數(shù)估計和最大似然估計法分析剔除共模誤差前后的時間序列噪聲特性, 為開展區(qū)域地殼現(xiàn)今構造運動特征研究提供高精度坐標速度場估值。

1 分析方法

1.1 主成分分析

主成分分析是一種廣義的空間濾波方法, 該方法把一組相關數(shù)據(jù)變量通過線性變換成一組不相關的變量[6,13]。假設一個區(qū)域GPS站網有n個測站, 每個站點觀測m天, 且滿足n≤m關系。設X矩陣的每一列表示站點某一方向時間序列的殘差值, 每一行則表示某一歷元所有站點的殘差值。X矩陣的協(xié)方差陣為B, 其元素定義為

(1)

協(xié)方差矩陣B是一個n×n對稱矩陣, 可由正交分解為

B=VΛVT,

(2)

式中:Λ是協(xié)方差陣B的非零特征值構成的對角矩陣, 其中特征值按對角線降序排序;VT是由特征向量構成一個n×n正交矩陣。根據(jù)特征值可得到主成分為

(3)

式中:ak是X矩陣的第k個主成分;Vk為對應第k個主成分的特征向量。以降序排列的前幾個特征值代表殘差時間序列的最大貢獻, 而特征向量代表各站點的空間響應, 則PCA最終提取的CME可定義為

(4)

式中:ε(Mi,j)表示由p個主成分計算的第j個站點第Mi天的共模誤差。

1.2 譜指數(shù)法

許多地球物理信號能用某種常見的統(tǒng)計模型來表征, 可描述為冪律過程[14-16]。 GPS坐標時間序列噪聲信號也屬于地球物理信號, 也可描述為冪律過程, 將該隨機過程的時間域用功率譜函數(shù)描述為

Px(f)=P0(f/f0)k,

(5)

式中:f是時間頻率;P0和f0為歸一化常數(shù);Px(f)表示功率譜密度;k是頻譜指數(shù), 其為雙對數(shù)坐標系中功率譜密度P0(f)對頻率圖像的斜率, 取值在[-3,1]。 當-3

1.3 最大似然估計法

借助CATS時間序列分析軟件[17]采用最大似然估計法計算最大似然值(MLE), 同時估計噪聲分量值。對于給定的一組觀測值x, 使其協(xié)方差的聯(lián)合概率值l達到最大

(6)

對式(6)兩邊取對數(shù)

MLE=ln[l(x,C)]

(7)

其絕對值是反映噪聲模型的一個關鍵指標, 可通過其估計出不同噪聲模型的最大似然值, 選取其中最大MLE的模型作為最佳模型。蒙特卡羅提出, 兩種噪聲模型MLE差值大于2.9可作為噪聲模型的區(qū)分閾值[10], 即假設WH+RWN與WH+FN的MLE差值大于閾值2.9, 則可判定MLE大的噪聲模型更優(yōu)。

2 數(shù)據(jù)來源與預處理

2.1 GPS數(shù)據(jù)

選取均勻分布于澳大利亞地區(qū)的19個GPS站(圖1)、 時間跨度為2010-06—2020-07的坐標時間序列數(shù)據(jù), 數(shù)據(jù)來源于內達華大地測量實驗室(http: //geodesy.unr.edu/gps_timeseries/)[18]。原始時間序列數(shù)據(jù)采用GIPSY/OASIS-II軟件進行處理, 同時已對電離層、 對流層延遲、 天線相位中心偏差進行校正, 并加入了極潮、 海潮及固體潮等模型改正, 最終處理得到IGS14框架下GPS站坐標時間序列。

圖1 澳大利亞區(qū)域19個GPS站點序號及分布Fig.1 Distribution of 19 GPS stations in Australia

2.2 GPS時間序列處理

為了合理估計GPS時間序列各項參數(shù)，需要對其進行預處理。一般GPS時間序列模型可表示為

y(ti)=a+bti+csin(2πti)+dcos(2πti)+esin(4πti)+

(8)

式中:ti表示以年為單位的時間; 待求系數(shù)a、b分別表示為測站的初始位置和速率; 系數(shù)c、d表示站點周年項運動; 系數(shù)e、f表示站點半周年項運動;gi為由各種原因(如儀器天線更換等)引起階躍式坐標突變;H為一階梯函數(shù);Tgj為發(fā)生坐標跳變時間;vi為觀測誤差項。

受到外界條件、 GPS接收機軟硬件等因素影響, 部分GPS站點某些年積日會出現(xiàn)觀測數(shù)據(jù)丟失或觀測數(shù)據(jù)質量不好等情況, 從而導致 GPS 時間序列存在缺失或粗差等問題。本文首先采用絕對中位數(shù)法對每個站點的時間序列進行粗差剔除; 然后用三次多項式進行插值與補齊缺失數(shù)據(jù), 得到等間隔的“干凈”時間序列; 最后對原始時間序列去均值項、 速度項以及階躍項, 得到各GPS站點的殘差時間序列。ALIC站點原始坐標時間序列經過預處理后得到的殘差時間序列如圖2所示。

圖2 ALIC站點殘差時間序列Fig.2 Residual time series of Station ALIC

3 GPS站坐標時間序列噪聲特性分析

3.1 共模誤差提取

對澳大利亞地區(qū)19個站點的殘差時間序列進行主成分分析, 分別得到N、 E、 U向的貢獻率, 結果如圖3所示?？芍? 采用PCA分析得到第一主成分貢獻率在N、 E、 U向為36.9%、 42.7%、 39.9%; 第二主成分貢獻率在N、 E、 U向為22.1%、 16.4%、 9.6%; 第三主成分貢獻率在N、 E、 U向為7.9%、 8.4%、 8.8%; 前3個主成分累計貢獻率在N、 E、 U向66.9%、 67.5%、 58.3%, 即前3個主成分綜合了大部分的空間響應, 因此后續(xù)僅對前3個主成分進行討論。

圖3 N、 E、 U方向的主成分貢獻率占比Fig.3 Contribution of principal components in N，E and U direction

研究表明, 當某一主分量模式中50%以上測站的標準化空間響應大于25%, 且該模式的特征值貢獻率大于1%, 則可作為測站間的共有模式[19]。為確定式(4)中的主分量個數(shù)p來計算CME, 將特征向量除以其絕對值最大的元素, 得到標準化的空間特征向量, 如圖4所示。第一主分量每個站點N、 E、 U方向上的標準化空間響應均大于25%, 主成分貢獻率超過1%, 而第二、 三主分量在3個方向上并不滿足該條件, 因此本文用第一主成分來計算區(qū)域各站的CME。以ALIC站為例, 扣除了CME前后的殘差時間序列和殘差RMS分別如圖5、 6所示。

圖4 各GPS站點標準化空間響應Fig.4 Standardized spatial responses of GPS stations

圖5 ALIC站濾波前后的殘差時間序列Fig.5 Residual time series before and after filtering at Station ALIC

圖6 各站點濾波前后的殘差RMS值Fig.6 Residual RMS of each station before and after filtering

從圖5可看出, 濾波前殘差時間序列N、E向存在較弱的周期波動,U向周期波動明顯,濾波后殘差時間序列在N、 E、 U向的振幅都有所減小。濾波后各站N、 E、 U向的殘差RMS都有減小, 平均減小約29.7%、 25.5%、 24.2%, 說明濾波后各站點的坐標分量估值精度得到有效提高(圖6)。此外, 部分站點濾波前后的殘差RMS值相差較小, 而部分站點則相差較大, 這可能與站點空間分布有一定關系。如U向站點序號5、 8的殘差RMS值相差較小, 對應站點DARW、 HOB2分布于澳大利亞邊緣。

3.2 最優(yōu)噪聲模型分析

用CATS軟件獲取空間濾波前后的時間序列在N、 E、 U向上的譜指數(shù), 如表1所示。各站點濾波前后的時間序列的譜指數(shù)在-2～0, 說明該區(qū)測站濾波前后的噪聲模型并非單一的白噪聲模型。

表1 空間濾波前后測站的譜指數(shù)

為確定最優(yōu)噪聲模型, 用CATS軟件獲取濾波后各測站在WN、 WN+FN、 WN+RWN、 WN+FN+RWN噪聲模型下的最大似然值, 并將其他3種模型與WN模型的最大似然值作差, 結果如圖7所示。WN+FN和WN+FN+RWN模型MLE差值基本相等, 且兩個模型都比WN+RWN模型的MLE差值大。根據(jù)蒙特卡羅準則判定WN+FN和WN+FN+RWN模型更優(yōu), 但兩個模型不具有可分性。現(xiàn)假設WN+FN+RWN為區(qū)域最優(yōu)噪聲, 且在該噪聲模型下估計濾波前后各站點的噪聲分量。

圖7 19個GPS測站N、 E、 U向MLE差值Fig.7 MLE difference of 19 stations in N，E，U direction

限于篇幅, 表2僅列出濾波前后的9個GPS站基于WN+FN+RWN模型的噪聲分量?？芍? 1)澳大利亞地區(qū)GPS站N、 E、 U向最優(yōu)噪聲模型為WN+FN+RWN, 9個站的坐標分量出現(xiàn)有白噪聲、 閃爍噪聲及隨機游走噪聲, 因此該區(qū)域GPS測站在N、 E、 U向不僅需要考慮WN+FN, 還應加入RWN噪聲; 2)空間濾波后測站N、 E、 U向白噪聲和閃爍噪聲的估值都有減小, 但U向有部分站點(ALIC、 TBOB)的白噪聲有增加情況, 該原因有待進一步研究; 3)部分站點(ALIC、 KARR等)空間濾波后還原出了被掩蓋的隨機游走噪聲, 這表明空間濾波能夠將被高頻信號掩蓋的隨機游走噪聲顯示出來。

表2 基于WN+FN+RWN模型的噪聲分量估計值

4 區(qū)域坐標速度場分析

4.1 基于WN+FN+RWN模型的坐標速度場估計

根據(jù)GPS站時間序列分析結果, 采用WN、 WN+FN+RWN模型分別估計GPS站N、 E、 U向的速度及其精度, 如表3所示?？芍? 模型WN+FN+RWN和WN在N、 E、 U向速率平均偏差為0.12、 0.06、 0.19 mm/a, 最大差值分別為0.42、 0.18、 0.93 mm/a。在速度不確定性方面, 最優(yōu)噪聲模型估計下N、 E、 U向速度不確定性為WN模型下的十至幾十倍。由此說明, 僅考慮白噪聲的影響并不能真實地反映測站速度及其不確定性, 尤其是對于速度不確定性的估計。

從內達華大地測量實驗室收集到的站點速率與WN+FN+RWN模型估計站點速率作差, 得到GPS站的速率差值(圖8)。由19個GPS測站點N、 E、 U向速率差異可知, GPS站點N、 E、 U向速率差異范圍-0.99～0.19 mm/a、 -0.39～0.73 mm/a、 -1.74～1.25 mm/a, 差異絕對值的均值分別為0.28、 0.17、 0.45 mm/a?；赪N+FN+RWN模型估計的GPS站坐標速率與內達華大地測量實驗室處理的結果相比較, 從整體上具有較好的一致性, 其差異主要表現(xiàn)在HIL1、 RHPT、 TOW2站點的垂向速率估值; 同時, 這3個站點速率估值的不確定性相對較大, 分別達到了±1.24、 ±1.09、 ±0.84 mm/a(表3)。

表3 WN+FN+RWN模型與WN模型下的GPS站速度及其精度

圖8 GPS站點N、 E、 U方向速率差Fig.8 Rate difference of GPS stations in N，E and U direction

4.2 區(qū)域構造運動分析

以WN+FN+RWN組合模型下N、 E、 U向的速率繪制澳大利亞板塊坐標速度場, 同時繪制輔助線上下10°范圍的GPS站速率子圖。由圖9a和表3分析可知, 測站N向速率在54.3～60.5 mm/a, E向速率在13.9～39.0 mm/a, 測站水平方向速率的差異較大, 說明澳大利亞板塊內部存在一定形變。IGS14框架下GPS站整體水平運動方向為北東方向, 平均速率約為65.1 mm/a; 在太平洋海板塊西南端的北西向推擠力、 東南印度洋中脊的北東向分離推擠力作用下, 測站運動方向由北北東逐漸向北方向過渡; 自澳大利亞地區(qū)的北西向南東方向, 水平速率估值逐漸減小。

圖9 WN+FN+RWN模型下的澳大利亞水平(a)與垂直(b)速度場Fig.9 Horizontal(a) and vertical(b) velocity fields based on WN+FN+RWN model in Australia

除TBOB、 CEDU站點近零速率的正值外(圖9b、 表3), 其余站點速率為-2.15～0 mm/a, 各GPS站點垂向速率差異較小。GPS站點監(jiān)測到澳大利亞板塊現(xiàn)今地殼垂向運動整體處于緩慢下沉狀態(tài), 平均速率約為-0.95 mm/a; 澳大利亞地區(qū)的東、 南面下沉速度相對于西、 北面更快。

5 結 論

通過分析濾波前后和不同噪聲模型下澳大利亞地區(qū)19個GPS站2010-06—2020-07的長時間序列, 得到以下結論:

(1)對區(qū)域GPS站開展空間濾波能夠有效降低白噪聲和閃爍噪聲。此外, 空間濾波能夠將被掩蓋的隨機游走噪聲還原出來, 因此應考慮對該區(qū)域進行空間濾波。

(2)澳大利亞地區(qū)GPS站坐標時間序列應考慮采用WN+FN+RWN的組合。雖然白噪聲模型與WN+FN+RWN模型的速率估值非常接近, 但只考慮白噪聲影響下的速度不確定性將被低估。因此, 在該區(qū)域開展GPS站時間序列分析時, 應考慮采用多噪聲組合模型代替白噪聲模型。

(3)在IGS14框架下, 澳大利亞區(qū)域19個GPS站現(xiàn)今地殼水平運動方向整體表現(xiàn)為北北東向, 平均速率約為65.1 mm/a, 且速率自北西向南東逐漸減小; 垂直運動處于緩慢下沉狀態(tài), 平均速率約為-0.95 mm/a, 且東、 南面下沉速率相對于西、 北面更快。