近20 a来，全球IGS基准站不断积累的时间序列数据为大地测量和地球动力学研究提供了宝贵的数据基础^[1-3]。这些数据可有效反映由地球物理效应引起的长期变化趋势和非线性变化^[4]，因此对GNSS坐标时间序列进行分析有助于地壳板块运动监测^[5]、大坝或桥梁变形监测^[6-7]、全球或区域坐标系维护^[8-9]等领域的发展。通过分析GNSS坐标时间序列，可以预测连续时间点坐标，为判断运动趋势提供重要依据。研究发现，GNSS坐标时间序列中高程方向的噪声分量通常大于水平方向，且噪声组合模型丰富，因此构建高精度的GNSS高程时间序列预测模型具有较大困难^[10-12]。

目前，已有许多学者在GNSS高程坐标时间序列预测研究中基于机器学习、深度学习等算法构建相应的预测模型，并针对GNSS高程坐标时间序列的信号特征进行高精度的预测研究工作^[11]。但是相关研究仍存在一定不足：1)模型预测精度容易受到异常值以及缺失值影响，抗干扰能力较弱，导致使用模型时对数据质量具有较高的要求。2)模型单次预测精度无法保证随着预测时间跨度的增长持续保持在较高水平上，在进行长期预测时需要进行滚动预测以保持模型精度，从而会增大模型的运行载荷。3)在相同参数设置和数据集下，部分模型的预测结果存在差异，并且无法进行动态调整。

针对这些问题，本文提出一种Prophet-XGBoost组合预测模型用于GNSS高程时间序列预测，并通过SOPAC公布的实测高程数据对本文方法的适用性和有效性进行检验。首先检验XGBoost模型和Prophet模型在GNSS时间序列中的适应度；再基于单模型在长期预测中表现出的特点，针对XGBoost模型对于趋势项部分预测能力较弱和Prophet模型对于非线性部分预测能力较弱的问题，构建Prophet-XGBoost组合预测模型；最后通过ALGO、ALRT和BRST基准站U分量坐标时间序列数据对模型预测精度进行检验。

1 模型原理与评价指标 1.1 XGBoost模型

XGBoost作为一种集成学习算法，可通过构建多个CART树对数据集进行预测，然后集成多个树模型的预测结果，从而输出最终预测数据。与同属于树模型算法的GBDT算法相比，XGBoost算法更为高效，并且在算法方面也进行过改进。XGBoost使用二阶泰勒展开式逼近目标函数的泛化误差部分，可有效简化目标函数的计算；XGBoost还可通过在目标函数中加入正则项来降低模型预测的波动性以及改善模型过拟合现象^[13]。

假设存在数据集A={(x_i, y_i): i=1，…，n, x_i∈$\mathscr{R} $^m, y_i∈$\mathscr{R} $}，其中含有n个观测值，每个观测值有m个特征和一个对应的变量y。定义数值${\hat y_i}$并通过广义模型表示为：

$ \hat{y}_i=\varphi\left(x_i\right)=\sum\limits_{k=1}^K f_k\left(x_i\right) $

(1)

式中，f_k表示决策树，f_k(x_i) 表示第k棵树赋予第i个观测值的分数。同时，使用函数f_k时，下述正则化目标函数应被最小化：

$ L(\varphi)=\sum\limits_i l\left(y_i, \hat{y}_i\right)+\sum\limits_k \varOmega\left(f_k\right) $

(2)

式中，l为损失函数。为防止模型过于复杂，在模型中设置惩罚项Ω：

$ \varOmega\left(f_k\right)=\gamma T+\frac{1}{2} \lambda\|w\|^2 $

(3)

式中，γ可控制惩罚项中枝叶数量T，λ可控制惩罚项中枝叶重量w。设置Ω(f_k) 项不仅可以简化算法生成模型，同时可以防止模型过拟合。

XGBoost算法采用迭代法最小化目标函数。模型通过增加f_j项在第j次迭代后得到减小的目标函数：

$ L^j=\sum\limits_{i=1}^n l\left(y_i, \hat{y}_i^{(j-1)}+f_j\left(x_i\right)\right)+\varOmega\left(f_j\right) $

(4)

式(4)可以通过泰勒展开进行简化，并且可以推导出决策树从给定节点分裂后的损失减少公式：

$ \begin{gathered} L_{\text {split }}= \\ \frac{1}{2}\left[\frac{\left(\sum\limits_{i \in I_L} g_i\right)^2}{\sum\limits_{i \in I_L} h_i+\lambda}+\frac{\left(\sum\limits_{i \in I_R} g_i\right)^2}{\sum\limits_{i \in I_R} h_i+\lambda}-\frac{\left(\sum\limits_{i \in I} g_i\right)^2}{\sum\limits_{i \in I} h_i+\lambda}\right]-\gamma \end{gathered} $

(5)

式中，I为当前节点中可用观测数据的一个子集，I_L和I_R分别为分裂后左右节点中可用观测数据的一个子集。函数g_i和h_i可表示为：

