顧及空間相關(guān)性的GNSS 坐標(biāo)序列插值比較

謝春橋，匡翠林

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083）

0 引言

在全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system, GNSS）連續(xù)運(yùn)行參考站的長(zhǎng)期觀測(cè)過(guò)程中，由于接收機(jī)與天線的故障、系統(tǒng)供電中斷、粗差的剔除以及其他未知原因，導(dǎo)致GNSS 坐標(biāo)序列易存在數(shù)據(jù)缺失。數(shù)據(jù)的缺失會(huì)帶來(lái)諸多不利影響，例如在空域分析中，由于區(qū)域網(wǎng)GNSS坐標(biāo)序列普遍受到共模誤差（common mode error,CME）的影響，而共模誤差的計(jì)算一般要求所有站點(diǎn)的數(shù)據(jù)必須連續(xù)，數(shù)據(jù)缺失將導(dǎo)致無(wú)法計(jì)算并剔除CME；在時(shí)域分析中，若GNSS 坐標(biāo)序列數(shù)據(jù)缺失過(guò)多，會(huì)影響模型參數(shù)估計(jì)值的精度，如速度不確定性和噪聲幅度等參數(shù)[1]；在頻域分析中，數(shù)據(jù)缺失即采樣率降低，會(huì)出現(xiàn)混頻現(xiàn)象。對(duì)于數(shù)據(jù)缺失問(wèn)題，國(guó)內(nèi)外學(xué)者通常選取數(shù)據(jù)缺失比例較小的站點(diǎn)進(jìn)行研究，而對(duì)于數(shù)據(jù)缺失比例高，尤其是數(shù)據(jù)連續(xù)缺失的站點(diǎn)，一般采用舍棄的策略，或截取其中缺失率低的部分?jǐn)?shù)據(jù)，這將導(dǎo)致GNSS坐標(biāo)序列數(shù)據(jù)不能得到充分利用。為解決數(shù)據(jù)缺失所引起的問(wèn)題，必須對(duì)缺失數(shù)據(jù)進(jìn)行插值，得到完整連續(xù)的數(shù)據(jù)，并盡可能恢復(fù)數(shù)據(jù)原有的噪聲特性。

目前，單站坐標(biāo)序列插值的方法有拉格朗日插值、多點(diǎn) 3 次樣條插值[2]和正交多項(xiàng)式插值[3]等，但這些插值方法只適用于缺失較少的情形，且未顧及站點(diǎn)之間的空間相關(guān)性?；诙嗾镜牟逯捣椒櫦傲苏军c(diǎn)間的空間相關(guān)性，如常見(jiàn)的多通道奇異譜（multi-channel singular spectrum analysis,MSSA）插值方法，MSSA 插值方法對(duì)多站坐標(biāo)序列通過(guò)構(gòu)建時(shí)滯軌跡矩陣進(jìn)行分解和重構(gòu)，提取出原序列區(qū)域性的周期信號(hào)和趨勢(shì)信號(hào)等，再通過(guò)迭代和交叉驗(yàn)證進(jìn)行插值，該方法無(wú)需先驗(yàn)信息，得到了廣泛應(yīng)用[4]。文獻(xiàn)[5-6]將 MSSA 應(yīng)用到 GNSS 坐標(biāo)序列缺失數(shù)據(jù)的插值中，取得了較好的插值效果，但其插值后的數(shù)據(jù)較為光滑，無(wú)法保留原序列的高頻信息。而克里金卡爾曼濾波（Kriged Kalman filter, KKF）方法能夠保留序列的高頻信息，該方法是1 種時(shí)空插值方法，通過(guò)構(gòu)建空間場(chǎng)，利用卡爾曼（Kalman）濾波，采用期望最大化（expectation maximization, EM）算法，估計(jì)Kalman 濾波初值和相關(guān)參數(shù)，基于克里金模型插補(bǔ)缺失數(shù)據(jù)[7]。文獻(xiàn)[8]在此基礎(chǔ)上，提出了抗差KKF 插值方法，并將其應(yīng)用到 GNSS 坐標(biāo)序列插值中，結(jié)果表明，對(duì)于空間測(cè)站分布稀疏、數(shù)據(jù)連續(xù)缺失較多的坐標(biāo)序列，KKF 仍能較好地插補(bǔ)缺失數(shù)據(jù)，且較好地保留了細(xì)節(jié)信息[8-9]。此外，文獻(xiàn)[10]提出 1 種正則期望最大化（regularized EM,RegEM）算法的插值方法，該方法利用不完整數(shù)據(jù)集，采用嶺估計(jì)和廣義交叉核實(shí)等方法，通過(guò)迭代估計(jì)均值和協(xié)方差矩陣，用得到的均值和協(xié)方差陣對(duì)缺失數(shù)據(jù)進(jìn)行插值，該方法也顧及了空間相關(guān)性，插值精度高，在氣象領(lǐng)域應(yīng)用廣泛。文獻(xiàn)[11]將 RegEM 算法引入到 GNSS 坐標(biāo)序列的插值中，證實(shí)了該方法插值精度高，且對(duì)于缺失比例較大的GNSS 坐標(biāo)序列，仍能較好地復(fù)現(xiàn)其細(xì)節(jié)信息。本文利用基于多站的 KKF 和 RegEM 插值方法，對(duì)GNSS 殘差坐標(biāo)序列進(jìn)行插補(bǔ)，同時(shí)與常用的 MSSA 插值方法進(jìn)行對(duì)比，分析這 3 種顧及空間相關(guān)性的插值方法的效果，比較插值結(jié)果對(duì)GNSS 坐標(biāo)序列噪聲特性的影響。

