基于獨立分量法的新疆GNSS時間序列共模誤差分析

雷傳金，魏冠軍，高茂寧，張沛

( 1. 蘭州交通大學 測繪與地理信息學院, 蘭州 730070；2. 地理國情監(jiān)測技術應用國家地方聯(lián)合工程研究中心, 蘭州 730070；3. 甘肅省地理國情監(jiān)測工程實驗室, 蘭州 730070 )

0 引 言

全球衛(wèi)星導航系統(tǒng)(GNSS)能提供高精度、可靠的坐標時間序列，現(xiàn)已廣泛應用于地球物理研究中，如板塊運動變化、參考框架構建與維持、區(qū)域地殼形變等[1-3]. 然而，GNSS觀測數(shù)據(jù)中存在一種與空間相關的噪聲，即共模誤差(CME)[4]. CME的存在會影響測站的參數(shù)估計，增加了測站參數(shù)信息的不確定性，這將導致解譯出不準確的地球物理信息，故需采用濾波方法來削弱其干擾，以便解譯出真實的地球物理信息. 針對這一問題，國內外學者們引入了堆棧法[4]、奇異譜分析[5]、主成分分析(PCA)[6-7]等濾波方法成功分離出CME. 在這些方法中，PCA濾波應用最為廣泛，該方法是基于二階統(tǒng)計量將輸入信號分解成方差最大分量來提取CME[6]，但未充分利用CME的高階統(tǒng)計信號，使得某一成分包含不同的地球物理信息，在一定程度上影響了CME分離的準確性[8].

相較于PCA，獨立分量分析(ICA)是一種基于高階統(tǒng)計量對盲源分離的信號進行分解的方法，對分離出獨立的非高斯信號方面展現(xiàn)出優(yōu)越的性能，能有效地分離出獨立的信號[9]，現(xiàn)已廣泛用于大地測量數(shù)據(jù)處理，如GNSS時間序列分析[8-9]、全球時變重力信號提取[10]、InSAR時序分析[11]等領域. 在GNSS網(wǎng)空間濾波應用中，文獻[8]使用ICA方法對中國區(qū)域的259個GPS站進行濾波分析，取得較好地濾波效果，但ICA分離后信號顯著性問題未得到較好地解釋，未討論濾波后對測站速度場、周期項的影響.

新疆地區(qū)是盆地、山地地震構造區(qū)，是研究地殼運動的良好試驗場，但一些學者在提取區(qū)域地球物理信號時，僅考慮時域上的噪聲項影響，未考慮空域上的CME. 因此，本文引入ICA對新疆地區(qū)GNSS坐標時間序列CME進行提取，并與PCA濾波效果進行對比分析，驗證了ICA在提取區(qū)域CME時的可行性，為CME物理成分的科學解釋提供參考. 此外，高精度地提取線性速率等地球物理信號，為更好地解釋這一區(qū)域的構造運動提供參考價值.

1 空間濾波

1.1 PCA空間濾波

設某一區(qū)域有m個時間跨度為n天的觀測站，且n＞m，數(shù)據(jù)矩陣X(n×m) 為測站的殘差序列，則可計算X的協(xié)方差陣CX

對CX進行正交分解

式中：Λ為具有r個非零特征值的對角矩陣；V為m×m維特征向量矩陣. 采用正交基V來擴充X矩陣，可得

即

式中：ak(t) 為X的第j個主分量；vk(s) 為對應主分量ak(t) 的空間響應(SR). 將得到的特征值 λi降序排列，以便提取X中的CME信號. 假設前p個方差較大的主分量代表CME，則CME可表示為[6]

1.2 ICA空間濾波

設有m個時間序列為X的觀測站，每個觀測序列Xi(t)=[x1(t),x2(t),···,xn(t)]T可由n個相互獨立的未知信號y進行疊加，某一歷元t的線性表達式為

