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

劉曉祥 高二濤 羅 益 付波霖

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

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

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

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

1 陸態(tài)網(wǎng)坐標時間序列預處理

1.1 數(shù)據(jù)選取

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

1.2 GNSS站點時間序列預處理

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

圖1為BJFS站在去除突跳項前后N、E、U方向的對比，可以看到，通過擬合處理后，突跳位置得到有效補充。此外對于因儀器設(shè)備更換而引起的突變，可以通過各個站點的log文件獲??；對于其他未知原因造成的突變，為了確保該類突變的擬合效果，通過目視排查獲取突變位置，從而提高突變探測的準確度。在本文數(shù)據(jù)獲取時段內(nèi)，發(fā)生過4次對觀測數(shù)據(jù)產(chǎn)生較大影響的地震，分別為2011-03-11日本東北MW9.0地震、2013-04-20蘆山MW7.0地震、2015-04-25和05-12尼泊爾MW7.9和MW7.3地震。參考各地震的發(fā)震時間及基準站站點的歷史觀測文件確定起始突變時刻。經(jīng)過以上處理，本文在GNSS坐標時間序列預處理中設(shè)置了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坐標序列中還會包含長時期趨勢信息，并存在粗差及缺失數(shù)據(jù)等問題，因此還要對GNSS坐標時間序列進行去除趨勢項、剔除粗差、缺失數(shù)據(jù)插值補齊等處理。本文利用四分位數(shù)粗差探測法對原始坐標時間序列進行粗差探測，將含有粗差的值從原始數(shù)據(jù)中剔除。對GNSS序列進行PCA分析時，須確保坐標序列是等時間間隔采樣，因此需要根據(jù)缺失數(shù)據(jù)的實際情況選擇合適的插值方法。三次樣條插值方法對于缺失數(shù)據(jù)較小的時段，其插值效果比較好，但在缺失數(shù)據(jù)量較大時效果不佳。鑒于此，實驗對缺失小于3 d的序列，利用三次樣條插值法進行補齊；對于缺失3 d以上的序列，采用基于PCA迭代的方法予以補充。具體步驟如下：1)首先用全部有效基準站時間序列的平均值代替缺失部分的數(shù)據(jù)，組建GNSS坐標時間序列矩陣；2)對組建的矩陣進行PCA分析，選取前3個主分量來恢復缺失的時間序列；3)設(shè)置缺失時間序列兩次迭代之差的閾值為10-6，直至全部的缺失時間序列的迭代之差均小于閾值為止。

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

圖2 AHBB站N、E、U方向原始數(shù)據(jù)、粗差剔除、缺失數(shù)據(jù)插值補齊時間序列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 陸態(tài)網(wǎng)坐標時間序列主成分分析

2.1 數(shù)據(jù)解算及結(jié)果分析

對經(jīng)過預處理后的陸態(tài)網(wǎng)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個主分量在時域上沒有表現(xiàn)出單純的隨機性或系統(tǒng)性變化，第1主分量波動較小，第2、3主分量則表現(xiàn)出較為顯著的非線性變化趨勢，第1、2主分量的振幅隨時間推移顯現(xiàn)出增大的趨勢。第1主分量的空間響應分布較為一致，且為正響應，響應大小在區(qū)域分布上存在差異；第2、3主分量的空間響應則相對雜亂，響應大小在區(qū)域上不存在一致性分布的特征，且正、負響應均存在，正響應主要分布在內(nèi)蒙和西北地區(qū)，負響應主要集中在華中、華南以及西南地區(qū)。

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

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

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

2.2 共模誤差分析

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

空間響應分布一致的區(qū)域，有時候會存在少量異常站點，如華南沿海的GDST站、云南地區(qū)的XIAG站、海南的HISY站、江西的JXHK站等，這些站點在區(qū)域上會表現(xiàn)出更明顯的空間響應，并且遠大于該區(qū)域正常響應的大小，其時間序列也存在異常情況。對于這些異常站點，在沒有剔除的情況下會對PCA的結(jié)果產(chǎn)生負面影響。本文剔除絕對值化的平均空間響應大于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%，說明異常站點的剔除有助于進一步研究區(qū)域站點的空間分布特征。

3 結(jié) 語

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