CN115204449A - 一种基于自适应勒让德皮卡迭代法的轨道预测方法 - Google Patents

Abstract

本发明公开了一种基于自适应勒让德皮卡迭代法的轨道预测方法，包括如下步骤：步骤1，根据航天器的初始位置与初始速度选择航天器的初始预估轨道；步骤2，设置预测时间步长以及预测时间步长内的节点数目；步骤3，基于带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息，并通过插值技术获得航天器在任意时刻的状态信息与速度信息。本发明应用于航天轨道设计领域，通过基于多次迭代后的近似轨迹获得理想的轨迹预测，不仅能提供足够的精度使预测迭代收敛加速，同时还大大减少了计算量。

Description

一种基于自适应勒让德皮卡迭代法的轨道预测方法

技术领域

本发明涉及航天轨道设计技术领域，具体是一种基于自适应勒让德皮卡迭代法的轨道预测方法。

背景技术

卫星轨道预测问题是航天领域一个基本且重要的实际问题，通常只能通过使用数值积分方法来解决。从历史上看，解决此问题的数值积分器可分为单步法或多步法。而这些只有使用较小的步长才能保证计算不发散，且必须采用高阶方法才能保证一定的计算精度，这导致轨迹预测过程需要大量的计算节点，且各个节点的计算难以并行处理。因此在应对长时间轨道计算和复杂动力学模型时，单步法、多步法需要消耗大量计算资源和在线计算时间，在卫星轨道预测方面面临精度差与效率低的问题。

发明内容

针对上述现有技术中的不足，本发明提供一种基于自适应勒让德皮卡迭代法的轨道预测方法，该方法不仅能提供足够的精度使预测迭代收敛加速，同时还大大减少了计算量。

为实现上述目的，本发明提供一种基于自适应勒让德皮卡迭代法的轨道预测方法，包括如下步骤：

步骤1，获得航天器的初始位置与初始速度并选择航天器的初始预估轨道；

步骤2，自适应设置预测时间步长以及预测时间步长内的节点数目；

步骤3，基于带积分反馈项的勒让德皮卡迭代法得到航天器在各时间步长内各节点上的状态信息与速度信息，并通过插值技术获得航天器在任意时刻的状态信息与速度信息。

在其中一个实施例，步骤2中，设置预测时间步长以及预测时间步长内的节点数目具体包括：

步骤2.1，预测时间步长根据真近点角的均匀划分来确定，从近地点开始将每个轨道周期分为X段，并设置每个预测时间步长上节点数目的最大数量为Y₀，其中，X为奇数；

步骤2.2，令X＝3；

步骤2.3，根据指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度；

步骤2.4，以近地点为起始点，沿着开普勒轨道向前传播轨道的1/X，利用Y(Y<Y₀)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数：

若切比雪夫多项式的最后三个系数均小于给定误差的百分之一，则采用当前1/X的轨道作为预测时间步长，每个预测时间步长上节点的数目设为Y；

若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一，则令节点数目Y＝Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y₀，其中，若Y≥Y₀则进行步骤2.5；

步骤2.5，令X＝X+2后重复步骤2.3-2.4，直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一；

其中，若初始和/或最终时间不在预先计算的分段点上，则缩短第一个预测时间步长和/或最后一个预测时间步长。

在其中一个实施例，步骤3中，所述基于带积分反馈项的勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息，具体包括：

步骤3.1，获取航天器的初始状态信息x₀和初始速度信息v₀；

步骤3.2，将轨道预测问题转换为二阶微分方程组的求解问题，采用勒让德皮卡迭代法得到各节点的速度信息近似解以及状态信息近似解；

步骤3.3，采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值，再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及状态信息近似解，得到带积分误差反馈的速度信息近似解与状态信息近似解，反复迭代该过程直到获得满足给定精度的速度信息近似解与状态信息近似解，并将其作为航天器在各个离散的节点上的速度信息与状态信息。

在其中一个实施例，步骤3.2的具体包括：

将轨道预测问题转换为二阶微分方程组的求解的问题，即：

x″(t)＝f(t,x(t),x′(t)),t₀≤t≤t_f (1)

式中，x″(t)表示位置状态函数的二阶导数即加速度函数，x(t)表示位置状态函数，x′(t)表示位置状态函数的一阶导数即速度函数，t表示时间节点，t₀表示预测时间步长的初始时间，t_f表示预测时间步长的终止时间，

