一种大偏心率轨道的星上快速高精度外推方法

摘要

一种大偏心率轨道的星上快速高精度外推方法，火星探测器任务轨道为大偏心率(e＞0.6)椭圆轨道，采用数值积分方法，要做到较高精度的轨道计算，需要考虑较高阶次的火星形状摄动以及其他各种摄动因素影响，从而导致计算量较大，若星上计算机复位或切机，数值积分方法将会中断；解析法只适用于偏心率较小(e＜0.6)的椭圆轨道；利用精密轨道提供的探测器位置拟合切比雪夫多项式系数，将会带来较大的拟合误差。本发明通过引入精密轨道与二体轨道的位置差，利用位置差拟合切比雪夫多项式系数，可以很好地解决大偏心率椭圆轨道的星上轨道计算问题，计算量小，且精度高。

说明书

技术领域

本发明涉及一种大偏心率轨道的星上方法，属于卫星轨道技术领域。

背景技术

航天器在空间的运动轨迹称为航天器的轨道，通过高精度的外推计算可以 得到航天器在过去、当前和未来一段时间内任一时刻的运动状态。火星探测器 处于星际飞行，飞行距离最大达到4亿公里，相比地球卫星和月球探测器，火 星探测器面临星地距离远、时延大和长时间日凌等问题，火星探测器对GNC 的自主性提出了更高要求。因此，星上高精度轨道计算和轨道外推非常重要。 求解卫星轨道动力学方程主要有数值法和解析法两种方法。数值方法解算精度 高，但存在计算量大以及星上计算机复位或切机时数值积分将会中断的弊端； 由于火星探测器任务轨道为大偏心率(e＞0.6)椭圆轨道，解析方法也不适用。 目前有文献利用精密轨道提供的探测器位置拟合切比雪夫多项式系数，进而实 现圆轨道的星上轨道计算，但该方法用于大偏心率椭圆轨道时将会带来较大的 拟合误差。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提出一种大偏心率轨道 的星上快速高精度外推方法，解决原来星上轨道外推算法拟合误差大的问题， 从而提高星上轨道外推计算的精度。

本发明的技术解决方案是：一种大偏心率轨道的星上快速高精度外推方法， 所述大偏心率指e＞0.6的椭圆轨道，包括以下步骤：

(1)在地面求解大偏心率椭圆轨道的精密轨道与二体轨道的X，Y，Z三 轴位置差，之后利用切比雪夫多项式拟合法拟合得到切比雪夫多项式系数；所 述精密轨道是指卫星实际运行的轨道，由地面测轨提供，二体轨道是指卫星只 受目标天体中心引力作用下的理想轨道，由理论计算得到；所述拟合过程求解 如下：

(1.1)首先将时间t∈[t₀,t₀+Δt]变换成τ∈[-1,1]，变换公式为 $\tau =\frac{2}{\mathrm{\Delta t}}(t-t_0)-1,$t∈[t₀,t₀+Δt]                (1)

其中t₀为开始历元，Δt为拟合时间区间的长度，τ为标准化的时间变量；

(1.2)将大偏心率椭圆轨道的精密轨道与二体轨道的三轴位置差的X，Y， Z分量采用如下切比雪夫多项式表示

$\left(\begin{array}{c}X\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{X_i}{T}_i\left(\tau \right)\\ Y\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{Y_i}{T}_i\left(\tau \right)\\ Z\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{Z_i}{T}_i\left(\tau \right)\end{array}\right)---\left(2\right)$

三轴速度差分量表示为

$\left(\begin{array}{c}\stackrel{.}{X}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{X_i}{\stackrel{.}{T}}_i\left(\tau \right)\\ \stackrel{.}{Y}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{Y_i}{\stackrel{.}{T}}_i\left(\tau \right)\\ \stackrel{.}{Z}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{Z_i}{\stackrel{.}{T}}_i\left(\tau \right)\end{array}\right)---\left(3\right)$

式中，i＝0,1,2,…,n为切比雪夫多项式的阶数，分别为X坐标分量、 Y坐标分量、Z坐标分量切比雪夫多项式的系数，其中切比雪夫多项式T_i和分 别由以下递推公式确定

