CN102306302B - 基于emd与garch模型的卫星时钟误差预报方法 - Google Patents

Abstract

基于EMD与GARCH模型的卫星时钟误差预报方法，涉及卫星时钟的钟差预测领域，解决了现有的卫星时钟误差预报方法因缺少非平稳随机项的预测，而导致预报精度难以提高的问题，它包括步骤一、获取时钟误差历史数据，并进行数据修正预处理得到卫星钟差数据；步骤二、对卫星钟差数据的经验模态进行分解，得到卫星钟差数据的随机项部分；步骤三、对卫星钟差数据的趋势项进行预测，通过建立卡尔曼预报模型进行卫星钟差数据趋势项的预测；步骤四、对卫星钟差数据的随机项进行预测，去除趋势项后得到的随机项，采用ARMA和GARCH模型对随机项进行预测，提高卫星钟差预测精度。用于卫星导航系统的高精度时间同步。

Description

基于EMD与GARCH模型的卫星时钟误差预报方法

技术领域

本发明涉及卫星时钟的钟差预测领域。

背景技术

时间同步是卫星导航系统的基础和关键，时间同步的最终精度取决于卫星钟差预报的精度，即卫星钟差预报误差是最大的误差源之一，因而，卫星钟差预报是卫星导航系统关键技术之一。

钟差预报就是利用当前已有的钟差观测数据，通过一定算法，得到未来所需某个时刻的钟差预报值。基于经验模态分解(EMD)与广义自回归条件异方差(GARCH)模型的卫星时钟误差预报方法是一种高精度、快速的钟差预测技术，这种方法利用EMD将卫星钟差数据分解为趋势项和随机项，根据卡尔曼滤波器对趋势项数据进行预测，根据自回归滑动平均(ARMA)模型和GARCH模型分别对随机项的平稳部分和非平稳部分进行预测，从而提高卫星时钟的钟差预报精度。

现有的卫星时钟误差预报方法因缺少非平稳随机项的预测，而导致预报精度难以提高。

发明内容

本发明的目的是为了解决现有的卫星时钟误差预报方法因缺少非平稳随机项的预测，而导致预报精度难以提高的问题，提供一种基于EMD与GARCH模型的卫星时钟误差预报方法。

基于EMD与GARCH模型的卫星时钟误差预报方法，它包括以下具体步骤：

步骤一、获取卫星时钟误差历史数据，并进行数据修正预处理得到卫星钟差数据；

步骤二、对卫星钟差数据的经验模态进行分解；通过EMD将卫星钟差数据分解为一系列具有不同频率成分的本征模态函数和残差函数分量；去除高频部分后，对本征模态函数取和得到卫星钟差数据的随机项部分，通过残差函数获得卫星钟差数据的趋势项；

步骤三、对卫星钟差数据的趋势项进行预测，根据卡尔曼滤波理论，通过建立卡尔曼预报模型进行卫星钟差数据趋势项的预测；完成卫星钟差数据趋势项的预报；

步骤四、对卫星钟差数据的随机项进行预测，卫星钟差数据去除趋势项后得到随机项，采用ARMA和GARCH模型对随机项进行预测，完成卫星钟差数据随机项的预报。

本发明的方法增加了非平稳随机项的预测，是一种新的高精度的钟差预测方法。

本发明与现有钟差预测技术相比优点在于：

1、目前钟差的预测技术大多数都是仅在趋势项上建立钟差函数模型，而不考虑钟差的随机项，少数考虑随机项建模也是主要采用ARMA模型，且存在ARMA模式识别困难的问题。本发明采用EMD与随机项建模方法相结合，解决模式识别困难的问题，并进一步提高了钟差预报精度；

2、相比当前主流的基于最小二乘和灰色系统模型等进行钟差趋势项长期预报技术而言，本发明提出的卡尔曼滤波器预测方法是一种递推算法，不需要保存过去的观测数据，因此，只要卫星时钟不进行调整，就保证卡尔曼预报有充足的数据，适合于预报时间小于一天的短期预报；

3、通常情况下仅采用ARMA模型对钟差随机项进行预测，本发明中对随机项的预测过程中考虑了实际问题中随着时间变化，随机扰动项的条件方差的变化，不仅采用ARMA模型对钟差随机项进行预测，还采用GARCH模型对随机扰动项进行预测，对比图5、图6，本发明中提出的钟差预测方法的预测精度比传统方法提高70％。

附图说明

图1为本发明的流程示意图，图2为本发明对获取的15天的时钟误差历史数据进行数据修正预处理后得到的卫星钟差数据，图3为本发明的时钟误差的趋势项，图4为本发明的时钟误差的随机项，图5为传统方法的钟差预测误差图，图6为本发明方法的钟差预测误差图。

