顾及空间相关性的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 方法，尤其是连续缺失较多的情形。