$\left(\begin{array}{c}{T}_0\left(\tau \right)=1\\ T_1\left(\tau \right)=\tau \\ ...\\ T_n\left(\tau \right)=2\tau T_{n-1}\left(\tau \right)-T_{n-2}\left(\tau \right)\end{array}\right)---\left(4\right)$

其中公式(4)中n≥2；

$\left(\begin{array}{c}{\stackrel{.}{T}}_1\left(\tau \right)=\frac{2}{\mathrm{\Delta t}}\\ {\stackrel{.}{T}}_2\left(\tau \right)=\frac{8}{\mathrm{\Delta t}}\tau \\ ...\\ {\stackrel{.}{T}}_n\left(\tau \right)=\frac{2n}{n-1}\tau {\stackrel{.}{T}}_{n-1}\left(\tau \right)-\frac{n}{n-2}{\stackrel{.}{T}}_{n-2}\left(\tau \right)\end{array}\right)---\left(5\right)$

其中公式(5)中n≥3；

(1.3)将卫星精密轨道和二体轨道的X轴位置差为观测值，则误差方程为

$V_{X_k}=\underset{i=0}{\overset{n}{\Sigma }}{C}_{X_i}{T}_i\left({\tau }_k\right)-X_k---\left(6\right)$其中k＝1,2,…,m为拟合时间区间内的观测值个数，X_k为对应k时刻的X轴位置 差，V_x为相应的误差向量；

误差方程的矩阵展开式为