具体实施方式

具体实施方式一：结合图1说明本实施方式，本实施方式包括具体步骤如下：

步骤一、获取时钟误差历史数据，并进行数据修正预处理得到卫星钟差数据；

步骤二、对卫星钟差数据的经验模态进行分解；通过EMD将卫星钟差数据分解为一系列具有不同频率成分的本征模态函数和残差函数分量；本征模态函数的频率具有从高到低的特点，去除高频部分后，对本征模态函数取和得到卫星钟差数据的随机项部分，通过残差函数获得卫星钟差数据的趋势项，而残差函数反映了卫星钟差数据的趋势项部分；

步骤三、对卫星钟差数据的趋势项进行预测，由EMD中残差函数得到的趋势项的走势具有伪线性特性，根据卡尔曼滤波理论，通过建立卡尔曼预报模型进行卫星钟差数据趋势项的预测；完成卫星钟差数据趋势项的预报；

步骤四、对卫星钟差数据的随机项进行预测，卫星钟差数据去除趋势项后得到随机项，通常情况下，该随机项为非平稳时间序列，采用ARMA和GARCH模型对随机项进行预测，提高卫星钟差数据的预测精度，完成卫星钟差数据随机项的预报。

具体实施方式二：具体实施方式一中的步骤一获取的时钟误差历史数据中如果存在异常点、无数据段或数据跳变的数据异常，对所述异常进行修正，对于异常点，采取删除异常点，然后合并为无数据段；对于无数据段采用多项式插值的方法得到无数据段的数据；对于数据跳变，采用滑动窗口探测出跳变，对于跳变数据采取舍去前端数据、用后段数据做钟差预测。其它组成和连接关系与实施方式一相同。

步骤二在卫星钟差数据的EMD分析前判断是否同时满足以下条件：a、被分析的数据至少包含极大值和极小值两个极值点；b、根据两个相邻极值点的时间距离能定义特征时间尺度；c、如果数据无极值点但包含拐点，则所述数据能通过一次或多次求导得到极值点。

从上述三个判断条件可以看出，EMD分析需要确定卫星钟差数据的极值点；卫星钟差数据在常值漂移和随机漂移的影响下，可以满足判断条件中对极值点的要求；为了保证EMD分析过程中所采用的极值点是卫星钟差数据的真实极值点，采用足够高的采样频率能够满足EMD分析的判断条件，能够保证EMD分析的准确性。

具体实施方式三：具体实施方式一中的步骤二包括的子步骤如下：

子步骤1、对卫星钟差数据y(t)求取极大值y(t_u)和极小值y(t_v)，其中u＝1，2，...，N_u，v＝1，2，...，N_v，N_u为极大值的个数，N_v为极小值的个数；采用三次样条函数分别构造极大值点和极小值点的上下包络线y_u(t)和y_v(t)，计算两个包络线的均值

$m_1=\frac{1}{2}(y_u\left(t\right)+y_v\left(t\right));---\left(1\right)$

子步骤2、判断h₁＝y(t)-m₁是否同时满足如下的作为IMF的两个条件：

a、数据极值点的数量与零点数量相等或相差一个；

b、数据由极大值定义的上包络和由极小值定义的下包络的局部均值为零；

如果h₁满足上述IMF条件，那么h₁为y(t)的第一个IMF分量c₁，并得到第一个残差r₁＝y(t)-c₁；

如果h₁不满足上述IMF条件，则将h₁作为新的数据，重复子步骤1求取均值m₁₁，进而求数据h₁与m₁₁的差值h₁₁；对h₁₁重复上述过程k次，直到h_1k满足IMF条件，那么h_1k为y(t)的第一个IMF分量c₁，并得到第一个残差r₁＝y(t)-h_1k；

子步骤3、把r₁作为新的数据，重复上面的步骤；依次得到IMF分量c₂，...，c_n和残差r₂，...，r_n；直到r_n为单调数据或者只存在一个极点为止；

综上，经过EMD处理后的卫星钟差数据y(t)表示为

$y\left(t\right)=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{c}_k+r_n=C_n+r_n---\left(2\right)$

本征模态函数c_k的频率具有从高到低的特点，对本征模态函数c_k取和得到C_n，C_n为卫星钟差数据的随机项，而残差函数r_n反映了卫星钟差数据的趋势项部分。

具体实施方式四：具体实施方式一中的步骤三中采用卡尔曼滤波理论对卫星钟差数据的趋势项进行预测的过程为：

