利用主成分分析法分析GNSS坐標時間序列

劉曉祥 高二濤 羅 益 付波霖

1 重慶市永川區規劃和自然資源局，重慶市人民北路6號，402160 2 桂林理工大學測繪地理信息學院，桂林市雁山街319號，541006

陸態網自運行以來，為我國的地震趨勢報告、地殼形變規律研究、邊界勘查、測繪基準建設以及高精度的坐標框架維持提供了大量可靠的GNSS觀測數據[1-4]。但是受外界各種因素的干擾，GNSS站點存在數據缺失、點位突變、粗差以及異常站點等現象[5]，對觀測站高精度坐標的提取具有較大影響。

此外，想要系統地研究基準站坐標時間序列的變化規律，更深層次地了解驅使基準站運動的內部機制，需要對GNSS基準站連續觀測坐標時間序列的噪聲特性及運動規律進行分析。基于主成分分析(principal component analysis, PCA)的GNSS時間序列研究已經很成熟，并取得了較理想的研究成果[6]。Dong等[7]充分考慮GNSS觀測網的時空相關性因素，綜合PCA和Karhunen-Loeve分解方法來扣除共模誤差。袁林果等[8]利用主成分空間濾波方法有效提高了香港地區GPS坐標序列信號的信噪比。隨著GNSS基準站坐標時間序列觀測長度的不斷增加，以及軟硬件分析工具和方法的不斷優化，對GNSS基準站坐標時間序列的研究也越來越方便。

本文利用主成分分析方法對陸態網224個GNSS基準站2010～2016年的坐標時間序列進行分析。首先分析各站點原始坐標序列，并進行突變項擬合、粗差剔除、缺失時間序列插值等預處理；然后對預處理后的結果分N、E、U方向組建坐標時間序列矩陣并進行主成分分析；最后對各主分量及其對應的空間特征向量進行分布特征、共模誤差、站點響應區域等分析，并討論異常站點對3個方向PCA結果的影響。

1 陸態網坐標時間序列預處理

1.1 數據選取

GNSS基準站原始坐標時間序列中不僅包含由于觀測設備改變、更換天線、板塊位移、遠場地震等引起的趨勢變化，還包括地震引發的同震位移、點位突變以及數據缺失等造成的影響[5]。這些因素易導致主成分分析產生較大的負面影響，因此必須采取有效措施對原始坐標時間序列進行預處理。為了提升原始坐標序列的可靠性，本文將扣除缺失數據大于30%的基準站點，共剔除26個不達標的站點；選取剔除后的224個基準站時序坐標作為分析對象，截取這些站點2010～2016年的時序觀測結果作為實驗數據。

1.2 GNSS站點時間序列預處理

本文首先對實驗選取的GNSS坐標時間序列中的突變項進行消除，對突跳歷元出現之后的數據統一進行突變變化值去除，依據GNSS觀測文件聯合各個方向的時序散點圖來確定突跳時刻。假設站點N、E、U方向均有突變信號，在采用最小二乘法進行初始擬合時，盡量大范圍地篩選出突跳的歷元，然后比較擬合前后的時序散點圖，對可能的突跳序列在N、E、U方向進行判別。

圖1為BJFS站在去除突跳項前后N、E、U方向的對比，可以看到，通過擬合處理后，突跳位置得到有效補充。此外對于因儀器設備更換而引起的突變，可以通過各個站點的log文件獲取；對于其他未知原因造成的突變，為了確保該類突變的擬合效果，通過目視排查獲取突變位置，從而提高突變探測的準確度。在本文數據獲取時段內，發生過4次對觀測數據產生較大影響的地震，分別為2011-03-11日本東北MW9.0地震、2013-04-20蘆山MW7.0地震、2015-04-25和05-12尼泊爾MW7.9和MW7.3地震。參考各地震的發震時間及基準站站點的歷史觀測文件確定起始突變時刻。經過以上處理，本文在GNSS坐標時間序列預處理中設置了289個突跳歷元，共擬合867個突跳項。

圖1 BJFS站突變項擬合前后對比Fig.1 Comparison of time series before and after fitting of mutations in N, E, and U directions at BJFS station