$ g_i=\partial_{\hat{y}_i^{(j-1)}} l\left(y_i, \hat{y}_i^{(j-1)}\right) $

(6)

$ h_i=\partial_{\hat{y}_i^{(j-1)}}^2 l\left(y_i, \hat{y}_i^{(j-1)}\right) $

(7)

从推导L_split公式中可找到任意给定节点的最佳分裂，该函数只依赖于损失函数和正则化参数γ。同时，基于式(5)，XGBoost算法可以优化任何损失函数，并且提供一阶和二阶梯度。

1.2 Prophet模型

Prophet模型是一种基于加性模型的时间序列预测模型，由趋势项、周期项和假日项组成^[14-15]，表达式为：

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

(8)

式中，g(t) 为趋势项；s(t) 为周期项；h(t) 为假日项，该项在GNSS时间序列中通常不存在，因此建模过程中不予考虑；ε_t为残差项，即模型未能预测到的随机趋势。

在Prophet模型中，趋势项包含逻辑回归函数和分段线性函数两个重要函数。根据GNSS高程时间序列的趋势变化特点，选用逻辑回归函数来计算趋势项。同时，Prophet模型使用傅里叶级数表示时间序列的周期性：

$ s(t)=\sum\limits_{n=1}^N\left(a_n \cos \left(\frac{2 \pi n t}{P}\right)+b_n \sin \left(\frac{2 \pi n t}{P}\right)\right) $

(9)

式中，P为时间序列周期，P=365.25时为年周期项，P=7时为周周期项。因此，Prophet模型在进行GNSS日坐标时间序列分解时会分别输出年周期项和周周期项。

1.3 精度评价指标

使用平均绝对误差(MAE)和均方根误差(RMSE)作为模型预测精度的评价指标：

$ {\mathrm{MAE}}=\frac{1}{n} \sum\limits_{i=1}^n\left|y_i-\hat{y}_i\right| $

(10)

$ {\mathrm{RMSE}}=\sqrt{\frac{1}{n} \sum\limits_{i=1}^n\left(y_i-\hat{y}_i\right)^2} $

(11)

式中，y_i为原始值，${\hat y_i}$为预测值。MAE和RMSE数值越小表示模型预测精度越高，越适用于该时间序列。反之，则代表模型预测精度越低，在该时间序列中适用性较差。

通过引入Pearson相关系数判断预测时间序列和原始时间序列之间的相关性。Pearson相关系数将时间序列视为变量，从而计算两个时间序列之间的相关性，其计算公式为：

$ \rho_{Y \hat{Y}}=\frac{E(Y \hat{Y})-E(Y) E(\hat{Y})}{\sigma_Y \sigma_{\hat{Y}}} $

(12)

式中，Y和$\hat{Y}$分别表示原始时间序列和预测时间序列。该公式值域为[-1, 1]，$\rho_{Y \hat{Y}}$>0表示两时间序列正相关，$\rho_{Y \hat{Y}}$ < 0表示两时间序列负相关，$\rho_{Y \hat{Y}}$=0表示两时间序列不相关。同时，$\rho_{Y \hat{Y}}$绝对值越接近1，两时间序列相关性越强^[16]。

2 数据与模型构建 2.1 实验数据

为验证Prophet-XGBoost模型的有效性，选取IGS基准站实测数据作为数据集。选用ALGO、ALRT和BRST站2010~2015年U分量日坐标时间序列作为实验数据，该数据来源于SOPAC(scripps orbit and permanent array center)。数据基本情况如表 1(单位mm)所示。

2.2 模型适用性检验

为检验XGBoost模型和Prophet模型对于GNSS坐标时间序列预测的适用性，引入经典时间序列预测模型ARIMA模型作为对比。选用ALGO站数据作为实验数据，分别用ARIMA模型、XGBoost模型和Prophet模型对数据集进行预测，将预测值与原始数据值进行对比，得到相应的MAE值和RMSE值。

该实验数据集中2010~2014年数据为训练集，2015年数据为验证集，共包含2 149个日坐标数据，缺失数据天数为42，数据完整率为98.08%，训练中最大连续缺失数据时间跨度为7 d。

由表 1和四分位距值(interquartile range，IQR)计算可知，该实验数据集中含有少量异常值。因此，该数据训练集中包含异常值、缺失值、连续缺失值等条件，可以更好地检验模型对于GNSS实测数据预测的适用性。各模型精度见表 2(单位mm)。

由表 2可知，相较于ARIMA模型，XGBoost模型和Prophet模型在MAE和RMSE两个精度评价指标上均有优势。基于3个模型的预测原理可知，在长时间跨度的GNSS坐标时间序列预测研究中，ARIMA模型精度更容易受到预测时间长度的影响，而XGBoost模型和Prophet模型的预测精度随预测跨度变化而变化的幅度较小，且XGBoost模型的算法实现更为简便，具有更高的预测效率。因此，XGBoost模型和Prophet模型均可应用于GNSS坐标时间序列预测，并在实测数据实验中具有适用性。图 1为ALGO站Prophet模型预测结果，其中离散点为原始时间序列值，曲线为拟合预测值。

