在GNSS坐标数据处理中，数据缺失较为常见，不同的缺失情况对GNSS坐标时间序列的分析会产生不同的影响^[1]。拉格朗日方法适用于连续缺失数据在3个以内的插值^[2], 但随着测量精度的提升，可能会出现插值效果不如简单线性插值的情况^[3]。武艳强等^[4]提出的三次样条法能够在一定程度上解决数据缺失较多的问题，但在数据连续缺失过多时，插值精度会大幅下降；邱荣海等^[5]使用奇异谱迭代的插值方法提高了插值效率，但处理较为复杂，需要对嵌入窗口的插值阶数进行大量验证；尹玲等^[6]使用神经网络进行插值，该方法在观测值缺失较多的情况下依然能够得到较好的插值数据，但操作较为繁琐；苏利娜等^[7]提出的基于模型和噪声间相关性的插值方法对长空缺数据有良好的插值效果，并发现数据空缺对周期性的影响较大，但该方法需要对时间序列具备大量的先验知识；蔡晓军等^[8]提出的多通道奇异谱插值方法不需要坐标时间序列先验信息也能得到较高精度的插值数据，但交叉验证模型的参数且缺失数据需要用0填充可能带来一定偏差。

本文针对GNSS坐标时间序列中的插值问题，提出使用Prophet模型对GNSS坐标时间序列进行插值。该方法不需要关于时间序列的大量先验知识，并能在各种缺失情况下快速得到高精度的插值数据。

1 模型原理与数据处理 1.1 Prophet模型

Prophet模型插值的本质是对缺失数据进行预测，即在拟合的过程中自动填补缺失值从而达到插值的效果^[9]。Prophet模型克服了传统插值方法的缺陷与限制，其在处理缺失值与异常值时鲁棒性极强^[10]。当坐标时间序列中存在缺失值时，Prophet模型能够较好地确定GNSS坐标时间序列的组成成分，并同时对具有多个季节性的周期数据进行模拟^[11], 使得到的插值数据能够较好地契合时间序列的变化规律。

Prophet模型基于贝叶斯方法的曲线拟合来实现预测^[11], 其可用广义加法模型表示：

$y(t)=g(t)+s(t)+h(t)+\varepsilon(t)$

(1)

式中，g(t)为趋势项，表示时间序列中非周期部分的变化趋势；s(t)为周期项或季节项，通常以周或年为单位；h(t)为假日项，表示假日对1 d或多天内不规则变化的影响；ε(t)为残差项，表示模型未预测到的趋势。

g(t)由逻辑回归函数(logistic function)和分段线性函数(piecewise linear function)组成。由逻辑回归函数表示g(t)为：

$g(t)=C /\left(1+\mathrm{e}^{-k(t-m)}\right)$

(2)

式中，C为曲线的最大渐近值，即g(t)随着t的增加趋于C; k为曲线的增长率；m为曲线的偏移量。

s(t)使用傅里叶级数对其进行构造，表达形式为：

$s(t)=\sum\limits_{n=1}^{N}\left(a_{n} \cos \left(\frac{2 \pi n t}{T}\right)+b_{n} \sin \left(\frac{2 \pi n t}{T}\right)\right)$

(3)

式中，T为期望时间序列具有的规则周期，当T=365.25、N=10时，表示以年为周期；当T=7、N=3时，表示以周为周期。

h(t)可用来反映时间序列中某时刻的特殊变动，Prophet模型根据每个假日项在不同时刻下产生的影响构建独立的模型，并为各个假日项设置不同的前后窗口期，以及产生相应的虚拟变量^[10]。h(t)的表达形式为：

$h(t)=\sum\limits_{i=1}^{L} K_{i} 1\left(t \in D_{i}\right) \\ Z(t)=\left[1\left(t \in D_{1}\right), \cdots, 1\left(t \in D_{L}\right)\right] \\ h(t)=Z(t)_{k}, k \sim \operatorname{normal}(0, \gamma)$

(4)

式中，i为节假日，K_i为节假日的影响范围；D_i为i对应的虚拟变量，表示第i个节假日的前后一段时间；L为时间序列中含有的节假日个数。

IGS基准站的坐标时间序列包含着测站的长期线性变化趋势，以及测站受地球物理效应等外界因素影响造成的非线性变化，使测站位置产生周期性的震荡变化^[12-13]。所以当坐标时间序列中存在缺失值时，传统的插值方法无法兼顾时间序列的周期变化。而Prophet模型能够对数据中大量潜在的突变点进行识别与监测，再对趋势变化的幅度进行稀疏先验(与L1正则化效果相同), 并且Prophet模型采用的傅里叶级数可构造适应周期性变化的模型^[14], 这使得Prophet模型在拟合时间序列的同时，能够得到适应时间序列周期变化的插值数据。图 1为Prophet模型对BJFS站坐标数据的拟合图，散点表示原始数据，曲线为拟合数据，阴影区域表示置信区间。由图 1可以看出，Prophet模型能够较好地表现数据的变化趋势以及坐标时间序列的周期性，且没有过度拟合原始数据，使得模型对粗差与离散点有较好的抵抗能力，这为得到高精度的插值数据提供了可能。

