在全球卫星导航系统中，星载原子钟的钟差预报在优化导航电文中的钟差参数、满足实时动态精密单点定位的需求和提供卫星自主导航所需的先验信息方面具有重要作用^[1-4]。但是，由于星载原子钟非常敏感，极易受到外界和自身因素变化的影响而很难掌握其复杂细致的变化规律。因此，建立精确的原子钟运行模型变得非常困难，相应的准确预报卫星钟差也非常困难^[5-7]。

目前，国内外许多学者在钟差预报方面开展了广泛深入的研究，取得了一系列成果。其中，应用较多的卫星钟差预报模型主要有：二次多项式模型(quadratic polynomial model，QPM)、灰色预报模型(gray model，GM(1, 1))、最小二乘支持向量机模型、人工神经网络模型、卡尔曼滤波模型等^[8-13]。这些预报方法分别适用于不同条件下导航卫星原子钟钟差的短期、中长期和长期的预报，但也均有各自的适用范围和局限性。

本文考虑卫星钟差的随机项呈现非平稳变化趋势，提出了GM(1, 1)与ARIMA的组合预报模型。该组合模型首先采用GM(1, 1)模型预报钟差的趋势项，然后利用ARIMA模型对GM(1, 1)模型残差序列进行建模和预报，最后将GM(1, 1)和ARIMA模型的预报结果对应相加即得到钟差的最终预报值。此外，采用IGS公布的事后精密钟差产品进行预报试验，通过与卫星钟差预报中常用的二次多项式模型和修正指数曲线法模型(modified exponential curve method，MECM)预报结果的对比分析，结果表明：该方法可以对GPS卫星钟差进行高精度的中短期预报，验证了该方法的优越性和有效性。

1 GM(1, 1)预报模型

在灰色系统理论中，GM(1, 1)预报模型是最常用的一种灰色系统模型，它是由一个仅包含单变量的一阶微分方程所构成的预报模型，适合对自身数据的预报^[14]。

设序列为，通过一次累加生成新的序列为

对累加序列X⁽¹⁾建立一阶微分方程

式中，参数a为发展系数；参数u为灰作用量。对式(2)离散化，可得矩阵方程

根据最小二乘法可得式(3)的最小二乘解为

将式(4)代入式(2)即可得预报模型为

式中，k为参与预报序列的个数。由以上模型即可求得未来任意时刻序列的预报值。

2 ARIMA预报模型

ARIMA(p, d, q)模型常用来做时间序列的预报，适用于平稳时间序列。对于非平稳时间序列需要进行d次差分将非平稳时间序列转化为平稳时间序列，然后对平稳时间序列进行模型定阶和模型参数的估计，最后选择合适的模型阶数建立ARIMA模型进行预报。ARIMA(p, d, q)模型的实质是差分运算与ARMA(p, q)模型的组合^[15]。

设有一非平稳时间序列z_t，经过d次差分运算已转化为平稳时间序列x_t，然后对序列x_t建立ARMA(p, q)模型，其定义为

式中，x_t为样本；φ_i(i=1, 2, …, p)为自回归系数；θ_i(i=1, 2, …, q)为移动平均系数；ε_t为零均值，方差为σ_pq²的平稳白噪声序列；p为自回归阶数；q为移动平均阶数。当p=0时，式(6)为MA(q)模型；当q=0时，式(6)为AR(p)模型。

ARIMA模型建立与预报的步骤可归结为如下过程：

(1) 根据时间序列散点图、自相关函数、偏相关函数和ADF单位根检验其方差、趋势变化规律，从而对序列的平稳性进行识别。

(2) 若是非平稳时间序列，需对时间序列进行d次差分处理转化为平稳序列。

(3) 利用平稳时间序列的自相关函数和偏相关函数图确定模型阶数。

(4) 进行模型参数的估计和假设检验。

(5) 选择合适的参数建立ARIMA模型并进行外推预报。

3 组合预报过程

卫星钟差数据序列可以认为是由趋势项和随机项两部分叠加而成，针对它的这一特点，提出了一种基于GM(1, 1)与ARIMA的组合预报模型。运用此组合模型去预报卫星钟差的预报流程如图 1所示，具体步骤如下所述：

(1) 卫星钟差观测数据预处理。首先剔除卫星钟差观测序列中的异常钟差数据，然后采用拉格朗日插值法将缺失的钟差数据补齐。

(2) 利用GM(1, 1)预报模型提取卫星钟差序列的趋势项，并且对其进行外推预报。

(3) 利用ARIMA预报模型对GM(1, 1)模型的模型残差建模，即对卫星钟差的随机项部分进行建模并外推预报。