将卫星钟差数据的趋势项r_n表示为

r_n(t)＝a+bt+ct²+Δε₁(t)+Δε₂(t)    (3)

其中，a，b，c为二次多项式系统误差参数，Δε₁(t)为相位噪声，Δε₂(t)为测量噪声；

建立如下的卡尔曼预测模型

r_n(t)＝x₁(t)+n₁(t)+n₂(t)

${\stackrel{\·}{x}}_1\left(t\right)=x_2\left(t\right)+n_3\left(t\right)$ (4)

${\stackrel{\·}{x}}_2\left(t\right)=x_3\left(t\right)+n_4\left(t\right)$

${\stackrel{\·}{x}}_3\left(t\right)=0$

其中，n_i(t)，i＝1，2，3，4；为相互独立的白噪声，并设噪声的方差为

当采样时间为τ时，离散化后的状态转移模型和观测模型为：

$\left[\begin{array}{c}{x}_1(k+1)\\ x_2(k+1)\\ x_3(k+1)\end{array}\right]=\left[\begin{array}{ccc}1& \mathrm{\τ}& \frac{1}{2}{\mathrm{\τ}}^2\\ 0& 1& \mathrm{\τ}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_1\left(k\right)\\ x_2\left(k\right)\\ x_3\left(k\right)\end{array}\right]+\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ 0\end{array}\right]---\left(5\right)$

r_n(k)＝x₁(k)+n₁(k)+n₄(k)    (6)

其中，噪声方差阵为

$Q=\left[\begin{array}{ccc}{\mathrm{\τ\σ}}_3^2+\frac{1}{3}{\mathrm{\τ}}^2{\mathrm{\σ}}_4^2& \frac{1}{2}{\mathrm{\τ}}^2{\mathrm{\σ}}_4^2& 0\\ \frac{1}{2}{\mathrm{\τ}}^2{\mathrm{\σ}}_4^2& {\mathrm{\τ\σ}}_4^2& 0\\ 0& 0& 0\end{array}\right]---\left(7\right)$

测量方差阵为

$R={\mathrm{\σ}}_1^2\left(k\right)+{\mathrm{\σ}}_2^2\left(k\right)---\left(8\right)$

设预测时刻距离最后一个滤波值的时间间隔为nτ，则预测值为

$r_n\left((k+n)\mathrm{\τ}\right)=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]{\left[\begin{array}{ccc}1& \mathrm{\τ}& \frac{1}{2}{\mathrm{\τ}}^2\\ 0& 1& \mathrm{\τ}\\ 0& 0& 1\end{array}\right]}^n\left[\begin{array}{c}{x}_1\left(k\right)\\ x_2\left(k\right)\\ x_3\left(k\right)\end{array}\right]---\left(9\right)$

式(5)至(9)就是完整的卫星钟差数据趋势项卡尔曼预测算法模型，根据式(9)得到卫星钟差数据趋势项的实时预测值。

具体实施方式五：具体实施方式一中的步骤四的具体过程为：

建立卫星钟差数据随机项C_n的ARMA模型

其中，

是待估的自回归参数；θ_j是移动平均参数；r是自回归阶数；m是移动平均的阶数；ε(t)是第t时刻的残差；c是常数；

ARMA模型的前提假设是模型扰动均值为零，方差为常数；而随机项的扰动部分方差是随时间变化的时间序列，采用广义自回归条件异方差模型进行处理；标准的GARCH(p，q)为：

${\mathrm{\σ}}_t^2=k+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{G}_i{\mathrm{\σ}}_{t-i}^2+\underset{j=1}{\overset{q}{\mathrm{\Σ}}}{A}_j{\mathrm{\ϵ}}_{t-j}^2---\left(11\right)$

其中，G_i＞0是GARCH影响的待估参数；A_j＞0是ARCH影响的待估参数；p和q是模型阶次；

是以前面数据为基础的一期向前预测方差，即条件方差；k是常数；

在选择模型的阶次时，采用赤池信息准则(AIC)来衡量统计模型拟合的优良性；通过选择较小的AIC值来控制模型的复杂程度，AIC表示为

x_AIC＝-2y_LLF+2N    (12)

其中，N是GARCH模型需要估计参数的个数，它影响模型的自由度，y_LLF为极大似然函数；

结合式(10)和(11)构成了考虑异方差性的卫星钟差数据随机项预测模型，在一定程度上提高卫星钟差数据随机项的建模精度。

基于EMD和GARCH模型的卫星时钟误差预报方法能实现卫星的高精度钟差预测。

Claims (4)

1.基于EMD与GARCH模型的卫星时钟误差预报方法，其特征是它包括具体步骤如下：

