CN103453906B - 卫星轨道的预测方法 - Google Patents

Abstract

本发明提供了一种卫星轨道的预测方法，包括以下步骤：获得过去N小时的卫星事后精密轨道，获得卫星初始状态；计算卫星受力参数；根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。本发明的技术可用于预测卫星轨道，满足导航定位类卫星的基于预测的方法实现实时定轨需求。

Description

卫星轨道的预测方法

技术领域

本发明涉及一种卫星轨道预测方法，尤其是一种适用于低轨卫星的精密轨道预测方法，属于航天轨道动力学技术领域。

背景技术

近年来，低轨卫星导航正引起越来越多地关注。与现有的全球导航卫星相比，低轨卫星导航具有三大主要优势：相同发射功率下，较低的路径损耗给地面用户带来高于20dB的接收功率增益；卫星视角变化快，有利于快速精确的解算载波相位整周模糊度；适应导航能力泛在化发展趋势。

低轨卫星导航的一个基本问题，是实时定轨问题。在低轨卫星实时定轨领域，最常见的手段是在低轨卫星上安装全球导航卫星***(GlobalNavigationSatelliteSystem,GNSS)接收机。全球导航卫星***主要包括中国的北斗(BeiDou-2)、美国的全球定位***(GlobalPositioningSystem,GPS)、俄罗斯的全球导航卫星***(GlobalNavigationSatelliteSystem,GLONASS)和欧洲的伽利略***(Galileo)。全球导航卫星***接收机工作的基本原理是：接收到导航卫星发送的无线电信号并提取伪距，并根据4个以上伪距计算自身在地理坐标系中的位置，常见的解算算法有最小二乘法和卡尔曼滤波法。

然而，这种方法的缺点是：全球导航卫星***自身的实时轨道确定存在误差，时钟也存在误差，这导致低轨卫星星载接收机的实时定位精度大于6米，而低轨导航增强通常要求低轨卫星的实时轨道确定精度达到甚至优于1米，从而导致卫星轨道的预测不够精确。

发明内容

综上所述，确有必要提供一种具有较高预测精度的卫星精密轨道的预测方法。

一种卫星轨道的预测方法，主要包括以下步骤：步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。

本发明提供的卫星精密轨道的预测方法，利用事后精密轨道数据来预测未来轨道，综合利用了轨道动力学、事后精密轨道数据和低轨卫星轨道根数特征，其预测位置精度能达到米级。

附图说明

图1是本发明第一实施例提供的低轨卫星导航***示意图。

图2是本发明提供的卫星轨道预测方法流程图。

图3本发明提供的卫星轨道预测方法应用于GRACE卫星的轨道预测位置误差。

图4本发明提供的卫星轨道预测方法应用于GRACE卫星的轨道预测速度误差。

图5是本发明第二实施例提供的卫星轨道预测方法。

具体实施方式

下面根据说明书附图并结合具体实施例对本发明的技术方案进一步详细表述。

请参阅图1，S1～S5为传统导航卫星，它们为中轨导航卫星，工作轨道高度一般为2万公里。现有一个低轨导航星L1，工作轨道不高于3000公里，能够发射导航信号。一用户终端U1通过接收低轨导航星L1和传统导航卫星S1～S5的信号实现定位。在这个***中，为了提高导航精度，要求低轨导航星L1的精度越高越好。

请一并参阅图2，图2为本发明提供的卫星轨道预测方法流程图，主要包括以下步骤：

步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；

步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；

步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；

步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；

步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；

步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；

步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及

步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。

在步骤S10中，所述“事后精密轨道”，指的是过去N小时的精密卫星轨道。所述卫星事后精密轨道可通过常规技术获得，典型精度为1～10厘米。具体的，在卫星过顶时，将过去N小时星上的全球导航卫星***接收机的初始观测量下载到地面，使用地面根据全球导航卫星***的轨道、钟差和电离层修正数据产品经过约化动力学等方法进行计算，得到事后精密轨道。

在步骤S20中，根据N小时的卫星事后精密轨道计算卫星受力参数的方法，可进一步包括以下步骤：

