基于级别建模的灰色GM(1,1)模型在卫星钟差中长期预报中的应用

Application of the Grey GM(1,1) Based on Stepwise-Ratio Modelling to Medium and Long-Term Prediction of Satellite Clock Bias

摘要: 针对常规灰色GM(1,1)模型无法反映卫星钟差级比动态变化的问题，提出以钟差相邻历元的级比为建模序列，建立能够体现钟差级比序列变化趋势的级比离散灰色GM(1,1) (stepwise-ratio discrete GM(1,1)，SDGM(1,1))模型。首先，在卫星钟差相邻历元之间作级比生成级比序列；其次，以钟差级比序列为建模对象，建立SDGM(1,1))模型对钟差级比进行外推；最后，根据钟差序列与级比序列之间的关系，将级比外推值还原得到卫星钟差预报值。采用国际GNSS服务组织(International GNSS Service, IGS)提供的最终精密GPS卫星钟差数据进行1 d短期和30 d中长期预报试验，并和常规GM(1,1)模型、二次多项式模型的预报结果进行对比。结果表明，SDGM(1,1))模型可以较好地预报出卫星钟差级比的变化趋势，1 d钟差平均预报精度相比于GM(1,1)和二次多项式模型分别提升62.38%和37.96%，平均预报稳定度分别提升26.6%和10.87%；对于30 d钟差预报，与二次多项式模型相比，SDGM(1,1))模型的平均预报精度和预报稳定度分别提升75.18%和76.74%，比GM(1,1)模型的平均预报精度和预报稳定度高2个数量级。

Abstract: Aiming at the problem that the conventional grey GM(1,1) model cannot reflect the dynamic change of satellite clock bias stepwise-ratio, this paper proposes to take the stepwise-ratio of adjacent epochs of clock bias as the modeling sequence, and establishes the stepwise-ratio discrete grey GM(1,1) model, denoted as SDGM(1,1), which can reflect the change trend of clock bias stepwise-ratio sequence. The stepwise-ratio sequence is first generated by making stepwise-ratio calculation between adjacent epochs of satellite clock bias. Secondly, taking the clock bias stepwise-ratio sequence as the modeling object, a SDGM(1,1) model is established to extrapolate the clock bias stepwise-ratio. Finally, according to the relationship between the clock bias and stepwise-ratio sequence, the extrapolated stepwise-ratio values are restored to obtain the satellite clock bias predictions. The final products of precise GPS satellite clock bias provide by the International GNSS Service (IGS) are taken as data basis to carry out satellite clock bias prediction up into one day and 30 days in future. The results show that the SDGM(1,1) model can better predict the change trend of satellite clock bias stepwise-ratio. Compared with the GM(1,1) and quadratic polynomial model, the average prediction accuracy of one-day clock bias is improved by 62.38% and 37.96% respectively, and the average prediction stability is improved by 26.6% and 10.87% respectively. For the 30-day clock bias prediction, in contrast to the quadratic polynomial model, the average prediction accuracy and prediction stability of the SDGM(1,1) model are improved by 75.18% and 76.74% respectively, and the average accuracy and stability are two orders of magnitude higher than those of the conventional GM(1,1) model.

1. 引言

卫星钟是全球卫星导航系统(Global Navigation Satellite System, GNSS)的核心载荷，它直接决定着GNSS定位、导航和授时的准确性，是整个导航卫星的“心脏” [1]。卫星钟差信息可以从GNSS导航电文发播的广播星历中实时获取，但GPS广播钟差和BDS-3广播钟差的精度分别约为1.4 ns和1.2 ns [2]，无法满足厘米级精密单点定位的应用需求。当前，国际GNSS服务组织(International GNSS Service, IGS)提供的最终精密卫星钟差产品的精度优于0.1 ns，可以满足厘米级精密单点定位的应用需求，但最终产品存在13 d的延时 [3]，无法满足用户对实时定位的需求，因此，通常需要预报12~24 h的卫星钟差。另一方面，在卫星自主导航阶段，导航系统面临地面运控站无法上传最新的卫星钟差数据而不得不使用钟差中长期预报值的问题，地面运控站需要将卫星钟差中长期预报值注入卫星作为卫星自主导航的先验信息 [4]，所以卫星中长期预报是卫星自主导航中必须要解决的一个关键问题。

