CN104501804A - 一种基于gps测量数据的卫星在轨轨道预报方法 - Google Patents

Abstract

本发明提出了一种基于GPS测量数据的卫星在轨轨道预报方法，主要针对低轨近圆卫星轨道给出了一种卫星解析星历参数模型，利用GPS测量信息对所述星历模型参数进行在轨估计，可对任意时间间隔的卫星星历进行预报。仿真结果表明，对于450km以上的近地圆轨道，预报精度可达20km以内，具有很好的实用性。

Description

一种基于GPS测量数据的卫星在轨轨道预报方法

技术领域

本发明涉及天体力学中的卫星轨道预报技术领域，具体涉及一种基于GPS测量数据的卫星在轨轨道预报方法。

背景技术

随着目前国内外卫星技术的不断发展尤其是卫星组网的发展，对于卫星在轨自主能力的需求不断增加，在轨实时轨道确定成为判断卫星是否具有自主能力的首要条件。随着低成本全球导航***接收机(包括美国的GPS及中国的北斗导航***)的应用，在微小卫星上进行实时轨道确定进而提高小卫星的自主能力，已成为目前的一种发展趋势。

目前基于GPS测量信息的卫星实时在轨轨道确定主要采用Kalman滤波及基于轨道动力学模型的轨道确定技术，对卫星的瞬时轨道状态(位置和速度矢量)或密切轨道根数进行估计，以实时提供卫星的高精度轨道确定信息。其结果主要用于卫星各种实时应用，例如图像的地理定位编码、敏感器与可驱动天线的指向以及卫星三轴姿态控制。但对于采用这种轨道确定技术的导航***，其结果很难用于长期的在轨星历预报(例如几个轨道周期后甚至几天后的轨道预报)，而这一点恰恰是卫星自主任务规划和自主管理所需要的。其主要原因是：基于瞬时轨道参数或密切轨道根数进行轨道预报，需要进行复杂的、计算量较大的轨道动力学数值积分运算过程，往往需要占用大量的星载计算机机时，因而不适于进行长期轨道预报。

发明内容

为了解决现有技术中存在的缺陷，本发明提供了一种基于GPS测量数据的卫星在轨轨道预报方法，该方法针对目前常用的低轨近圆卫星轨道，基于解析的轨道动力学模型，采用Kalman滤波技术利用GPS测量信息对相关星历模型参数进行在轨估计。可对任意时间间隔的卫星星历进行预报，而不用按步长对轨道进行积分，其计算量相对较小，可对卫星进行中期或长期的轨道预报(多个轨道周期或一周以上)。

为了达到上述目的，本发明采取了以下技术方案：

一种基于GPS测量数据的卫星在轨轨道预报方法，包括以下步骤：

S1：根据初始化标志判断当前计算状态，是否需要进行初始化，若需要进行初始化，则执行S2，否则执行S4；

S2：根据初始化方法标志进行滤波器初始化，其中初始化方法分为三种：a)根据GPS测量数据进行滤波器初始化，即根据GPS提供的测量数据，计算所需测量时刻对应的轨道平根数集的初始猜测；b)根据地面上注的星历模型参数进行初始化；c)用上一次滤波器计算收敛得到的星历参数进行初始化；

S3：初始化完成后，将初始化标志置位，回到S1；

S4：进行Kalman滤波计算，对星历模型参数进行滤波，得到更新的星历参数。

进一步地，所述Kalman滤波计算为Bierman-UD分解形式的Kalman滤波。

进一步地，所述步骤S1中的初始化标志为外部给出的是否进行滤波器初始化的标志。

进一步地，所述步骤S2中的初始化方法标志预先设置，在轨道的某一段时间内没有GPS测量数据情况下，初始化方法标志设置为a。