$\left(\begin{array}{c}{V}_{X_1}\\ V_{X_2}\\ .\\ .\\ .\\ V_{X_m}\end{array}\right)=\left(\begin{array}{cccc}{T}_0\left({\tau }_1\right)& T_1\left({\tau }_1\right)& ...& T_n\left({\tau }_1\right)\\ T_0\left({\tau }_2\right)& T_1\left({\tau }_2\right)& ...& T_n\left({\tau }_2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ T_0\left({\tau }_m\right)& T_1\left({\tau }_m\right)& ...& T_n\left({\tau }_m\right)\end{array}\right)\left(\begin{array}{c}{C}_{X_0}\\ C_{X_1}\\ .\\ .\\ .\\ C_{X_m}\end{array}\right)-\left(\begin{array}{c}{X}_1\\ X_2\\ .\\ .\\ .\\ X_m\end{array}\right)---\left(7\right)$

令

$V_X={\left(\begin{array}{cccc}{V}_{X_1}& V_{X_2}& ...& V_{X_m}\end{array}\right)}^T,$f_x＝[X₁ X₂…X_m]^T，$M={\left(\begin{array}{cccc}{C}_{X_0}& C_{X_1}& ...& C_{X_m}\end{array}\right)}^T$

$B=\left(\begin{array}{cccc}{T}_0\left({\tau }_1\right)& T_1\left({\tau }_1\right)& ...& T_n\left({\tau }_1\right)\\ T_0\left({\tau }_2\right)& T_1\left({\tau }_2\right)& ...& T_n\left({\tau }_2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ T_0\left({\tau }_m\right)& T_1\left({\tau }_m\right)& ...& T_n\left({\tau }_m\right)\end{array}\right)$

则有向量表达式

V_x＝BM-f_x                    (8)

为了求得唯一解，采用最小二乘法，使所求的切比雪夫多项式系数阵M满 足V_x^TV_x＝min，其中V_x^T为V_x的转置矩阵，求解得到

M＝(B^TB)^-1B^Tf_x                        (9)

即得到X坐标分量切比雪夫多项式的拟合系数

同理，对Y、Z位置差重复步骤(1.1)～(1.3)的计算得到Y坐标分量和 Z坐标分量切比雪夫多项式的拟合系数知

(2)将大偏心率椭圆轨道初值和拟合得到的切比雪夫多项式系数和注入卫星；

(3)卫星利用轨道初值计算二体轨道的三轴位置和三轴速度，利用切比雪 夫多项式系数计算三轴位置差和三轴速度差，将其与对应时刻的二体轨道位置 和速度相加，得到卫星的实际位置和速度。

本发明与现有技术相比的优点在于：

(1)本发明引入精密轨道与二体轨道的位置差，利用位置差拟合切比雪夫 多项式系数，可以很好地解决大偏心率椭圆轨道的星上轨道计算问题，计算量 小，且精度高；

(2)本发明本用于星上快速轨道计算，适用于大偏心率的轨道，算法简单， 易操作，精度高，具有很强的实用性。

附图说明

图1为本发明的实现流程图；

图2为本发明实施例中利用精密轨道提供的探测器位置拟合切比雪夫多项 式系数，进行环火大椭圆轨道的拟合结果；左图代表拟合的位置误差，右图代 表拟合的速度误差；

图3为本发明实施例中利用精密轨道与二体轨道的位置差拟合切比雪夫多 项式系数，进行环火大椭圆轨道的拟合结果；左图代表拟合的位置误差，右图 代表拟合的速度误差。

具体实施方式

如图1所示，本发明具体实现包括以下步骤：

(1)在地面求解大偏心率椭圆轨道的精密轨道与二体轨道的X，Y，Z三 轴位置差，之后利用切比雪夫多项式拟合法拟合得到切比雪夫多项式系数；

切比雪夫多项式拟合就是根据给定的数据拟合出一个函数，使其在给定点 的函数值与给定值之间的方差和最小，且该函数是以切比雪夫多项式为基函数 构成的函数。

所述精密轨道是指卫星实际运行的轨道，由地面测轨提供，二体轨道是指 卫星只受目标天体中心引力作用下的理想轨道，由理论计算得到；二体轨道的 初值与精密轨道的初值一致，而平近点角M的取值从精密轨道初始平近点角M₀开始，以的步长增加，其中μ_m为火星的引力常数，a为精密轨道的 半长轴，n₀为轨道角速度，h为轨道角速度的整数倍。

所述拟合过程求解如下：

(1.1)首先将时间t∈[t₀,t₀+Δt]变换成τ∈[-1,1]，变换公式为

$\tau =\frac{2}{\mathrm{\Delta t}}(t-t_0)-1,$t∈[t₀,t₀+Δt]

其中t₀为开始历元，Δt为拟合时间区间的长度，τ为标准化的时间变量；

(1.2)将大偏心率椭圆轨道的精密轨道与二体轨道的三轴位置差的X，Y，

Z分量采用如下切比雪夫多项式表示

$\left(\begin{array}{c}X\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{X_i}{T}_i\left(\tau \right)\\ Y\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{Y_i}{T}_i\left(\tau \right)\\ Z\left(t\right)=\underset{i=0}{\overset{n}{\Sigma }}{C}_{Z_i}{T}_i\left(\tau \right)\end{array}\right)$

三轴速度差分量表示为

$\left(\begin{array}{c}\stackrel{.}{X}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{X_i}{\stackrel{.}{T}}_i\left(\tau \right)\\ \stackrel{.}{Y}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{Y_i}{\stackrel{.}{T}}_i\left(\tau \right)\\ \stackrel{.}{Z}\left(t\right)=\underset{i=1}{\overset{n}{\Sigma }}{C}_{Z_i}{\stackrel{.}{T}}_i\left(\tau \right)\end{array}\right)$

式中，i＝0,1,2,…,n为切比雪夫多项式的阶数，分别为X坐标分量、 Y坐标分量、Z坐标分量切比雪夫多项式的系数，其中切比雪夫多项式T_i和分 别由以下递推公式确定

$\left(\begin{array}{c}{T}_0\left(\tau \right)=1\\ T_1\left(\tau \right)=\tau \\ ...\\ T_n\left(\tau \right)=2\tau T_{n-1}\left(\tau \right)-T_{n-2}\left(\tau \right)\end{array}\right)$

其中n≥2；

$\left(\begin{array}{c}{\stackrel{.}{T}}_1\left(\tau \right)=\frac{2}{\mathrm{\Delta t}}\\ {\stackrel{.}{T}}_2\left(\tau \right)=\frac{8}{\mathrm{\Delta t}}\tau \\ ...\\ {\stackrel{.}{T}}_n\left(\tau \right)=\frac{2n}{n-1}\tau {\stackrel{.}{T}}_{n-1}\left(\tau \right)-\frac{n}{n-2}{\stackrel{.}{T}}_{n-2}\left(\tau \right)\end{array}\right)$

其中n≥3；

(1.3)将卫星精密轨道和二体轨道的X轴位置差为观测值，则误差方程为

$V_{X_k}=\underset{i=0}{\overset{n}{\Sigma }}{C}_{X_i}{T}_i\left({\tau }_k\right)-X_k$

其中k＝1,2,…,m为拟合时间区间内的观测值个数，X_k为对应k时刻的X轴位置 差，V_x为相应的误差向量；

误差方程的矩阵展开式为

$\left(\begin{array}{c}{V}_{X_1}\\ V_{X_2}\\ .\\ .\\ .\\ V_{X_m}\end{array}\right)=\left(\begin{array}{cccc}{T}_0\left({\tau }_1\right)& T_1\left({\tau }_1\right)& ...& T_n\left({\tau }_1\right)\\ T_0\left({\tau }_2\right)& T_1\left({\tau }_2\right)& ...& T_n\left({\tau }_2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ T_0\left({\tau }_m\right)& T_1\left({\tau }_m\right)& ...& T_n\left({\tau }_m\right)\end{array}\right)\left(\begin{array}{c}{C}_{X_0}\\ C_{X_1}\\ .\\ .\\ .\\ C_{X_m}\end{array}\right)-\left(\begin{array}{c}{X}_1\\ X_2\\ .\\ .\\ .\\ X_m\end{array}\right)$

令

$V_X={\left(\begin{array}{cccc}{V}_{X_1}& V_{X_2}& ...& V_{X_m}\end{array}\right)}^T,$f_x＝[X₁X₂…X_m]^T，$M={\left(\begin{array}{cccc}{C}_{X_0}& C_{X_1}& ...& C_{X_m}\end{array}\right)}^T$

$B=\left(\begin{array}{cccc}{T}_0\left({\tau }_1\right)& T_1\left({\tau }_1\right)& ...& T_n\left({\tau }_1\right)\\ T_0\left({\tau }_2\right)& T_1\left({\tau }_2\right)& ...& T_n\left({\tau }_2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ T_0\left({\tau }_m\right)& T_1\left({\tau }_m\right)& ...& T_n\left({\tau }_m\right)\end{array}\right)$

则有向量表达式

V_x＝BM-f_x

为了求得唯一解，采用最小二乘法，使所求的切比雪夫多项式系数阵M满 足V_x^TV_x＝min，其中V_x^T为V_x的转置矩阵，求得法方程为

(B^TB)M-B^T f_x＝0

解得

M＝(B^TB)^-1B^T f_x

即得到X坐标分量切比雪夫多项式的拟合系数

同理，对Y、Z位置差重复步骤(1.1)～(1.3)的计算得到Y坐标分量和 Z坐标分量切比雪夫多项式的拟合系数和

(2)将大偏心率椭圆轨道初值和拟合得到的切比雪夫多项式系数和注入卫星；

(3)卫星利用轨道初值计算二体轨道的三轴位置和三轴速度，利用切比雪 夫多项式系数计算三轴位置差和三轴速度差，将其与对应时刻的二体轨道位置 和速度相加，得到卫星的实际位置和速度。

为了更清楚的表明本发明的优点，在此进行数学仿真，仿真条件：火星探 测器半长轴9296.7公里，倾角87.6度，偏心率0.607403，升交点赤经73.5 度，近地点幅角321.0度，真近点角0.0度。

图2为利用精密轨道提供的探测器位置拟合切比雪夫多项式系数，进行环 火大椭圆轨道的拟合结果；图3为利用精密轨道与二体轨道的位置差拟合切比 雪夫多项式系数，进行环火大椭圆轨道的拟合结果。可以看出，采用本发明方 法后拟合精度得到显著提高。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。