利用径向力平衡飞行控制的航天器高精度轨道捕获方法

孟云鹤; 侯佳睿; 罗宇飞; MENG Yunhe; HOU Jiarui; LUO Yufei

doi: 10.11887/j.cn.202402003

孟云鹤¹ ，侯佳睿² ，罗宇飞¹

1. 中山大学人工智能学院，广东珠海 519082

2. 中山大学物理与天文学院，广东珠海 519082

基金项目: 国家自然科学基金资助项目（61673390）；广东省基础与应用基础研究重大资助项目（2019B030302001）

详细信息

作者简介

孟云鹤（1978—），男，河北永清人，教授，博士，博士生导师，E-mail: mengyh7@mail.sysu.edu.cn

通讯作者

侯佳睿（1996—），女，陕西宝鸡人，硕士，E-mail: 503599729@qq.com

中图分类号: V412.4

文献标识码: A

文章编号: 1001-2486(2024)02-027-09

High-precision orbit capture method of spacecraft based on radial equilibrium flight control

MENG Yunhe¹ ， HOU Jiarui² ， LUO Yufei¹

1. School of Artificial Intelligence, Sun Yat-sen University, Zhuhai 519082 ， China

2. School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519082 , China

摘要

为实现对探测器轨道形状与高度的精准调整，提出一种径向力平衡飞行的航天器连续推力控制新方法。建立连续推力平衡飞行的动力学极坐标模型，并推导出特殊条件下的解析轨道解，进一步分析边值条件，给出连续推力的控制律。利用这一平衡飞行控制理论，构建轨道捕获的最优控制策略。考虑推力器的推力水平，通过一次或多次的控制过程，实现对轨道形状、轨道高度及轨道相位的综合调整。数值仿真表明：利用平衡飞行的轨道控制方法，配置微小推力器的空间引力波探测器可以实现高精度的轨道捕获；该方法具有控制过程可解析、计算量小、简便、实用等特点。

关键词

空间引力波探测器 / 连续推力 / 平衡飞行 / 轨道捕获 / 最优控制

Abstract

In order to acquire the accurate adjustment of the orbit shape and height of the detector, a new method of spacecraft continuous thrust control based on radial force equilibrium flight was proposed. The dynamic polar coordinate model of continuous thrust equilibrium flight was established, and the analytical orbit solution under special conditions was derived, the boundary conditions were further analyzed, and the control law of continuous thrust was given. Using this equilibrium flight control theory, the optimal control strategy for orbit capture can be constructed. Considering the thrust level of thrusters, the integrated adjustment of orbit shape, orbit height and orbit phase can be adjusted through one or more control processes. Numerical simulations show that the space gravitational wave detector with micro thruster can achieve high-precision orbit capture by using the orbit control method of equilibrium flight. This method has the advantages of analytical control process, small calculation, simple and practical.

Keywords

space gravitational wave detector / continuous thrust / equilibrium flight / orbit capture / optimal control

1 轨道捕获任务建模 1.1 轨道机动任务描述 1.2 坐标系建立 1.3 平衡飞行轨道控制方法 2 动力学方程的解析求解 2.1 化简为椭圆积分进行求解 2.2 边值条件分析 2.3 控制律——推力加速度矢量的求解 2.3.1 求解推力加速度倾角Θ_P的步骤 2.3.2 将推力矢量在惯性系表达 3 平衡飞行轨道控制方法的应用 3.1 航天器角速度（相位角）调整 3.2 轨道形状调整 4 仿真验证 4.1 仿真初值设置 4.2 仿真结果 4.2.1 由椭圆轨道转移到圆轨道 4.2.2 由椭圆轨道经多次控制转移到圆轨道 5 结论

引力波是爱因斯坦相对论的重要预言之一，开展引力波探测是人类认识宇宙的重要手段。2016年初，美国LIGO合作组宣布探测到引力波，人类进入利用引力波探测宇宙的新时期。相比地面探测装置，空间引力波探测具有波源类型丰富、时空尺度广、探测精度高等特殊优势，将打开毫赫兹频段的引力波探测窗口^[1-3]。当前国际知名的空间引力波探测计划有欧美的LISA计划^[3-5]，中国的天琴计划^[6]、太极计划^[7]等。其中，LISA计划、太极计划采用日心轨道编队方案，而天琴计划采取地心轨道星座方案。不管哪种方案，探测系统都对探测器的空间分布构形精度提出了极高要求^[8]。根据空间引力波探测器特点，为实现高精度控制，需要配置微推力器^[9]。因而，基于微推力的轨道精确控制是其中的关键技术之一。