进一步地，所述星历模型参数为星历模型参数对应的星上时间t_k、平均轨道半长轴平均轨道倾角平均升交点经度偏心率向量长周期变化的幅值偏心率向量长周期变化的相位平均纬度角升交点赤经的一次长期项摄动系数Ω₁、近地点幅角的一次长期项摄动系数ω₁、平均纬度角一次长期摄动系数λ₁、平均半长轴一次长期摄动系数a_d、偏心率向量x分量ξ的一次长期摄动系数ξ_d、偏心率向量y分量η的一次长期摄动系数η_d、偏心率向量固定偏移e_f、平均纬度角的二次长期摄动系数λ₂。

进一步地，所述步骤S4后还包括步骤S5:根据更新的星历参数基于无奇点变量的拟平均根数法计算出任意时刻卫星的瞬时轨道根数，进而得到卫星在瞬时真赤道TOD坐标系中的位置和速度。

本发明针对卫星长期自主运行需求，提出了一种基于GPS测量数据的卫星自主在轨轨道预报方法，该方法基于解析的轨道动力学模型，采用Kalman滤波技术，可对卫星进行中期或长期的轨道预报(多个轨道周期或一周以上)。具有以下优点：(1)采用基于无奇点变量的拟平均根数法，具有自启动(自初始化)、收敛性好，且大大简化了星上计算量；(2)对测量数据的采样要求不严格(允许测量数据以非均匀间隔时间给出，采样时间可在几秒～10多分钟内变化)的优点；(3)对于450km以上的近地近圆轨道，在滤波2.5天的情况下，预报精度在20km以内，能够满足一般的任务需求；(4)允许在轨道的某一段时间内没有测量数据情况下(例如GPS天线被地球遮挡，或接收机出现临时性故障)，滤波器仍然能正常运行。

附图说明

图1是本发明的方法的流程图；

图2是本发明的方法在轨道高度450km的预报精度仿真结果图；

图3是本发明的方法在轨道高度600km的预报精度仿真结果图；

图4是本发明的方法在轨道高度700km的预报精度仿真结果图；

图5是本发明的方法在轨道高度800km的预报精度仿真结果图。

具体实施方式

下面结合附图说明及具体实施方式对本发明进一步说明。

本发明采用基于无奇点变量的拟平均根数法和Kalman滤波技术，利用GPS测量信息对相关星历模型参数进行在轨估计，实现卫星的在轨轨道预报。为减少滤波过程中计算误差和舍入误差对滤波器状态误差协方差矩阵正定性的影响，避免长时间计算过程滤波器发散，可采用Bierman-UD分解形式的Kalman滤波技术。

根据定义的不同，平均根数又可以分为平均根数和拟平均根数，其中平均根数只包含了长期变化项，拟平均根数包含长期变化项和长周期变化项，长周期变化项的周期能达到数个月，因而在轨道问题研究中其影响是不可忽略的，同时，拟平均根数能消除通约奇点问题，因而本发明采用基于无奇点变量的拟平均根数法，为简化星上计算量，计算过程仅考虑了一阶精度。

近地卫星主要受到地球引力、大气阻力、太阳光压以及日月引力等作用，因而轨道变化包括长期变化项、长周期变化项和短周期变化项。长期变化即轨道参数随时间的线性或二阶及更高阶的变化，长周期变化主要是由地球引力场的带谐调和项引起的，短周期项一类是由地球扁率引起的，一类是由地球引力场的田谐调和项引起的。为了简化计算，在模型中仅考虑长期和长周期变化项。

在星历模型的设计中，通常是给出能够描述描述轨道参数变化的数学模型，根据这个数学模型可以计算出任意时刻卫星的瞬时轨道根数，进而可得到卫星在瞬时真赤道(TOD)坐标系中的位置和速度。本发明的方法针对的是近圆轨道，因而选用以下无奇点变量：

a,i,ξ＝ecosω,η＝esinω,Ω,λ＝M+ω

上式中a,i,e,ω,Ω是五个经典Kepler要素，M是平近点角，为了精确描述轨道变化，再引入表征轨道长期变化和长周期变化的固定参数，同时为了与国内采用的卫星解析星历模型相匹配，本发明采用了目前国内卫星普遍使用的星历模型参数形式(15参数)，其定义为:

t_k——星历模型参数对应的星上时间；

——平均轨道半长轴；