步骤S21，将卫星所受力假想为保守力，并使用重力模型表示该保守力；

步骤S22，根据N小时的卫星事后精密轨道计算重力模型参数修正量；

步骤S23，利用所述重力模型参数修正量修正重力模型，所获得的重力模型即所求卫星受力参数。

在步骤S22中，所述重力模型参数修正量包括余弦规范化球面谐波系数C_nm的修正量ΔC_nm以及正弦规范化球面谐波系数S_nm的修正量ΔS_nm。其中，C_nm和S_nm来源于重力场模型，在地心大地坐标系中的具体表达式为：

$U=\frac{\mathrm{\μ}}{r}+\frac{\mathrm{\μ}}{r}\underset{n=1}{\overset{\mathrm{\∞}}{\Σ}}{\underset{m=0}{\overset{n}{\Σ}}\left(\frac{R_e}{r}\right)}^n{P}_{nm}\left(\sin \mathrm{\φ}\right)\left(C_{nm}\cos \right(m\mathrm{\λ})+S_{nm}\sin (m\mathrm{\λ}\left)\right)=\frac{\mathrm{\μ}}{r}+U_{per};$

其中，U_per为地球摄动势能；μ为地球引力常数；R_e为地球平均赤道半径；r为卫星事后精密轨道径矢；P_nm(sinφ)为sinφ的n阶m次伴随勒让德多项式；φ代表地心大地坐标系中的纬度，λ代表地心大地坐标系中的经度。在一般应用中，可使用重力模型获得规范化球面谐波系数的值C_nm、S_nm，如采用GRACE重力模型(GRACEGravityModel03,GGM03)参数。所述C_nm和S_nm的初始化方法可采用GGM03重力模型中的C_nm和S_nm作为卫星受力的初始参数，但并不限于所述方法。所述重力模型参数修正量为实际重力场规范化球面谐波系数与已知的规范化球面谐波系数的差值，包括前述余弦规范化球面谐波系数修正量ΔC_nm和正弦规范化球面谐波系数修正量ΔS_nm。

所述重力模型参数修正量ΔC_nm和ΔS_nm的计算可包括以下步骤：

步骤S221，将过去N小时中根据事后精密定轨结果所获得的任意时刻轨道状态定义为过去N小时中任意时刻的真实轨道状态；

步骤S222，根据所述真实轨道状态进行重力场反演，得到重力模型参数修正量ΔC_nm和ΔS_nm。

在步骤S222中，重力场反演的具体方法可包括如下步骤：

将过去N小时的真实轨道状态在动力学积分轨道处进行下述泰勒展开：

$x\left(t\right)=\hat{x}\left(t\right)+\mathrm{\Φ}(t,t_0)\mathrm{\Δ}x\left(t_0\right)+S\left(t\right)\mathrm{\Δ}h+\mathrm{\ϵ}$

其中，＝给定初始状态及初始重力场模型里预测的轨道状态；

Δx(t₀)＝初始轨道状态误差；

为状态向量对初始状态的偏导数，即状态转移矩阵；

为状态对重力场球面谐波系数的偏导数，即参数敏感矩阵。

其中，h即为待求重力场球面谐波系数构成的列向量，元素个数为nh个，C_nm,S_nm的排列顺序不限，举例如下：

h＝[C₁₀,C₁₁,C₂₀,C₂₁,C₂₂,...C_nn,S₁₁,S₂₁,S₂₂,...,S_nn]^T。

h的初值h₀可以采用GRACE重力模型(GRACEGravityModel03,GGM03)参数。

所述状态转移矩阵、参数敏感矩阵满足如下关系

$\left\{\begin{array}{c}\frac{d\mathrm{\Φ}(t,t_0)}{dt}=\frac{\∂}{\∂x\left(t_0\right)}\frac{dx\left(t\right)}{dt}=\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂x\left(t\right)}\frac{\∂x\left(t\right)}{\∂x\left(t_0\right)}=\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂x\left(t\right)}\mathrm{\Φ}(t,t_0)\\ \frac{dS\left(t\right)}{dt}=\frac{\∂}{\∂h}\frac{dx\left(t\right)}{dt}=\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂h}+\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂x\left(t\right)}S\left(t\right)\\ \mathrm{\Φ}(t_0,t_0)=I_{6\×6},S\left(t_0\right)=0_{6\×nh}\end{array}\right.$

