CN109255096A - 一种基于微分代数的地球同步卫星轨道不确定演化方法 - Google Patents

Abstract

本发明公开了一种基于微分代数技术的轨道不确定性分析方法，基于多元函数泰勒展开和多项式运算框架，基于地球同步卫星轨道要素描述的动力学模型，在动力学模型上加入太阳光压、第三引力摄动和地球扁率三个摄动力项，将动力学模型的右边项沿着标称轨道进行泰勒展开，得到以初始偏差为变量的展开多项式，在微分代数框架下，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，将初始偏差的具体数值带入多项式结果，即可得到最终航天器的状态，本发明针对不同摄动力，分析了最优的展开阶，平衡计算时间和计算精度；能够应用于快速分析同步卫星在存在初始状态偏差和参数不确定时的轨道演化问题，也可用于其他航天器轨道演化和姿态演化任务中。

Description

一种基于微分代数的地球同步卫星轨道不确定演化方法

技术领域

本发明涉及航空航天领域，特别涉及一种基于微分代数技术的地球同步卫星轨道不确定性演化方法。

背景技术

近年来，随着航天技术的不断发展和对地球同步轨道卫星需求的增加，位于地球同步轨道区域的物体数量不断增加。欧空间最新的空间物体监测报告显示当前大约有1533颗非涉密物体运行在同步轨道区域，其中502颗具有自控制能力，其余的为不可控卫星或空间碎片。为了避免地球同步轨道卫星之间的碰撞，保证它们的安全，监测它们的当前状态和预测它们在未来一段时间内的状态演化具有十分重要的实际价值。

目前，关于地球同步卫星轨道的描述方法主要包括位置和速度矢量描述，经典的轨道六要素描述以及地球同步轨道要素描述。位置和速度矢量描述方法对航天器运动的描述具有明确的物理意义，然而它忽略了同步轨道卫星相对于地球固连坐标系几乎静止的特点。因此，在动力学模型积分过程中需要较小的积分步长，增加了计算成本。传统的开普勒轨道要素考虑了相对地球静止这一特点，解决了小积分步长引起计算成本增加的问题，然而，该轨道描述方法会在轨道倾角和偏心率为0时产生歧义。为了克服上面的缺点，研究人员发展了地球同步轨道要素。一般而言，为了研究航天器轨道的不确定性演化问题，线性演化模型是最常用的方法。但是在强非线性动力学系统中，线性演化模型很容易丢失精度。在不考虑大量计算成本的情形中，可选的解决该问题的方法是非线性的蒙特卡洛打靶方法，该方法能够提供较为准确的结果。随着动力学系统的复杂性增加，计算成本将会限制该方法的应用。今年来，非线性的演化模型得到了快速的发展，企图通过考虑动力学系统的高阶非线性项来提高方法的精度。该方法的缺点是在任意积分时刻都需要计算动力学方程右边项的高阶展开，这是一个复杂而又容易出错的过程。由于地球同步轨道航天器的初始状态和航天器参数存在不可预测的偏差(常假设为高斯分布)，致使无法根据航天器的动力学模型预测未来某一时刻的航天器状态。

发明内容

本发明的目的在于提供一种基于微分代数技术的地球同步卫星轨道不确定性演化方法，以克服现有技术的不足。

为达到上述目的，本发明采用如下技术方案：

一种基于微分代数的地球同步卫星轨道不确定演化方法，包括以下步骤：首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，在动力学模型上加入太阳光压、第三引力摄动和地球扁率三个摄动力项，将加入摄动力项的动力学模型的多项式沿着标称轨道进行泰勒展开，得到以初始偏差为变量的展开多项式，计算展开阶与精度之间的关系，得到各摄动力的最优展开阶，最后，采用微分代数对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，将初始偏差的具体数值带入多项式，即可得到最终航天器的状态。

进一步的，采用地球同步轨道要素描述的动力学模型，地球同步轨道要素根据经典的轨道六要素进行定义：

其中λ表示相对于本初子午线的恒星时角，径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t₀)+ω_e(t-t₀)表示格林尼治恒星时角，ω_e表示地球自转角速度；e_x、e_y表示轨道偏心率e在x，y坐标轴上的投影；Q₁、Q₂是地球同步轨道要素集合的第五第六个要素，定义如上。使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：

其中，a＝(a_r,a_θ,a_h)表示摄动加速度沿着轨道径向、横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。