——平均轨道倾角；

——平均升交点经度；

——偏心率向量长周期变化的幅值；

——偏心率向量长周期变化的相位；

——平均纬度角；

Ω₁——升交点赤经的一次长期项摄动系数；

ω₁——近地点幅角的一次长期项摄动系数；

λ₁——平均纬度角一次长期摄动系数；

a_d——平均半长轴一次长期摄动系数；

ξ_d——偏心率向量x分量ξ的一次长期摄动系数；

η_d——偏心率向量y分量η的一次长期摄动系数；

e_f——偏心率向量固定偏移；

λ₂——平均纬度角的二次长期摄动系数。

本发明的方法的流程如附图1所示。

S1：根据初始化标志判断当前计算状态，是否需要进行初始化，若需要进行初始化，则执行S2，否则执行S4；

S2：根据初始化方法标志进行滤波器初始化，其中初始化方法分为三种：a)根据GPS测量数据进行滤波器初始化，即根据GPS提供的测量数据，计算所需测量时刻对应的轨道平根数集的初始猜测；b)根据地面上注的星历模型参数进行初始化；c)用上一次滤波器计算收敛得到的星历参数进行初始化；

S3：初始化完成后，将初始化标志置位，回到S1；

S4：进行Kalman滤波计算，对星历模型参数进行滤波，得到更新的星历参数。

利用以上选取的15个参数可以计算出任意时刻卫星在TOD坐标系中的位置R和速度V，具体计算过程如下所示：

1)t₁时刻平根数计算

$\stackrel{\‾}{a}={\stackrel{\‾}{a}}_0+{\stackrel{\‾}{a}}_d(t_1-t_0)$

$\stackrel{\‾}{\mathrm{\Ω}}={\stackrel{\‾}{\mathrm{\Ω}}}_0+{\stackrel{\‾}{\mathrm{\Ω}}}_1(t_1-t_0)$

$\stackrel{\‾}{\mathrm{\ξ}}={\stackrel{\‾}{l}}_a\cos [{\mathrm{\ω}}_1(t_1-t_0)-{\stackrel{\‾}{l}}_p]+{\mathrm{\ξ}}_d(t_1-t_0)$

$\stackrel{\‾}{\mathrm{\η}}=-{\stackrel{\‾}{l}}_a\sin [{\mathrm{\ω}}_1(t_1-t_0)-{\stackrel{\‾}{l}}_p]+{\mathrm{\η}}_d(t_1-t_0)+e_f$

$\stackrel{\‾}{\mathrm{\λ}}={\stackrel{\‾}{\mathrm{\λ}}}_0+{\mathrm{\λ}}_1(t_1-t_0)+{\mathrm{\λ}}_2{(t_1-t_0)}^2$

$\stackrel{\‾}{p}=\stackrel{\‾}{a}(1-{\stackrel{\‾}{\mathrm{\ξ}}}^2-{\stackrel{\‾}{\mathrm{\η}}}^2)$

2)计算纬度幅角的平均值

${\stackrel{\‾}{\mathrm{\ξ}}}_0={\stackrel{\‾}{l}}_a\cos \left({\stackrel{\‾}{l}}_p\right),{\stackrel{\‾}{\mathrm{\η}}}_0={\stackrel{\‾}{l}}_a\sin \left({\stackrel{\‾}{l}}_p\right)+e_f$

${\stackrel{\‾}{e}}_0=\sqrt{{\stackrel{\‾}{\mathrm{\ξ}}}_0^2+{\stackrel{\‾}{\mathrm{\η}}}_0^2}$

$\tilde{u}=\stackrel{\‾}{\mathrm{\λ}}$

for i＝1:6,

$\tilde{u}=\stackrel{\‾}{\mathrm{\λ}}+\stackrel{\‾}{\mathrm{\ξ}}\sin \tilde{u}+\stackrel{\‾}{\mathrm{\η}}\cos \tilde{u}$

end