除了突跳項以外，GNSS坐標序列中還會包含長時期趨勢信息，并存在粗差及缺失數據等問題，因此還要對GNSS坐標時間序列進行去除趨勢項、剔除粗差、缺失數據插值補齊等處理。本文利用四分位數粗差探測法對原始坐標時間序列進行粗差探測，將含有粗差的值從原始數據中剔除。對GNSS序列進行PCA分析時，須確保坐標序列是等時間間隔采樣，因此需要根據缺失數據的實際情況選擇合適的插值方法。三次樣條插值方法對于缺失數據較小的時段，其插值效果比較好，但在缺失數據量較大時效果不佳。鑒于此，實驗對缺失小于3 d的序列，利用三次樣條插值法進行補齊；對于缺失3 d以上的序列，采用基于PCA迭代的方法予以補充。具體步驟如下：1)首先用全部有效基準站時間序列的平均值代替缺失部分的數據，組建GNSS坐標時間序列矩陣；2)對組建的矩陣進行PCA分析，選取前3個主分量來恢復缺失的時間序列；3)設置缺失時間序列兩次迭代之差的閾值為10-6，直至全部的缺失時間序列的迭代之差均小于閾值為止。

本文選取224個站點的坐標單日解共538 048個，剔除N、E、U方向粗差值4 925個，缺失數據共有65 759個。圖2為AHBB站N、E、U方向在粗差剔除及缺失時間序列補齊前后的對比，可以看到，未經過處理的序列中分布著較明顯的粗差點和缺失點(矩形標記處)，經過粗差剔除、數據插值等處理后的坐標序列得到明顯改善。

圖2 AHBB站N、E、U方向原始數據、粗差剔除、缺失數據插值補齊時間序列Fig.2 Time series of original data, gross error elimination, and interpolation and completion of missing data in N,E and U directions at AHBB station

2 陸態網坐標時間序列主成分分析

2.1 數據解算及結果分析

對經過預處理后的陸態網GNSS站點分N、E、U方向組建坐標時間序列矩陣并進行PCA分析。圖3、4表示3個方向的主分量特征值分布情況及各方向主分量累積貢獻率。由圖3可見，3個方向特征值的變化曲線相似，其中，N、E方向的特征值差別較小，U方向特征值明顯比N、E方向大。

圖3 N、E、U方向特征值分布Fig.3 Distribution of eigenvalues in N,E and U directions

由圖4可見，N、E、U方向第1、2、3主分量貢獻率分別為31.1%、30.3%、34.8%，12.4%、10.9%、13.7%，3%、6.8%、5.8%，前3個主分量累積貢獻率分別為50.8%、48.0%、54.3%。除了前3個主分量外，其余主分量的貢獻率逐步降低，N、E、U方向貢獻率大于1%的主分量分別有11、12、9個。

圖4 前30個主分量累積貢獻率Fig.4 Cumulative contribution rate of the top 30 principal components in N,E, and U directions

主分量表示的是時間域上的變化，空間特征向量則表示各GNSS基準站的空間響應。為了對N、E、U方向各主分量的空間特征向量進行對比，本文分別對每個方向的前3個主分量對應的空間特征向量進行正則化處理，使每個站點的空間響應大小在-1～1之間(圖5)。

N方向前3個主分量及其對應的空間特征向量如圖5(a)所示，箭頭向上表示正響應，向下表示負響應。3個主分量在時域上沒有表現出單純的隨機性或系統性變化，第1主分量波動較小，第2、3主分量則表現出較為顯著的非線性變化趨勢，第1、2主分量的振幅隨時間推移顯現出增大的趨勢。第1主分量的空間響應分布較為一致，且為正響應，響應大小在區域分布上存在差異；第2、3主分量的空間響應則相對雜亂，響應大小在區域上不存在一致性分布的特征，且正、負響應均存在，正響應主要分布在內蒙和西北地區，負響應主要集中在華中、華南以及西南地區。

E方向前3個主分量及其對應的空間特征向量如圖5(b)所示，時域特征與N方向相似，沒有表現出隨機性或系統性變化。基準站點的第1主分量空間響應上與N方向類似，分布較為一致，西北和東北地區比其他地區的響應要小；第2主分量相比第1主分量要混亂許多，東部地區沒有表現出明顯的空間響應，西部地區的響應也比較弱，但分布較一致，少部分站點表現出較強的響應；第3主分量在西北地區表現出較強的負響應，而在華北、華南和華中地區則表現出正響應分布。

