CN112141366B - 一种地球轨道中航天器的轨迹预测方法及*** - Google Patents

一种地球轨道中航天器的轨迹预测方法及*** Download PDF

Info

Publication number
CN112141366B
CN112141366B CN202011009507.5A CN202011009507A CN112141366B CN 112141366 B CN112141366 B CN 112141366B CN 202011009507 A CN202011009507 A CN 202011009507A CN 112141366 B CN112141366 B CN 112141366B
Authority
CN
China
Prior art keywords
spacecraft
state information
matching point
time
prediction
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011009507.5A
Other languages
English (en)
Other versions
CN112141366A (zh
Inventor
汪雪川
冯浩阳
刘绮帆
岳晓奎
代洪华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202011009507.5A priority Critical patent/CN112141366B/zh
Publication of CN112141366A publication Critical patent/CN112141366A/zh
Application granted granted Critical
Publication of CN112141366B publication Critical patent/CN112141366B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开的一种地球轨道中航天器的轨迹预测方法，包括：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。本发明采用高精度的反馈加速Picard迭代方法对初始预测状态信息进行修正，无需对矩阵求逆，在提高轨迹预测精度的同时，能够降低轨道计算中的复杂加速度项的计算耗时，并且结合配点法的思想，进一步减少了计算量，本发明为航天器提供了精确、实时的轨道预测方法。

Description

一种地球轨道中航天器的轨迹预测方法及***
技术领域
本发明涉及航天动力学技术领域，特别涉及一种地球轨道中航天器的轨迹预测方法及***。
背景技术
在未来空间任务中，航天器自主程度更高，协同和交互操作更加频繁，研究新型高性能空间轨道计算方法对于航天工程具有重大应用价值。对于轨道递推，研究人员关注的焦点在于计算精度和计算效率。过去的几十年中，已经有大量研究人员对轨道递推计算方法进行研究。传统的数值方法中，为了得到高精度的计算结果，需要选择非常小的计算步长，并构造大规模的代数矩阵方程，这在实际操作中往往难以实现。
发明内容
本发明的目的是提供一种地球轨道中航天器的轨迹预测方法及***，以提高轨迹预测的精度和效果。
为实现上述目的，本发明提供了如下方案：
一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：
从预测时间步长内选取多个配点；
根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；
采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；
对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
可选的，所述从预测时间步长内选取多个配点，具体包括：
利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；
分别令ξ＝ξ1,…,ξM，利用公式
Figure BDA0002697108870000021
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；
其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
可选的，所述根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息，具体包括：
判断预测时间步长是否小于一个轨道周期，获得第一判断结果；
若所述第一判断结果表示是，则，分别令t＝t1,t2,…,tM，根据航天器的初始位置和初始速度，利用公式
Figure BDA0002697108870000022
确定航天器在每个配点时的初始预测状态信息；
其中，x(t)表示实际时间t的初始预测状态信息，r(t)和v(t)分别表示实际时间t的初始预测位置矢量和初始预测速度矢量，r0和v0为航天器的初始位置和初始速度，t1、t2和tM分别表示第1个、第2个和第M个配点；
若所述第一判断结果表示否，则分别令t＝t1,t2,…,tM，根据航天器的初始位置r0和初始速度v0，利用公式
Figure BDA0002697108870000023
确定航天器在每个配点时的初始预测状态信息；
其中，a表示轨道长半轴，a＝-μ/(||v0||2-2μ/||r0||)，e为偏心率，
Figure BDA0002697108870000024
E为偏近点角，通过求解关系式
Figure BDA0002697108870000025
获得，
Figure BDA0002697108870000028
为偏近点角的时间导数，
Figure BDA0002697108870000026
μ为地球引力常数。
可选的，所述采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：
利用公式
Figure BDA0002697108870000027
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；
其中，Un和Un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息；T0、T1和G分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，T0＝diag[-1,-1,...]，
Figure BDA0002697108870000031
G＝EUn-f(Un,t)，J为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tM]，t1、t2和tM分别表示第1个、第2个和第M个配点，f为一阶形式的轨道动力学方程的右端函数，E和
Figure BDA0002697108870000032
分别为第一常值矩阵和第二常值矩阵，E＝diag[(LB)B-1,(LB)B-1,...,(LB)B-1,...,(LB)B-1]，
Figure BDA0002697108870000033
L为微分算子，L-1为积分算子，diag[·]表示对角矩阵，B表示变换矩阵，B＝[Φ(t1)T,Φ(t2)T,...,Φ(tM)T]T，Φ(t)为横向量，Φ(t)＝[T1(t),T2(t)...TM(t)]，t为实际时间，t＝t1,t2,…,tM,，T1(t)、T2(t)和TM(t)分别为横向量Φ(t)中的第1个、第2个和第M个元素，N为插值函数阶数，M为配点个数；
判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；
若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式
Figure BDA0002697108870000035
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；
若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的的状态信息作为航天器在每个配点时的修正后的状态信息。
可选的，所述利用公式
Figure BDA0002697108870000036
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：
利用公式Tm(t)＝cos[m arccos(ξ)]，计算横向量Φ(t)中的每个元素；
其中，Tm(t)表示横向量Φ(t)中的第m个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure BDA0002697108870000034
ti和ti+1为第i个预测时间步长的左端点和右端点。
可选的，所述利用公式
Figure BDA0002697108870000041
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：
利用公式T0(t)＝1计算横向量Φ(t)中的第0个元素T0(t)；
利用公式T1(t)＝ξ计算横向量Φ(t)中的第1个元素T1(t)；
根据横向量Φ(t)中的第0个元素和第1个元素，利用公式Tm+1(t)＝2ξTm(t)-Tm-1(t)，采用迭代的方式计算横向量Φ(t)中的第3个到第M个元素；
其中，Tm-1(t)、Tm(t)和Tm+1(t)分别表示横向量Φ(t)中的第m-1个、第m个和第m+1个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure BDA0002697108870000042
ti和ti+1为第i个预测时间步长的左端点和右端点。
可选的，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，具体包括：
采用Chebyshev多项式作为差值函数，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
可选的，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，之后还包括：
将当前的预测时间步长输出的航天器在每个时刻的修正后的状态信息，作为下一个预测时间步长的初始状态信息，对下个预测时间步长的航天器的轨迹进行预测；所述初始状态信息为航天器的初始位置和初始速度。
一种地球轨道中航天器的轨迹预测***，所述预测***包括：
配点选取模块，用于从预测时间步长内选取多个配点；
初始预测模块，用于根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；
修正模块，用于采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；
插值模块，用于对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
可选的，所述配点选取模块，具体包括：
节点选取子模块，用于利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；
配点选取子模块，用于分别令ξ＝ξ1,…,ξM，利用公式
Figure BDA0002697108870000051
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；
其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。本发明采用高精度的反馈加速Picard迭代方法对初始预测状态信息进行修正，无需对矩阵求逆，在提高轨迹预测精度的同时，能够降低轨道计算中的复杂加速度项的计算耗时，并且结合配点法的思想，进一步减少了计算量，本发明为航天器提供了精确、实时的轨道预测方法。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明提供的一种地球轨道中航天器的轨迹预测方法的流程图；
图2为本发明提供的一种地球轨道中航天器的轨迹预测方法的原理图；
图3为本发明提供的一个预测时间步长的轨道递推结果；
图4为本发明提供的本发明获取的计算结果的Hamilton(哈密顿)量的相对误差；
图5为本发明提供的反馈加速Picard迭代方法与RK12(10)计算结果的位置误差对比图。
具体实施方式
本发明的目的是提供一种地球轨道中航天器的轨迹预测方法及***，以提高轨迹预测的精度和效果。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。
如图1和2所示，本发明提供一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：
步骤101，从预测时间步长内选取多个配点。
配点法也是求解轨道递推问题的一个重要的工具。在配点法中，问题的解是通过多项式、谐波函数、或者其他的基函数估计得到的。相比于其他纯粹的数值方法，配点法的优势在于能够在相对更大的时间尺度上得到半解析的估计解，同时只需要存储较少的基函数系数。Newton-Kantorovich/Chebyshev(牛顿·坎托罗维奇/切比雪夫)伪谱法和RBFs(RadialBasisFunctions，径向基插值方法)方法分别使用Chebyshev(切比雪夫)多项式和径向基函数作为基函数对配点插值，都能够精确高效的求解轨道动力学问题。
1-1选择初始化轨道方案：
根据用户或任务需求，选择单次轨道预测时间(预测时间步长)Δt＝ti+1-ti，其中ti为单次轨道预测初始时刻，ti+1为单次轨道预测末端时刻。
当单次轨道预测时间Δt小于一个轨道周期时，采用冷启动方案，以降低初始化过程计算量，跳转至步骤1-2；
当单次轨道预测时间Δt大于一个轨道周期时，采用热启动方案，以保证后期迭代修正过程能够收敛，跳转至步骤1-3。
1-2冷启动方案：
冷启动方案是利用初始位置r0及初始速度v0得到直线型估计轨道及其整体状态信息x(t)，其中包括轨道位置信息r(t)及速度信息v(t)。实现方法如下：
在已知航天器初始位置r0和速度v0的条件下，假设航天器从起点出发做匀速直线运动，则航天器在t时刻的状态信息x(t)为：
Figure BDA0002697108870000071
x(t)由t时刻的位置矢量r(t)和速度矢量v(t)组成，r0和v0为初始位置和速度，t0为初始时刻。
1-3热启动方案：
热启动方案是利用初始位置r0及初始速度v0得到开普勒模型的估计轨道及其整体状态信息x(t)，其中包括轨道位置信息r(t)及速度信息v(t)。实现方法如下：
基于线性化中心引力场的开普勒轨道为周期性椭圆轨道，具有解析解如下：
||r||＝a(1-e cos E) (2)
E-e sin E＝2π(t/T) (3)
联立以上二式可得轨道位置信息r(t)和速度信息v(t)，从而得到线性化估计轨道整体状态信息x(t)。以上两式中E为偏近点角，a为轨道半长轴，e为偏心率，这些参数均可由经典轨道动力学理论得到。
2-1在估计轨道上选取Chebyshev-Gauss-Lobatto(切比雪夫·高斯·洛巴托，CGL)配点：
基于轨道状态信息x(t)，通过在时间域上[ti,ti+1]选取离散时间点，得到M个离散化的CGL配点状态信息x(tm)，其中m＝1,2,…,M。实现方法如下：
将局部时间区间ti≤t≤ti+1缩放至[-1,1]区间，缩放规则为：
Figure BDA0002697108870000081
其中ti和ti+1为第i个步长的左右端点，t为实际时间，ξ为缩放后的时间。在[-1,1]区间内通过如下公式选取CGL节点
ξm＝cos[(m-1)π/(M-1)],m＝1,...,M(5)
并通过缩放公式反推出时间节点tm，从而得到CGL配点状态信息x(tm)。
步骤101，所述从预测时间步长内选取多个配点，具体包括：利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；分别令ξ＝ξ1,…,ξM，利用公式
Figure BDA0002697108870000082
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
步骤102，根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息。
步骤102所述根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息，具体包括：判断预测时间步长是否小于一个轨道周期，获得第一判断结果；若所述第一判断结果表示是，则，分别令t＝t1,t2,…,tM，根据航天器的初始位置和初始速度，利用公式
Figure BDA0002697108870000083
确定航天器在每个配点时的初始预测状态信息；其中，x(t)表示实际时间t的初始预测状态信息，r(t)和v(t)分别表示实际时间t的初始预测位置矢量和初始预测速度矢量，r0和v0为航天器的初始位置和初始速度，t1、t2和tM分别表示第1个、第2个和第M个配点；若所述第一判断结果表示否，则分别令t＝t1,t2,…,tM，根据航天器的初始位置r0和初始速度v0，利用公式
Figure BDA0002697108870000091
确定航天器在每个配点时的初始预测状态信息；
其中，a表示轨道长半轴，a＝-μ/(||v0||2-2μ/||r0||)，e为偏心率，
Figure BDA0002697108870000092
E为偏近点角，通过求解关系式
Figure BDA0002697108870000093
获得，
Figure BDA0002697108870000094
为偏近点角的时间导数，
Figure BDA0002697108870000095
μ为地球引力常数
步骤103，采用反馈加速Picard(比卡)迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息。
基于步骤102中航天器在每个配点时的初始预测状态信息x(tm)，通过反馈Picard迭代方程对其进行修正，得到精确的C航天器在每个配点时的初始预测状态信息
Figure BDA0002697108870000098
实现方法如下：
选取第一类Chebyshev(切比雪夫)正交多项式作为插值基函数，CGL节点作为时间域节点。通过试算确定基函数个数和配点个数，一般情况下二者相同，从而为基函数系数的确定提供足够多的约束，同时也便于迭代计算并得到高精度解。记基函数个数和配点个数均为M，CGL配点的确定规则为：
ξm＝cos[(m-1)π/(M-1)],m＝1,...,M(6)
(i)构造反馈Picard迭代方程中的常值矩阵E和
Figure BDA0002697108870000099
常值矩阵E的构造表达式为E＝diag[(LB)B-1,(LB)B-1,...,(LB)B-1,...,(LB)B-1]，其中包括M个(LB)B-1；
常值矩阵
Figure BDA0002697108870000096
的构造表达式为
Figure BDA0002697108870000097
其中包括M个(L-1B)B-1。
上式中B＝[Φ(t1)T,Φ(t2)T,...,Φ(tM)T]T，L为微分算子，L-1为积分算子。diag[·]表示以括号中元素构成的对角矩阵。Φ(t)为横向量，表示第一类Chebyshev正交基函数系在t时刻的值，Φ(t)＝[T1(t),T2(t)...TM(t)]，由函数Tn(ξ)＝cos[n arccos(ξ)]或者如下迭代关系确定：
T0(ξ)＝1，T1(ξ)＝ξ，Tn+1(ξ)＝2ξTn(ξ)-Tn-1(ξ) (7)
ξ为缩放后的时间尺度。由于第一类Chebyshev多项式和CGL配点的定义域均为[-1,1],在实际计算中，将局部时间区间ti≤t≤ti+1缩放至[-1,1]区间进行迭代，缩放规则为：
Figure BDA0002697108870000101
其中ti和ti+1为第i个步长的左右端点，t为实际时间，ξ为缩放后的时间，ξ＝t1,t2,...,tM∈[-1,1]为CGL配点。
(ii)构造反馈Picard迭代方程中的变量矩阵T0,T1,G：
变量矩阵T0的表达式为T0＝diag[-1,-1,...]，为MN×MN维单位阵，
变量矩阵T1的表达式为
Figure BDA0002697108870000106
其中J为雅可比阵，
变量矩阵G的表达式为G＝EUn-f(Un,t)，其中t＝[t1,t2,...,tM]。
以上各式中N为插值函数阶数，M为配点个数，一般使M＝N。另外t＝[t1,t2,...,tM]，f为如下一阶形式的轨道动力学方程右端函数，
Figure BDA0002697108870000102
Un为第n次迭代修正的状态矢量，其初始估计U0由步骤2-1中CGL配点状态信息x(tm)，m＝1,2,…,M，合并得到，其中U0＝[x(t1)T,x(t2)T,...x(tM)T]T，上标T代表矩阵转置。
(iii)反馈Picard修正：
将常值矩阵E、
Figure BDA0002697108870000103
变量矩阵T0,T1,G代入下式，对状态矢量Un修正，得到修正状态矢量Un+1
Figure BDA0002697108870000104
上式的推导过程如下：
对于一阶形式的轨道动力学方程组
Figure BDA0002697108870000105
其中x＝x(t),x(t)＝[x1(t),x2(t),...,xd(t),...xN(t)]。根据变分迭代法可得到如下迭代公式
Figure BDA0002697108870000111
其中λ(τ)为待求的拉格朗日乘子。令上式右端项变分为零,得到关于λ(τ)的约束条件：
Figure BDA0002697108870000112
接下来，根据(13)式给出变分迭代法的变型。将λ(τ)估计成Taylor级数形式，其估计式如下：
λ(τ)≈T0[λ]+T1[λ](τ-t)+...+TK[λ](τ-t)K (14)
TK[λ]是λ(τ)的k阶微分变换，即
Figure BDA0002697108870000113
根据约束条件，TK[λ]可以由如下迭代过程得到
Figure BDA0002697108870000114
其中
Figure BDA0002697108870000117
令
Figure BDA0002697108870000118
将λ(τ)的估计式代入原迭代方程，可以得到
Figure BDA0002697108870000115
截取一阶估计式得到
Figure BDA0002697108870000116
通过CGL配点及第一类Chebyshev多项式对其进行离散化后即可得到反馈Picard修正公式。
(iv)判断迭代修正是否收敛：
计算当前修正状态矢量Un+1与上一步状态矢量Un之间的范数误差εn+1＝||Un+1-Un||。设定计算误差限为δ，若εn+1＜δ，结束迭代修正计算，进入步骤(v)；若εn+1≥δ，则转到步骤(ii)继续计算。
(v)提取精确的CGL配点：
从迭代修正收敛的状态矢量U提取精确的CGL配点信息
Figure BDA0002697108870000121
其提取规则为构造规则的逆过程，如下式所示：
Figure BDA0002697108870000122
步骤104所述采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：利用公式
Figure BDA0002697108870000123
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；其中，Un和Un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息；T0、T1和G分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，T0＝diag[-1,-1,...]，
Figure BDA0002697108870000124
G＝EUn-f(Un,t)，J为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tM]，t1、t2和tM分别表示第1个、第2个和第M个配点，f为一阶形式的轨道动力学方程的右端函数，E和
Figure BDA0002697108870000125
分别为第一常值矩阵和第二常值矩阵，E＝diag[(LB)B-1,(LB)B-1,...,(LB)B-1,...,(LB)B-1]，
Figure BDA0002697108870000126
L为微分算子，L-1为积分算子，diag[·]表示对角矩阵，B表示变换矩阵，B＝[Φ(t1)T,Φ(t2)T,...,Φ(tM)T]T，Φ(t)为横向量，Φ(t)＝[T1(t),T2(t)...TM(t)]，t为实际时间，t＝t1,t2,…,tM,，T1(t)、T2(t)和TM(t)分别为横向量Φ(t)中的第1个、第2个和第M个元素，N为插值函数阶数，M为配点个数；判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式
Figure BDA0002697108870000127
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的的状态信息作为航天器在每个配点时的修正后的状态信息。
其中，所述利用公式
Figure BDA0002697108870000128
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：利用公式Tm(t)＝cos[m arccos(ξ)]，计算横向量Φ(t)中的每个元素；其中，Tm(t)表示横向量Φ(t)中的第m个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure BDA0002697108870000131
ti和ti+1为第i个预测时间步长的左端点和右端点。或利用公式T0(t)＝1计算横向量Φ(t)中的第0个元素T0(t)；利用公式T1(t)＝ξ计算横向量Φ(t)中的第1个元素T1(t)；根据横向量Φ(t)中的第0个元素和第1个元素，利用公式Tm+1(t)＝2ξTm(t)-Tm-1(t)，采用迭代的方式计算横向量Φ(t)中的第3个到第M个元素；其中，Tm-1(t)、Tm(t)和Tm+1(t)分别表示横向量Φ(t)中的第m-1个、第m个和第m+1个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure BDA0002697108870000132
ti和ti+1为第i个预测时间步长的左端点和右端点。
步骤104，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
基于步骤103中精确的CGL配点的状态信息
Figure BDA0002697108870000133
以Chebyshev多项式为插值函数，生成精确的预测轨道及其整体状态信息
Figure BDA0002697108870000134
实现方法如下：
Figure BDA0002697108870000135
步骤104所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，具体包括：采用Chebyshev多项式作为差值函数，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
步骤104之后还包括：判断是否结束轨道预测：
若当前轨道预测时间长度满足用户或任务需求，结束轨道预测。若需要做更长时间的轨道预测，基于步骤103中预测轨道整体状态信息
Figure BDA0002697108870000136
获得末端时间节点的位置信息r(tf)和速度信息v(tf)。转到步骤101，并使用末端时间节点的位置信息r(tf)和速度信息v(tf)更新步骤1中的初始位置r0及初始速度v0，重复步骤101至步骤104。
步骤104所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，之后还包括：将当前的预测时间步长输出的航天器在每个时刻的修正后的状态信息，作为下一个预测时间步长的初始状态信息，对下个预测时间步长的航天器的轨迹进行预测；所述初始状态信息为航天器的初始位置和初始速度。
下面结合附图对本发明方法的实施方式作详细说明。
本发明针对低轨情形下的航天器轨道递推问题，提出了一种反馈加速Picard迭代方法。其具体实施方式包括如下步骤：
第一步：对于LEO情形下的轨道递推问题，设计修正变分迭代法，得到轨道递推公式。首先推导航天器运动方程的变分迭代法递推公式，根据泛函取极值条件和泰勒展开，得到拉格朗日乘子的估计形式，将其代入原递推公式，得到修正变分迭代法的递推方程。计算E,
Figure BDA0002697108870000141
T0,T1常值矩阵，代入迭代方程(11)。根据工程需要，配置计算步长和迭代终止条件参数。
第二步：根据初值条件和修正迭代方程在第一个时间步长内进行轨道计算。由初始位置r0和初始速度
Figure BDA0002697108870000142
得到U0，代入迭代方程(11)式进行轨道计算。
第三步：计算状态变量修正量的模，和迭代终止条件对比，判定是否继续进行迭代。将修正量与迭代终止条件比对，若满足迭代终止条件ε≤δ，记录Un+1，转入下一区间进行迭代。若ε＞δ，回到第二步继续迭代。
第四步：根据上一步长终点信息，估计下一步长的初始状态信息，计算下一步长内的运行轨道。根据上一步长终点处位置、速度信息估计下一步长起点处状态信息，从而估计出U0，转入第二步继续计算，直至终点，从而得到航天器在给定时间区间的完整运行轨道。
实例1：在LEO情形下，不考虑空气阻力，利用球谐重力场40阶模型(EMG2008)进行轨道递推，将常值矩阵E,
Figure BDA0002697108870000143
T0,T1及表1中的计算参数代入公式(11)进行计算，得到一个周期内的轨道，如图3所示。为评估本方法的计算精度，对数值结果Hamilton量的变化进行分析。数值结果中Hamilton量的相对误差为εH＝ΔH/H0，其中H0是***在初始时刻的Hamilton量，ΔH为计算过程中Hamilton量的偏差。εH的变化如图4所示。Hamilton量的相对计算误差均被限制在10-13附近，这几乎达到了机器计算精度。
表1 实例1采用的参数值
Figure BDA0002697108870000144
Figure BDA0002697108870000151
实例2：在实例1的重力场模型中加入大气阻力模型，并将计算结果与Runge-Kutta12(10)，即龙格库塔12(10)阶方法进行对比。此处采用的大气阻力模型为简单的球体阻力模型，相关参数如表2所示。其中大气密度通过MSIS-E-90大气模型得到。计算结果的位置误差如图5所示，可以看到在近地轨道，本发明所提方法的最大位置误差为10-6米。计算成本方面，本发明的方法相比于Runge-Kutta12(10)方法可以节约大约95％的计算时间。其中位置误差的计算方法如下：
Figure BDA0002697108870000152
其中Δri是采样点处的位置误差，rA是远地点相对于地心的距离，Norbits为轨道个数。
表2 实例2采用的参数值
参数 数值
RK12(10)的计算步长 50s
节点个数 25
计算步长(T<sub>p</sub>)(s) 500
迭代终止条件 1e<sup>-7</sup>
弹道系数 0.01m<sup>2</sup>/kg
本发明的方法不仅适用于低地球轨道，对其他的轨道类型也是适用的。本发明的前2阶形式，如表3所示。当用户允许的计算时间较长时，可采用0阶估计方法，反之，可采用1阶估计方法。
表3 低阶修正公式
Figure BDA0002697108870000161
本发明还提供一种地球轨道中航天器的轨迹预测***，所述预测***包括：
配点选取模块，用于从预测时间步长内选取多个配点。
所述配点选取模块，具体包括：节点选取子模块，用于利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；配点选取子模块，用于分别令ξ＝ξ1,…,ξM，利用公式
Figure BDA0002697108870000162
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
初始预测模块，用于根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息。
修正模块，用于采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息。
插值模块，用于对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。本发明采用高精度的反馈加速Picard迭代方法对初始预测状态信息进行修正，无需对矩阵求逆，在提高轨迹预测精度的同时，能够降低轨道计算中的复杂加速度项的计算耗时，并且结合配点法的思想，进一步减少了计算量，本发明为航天器提供了精确、实时的轨道预测方法。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Claims (9)