其中， $\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂x\left(t\right)}=\left(\begin{array}{cc}{0}_{3\×3}& I_{3\×3}\\ \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂r\left(t\right)}& \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂\stackrel{\·}{r}\left(t\right)}\end{array}\right),\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂h}=\left(\begin{array}{c}{0}_{3\×n_h}\\ \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂h}\end{array}\right).$

若令 $\left(\begin{array}{cc}\mathrm{\Φ}\left(t,t_0\right)& S\left(t\right)\end{array}\right)=\left(\frac{\∂x\left(t\right)}{\∂\left(x\right(t_0),h)}\right)=\left(\begin{array}{cc}\frac{\∂r\left(t\right)}{\∂\left(r\right(t_0),\stackrel{\·}{r}(t_0\left)\right)}& \frac{\∂r\left(t\right)}{\∂h}\\ \frac{\∂\stackrel{\·}{r}\left(t\right)}{\∂\left(r\right(t_0),\stackrel{\·}{r}(t_0\left)\right)}& \frac{\∂\stackrel{\·}{r}\left(t\right)}{\∂h}\end{array}\right)=\left(\begin{array}{c}\mathrm{\Ψ}\left(t\right)\\ \stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right),$ 则状态转移矩阵和参数敏感矩阵的联合微分可写为：

$\begin{array}{c}\frac{s\left(\begin{array}{cc}\mathrm{\Φ}& S\end{array}\right)}{dt}=\frac{d}{dt}\frac{\∂x\left(t\right)}{\∂\left(x\right(t_0),h)}=\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂x\left(t\right)}\frac{\∂x\left(t\right)}{\∂\left(x\right(t_0),h)}+\frac{\∂\stackrel{\·}{x}\left(t\right)}{\∂\left(x\right(t_0),h)}\\ =\left(\begin{array}{cc}{0}_{3\×3}& I_{3\×3}\\ \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂r\left(t\right)}& \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂\stackrel{\·}{r}\left(t\right)}\end{array}\right)\frac{\∂x\left(t\right)}{\∂\left(x\right(t_0),h)}+\left(\begin{array}{cc}{0}_{3\×6}& 0_{3\×nh}\\ 0_{3\×6}& \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂h}\end{array}\right)\\ =\left(\begin{array}{c}\stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\\ \stackrel{\·\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right)\end{array}---\left(1\right)$

考虑到 $\mathrm{\Ψ}\left(t\right)=\frac{\∂r\left(t\right)}{\∂\left(r\right(t_0),\stackrel{\·}{r}(t_0),h)}$ 中时间参数t连续且与x(t₀)和h相关性小或者无关，可以交换微分的先后次序：

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\Ψ}}=\frac{d\mathrm{\Ψ}}{dt}=\frac{\∂}{\∂\left(r\right(t_0),r(t_0),h)}\frac{dr\left(t\right)}{dt}\\ \stackrel{\·\·}{\mathrm{\Ψ}}=\frac{d^2\mathrm{\Ψ}}{{\mathrm{dt}}^2}\frac{\∂}{\∂\left(r\right(t_0),r(t_0),h)}\frac{d^{2}r\left(t\right)}{{\mathrm{dt}}^2}\end{array}\right.$

令 $A\left(t\right)=\frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂r\left(t\right)},B\left(t\right)=\frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂\stackrel{\·}{r}\left(t\right)},C\left(t\right)=\left(\begin{array}{cc}{0}_{3\×6}& \frac{\∂\stackrel{\·\·}{r}\left(t\right)}{\∂h}\end{array}\right),$ 上述变分方程可写为