圖5 N、E、U方向第1、2、3主成分分析結果Fig.5 The first, second, and third principal components in N,E, and U directions

U方向前3個主分量及其對應的空間特征向量如圖5(c)所示，時域特征波動性變化相比N、E方向而言，變化的系統性和周期性較清晰。第1主分量的空間響應也表現出比較一致的空間分布；第2主分量的空間響應在西北和西南地區表現出一致的負響應分布，云南地區表現尤為顯著，東北地區為正響應分布；第3主分量的空間特征向量在西北地區表現出較強的負響應分布，華南及東南沿海地區表現出較強的正響應分布。

2.2 共模誤差分析

共模誤差(common mode error, CME)是GNSS坐標時間序列的主要誤差之一，嚴重影響GNSS的解算精度[9]。本文根據各主分量的貢獻率、變化特征及其對應的絕對值化的平均空間響應對共模誤差進行分析。N、E、U方向前10個主分量的貢獻率及對應的平均空間響應情況如圖6、7所示，3個方向第1主分量的空間向量都表現出較為一致的空間響應，第1主分量的貢獻率分別為31.1%、30.3%、34.8%，平均空間響應分別為0.495、0.576、0.537；第2、3主分量的空間響應在局部區域也表現出了一致性分布的特征，兩者的貢獻率也分別達到了12.4%、10.9%、13.7%和7.3%、6.8%、5.8%。從統計結果可以看出，如果僅將第1主分量當作共模誤差或套用其他區域的研究標準不能準確地反映共模誤差的時空特點。鑒于此，可對第1、2、3主分量共同進行共模誤差分析。

圖6 N、E、U方向前10個主分量貢獻率Fig.6 Distribution of the top ten principal components contribution rates in N,E, and U directions

圖7 N、E、U方向前10個主分量正則化后平均空間響應對比Fig.7 Regularized spatial average response of the top ten principal components in N,E, and U directions

2.3 異常站點影響分析

空間響應分布一致的區域，有時候會存在少量異常站點，如華南沿海的GDST站、云南地區的XIAG站、海南的HISY站、江西的JXHK站等，這些站點在區域上會表現出更明顯的空間響應，并且遠大于該區域正常響應的大小，其時間序列也存在異常情況。對于這些異常站點，在沒有剔除的情況下會對PCA的結果產生負面影響。本文剔除絕對值化的平均空間響應大于20%的站點，剔除后，對剩余的196個站點再次進行PCA計算，得到3個方向的主分量貢獻率及站點絕對值化的平均空間響應情況(表1)。

表1 N、E、U方向異常站點剔除后站點絕對值平均空間響應及貢獻率變化(“+”代表提高，“-”代表降低)Tab.1 Changes in the average spatial response and contribution rate of the absolute value of the stations after removing abnormal sites in N,E, and U directions (+ means increase, -means decrease)

由表1可見，剔除異常站點后，N、E、U方向的第1主分量貢獻率分別提升2.0%、3.9%、5.7%，而第2主分量分別下降1.1%、1.9%、6.7%。除前3個主分量以外，其他主分量貢獻率總體變化不大。此外，N、E、U方向的平均空間響應也得到較為明顯的提升，其中第1主分量最為顯著，在3個方向分別提高0.22、0.10、0.24。N、E、U方向的前3個主分量累積空間響應分別提高0.35、0.22、0.34，累積貢獻率分別為提升0.4%、3.2%和下降2.6%，說明異常站點的剔除有助于進一步研究區域站點的空間分布特征。

3 結 語

本文利用主成分分析法對陸態網GNSS坐標時間序進行分析。首先對各站點初始坐標序列進行突變項擬合、粗差剔除、缺失數據插值補齊等預處理，然后分N、E、U方向組建矩陣進行主成分分析。依據主分量時域上的變化特征、空間響應特征以及各個主分量貢獻率等，建議將前3個主分量納入共模誤差分析；經過分析，水儲量變化是引起華北、西北和云南地區空間響應在U方向一致性分布的主要原因；通過異常站點剔除前后的對比發現，異常站點對主成分分析結果影響顯著，剔除異常站點后各方向的平均空間響應和貢獻率在第1主分量上都有明顯提高。