顧及空間相關性的GNSS 坐標序列插值比較

謝春橋，匡翠林

（中南大學 地球科學與信息物理學院，長沙 410083）

0 引言

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

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

1 RegEM 和 KKF 插值算法

由于GNSS 坐標序列中含有線性項和周期項，同時還存在粗差以及階躍，粗差與階躍會影響插值結果。因此，本文采用經過四分位距（interquartile range, IQR）方法剔除粗差并修復階躍后的殘差坐標序列，再進行插值[12-13]。

1.1 KKF 插值算法

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

在計算過程中，首先對Xm|t向量和Hm矩陣賦初值0，則設Q為觀測噪聲方差－協方差矩陣；R為狀態噪聲的等價觀測方差－協方差陣；Φ為狀態轉移矩陣；P為估計的協方差陣，利用卡爾曼濾波（式（2）～式（6））可求得tα：

1）狀態一步預測為

2）一步預測協方差陣為

3）濾波增益為

4）狀態估計為

5）估計的協方差陣為

再利用 Kriging 模型求當前歷元缺失觀測值的空間場Hm，最后計算當前歷元缺失的觀測值為

式中，通過EM 算法確定Kalman 濾波初值及相關參數與Q等。關于 KKF 插值更詳細的描述可參考文獻[7-9]。

1.2 RegEM 插值算法

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

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

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

求出缺失數據估值的公式為

插補完成后，需對插值后的數據進行更新，重新計算均值和協方差陣，再次進行迭代計算，當估計的缺失值的均值和協方差達到指定的閾值時，便停止迭代。關于RegEM 算法更詳細的推導見參考文獻[10]。

1.3 插值精度評價指標

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

在這些指標中，MAE 為絕對指標，Pearson 系數和Spearman 系數為相對指標。若插值算法效果越好，則MAE 值越小，Pearson 系數和Spearman系數越趨近于1。

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

2 GNSS 坐標序列插值

2.1 模擬插值實驗

本實驗采用中國陸態網數據缺失率較低的NMDW 站及其周邊9 個站點的坐標序列，站點分布如圖 1 所示，其中 NMDW 站數據缺失率為3.06 %，該站有1 段800 多天完全連續數據。為研究不同比例連續缺失下，各插值方法的效果，針對 NMDW 站各方向連續段第 20～800 個數據，以20 個數據為步長，人為地移除數據，從而構成40 組插值實驗數據，再與周邊9 個站點的數據進行插值。

圖1 GNSS 實驗站點分布

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

圖2 NMDW 站GNSS 殘差坐標序列

圖3 原殘差序列與MSSA 插值結果比較

圖6 為NMDW 站實驗中，不同缺失比例下插值后求得的MAE、Pearson 系數及Spearman 系數，其中Pearson 系數及Spearman 系數為相對指標。從圖6 可以得出：

圖4 原殘差序列與KKF 插值結果比較

圖5 原殘差序列與RegEM 插值結果比較

對于不同比例的連續缺失情形，RegEM 和KKF插值后的各統計指標均優于 MSSA 插值的結果，且隨著連續缺失比例的增大，該特點越凸顯，水平方向上，RegEM 插值后求得的 Pearson 系數和Spearman 系數均大于KKF 插值后求得的相應的系數，垂向上2 者相當。

RegEM 和 KKF 插值后的 3 個統計指標隨連續缺失比例的變化表現得較平穩，即插值結果與連續缺失比例基本無關，其中各方向插值結果與原殘差序列的 Pearson 系數及 Spearman 系數均在 0.63～0.93 之間，表現出強相關性。而 MSSA方法的 3 個統計指標隨數據缺失比例變化呈現出較大波動，總體表現為缺失比例越大，插值效果越差。

圖6 不同比例連續缺失數據的插值效果比較

在連續缺失比例較小時，RegEM 和KKF 插值后的 3 個統計指標均表現出一定的波動性，其原因可能是短期內站點噪聲隨機性大，使插值結果偏差大，導致統計結果呈現較大偏差，而且本文所采用的最小二乘求殘差坐標序列的前提，是將噪聲當做高斯白噪聲處理，但 GNSS 坐標序列的噪聲并非高斯白噪聲，這也導致結果有偏差，尤其對于短時段序列更是如此，此外，GNSS 站點自身的噪聲特性和原始坐標序列中的未探測到的階躍也將對插值結果產生影響。

整體上，RegEM 插值效果在垂直方向上與KKF 相當，在水平方向上優于KKF 方法，而MSSA插值效果較差。

2.2 陸態網實測數據插值實驗

本實驗利用RegEM、KKF 和MSSA 3 種插值方法，對中國陸態網華東地區 10 個站點的 GNSS坐標序列進行插值，其中JSLY 和JSYC 站缺失比例較大，且存在長時間缺失，連續缺失數據個數最大值分別為 120 和 184，其他站點缺失比例較小，對插值后的連續GNSS 坐標序列進行主成分分析，各坐標方向前j（j=1,2,3）個主成分累積貢獻率如表1 所示。

表1 不同插值方法主成分累積貢獻率 %

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

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

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

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

表2 各種插值方法下的GNSS 坐標序列最優噪聲模型

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

表4 各種插值方法下求得最優噪聲模型的MLE 值

3 結束語

本文利用中國陸態網NMDW 站及周邊9 個站點、華東10 個站點的GNSS 坐標序列進行了插值實驗，比較分析了顧及空間相關性的RegEM、KKF和常用的 MSSA 方法的插值效果，從3 種方法插補連續缺失數據的結果中可得出以下結論：

1）整體而言，RegEM 和 KKF 插值方法優于MSSA 插值方法，而RegEM 插值效果最優。MSSA插值結果嚴重依賴于窗口長度與插值階數。RegEM方法的缺陷是，當所有站點觀測數據在同一歷元全部缺失時將無法進行插值。

2）RegEM、KKF 插值效果較為穩定，與連續缺失比例關系不大，但其插值精度依賴于周圍站點的密集程度、各站點之間的相關性以及數據質量。當參與插值的站點達到一定數量時，插值效果趨于穩定。

3）MSSA 插值后的數據較為光滑，使得完全通過插值產生的序列對白噪聲幾乎無貢獻，而僅對有色噪聲有貢獻且極小；KKF 插值對有色噪聲有貢獻而對白噪聲貢獻小；RegEM 插值對2 者均有貢獻，并與原始序列的較為接近。對于各插值方法插補后的完整序列，當缺失率在一定范圍內時，其最優噪聲模型及相應的速度及不確定性與原始序列的一致。本文建議在對 GNSS 坐標序列進行插值時，應盡可能采用RegEM 方法，尤其是連續缺失較多的情形。