进一步的，通过地球同步轨道要素得到：

其中，s＝ω+Ω+ν表示航天器的恒星时角，μ表示地球引力常数，表示地球同步轨道要素。

进一步的，使用经典轨道六要素作为桥梁，建立笛卡尔坐标与地球同步轨道要素之间的显式关系：

其中，s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。

进一步的，将加入摄动力项的动力学模型的多项式沿着标称轨道进行泰勒展开，分析太阳光压加速度，其动力学模型为：

其中，a_SRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，r_r表示太阳相对于航天器的相对位置，AU表示一个天文单位，P表示在距太阳1个天文单位处单位面积上的太阳压力，C_r表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量。

进一步的，对方程(4)进行积分过程中，需要将a^SRP,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为

其中，a^SRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，通过以下公式计算：

其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为：

I＝(1,0,0)，J＝(0,1,0)，K＝(0,0,1)

其中，r＝(x,y,z)表示航天器的位置，表示航天器的速度，由方程(8)计算得到。

进一步的，将太阳光压加速度模型在标称轨道上以航天器位置偏差进行展开得到k阶近似多项式，在标称轨道(±δx,±δy,±δz)的邻域内进行均匀分布抽样，对抽样点分别使用k阶近似多项式和精确的太阳光压加速度模型计算各点处太阳光压，并计算它们之间的误差，进而得到平均误差，然后对于不同大小的位置偏差，寻找最小的展开阶，使得平均近似误差小于精度阈值。

进一步的，质量为M的第三引力体对近地航天器的引力摄动加速度表示为：

其中，r_M和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数，a_M,ECI表示第三引力体对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影，将a_M,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：

其中，a_M,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(11)计算得到。

进一步的，在微分代数框架下，对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，并对该结果与准确的数值积分值进行比较，评估其性能。

进一步的，采用微分代数对动力学模型进行多项式形式的积分的详细过程如下，基于微分代数技术，能够得到常微分方程解的任意高阶近似多项式解；其具体过程如下：设一个n维变量组成的常微分动力学系统为：

该微分方程的解表示为x(t)＝Φ(t；t₀,x₀)；假设t₀时刻航天器的标称状态为初始偏差为δx₀，采用该微分方程对某一数值初值的邻域进行积分；其中[x₀]表示多项式，为微分代数变量；在区间[t₀,t₁]上进行积分后得到微分方程在t₁时刻的解Φ(t₁；t₀,x₀+δx₀)的高阶展开式其中将上述积分过程重复执行直到最终时刻t_d，并得到微分方程解的k阶多项式近似解：

其中表示标称轨迹状态；δx₀＝[δx_0,1，…,δx_0,n]^T表示初始偏差，表示展开多项式的系数，在上述过程中，微分代数技术会自动将微分方程的右边项f(x,t)沿着标称轨迹进行高阶泰勒展开，同时，将传统数值积分器中的浮点数运算替换为相应多项式运算。

与现有技术相比，本发明具有以下有益的技术效果：

本发明一种基于微分代数的地球同步卫星轨道不确定演化方法，基于多元函数泰勒展开和多项式运算框架，基于地球同步卫星轨道要素描述的动力学模型，在动力学模型上加入太阳光压、第三引力摄动和地球扁率三个摄动力项，将动力学模型的右边项沿着标称轨道进行泰勒展开，得到以初始偏差为变量的展开多项式，在微分代数框架下，积分器中的所有浮点数运算被多项式运算替代，最终可得到任意时刻以初始偏差为变量的多项式表示的轨道状态，将初始偏差的具体数值带入多项式结果，即可得到最终航天器的状态，本发明针对不同摄动力，分析了最优的展开阶，减少了计算时间，平衡计算时间和计算精度；本发明能够应用于快速分析同步卫星在存在初始状态偏差和参数不确定时的轨道演化问题，也可应用于其他航天器轨道演化和姿态演化任务中。

进一步的，基于微分代数技术的地球同步轨道卫星轨道不确定性演化方法，结合轨道观测数据，可提出高阶卡尔曼滤波方法，提高轨道估计的精度。

附图说明

图1为航天器与第三体最小距离结构示意图。

图2为针对太阳光压加速度模型，k阶近似多项式的平均误差与位置偏差之间的关系示意图。

图3为针对太阳引力摄动加速度，k阶近似多项式的平均误差与位置偏差之间的关系示意图。

