顧及共模誤差的西南極GPS坐標時間序列噪聲分析

王朝陽 邢 喆 張 峰 楊曉彤 王力彥

1 國家海洋信息中心，天津市六緯路93號，300171

區域GPS網因受到地殼運動、氣候變化、海潮及大氣潮負荷、解算策略等因素的影響，導致連續GPS坐標時間序列中存在時空相關誤差，稱為共模誤差(common mode error, CME)，準確提取CME對于研究噪聲特征具有重要意義[1]。國內外學者對CME空間濾波方法進行了研究，先后提出帶權均值法、主分量分析法(PCA)和Karhunen-Loeve展開法，通過比較這3種方法對GPS坐標時間序列的濾波處理效果發現，PCA提取共模誤差的效果更優[2-3]。另外，明鋒等[4]通過對比PCA和獨立分量分析法(ICA)發現，ICA能夠有效提取區域GPS觀測網中的CME。

南極地殼運動特征是當今地球動力學研究的熱點，南極地區GPS基準站的建設和連續運行可為研究坐標時間序列噪聲特征提供數據保障。目前對于GPS坐標時間序列中CME和噪聲特征的研究主要集中在南極半島等區域，而對西南極區域的研究較少。本文利用POLENET計劃和南極IGS站的GPS觀測數據，對比分析PCA和ICA提取西南極區域GPS坐標時間序列中 CME的效果，再采用極大似然估計法定量估計濾波前后各站坐標時間序列在不同組合噪聲模型下的噪聲含量，并確定最優噪聲模型。

1 GPS站坐標殘差時間序列獲取

1.1 GPS數據及處理方法

選取西南極區域24個GPS連續觀測站，測站分布如圖1所示，包括10個IGS站、13個POLENET網站及中國南極長城站(GRW1)，數據時間跨度為2010～2018年，測站最短時間跨度為6 a，最長時間跨度為9 a。

圖1 西南極24個GPS觀測站分布Fig.1 Distribution of 24 GPS stations inwest Antarctica region

采用Bernese 5.2軟件雙差定位程序對GPS數據進行處理，獲得各測站單日解坐標。為將GPS基準站坐標精確統一到ITRF框架下，在解算過程中引入南極大陸其他5個IGS站(SYOG、MAW1、DAV1 、CAS1和DUM1)同時段的觀測數據。在數據處理時，采用IGS歐洲定軌中心(CODE)提供的精密星歷和地球定向參數等產品，整周模糊度采用無電離層組合模糊度固定(QIF)策略，基線組成方法采用OBS-MAX，衛星截止高度角為7°，采樣間隔為30 s，對流層改正采用Saastamoinen模型，映射函數選用VMF3模型，潮汐改正選用FES2014b模型，采用ITRF2014框架。

1.2 殘差時間序列預處理

受儀器設備故障及外界環境因素影響，GPS坐標時間序列中不可避免地存在數據缺失和粗差，導致分析結果不準確，因此有必要對單日坐標時間序列進行粗差剔除和數據插值。首先采用多通道奇異譜分析和四分位距相結合的方法(MSSA-IQR)[5]進行異常值探測和剔除，然后從各站坐標時間序列中剔除線性趨勢和周期項，得到殘差時間序列，計算公式為：

