顧及共模誤差的西南極GPS坐標(biāo)時(shí)間序列噪聲分析

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

1 國(guó)家海洋信息中心，天津市六緯路93號(hào)，300171

區(qū)域GPS網(wǎng)因受到地殼運(yùn)動(dòng)、氣候變化、海潮及大氣潮負(fù)荷、解算策略等因素的影響，導(dǎo)致連續(xù)GPS坐標(biāo)時(shí)間序列中存在時(shí)空相關(guān)誤差，稱為共模誤差(common mode error, CME)，準(zhǔn)確提取CME對(duì)于研究噪聲特征具有重要意義[1]。國(guó)內(nèi)外學(xué)者對(duì)CME空間濾波方法進(jìn)行了研究，先后提出帶權(quán)均值法、主分量分析法(PCA)和Karhunen-Loeve展開法，通過比較這3種方法對(duì)GPS坐標(biāo)時(shí)間序列的濾波處理效果發(fā)現(xiàn)，PCA提取共模誤差的效果更優(yōu)[2-3]。另外，明鋒等[4]通過對(duì)比PCA和獨(dú)立分量分析法(ICA)發(fā)現(xiàn)，ICA能夠有效提取區(qū)域GPS觀測(cè)網(wǎng)中的CME。

南極地殼運(yùn)動(dòng)特征是當(dāng)今地球動(dòng)力學(xué)研究的熱點(diǎn)，南極地區(qū)GPS基準(zhǔn)站的建設(shè)和連續(xù)運(yùn)行可為研究坐標(biāo)時(shí)間序列噪聲特征提供數(shù)據(jù)保障。目前對(duì)于GPS坐標(biāo)時(shí)間序列中CME和噪聲特征的研究主要集中在南極半島等區(qū)域，而對(duì)西南極區(qū)域的研究較少。本文利用POLENET計(jì)劃和南極IGS站的GPS觀測(cè)數(shù)據(jù)，對(duì)比分析PCA和ICA提取西南極區(qū)域GPS坐標(biāo)時(shí)間序列中 CME的效果，再采用極大似然估計(jì)法定量估計(jì)濾波前后各站坐標(biāo)時(shí)間序列在不同組合噪聲模型下的噪聲含量，并確定最優(yōu)噪聲模型。

1 GPS站坐標(biāo)殘差時(shí)間序列獲取

1.1 GPS數(shù)據(jù)及處理方法

選取西南極區(qū)域24個(gè)GPS連續(xù)觀測(cè)站，測(cè)站分布如圖1所示，包括10個(gè)IGS站、13個(gè)POLENET網(wǎng)站及中國(guó)南極長(zhǎng)城站(GRW1)，數(shù)據(jù)時(shí)間跨度為2010～2018年，測(cè)站最短時(shí)間跨度為6 a，最長(zhǎng)時(shí)間跨度為9 a。

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

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

1.2 殘差時(shí)間序列預(yù)處理

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

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

(1)

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

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

2 共模誤差提取

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

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

圖3 各PC特征值所占百分比和標(biāo)準(zhǔn)化空間響應(yīng)超過25%的測(cè)站數(shù)Fig.3 The percentage of variance and number of stations with standardized spatial response exceeding 25% derived from PC

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

圖4 各IC標(biāo)準(zhǔn)化空間響應(yīng)超過25%的測(cè)站數(shù)Fig.4 Number of stations with standardized spatialresponse exceeding 25% derived from IC