图4为针对月球引力摄动加速度，k阶近似多项式的平均误差与位置偏差之间的关系示意图。

图5为针对地球扁率摄动加速度，k阶近似多项式的平均误差与位置偏差之间的关系示意图。

图6为基于微分代数方法的轨道不确定演化流程图。

图7为初始位置偏差或速度偏差存在时，航天器状态的3阶近似多项式的平均位置误差随着时间的演化示意图。

图8为初始位置偏差或速度偏差存在时，航天器状态的3阶近似多项式的平均速度误差随着时间的演化示意图示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，在此基础上，为了更加准确的描述同步轨道卫星的轨道演化，在动力学模型上加入太阳光压、日引力摄动、月引力摄动和地球扁率四个摄动力项，将动力学模型多项式沿着标称轨道进行泰勒展开，得到以初始偏差为变量的展开多项式，并分析展开阶与精度之间的关系，得到各摄动力的最优展开阶，最后，在微分代数框架下，对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，将初始偏差的具体数值带入多项式结果，即可得到最终航天器的状态。将一系列初始偏差数值代数该多项式结果，可以替换非线性蒙特卡洛打靶中的一系列数值积分，实现地球同步轨道状态预测。

具体的包括以下步骤：

步骤1)、动力学模型选择：地球同步卫星在地球固连坐标系中运动缓慢，地球同步卫星轨道要素考虑了这个特征，因此可增加积分步长而不损失精度；

步骤2)、摄动力建模：为了增加轨道演化的精度，按照地球同步轨道的高度，分析得到太阳光压、日引力摄动、月引力摄动和地球扁率摄动是影响航天器轨道最为主要的四个摄动力，因此，该方法考虑这四种摄动力影响；

步骤3)、摄动力的最优展开阶分析：在微分代数框架下，对摄动力模型进行泰勒展开，随着展开阶的增加，近似精度会增加，同时计算成本也会增加，满足某一精度阈值的最低展开阶被称为最优展开阶；

步骤4)、最优积分步长分析：针对某一积分器(固定步长四阶龙格库塔积分方法或者八阶龙格库塔方法)，分析最优积分步长，即满足某一精度阈值的最大积分步长；

步骤5)、方法性能评估：根据上述最优展开阶分析和最优积分步长分析，在微分代数框架下，对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，并对该结果与准确的数值积分值进行比较，评估其性能。

选择动力学模型：

根据任务要求，为了在数值积分过程中使用较大的积分步长且避免动力学模型出现奇异点，采用地球同步轨道要素描述的动力学模型，地球同步轨道要素可根据经典的轨道六要素进行定义：

其中λ表示相对于本初子午线的恒星时角，径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t₀)+ω_e(t-t₀)表示格林尼治恒星时角，ω_e表示地球自转角速度；e_x、e_y表示轨道偏心率e在x，y坐标轴上的投影；Q₁、Q₂是地球同步轨道要素集合的第五第六个要素，定义如上。使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：

其中，a＝(a_r,a_θ,a_h)表示摄动加速度沿着轨道径向，横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。通过地球同步轨道要素得到：

其中，s＝ω+Ω+ν表示航天器的恒星时角，μ表示地球引力常数，表示地球同步轨道要素。。

笛卡尔坐标和地球同步轨道要素之间的关系推导：

在对方程(4)进行积分过程中，需要计算航天器在时刻t的摄动力，然而，摄动力模型是航天器位置矢量的函数，因此，推导笛卡尔坐标和地球同步轨道要素之间显得尤为重要，使用经典轨道六要素作为桥梁，能够得到笛卡尔坐标与地球同步轨道要素之间的显式关系为：

其中，s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。

摄动力最优展开阶的分析方法：

在该框架中，针对地球同步卫星轨道的演化情形，主要考虑太阳光压，日引力摄动、月引力摄动和地球扁率这四个显著影响航天器轨道的摄动力。

首先分析太阳光压加速度，其动力学模型为：

其中，a_SRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，r_r表示太阳相对于航天器的相对位置，AU表示一个天文单位，P表示在距太阳1个天文单位(1AU)处单位面积上的太阳压力，C_r表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量；值得注意的是，在对方程(4)进行积分过程中，需要将a_SRP,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为

其中，a_SRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，可通过以下公式计算

其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为

I＝(1,0,0)，J＝(0,1,0)，K＝(0,0,1)