$\stackrel{\‾}{u}={\sin }^{-1}\left[\frac{\tilde{u}-\stackrel{\‾}{\mathrm{\λ}}}{1-\stackrel{\‾}{\mathrm{\ξ}}\cos \tilde{u}+\stackrel{\‾}{\mathrm{\η}}\sin \tilde{u}}(1+\frac{\stackrel{\‾}{\mathrm{\η}}\sin \tilde{u}-\stackrel{\‾}{\mathrm{\ξ}}\cos \tilde{u}}{1+\sqrt{1-{{\stackrel{\‾}{e}}_0}^2}})\right]+\tilde{u}$

3)t₁时刻短周期项的计算

$a_s=\frac{A_2}{\stackrel{\‾}{a}}{\sin }^2{\stackrel{\‾}{i}}_0\cos 2\stackrel{\‾}{u}$

$i_s=\frac{A_2}{2{{}_{}\stackrel{\‾}{p}}^2}\sin {\stackrel{\‾}{i}}_0\cos {\stackrel{\‾}{i}}_0\cos 2\stackrel{\‾}{u}$

${\mathrm{\Ω}}_s=\frac{A_2}{2{\stackrel{\‾}{p}}^2}\cos {\stackrel{\‾}{i}}_0\sin 2\stackrel{\‾}{u}$

${\mathrm{\ξ}}_s=\frac{A_2}{{\stackrel{\‾}{p}}^2}[(1-\frac{5}{4}s{\mathrm{in}}^2{\stackrel{\‾}{i}}_0)\cos \stackrel{\‾}{u}+\frac{7}{12}{\sin }^2{\stackrel{\‾}{i}}_0\cos 3\stackrel{\‾}{u}]$

${\mathrm{\η}}_s=\frac{A_2}{{\stackrel{\‾}{p}}^2}[(-1+\frac{7}{4}s{\mathrm{in}}^2{\stackrel{\‾}{i}}_0)\sin \stackrel{\‾}{u}-\frac{7}{12}{\sin }^2{\stackrel{\‾}{i}}_0\sin 3\stackrel{\‾}{u}]$

${\mathrm{\λ}}_s=\frac{A_2}{{\stackrel{\‾}{p}}^2}(-\frac{1}{2}+\frac{5}{4}s{\mathrm{in}}^2{\stackrel{\‾}{i}}_0)\sin \stackrel{}{}2\stackrel{\‾}{u}$

4)t₁时刻瞬时根数的计算

$a=\stackrel{\‾}{a}+a_s$

$i=\stackrel{\‾}{i}+i_s$

$\mathrm{\Ω}=\stackrel{\‾}{\mathrm{\Ω}}+{\mathrm{\Ω}}_s$

$\mathrm{\ξ}=\stackrel{\‾}{\mathrm{\ξ}}+{\mathrm{\ξ}}_s$

$\mathrm{\η}=\stackrel{\‾}{\mathrm{\η}}+{\mathrm{\η}}_s$

$\mathrm{\λ}=\stackrel{\‾}{\mathrm{\λ}}+{\mathrm{\λ}}_s$

$e=\sqrt{{\mathrm{\ξ}}^2+{\mathrm{\η}}^2}$

ω＝atan2(-η,ξ)

M＝λ-ω

u＝λ+2esinM+1.25e²sin2M

f＝u-ω

5)t₁时刻卫星位置和速度的计算

$r=\frac{a(1-e^2)}{1+e\cos f}$

$\underset{\‾}{P}=\left[\begin{array}{c}\cos \mathrm{\ω}\cos \mathrm{\Ω}-\sin \mathrm{\ω}\cos i\sin \mathrm{\Ω}\\ \cos \mathrm{\ω}\sin \mathrm{\Ω}+\sin \mathrm{\ω}\cos i\cos \mathrm{\Ω}\\ \sin \mathrm{\ω}\sin i\end{array}\right]$

$\underset{\‾}{Q}=\left[\begin{array}{c}-\sin \mathrm{\ω}\cos \mathrm{\Ω}-\cos \mathrm{\ω}\cos i\sin \mathrm{\Ω}\\ -\sin \mathrm{\ω}\sin \mathrm{\Ω}+\cos \mathrm{\ω}\cos i\cos \mathrm{\Ω}\\ \cos \mathrm{\ω}\sin i\end{array}\right]$