轨道的捕获控制是星座构形控制的基础。根据推力器配置的不同，实现航天器轨道捕获的方式有冲量式机动和有限推力机动等方式。冲量式轨道机动计算简便，学者们对冲量式轨道机动设计已经进行了大量研究^[10]。冲量式机动要求航天器的速度变化在瞬间完成^[11]，实际工程中很难实现理想的冲量推力，冲量式轨道机动存在着燃耗较大，精度不高的缺点。随着航天任务对精度要求的不断提高以及小推力技术的进步，近年来许多学者对小推力轨道机动技术进行了很多研究^[12-13]。小推力轨道机动技术拓展了轨道动力学理论研究范畴：Yang等使用连续推力在小行星附近生成人工平动点，并开展相关研究^[14]；EL Mabsout等通过优化连续小推力方向给出了卫星交会的新方法^[15]；Kuettel等提出了基于连续推力的机动瞄准方法^[16]；Dachwald等使用进化神经控制研究了连续推力最优轨迹问题^[17]；Maestrini等对发动机施加推力时的不确定性因素进行分析，提出了一些应对方法^[18]；Izzo等在地球到金星的连续小推力星际转移背景中，用深度神经网络方法表示最优转移问题，得到优化的结果^[19]；Zavoli等使用强化学习方法对存在强不确定性和干扰的小推力星际转移问题开展研究^[20]。小推力轨道控制设计方法主要分为直接法和间接法。直接法将轨道设计问题转化为参数优化问题，可采用非线性规划算法进行求解^[21-22]；Pontani等提出了一种基于李雅普诺夫稳定性的连续小推力反馈控制律，对非线性地球轨道控制进行研究^[23]；此外，他们还采用小推力与非线性轨道控制结合的方法对低月球轨道捕获任务开展研究，结果表明采用上述控制方法可以有效实现航天器的月球轨道捕获^[24]。间接法直接利用庞特里亚金极大值原理将轨道设计问题转化为两点边值问题，而后者一般采用数值打靶进行求解^[25-26]，Hernando-Ayuso等在庞特里亚金极大值原理的基础上研究了圆轨道航天器最优小推力避碰问题，给出了较为准确的解析解^[27]。无论是直接法、间接法，都离不开数值求解，几乎都需要较为准确的初值猜测以及多次迭代计算。初值猜测不准确又会导致结果不收敛，多次迭代计算会带来巨大的计算量。因此，在目前的小推力轨道控制研究领域，可解析求解的轨道控制策略有重要意义。

本文面向地心轨道空间引力波探测器，研究轨道精确捕获问题的解析算法。提出了平衡飞行轨道控制方法，建立平衡飞行条件下基于连续小推力的轨道运动方程；进行轨道积分，并将轨道积分转化为椭圆积分进行求解，进而可以解析计算连续小推力的轨道控制律，该控制律可以在航天器轨道高度不变的情况下，实现瞬时轨道大小与形状（长半轴、偏心率等）的调整；通过数值仿真验证了本方法的可行性。

1 轨道捕获任务建模

1.1 轨道机动任务描述

本文的研究对象是，地心轨道的空间引力波探测器，以天琴计划为例，三星均运行在地心距约10⁵ km的地心圆轨道上，呈等边三角形排列。在科学探测期间，探测精度与分辨率对星座构形稳定性提出了极高要求^[8]，具体的构形稳定性指标见表1。

表1表征等边三角形稳定性的指标要求

Tab.1 Index requirements of the stability of the equilateral triangle

由于存在入轨初始偏差以及在轨运行中的摄动，航天器的瞬时轨道可能与其标称轨道不一致，若不施加轨道控制，卫星构形会很快偏离期望构形，需要进行单星轨道捕获和多星构形捕获，使其进入期望构形并能够稳定保持。欲使三颗卫星构成的星座构形长期稳定，首先需要三星各自的轨道都为高近圆轨道，即偏心率接近为零。

本文重在解决轨道形状调整的问题，天琴三星星座如图1所示，目标是在大小确定但方向可调的连续微推力作用下，将卫星的瞬时轨道由椭圆形调整为高近圆形。

图1天琴三星星座示意图

Fig.1Schematic diagram of TianQin 3 satellites constellation

1.2 坐标系建立

地心惯性坐标系OXYZ：原点O位于地球中心，OX轴沿地球赤道面和黄道面交线指向某历元的春分点方向，OZ轴指向地球北极，OY轴由右手定则确定。

极坐标系：如图2所示建立极坐标系，以地心为极点，从地球指向航天器初始位置的射线为极轴，以轨道所在平面为极平面，规定航天器飞行方向为角度变化的正方向。航天器的位置及速度可用极径与极角表示。

图2轨道与极坐标示意图

Fig.2Illustration of orbit and polar coordinates

1.3 平衡飞行轨道控制方法