1.一种地球轨道中航天器的轨迹预测方法，其特征在于，所述预测方法包括如下步骤：
从预测时间步长内选取多个配点；
根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；
采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；
所述采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：
利用公式
Figure FDA0003335454650000011
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；
其中，Un和Un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息，T0、T1和G分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，T0＝diag[-1,-1,...]，
Figure FDA0003335454650000012
G＝EUn-f(Un,t)，J为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tM]，t1、t2和tM分别表示第1个、第2个和第M个配点，f为一阶形式的轨道动力学方程的右端函数，E和
Figure FDA0003335454650000013
分别为第一常值矩阵和第二常值矩阵，E＝diag[(LB)B-1,(LB)B-1,...,(LB)B-1,...,(LB)B-1]，
Figure FDA0003335454650000014
L为微分算子，L-1为积分算子，diag[·]表示对角矩阵，B表示变换矩阵，B＝[Φ(t1)T,Φ(t2)T,...,Φ(tM)T]T，Φ(t)为横向量，Φ(t)＝[T1(t),T2(t)...TM(t)]，t为实际时间，t＝t1,t2,…,tM，T1(t)、T2(t)和TM(t)分别为横向量Φ(t)中的第1个、第2个和第M个元素，N为插值函数阶数，M为配点个数；
判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；
若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式
Figure FDA0003335454650000021
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；
若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的状态信息作为航天器在每个配点时的修正后的状态信息；
对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
2.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述从预测时间步长内选取多个配点，具体包括：
利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；
分别令ξ＝ξ1,…,ξM，利用公式
Figure FDA0003335454650000022
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；
其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
3.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息，具体包括：
判断预测时间步长是否小于一个轨道周期，获得第一判断结果；
若所述第一判断结果表示是，则分别令t＝t1,t2,…,tM，根据航天器的初始位置和初始速度，利用公式
Figure FDA0003335454650000023
确定航天器在每个配点时的初始预测状态信息；
其中，x(t)表示实际时间t的初始预测状态信息，r(t)和v(t)分别表示实际时间t的初始预测位置矢量和初始预测速度矢量，r0和v0为航天器的初始位置和初始速度，t1、t2和tM分别表示第1个、第2个和第M个配点，t0为初始时刻；
若所述第一判断结果表示否，则分别令t＝t1,t2,…,tM，根据航天器的初始位置r0和初始速度v0，利用公式
Figure FDA0003335454650000031
确定航天器在每个配点时的初始预测状态信息；
其中，a表示轨道长半轴，a＝-μ/(||v0||2-2μ/||r0||)，e为偏心率，
Figure FDA0003335454650000032
E为偏近点角，通过求解关系式
Figure FDA0003335454650000033
获得，
Figure FDA0003335454650000034
为偏近点角的时间导数，
Figure FDA0003335454650000035
μ为地球引力常数。
4.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述利用公式
Figure FDA0003335454650000036
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：
利用公式Tm(t)＝cos[marccos(ξ)]，计算横向量Φ(t)中的每个元素；
其中，Tm(t)表示横向量Φ(t)中的第m个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure FDA0003335454650000037
ti和ti+1为第i个预测时间步长的左端点和右端点。
5.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述利用公式
Figure FDA0003335454650000038
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：
利用公式T0(t)＝1计算横向量Φ(t)中的第0个元素T0(t)；
利用公式T1(t)＝ξ计算横向量Φ(t)中的第1个元素T1(t)；
根据横向量Φ(t)中的第0个元素和第1个元素，利用公式Tm+1(t)＝2ξTm(t)-Tm-1(t)，采用迭代的方式计算横向量Φ(t)中的第3个到第M个元素；
其中，Tm-1(t)、Tm(t)和Tm+1(t)分别表示横向量Φ(t)中的第m-1个、第m个和第m+1个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，
Figure FDA0003335454650000041
ti和ti+1为第i个预测时间步长的左端点和右端点。
6.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，具体包括：
采用Chebyshev多项式作为差值函数，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
7.根据权利要求1所述的地球轨道中航天器的轨迹预测方法，其特征在于，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，之后还包括：
将当前的预测时间步长输出的航天器在每个时刻的修正后的状态信息，作为下一个预测时间步长的初始状态信息，对下个预测时间步长的航天器的轨迹进行预测；所述初始状态信息为航天器的初始位置和初始速度。
8.一种地球轨道中航天器的轨迹预测***，其特征在于，所述预测***包括：
配点选取模块，用于从预测时间步长内选取多个配点；
初始预测模块，用于根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；
修正模块，用于采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；
所述采用反馈加速Picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：
利用公式
Figure FDA0003335454650000051
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；
其中，Un和Un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息，T0、T1和G分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，T0＝diag[-1,-1,...]，
Figure FDA0003335454650000052
G＝EUn-f(Un,t)，J为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tM]，t1、t2和tM分别表示第1个、第2个和第M个配点，f为一阶形式的轨道动力学方程的右端函数，E和
Figure FDA0003335454650000053
分别为第一常值矩阵和第二常值矩阵，E＝diag[(LB)B-1,(LB)B-1,...,(LB)B-1,...,(LB)B-1]，
Figure FDA0003335454650000054
L为微分算子，L-1为积分算子，diag[·]表示对角矩阵，B表示变换矩阵，B＝[Φ(t1)T,Φ(t2)T,...,Φ(tM)T]T，Φ(t)为横向量，Φ(t)＝[T1(t),T2(t)...TM(t)]，t为实际时间，t＝t1,t2,…,tM，T1(t)、T2(t)和TM(t)分别为横向量Φ(t)中的第1个、第2个和第M个元素，N为插值函数阶数，M为配点个数；
判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；
若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式
Figure FDA0003335454650000055
对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；
若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的状态信息作为航天器在每个配点时的修正后的状态信息；
插值模块，用于对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。
9.根据权利要求8所述的地球轨道中航天器的轨迹预测***，其特征在于，所述配点选取模块，具体包括：
节点选取子模块，用于利用公式ξm＝cos[(m-1)π/(M-1)],m＝1,...,M，在[-1,1]区间中选取M个节点；
配点选取子模块，用于分别令ξ＝ξ1,…,ξM，利用公式
Figure FDA0003335454650000061
确定每个节点对应的第i个预测时间步长内的实际时间作为配点；
其中，ξ1、ξm和ξM表示分别表示第1个、第m个和第M个节点，M表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。
CN202011009507.5A 2020-09-23 2020-09-23 一种地球轨道中航天器的轨迹预测方法及*** Active CN112141366B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011009507.5A CN112141366B (zh) 2020-09-23 2020-09-23 一种地球轨道中航天器的轨迹预测方法及***

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011009507.5A CN112141366B (zh) 2020-09-23 2020-09-23 一种地球轨道中航天器的轨迹预测方法及***