$\left\{\begin{array}{c}\stackrel{\·\·}{\mathrm{\Ψ}}\left(t\right)=A\left(t\right)\mathrm{\Ψ}\left(t\right)+B\left(t\right)\stackrel{\·}{\mathrm{\Ψ}}\left(t\right)+C\left(t\right)\\ \left(\begin{array}{c}\mathrm{\Ψ}\left(t\right)\\ \stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{I}_{6\×6}& 0_{6\×n_h}\end{array}\right)\end{array}\right.---\left(2\right)$

其中n_h＝(N+1)²(N为重力场球谐系数的阶数)为球谐系数个数，由于卫星加速度由受力决定，可令B(t)＝0。

将式 $x\left(t\right)=\hat{x}\left(t\right)+\mathrm{\Φ}(t,t_0)\mathrm{\Δ}x\left(t_0\right)+S\left(t\right)\mathrm{\Δ}h+\mathrm{\ϵ}$ 写成矩阵形式为

$\mathrm{\Δ}x\left(t\right)=\left(\begin{array}{c}\mathrm{\Ψ}\left(t\right)\\ \stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right)\left(\begin{array}{c}\mathrm{\Δ}x\left(t_0\right)\\ \mathrm{\Δ}h\end{array}\right)+\mathrm{\ϵ}$

计算出Ψ(t)和后，即根据所述反演变分方程形式(2)并利用最小二乘法解算地球重力场球面谐波系数修正量Dh＝[DC₁₀,DC₁₁,DC₂₀,DC₂₁,DC₂₂,...DC_nn,DS₁₁,DS₂₁,DS₂₂,...,DS_nn]^T。

所述反演变分方程形式(2)可通过微分法求解。

在步骤S23中，获得了Dh后，就可以获得反演后的h＝Dh+h₀,即修正后的重力模型参数，即所述的卫星受力参数。

在步骤S30中，选择一个初始时间和历史精密轨道得到初始卫星状态，根据卫星初始状态和卫星受力参数来计算未来M小时卫星轨道根数(即第一预测轨道根数)。

所述第一预测轨道根数包括第一预测轨道轨道倾角i₁(t)和第一预测轨道升交点赤经Ω₁(t)等的计算方法可通过对下式进行积分获得：

$\frac{da}{dt}=\frac{2esin\mathrm{\θ}}{n\mathrm{\χ}}{F}_r+\frac{2a\mathrm{\χ}}{nr}{F}_t$

$\frac{de}{dt}=\frac{\mathrm{\χ}sin\mathrm{\θ}}{na}{F}_r+\frac{\mathrm{\χ}}{{\mathrm{na}}^{2}e}(\frac{a^2{\mathrm{\χ}}^2}{r}-r)F_t$

$\frac{dM}{dt}=n-\frac{1}{na}(\frac{2r}{a}-\frac{{\mathrm{\χ}}^2}{e}cos\mathrm{\θ})F_r-\frac{{\mathrm{\χ}}^2}{nae}(1+\frac{r}{{\mathrm{a\χ}}^2}){\mathrm{sin\θF}}_t$

$\frac{di}{dt}=\frac{rcos\mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}}{F}_w$

$\frac{d\mathrm{\Ω}}{dt}=\frac{rsin\mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}\sin i}{F}_w$

$\frac{d\mathrm{\ω}}{dt}=\frac{\mathrm{\χ}\cos \mathrm{\θ}}{nae}{F}_r+\frac{p}{eh}\left(\sin \mathrm{\θ}\right(1+\frac{1}{1+e\cos \mathrm{\θ}}\left)\right)F_t-\frac{r\cot i\sin \mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}}{F}_w$

其中，F_r,F_t,F_w＝摄动加速度分量，方向分别为径向e_r＝r/r背向地心为正，横向e_t＝e_w×e_r速度增加方向为正，轨道面法线方向e_w＝r×v/|r×v|；

θ＝真近点角；

n＝平均角速度；

$\mathrm{\χ}=\sqrt{1-e^2};$

Φ＝θ+ω，升交点角距；

p＝a(1-e²)，其中a为卫星轨道半长轴，e为卫星离心率；

$h=\sqrt{\mathrm{\μ}p};$

i为卫星轨道倾角；Ω为卫星轨道升交点赤经；M为卫星轨道平近点角；ω为卫星轨道近地点角距；μ为地心引力常数。

摄动源为非球形引力摄动，是为保守力，则摄动加速度为