1.2 模拟数据的生成

几乎所有的GNSS坐标时间序列都表现出一个季节性的位移周期，它可以被模拟成一个周期项为1 a或0.5 a的4项傅里叶级数序列。最常见的轨迹模型是使用恒定的速度趋势，在这种情况下，模拟数据遵循式(5):

$\begin{aligned} \boldsymbol{x}(t)=& \boldsymbol{x}_{R}+\boldsymbol{v}\left(t-t_{R}\right)+\sum\limits_{j=1}^{n_{J}} \boldsymbol{b}_{j} H\left(t-t_{j}\right)+\\ & \sum\limits_{k=1}^{n_{F}}\left[\boldsymbol{s}_{k} \sin \left(\omega_{k} t\right)+\boldsymbol{c}_{k} \cos \left(\omega_{k} t\right)\right] \end{aligned}$

(5)

式中，x_R为参考位置；t_R为任意的参考时间，通常设置为平均观测时间；v为测站的速度向量，可假定为常数；H为单位步长函数；向量b_j表示在时间t_j上产生跳跃的方向和幅度，n_J为跳跃的次数；s_k与c_k为角频率ω_k的谐波傅里叶系数(每个位置向量对应一个); n_F为不同频率的数量；角频率ω_k=2πτ_k, τ_k为角频率对应的周期。为了模拟年位移周期，选择基本周期τ₁=1 a, 高次谐波周期τ_k=1/k a, 这样可确保由n_F个正弦和n_F个余弦构成的循环每年仅重复一次。根据式(5), 分别模拟东、北、高程3个方向上的坐标时间序列。通常GNSS坐标时间序列的噪声组成不能由单一的噪声模型描述^[15], 故模拟数据中包含线性趋势+季节性信号，噪声选用白噪声与闪烁噪声(WN+FN)的组合，且垂直分量大于水平分量。

2 Prophet模型插值实验

由于实测的GNSS坐标数据包含各种误差以及复杂的噪声模型组合，需要通过模拟数据来检验Prophet模型的插值能力。针对短周期和长周期的坐标时间序列，以及随机缺失和连续缺失等情形，设计3组模拟实验。模拟实验1、2分别是在短周期和长周期的坐标时间序列中对随机剔除的数据进行插值，模拟实验3则对不同长度的连续缺失数据进行插值。在实测数据实验中，分别进行数据随机缺失和连续缺失的插值实验。各实验均使用均方根误差(RMSE)与平均绝对误差(MAE)作为插值精度的评价指标。

2.1 模拟实验1

模拟实验1主要针对短周期坐标时间序列中随机缺失的数据插值。根据式(5)得到模拟测站S1东、北、高程3个方向连续4 a共1 461 d的模拟数据。图 2为S1站3个方向的坐标时间序列图，将各方向时间序列随机剔除5%、10%、15%的数据，再分别用3种方法对缺失数据进行插值，并将得到的插值数据分别与缺失的原始数据进行对比。

由表 1中的精度指标可以看出，3种方法在东向与北向的坐标时间序列插值中都能保持较高的精度；而在高程方向上，尽管各方法的精度都有所下降，但仍能保持良好的插值效果。3种方法在3个方向上的插值精度并没有随着缺失比例的增加而出现明显的降低，表明3种方法在短周期的时间序列中都有较好的插值效果。其中Prophet模型的插值精度最高，三次样条法次之，拉格朗日法最差。

2.2 模拟实验2

模拟实验2主要针对中长周期坐标时间序列随机缺失的数据插值。由式(5)生成模拟测站S2东、北、高程3个方向连续10 a共3 653 d的模拟数据。图 3为S2站在3个方向上的时间序列图，将各方向时间序列随机剔除5%、10%、15%的数据，再分别采用拉格朗日法、三次样条法与Prophet模型对缺失数据进行插值，并将3种方法得到的插值数据与缺失的原始数据进行对比。

从表 2中的精度指标可以看出，针对较长周期的坐标时间序列，3种模型都能较好地对缺失数据进行插值。但拉格朗日法的插值精度波动较大，这可能是由于当缺失点较多时，拉格朗日插值多项式的次数增高，导致插值数据不稳定；而Prophet模型与三次样条法依然能够提供高精度的插值数据。与模拟数据实验1的结果相同，Prophet模型在3种方法中插值精度最高。

2.3 模拟实验3

模拟实验3主要检验Prophet模型在数据连续缺失时的插值效果。由式(5)得到模拟测站S3东、北、高程3个方向连续6 a共2 192 d的模拟数据。图 4为S3站在3个方向上的时间序列图，由于拉格朗日法对连续缺失数据的插值效果较差，故实验中不使用该方法。以2003-01-01为起点，分别剔除连续7 d、30 d、60 d、0.5 a、1 a、2 a共6种时间长度的数据，再采用三次样条法与Prophet模型对连续缺失数据进行插值，并将2种方法得到的插值数据与缺失的原始数据进行对比。由于三次样条法在60 d时的插值精度已经较低，故在缺失时间为0.5 a、1 a、2 a时只检验Prophet模型的插值精度。