针对卫星钟差预报，学者们建立了多种预报模型，如二次多项式模型 [5]、灰色GM(1,1)模型 [6]、自回归滑动平均(autoregressive moving average, ARMA)模型 [7]、卡尔曼滤波模型 [8] 和神经网络模型 [9] 等。其中，二次多项式模型物理意义明确，钟差短期预报效果较好，但该模型是时间的函数，预报误差随着时间累计而增大；ARMA对线性变化的卫星钟差预报效果较好，但对非线性变化的钟差预报误差较大；神经网络的钟差预报精度较高，但存在预报结果不稳定的现象；卡尔曼滤波模型需要精确确定钟差噪声方差，但噪声方差难以精确确定；灰色GM(1,1)模型所需样本数据少，常用于钟差中长期预报，但其要求钟差数据呈现严格的指数变化，虽然通过钟差数据累加处理可以使累加生产序列呈现单调变化，但不能保证生成序列具有指数变化特性，这限制了GM(1,1)模型的适应性。

为提升灰色GM(1,1)模型的适应性，一些学者通过提高数据光滑性和优化模型灰参数改进GM(1,1)模型，在不同程度上提高了钟差预报精度 [10] [11]。然而，GM(1,1)模型根据灰参数确定的恒定公比的等比数列对钟差进行拟合和外推，无法反映钟差级比序列的动态变化特性，只对纯指数变化的钟差序列拟合具有无偏性，易造成较大的预报误差 [12]。针对GM(1,1)模型无法反映钟差级比动态变化的问题，本文以卫星钟差级比序列为建模对象，建立基于级比建模的指数函数模型来预报卫星钟差。首先，在相邻历元的卫星钟差间作一次级比生成级比序列，然后将级比序列作为建模序列，建立级比离散GM(1,1) (stepwise-ratio discrete GM(1,1), SDGM(1,1))模型对钟差级比进行外推，在此基础上，根据钟差序列与级比序列之间的关系，将级比外推值还原得到卫星钟差预报值。选取IGS提供的最终精密GPS卫星钟差产品进行预报试验，结果表明，SDGM(1,1)模型能很好地反映卫星钟差级比序列的总体变化趋势，钟差1 d短期和30 d中长期的预报精度与预报稳定度明显好于常规灰色GM(1,1)模型和常用的二次多项式模型。

2. GM(1,1)模型拟合序列级比分析

设 $X^{\left(0\right)}=\left\{x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\cdots ,x^{\left(0\right)}\left(n\right)\right\}$ 是一非负钟差序列，通过对 $X^{\left(0\right)}$ 一次累加处理构造生成序列建立离散GM(1,1)模型，再进行差分还原处理获得原始 $X^{\left(0\right)}$ 序列的离散GM(1,1)预报模型：

${\hat{x}}^{\left(0\right)}\left(k\right)=\{\begin{array}{l}{\hat{x}}^{\left(0\right)}\left(1\right)=x^{\left(0\right)}\left(1\right)k=1\\ \left(1-{\text{e}}^{\hat{a}}\right)\left[x^{\left(0\right)}\left(1\right)-\frac{\hat{u}}{\hat{a}}\right]{\text{e}}^{-\hat{a}\left(k-1\right)}k>1\end{array}$ (1)

式中， $\hat{a}$ 为发展系数， $\hat{u}$ 为灰作用量，可以通过最小二乘法求解。

根据式(1)可以得到GM(1,1)模型拟合序列的级比 ${\hat{c}}^{\left(0\right)}$ ：

${\hat{c}}^{\left(0\right)}\left(k\right)=\frac{{\hat{x}}^{\left(0\right)}\left(k+1\right)}{{\hat{x}}^{\left(0\right)}\left(k\right))}={\text{e}}^{-\hat{a}}$， $1\le k\le n-1$ (2)

由于通过最小二乘法求解的发展系数 $\hat{a}$ 为恒定值，所以由式(2)确定的GM(1,1)模型拟合序列的级比也为恒定值 ${\text{e}}^{-\hat{a}}$。级别恒定的GM(1,1)模型仅对纯指数增长序列的拟合具有无偏性，但卫星钟差序列并非纯指数增长序列，级别并不是恒定值，用恒定的级比代替的变化拟合卫星钟差无法体现钟差的动态变化特性，必然存在预报误差，特别是在钟差序列级比变化幅度较大的情况下，容易导致较大的预报误差。

3. 基于级比建模的卫星钟差预报模型

为了更好地体现卫星钟差的动态变化特性，本文以卫星钟差级比序列作为建模对象建立GM(1,1)模型。本节首先给出了基于级比序列的离散GM(1,1)建模原理，然后阐述了级比离散GM(1,1)模型的建模步骤。

3.1. 级比离散GM(1,1)模型

1) 级比一次生成序列

在相邻历元的卫星钟差之间作级比构造一次生成序列 $C^{\left(0\right)}=\left\{c^{\left(0\right)}\left(1\right),c^{\left(0\right)}\left(2\right),\cdots ,c^{\left(0\right)}\left(n-1\right)\right\}$，其中，级比 $c^{\left(0\right)}\left(k\right)$ 可表示为

$c^{\left(0\right)}\left(k\right)=x^{\left(0\right)}\left(k+1\right)/c^{\left(0\right)}\left(k\right)$， $1\le k\le n-1$ (3)

2) 累加二次生成序列

对级比一次生成序列做累加处理构造二次生成序列 $1\le k\le n-1$，其中， $c^{\left(1\right)}\left(k\right)$ 可表示为

$c^{\left(1\right)}\left(k\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}{c}^{\left(0\right)}\left(k\right)}$ (4)

3) 模型建立

以级比累加二次生成序列 $C^{\left(1\right)}$ 为灰信息建立离散GM(1,1)模型，其数学表达式为

$c^{\left(1\right)}\left(k+1\right)={\beta }_1{c}^{\left(1\right)}\left(k\right)+{\beta }_2$ (5)

式中， ${\beta }_1$ 和 ${\beta }_2$ 为待估参数，可以通过最小二乘法求解：

${\left[{\hat{\beta }}_1{\hat{\beta }}_2\right]}^{\text{T}}={\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}Y$ (6)

其中， $B={\left[\begin{array}{cccc}{c}^{\left(1\right)}\left(1\right)& c^{\left(1\right)}\left(2\right)& \cdots & c^{\left(1\right)}\left(n-2\right)\\ 1& 1& \cdots & 1\end{array}\right]}^{\text{T}}$， $Y={\left[\begin{array}{cccc}{c}^{\left(1\right)}\left(2\right)& c^{\left(1\right)}\left(3\right)& \cdots & c^{\left(1\right)}\left(n-1\right)\end{array}\right]}^{\text{T}}$。

根据初始条件 ${\hat{c}}^{\left(1\right)}\left(1\right)=c^{\left(1\right)}\left(1\right)=c^{\left(0\right)}\left(1\right)$，可以得到级比累加二次生成序列的离散GM(1,1)预报模型：

${\hat{c}}^{\left(1\right)}\left(k\right)=\{\begin{array}{l}{c}^{\left(0\right)}\left(1\right)k=1\\ {\hat{\beta }}_1^{k-1}{c}^{\left(0\right)}\left(1\right)+\frac{1-{\hat{\beta }}_1^{k-1}}{{\hat{\beta }}_1}{\hat{\beta }}_{2}k>1\end{array}$ (7)

4) 模型还原

首先，通过一次累减还原运算，得到级比序列 $C^{\left(0\right)}$ 的外推公式：