式中：xt為t時刻m×1 維 觀測向量；為t時刻第j個信號源；aij為第j個信 號 源對i個觀 測值的空間響 應 情況，且滿足1 ≤i≤m,1≤j≤n. 將所有觀測數(shù)據(jù)寫成矩陣形式

引入混合矩陣W來恢復A的 原始信號，即W=A-1. 為了解算W，本文利用基于負熵峭度的魯棒方法計算非高斯性. 其近似負熵為

式中：Y=WTZ代表獨立分量，Z為觀測向量中心化、白化處理后的值；g函數(shù)的形式為g(y)=yexp(-y2/2)，通過式(8)便可得到W的最優(yōu)解及每個獨立分量. 本文采用文獻[9]中時域上相互獨立的tICA進行解算混合矩陣W，并計算出第k個獨立信號的貢獻率為

式中： ‖ ·‖ 為F-范數(shù)；rk為第k成分的貢獻率. 將所有獨立分量按貢獻比降序排列，用前幾個顯著的獨立分量表示最顯著信號，則CME可表示為[8]

2 GNSS數(shù)據(jù)處理

本文對我國新疆地區(qū)連續(xù)站的原始GNSS時間序列進行CME分析，分析站點的時間跨度為2011—2018年，數(shù)據(jù)缺失率較低，平均缺失率為5%，能夠有效地避免數(shù)據(jù)缺失對主頻項分析的影響. 數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產品服務平臺，GNSS站點分布如圖1所示.

圖 1 新疆地區(qū)GNSS站點分布

2.1 數(shù)據(jù)預處理

原始時間序列通過粗差剔除、階躍項改正、數(shù)據(jù)插值等預處理工作，獲得“干凈”的GNSS殘差時序.本文采用三倍四分位法進行粗差剔除，選用克里金卡爾曼濾波算法(KKF)對缺失數(shù)據(jù)進行插值[12]，利用最小二乘擬合去除趨勢項，得到殘差時間序列[13]. 圖2給出了XJQH站預處理后的殘差時間序列，黑點為插值效果，藍點為原始殘差序列. 由圖2可知，經預處理后，測站坐標時序的粗差、周期項得以有效消除，較真實地還原出缺失數(shù)據(jù).

圖 2 XJQH站預處理后的殘差序列

2.2 共模誤差的提取

2.2.1 PCA濾波

對經預處理后的GNSS殘差時間序列進行PCA濾波，計算出N、E、U三個方向分量的PC1的貢獻率分別為59.91%、64.34%、69.95%，前三個分量的累計貢獻分別為79.11%、77.92%、83.52%. 各站點的N、E、U三個方向分量均在PC1所對應的空間響應上表現(xiàn)出一致性，故本文采用PC1提取區(qū)域網(wǎng)的CME.

根據(jù)式(5)對各測站進行CME提取并剔除，圖3(a)給出了經PCA濾波剔除CME前后XJSS站殘差時序變化，綠色為原始GNSS殘差時間序列，藍色為經PCA濾波后的殘差序列. 由圖3(a)可知，經PCA剔除CME后，測站殘差時序的離散程度明顯降低；XJSS站的E、N、U三個方向分量的殘差時序均方根(RMS)值分別降低了58.04%、58.87%、59.56%，PCA濾波后的N、E、U三個方向的殘差序列的平均RMS分別減少了42.44%、46.14%、48.94%，顯著提高了測站時序信號的信噪比.

圖3(b)為經ICA濾波剔除CME前后XJSS站殘差時序變化，綠色為原始GNSS殘差時間序列，藍色為經過ICA濾波后的殘差序列，與PCA濾波后得到的結果大致相同，測站殘差時序的離散程度也明顯降低，XJSS站的E、N、U三個方向分量的殘差時序RMS值分別降低了45.55%、43.28、47.42%，ICA濾波后殘差序列的平均RMS分別減少31.83%、32.29%、35.49%.

圖 3 XJSS站殘差時序

2.2.2 ICA濾波