Publications (2)

Publication Number Publication Date
CN112141366A CN112141366A (zh) 2020-12-29
CN112141366B true CN112141366B (zh) 2022-03-25

Family

ID=73898000

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011009507.5A Active CN112141366B (zh) 2020-09-23 2020-09-23 一种地球轨道中航天器的轨迹预测方法及***

Country Status (1)

Country Link
CN (1) CN112141366B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113139233B (zh) * 2021-03-19 2022-02-25 徐州九鼎机电总厂 一种基于沉浸式人机交互的武器弹道仿真方法
CN115204449B (zh) * 2022-05-26 2023-04-25 中国人民解放军国防科技大学 一种基于自适应勒让德皮卡迭代法的轨道预测方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103708045A (zh) * 2014-01-16 2014-04-09 中国人民解放军国防科学技术大学 一种探月飞船跳跃式再入的在线参数辨识方法
CN103863579A (zh) * 2014-03-31 2014-06-18 北京控制工程研究所 一种深空探测返回过程的预测校正制导方法
CN109032176A (zh) * 2018-07-25 2018-12-18 西北工业大学 一种基于微分代数的地球同步轨道确定和参数确定方法
CN109927941A (zh) * 2019-04-08 2019-06-25 北京电子工程总体研究所 一种基于预测离轨点精度的自主允许离轨判断方法
CN111678525A (zh) * 2020-08-11 2020-09-18 北京控制与电子技术研究所 基于互测信息的多航天器自主导航方法、***及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103708045A (zh) * 2014-01-16 2014-04-09 中国人民解放军国防科学技术大学 一种探月飞船跳跃式再入的在线参数辨识方法
CN103863579A (zh) * 2014-03-31 2014-06-18 北京控制工程研究所 一种深空探测返回过程的预测校正制导方法
CN109032176A (zh) * 2018-07-25 2018-12-18 西北工业大学 一种基于微分代数的地球同步轨道确定和参数确定方法
CN109927941A (zh) * 2019-04-08 2019-06-25 北京电子工程总体研究所 一种基于预测离轨点精度的自主允许离轨判断方法
CN111678525A (zh) * 2020-08-11 2020-09-18 北京控制与电子技术研究所 基于互测信息的多航天器自主导航方法、***及装置