$F_r=\frac{\∂U_{per}}{\∂r},F_t=\frac{\∂U_{per}}{r\∂\mathrm{\θ}},F_w=\frac{\∂U_{per}}{rsin\mathrm{\Φ}\∂i}$

至此，我们获得了未来M小时卫星第一预测轨道根数，包括第一预测轨道根数中的轨道倾角i₁(t)和第一预测轨道根数升交点赤经Ω₁(t)等。

在步骤S40中，将卫星所有受力归结为重力进行轨道预测仍然存在误差，主要是轨道倾角误差和升交点赤经误差。因此，我们需要进一步根据所述过去N小时的卫星事后精密轨道计算轨道根数修正量，并根据所述轨道根数修正系数修正轨道根数，获得修正后的修正轨道根数。具体地，轨道根数修正量包括轨道倾角修正量和升交点赤经修正量。

所述轨道倾角修正量可通过以下方法计算：

步骤S411，计算过去N小时内事后精密轨道状态对应的轨道根数的轨道倾角i(t)；

步骤S412，将过去N小时的i(t)与过去N小时第一预测轨道根数的轨道倾角相减，获得i'(t)；

步骤S413，用多项式拟合i'(t)的长期缓慢变化：其中，I的具体值可以根据精度要求给定，是一个设计参数，t₀是过去N小时的起始时刻，α、w和γ为Ω’(t)拟合多项式的系数；

步骤S414，根据将过去N小时的i'(t)的值使用最小二乘方法计算上式中的a_i；以及

步骤S415，根据a_i计算未来M小时的轨道倾角修正量。

所述升交点赤经修正量可通过以下方法计算：

步骤S421，计算过去N小时内事后精密轨道状态对应的轨道根数的升交点赤经Ω(t)；

步骤S422，将过去N小时的Ω(t)与过去N小时第一预测轨道根数的升交点赤经相减，获得Ω'(t)；

步骤S423，用多项式拟合Ω(t)的长期缓慢变化：其中，J的具体值可以根据精度要求给定，是一个设计参数，t₀是过去N小时的起始时刻；

步骤S424，根据过去N小时的Ω'(t)的值使用最小二乘方法计算上式中的α,w,γ；

步骤S425，根据α,w,γ计算未来M小时的升交点赤经修正量。

在步骤S50中，所述修正后的修正轨道根数包括修正后的轨道倾角和修正后的升交点赤经，其中所述修正后的轨道倾角为i₂(t)＝i₁(t)+i'(t)，修正后的轨道根数中的升交点赤经为Ω₂(t)＝Ω₁(t)+Ω'(t)。

在步骤S60中，所述轨道状态包括卫星在地心坐标系下的三维位置和速度，所述第一预测轨道状态就可以通过上述修正后的轨道倾角和修正后的升交点赤经计算。

在步骤S70中，初始状态误差也是影响轨道预测的重要误差项。因此，我们需要进一步根据所述过去N小时的卫星事后精密轨道计算初始状态修正量，并根据初始状态修正量修正第一预测轨道状态，获得修正轨道状态即最终轨道状态。

所述初始误差修正量的计算方法如下：

根据位置与加速度的关系，由卫星初始状态造成的位置误差可以表示为δ(x)＝Ax+B，其中，x为任意一轴的卫星实际轨道位置分量。A、B由卫星初始状态误差决定，这说明在初始状态误差一定的情况下，某时刻有初始状态引起的位置预测误差与当前的位置为线性关系。若知道K(K大于等于2)个时刻对应的位置分量x对应预测位置误差δ(x)，则上式可写为

δ(x)＝Ax+B

即可通过最小二乘法计算出A和B，进而可以由上式进而根据A、B计算未来M小时由初始状态引起的位置误差δ(x)，即所述初始误差修正量。

步骤S80，将初始误差修正量加到第一预测轨道状态，就获得了修正轨道状态即最终轨道状态。