步骤一、获取卫星时钟误差历史数据，并进行数据修正预处理得到卫星钟差数据；

步骤二、对卫星钟差数据的经验模态进行分解；通过经验模态分解EMD将卫星钟差数据分解为一系列具有不同频率成分的本征模态函数和残差函数分量；去除高频部分后，对本征模态函数取和得到卫星钟差数据的随机项部分，通过残差函数获得卫星钟差数据的趋势项；

步骤三、对卫星钟差数据的趋势项进行预测，根据卡尔曼滤波理论，通过建立卡尔曼预报模型进行卫星钟差数据趋势项的预测，完成卫星钟差数据趋势项的预报；

步骤四、对卫星钟差数据的随机项进行预测，卫星钟差数据去除趋势项后得到随机项，采用自回归滑动平均ARMA和广义自回归条件异方差GARCH模型对随机项进行预测，完成卫星钟差数据随机项的预报。

2.根据权利要求1所述基于EMD与GARCH模型的卫星时钟误差预报方法，其特征在于步骤一获取的时钟误差历史数据中如果存在异常点、无数据段或数据跳变的数据异常，对所述异常进行修正，对于异常点，采取删除异常点，然后合并为无数据段；对于无数据段采用多项式插值的方法得到无数据段的数据；对于数据跳变，采用滑动窗口探测出跳变，对于跳变数据采取舍去前端数据、用后段数据做钟差预测。

3.根据权利要求1所述基于EMD与GARCH模型的卫星时钟误差预报方法，其特征在于步骤二包括的子步骤如下：

子步骤1、对卫星钟差数据y(t)求取极大值y(t_u)和极小值y(t_v)，其中u=1,2,...,N_u；v=1,2,...,N_v；N_u为极大值的个数，N_v为极小值的个数；采用三次样条函数分别构造极大值点和极小值点的上下包络线y_u(t)和y_v(t)，计算两个包络线的均值

$m_1=\frac{1}{2}(y_u\left(t\right)+y_v\left(t\right));$

子步骤2、判断h₁=y(t)-m₁是否同时满足如下的作为IMF的两个条件：

a、数据极值点的数量与零点数量相等或相差一个；

b、数据由极大值定义的上包络和由极小值定义的下包络的局部均值为零；

如果h₁满足上述IMF条件，那么h₁为y(t)的第一个IMF分量c₁，并得到第一个残差r₁=y(t)-c₁；

如果h₁不满足上述IMF条件，则将h₁作为新的数据，重复子步骤1求取均值m₁₁，进而求数据h₁与m₁₁的差值h₁₁；对h₁₁重复上述过程k次，直到h_1k满足IMF条件，那么h_1k为y(t)的第一个IMF分量c₁，并得到第一个残差r₁=y(t)-h_1k；

子步骤3、把r₁作为新的数据，重复上面的步骤；依次得到IMF分量c₂,...,c_n和残差r₂,...,r_n，直到r_n为单调数据或者只存在一个极点为止；

综上，经过EMD处理后的卫星钟差数据y(t)表示为

$y\left(t\right)=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{c}_k+r_n=C_n+r_n$

C_n是卫星钟差数据的随机项，r_n是卫星钟差数据的趋势项。

4.根据权利要求1所述基于EMD与GARCH模型的卫星时钟误差预报方法，其特征在于步骤四的具体过程为：

建立卫星钟差数据随机项C_n的ARMA模型

其中，

是待估的自回归参数；θ_j是移动平均参数；r是自回归阶数；m是移动平均的阶数；

ε(t)是第t时刻的残差；c是常数；

ARMA模型的前提假设是模型扰动均值为零，方差为常数；而随机项的扰动部分方差是随时间变化的时间序列，采用广义自回归条件异方差模型进行处理；标准的GARCH(p,q)为：

${\mathrm{\σ}}_t^2=k+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{G}_i{\mathrm{\σ}}_{t-i}^2+\underset{j=1}{\overset{q}{\mathrm{\Σ}}}{A}_j{\mathrm{\ϵ}}_{t-j}^2$

其中，G_i>0是GARCH影响的待估参数；A_j>0是ARCH影响的待估参数；p和q是模型阶次；

是以前面数据为基础的一期向前预测方差，即条件方差；k是常数；

在选择模型的阶次时，采用赤池信息准则AIC来衡量统计模型拟合的优良性；通过选择较小的AIC值来控制模型的复杂程度，AIC表示为

x_AIC=-2y_LLF+2N

其中，N是GARCH模型需要估计参数的个数，它影响模型的自由度，y_LLF为极大似然函数。

Images