Also Published As

Publication number Publication date
CN112141366A (zh) 2020-12-29

Similar Documents

Publication Publication Date Title
CN112141366B (zh) 一种地球轨道中航天器的轨迹预测方法及***
Mueller et al. Covariance correction step for kalman filtering with an attitude
Ainscough et al. Q-method extended Kalman filter
Taheri et al. Co-state initialization for the minimum-time low-thrust trajectory optimization
Mikkola Practical symplectic methods with time transformation for the few-body problem
Oberle et al. Existence and multiple solutions of the minimum-fuel orbit transfer problem
Wang et al. Feedback-accelerated Picard iteration for orbit propagation and lambert’s problem
Wu et al. Adaptive control for spacecraft relative translation with parametric uncertainty
Haug et al. Generalized coordinate partitioning methods for numerical integration of differential-algebraic equations of dynamics
CN106248300A (zh) 基于成对推力器连续工作的卫星质心位置测量方法
Yan et al. High-accuracy trajectory optimization for a trans-earth lunar mission
Peng et al. Multi-objective transfer to libration-point orbits via the mixed low-thrust and invariant-manifold approach
Bradley et al. A continuation method for converting trajectories from patched conics to full gravity models
Wu et al. Approximate time-optimal low-thrust rendezvous solutions between circular orbits
Ellison et al. Analytical partial derivative calculation of the sims-flanagan transcription match point constraints
Stacey et al. Analytical process noise covariance modeling for absolute and relative orbits
Yin et al. Midcourse correction of Earth-Moon distant retrograde orbit transfer trajectories based on high-order state transition tensors
CN110955255A (zh) 基于cmg的高精度轨控姿态维持方法、***及介质
Kluever Efficient computation of optimal interplanetary trajectories using solar electric propulsion
Margolis et al. Nonlinear model predictive control of reentry vehicles based on Takagi-Sugeno fuzzy models
Arya et al. Hyperbolic tangent-based smoothing with state transition matrix implementation for generating fuel-optimal trajectories
CN111209523B (zh) 适用于大偏心率轨道密集星历精密计算的快速处理方法
Lee et al. A combinatorial optimal control problem for spacecraft formation reconfiguration
CN109581356B (zh) 一种常值机动空间目标的约束滤波追踪方法
Ma et al. Velocity corrections to Kepler energy and Laplace integral

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant