劉曉祥 高二濤 羅 益 付波霖

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

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

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

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

1 陸態(tài)網(wǎng)坐標(biāo)時(shí)間序列預(yù)處理

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

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

1.2 GNSS站點(diǎn)時(shí)間序列預(yù)處理

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

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

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

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

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

圖2 AHBB站N、E、U方向原始數(shù)據(jù)、粗差剔除、缺失數(shù)據(jù)插值補(bǔ)齊時(shí)間序列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)坐標(biāo)時(shí)間序列主成分分析

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

對(duì)經(jīng)過(guò)預(yù)處理后的陸態(tài)網(wǎng)GNSS站點(diǎn)分N、E、U方向組建坐標(biāo)時(shí)間序列矩陣并進(jìn)行PCA分析。圖3、4表示3個(gè)方向的主分量特征值分布情況及各方向主分量累積貢獻(xiàn)率。由圖3可見(jiàn)，3個(gè)方向特征值的變化曲線相似，其中，N、E方向的特征值差別較小，U方向特征值明顯比N、E方向大。

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

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

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

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

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

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

圖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個(gè)主分量及其對(duì)應(yīng)的空間特征向量如圖5(c)所示，時(shí)域特征波動(dòng)性變化相比N、E方向而言，變化的系統(tǒng)性和周期性較清晰。第1主分量的空間響應(yīng)也表現(xiàn)出比較一致的空間分布；第2主分量的空間響應(yīng)在西北和西南地區(qū)表現(xiàn)出一致的負(fù)響應(yīng)分布，云南地區(qū)表現(xiàn)尤為顯著，東北地區(qū)為正響應(yīng)分布；第3主分量的空間特征向量在西北地區(qū)表現(xiàn)出較強(qiáng)的負(fù)響應(yīng)分布，華南及東南沿海地區(qū)表現(xiàn)出較強(qiáng)的正響應(yīng)分布。

2.2 共模誤差分析

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

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

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

2.3 異常站點(diǎn)影響分析

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

表1 N、E、U方向異常站點(diǎn)剔除后站點(diǎn)絕對(duì)值平均空間響應(yīng)及貢獻(xiàn)率變化(“+”代表提高，“-”代表降低)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可見(jiàn)，剔除異常站點(diǎn)后，N、E、U方向的第1主分量貢獻(xiàn)率分別提升2.0%、3.9%、5.7%，而第2主分量分別下降1.1%、1.9%、6.7%。除前3個(gè)主分量以外，其他主分量貢獻(xiàn)率總體變化不大。此外，N、E、U方向的平均空間響應(yīng)也得到較為明顯的提升，其中第1主分量最為顯著，在3個(gè)方向分別提高0.22、0.10、0.24。N、E、U方向的前3個(gè)主分量累積空間響應(yīng)分別提高0.35、0.22、0.34，累積貢獻(xiàn)率分別為提升0.4%、3.2%和下降2.6%，說(shuō)明異常站點(diǎn)的剔除有助于進(jìn)一步研究區(qū)域站點(diǎn)的空間分布特征。

3 結(jié) 語(yǔ)

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