平衡飞行指航天器在飞行过程中保持在径向方向的受力平衡，如果径向初速度也为零，则航天器的轨道高度（或地心距离）不变，势能稳定，动能的增减仅来自推力加速度的施加。基于平衡飞行条件，可使瞬时轨道为椭圆的航天器，在轨道高度不变的情况下实现轨道形状控制，即进入高近圆轨道，而这正是空间引力波探测器对轨道形状的基本要求。

在极坐标系下描述航天器的质心运动，列写径向运动方程以及周向运动方程：

\{\begin{matrix} \frac{d^{2} r}{d t^{2}} = r {(\frac{d β}{d t})}^{2} - \frac{μ}{r^{2}} + a_{T r} \\ \frac{d}{d t} (r^{2} \frac{d β}{d t}) = r a_{T θ} \end{matrix}

(1)

其中：r为航天器地心矢径，β为t时刻航天器转过的地心角，μ为地球引力常数，a_T_r为推力加速度矢量的径向分量，a_T_θ为推力加速度矢量的周向分量。

连续推力加速度大小恒定，方向可以调整，定义推力加速度矢量（大小为a_T）与当地水平面夹角为当地推力加速度倾角，记为Θ_P。则推力加速度可分解为：

\{\begin{matrix} a_{T r} = a_{T} s i n Θ_{P} \\ a_{T θ} = a_{T} c o s Θ_{P} \end{matrix}

(2)

将式（2）代入式（1）可得：

\{\begin{matrix} \frac{d^{2} r}{d t^{2}} = r {(\frac{d β}{d t})}^{2} - \frac{μ}{r^{2}} + a_{T} s i n Θ_{P} \\ \frac{d}{d t} (r^{2} \frac{d β}{d t}) = r a_{T} c o s Θ_{P} \end{matrix}

(3)

定义平衡飞行条件，在极坐标系下航天器径向受力平衡，径向速度和径向加速度都为零，即：

\{\begin{matrix} \frac{d^{2} r}{d t^{2}} = 0 \\ \frac{d r}{d t} = 0 \end{matrix}

(4)

将式（4）代入式（3）并化简，可得：

\{\begin{matrix} \frac{a_{T}}{r} s i n Θ_{P} = \frac{μ}{r^{3}} - {(\frac{d β}{d t})}^{2} \\ \frac{a_{T}}{r} c o s Θ_{P} = \frac{d^{2} β}{d t^{2}} \end{matrix}

(5)

根据平衡飞行条件，认为地心距离恒定，在轨道机动过程中不变，这也限定了这种理想的转移应该满足近地点与远地点的时机约束。

如果飞行器从椭圆轨道向圆轨道过渡，r_f为目标圆轨道的地心距，圆轨道的轨道角速度定义为

\sqrt{μ / r_{f}^{3}} ≜ n

，将椭圆轨道的瞬时角速度定义为

ω ≜ d β / d t

，将式（5）两边取平方和，可以得到：

{(\frac{a_{T}}{r})}^{2} = {(n^{2} - ω^{2})}^{2} + {(\frac{d ω}{d t})}^{2}

(6)

因此

\pm \frac{d ω}{\sqrt{{(\frac{a_{T}}{r})}^{2} - {(n^{2} - ω^{2})}^{2}}} = d t

(7)

对两边积分

t = \int_{0}^{t} d t = \pm \int_{ω}^{n} \frac{d ω}{\sqrt{{(\frac{a_{T}}{r})}^{2} - {(n^{2} - ω^{2})}^{2}}}

(8)

如果式（8）可以计算，则在大小为a_T、倾角为Θ_P的推力加速度作用下，均可以求解从椭圆轨道向圆轨道的平衡飞行的轨道机动时间Δt与过程参数。

2 动力学方程的解析求解

2.1 化简为椭圆积分进行求解

对式（8）中积分项的分母进行变换，定义航天器的引力加速度为

g_{0} ≜ μ / r^{2}

，并令推力加速度与引力加速度之比为

ρ ≜ a_{T} / g_{0}

，可以做如下变换：

\sqrt{{(\frac{a_{T}}{r})}^{2} - {(n^{2} - ω^{2})}^{2}} = n^{2} \sqrt{ρ^{2} - {(1 - \frac{ω^{2}}{n^{2}})}^{2}}

(9)

基于此可以进行一系列推导。

1）当ω＜n时，设推力加速度与圆轨道引力加速度之比定义为

ρ ≜ \frac{a_{C}}{g_{0}} = \frac{1}{k}

，并令

τ^{2} = 1 - \frac{ω^{2}}{n^{2}}

，经一系列变换推导，当ω＜n时，可得：

t = \int_{0}^{t} d t = \frac{k}{2 n} \int_{0}^{τ_{0}^{2}} \frac{d α}{\sqrt{(1 - α) (1 - k^{2} α^{2})}}

(10)

式中：

τ_{0}^{2} = 1 - \frac{ω^{2}}{n^{2}}; k = \frac{1}{ρ} = \frac{g_{0}}{a_{T}} ， g_{0} = \frac{μ}{r^{2}}

；α为积分变量。经分析可知，式（10）属于第一类椭圆积分^[28]。其定积分可以求解，可用函数EllipticF函数表示为：

t = \frac{1}{n} \sqrt{\frac{k}{1 + k}} EllipticF {(a r c s i n \sqrt{\frac{1 + k α}{2}}, \frac{2}{1 + k})|}_{α = 0}^{α = 1 - \frac{ω^{2}}{ω_{n}^{2}}}

(11)

2）当ω＞n时，可得

t = \int_{0}^{t} d t = \frac{k}{2 n} \int_{0}^{τ_{0}^{2}} \frac{d α}{\sqrt{(1 + α) (1 - k^{2} α^{2})}}

(12)

式中，

τ_{0}^{2} = \frac{ω^{2}}{n^{2}} - 1

。对此进行积分，可用函数EllipticF表示为：

t = \frac{1}{n} \sqrt{\frac{k}{k + 1}} EllipticF {(a r c s i n \sqrt{\frac{1 + k α}{2}}, \frac{2}{k + 1})|}_{α = 1 - \frac{ω^{2}}{ω_{n}^{2}}}^{α = 0}

(13)

2.2 边值条件分析

分析式有解的条件是：

{(\frac{a_{T}}{r})}^{2} - {(n^{2} - ω^{2})}^{2} ⩾ 0

(14)

进而有

\frac{a_{T}}{r} ⩾ |n^{2} - ω^{2}|

(15)

其代表的物理意义是：在平衡飞行控制过程中，需要施加的推力加速度与可实现调整的角速度范围存在约束关系。

比如，对于标称圆轨道半径为r、角速度为n并且推力值大小为a_T的控制过程来说，航天器初始轨道角速度需要满足一定范围才能够实现轨道机动任务，由式（15）可知，初始角速度范围为：

ω \in [\sqrt{n^{2} - \frac{a_{T}}{r}}, \sqrt{n^{2} + \frac{a_{T}}{r}}]

(16)

对式（16）进行变换，式（16）两边同除以n²，注意到前文定义过

n ≜ \sqrt{μ / r_{f}^{3}}

，以及

g_{0} ≜ μ / r^{2}

，可以得到：

\frac{a_{T}}{g_{0}} ⩾ |1 - \frac{ω^{2}}{n^{2}}|

(17)

定义椭圆轨道瞬时速度为v_E，标称圆轨道速度为v_C。假设飞行过程的轨道地心距不变，故而有：

\frac{ω^{2}}{n^{2}} = \frac{v_{E}^{2}}{v_{C}^{2}}

(18)

活力公式给出了航天器瞬时速度、轨道地心距与瞬时轨道半长轴的关系：

v^{2} = μ (\frac{2}{r} - \frac{1}{a})

(19)

将式（19）代入式（18），化简可得：

\frac{ω^{2}}{n^{2}} = \frac{μ (\frac{2}{r} - \frac{1}{a_{E}})}{μ (\frac{2}{r} - \frac{1}{a_{C}})} = 2 - \frac{r}{a_{E}}

(20)

将式（20）代入式（16），则实现调整的半长轴范围为：

a \in [\frac{r}{1 + a_{T} / g_{0}}, \frac{r}{1 - a_{T} / g_{0}}]

(21)

式（21）的含义是：半长轴在这个范围内，且近地点地心距或远地点地心距为r的椭圆轨道，可以用上述方法调整到半径为r的圆轨道。

式（16）与式（21）本质相同，表达的都是某确定推力加速度的调整能力。若不满足这一条件，则说明调整任务的目标范围过大或发动机推力过小，此时可以通过增大推力值或者采取多次控制的方法，实现控制目标。

2.3 控制律——推力加速度矢量的求解

采取的是方向可调的定值连续推力，故而轨道控制律为整个控制过程中推力方向关于时间的函数。本模型是平面模型，且只以调整轨道形状为目标，暂不涉及轨道面方向的调整，控制律可表达为推力倾角关于时间的函数。

2.3.1 求解推力加速度倾角Θ_P的步骤

1）对于从ω₀→n的控制过程，等间距给出角速度变化序列，ω=ω₀∶Δω ∶n。

2）根据式（11）（ω＜n时）和式（13）（ω＞n时），计算时间序列[t₀，t₁，···，t_n]。

3）前文定义

d β / d t ≜ ω

，将式（5）改写为：

\{\begin{matrix} \frac{a_{T}}{r} s i n Θ_{P} = \frac{μ}{r^{3}} - ω^{2} \\ \frac{a_{T}}{r} c o s Θ_{P} = \frac{d ω}{d t} \end{matrix}

(22)

由式（22）可确定推力加速度倾角Θ_P的大小与象限。

4）容易计算出t-ω-Θ_P序列。

2.3.2 将推力矢量在惯性系表达

定义如下物理量：H表示航天器所在位置的当地水平线，方向沿航天器飞行方向（单位矢量）；n表示轨道平面的法向量（单位矢量）。

使H绕n在轨道平面内旋转Θ_P，可得推力加速度矢量的方向，进一步利用坐标变换可计算得出推力加速度矢量在地心惯性系下的表达，计算方法为：

a_{T} = |a_{T}| \cdot [H c o s Θ_{P} + (n \cdot H) n (1 - c o s Θ_{P}) + (n \times H) s i n Θ_{P}]

(23)

3 平衡飞行轨道控制方法的应用

在平衡飞行条件及上述推力作用下，航天器可以在轨道高度不变的情况下，通过施加连续推力，连续地改变关键的轨道参数。

3.1 航天器角速度（相位角）调整

前文计算t-ω-Θ_P序列的过程即可直观表达航天器角速度的调整过程。当推力矢量与飞行方向夹角为锐角时，航天器的运动为加速过程；当推力矢量与飞行方向夹角为钝角时，航天器的运动为减速过程。

由于要满足平衡飞行条件，无论是加速过程还是减速过程，始终保持航天器径向受力平衡，所以航天器轨道高度不变，做变速圆周运动，从而可以实现相位角的调整。由二体运动理论可得出如下结论：在变速圆周运动过程中，航天器的瞬时轨道半长轴在连续地变化，加速过程轨道半长轴持续增加，减速过程轨道半长轴则持续减小；同时，偏心率参数也在发生变化。

3.2 轨道形状调整

上文介绍了平衡飞行的控制策略，这是个变速圆周运动的过程。航天器高度不变，瞬时轨道的半长轴持续变化，同时瞬时轨道偏心率在持续变化。经仿真分析可知，轨道偏心率调整策略如图3所示。

图3偏心率调整策略

Fig.3Strategies for eccentricity adjustment

需要说明的是，为保证在控制过程高度（或地心距）不变，控制开始的点应在瞬时轨道的近地点或远地点，且整个控制过程中，航天器总位于该时刻瞬时轨道的近地点或远地点，瞬时轨道的拱线随航天器控制机动过程转动。

4 仿真验证

为验证本文所提方法的有效性，在此以天琴计划卫星的轨道捕获控制为例开展仿真分析。

4.1 仿真初值设置

天琴卫星轨道^[8]半长轴取10⁵ km，轨道倾角为98.19°、初始近地点辐角为0°、初始升交点赤经为0°、初始真近点角为0°，历元选取2035年1月1日0点0分，卫星面质比设为0.02 m²/kg，地球重力场使用EGM2008引力模型，大气阻力使用Jacchia-Robert模型，太阳、月球引力根据星历生成，太阳光压常数设为1；卫星发动机提供连续推力的大小为1000 μN，推力方向可调；天琴卫星的初始质量取200 kg。

根据式（21）计算可得，在以上初始仿真参数的条件下，初始椭圆轨道半长轴在[10⁵-12.542 3，10⁵+12.545 4]km范围内时，可以在控制结束时刻将轨道形状调整为圆形。本文取初始椭圆轨道半长轴分别为（10⁵-12.5）km和（10⁵+12.5）km进行仿真验证。

4.2 仿真结果

4.2.1 由椭圆轨道转移到圆轨道

初始椭圆轨道半长轴分别为（10⁵-12.5）km和（10⁵+12.5）km，由式（22）可计算出控制律，即整个控制过程中推力加速度倾角的变化，由半长轴较小和较大的椭圆轨道转到圆轨道的推力倾角如图4所示。

图4控制律

Fig.4Control law

从较小椭圆开始轨道转移，控制时长为37 282 s（约10.3 h）；由较大椭圆轨道转移，控制时长为37 193 s（约10.3 h）。控制过程结束，航天器理论上在此时刻进入标称圆轨道，即完成了轨道捕获任务。

航天器角速度变化的解析解和数值解对比如图5所示，结果显示解析方法所求结果与数值验证吻合。其中数值积分参数设置如下：仿真步长为0.1 s，采取4阶龙格-库塔法进行数值积分，相对误差容限为10^-13，绝对误差容限为10^-13，变量的初值设置与上述解析仿真初值相同。

由半长轴较小和较大的椭圆轨道转到圆轨道的半长轴变化如图6所示，偏心率变化如图7所示。

图5角速度变化解析解与数值解对比

Fig.5Comparison of analytical and numerical solutions of angular velocity

图6半长轴变化

Fig.6Changes of semi-major axis

图7偏心率变化

Fig.7Changes of eccentricity

从图6、图7可看出，发动机关机时，长半轴达到标称值，而轨道偏心率已降为0，卫星进入目标圆轨道。

航天器高度变化如图8所示，仿真结果显示航天器高度几乎不变，与“平衡飞行”条件相符合。

上述是在CPU为AMD Ryzen 5 5600U、RAM为16 GB、操作系统为Windows10的计算机上使用MATLAB软件开展的仿真，应用平衡飞行理论的解析计算，由大椭圆到圆轨道和小椭圆到圆轨道的两次仿真时长分别为0.154 491 s与0.220 807 s，对比传统多冲量轨道机动控制律优化需要动辄几分钟甚至几小时的大规模优化计算相比，解析控制律大幅度缩短了计算时间，并且连续推力控制对比多冲量控制方法提升了控制精度。

图8高度变化

Fig.8Changes of height

4.2.2 由椭圆轨道经多次控制转移到圆轨道

初始椭圆轨道半长轴为99 985 km，初始偏心率为5×10^-5，轨道倾角为98.19°、近地点辐角为0°、升交点赤经为0°、真近点角为0°，此时远地点地心距为99 990 km，近地点地心距为99 980 km；目标轨道是半径为10⁵ km的圆轨道，其余轨道根数基本一致。

根据平衡飞行控制理论，轨道转移到少需要两次控制，第一次控制从轨道的远地点开始，近地点地心距从99 980 km抬升至10⁵ km，第一次轨道转移后，轨道半长轴为9.999 5 km，近地点地心距为99 990 km，远地点地心距为10⁵ km，第二次控制从转移后轨道的远地点开始，再次将近地点地心距从99 990 km抬升到10⁵ km后，完成轨道转移任务。由式（11）、式（13）计算，第一次控制时间为20 539 s，第二次控制时间为10 268 s。

在两段控制过程中，航天器的控制律如图9所示，半长轴变化如图10所示，偏心率变化如图11所示。在两段控制中，推力倾角均沿着时间下降，第二段控制后推力倾角降为0代表最终轨道为圆形。这里偏心率经历了一个阶梯式的提升过程，即两次提升到达目标；而偏心率第一次凑巧基本保持不变，在经过第二次控制时降到0，从而由椭圆轨道变成高近圆轨道。

图9轨道捕获过程控制律

Fig.9Control law in orbit acquisition

图10轨道捕获过程半长轴变化

Fig.10Changes of semi-major axis in orbit acquisition

图11轨道捕获过程偏心率变化

Fig.11Changes of eccentricity in orbit acquisition

由图10、图11可以看出，经过两段轨道转移后航天器轨道半长轴提升至10⁵ km，偏心率降低至接近0，完成任务指标，控制律的解析计算仿真总共用时0.285 08 s，发动机开机两次，关闭两次。当然在推力较大时，有一定裕量的条件之下，这两步控制过程还可以进一步优化，寻找出全局最优解。

5 结论

本文针对连续小推力机动条件下的地心轨道航天器的轨道捕获问题进行了具体研究，提出了考虑径向力的平衡飞行轨道控制方法，从航天器角速度调整、轨道形状调整、轨道高度调整三个方面展开分析，进一步提出了特殊飞行条件下轨道捕获问题的控制律，并开展了仿真分析，结果表明：

1）本文所提出的轨道控制策略可以精确实现轨道形状调整、航天器角速度及高度调整等目标；而且如果需要实现的轨道参数调整范围相对于推力过大，则可以使用多次、多圈控制方式灵活进行调整。

2）解析计算控制律极大减小了计算量和仿真时间，避免了传统优化算法复杂的迭代和初值猜测，且控制精度极高，可以达到引力波探测任务的精度要求。

3）本文所提的平衡飞行轨道控制方法给出了特殊情况下小推力轨道飞行的解析解，为小推力机动轨道的设计研究提供了新的思路，具有启发意义。

单星轨道的精确捕获可以为后续星座构形的精确捕获研究提供基础，下一步将针对空间引力波探测器星座构形的高精度捕获与保持控制问题继续开展研究。

图1天琴三星星座示意图

Fig.1Schematic diagram of TianQin 3 satellites constellation

下载: 全尺寸图片

图2轨道与极坐标示意图

Fig.2Illustration of orbit and polar coordinates

下载: 全尺寸图片

图3偏心率调整策略

Fig.3Strategies for eccentricity adjustment

下载: 全尺寸图片

图4控制律

Fig.4Control law

下载: 全尺寸图片

图5角速度变化解析解与数值解对比

Fig.5Comparison of analytical and numerical solutions of angular velocity

下载: 全尺寸图片

图6半长轴变化

Fig.6Changes of semi-major axis

下载: 全尺寸图片

图7偏心率变化

Fig.7Changes of eccentricity

下载: 全尺寸图片

图8高度变化

Fig.8Changes of height

下载: 全尺寸图片

图9轨道捕获过程控制律

Fig.9Control law in orbit acquisition

下载: 全尺寸图片

图10轨道捕获过程半长轴变化

Fig.10Changes of semi-major axis in orbit acquisition

下载: 全尺寸图片

图11轨道捕获过程偏心率变化

Fig.11Changes of eccentricity in orbit acquisition

下载: 全尺寸图片

表1表征等边三角形稳定性的指标要求

下载: 全尺寸图片

图1天琴三星星座示意图

Fig.1Schematic diagram of TianQin 3 satellites constellation

图2轨道与极坐标示意图

Fig.2Illustration of orbit and polar coordinates

图3偏心率调整策略

Fig.3Strategies for eccentricity adjustment

图4控制律

Fig.4Control law

图5角速度变化解析解与数值解对比

Fig.5Comparison of analytical and numerical solutions of angular velocity

图6半长轴变化

Fig.6Changes of semi-major axis

图7偏心率变化

Fig.7Changes of eccentricity

图8高度变化

Fig.8Changes of height

图9轨道捕获过程控制律

Fig.9Control law in orbit acquisition

图10轨道捕获过程半长轴变化

Fig.10Changes of semi-major axis in orbit acquisition

图11轨道捕获过程偏心率变化

Fig.11Changes of eccentricity in orbit acquisition

表1表征等边三角形稳定性的指标要求

图1天琴三星星座示意图

Fig.1Schematic diagram of TianQin 3 satellites constellation

图2轨道与极坐标示意图

Fig.2Illustration of orbit and polar coordinates

图3偏心率调整策略

Fig.3Strategies for eccentricity adjustment

图4控制律

Fig.4Control law

图5角速度变化解析解与数值解对比

Fig.5Comparison of analytical and numerical solutions of angular velocity

图6半长轴变化

Fig.6Changes of semi-major axis

图7偏心率变化

Fig.7Changes of eccentricity

图8高度变化

Fig.8Changes of height

图9轨道捕获过程控制律

Fig.9Control law in orbit acquisition

图10轨道捕获过程半长轴变化

Fig.10Changes of semi-major axis in orbit acquisition

图11轨道捕获过程偏心率变化

Fig.11Changes of eccentricity in orbit acquisition

表1表征等边三角形稳定性的指标要求

DECHER R, RANDALL J L, BENDER P L,et al. Design aspects of a laser gravitational wave detector in space[C]//Proceedings of Seminar on Active Optical Devices and Applications,1980.

Gravitational-Wave Community Science Team, Gravitational-Wave Core Team, Gravitational-Wave Science Task Force. Gravitational-wave mission concept study final report[R/OL].(2012-07-20)[2022-08-01].https://pcos.gsfc.nasa.gov/physpag/meetings/PhysPAG_Aug2012/presentations/GW_Study_Rev3_Aug2012-Final.pdf.

ISOYAMA S, NAKANO H, NAKAMURA T. Multiband gravitational-wave astronomy:observing binary inspirals with a decihertz detector, B-DECIGO[J]. Progress of Theoretical and Experimental Physics,2018,2018(7):073E01.

LISA Study Team. Laser interferometer space antenna for the detection and observation of gravitational waves[R/OL].[2022-08-01].https://www.lisa.nasa.gov/archive2011/Documentation/ppa2.08.pdf.

AMARO-SEOANE P, AOUDIA S, BABAK S,et al.eLISA:astrophysics and cosmology in the millihertz regime[EB/OL].(2012-01-17)[2022-08-01].https://arxiv.org/pdf/1201.3621.pdf.

胡一鸣, 梅健伟, 罗俊. 天琴计划与国际合作[J]. 科学通报,2019,64(24):2475-2483. HU Y M, MEI J W, LUO J. TianQin project and international collaboration[J]. Chinese Science Bulletin,2019,64(24):2475-2483.(in Chinese)

刘志国, 朴云松, 乔从丰. 多波段引力波宇宙研究和空间太极计划[J]. 现代物理知识,2016,28(5):28-33. LIU Z G, PIAO Y S, QIAO C F. Muti-band gravitational wave cosmic research and space Taiji program[J]. Modern Physics,2016,28(5):28-33.(in Chinese)

张雪峰, 叶伯兵, 檀庄斌, 等. 天琴轨道与星座设计进展[J]. 中山大学学报(自然科学版),2021,60(1/2):123-128. ZHANG X F, YE B B, TAN Z B,et al. Orbit and constellation design for TianQin:progress review[J]. Acta Scientiarum Naturalium Universitatis SunYatseni,2021,60(1/2):123-128.(in Chinese)

LEVCHENKO I, BAZAKA K, DING Y J,et al. Space micropropulsion systems for Cubesats and small satellites:from proximate targets to furthermost frontiers[J]. Applied Physics Reviews,2018,5(1):011104.

WANG J X, BAOYIN H X, LI J F,et al. Optimal four-impulse rendezvous between coplanar elliptical orbits[J]. Science China Physics, Mechanics and Astronomy,2011,54(4):792-802.

周军, 常燕. 考虑地球扁率J2摄动影响的异面椭圆轨道多冲量最优交会[J]. 宇航学报,2008,29(2):472-475. ZHOU J, CHANG Y. Optimal multiple-impulse rendezvous between non-coplanar elliptic orbits considering the J2 perturbation effects[J]. Journal of Astronautics,2008,29(2):472-475.(in Chinese)

CONWAY B A. Spacecraft trajectory optimization[M]. Cambridge: Cambridge University Press,2010.

BETTS J T. Survey of numerical methods for trajectory optimization[J]. Journal of Guidance Control Dynamics,1998,21(2):193-207.

YANG H W, BAI X L, LI S. Artificial equilibrium points near irregular-shaped asteroids with continuous thrust[J]. Journal of Guidance, Control,and Dynamics,2018,41(6):1308-1319.

EL MABSOUT B, CLAUDÉ P, DUDECK M. Rendezvous of a continuous low-thrust cubesat with a satellite[J]. Journal of Physics Communications,2019,3(5):055014.

KUETTEL Ⅲ D H, MCMAHON J W. Autonomous maneuver targeting around small bodies using continuous-thrust propulsion[J]. Journal of Guidance Control Dynamics,2022,45(3):499-516.

DACHWALD B, OHNDORF A. Global optimization of continuous-thrust trajectories using evolutionary neurocontrol[M]//Modeling and Optimization in Space Engineering. Switzerland: Springer,2019:33-57.

MAESTRINI M, DE VITTORI A, DI LIZIA P,et al. Dynamics-based uncertainty propagation with low-thrust[C]//Proceedings of AAS/AIAA Astrodynamics Specialist Conference,2022.

IZZO D,ÖZTÜRK E. Real-time guidance for low-thrust transfers using deep neural networks[J]. Journal of Guidance Control Dynamics,2021,44(2):315-327.

ZAVOLI A, FEDERICI L. Reinforcement learning for robust trajectory design of interplanetary missions[J]. Journal of Guidance Control Dynamics,2021,44(8):1440-1453.

YANG G, KLUEVER C A. An algorithm for computing near optimal,many revolution earth-orbit transfer[C]//Proceedings of Advances in the Astronautical Sciences,2006:1861-1880.

BETTS J T, ERB S O. Optimal low thrust trajectories to the moon[J]. SIAM Journal on Applied Dynamical Systems,2003,2(2):144-170.

PONTANI M, PUSTORINO M. Nonlinear earth orbit control using low-thrust propulsion[J]. Acta Astronautica,2021,179:296-310.

PONTANI M, PUSTORINO M. Low-thrust lunar capture leveraging nonlinear orbit control[J]. The Journal of the Astronautical Sciences,2023,70(5):28.

NAH R S, VADALI S R, BRADEN E. Fuel-optimal,low-thrust,three-dimensional earth-mars trajectories[J]. Journal of Guidance Control Dynamics,2001,24(6):1100-1107.

NAH R S. Fuel-optimal Earth-Mars trajectories using low-thrust exhaust-modulated plasma propulsion[D]. Texas: Texas A & M University,1999.

HERNANDO-AYUSO J, BOMBARDELLI C. Low-thrust collision avoidance in circular orbits[J]. Journal of Guidance Control Dynamics,2021,44(5):983-995.

理查德·H. 巴廷. 航天动力学的数学方法[M]. 倪彦硕, 蒋方华, 李俊峰, 译. 北京: 中国宇航出版社,2018. RICHARD H BATTIN. An introduction to the mathematics and methods of astrodynamics[M]. Translated by NI Y S, JIANG F H, LI J F. Beijing: China Aerospace Press,2018.(in Chinese)

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

Email订阅

Rss

English

1 轨道捕获任务建模

2 动力学方程的解析求解

3 平衡飞行轨道控制方法的应用

4 仿真验证

5 结论