其中，r＝(x,y,z)表示航天器的位置，表示航天器的速度，可由方程(8)计算得到。

特别需要强调的是类似于太阳光压加速度这类摄动力，它的大小随着航天器到太阳的距离增加而减小，因此，在距离最小时，航天器位置偏差引起的太阳光压大小变化最为显著，如图1所示，采用泰勒多项式近似摄动力模型，在分析针对太阳光压加速度模型的最优展开阶时，根据图1要求，仿真时刻设置为2015年1月5日0:0:0UTC，该时刻日地距离最小。加速度精度阈值设置为ε_a＝10^-11m/s。将太阳光压加速度模型在标称轨道上以航天器位置偏差进行展开得到k阶近似多项式。在标称轨道(±δx,±δy,±δz)的邻域内进行均匀分布抽样，对抽样点分别使用上述k阶近似多项式和精确的太阳光压加速度模型计算各点处太阳光压，并计算它们之间的误差，进而得到平均误差。为了降低计算成本，需要寻找在满足时的最小加速度精度阈值的展开阶，即对于不同大小的位置偏差，寻找最小的展开阶，使得平均近似误差小于精度阈值。图2给出了针对太阳光压模型，k阶近似多项式的平均误差与位置偏差之间的关系。图2表明，当航天器位置偏差小于500km时，太阳光压可被视为常值。

日引力摄动和月引力摄动统称为第三引力体对近地航天器的引力摄动，即第三引力摄动，质量为M的第三引力体(太阳或月球)对近地航天器的引力摄动加速度可以表示为

其中，r_M和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数。a_M,ECI表示第三引力体(太阳或月球)对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影。类似于方程(10)，需要将a_M,ECI投影到航天器轨道的径向，横向和法向方向，其投影关系为

其中，a_M,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(11)计算得到。

采用泰勒多项式近似摄动力模型，在分析针对日月摄动加速度模型的最优展开阶时，根据图1要求，仿真时刻分别设置为2015年1月5日0:0:0UTC和2015年1月21日21:36:0UTC，这两个时刻日地距离和日月距离最小。类似于上面对太阳光压k阶近似多项式的平均误差分析。针对日引力摄动、月引力摄动模型，图3和图4分别给出了k阶近似多项式的平均误差与位置偏差之间的关系。图3表明当航天器位置偏差小于500km时，可使用一阶泰勒展开多项式近似太阳引力摄动加速度模型。图4表明当航天器位置偏差小于100km时，可使用一阶泰勒展开多项式近似月球引力摄动加速度模型；当位置偏差大于100km小于500km时，必须使用二阶泰勒展开多项式近似月球引力摄动加速度模型。

在以地心为中心的地球固连坐标系中，5×5地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影a_{nonspherical,ECF}可表示为

其中，GM表示地球引力常数，系数C_nm,S_nm以及标称的地球半径分别由EGM96地球引力场模型给出。V_n,m和W_n,m可通过下面迭代公式求得

且

其中，r＝(x,y,z)表示航天器在地球固连坐标系中的位置矢量，r＝|r|表示航天器到地心的距离，标称的地球半径由EGM96地球引力场模型给出。类似于方程(10)，需要将a_{nonspherical,ECF}投影到航天器轨道的径向，横向和法向方向，其投影关系为

其中，a_{nonspherical,LVLH}表示地球非球形引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向；a_{nonspherical,ECF}表示地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影；表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(11)计算得到；表示从以地心为中心的地球固连坐标系到以地心为中心的惯性坐标系的转换矩阵，

其中，ω_e表示地球自转角速度，t表示转过的时间。

采用泰勒多项式近似摄动力模型，在分析针对地球非球形引起的引力摄动力模型的最优展开阶时，类似于上面对太阳光压k阶近似多项式的平均误差分析。针对地球非球形引起的引力摄动力模型，图5给出了k阶近似多项式的平均误差与位置偏差之间的关系。图5表明，当航天器位置偏差小于10km时，可使用一阶泰勒展开多项式近似地球非球形摄动加速度；当位置偏差大于10km小于200km时，必须使用二阶泰勒展开多项式近似地球非球形摄动加速度；当位置偏差大于200km小于500km时，必须使用三阶泰勒展开多项式近似地球非球形摄动加速度。

基于微分代数技术的常微分方程多项式解：