R＝r(P cos f+Q sin f)

$\underset{\‾}{V}=\sqrt{\frac{\mathrm{\μ}}{a(1-e^2)}}[-\underset{\‾}{P}\sin f+\underset{\‾}{Q}(e+\cos f)]$

其中地球赤道平均半径R_e＝6378.14km，地球引力场二阶带谐调和项系数J₂＝0.00108264。

通过数学仿真对本发明的方法进行验证，用卫星工具箱STK中的生成基准轨道，具体的仿真参数如下所示：

轨道积分：HPOP

引力模型：JGM-3；

太阳光压：Cr＝1.0；面质比：0.02m²/kg；阴影区模型：Dual Cone；

气动阻力：Cd＝2.2；面质比：0.02m²/kg；大气密度模型：Jachia-Roberts；

转动惯量：I_xx＝4500kg·m²，I_yy＝4500kg·m²，I_zz＝4500kg·m²；

三体引力：Sun，Earth；

输出的WGS84坐标系中卫星信息作为GPS测量数据，采样间隔为1min，滤波时间5天，对卫星的轨道进行预报，预报时间1～7天，与生成的TOD坐标系下轨道数据进行对比，对于不同的轨道高度，预报精度如附图2-5所示。仿真结果表明，对于450km以上的近地圆轨道，预报精度可达20km以内，具有很好的实用性。

本发明具有自启动(自初始化)、收敛性好、对测量数据的采样要求不严格(允许测量数据以非均匀间隔时间给出，采样时间可在几秒～10多分钟内变化)的优点，甚至允许在轨道的某一段时间内没有测量数据情况下(例如GPS天线被地球遮挡，或接收机出现临时性故障)，滤波器仍然能正常运行。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (7)

1.一种基于GPS测量数据的卫星在轨轨道预报方法，其特征在于：所述方法包括以下步骤：

S1：根据初始化标志判断当前计算状态，是否需要进行初始化，若需要进行初始化，则执行S2，否则执行S4；

S2：根据初始化方法标志进行滤波器初始化，其中初始化方法分为三种：a)根据GPS测量数据进行滤波器初始化，即根据GPS提供的测量数据，计算所需测量时刻对应的轨道平根数集的初始猜测；b)根据地面上注的星历模型参数进行初始化；c)用上一次滤波器计算收敛得到的星历参数进行初始化；

S3：初始化完成后，将初始化标志置位，回到S1；

S4：进行Kalman滤波计算，对星历模型参数进行滤波，得到更新的星历参数。

2.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述Kalman滤波计算为Bierman-UD分解形式的Kalman滤波。

3.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述步骤S1中的初始化标志为外部给出的是否进行滤波器初始化的标志。

4.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述步骤2中的初始化方法标志预先设置，在轨道的某一段时间内没有GPS测量数据情况下，初始化方法标志设置为a。

5.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述星历模型参数为星历模型参数对应的星上时间t_k、平均轨道半长轴平均轨道倾角平均升交点经度偏心率向量长周期变化的幅值偏心率向量长周期变化的相位平均纬度角升交点赤经的一次长期项摄动系数Ω₁、近地点幅角的一次长期项摄动系数ω₁、平均纬度角一次长期摄动系数λ₁、平均半长轴一次长期摄动系数a_d、偏心率向量x分量ξ的一次长期摄动系数ξ_d、偏心率向量y分量η的一次长期摄动系数η_d、偏心率向量固定偏移e_f、平均纬度角的二次长期摄动系数λ₂。

6.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述步骤S4后还包括步骤S5:根据更新的星历参数基于无奇点变量的拟平均根数法计算出任意时刻卫星的瞬时轨道根数，进而得到卫星在瞬时真赤道TOD坐标系中的位置和速度。

7.根据权利要求1所述的卫星在轨轨道预报方法，其特征在于：所述方法对测量数据的采样要求不严格，允许测量数据以非均匀间隔时间给出。