${\hat{c}}^{\left(0\right)}(k)=\{\begin{array}{l}{\hat{c}}^{\left(1\right)}\left(1\right)k=1\\ {\hat{c}}^{\left(1\right)}(k)-{\hat{c}}^{\left(1\right)}\left(k-1\right)=\left({\hat{\beta }}_1-1\right){\hat{\beta }}_1^{k-1}{c}^0\left(1\right)+{\hat{\beta }}_1^{k-1}{\hat{\beta }}_{2}k>1\end{array}$ (8)

其一般形式为

${\hat{c}}^{\left(0\right)}\left(n+p\right)=\left({\hat{\beta }}_1-1\right){\hat{\beta }}_1^{n+p-1}c\left(1\right)+{\hat{\beta }}_1^{n+p-1}{\hat{\beta }}_2$ (9)

式中，p表示预报点数，故由式(9)可以对未来任意历元的级比进行外推。将式(9)代入式(3)，获得基于级比序列建模的卫星钟差递推预报公式：

${\hat{x}}^{\left(0\right)}\left(k+1\right)={\hat{c}}^{\left(0\right)}\left(k\right){\hat{x}}^{\left(0\right)}\left(k\right)$ (10)

3.2. 建模步骤

基于级比离散GM(1,1)模型对卫星钟差进行外推预报的具体步骤为

1) 钟差异常处理。利用中位数(median absolute deviation, MAD)方法对卫星钟差序列钟的粗差和跳变等异常进行检测和识别，将钟差异常值视为数据缺失点，并采用线性插值方法将缺失数据补全。

2) 根据式(3)对相邻历元钟差作一次级比生成一组级比序列，然后根据式(4)对级比生成序列进行累加生成一组累加序列，利用级比累加二次生成序列建立离散GM(1,1)模型。

3) 根据式(9)外推未来任意历元的钟差级比值，结合钟差序列和级比序列之间的关系式(10)，得到基于级比预报序列的卫星钟差外推函数：

${\hat{x}}^{\left(0\right)}\left(n+p\right)=x^{\left(0\right)}\left(n\right){\displaystyle \underset{k=1}{\overset{p}{\prod }}{\hat{c}}^{\left(0\right)}}\left(n+k-1\right)$ (11)

4. 试验结果与分析

4.1. 钟差数据来源

选取IGS提供的2021-03-28至2021-04-29最终精密GPS卫星钟差数据进行预报分析，钟差数据采样间隔为5 min。2021-03-28至2021-04-29时段GPS星座由32颗卫星构成，卫星钟类型包括BLOCK IIR 铷钟、BLOCK IIR-M铷钟、BLOCKIIF铷钟、BLOCKIIF铯钟与BLOCKIIIA铷钟5种，从每种卫星钟当中各选一颗卫星进行试验，选择的卫星为GPSPRN08、PRN17、PRN22、PRN23与PRN32，这5颗卫星的发射日期、入轨日期及卫星钟类型如表1所示。

Table 1. Information on the five GPS satellites

表1. 5颗GPS卫星的相关信息

4.2. 预报误差评估

通过将卫星钟差预报结果和IGS提供的最终精密钟差数据作差对预报误差进行分析，选择均方根误差(root mean square, RMS)与极差(range)作为预报误差评估指标，其中，RMS反映预报精度，range反映预报稳定性。RMS和range的计算公式分别为

$\text{RMS}=\sqrt{\frac{{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\left[{\hat{x}}^{\left(0\right)}\left(n+k\right)-x^{\left(0\right)}\left(n+k\right)\right]}^2}}{p}}$ (11)

$\text{range}=\left|{\text{error}}_{\max }-{\text{error}}_{\min }\right|$ (12)

式中， ${\hat{x}}^{\left(0\right)}\left(n+k\right)$ 和 $x^{\left(0\right)}\left(n+k\right)$ 分别为历元 $n+k$ 的卫星钟差预报值与真实值； ${\text{error}}_{\max }$ 和 ${\text{error}}_{\min }$ 分别为卫星钟差预报误差序列中的最大值与最小值。

4.3. 钟差级比分析

利用式(3)提取2021-03-28至2021-04-29时段GPSPRN08、PRN17、PRN22、PRN23与PRN32卫星钟差的级比序列，如图1所示。由图1可以看出，这5颗卫星钟差级比并不是恒定值，而是呈现动态变化，且不同卫星钟差级比变化趋势与变化幅度各不相同，其中，PRN08、PRN17与PRN22卫星钟差级比序列变化比较平稳，而PRN23与PRN32卫星钟差级比序列表现明显的降低趋势，变化趋势较为明显，所以选择的这5颗卫星钟差的级比序列具有充分的代表性。卫星钟差序列主要由趋势项分量和随机项分量构成，图1中每颗卫星钟差级比序列的线性变化可视为由钟差趋势项变化引起，而钟差随机性变化则导致了级比序列的非线性变化。传统GM(1,1)模型通过恒定级比对钟差进行拟合和外推存在局限性，因此本文利用钟差级比序列建立离散GM(1,1)模型对卫星钟差进行拟合和外推。

Figure 1. Calculated the step-ratio time-series from the bias of the five GPS satellite clocks

图1. 5颗卫星钟差计算级比序列

在试验中级比离散GM(1,1)模型是利用3 d的历史钟差级比序列去预报中长期的级比序列，再通过式(10)还原得到卫星钟差外推序列。2021-03-31至2021-04-29 PRN08、PRN17、PRN22、PRN23与PRN32卫星钟差的预报级比序列和真实级比序列如图2所示。由图2可以看出，级比离散GM(1,1)模型可以很好地预报出级比序列的变化趋势，表明基于级比序列建模的离散GM(1,1)模型对卫星钟差进行外推是合理的、可行的。

Figure 2. Extrapolated the step-ratio time-series from the bias of the five GPS satellite clocks

图2. 5颗卫星钟差预报级比序列

4.4. 预报算例分析

1) 1 d预报

为了验证级比离散GM(1,1)模型的卫星钟差预报性能，本文将其卫星钟差预报结果和二次多项式模型、常规GM(1,1)模型的钟差预报结果作对比分析。选择2021-03-28至2021-03-30时段的卫星钟差序列作为基础序列，对2021-03-31的卫星钟差进行预报，级比离散GM(1,1)模型、二次多项式模型与常规GM(1,1)模型的预报误差(真实序列与预报序列之差)如图3所示，图中，SDGM(1,1)代表级比离散GM(1,1)模型，QP表示二次多项式模型，三种模型的RMS和range统计结果见表2。

Figure 3. Errors for the 1-day predictions of the GM(1,1), QP and SDGM(1,1) models

图3. GM(1,1)、QP与SDGM (1,1)模型的1 d预报误差

Table 2. RMS and range statistics for1-day predictions of the GM(1,1), QP and SDGM(1,1) models

表2. GM(1,1)、QP和SDGM (1,1)模型的1 d预报RMS和range统计

① 由图3可以发现，二次多项式模型、GM(1,1)模型和级比离散GM(1,1)模型的预报误差均随着预报跨度的增大而变大，其中，GM(1,1)模型的预报误差发散现象非常明显，特别是对于PRN08卫星钟差预报，GM(1,1)模型的最大预报误差超过30 ns，这是由于该卫星钟差的级别序列变化幅度较大(见图1)，GM(1,1)模型用恒定级比去预报容易引起较大的预报误差，而级比离散GM(1,1)模型能够很好地预报出钟差级比序列的变化趋势，预报误差显著降低，5颗卫星钟差的预报误差控制在4 ns以内，表明级比离散GM(1,1)模型比常规GM(1,1)模型拥有更好的适应性。从图3还可看到，级比离散GM (1,1)模型在预报开始历元的误差非常接近于零，这是因为级比离散GM(1,1)模型的预报值是通过原始钟差序列最后一个历元时刻的钟差观测值与级比预报值累乘获得。

② 由表2可以看到，级比离散GM(1,1)模型对PRN08、PRN17、PRN22、PRN23与PRN32卫星钟差的预报RMS分别为1.14 ns、1.33 ns、0.53 ns、0.35 ns和0.26 ns，5颗卫星钟差的平均预报RMS为0.72 ns，预报精度显著提高，相比于二次多项式模型和GM(1,1)模型，平均预报RMS分别降低了37.96%和62.38%；PRN08、PRN17、PRN22、PRN23与PRN32卫星钟差的预报range分别为4.94 ns、4.53 ns、1.69 ns、0.77 ns和0.66 ns，平均预报range为2.52 ns，相比于二次多项式模型和GM(1,1)模型，平均预报range分别降低了10.87%和26.6%，说明级比离散GM(1,1)模型不仅比另外两种模型的预报精度更高，并且预报稳定度也高于另外两种模型。

2) 30 d预报