基于微分代数技术，能够得到常微分方程解的任意高阶近似多项式解；其具体过程如下：设一个n维变量组成的常微分动力学系统为：

因此，该微分方程的解可表示为x(t)＝Φ(t；t₀,x₀)；假设t₀时刻航天器的标称状态为初始偏差为δx₀，不同于传统的数值积分方法只能从某一数值初值开始积分，微分代数技术能够对某一数值初值的邻域进行积分，其中[x₀]表示多项式，被称为微分代数变量；在区间[t₀,t₁]上进行积分后可得到微分方程在t₁时刻的解Φ(t₁；t₀,x₀+δx₀)的高阶展开式其中将上述积分过程重复执行直到最终时刻t_d，并得到微分方程解的k阶多项式近似解：

其中表示标称轨迹状态；δx₀＝[δx_0,1，…,δx_0,n]^T表示初始偏差，表示展开多项式的系数，在上述过程中，微分代数技术会自动将微分方程的右边项f(x,t)沿着标称轨迹进行高阶泰勒展开，同时，将传统数值积分器中的浮点数运算替换为相应多项式运算。

最优积分步长分析：

在微分代数框架下，针对固定步长的积分器，过大的积分步长使得积分精度下降，而过小的积分步长则容易引起截断误差的累积和增加计算成本。因此，寻找最优的积分步长不仅可以提高积分过程的精度，也能降低计算成本。最优的积分步长被定义为在满足给定的积分精度阈值内，所能使用的最大积分步长。为了分析8阶龙格库塔方法的最优积分步长，航天器初始位置偏差被设置为-10km≤δx,δy,δz≤10km，初始速度偏差被设置为轨道周期为2天，地球同步卫星轨道和相关的卫星参数如表1所示。根据上述针对不同摄动力模型的最优展开阶的分析，为了降低计算成本，应用相应的最优展开阶。表2，表3分别给出了只存在航天器初始位置偏差或者速度偏差时，使用8阶龙格库塔方法在微分代数框架下对地球同步卫星轨道动力学方程进行积分得到的k阶卫星状态近似多项式的平均误差。表2和表3表明，当使用8阶龙格库塔方法对地球同步卫星轨道动力学方程进行多项式积分时，且精度阈值小于1厘米时，2阶的近似多项式是足够精确的，最优的积分步长为13千秒左右。

表1.初始条件和航天器参数

表2.八阶龙格库塔方法演化初始位置偏差的k阶多项式的平均误差(米)

表3.八阶龙格库塔方法演化初始速度偏差的k阶多项式的平均误差(米)

实施例1：地球同步卫星轨道不确定性演化

基于微分代数技术的地球同步卫星轨道不确定性演化流程如图6所示，首先根据测量设备性能和观测精度，确定地球同步卫星的初始状态偏差，一般来说，该偏差可被假设为均值为0，标准差为σ的正态分布；然后根据该正态分布偏差的大小，分析使用多项式近似四类主要摄动力的最优阶，包括太阳光压，日引力摄动、月引力摄动和地球非球形摄动。在此基础上，应用八阶龙格库塔积分方法和最优步长对同步卫星动力学进行积分，而得到最终状态的k阶泰勒近似多项式。最后，对初始状态偏差进行采样，并将采样得到的初始偏差代入k阶泰勒近似多项式，得到最终的航天器状态。另一方面，通过蒙特卡洛打靶方法，对相同的初始条件进行蒙特卡洛打靶得到最终的航天器状态。比较上面两种结果的误差，当误差小于精度阈值时，该方法可以被使用，否则，需要增加展开阶的值，重复上述步骤，直到达到精度阈值。随着展开阶的增加，基于微分代数技术的轨道不确定演化方法可以达到任意精度。以地球同步卫星轨道不确定性演化为例，在本例中，航天器的初始轨道参数如表1所示，我们假设航天器只存在初始的位置偏差或者速度偏差，位置偏差我们假设为均值为0，标准差为σ_x＝σ_y＝σ_z＝1km；速度偏差我们假设为均值为0，标准差为图7和图8分别给出了初始位置偏差或速度偏差存在时，航天器状态的3阶近似多项式的平均位置误差和平均速度误差随着时间的演化示意图。图7图8显示3阶近似多项式在微分代数框架下的位置精度小于1mm，速度精度在10^-8m/s量级上，因此表明该方法是可以取代非线性蒙特卡洛打靶方法的。