υi=y(ti)-(a+bti+csin(2πti)+dcos(2πti)+

(1)

式中，υi為第i個歷元的觀測殘差，y(ti)為測站坐標時間序列，ti為以a為單位的歷元時間，a、b為測站坐標的初始歷元和線性速度，c、d為時間序列的年周期系數，e、f為半年周期系數，g為歷元Tg處的階躍式偏移量，H為Heaviside階梯函數。最后采用基于數據驅動的RegEM算法[6]進行殘差時間序列插值，RegEM數據插值法可考慮站點坐標時間序列的物理背景，并顧及各站之間的相關性。圖2為RegEM算法插值后的FALL站和WHN0站殘差時間序列，圖中藍色為原始數據，紅色為RegEM插值結果。

圖2 FALL站和WHN0站RegEM插值結果Fig.2 The RegEM interpolation results of FALL and WHN0 station

2 共模誤差提取

本文引入KMO檢驗[7]以驗證各站坐標殘差時間序列是否適用于PCA和ICA。KMO通過檢驗各站殘差時間序列之間的簡單相關系數和偏相關系數來確定數據是否適合進行因子分析，KMO值越接近1，說明變量間的相關性越強，數據越適合進行因子分析；KMO值接近0，則說明數據不適合進行因子分析。結果表明，各站坐標殘差時間序列在N、E、U分量上的KMO值分別為0.804、0.824、0.888，通過Barttest檢驗顯著性水平為0.05的結果，數據可用于主成分分析。

首先利用PCA方法對N、E、U分量的殘差時間序列進行濾波處理，得到各主分量(PC)特征值占總特征值的百分比，結果如圖3(a)所示。從圖中可以看出，N分量前3個PC的貢獻率分別為25.75%、17.50%、7.37%，E分量前3個PC的貢獻率分別為24.97%、12.85%、7.20%，U分量前3個PC的貢獻率分別為33.24%、11.92%、6.99%。根據Dong等[8]對共模誤差的定義，要求參與計算的主分量中至少有50%的測站標準化空間響應超過25%，且該分量的特征值超過所有特征值總和的1%，據此對各主分量的標準化空間響應進行顯著性檢驗，統計各PC的標準化空間響應大于25%的測站數，結果如圖3(b)所示。從圖中可以看出，N、E分量第1、2主分量的標準化空間響應大于25%的測站超過50%，即前2個PC可作為區域的CME進行計算，U分量前3個PC滿足共模誤差要求。

圖3 各PC特征值所占百分比和標準化空間響應超過25%的測站數Fig.3 The percentage of variance and number of stations with standardized spatial response exceeding 25% derived from PC

采用FastICA方法分別對N、E、U分量的殘差時間序列進行ICA分解，得到各獨立分量(IC)及其標準化空間響應。FastICA算法以各估計的信號源之間的互信息最小化為目標函數，采用基于牛頓迭代的固定點優化方案，其計算效率比傳統的梯度法快10～100倍，且穩健性更好[8]。與PCA方法不同的是，ICA方法無法確定各IC之間的次序及各IC的特征值，因此本文分別統計N、E、U分量各IC的標準化空間響應大于25%的測站，并按照測站數大小對IC進行排列，結果如圖4所示。從圖中可以看出，N、E分量前2個IC的標準化空間響應大于25%的測站數剛好超過50%，U分量前3個IC滿足共模誤差要求，因此水平分量選擇前2個IC計算CME，U分量選擇前3個IC計算CME。

圖4 各IC標準化空間響應超過25%的測站數Fig.4 Number of stations with standardized spatialresponse exceeding 25% derived from IC

為比較PCA和ICA方法提取共模誤差的效果，對剔除共模誤差前后殘差的絕對值進行求和，再除以測站數，得到單日L1范數離散度序列。圖5為2種方法剔除共模誤差前后殘差L1范數離散度，從圖中可以看出，PCA方法濾波后殘差的振幅略小于ICA方法，表明PCA方法對西南極區域GPS數據的適應性更好。此外，本文還統計了原始殘差時間序列分別扣除PCA、ICA共模誤差前后各測站殘差時間序列的標準差RMS，表1中括號內數值為RMS的降低率。從表1可以看出，PCA濾波后各站殘差RMS均值在N、E、U分量均比ICA方法小，尤其在U分量差異較大。綜上可知，PCA和ICA空間濾波方法均能有效降低殘差時間序列的標準差，提高GPS站坐標精度；PCA濾波法效果更優，是更適合在西南極區域提取共模誤差的濾波方法。

圖5 PCA、ICA剔除共模誤差后的L1范數離散度序列比較Fig.5 Comparison ofthe L1-norm dispersion series after CME elimination by PCA and ICA

表1 空間濾波前后測站殘差時間序列的平均RMS

3 GPS坐標時間序列噪聲特征研究

3.1 噪聲特征分析方法

GPS坐標時間序列噪聲分析方法主要有極大似然估計法和功率譜分析法，分別從時間域和頻率域對時間序列進行噪聲分析，以確定噪聲特征。其中，極大似然估計法不僅可同時估計噪聲特征、測站速度場和周期性振幅等，還可避開頻譜分析需要數據均勻采樣的局限性。

本文采用Hector軟件改進的極大似然估計法分析GPS坐標時間序列的噪聲特征，該軟件具有更高的數據處理速率，相比于CATS軟件可選擇更多的噪聲模型。在數據處理過程中，顧及噪聲模型參數對估計結果的影響，采用赤池信息量(Akaike information criteria, AIC)或貝葉斯信息量(Bayesian information criteria, BIC)最小的噪聲模型估計準則，以獲得更加穩健的噪聲模型估計結果。使用AIC和BIC標準對不同噪聲模型進行分析的公式為：

(2)

AIC=2k+2ln(L)

(3)

BIC=kln(N)+2ln(L)

(4)

式中，k為待估參數的個數。

3.2 共模誤差對測站噪聲的影響

利用PCA空間濾波法提取測站的共模誤差，分析濾波前后時間序列的最優噪聲模型，以確定共模誤差對較長時間序列噪聲的影響。時間序列中的噪聲譜可用指數定律來描述，通過計算譜指數u可以判斷噪聲的大概類型。當u=0時對應白噪聲，當u=-1時對應閃爍噪聲，當u=-2時對應隨機游走噪聲[9]。圖6為濾波后的西南極GPS坐標時間序列譜指數，從圖中可以看出，西南極24個GPS基準站各坐標分量的噪聲譜指數均介于-1～-0.5之間，說明這些測站的噪聲類型均不具有純白噪聲特征，而是白噪聲與有色噪聲的組合。

圖6 空間濾波后測站3個方向的譜指數Fig.6 Spectral index values in the three directions after spatial filtering

為進一步分析各測站在空間濾波前后不同坐標分量的噪聲特征，根據譜指數分析結果及有色噪聲的確定性，假設坐標時間序列除包含WN外，還含有FN、RWN、PL、廣義高斯-馬爾可夫噪聲(generalized Gauss Markov noise, GGM)及一階自回歸滑動平均模型噪聲(auto-regressive moving average model noise(1,1), ARMA(1))，選取WN、WN+FN、WN+RWN、WN+FN+RWN、WN+PL、WN+ARMA(1)、GGM+FN及GGM+RWN等8種噪聲模型，根據極大似然估計法計算不同噪聲模型下的AIC和BIC值。由于本文所選樣本較小，選用BIC值作為判斷模型可靠性的指標，即BIC值越小結果越可靠。

使用上述8種模型計算空間濾波前后24個GPS站N、E、U三分量的BIC值，然后確定每個測站不同坐標分量的最佳噪聲模型，結果見表2。從表中可以看出，西南極GPS站3個坐標分量時間序列的噪聲特征并不完全相同，濾波后部分測站的噪聲特征發生改變。空間濾波前，N分量的最優噪聲模型比例最大的為WN+FN(45.83%)，其次為WN+ARMA(1)(29.17%)；E分量的最優噪聲模型比例最大的為WN+FN(41.67%)，其次為WN+ARMA(1)(33.33%)；U分量的最優噪聲模型比例最大的為WN+FN(48.00%)，其次為GGM+FN(16.00%)和GGM+RWN(16.00%)。空間濾波后，E分量的最優噪聲模型WN+FN和WN+ARMA(1)均占33.33%；N分量的最優噪聲模型比例最大的為WN+ARMA(1)(41.76%)，其次為WN+FN(25.00%)；U分量的最優噪聲模型比例最大的為WN+FN(45.83%)，其次為WN+ARMA(1)(25.00%)和GGM+FN(16.67%)。通過對比濾波前后各測站時間序列的最優噪聲模型可知，共模誤差同時具有白噪聲和有色噪聲的特征。

表2 空間濾波前后各站坐標最優噪聲模型統計

此外，CME會影響時間序列最優噪聲模型的地域分布特征，空間濾波前GPS坐標時間序列最優噪聲模型具有較強的區域分布一致性，空間濾波后區域分布特征不明顯。以維多利亞地的MCMC、MCM4、SCTB、CRAR、COTE和WHN0站為例，空間濾波前6個測站三分量的最優噪聲模型以WN+FN為主，濾波后各站不同分量的最優噪聲模型按比例大小依次為WN+FN、GGM+RWN等，其原因為空間濾波會大幅降低GPS坐標時間序列的白噪聲和閃爍噪聲，當閃爍噪聲占主導地位時，容易掩蓋其他有色噪聲。這同時也表明，對GPS坐標時間序列進行分析時，利用空間濾波扣除CME具有必要性。

3.3 共模誤差對測站速度影響分析

根據前文分析得到的空間濾波前后西南極24個測站各分量時間序列的最優噪聲模型，采用極大似然估計法估計濾波前后各分量在ITRF2014框架下的速度及其不確定度，結果見表3(單位mm/a)。從表中可以看出，經過PCA濾波后，87.50%(21個)測站的水平速度變化在±0.2 mm/a，79.17%(19個)測站的垂直速度變化在±0.3 mm/a，說明空間濾波對時間序列速度估計值的影響較小；N、E、U三分量的速度不確定度分別平均降低14.56%、20.66%和36.40%，說明垂直方向通過空間濾波扣除的共模誤差量級大于水平方向，垂直方向不確定度的減小幅度更大。時間序列較短或非線性運動幅度較大的測站差異較大，因為這些測站的速度更易受系統誤差的影響。

表3 空間濾波前后最優噪聲模型下各站速度及不確定度

4 結 語

通過分析不同空間濾波方法和濾波前后西南極區域24個GPS站的連續坐標時間序列，得到以下結論：

1)PCA和ICA空間濾波方法均能有效提取區域GPS坐標時間序列中的共模誤差，但PCA方法的濾波效果略優于ICA方法，且在垂直方向的濾波效果優于水平方向。

2)利用PCA濾波后的GPS站各坐標分量的噪聲譜指數均介于-1～-0.5之間，表明這些測站的噪聲是包括白噪聲在內的多種噪聲模型的組合。

3)GPS站各坐標分量時間序列的最優噪聲模型以WN+FN、WN+ARMA(1)為主，其中多個測站的U分量還含有GGM。去除共模誤差后，部分測站的最優噪聲模型發生改變，但仍以WN+FN、WN+ARMA(1)為主，最優噪聲的區域分布特征不明顯。

4)空間濾波前后GPS站各坐標分量的速度平均變化均小于0.1 mm/a，速度不確定度分別平均降低14.56%、20.66%和36.40%，說明共模誤差對時間序列速度估計值的影響較小，但能有效降低白噪聲和有色噪聲的量級，從而大幅減小線性項的不確定性。