為比較PCA和ICA方法提取共模誤差的效果，對(duì)剔除共模誤差前后殘差的絕對(duì)值進(jìn)行求和，再除以測(cè)站數(shù)，得到單日L1范數(shù)離散度序列。圖5為2種方法剔除共模誤差前后殘差L1范數(shù)離散度，從圖中可以看出，PCA方法濾波后殘差的振幅略小于ICA方法，表明PCA方法對(duì)西南極區(qū)域GPS數(shù)據(jù)的適應(yīng)性更好。此外，本文還統(tǒng)計(jì)了原始?xì)埐顣r(shí)間序列分別扣除PCA、ICA共模誤差前后各測(cè)站殘差時(shí)間序列的標(biāo)準(zhǔn)差RMS，表1中括號(hào)內(nèi)數(shù)值為RMS的降低率。從表1可以看出，PCA濾波后各站殘差RMS均值在N、E、U分量均比ICA方法小，尤其在U分量差異較大。綜上可知，PCA和ICA空間濾波方法均能有效降低殘差時(shí)間序列的標(biāo)準(zhǔn)差，提高GPS站坐標(biāo)精度；PCA濾波法效果更優(yōu)，是更適合在西南極區(qū)域提取共模誤差的濾波方法。

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

表1 空間濾波前后測(cè)站殘差時(shí)間序列的平均RMS

3 GPS坐標(biāo)時(shí)間序列噪聲特征研究

3.1 噪聲特征分析方法

GPS坐標(biāo)時(shí)間序列噪聲分析方法主要有極大似然估計(jì)法和功率譜分析法，分別從時(shí)間域和頻率域?qū)r(shí)間序列進(jìn)行噪聲分析，以確定噪聲特征。其中，極大似然估計(jì)法不僅可同時(shí)估計(jì)噪聲特征、測(cè)站速度場(chǎng)和周期性振幅等，還可避開頻譜分析需要數(shù)據(jù)均勻采樣的局限性。

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

(2)

AIC=2k+2ln(L)

(3)

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

(4)

式中，k為待估參數(shù)的個(gè)數(shù)。

3.2 共模誤差對(duì)測(cè)站噪聲的影響

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

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

為進(jìn)一步分析各測(cè)站在空間濾波前后不同坐標(biāo)分量的噪聲特征，根據(jù)譜指數(shù)分析結(jié)果及有色噪聲的確定性，假設(shè)坐標(biāo)時(shí)間序列除包含WN外，還含有FN、RWN、PL、廣義高斯-馬爾可夫噪聲(generalized Gauss Markov noise, GGM)及一階自回歸滑動(dòng)平均模型噪聲(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種噪聲模型，根據(jù)極大似然估計(jì)法計(jì)算不同噪聲模型下的AIC和BIC值。由于本文所選樣本較小，選用BIC值作為判斷模型可靠性的指標(biāo)，即BIC值越小結(jié)果越可靠。

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

表2 空間濾波前后各站坐標(biāo)最優(yōu)噪聲模型統(tǒng)計(jì)

此外，CME會(huì)影響時(shí)間序列最優(yōu)噪聲模型的地域分布特征，空間濾波前GPS坐標(biāo)時(shí)間序列最優(yōu)噪聲模型具有較強(qiáng)的區(qū)域分布一致性，空間濾波后區(qū)域分布特征不明顯。以維多利亞地的MCMC、MCM4、SCTB、CRAR、COTE和WHN0站為例，空間濾波前6個(gè)測(cè)站三分量的最優(yōu)噪聲模型以WN+FN為主，濾波后各站不同分量的最優(yōu)噪聲模型按比例大小依次為WN+FN、GGM+RWN等，其原因?yàn)榭臻g濾波會(huì)大幅降低GPS坐標(biāo)時(shí)間序列的白噪聲和閃爍噪聲，當(dāng)閃爍噪聲占主導(dǎo)地位時(shí)，容易掩蓋其他有色噪聲。這同時(shí)也表明，對(duì)GPS坐標(biāo)時(shí)間序列進(jìn)行分析時(shí)，利用空間濾波扣除CME具有必要性。

3.3 共模誤差對(duì)測(cè)站速度影響分析

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

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

4 結(jié) 語

通過分析不同空間濾波方法和濾波前后西南極區(qū)域24個(gè)GPS站的連續(xù)坐標(biāo)時(shí)間序列，得到以下結(jié)論：

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

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

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

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