2.3 Prophet-XGBoost模型构建

基于XGBoost模型和Prophet模型在GNSS坐标时间序列预测中表现出的模型特性，构建Prophet-XGBoost组合预测模型，以更好地发挥两个模型的预测优势，削弱模型对于GNSS时间序列部分分解项预测能力不足的影响，进而实现更高精度的GNSS坐标时间序列预测。图 2为Prophet-XGBoost模型构建流程。

Prophet-XGBoost模型实现流程如下：

1) 通过Prophet模型对GNSS坐标时间序列U方向数据进行分解，将时间序列分解为趋势项、周期项和残差项，其中周期项由年周期项和周周期项等权相加构成。图 3为ALGO站Prophet模型分解结果。

2) 通过Prophet模型对分解出的趋势项和周期项进行拟合，然后使用XGBoost模型对拟合时间序列进行预测，得到预测值y_fit。图 4为ALGO站XGBoost模型拟合时间序列预测结果。

3) 使用XGBoost模型对Prophet模型分解时间序列过程中遗留的残差项进行预测，得到预测值y_residual。

4) 将两个预测结果等权相加，得到预测值y_pred：

$ y_{\text {pred }}=y_{\text {fit }}+y_{\text {residual }} $

(13)

2.4 实验结果分析

本文选用ALGO、ALRT、BRST三个IGS站U分量GNSS坐标时间序列作为实验数据，分别采用XGBoost模型、Prophet模型和Prophet-XGBoost模型进行预测。其中，ALGO和ALRT站纬度差异较大，ALGO和BRST站经度差异较大。预测精度对比结果如表 3所示。

由表 3可知，相较于XGBoost模型和Prophet模型，ALGO站Prophet-XGBoost模型的MAE值分别降低7.4%和11.5%，RMSE值分别降低16.6%和3.5%；ALRT站Prophet-XGBoost模型的MAE值分别降低7.1%和20.5%，RMSE值分别降低7.2%和6.5%；BRST站Prophet-XGBoost模型的MAE值分别降低10.5%和9.2%，RMSE值分别降低5.6%和3.5%。相较于Prophet模型，组合模型的RMSE值降低幅度较小，这是因为Prophet-XGBoost模型进行预测时的拟合数据是由Prophet模型拟合所得。对Prophet-XGBoost模型预测结果和原始数据值进行相关性检验，通过引入Pearson相关系数计算ALGO、ALRT和BRST三个站的ρ值分别为0.53、0.64和0.49，其中ALGO和BRST两个站预测结果和原始数据呈中等相关，ALRT站预测结果和原始数据呈强相关。

以ALGO站为例，对比分析XGBoost模型、Prophet模型和Prophet-XGBoost模型预测结果误差分布。图 5为各模型预测结果绝对误差分布情况。

由图 5可知，Prophet-XGBoost组合模型的较高精度预测结果占比最高，即相较于Prophet和XGBoost单模型预测结果，组合模型准确率得到提升。图 6为模型预测结果误差对比。

由图 6可知，相较于XGBoost模型，Prophet-XGBoost组合模型的误差变化幅度和凸点值更小；Prophet模型和Prophet-XGBoost组合模型的误差变化趋势较为一致，但Prophet-XGBoost组合模型误差更小。结合图 4结果可知，Prophet模型对于时间序列的趋势变化较为敏感。

3 结语

本文基于XGBoost模型和Prophet模型提出一种GNSS高程时间序列预测组合模型Prophet-XGBoost。在数据处理阶段，利用Prophet模型对时间序列进行分解，将时间序列分解成3个主要组成部分，有利于XGBoost模型在预测时提取时间序列中不同组成部分的特征，减弱预测模型提取特征时受到其他组成部分的影响。得益于此，XGBoost模型可以进行更高精度的GNSS坐标时间序列预测工作，也可弥补两个模型在单模型预测时存在的问题。此外，在计算预测结果时，也可依据XGBoost模型提供的特征分数进行定权相加。

XGBoost模型和Prophet模型在实验过程中也存在需要改进之处：XGBoost模型对于GNSS坐标时间序列的趋势变化不够敏感以及Prophet模型对于GNSS坐标时间序列非线性部分预测能力较弱。同时考虑到XGBoost模型的特征项构建具有很大灵活性，研究人员可以使用非GNSS坐标时间序列数据辅助进行特征构造，如水文载荷对GNSS坐标时间序列U方向振幅具有较为显著的影响^[17]。因此，在以后研究中可以将水文载荷等因素纳入特征构造条件中，从而通过更多的特征子集进行模型构建和预测。基于该特性，Prophet-XGBoost模型可以作为通用预测模型，将其推广到其他研究领域，如海平面时间序列预测、天气预报和发病率预测，可极大增强该模型的适用性。