Claims (10)

1.一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，包括以下步骤：首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，在动力学模型上加入太阳光压、第三引力摄动和地球扁率三个摄动力项，将加入摄动力项的动力学模型的多项式沿着标称轨道进行泰勒展开，得到以初始偏差为变量的展开多项式，计算展开阶与精度之间的关系，得到各摄动力的最优展开阶，最后，采用微分代数对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态解，将初始偏差的具体数值带入多项式，即可得到最终航天器的状态。

2.根据权利要求1所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，采用地球同步轨道要素描述的动力学模型，地球同步轨道要素根据经典的轨道六要素进行定义：

其中λ表示相对于本初子午线的恒星时角，径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t₀)+ω_e(t-t₀)表示格林尼治恒星时角，ω_e表示地球自转角速度，e_x、e_y表示轨道偏心率e在x，y坐标轴上的投影；Q₁、Q₂是地球同步轨道要素集合的第五第六个要素；使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：

其中，a＝(a_r,a_θ,a_h)表示摄动加速度沿着轨道径向、横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。

3.根据权利要求2所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，通过地球同步轨道要素得到：

其中，s＝ω+Ω+ν表示航天器的恒星时角，μ表示地球引力常数，表示地球同步轨道要素。

4.根据权利要求2所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，使用经典轨道六要素作为桥梁，建立笛卡尔坐标与地球同步轨道要素之间的显式关系：

其中，s＝ω+Ω+ν表示航天器的恒星时角，ω_e表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。

5.根据权利要求1所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，将加入摄动力项的动力学模型的多项式沿着标称轨道进行泰勒展开，分析太阳光压加速度，其动力学模型为：

其中，a_SRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，r_r表示太阳相对于航天器的相对位置，AU表示一个天文单位，P表示在距太阳1个天文单位处单位面积上的太阳压力，C_r表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量。

6.根据权利要求5所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，对方程(4)进行积分过程中，需要将a_SRP,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为

其中，a_SRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，通过以下公式计算：

其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为：

I＝(1,0,0),J＝(0,1,0),K＝(0,0,1)

其中，r＝(x,y,z)表示航天器的位置，表示航天器的速度，由方程(8)计算得到。

7.根据权利要求6所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，将太阳光压加速度模型在标称轨道上以航天器位置偏差进行展开得到k阶近似多项式，在标称轨道(±δx,±δy,±δz)的邻域内进行均匀分布抽样，对抽样点分别使用k阶近似多项式和精确的太阳光压加速度模型计算各点处太阳光压，并计算它们之间的误差，进而得到平均误差，然后对于不同大小的位置偏差，寻找最小的展开阶，使得平均近似误差小于精度阈值。

8.根据权利要求5所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，质量为M的第三引力体对近地航天器的引力摄动加速度表示为：

其中，r_M和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数，a_M,ECI表示第三引力体对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影，将a_M,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：

其中，a_M,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向，表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(11)计算得到。

9.根据权利要求5所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，在微分代数框架下，对动力学模型进行多项式形式的积分，得到任意时刻以初始偏差为变量的多项式表示的轨道状态，并对该结果与准确的数值积分值进行比较，评估其性能。

10.根据权利要求1所述的一种基于微分代数的地球同步卫星轨道不确定演化方法，其特征在于，采用微分代数对动力学模型进行多项式形式的积分，基于微分代数技术，能够得到常微分方程解的任意高阶近似多项式解；其具体过程如下：设一个n维变量组成的常微分动力学系统为：

该微分方程的解表示为x(t)＝Φ(t；t₀,x₀)；假设t₀时刻航天器的标称状态为初始偏差为δx₀，采用该微分方程对某一数值初值的邻域进行积分；其中[x₀]表示多项式，为微分代数变量；在区间[t₀,t₁]上进行积分后得到微分方程在t₁时刻的解Φ(t₁；t₀,x₀+δx₀)的高阶展开式其中将上述积分过程重复执行直到最终时刻t_d，并得到微分方程解的k阶多项式近似解：

其中表示标称轨迹状态；δx₀＝[δx_0,1，…,δx_0,n]^T表示初始偏差，表示展开多项式的系数，在上述过程中，微分代数技术会自动将微分方程的右边项f(x,t)沿着标称轨迹进行高阶泰勒展开，同时，将传统数值积分器中的浮点数运算替换为相应多项式运算。