由表 3的精度指标可知，2种模型在东向与北向的插值精度都要优于高程方向。其中，三次样条法在数据连续缺失7 d和30 d时，能够提供可用的插值数据；当连续缺失数据为60 d时，其插值精度会大幅度下降。而Prophet模型在3种不同的缺失长度下都能够保持极高的插值精度；在缺失时间为0.5 a、1 a、2 a时，Prophet模型的插值精度并没有随着缺失时间的增加而明显降低，仍能得到高精度的插值数据。

2.4 GNSS实测高程坐标时间序列插值实验

为进一步检验Prophet模型插值方法的可靠性，将处理后得到的BJFS站(2008~2014年)高程方向坐标时间序列(图 1)进行插值研究与分析。数据来源于中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn/), 数据采样间隔为1/365.25 a, 采样频率为1/365.25 Hz。该站的时间跨度为6 a共2 192 d。按照5%、10%、15%的比例随机剔除数据，然后对比Prophet模型、拉格朗日法以及三次样条法对实测数据的插值效果。

由表 4可知，3种插值方法在数据随机缺失比例为5%、10%、15%时，对高程方向实测数据的插值都能保持较高的精度，且插值效果十分稳定。其中，Prophet模型的插值精度最高。

当BJFS站高程方向的实测数据连续缺失时，与模拟实验3相同，只选用Prophet模型与三次样条法对缺失数据进行插值。以2010-01-01为起点，分别剔除7 d、30 d、60 d、0.5 a、1 a、2 a共6种时间长度数据，同样在缺失时间为0.5 a、1 a、2 a时只检验Prophet模型的精度(表 5)。

由表 5可知，Prophet模型与三次样条法在实测数据连续缺失7 d时，都能够得到较高的插值精度；而当实测数据连续缺失30 d和60 d时，2种方法的插值精度都有不同程度的下降，但Prophet模型的插值精度远高于三次样条法；在连续缺失时间为0.5 a、1 a、2 a时，Prophet模型的精度只有小幅度下降，且表现得十分稳定。

图 5是Prophet模型与三次样条法对BJFS站高程方向上连续缺失数据的插值对比图。可以看出，三次样条法在连续空缺较多时插值数据会偏离原始数据，从而产生较大的偏差；而Prophet模型得到的插值数据更加符合原始数据的变化。图 6(2010年)是使用Prophet模型在连续缺失0.5 a、1 a、2 a时的插值图。可以看出，在缺失量较大时，Prophet模型的插值效果依然稳定，并能够较好地体现原始数据的趋势变化与周期性。

2.5 Prophet模型的灵活度

使用Prophet模型插值时，时间序列趋势的变化灵活性不足(即拟合不足)或过于灵活(即过度拟合), 导致模型的默认参数(默认灵活度为0.05)有时不能得到最优的插值数据，这时需要手动调整模型的灵活性以调整稀疏先验的程度。

以模拟测站S1与BJFS站高程方向坐标时间序列为例分析灵活度对Prophet模型插值的影响。按照15%的比例随机剔除数据，再选择常用的灵活度参数得到Prophet模型在不同灵活度下的插值数据，最后将每次得到的插值数据与原始数据对比得到相应的RMSE。图 7中，x轴表示选取的灵活度参数，y轴为对应的RMSE。可以看出，通常在灵活度为0.001时模型的插值效果最差，当灵活度上升到0.05时插值效果会有较大提升。插值精度往往会随着灵活度的上升而提高；但达到一定的精度后，即使灵活度上升，精度也不再提高，而是逐渐趋于平稳甚至有小幅下降。

3 结语

1) 在建模过程中，由于Prophet模型不会过度拟合原始数据，所以能够降低粗差或离群值对模型的影响。Prophet模型在GNSS坐标时间序列中具有较好的适用性，这为Prophet模型的高精度插值提供了可能性。

2) Prophet不是单纯的数学模型，它能够兼顾数据的趋势变化，从而提供良好的插值效果。Prophet模型在短周期与中长周期的坐标时间序列中都能保持良好的插值效果。而当数据连续缺失时，Prophet模型较传统模型优势更明显。在小比例和非连续缺失的情况下，仍可选用拉格朗日法与三次样条法进行插值，但Prophet模型通常能够提供更高精度的插值数据。

3) Prophet模型的优势在于简便性与快捷性，其不需要太多关于时间序列的先验知识，尽管它可以根据专业需求来调整模型以达到更好的效果，但在默认设置下仍能给予较好的反馈。所以处理数据时不用过多地关注模型的选择，而可以将重点放在对数据的分析中。

4) 当对数据插值的精度要求较高时，Prophet模型有时不能一次得出理想的结果，这时需要不断调整模型的灵活度，而模型本身又无法自动给出最优解，所以尽管可以根据经验较快地锁定最优灵活度的范围，但仍会增加数据处理的时间。