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

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

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

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

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

1 分析方法

1.1 主成分分析

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

(1)

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

B=VΛVT,

(2)

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

(3)

式中:ak是X矩陣的第k個(gè)主成分;Vk為對(duì)應(yīng)第k個(gè)主成分的特征向量。以降序排列的前幾個(gè)特征值代表殘差時(shí)間序列的最大貢獻(xiàn), 而特征向量代表各站點(diǎn)的空間響應(yīng), 則PCA最終提取的CME可定義為

(4)

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

1.2 譜指數(shù)法

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

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

(5)

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

1.3 最大似然估計(jì)法

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

(6)

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

MLE=ln[l(x,C)]

(7)

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

2 數(shù)據(jù)來(lái)源與預(yù)處理

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

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

圖1 澳大利亞區(qū)域19個(gè)GPS站點(diǎn)序號(hào)及分布Fig.1 Distribution of 19 GPS stations in Australia

2.2 GPS時(shí)間序列處理

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

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

(8)

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

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

圖2 ALIC站點(diǎn)殘差時(shí)間序列Fig.2 Residual time series of Station ALIC

3 GPS站坐標(biāo)時(shí)間序列噪聲特性分析

3.1 共模誤差提取

對(duì)澳大利亞地區(qū)19個(gè)站點(diǎn)的殘差時(shí)間序列進(jìn)行主成分分析, 分別得到N、 E、 U向的貢獻(xiàn)率, 結(jié)果如圖3所示?？芍? 采用PCA分析得到第一主成分貢獻(xiàn)率在N、 E、 U向?yàn)?6.9%、 42.7%、 39.9%; 第二主成分貢獻(xiàn)率在N、 E、 U向?yàn)?2.1%、 16.4%、 9.6%; 第三主成分貢獻(xiàn)率在N、 E、 U向?yàn)?.9%、 8.4%、 8.8%; 前3個(gè)主成分累計(jì)貢獻(xiàn)率在N、 E、 U向66.9%、 67.5%、 58.3%, 即前3個(gè)主成分綜合了大部分的空間響應(yīng), 因此后續(xù)僅對(duì)前3個(gè)主成分進(jìn)行討論。

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

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

圖4 各GPS站點(diǎn)標(biāo)準(zhǔn)化空間響應(yīng)Fig.4 Standardized spatial responses of GPS stations

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

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

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

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

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

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

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

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

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

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

4 區(qū)域坐標(biāo)速度場(chǎng)分析

4.1 基于WN+FN+RWN模型的坐標(biāo)速度場(chǎng)估計(jì)

根據(jù)GPS站時(shí)間序列分析結(jié)果, 采用WN、 WN+FN+RWN模型分別估計(jì)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)噪聲模型估計(jì)下N、 E、 U向速度不確定性為WN模型下的十至幾十倍。由此說(shuō)明, 僅考慮白噪聲的影響并不能真實(shí)地反映測(cè)站速度及其不確定性, 尤其是對(duì)于速度不確定性的估計(jì)。

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

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

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

4.2 區(qū)域構(gòu)造運(yùn)動(dòng)分析

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

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

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

5 結(jié) 論

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

(1)對(duì)區(qū)域GPS站開展空間濾波能夠有效降低白噪聲和閃爍噪聲。此外, 空間濾波能夠?qū)⒈谎谏w的隨機(jī)游走噪聲還原出來(lái), 因此應(yīng)考慮對(duì)該區(qū)域進(jìn)行空間濾波。

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

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