準確選取獨立分量是進行ICA濾波的關鍵，根據(jù)式(9)計算出獨立分量的貢獻比并降序排序，再結合對應的SR來選取恰當?shù)莫毩⒎至? 從實驗計算結果可分析出，N方向的IC1～IC4、IC6～IC7對應的測站具有一致的SR，E方向的IC1～IC3、IC5～IC7對應的測站具有一致的SR，U方向的IC1～IC6對應的測站具有一致的SR，而其他獨立分量的SR無明顯的空間一致性. 限于篇幅，本文現(xiàn)給出N、E、U前3個標準化后的獨立分量及其對應的歸一化SR，如圖4所示，紅色箭頭向上表示空間正響應，藍色箭頭向下表示空間負響應.

由圖4可知，IC1在天山區(qū)域空間響應較強，在柴達木盆地周圍空間響應較弱. 相反，IC3在柴達木盆地周圍空間響應較強，在天山區(qū)域空間響應較弱.這可能與地表質量荷載有關. 而各測站在IC2上呈負響應，且響應程度較強，其他空間響應呈現(xiàn)出較均勻的強度，具體原因后續(xù)將進一步討論. 根據(jù)式(10)計算出各分量的CME，再將具有一致SR的ICA分量進行累加計算得到各測站的CME.

圖 4 N、E、U方向分量前3個獨立分量及其歸一化SR

表1統(tǒng)計了經PCA、ICA濾波后的測站殘差時 序RMS值，其經ICA濾波后的測站殘差時序的RMS值不如PCA濾波后的結果，這是因為ICA能夠最大化非高斯分量，把異常信號當成獨立的分量分離出來，使得異常信號能夠和CME分離. 而PCA是利用二階統(tǒng)計量將輸入信號分解成方差最大分量來提取CME，存在異常信號被作為CME成分剔除，一定程度上影響CME提取的準確度.

表 1 PCA、ICA濾波前后殘差序列的平均RMS變化 mm

3 討 論

3.1 CME對GNSS測站時間序列的影響

為了探討CME對GNSS測站時間序列的影響，經探測識別剔除異常信號，選用白噪聲+冪律噪聲(WN+PL)組合模型[14-15]進行時間序列分析. 圖5為經ICA濾波前后測站速度場不確定度變化. 經ICA濾波剔除CME后，N、E、U三個方向分量的速度不確定度分別減少44.14%、38.49%、43.32%，說明ICA能夠降低測站時間序列的不確定性，從而提高GNSS坐標序列的可靠性.

圖6～7為經ICA、PCA濾波后得到的新疆區(qū)域的水平和垂直速度場. 從圖6～7可以看出，在水平方向上，新疆地區(qū)整體運動趨勢以E方向為主，天山東部及準噶爾盆地的運動趨勢為E方向，新疆塔里木盆地為NE方向的運動趨勢，最大速度為WUSH站，達到33.7 mm·a-1，其總體運動速度由西向東逐漸減弱. 在垂直方向上，大多數(shù)測站運動速度小于2 mm·a-1，最大抬升為XJKL站，達到5.8 mm·a-1，經查閱資料，該地區(qū)主要是由于油田注水導致了該站出現(xiàn)異常抬升[16]. 天山地區(qū)整體上呈現(xiàn)出抬升的趨勢，塔里木盆地北部、中部及準噶爾盆地呈現(xiàn)下降趨勢. 但這兩種方法在天山地區(qū)的XJYN、XJXY測站、塔里木盆地南緣的XJHT、XJYT及盆地中部的XJQM測站呈現(xiàn)相反的運動趨勢. 通過查閱相關文獻，天山地區(qū)受盆地擠壓而產生抬升的運動趨勢，與ICA濾波后的結果一致.

圖 5 ICA濾波前后測站速度場不確定度變化

圖 6 經ICA濾波后新疆區(qū)域水平和垂直速度場

圖 7 經PCA濾波后新疆區(qū)域水平和垂直速度場

3.2 CME對GNSS測站周期項的影響