1 RegEM 和 KKF 插值算法

由于GNSS 坐標(biāo)序列中含有線性項(xiàng)和周期項(xiàng)，同時(shí)還存在粗差以及階躍，粗差與階躍會(huì)影響插值結(jié)果。因此，本文采用經(jīng)過(guò)四分位距（interquartile range, IQR）方法剔除粗差并修復(fù)階躍后的殘差坐標(biāo)序列，再進(jìn)行插值[12-13]。

1.1 KKF 插值算法

設(shè)歷元t的狀態(tài)量為tα，GNSS 坐標(biāo)序列非缺失情況下的觀測(cè)值和觀測(cè)噪聲分別為Xa|t、εa|t，而為缺失情況下的觀測(cè)值和觀測(cè)噪聲，Ha和Hm為非缺失和缺失情況下的空間場(chǎng)，可由克里金（Kriging）模型確定，則歷元t的觀測(cè)方程可表示為

在計(jì)算過(guò)程中，首先對(duì)Xm|t向量和Hm矩陣賦初值0，則設(shè)Q為觀測(cè)噪聲方差－協(xié)方差矩陣；R為狀態(tài)噪聲的等價(jià)觀測(cè)方差－協(xié)方差陣；Φ為狀態(tài)轉(zhuǎn)移矩陣；P為估計(jì)的協(xié)方差陣，利用卡爾曼濾波（式（2）～式（6））可求得tα：

1）狀態(tài)一步預(yù)測(cè)為

2）一步預(yù)測(cè)協(xié)方差陣為

3）濾波增益為

4）狀態(tài)估計(jì)為

5）估計(jì)的協(xié)方差陣為

再利用 Kriging 模型求當(dāng)前歷元缺失觀測(cè)值的空間場(chǎng)Hm，最后計(jì)算當(dāng)前歷元缺失的觀測(cè)值為

式中，通過(guò)EM 算法確定Kalman 濾波初值及相關(guān)參數(shù)與Q等。關(guān)于 KKF 插值更詳細(xì)的描述可參考文獻(xiàn)[7-9]。

1.2 RegEM 插值算法

設(shè)1 個(gè)含有n個(gè)歷元、p個(gè)測(cè)站（n＞p）的觀測(cè)值矩陣Z，其元素zij表示測(cè)站j（j=1,2,…,p）在歷元ti（i=1,2,…,n）時(shí)的觀測(cè)值，設(shè)Z的待估均值為μ（1×p維），協(xié)方差陣為Σ（p·p維），記某1 歷元ti由a個(gè)非缺失數(shù)據(jù)構(gòu)成向量，其均值和協(xié)方差分別為和而其余m個(gè)缺失數(shù)據(jù)構(gòu)成向量Zm（1×m維），其均值和協(xié)方差分別為μm（1×m維）和Σm（m·m維），Za和Zm的交叉協(xié)方差矩陣為或?qū)τ赯中每1 個(gè)歷元的觀測(cè)值，缺失觀測(cè)值向量可通過(guò)非缺失觀測(cè)值的線性回歸模型表示為

式中：B（a·m維）為回歸系數(shù)，可通過(guò)嶺回歸求得；殘差e為一均值為0、協(xié)方差陣為C（m·m維）的未知隨機(jī)變量。

在RegEM 算法中，通常先計(jì)算均值μ（忽略缺失數(shù)據(jù)），再將缺失數(shù)據(jù)初值置為 0，并求協(xié)方差陣Σ，通過(guò)條件極大似然估計(jì)求出Z中每 1 歷元包含數(shù)據(jù)缺失的觀測(cè)值的回歸系數(shù)的估值?B和殘差協(xié)方差陣的估值等，則有：

求出缺失數(shù)據(jù)估值的公式為

插補(bǔ)完成后，需對(duì)插值后的數(shù)據(jù)進(jìn)行更新，重新計(jì)算均值和協(xié)方差陣，再次進(jìn)行迭代計(jì)算，當(dāng)估計(jì)的缺失值的均值和協(xié)方差達(dá)到指定的閾值時(shí)，便停止迭代。關(guān)于RegEM 算法更詳細(xì)的推導(dǎo)見(jiàn)參考文獻(xiàn)[10]。

1.3 插值精度評(píng)價(jià)指標(biāo)

若人為地將 GNSS 坐標(biāo)序列數(shù)據(jù)移除后再進(jìn)行插值，由于待插值點(diǎn)的真值已知，可采用平均絕對(duì)值誤差（mean absolute error, MAE）、皮爾遜（Pearson）系數(shù)以及斯皮爾曼（Spearman）秩相關(guān)系數(shù)等指標(biāo)評(píng)價(jià)插值精度。記某序列插補(bǔ)值總數(shù)為s，fk、rk分別為第k（k=1,2,…,s）個(gè)插補(bǔ)值和真實(shí)值，fk、r k按由大到小的順序排列后的秩次為則MAE、Pearson 系數(shù)ρ1和Spearman 系數(shù)ρ2的計(jì)算公式為：

在這些指標(biāo)中，MAE 為絕對(duì)指標(biāo)，Pearson 系數(shù)和Spearman 系數(shù)為相對(duì)指標(biāo)。若插值算法效果越好，則MAE 值越小，Pearson 系數(shù)和Spearman系數(shù)越趨近于1。

在對(duì)實(shí)際缺失數(shù)據(jù)進(jìn)行插值時(shí)，由于缺失數(shù)據(jù)的真值未知，不能采用上述指標(biāo)。由于插值效果好的算法應(yīng)盡可能還原出原序列的噪聲水平，所保留的方差也較大。通常在利用主成分分析（principle components analysis, PCA）提取共模誤差時(shí)，將前幾個(gè)方差較大的主成分代表區(qū)域網(wǎng)的共模誤差[14]。故本文利用前幾個(gè)主成分貢獻(xiàn)率之和這一指標(biāo)評(píng)價(jià)插值效果[11]。此外，插值會(huì)使得GNSS 坐標(biāo)序列噪聲受到影響，時(shí)間序列噪聲分析通常利用極大似然估計(jì)（maximum likelihood estimate, MLE）方法，該方法可同時(shí)估計(jì)各噪聲模型下的噪聲幅度、周期項(xiàng)、站點(diǎn)速度及不確定性等。由該方法計(jì)算得到的 MLE 值是噪聲水平和速度不確定性等的綜合反映，也是評(píng)價(jià)最優(yōu)噪聲模型的依據(jù)[15]，因此也可通過(guò)插值前后 GNSS 坐標(biāo)序列的 MLE 值來(lái)評(píng)價(jià)插值效果。MLE 值的計(jì)算式為

2 GNSS 坐標(biāo)序列插值

2.1 模擬插值實(shí)驗(yàn)

本實(shí)驗(yàn)采用中國(guó)陸態(tài)網(wǎng)數(shù)據(jù)缺失率較低的NMDW 站及其周邊9 個(gè)站點(diǎn)的坐標(biāo)序列，站點(diǎn)分布如圖 1 所示，其中 NMDW 站數(shù)據(jù)缺失率為3.06 %，該站有1 段800 多天完全連續(xù)數(shù)據(jù)。為研究不同比例連續(xù)缺失下，各插值方法的效果，針對(duì) NMDW 站各方向連續(xù)段第 20～800 個(gè)數(shù)據(jù)，以20 個(gè)數(shù)據(jù)為步長(zhǎng)，人為地移除數(shù)據(jù)，從而構(gòu)成40 組插值實(shí)驗(yàn)數(shù)據(jù)，再與周邊9 個(gè)站點(diǎn)的數(shù)據(jù)進(jìn)行插值。

圖1 GNSS 實(shí)驗(yàn)站點(diǎn)分布

限于篇幅，此處僅給出連續(xù)移除100 個(gè)數(shù)據(jù)情形下的插值時(shí)序結(jié)果，圖2 表示NMDW 站N、E、U方向上的 GNSS 殘差坐標(biāo)序列，其中菱形點(diǎn)表示為模擬數(shù)據(jù)缺失而連續(xù)移除的數(shù)據(jù)點(diǎn)，實(shí)心圓點(diǎn)表示用于后續(xù)插值而保留的數(shù)據(jù)點(diǎn)。圖 3～圖 5為 3 種插值方法插補(bǔ)的結(jié)果，并與原殘差坐標(biāo)序列進(jìn)行對(duì)比，由此可以看出，MSSA 插值結(jié)果能保證數(shù)據(jù)的連續(xù)性和周期特性，但過(guò)于光滑，無(wú)法復(fù)現(xiàn)缺失數(shù)據(jù)的噪聲特性。而 RegEM 和 KKF 插值均能保持 GNSS 坐標(biāo)序列總體趨勢(shì)，能夠較好地還原坐標(biāo)序列的噪聲水平，但RegEM 插值在細(xì)節(jié)方面處理得更好。

圖2 NMDW 站GNSS 殘差坐標(biāo)序列

圖3 原殘差序列與MSSA 插值結(jié)果比較

圖6 為NMDW 站實(shí)驗(yàn)中，不同缺失比例下插值后求得的MAE、Pearson 系數(shù)及Spearman 系數(shù)，其中Pearson 系數(shù)及Spearman 系數(shù)為相對(duì)指標(biāo)。從圖6 可以得出：

圖4 原殘差序列與KKF 插值結(jié)果比較

圖5 原殘差序列與RegEM 插值結(jié)果比較

對(duì)于不同比例的連續(xù)缺失情形，RegEM 和KKF插值后的各統(tǒng)計(jì)指標(biāo)均優(yōu)于 MSSA 插值的結(jié)果，且隨著連續(xù)缺失比例的增大，該特點(diǎn)越凸顯，水平方向上，RegEM 插值后求得的 Pearson 系數(shù)和Spearman 系數(shù)均大于KKF 插值后求得的相應(yīng)的系數(shù)，垂向上2 者相當(dāng)。

RegEM 和 KKF 插值后的 3 個(gè)統(tǒng)計(jì)指標(biāo)隨連續(xù)缺失比例的變化表現(xiàn)得較平穩(wěn)，即插值結(jié)果與連續(xù)缺失比例基本無(wú)關(guān)，其中各方向插值結(jié)果與原殘差序列的 Pearson 系數(shù)及 Spearman 系數(shù)均在 0.63～0.93 之間，表現(xiàn)出強(qiáng)相關(guān)性。而 MSSA方法的 3 個(gè)統(tǒng)計(jì)指標(biāo)隨數(shù)據(jù)缺失比例變化呈現(xiàn)出較大波動(dòng)，總體表現(xiàn)為缺失比例越大，插值效果越差。

圖6 不同比例連續(xù)缺失數(shù)據(jù)的插值效果比較

在連續(xù)缺失比例較小時(shí)，RegEM 和KKF 插值后的 3 個(gè)統(tǒng)計(jì)指標(biāo)均表現(xiàn)出一定的波動(dòng)性，其原因可能是短期內(nèi)站點(diǎn)噪聲隨機(jī)性大，使插值結(jié)果偏差大，導(dǎo)致統(tǒng)計(jì)結(jié)果呈現(xiàn)較大偏差，而且本文所采用的最小二乘求殘差坐標(biāo)序列的前提，是將噪聲當(dāng)做高斯白噪聲處理，但 GNSS 坐標(biāo)序列的噪聲并非高斯白噪聲，這也導(dǎo)致結(jié)果有偏差，尤其對(duì)于短時(shí)段序列更是如此，此外，GNSS 站點(diǎn)自身的噪聲特性和原始坐標(biāo)序列中的未探測(cè)到的階躍也將對(duì)插值結(jié)果產(chǎn)生影響。

整體上，RegEM 插值效果在垂直方向上與KKF 相當(dāng)，在水平方向上優(yōu)于KKF 方法，而MSSA插值效果較差。

2.2 陸態(tài)網(wǎng)實(shí)測(cè)數(shù)據(jù)插值實(shí)驗(yàn)

本實(shí)驗(yàn)利用RegEM、KKF 和MSSA 3 種插值方法，對(duì)中國(guó)陸態(tài)網(wǎng)華東地區(qū) 10 個(gè)站點(diǎn)的 GNSS坐標(biāo)序列進(jìn)行插值，其中JSLY 和JSYC 站缺失比例較大，且存在長(zhǎng)時(shí)間缺失，連續(xù)缺失數(shù)據(jù)個(gè)數(shù)最大值分別為 120 和 184，其他站點(diǎn)缺失比例較小，對(duì)插值后的連續(xù)GNSS 坐標(biāo)序列進(jìn)行主成分分析，各坐標(biāo)方向前j（j=1,2,3）個(gè)主成分累積貢獻(xiàn)率如表1 所示。

表1 不同插值方法主成分累積貢獻(xiàn)率 %

從表 1 中可以看出，N、E、U3 個(gè)方向前 1～3 個(gè)主成分累積貢獻(xiàn)率從大到小依次是 RegEM、KKF和 MSSA，即 3 種方法中，RegEM 插值結(jié)果所保留的方差最大，插值效果最好，其次為KKF 插值，MSSA 插值效果最差。

為分析插值前后坐標(biāo)序列的噪聲特性，本文利用GNSS 坐標(biāo)序列分析軟件est_noise 對(duì)上述2 個(gè)插補(bǔ)實(shí)驗(yàn)結(jié)果進(jìn)行分析[15-16]，采用的噪聲模型有白噪聲（white noise, WN）+閃爍噪聲（flicker noise,FN）、白噪聲+隨機(jī)游走噪聲（random walk noise,RWN）、白噪聲+閃爍噪聲+隨機(jī)游走噪聲、白噪聲+冪律噪聲（power law noise, PL）和白噪聲+帶通冪律噪聲（band-pass filtered noise+ power law noise,BPPL）等。

通過(guò) RegEM、KKF 和 MSSA 方法插值后的NMDW 站連續(xù)段的GNSS 坐標(biāo)序列（完全通過(guò)插補(bǔ)產(chǎn)生的序列），其MLE 值均比原始序列的大，其中，MSSA 的最大，其次為KKF；而RegEM 插值后的 MLE 值，更接近于原始序列的 MLE 值，所有噪聲模型均表現(xiàn)出此規(guī)律。在白噪聲（white noise, WN）+閃爍噪聲（flicker noise, FN）模型（WN+FN）中，3 種插值方法得到的速度不確定性均比原始序列的小，RegEM 插值后求得的速度不確定性更接近于原始序列的速度不確定性。在所有噪聲模型中，MSSA 插值結(jié)果不僅使得速度不確定性變小，而且對(duì)白噪聲幾乎無(wú)貢獻(xiàn)（WN 模型除外），對(duì)有色噪聲貢獻(xiàn)也較小。

本文同樣比較分析了華東10 個(gè)站點(diǎn)中，含有隨機(jī)小間隔缺失的 AHAQ 站（各方向缺失率約5 %）、連續(xù)缺失的JSLY 站和JSYC 站的GNSS 坐標(biāo)序列插值前后的噪聲特性，基于MLE 方法，利用 est_noise 估計(jì)了 3 個(gè)站點(diǎn) GNSS 坐標(biāo)序列在各噪聲模型下的相關(guān)參數(shù)，求得其最優(yōu)噪聲模型。各插值方法下求得的最優(yōu)噪聲模型結(jié)果如表2 所示，表 3 為最優(yōu)噪聲模型下求得的速度及不確定性，表 4 為最優(yōu)噪聲模型下的 MLE 值。從表 2～表 4 可以看出，各插值方法下求得的最優(yōu)噪聲模型幾乎一致，1 個(gè)可能原因是本文選取的站點(diǎn)連續(xù)缺失的個(gè)數(shù)相對(duì)于整個(gè)序列長(zhǎng)度仍然較少，其缺失率低于10 %。一般而言，由于插值后的序列相比于原始序列增加了新數(shù)據(jù)，故插值后的MLE 值均變小，但RegEM 插值后的MLE 值最小，MSSA的最大。當(dāng)連續(xù)缺失比例較小時(shí)，3 種插值方法求得的整個(gè) GNSS 坐標(biāo)序列，在各噪聲模型下的噪聲幅度和速度不確定性無(wú)明顯差異，隨著連續(xù)缺失比例增大，但隨著連續(xù)缺失比例增大，RegEM 插值結(jié)果更優(yōu)。

表2 各種插值方法下的GNSS 坐標(biāo)序列最優(yōu)噪聲模型

表3 各種插值方法下求得最優(yōu)噪聲模型的速度 單位：mm/a

表4 各種插值方法下求得最優(yōu)噪聲模型的MLE 值

3 結(jié)束語(yǔ)

本文利用中國(guó)陸態(tài)網(wǎng)NMDW 站及周邊9 個(gè)站點(diǎn)、華東10 個(gè)站點(diǎn)的GNSS 坐標(biāo)序列進(jìn)行了插值實(shí)驗(yàn)，比較分析了顧及空間相關(guān)性的RegEM、KKF和常用的 MSSA 方法的插值效果，從3 種方法插補(bǔ)連續(xù)缺失數(shù)據(jù)的結(jié)果中可得出以下結(jié)論：

1）整體而言，RegEM 和 KKF 插值方法優(yōu)于MSSA 插值方法，而RegEM 插值效果最優(yōu)。MSSA插值結(jié)果嚴(yán)重依賴于窗口長(zhǎng)度與插值階數(shù)。RegEM方法的缺陷是，當(dāng)所有站點(diǎn)觀測(cè)數(shù)據(jù)在同一歷元全部缺失時(shí)將無(wú)法進(jìn)行插值。

2）RegEM、KKF 插值效果較為穩(wěn)定，與連續(xù)缺失比例關(guān)系不大，但其插值精度依賴于周圍站點(diǎn)的密集程度、各站點(diǎn)之間的相關(guān)性以及數(shù)據(jù)質(zhì)量。當(dāng)參與插值的站點(diǎn)達(dá)到一定數(shù)量時(shí)，插值效果趨于穩(wěn)定。

3）MSSA 插值后的數(shù)據(jù)較為光滑，使得完全通過(guò)插值產(chǎn)生的序列對(duì)白噪聲幾乎無(wú)貢獻(xiàn)，而僅對(duì)有色噪聲有貢獻(xiàn)且極??；KKF 插值對(duì)有色噪聲有貢獻(xiàn)而對(duì)白噪聲貢獻(xiàn)??；RegEM 插值對(duì)2 者均有貢獻(xiàn)，并與原始序列的較為接近。對(duì)于各插值方法插補(bǔ)后的完整序列，當(dāng)缺失率在一定范圍內(nèi)時(shí)，其最優(yōu)噪聲模型及相應(yīng)的速度及不確定性與原始序列的一致。本文建議在對(duì) GNSS 坐標(biāo)序列進(jìn)行插值時(shí)，應(yīng)盡可能采用RegEM 方法，尤其是連續(xù)缺失較多的情形。