(4) 将利用GM(1, 1)模型预报的趋势项序列和利用ARIMA模型预报的随机项序列对应相加，即得到卫星钟差序列最终的预报值。

4 算例与分析 4.1 数据来源

为了验证该组合模型的有效性，从IGS服务器(ftp://cddis.gsfc.nasa.gov/)上下载了2017年3月5日至7日共3天的IGS事后精密钟差产品，其采样间隔为15 min。考虑目前在轨的GPS卫星有30多颗且我国北斗二代系统均搭载的是铷原子钟，因此随机选3颗装载铷原子钟的卫星进行预报试验，它们分别为PRN17(Block ⅡF-M-Rb)、PRN03(Block ⅡF-Rb)和PRN16(Block ⅡR-Rb)号卫星。

4.2 建模方案与结果分析

算例中采用12 h钟差数据(2017年3月5日00:00至12:00)共49个历元，分别建立二次多项式模型(QPM)、修正指数曲线法(MECM)预报模型、GM(1, 1)与ARIMA的组合模型(GM(1, 1)+ARIMA)，对未来6、12、24和48 h的卫星钟差进行中短期预报。

由于本试验使用的是IGS服务器上公布的事后精密钟差产品，其自身的误差小于0.1 ns，故可以作为“真值”，使用均方根误差(RMS)作为统计量，用于检验二次多项式模型、修正指数曲线法预报模型、GM(1, 1)与ARIMA的组合模型所预报结果的准确程度。

将接下来6、12、24和48 h的实际观测钟差数据与各模型预报的钟差数据对应相减即得到预报误差，取各模型预报误差序列的均方根误差(RMS)为衡量预报精度的标准，结果如图 2—图 5和表 1—表 4所示。其中图 2、图 3和图 4分别为用12 h钟差数据建立二次多项式模型、修正指数曲线法预报模型、GM(1, 1)与ARIMA的组合模型预报未来6、12、24和48 h钟差数据的预报误差变化图，图 5为各模型预报钟差的平均预报精度对比图。表 1、表 2和表 3分别为用12 h钟差数据建立二次多项式模型、修正指数曲线法预报模型、GM(1, 1)与ARIMA的组合模型预报未来6、12、24和48 h钟差数据预报误差的最大值(Max)、平均值(Mean)和均方根误差(RMS)统计表。

结合图 2—图 5和分析表 1—表 4的统计量可以看出：在所选择的预报时间长度内，GM(1, 1)与ARIMA组合模型对于PRN03卫星钟差预报的精度最高，其次是PRN17卫星，再次是PRN16卫星，且该组合模型预报的均方根误差均小于二次多项式模型和修正指数曲线法模型预报对应的均方根误差。这说明对于卫星钟差预报，该组合模型较常用的二次多项式模型和修正指数曲线法模型有更高的预报精度和稳定性。在短期预报中，6 h预报的平均预报精度控制在1 ns以内，相比于常用的二次多项式模型，平均预报精度提高了29.70%，相比于修正指数曲线法，平均预报精度提高了18.39%；12 h预报的平均预报精度控制在2 ns以内，相比于常用的二次多项式模型，平均预报精度提高了43.75%，相比于修正指数曲线法，平均预报精度提高了33.90%；24 h预报的平均预报精度控制在2 ns以内，相比于常用的二次多项式模型，平均预报精度提高了67.62%，相比于修正指数曲线法，平均预报精度提高了61.40%。在中长期预报中，48 h预报的平均预报精度控制在5 ns以内，相比于常用的二次多项式模型，平均预报精度提高了76.21%，相比于修正指数曲线法，平均预报精度提高了70.49%。

5 结语

针对卫星钟差呈趋势项和随机项变化的特点，本文将GM(1, 1)模型与ARIMA模型有机地结合起来，建立了一种基于GM(1, 1)与ARIMA的组合模型，并将其应用于卫星钟差的中短期预报。该组合预报模型首先采用GM(1, 1)模型预报卫星钟差序列中的趋势项成分，然后提取出卫星钟差序列中的随机项成分并建立ARIMA模型，最后将两部分预报的结果对应相加即得到卫星钟差的最终预报值。该方法不仅顾及了卫星钟差序列变化的特点，而且充分发挥了GM(1, 1)模型和ARIMA模型各自的优势，具有所需建模数据少、易于实现和计算速度快等优点。从上述预报试验的结果可以看出，该组合模型只需要使用较少的钟差数据建模即可实现高精度的预报，尤其是在中长期的预报中。该方法相比于各单一预报方法，突出了其预报的有效性和优越性，从而验证了该组合方法的有效性和可行性。