为进一步评估分析级比离散GM(1,1)模型的钟差中长期预报效果，利用2021-03-28至2021-03-30共3的历史钟差序列来预报2021-03-31至2021-04-29时段的卫星钟差，级比离散GM(1,1)模型、二次多项式模型与常规GM(1,1)模型的预报误差如图4所示，三种模型对5颗卫星钟差的平均预报RMS和平均预报range统计结果如图5，图6所示。从图4~6可以看出：

① GM(1,1)模型的预报误差随着预报时长的增加而快速增大，预报误差发散现象严重，对级比变化幅度较大的PRN08卫星钟差，GM(1,1)模型1~30 d的最大预报误差超过1500 ns，对于另外4颗卫星钟差预报，其1~30d的最大预报误差可以达到800 ns，说明GM(1,1)模型利用恒定的级比取代变化的级比建立指数函数模型来预报卫星钟差，无法体现钟差级比序列的动态变化特性，在钟差中长期预报中可能预报失效。基于级比序列建模的离散GM(1,1)能够较好地预报出钟差级比的变化趋势，预报误差明显降低，对PRN08、PRN17、PRN22、PRN23与PRN32卫星钟差的最大绝对预报误差分别控制在430 ns、390 ns、106 ns、20 ns和46 ns，显著小于GM(1,1)模型和二次多项式模型。

② 级比离散GM(1,1)模型的1~10 d、1~20 d和1~30 d的平均预报RMS和平均预报range均明显低于GM(1,1)模型和二次多项式模型，1~10 d、1~20 d和1~30 d预报的平均RMS相比于二次多项式模型分别降低68.15%、76.04%和75.18%，1~10 d和1~30 d预报的平均RMS比GM(1,1)模型低1个数量级，1~20 d预报的平均RMS比GM(1,1)模型低2个数量级；对于预报稳定度，级比离散GM(1,1)模型的1~10 d、1~20 d和1~30预报的平均range相比于二次多项式模型分别降低68.35%、75.12%和76.74%，比GM(1,1)模型的平均预报range高1~2个数量级。从图4还可以发现，级比离散GM(1,1)模型对PRN08和PRN17卫星钟差的预报效果不如PRN22、PRN23与PRN32卫星钟差，原因在于PRN08和PRN17卫星钟差的级别变化幅度大于PRN22、PRN23与PRN32卫星钟差，导致预报误差相对较大。总体而言，级比离散GM(1,1)模型对于卫星钟差中长期预报的适应性由于二次多项式和常规GM(1,1)模型。

Figure 4. Errors for the 30-day predictions of the GM(1,1), QP and SDGM(1,1) models

图4. GM(1,1)、QP与SDGM (1,1)模型的30 d预报误差

Figure 5. Mean RMS for the predictions of the GM(1,1), QP and SDGM(1,1) models

图5. GM(1,1)、QP与SDGM (1,1)模型的平均预报RMS

Figure 6. Mean range for the predictions of the GM(1,1), QP and SDGM(1,1) models

图6. GM(1,1)、QP与SDGM (1,1)模型的平均预报range

5. 结论

为提升卫星钟差的中长期预报效果，提出基于钟差级比序列建模的卫星钟差预报方法，该方法利用变化的钟差级比代替恒定的级比建立离散GM(1,1)模型对级比进行外推，再结合钟差序列与级比序列之间的关系得到钟差预报值，这种方法考虑了钟差级比的动态变化特性。卫星钟差预报试验结果表明，级比离散GM(1,1)模型能够很好地反映卫星钟差级比序列的变化趋势，其预报精度和预报稳定度优于常用的二次多项式模型和常规GM(1,1)模型。特别是对于卫星钟差中长期预报，级比离散GM(1,1)模型的预报性能显著优于另外两种常用模型。此外，级比离散GM(1,1)模型解决了当钟差级比波动幅度较大时常规GM(1,1)模型预报失效的问题，比常规GM(1,1)模型具有更好的普适性，为卫星钟差中长期预报提供一种新方法。

致谢

感谢国际GNSS服务组织(IGS)提供的GPS卫星钟差数据，感谢审稿人提出的宝贵意见。

基金项目

本文受徐州市重点研发计划项目(KC18079)资助。

NOTES

^*通讯作者。