GNSS坐標時序中存在明顯的周期特性[17]，通過計算得到剔除CME前后測站時間序列的基本參數(shù)，如圖8所示. 由圖8可知，除個別測站外，經ICA濾波后其測站N、E、U三個方向分量的年振幅基本約在1 mm、0.5 mm、3 mm，且測站垂直方向的年周期項振幅高于測站水平方向. 對ICA濾波前后測站的振幅進行對比，可發(fā)現(xiàn)絕大多數(shù)測站經濾波后其振幅較濾波前更具有一致性，表明CME會影響季節(jié)項振幅的估算，得到不準確的物理參數(shù)信息. 表2統(tǒng)計了濾波前后測站的年振幅、半年振幅的標準差變化，可以發(fā)現(xiàn)剔除CME后測站坐標時間序列參數(shù)估計的不確定度有了顯著的下降，且N、E、U三個方向坐標分量的年振幅標準差分別減少44.49%、32.89%、43.12%，半年振幅標準差分別減少45.33%、33.07%、43.18%，表明CME對測站坐標時序振幅有較大的影響.

圖 8 ICA濾波前后年振幅變化統(tǒng)計圖

表 2 測站振幅標準差統(tǒng)計

3.3 CME的周期分析

CME是區(qū)域內大部分測站都存在的多種時空相關的誤差集. 為了對CME內的各種誤差信號進行分離，本文使用集合經驗模態(tài)分解(EEMD)對CME序列進行分解，再采用快速傅里葉變換(FFT)對分離的序列進行周期探測[18]. 通過EEMD計算出IMF1～IMF3為噪聲序列， I MF4～IMF7為有效成分，其中IMF6分量具有明顯季節(jié)性周期變化規(guī)律， I MF7分量呈現(xiàn)出明顯年周期變化規(guī)律，說明CME是一種與時間相關的誤差. 圖9給出U方向 I MF4～IMF7分量分析結果. 為了直觀準確分析周期變化，對EEMD分解后的IMF分量進行了FFT周期探測，如圖10所示.功率曲線表現(xiàn)出先增后減的變化趨勢，最大值分別在56.3天、83.6天、276天和345天. 說明CME具有半季節(jié)性、季節(jié)性、半年及年的周期變化規(guī)律，這些周期可能是由于參考框架、地表質量荷載、鐘差以及未建模的衛(wèi)星軌道等不確定性因素導致的[19-21]，具體成因還需要進一步的研究和探討.

圖 9 U方向CME序列的EEMD分解

圖 10 經EEMD的CME序列的FFT周期探測

針對PCA提取CME時，存在異常信號被作為CME成分的問題，采用ICA方法對新疆地區(qū)的GNSS坐標時間序列進行空間濾波，并與PCA的濾波結果進行驗證. 發(fā)現(xiàn)RMS值均有大幅下降，說明兩種方法均能有效提取CME，而ICA能從觀測數(shù)據(jù)中分離出混合的潛在成分和來源，能提取出更準確的CME信號.

4 結束語

通過采用ICA方法進行空間濾波，測站時間序列的速度不確定度、年周期、半年周期振幅標準差均有明顯下降，濾波后絕大多數(shù)測站的振幅較濾波前更具有一致性，說明ICA濾波能夠有助于提高測站時間序列的參數(shù)估計的精度.

經ICA濾波后，新疆地區(qū)在水平方向整體以E方向的運動趨勢為主，總體的運動速度由西向東逐漸減弱. 在垂直方向上，天山地區(qū)為抬升趨勢，塔里木盆地北部、中部及準噶爾盆地呈現(xiàn)下降趨勢. ICA濾波后同時在對CME進行周期性分析時發(fā)現(xiàn)，CME具有半季節(jié)性、季節(jié)性、半年周期及年周期變化規(guī)律，產生這些周期的具體成因還需要進一步的研究和探討.

致謝：感謝中國地殼運動觀測網(wǎng)絡(CMONOC)提供的GNSS坐標時間序列數(shù)據(jù).