根据本发明的方法，可以利用前一天精密轨道数据，预测下一天的低轨卫星轨道。作为具体的应用，我们针对GRACE卫星的数据进行理论计算，I＝J＝6，得到的结果如图3和图4所示。其中图3的横坐标单位为分钟，纵坐标单位为米。可以看到，实际单轴24小时的预测误差最大约为1米。图4是速度预测误差，其中横坐标单位为分钟，纵坐标为单位为米/秒，可以看到，单轴速度预测速度误差最大约为0.06米/秒。

进一步，请参阅图5，图5为本发明第二实施例提供的卫星轨道预测方法的示意图，图5的情况出现在本发明的方法首次运行时。不失一般性，我们假设N＝M。在时刻1，我们可以使用第一个N小时的事后精密轨道和重力模型做轨道状态预测；由于不存在第一个N小时的预测轨道，因此不能做根据本发明的修正计算，包括轨道根数修正和初始轨道状态。

根据本发明的轨道预测方法进行实时定轨的方法，包括以下步骤：

(1)在卫星过顶时下传过去N小时的全球导航卫星接收机观测量；

(2)根据本发明的轨道预测方法预测未来M小时的轨道状态,并将所预测的未来M小时的轨道状态在卫星过顶同时上传的卫星；

(3)在未来M小时使用根据当前时间所上传的预测的轨道状态，通过内插计算卫星的实时轨道位置和速度。

在实际应用中，卫星的平均过顶间隔约为12小时，考虑事后精密定轨的延迟，也可以考虑将时间分为若干12小时，即M＝N＝12小时，使用根据过去[2N,N]小时的事后精密轨道和本发明的轨道预测方法预测未来M小时的轨道状态(也可以是轨道根数)也能得到较好的预测结果，计算过程可以发生在两次卫星过顶之间。因此，本发明的上述实施例，是为了说明而不是限制本发明的方法。

本发明提供的卫星精密轨道的预测方法，利用事后精密轨道数据来预测未来轨道，综合利用了轨道动力学、事后精密轨道数据和低轨卫星轨道根数特征，其预测位置精度能达到米级，实现了24小时预测精度优于1米。本发明的技术可用于预测卫星轨道，满足导航定位类卫星的基于预测的方法实现实时定轨需求，从而能够解决低轨卫星的实时轨道计算问题。另外，所述卫星精密轨道的预测方法不仅可以用于低轨导航星，也可以用于其它需要实时精密确定轨道的卫星。

另外，本领域技术人员还可在本发明精神内作其它变化，当然这些依据本发明精神所作的变化，都应包含在本发明所要求保护的范围内。

Claims (13)

1.一种卫星轨道的预测方法，主要包括以下步骤：

步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；

步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；

步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；

步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；

步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；

步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；

步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及

步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。

2.如权利要求1所述的卫星轨道的预测方法，其特征在于，所述卫星受力参数的计算方法包括以下步骤：

将卫星所受力假想为保守力，并使用重力模型表示该保守力；

根据N小时的卫星事后精密轨道计算重力模型参数修正量；

根据所述重力模型参数修正量修正重力模型，修正后的重力模型即所述卫星受力参数。

3.如权利要求2所述的卫星轨道的预测方法，其特征在于，所述重力模型参数修正量包括余弦规范化球面谐波系数C_nm的修正量ΔC_nm以及正弦规范化球面谐波系数S_nm的修正量ΔS_nm。

4.如权利要求3为所述的卫星轨道的预测方法，其特征在于，C_nm和S_nm来源于重力场模型，在地心大地坐标系中的具体表达式为：

$U=\frac{\mathrm{\μ}}{r}+\frac{\mathrm{\μ}}{r}\underset{n=1}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=0}{\overset{n}{\Σ}}{\left(\frac{R_e}{r}\right)}^n{P}_{nm}\left(sin\mathrm{\φ}\right)\left(C_{nm}cos\right(m\mathrm{\λ})+S_{nm}sin(m\mathrm{\λ}\left)\right)=\frac{\mathrm{\μ}}{r}+U_{per};$

其中，U_per为地球摄动势能；μ为地球引力常数；R_e为地球平均赤道半径；r为卫星事后精密轨道径矢；P_nm(sinφ)为sinφ的n阶m次伴随勒让德多项式；φ代表地心大地坐标系中的纬度，λ代表地心大地坐标系中的经度。

5.如权利要求4所述的卫星轨道的预测方法，其特征在于，所述重力模型参数修正量ΔC_nm和ΔS_nm的计算包括以下步骤：

将过去N小时中根据事后精密定轨结果所获得的任意时刻轨道状态定义为过去N小时中任意时刻的真实轨道状态；

根据所述真实轨道状态进行重力场反演，得到重力模型参数修正量ΔC_nm和ΔS_nm。

6.如权利要求5所述的卫星轨道的预测方法，其特征在于，所述重力场反演包括如下步骤：

将过去N小时的真实轨道状态在动力学积分轨道处进行下述泰勒展开：

$x\left(t\right)=\hat{x}\left(t\right)+\mathrm{\Φ}(t,t_0)\mathrm{\Δ}x\left(t_0\right)+S\left(t\right)\mathrm{\Δ}h+\mathrm{\ϵ},$

其中，

Δx(t₀)＝初始轨道状态误差；

为状态向量对初始状态的偏导数，即状态转移矩阵；

为状态对重力场球面谐波系数的偏导数，即参数敏感矩阵；

其中，h即为重力场球面谐波系数Cnm,Snm构成的列向量，元素个数为nh个，h的初值h0采用GRACE重力模型参数；

将式 $x\left(t\right)=\hat{x}\left(t\right)+\mathrm{\Φ}(t,t_0)\mathrm{\Δ}x\left(t_0\right)+S\left(t\right)\mathrm{\Δ}h+\mathrm{\ϵ}$ 写成矩阵形式：

$\mathrm{\Δ}x\left(t\right)=\left(\begin{array}{c}\mathrm{\Ψ}\left(t\right)\\ \stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right)\left(\begin{array}{c}\mathrm{\Δ}x\left(t_0\right)\\ \mathrm{\Δ}h\end{array}\right)+\mathrm{\ϵ};$

计算Ψ(t)和根据反演变分方程式：

$\{\begin{array}{c}\stackrel{\·\·}{\mathrm{\Ψ}}\left(t\right)=A\left(t\right)\mathrm{\Ψ}\left(t\right)+B\left(t\right)\stackrel{\·}{\mathrm{\Ψ}}\left(t\right)+C\left(t\right)\\ \left(\begin{array}{c}\mathrm{\Ψ}\left(t\right)\\ \stackrel{\·}{\mathrm{\Ψ}}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{I}_{6\×6}& 0_{6\×n_h}\end{array}\right)\end{array};$

并利用最小二乘法解算地球重力场球面谐波系数修正量ΔC_nm和ΔS_nm，其中n_h＝(N+1)²(N为重力场球谐系数的阶数)为球谐系数个数。

7.如权利要求6所述的卫星轨道的预测方法，其特征在于，所述第一预测轨道根数包括第一预测轨道轨道倾角i1(t)和第一预测轨道升交点赤经Ω₁(t)，所述第一预测轨道轨道倾角i₁(t)和第一预测轨道升交点赤经Ω₁(t)的计算方法通过对下式进行积分获得：

$\frac{da}{dt}=\frac{2esin\mathrm{\θ}}{n\mathrm{\χ}}{F}_r+\frac{2a\mathrm{\χ}}{nr}{F}_t$

$\frac{de}{dt}=\frac{\mathrm{\χ}sin\mathrm{\θ}}{na}{F}_r+\frac{\mathrm{\χ}}{{\mathrm{na}}^{2}e}(\frac{a^2{\mathrm{\χ}}^2}{r}-r)F_t$

$\frac{dM}{dt}=n-\frac{1}{na}(\frac{2r}{a}-\frac{{\mathrm{\χ}}^2}{e}cos\mathrm{\θ})F_r-\frac{{\mathrm{\χ}}^2}{nae}(1+\frac{r}{{\mathrm{a\χ}}^2}){\mathrm{sin\θF}}_t$

$\frac{di}{dt}=\frac{rcos\mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}}{F}_w$

$\frac{d\mathrm{\Ω}}{dt}=\frac{rsin\mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}\sin i}{F}_w$

$\frac{d\mathrm{\ω}}{dt}=-\frac{\mathrm{\χ}\cos \mathrm{\θ}}{nae}{F}_r+\frac{p}{eh}\left(\sin \mathrm{\θ}\right(1+\frac{1}{1+e\cos \mathrm{\θ}}\left)\right)F_t-\frac{r\cot isin\mathrm{\Φ}}{{\mathrm{na}}^2\mathrm{\χ}}{F}_w;$

其中，F_r,F_t,F_w＝摄动加速度分量，方向分别为径向e_r＝r/r背向地心为正，横向e_t＝e_w×e_r速度增加方向为正，轨道面法线方向e_w＝r×v/|r×v|；

θ＝真近点角；

n＝平均角速度；

$\mathrm{\χ}=\sqrt{1-e^2};$

Φ＝θ+ω，升交点角距；

p＝a(1-e²)，其中a为卫星轨道半长轴，e为卫星离心率；

$h=\sqrt{\mathrm{\μ}p};$

其中，i为卫星轨道倾角；Ω为卫星轨道升交点赤经；M为卫星轨道平近点角；ω为卫星轨道近地点角距；μ为地心引力常数。

8.如权利要求7所述的卫星轨道的预测方法，其特征在于，摄动源为非球形引力摄动，是为保守力，摄动加速度为：

$\begin{array}{ccc}{F}_r=\frac{\∂U_{per}}{\∂r},& F_t=\frac{\∂U_{per}}{r\∂\mathrm{\θ}},& F_w=\frac{\∂U_{per}}{r\sin \mathrm{\Φ}\∂i}\end{array}.$

9.如权利要求8所述的卫星轨道的预测方法，其特征在于，所述轨道根数修正量包括轨道倾角修正量和升交点赤经修正量。

10.如权利要求9所述的卫星轨道的预测方法，其特征在于，所述轨道倾角修正量的计算包括以下步骤：

计算过去N小时内事后精密轨道状态对应的轨道根数的轨道倾角i(t)；

将过去N小时的i(t)与过去N小时第一预测轨道根数的轨道倾角相减，获得i'(t)；

用多项式拟合i'(t)的长期缓慢变化：其中，I的具体值根据精度要求给定，是一个设计参数，t₀是过去N小时的起始时刻，a_i为i'(t)拟合多项式的第i阶系数；

根据将过去N小时的i'(t)的值使用最小二乘方法计算上式中的a_i；以及根据a_i计算未来M小时的轨道倾角修正量。

11.如权利要求10所述的卫星轨道的预测方法，其特征在于，所述升交点赤经修正量的计算包括以下步骤：

计算过去N小时内事后精密轨道状态对应的轨道根数的升交点赤经Ω(t)；

将过去N小时的Ω(t)与过去N小时第一预测轨道根数的升交点赤经相减，获得Ω'(t)；

用多项式拟合Ω(t)的长期缓慢变化：其中，J的具体值根据精度要求给定，是一个设计参数，t₀是过去N小时的起始时刻，α、w和γ为Ω’(t)拟合多项式的系数；

根据过去N小时的Ω'(t)的值使用最小二乘方法计算上式中的α,w,γ；

根据α,w,γ计算未来M小时的升交点赤经修正量。

12.如权利要求11所述的卫星轨道的预测方法，其特征在于，所述初始状态误差修正量由以下方法计算：由卫星初始状态造成的位置误差表示为δ(x)＝Ax+B，其中x为任意一轴的卫星实际轨道位置分量；A、B由卫星初始状态误差决定；通过最小二乘法计算出A和B，进而上式计算未来M小时由初始状态引起的位置误差δ(x)，即所述初始状态误差修正量。

13.如权利要求1所述的卫星轨道的预测方法，其特征在于，进一步包括以下步骤：

在卫星过顶时下传过去N小时的全球导航卫星接收机观测量；

预测卫星未来M小时的轨道状态,并将所预测的未来M小时的轨道状态在卫星过顶同时上传；

在未来M小时使用根据当前时间所上传的预测的轨道状态，通过内插计算卫星的实时轨道位置和速度。