表示一个四元向量值函数；

在初始条件x(t₀)＝x₀，x′(t₀)＝v₀下，可以将式(1)改写成为一阶状态空间方程组，为：

将其转换为积分方程，为：

给定初始的状态信息x₀、速度信息v₀和以开普勒轨道作为初始轨道可通过勒让德皮卡迭代产生一列近似解xⁱ(t)、vⁱ(t)，i＝1,2，···，其中，二阶微分方程组的勒让德皮卡迭代公式为：

式中，vⁱ(t)表示第i次勒让德皮卡迭代的速度信息近似解，x^i-1(τ)表示第i-1次勒让德皮卡迭代的状态信息近似解，v^i-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解，xⁱ(t)表示第i次勒让德皮卡迭代的状态信息近似解，vⁱ(τ)表示第i次勒让德皮卡迭代的速度信息近似解；

令

将[t₀,t_f]变换到标准区间[-1,1]，基于该变换，可以将式(4)重写为：

式中，s表示积分变量，g(·)为力函数，g(τ,x(τ))＝f(t(τ),x(t(τ))以及x(τ)＝x(t(τ)；

采用勒让德多项式来逼近未知轨迹xⁱ和沿轨迹x^i-1的力函数g并使用勒让德多项式L_N+1在(-1,1)上的零点来表示状态轨迹，其中，将k次勒让德多项式记为L_k，零点记为τ_j，j＝0,1,···,N，通过勒让德多项式插值来近似力函数，即：

式中，L_k(τ)表示k次勒让德多项式；

使得

给定x^i-1(τ_j)的值，使得系统向量

可以借助前向离散勒让德变换立即计算，为：

式中，

L_k(τ_j)表示k次勒让德多项式在点τ_j上的函数值，ω_j是针对{τ_j}的勒让德-高斯求积公式的权重系数，为：

式中，L′_N+1(τ_j)表示(N+1)次勒让德多项式在点τ_j处的导数值；

通过结合g(τ,x^i-1(τ))的插值及式(5)，得到xⁱ(t)、vⁱ(t)的近似值，即式(5)可以重写为：

因此xⁱ(τ_j)、vⁱ(τ_j)的值可以通过后向离散勒让德变换来得到，为：

式中，

为系数向量；

获取积分矩阵S，将式(10)转换为简化的矩阵向量形式，为：

基于勒让德皮卡迭代，得到状态信息与速度信息第i次迭代的近似解，为：

式中，

表示第i次勒让德皮卡迭代后速度信息的近似解，

表示第i次勒让德皮卡迭代后状态信息的近似解。

在其中一个实施例，通过勒让德多项式的导数递推关系得到系数向量

具体过程为：

勒让德多项式的导数递推关系为：

(2k+1)L_k(s)＝L′_k+1(s)-L′_k-1(s),k≥1

则有：

式中，k表示自然数，L_k(s)表示k次勒让德多项式，L′_k+1(s)表示(k+1)次勒让德多项式的导数，L′_k-1(s)表示(k-1)次勒让德多项式的导数，L₀(s)表示零次勒让德多项式，L₁(τ)表示一次勒让德多项式，L₀(τ)表示零次勒让德多项式，L_k+1(τ)表示(k+1)次勒让德多项式，L_k-1(τ)表示(k-1)次勒让德多项式；

通过逐项解析计算式(9)的积分，即能得到系数向量

为：

式中，

表示公式(6)中的第k个系数，

表示公式(6)中的第(k+1)个系数，

表示公式(6)中的第(k-1)个系数。

在其中一个实施例，积分矩阵S的确定过程具体为：

首先，令：

其中，ω_j高斯积分的权重系数，

根据勒让德多项式和点{τ_j}定义两个矩阵，为：

其中，由于积分会使勒让德多项式的次数增加1次，因此L矩阵与

矩阵的大小不同，L为(N+1)×(N+1)的矩阵，而

为(N+1)×(N+2)的矩阵；

再定义矩阵T，为：

用[v_k]表示一个由向量v_k组成的矩阵，因此，可以将式(7)转换为矩阵形式，为：

其中，

而系数向量

可通过矩阵形式计算，为：

因此可以得到式(10)的矩阵形式，为：

最终得到积分矩阵，为：

其中，积分矩阵S仅取决于N。

在其中一个实施例，在步骤3.3中，得到带积分误差反馈的速度信息近似解与状态信息近似解的具体程为：

状态信息近似解与速度信息近似解的误差为：

又由多元函数的泰勒公式有：

式中，O(·)表示同阶无穷小；

利用式(14)，同时用

近似真实解，得到速度信息的误差近似值，为：

通过速度信息的误差近似值修正勒让德皮卡迭代的近似解，从而得到带积分误差反馈的近似解，为：

在迭代过程中判断是否有max{||xⁱ-x^i-1||,||vⁱ-v^i-1||}＜ε成立，ε为容差，若成立则退出迭代循环，否则继续迭代；

退出迭代时输出状态信息与速度信息的近似解x^J、v^J，即得到航天器在各个离散的节点上的状态信息与速度信息，最终通过插值技术获得航天器在任意时刻的状态信息与速度信息即能完成轨道预测；其中，J为退出迭代循环时的迭代次数。

在其中一个实施例，步骤1中，采用开普勒轨道作为初始预估轨道。

本发明提供的一种基于自适应勒让德皮卡迭代法的轨道预测方法，采用自适应策略确定时间步长及节点数目、基于勒让德高斯积分点利用勒让德正交多项式逼近力函数、结合积分误差反馈项等实现轨道预测的高精度，同时该方法易于实现并行计算，可大大降低轨道预测计算中的计算耗时，提高计算效率。本发明为卫星轨道预测提供了高精度、高效的轨道预测方法。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。

图1为本发明实施例中轨道预测方法的流程图；

图2为本发明实施例中置预测时间步长以及预测时间步长内的节点数目的设置流程图；

图3为本发明实施例中状态信息与速度信息近似解的获取流程图；

图4为本发明实施例中Molniya轨道模拟示意图；

图5为本发明实施例中Hamilton量误差随递推时间的变化示意图。

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。

如图1所示为本实施例公开的一种基于自适应勒让德皮卡迭代法的轨道预测方法，具体包括如下步骤：

步骤1，根据航天器的初始位置与初始速度选择航天器的初始预估轨道，本实施例中，采用开普勒轨道作为预估轨道。

步骤2，设置预测时间步长以及预测时间步长内的节点数目；

步骤3，基于带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息，并通过插值技术获得航天器在任意时刻的状态信息与速度信息。

参考图2，在步骤2中，本实施例中采用一种自适应技术确定预测时间步长和每个预测时间步长上的节点数目，其具体实施过程为：

步骤2.1，预测时间步长根据真近点角的均匀划分来确定，从近地点开始将每个轨道周期分为X段，并设置每个预测时间步长上节点数目的最大数量为Y₀，其中，X为奇数；

步骤2.2，令X＝3；

步骤2.3，根据指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度；

步骤2.4，以近地点为起始点，沿着开普勒轨道向前传播轨道的1/X，利用Y(Y<Y₀)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数：

若切比雪夫多项式的最后三个系数均小于给定误差的百分之一，则采用当前1/X的轨道作为预测时间步长，每个预测时间步长上节点的数目设为Y；

若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一，则令节点数目Y＝Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y₀，其中，若Y≥Y₀则进行步骤2.5；

步骤2.5，令X＝X+2后重复步骤2.3-2.4，直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一；

其中，若初始和/或最终时间不在预先计算的分段点上，则缩短第一个预测时间步长和/或最后一个预测时间步长。

下面结合具体的示例对预测时间步长以及节点数目的设置进行进一步的说明。

例如：从近地点开始将每个轨道周期分为奇数段，从三段开始，最开始将轨道周期划分为三段。每个预测时间步长的最大节点数40。根据用户指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度。以近地点为起始点，沿着开普勒轨道向前传播轨道的三分之一，利用10个切比雪夫点的切比雪夫多项式逼近该段轨道上的引力函数，如果切比雪夫多项式的最后三个系数均小于给定误差的百分之一，则采用三分之一的轨道作为预测时间步长，并将预测时间步长上节点的数目设为10。如果达不到误差要求，则将节点数目加倍。例如20、40个节点，直到达到误差要求或节点数目达到最大值40。若节点数目达到最大值还未满足误差要求，则减少预测时间步长为五分之一的轨道，节点数目再从10个节点开始计算，预测时间步长按下一个奇数分之一的轨道进行减少，依次类推，直至达到误差要求。一旦确定了预测时间步长和每个预测时间步长上的节点数，就可以开始逐段利用勒让德皮卡迭代进行轨道传播计算。

在步骤3，首先采用带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的速度信息与状态信息，在通过插值技术获得航天器在任意时刻的速度信息与状态信息，即完成了航天器的轨道预测。参考图3，其中采用带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息的具有实施过程为：

步骤3.1，获取航天器的初始状态信息x₀和初始速度信息v₀，可以直接获取或通过包含六要素的TLE两行轨道数据转换得到；

步骤3.2，将轨道预测问题转换为二阶微分方程组的求解的问题，采用勒让德皮卡迭代法得到各节点的速度信息近似解以及状态信息近似解；

步骤3.3，采用速度信息的近似解近似真实解，得到速度信息的误差近似值，再通过速度信息的误差近似值修正勒让德皮卡迭代的信息近似解以及状态信息近似解，得到带积分误差反馈的速度信息近似解与状态信息近似解，并将其作为航天器在各个离散的节点上的速度信息与状态信息。

本实施例中，步骤3.2的具体实施过程包括：

将轨道预测问题转换为二阶微分方程组的求解的问题，即：

x″(t)＝f(t,x(t),x′(t)),t₀≤t≤t_f (1)

式中，x″(t)表示位置状态函数的二阶导数即加速度函数，x(t)表示位置状态函数，x′(t)表示位置状态函数的一阶导数即速度函数，t表示时间节点，t₀表示预测时间步长的初始时间，t_f表示预测时间步长的终止时间，

表示一个四元向量值函数；

在初始条件x(t₀)＝x₀，x′(t₀)＝v₀下，可以将式(1)改写成为一阶状态空间方程组，为：

将其转换为积分方程，为：

给定初始的状态信息x₀,速度信息v₀和以开普勒轨道作为初始轨道，可通过勒让德皮卡迭代产生一列近似解xⁱ(t)、vⁱ(t)，i＝1,2，···，其中，二阶微分方程组的勒让德皮卡迭代公式为：

式中，vⁱ(t)表示第i次勒让德皮卡迭代的速度信息近似解，x^i-1(τ)表示第i-1次勒让德皮卡迭代的状态信息近似解，v^i-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解，xⁱ(t)表示第i次勒让德皮卡迭代的状态信息近似解，vⁱ(τ)表示第i次勒让德皮卡迭代的速度信息近似解；

当f是光滑、可积的单值函数且区间[t₀,t_f]的长度足够小时，这一列近似解收敛到精确解，因此，令

将[t₀,t_f]变换到标准区间[-1,1]，基于该变换，可以将式(4)重写为：

式中，s表示积分变量，g(·)为力函数，g(τ,x(τ))＝f(t(τ),x(t(τ))以及x(τ)＝x(t(τ)；

采用勒让德多项式来逼近未知轨迹xⁱ和沿轨迹x^i-1的力函数g并使用勒让德多项式L_N+1在(-1,1)上的零点来表示状态轨迹，其中，将k次勒让德多项式记为L_k，零点记为τ_j，j＝0,1,···,N，通过勒让德多项式插值来近似力函数，即：

式中，L_k(τ)表示k次勒让德多项式；

使得

给定x^i-1(τ_j)的值，使得系统向量

可以借助前向离散勒让德变换立即计算，为：

式中，

L_k(τ_j)表示k次勒让德多项式在点τ_j上的函数值，ω_j是针对{τ_j}的勒让德-高斯求积公式的权重系数，为：

式中，L′_N+1(x_j)表示(N+1)次勒让德多项式在点τ_j处的导数值；

通过结合g(τ,x^i-1(τ))的插值及式(5)，得到xⁱ(t)、vⁱ(t)的近似值，即式(5)可以重写为：

因此xⁱ(τ_j)、vⁱ(τ_j)的值可以通过后向离散勒让德变换来得到，为：

式中，

为系数向量，在具体实施过程中，可以通过勒让德多项式的导数递推关系得到系数向量

其具体实施过程为：

勒让德多项式的导数递推关系为：

(2k+1)L_k(s)＝L′_k+1(s)-L′_k-1(s),k≥1

则有：

式中，k表示自然数，L_k(s)表示k次勒让德多项式，L′_k+1(s)表示(k+1)次勒让德多项式的导数，L′_k-1(s)表示(k-1)次勒让德多项式的导数，L₀(s)表示零次勒让德多项式，L₁(τ)表示一次勒让德多项式，L₀(τ)表示零次勒让德多项式，L_k+1(τ)表示(k+1)次勒让德多项式，L_k-1(τ)表示(k-1)次勒让德多项式；

通过逐项解析计算式(9)的积分，即能得到系数向量

为：

式中，

表示公式(6)中的第k个系数，

表示公式(6)中的第(k+1)个系数，

表示公式(6)中的第(k-1)个系数。

现有的勒让德皮卡迭代中都是采用两个常数矩阵来形成矩阵向量形式，其具有计算量大的弊端，基于此，本实施例提出了仅采用一个常数矩阵来形成矩阵向量形式。其具体实施过程为：

通过获取积分矩阵S，将式(10)转换为简化的矩阵向量形式，为：

基于勒让德皮卡迭代，得到状态信息与速度信息第i次迭代的近似解，为：

式中，

表示第i次勒让德皮卡迭代后速度信息的近似解，

表示第i次勒让德皮卡迭代后状态信息的近似解。

在具体实施过程中，积分矩阵S的确定过程具体为：

首先，令：

其中，ω_j高斯积分的权重系数，

根据勒让德多项式和点{τ_j}定义两个矩阵，为：

其中，由于积分会使勒让德多项式的次数增加1次，因此L矩阵与

矩阵的大小不同，L为(N+1)×(N+1)的矩阵，而

为(N+1)×(N+2)的矩阵；

再定义矩阵T，为：

用[v_k]表示一个由向量v_k组成的矩阵，因此，可以将式(7)转换为矩阵形式，为：

其中，

而系数向量

可通过矩阵形式计算，为：

因此可以得到式(10)的矩阵形式，为：

最终得到积分矩阵，为：

其中，积分矩阵S仅取决于N，即本实施例中可以仅通过一个常值矩阵来实现，使得LPI方法更容易理解，并且节省了一半矩阵向量乘积的计算量。

在步骤3.3中，采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值，再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及状态信息近似解，得到带积分误差反馈的速度信息近似解与状态信息近似解的具体实施过程为：

在步骤3.2中得到状态信息与速度信息第i次迭代的近似解的前提下，此时状态信息近似解与速度信息近似解的误差为：

又由多元函数的泰勒公式有：

式中，O(·)表示同阶无穷小；

利用式(14)，同时用

近似真实解，得到速度信息的误差近似值，为：

通过速度信息的误差近似值修正勒让德皮卡迭代的近似解，从而得到带积分误差反馈的近似解，为：

迭代过程中判断是否有max{||xⁱ-x^i-1||,||vⁱ-v^i-1||}＜ε成立，ε为容差，若成立则退出迭代循环，否则继续迭代；

退出迭代时输出状态信息与速度信息的近似解x^J、v^J，即得到航天器在各个离散的节点上的状态信息与速度信息，最终通过插值技术获得航天器在任意时刻的状态信息与速度信息即能完成轨道预测。其中，J为退出迭代循环时的迭代次数。

对于轨道问题，当使用复杂的引力模型时，雅可比矩阵的精确计算是非常昂贵的。通常，不需要超过4位精度的J，因为积分误差项本身就是随着迭代而衰减得很快。在近地球轨道的情况下，我们发现使用基于前次迭代的近似轨迹获得的理想二体问题的雅可比矩阵就可提供足够的精度使得迭代解收敛加速，使得计算大大减少。

下面结合具体的实例对本发明作出进一步的说明。

参考图4，模拟一个Molniya轨道。在实例中不考虑空气阻力，考虑球谐重力场70阶模型(EMG2008)进行轨道递推，向前递推预测10个开普勒周期，约5.9×10⁵s。为评估本方法的计算精度，对数值结果Hamilton量的变化进行分析。数值结果中Hamilton量的相对误差为

其中H₀是系统在初始时刻的Hamilton量，H(t)为t时刻Hamilton量，e(t)随递推时间的变化如图5所示。Hamilton量的相对计算误差控制在3×10^-14内，达到了机器计算精度。该实例如表1所示。

表1实例1采用的轨道参数及自适应确定的时间步长与节点数目参数

参数	数值
		初始位置r<sub>0</sub>(km)	[9.2210×10<sup>3</sup>,-2.5413×10<sup>-12</sup>,-4.9875×10<sup>-12</sup>]<sup>T</sup>
初始速度v<sub>0</sub>(km/s)	[3.0433×10<sup>-15</sup>,3.9146,7.6829]<sup>T</sup>
		偏心率	0.72
倾斜角	63°
		一个周期分段数	5
节点数目	40
		迭代中止误差	10<sup>-15</sup>

以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。

Claims (8)

1.一种基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，包括如下步骤：

步骤1，获得航天器的初始位置与初始速度并选择航天器的初始预估轨道；

步骤2，自适应设置预测时间步长以及预测时间步长内的节点数目；

步骤3，基于带积分反馈项的勒让德皮卡迭代法得到航天器在各时间步长内各节点上的状态信息与速度信息，并通过插值技术获得航天器在任意时刻的状态信息与速度信息。

2.根据权利要求1所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，步骤2中，设置预测时间步长以及预测时间步长内的节点数目具体包括：

步骤2.1，预测时间步长根据真近点角的均匀划分来确定，从近地点开始将每个轨道周期分为X段，并设置每个预测时间步长上节点数目的最大数量为Y₀，其中，X为奇数；

步骤2.2，令X＝3；

步骤2.3，根据指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度；

步骤2.4，以近地点为起始点，沿着开普勒轨道向前传播轨道的1/X，利用Y(Y<Y₀)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数：

若切比雪夫多项式的最后三个系数均小于给定误差的百分之一，则采用当前1/X的轨道作为预测时间步长，每个预测时间步长上节点的数目设为Y；

若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一，则令节点数目Y＝Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y₀，其中，若Y≥Y₀则进行步骤2.5；

步骤2.5，令X＝X+2后重复步骤2.3-2.4，直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一；

其中，若初始和/或最终时间不在预先计算的分段点上，则缩短第一个预测时间步长和/或最后一个预测时间步长。

3.根据权利要求1所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，步骤3中，所述基于带积分反馈项的勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息，具体包括：

步骤3.1，获取航天器的初始状态信息x₀和初始速度信息v₀；

步骤3.2，将轨道预测问题转换为二阶微分方程组的求解问题，采用勒让德皮卡迭代法得到各节点的速度信息近似解以及状态信息近似解；

步骤3.3，采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值，再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及状态信息近似解，得到带积分误差反馈的速度信息近似解与状态信息近似解，反复迭代该过程直到获得满足给定精度的速度信息近似解与状态信息近似解，并将其作为航天器在各个离散的节点上的速度信息与状态信息。

4.根据权利要求3所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，步骤3.2的具体包括：

将轨道预测问题转换为二阶微分方程组的求解的问题，即：

x″(t)＝f(t,x(t),x′(t)),t₀≤t≤t_f (1)

式中，x″(t)表示位置状态函数的二阶导数即加速度函数，x(t)表示位置状态函数，x′(t)表示位置状态函数的一阶导数即速度函数，t表示时间节点，t₀表示预测时间步长的初始时间，t_f表示预测时间步长的终止时间，

表示一个四元向量值函数；

在初始条件x(t₀)＝x₀，x′(t₀)＝v₀下，可以将式(1)改写成为一阶状态空间方程组，为：

将其转换为积分方程，为：

给定初始的状态信息x₀、速度信息v₀和以开普勒轨道作为初始轨道可通过勒让德皮卡迭代产生一列近似解xⁱ(t)、vⁱ(t)，i＝1,2，…，其中，二阶微分方程组的勒让德皮卡迭代公式为：

式中，vⁱ(t)表示第i次勒让德皮卡迭代的速度信息近似解，x^i-1(τ)表示第i-1次勒让德皮卡迭代的状态信息近似解，v^i-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解，xⁱ(t)表示第i次勒让德皮卡迭代的状态信息近似解，vⁱ(τ)表示第i次勒让德皮卡迭代的速度信息近似解；

令

将[t₀,t_f]变换到标准区间[-1,1]，基于该变换，可以将式(4)重写为：

式中，s表示积分变量，g(·)为力函数，g(τ,x(τ))＝f(t(τ),x(t(τ))以及x(τ)＝x(t(τ)；

采用勒让德多项式来逼近未知轨迹xⁱ和沿轨迹x^i-1的力函数g并使用勒让德多项式L_N+1在(-1,1)上的零点来表示状态轨迹，其中，将k次勒让德多项式记为L_k，零点记为τ_j，j＝0,1,…,N，通过勒让德多项式插值来近似力函数，即：

式中，L_k(τ)表示k次勒让德多项式；

使得

给定x^i-1(τ_j)的值，使得系统向量

可以借助前向离散勒让德变换立即计算，为：

式中，

L_k(τ_j)表示k次勒让德多项式在点τ_j上的函数值，ω_j是针对{τ_j}的勒让德-高斯求积公式的权重系数，为：

式中，L′_N+1(τ_j)表示(N+1)次勒让德多项式在点τ_j处的导数值；

通过结合g(τ,x^i-1(τ))的插值及式(5)，得到xⁱ(t)、vⁱ(t)的近似值，即式(5)可以重写为：

因此xⁱ(τ_j)、vⁱ(τ_j)的值可以通过后向离散勒让德变换来得到，为：

式中，

为系数向量；

获取积分矩阵S，将式(10)转换为简化的矩阵向量形式，为：

基于勒让德皮卡迭代，得到状态信息与速度信息第i次迭代的近似解，为：

式中，

表示第i次勒让德皮卡迭代后速度信息的近似解，

表示第i次勒让德皮卡迭代后状态信息的近似解。

5.根据权利要求4所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，通过勒让德多项式的导数递推关系得到系数向量

具体过程为：

勒让德多项式的导数递推关系为：

(2k+1)L_k(s)＝L′_k+1(s)-L′_k-1(s),k≥1

则有：

式中，k表示自然数，L_k(s)表示k次勒让德多项式，L′_k+1(s)表示(k+1)次勒让德多项式的导数，L′_k-1(s)表示(k-1)次勒让德多项式的导数，L₀(s)表示零次勒让德多项式，L₁(τ)表示一次勒让德多项式，L₀(τ)表示零次勒让德多项式，L_k+1(τ)表示(k+1)次勒让德多项式，L_k-1(τ)表示(k-1)次勒让德多项式；

通过逐项解析计算式(9)的积分，即能得到系数向量

为：

式中，

表示公式(6)中的第k个系数，

表示公式(6)中的第(k+1)个系数，

表示公式(6)中的第(k-1)个系数。

6.根据权利要求5所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，积分矩阵S的确定过程具体为：

首先，令：

其中，ω_j高斯积分的权重系数，

根据勒让德多项式和点{τ_j}定义两个矩阵，为：

其中，由于积分会使勒让德多项式的次数增加1次，因此L矩阵与

矩阵的大小不同，L为(N+1)×(N+1)的矩阵，而

为(N+1)×(N+2)的矩阵；

再定义矩阵T，为：

用[v_k]表示一个由向量v_k组成的矩阵，因此，可以将式(7)转换为矩阵形式，为：

其中，

而系数向量

可通过矩阵形式计算，为：

因此可以得到式(10)的矩阵形式，为：

最终得到积分矩阵，为：

其中，积分矩阵S仅取决于N。

7.根据权利要求4或5或6所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，在步骤3.3中，得到带积分误差反馈的速度信息近似解与状态信息近似解的具体程为：

状态信息近似解与速度信息近似解的误差为：

又由多元函数的泰勒公式有：

式中，O(·)表示同阶无穷小；

利用式(14)，同时用

近似真实解，得到速度信息的误差近似值，为：

通过速度信息的误差近似值修正勒让德皮卡迭代的近似解，从而得到带积分误差反馈的近似解，为：

在迭代过程中判断是否有max{||xⁱ-x^i-1||,||vⁱ-v^i-1||}＜ε成立，ε为容差，若成立则退出迭代循环，否则继续迭代；

退出迭代时输出状态信息与速度信息的近似解x^J、v^J，即得到航天器在各个离散的节点上的状态信息与速度信息，最终通过插值技术获得航天器在任意时刻的状态信息与速度信息即能完成轨道预测；其中，J为退出迭代循环时的迭代次数。

8.根据权利要求1至6任一项所述基于自适应勒让德皮卡迭代法的轨道预测方法，其特征在于，步骤1中，采用开普勒轨道作为初始预估轨道。

Images