高精度的轨道机动决策是空间交会、在轨监视与捕获等空间任务的关键技术^[1-4]。空间中一点到另一点固定时间的机动是轨道机动的基础问题。在二体模型条件下，该问题可通过Gauss法^[5]、Bate法^[5]等快速求解，实际在轨机动要考虑摄动因素，直接应用二体解会产生较大的末端位置偏差，通过迭代优化的方式可以获得考虑摄动影响的高精度解^[6-10]。然而，迭代优化运算量大，不利于在轨自主实时解算。
近年来，机器学习不断发展，在实时控制方面展现出优势^[11-13]。Sánchez-Sánchez和Izzo^[14]的研究表明深度学习技术可以应用于航天器导航和控制的星上决策。在行星际任务^[15]和月球着陆^[16]方面，神经网络也得到应用，通过优化等间接方法获得精确解，训练网络用于实时决策。
对于高精度在轨机动问题，直接对问题的精确解进行训练得到的神经网络决策解精度不够。考虑将二体运动规律滤除，只对摄动引起的复杂规律进行训练学习。据此，本文提出了一种基于机器学习的在轨实时轨道机动决策方法，使得决策精度与优化算法精度一致，且耗时大大缩短。
1 轨道机动问题 本文的轨道机动问题可描述为：给定卫星初始时刻的位置r₁和速度v₁, 要求经过Δt飞行时间，卫星达到终端时刻的位置r₂和速度v₂, 确定出卫星在初始时刻所需的速度增量Δv₁。
在航天工程应用中，往往会根据复杂约束确定标称工况，对于上述问题的标称工况，考虑摄动的影响，上述问题通常采用迭代优化的方法求解^[6]。
已知卫星的运动满足如下微分方程^[17]：

(1)

式中: μ为地心引力常数；F_p为各项摄动力，包括地球非球型引力、日月引力、太阳光压、大气阻力等^[17]；m为卫星质量。
在初始时刻，对卫星施加速度增量Δv₁，则上述问题的边界条件为

(2)

卫星飞行Δt时间后到达位置r(Δt)，r(Δt)可通过求解微分方程式(1)得到，则终端位置偏差为

(3)

于是，Δv₁可通过求解下述优化问题得到：Find: Δv₁

2 机器学习方法 在实际工程应用中，由于不确定性因素的存在，会在标称工况周围形成一个可行解域，即实际初始位置速度和终端位置速度会在一定范围内变动。在该范围内采样并经优化算法得到每个样本点的解，训练神经网络，实现该范围内任一点的轨道机动问题的实时高精度求解。
2.1 初始采样 给定样本中心点的输入向量为

给定样本点输入向量到样本中心点输入向量各维度的最远距离δr_1x, δr_1y, δr_1z, δv_1x, δv_1y, δv_1z, δr_2x, δr_2y, δr_2z, δv_2x, δv_2y, δv_2z，可以得到初始位置速度和终端位置速度各分量的变动范围，即样本输入向量的取值范围X_min~X_max。


在上述范围中，对于每个维度，按照均匀分布采样一次得到一个样本输入向量X=(X₁, X₂, X₃, X₄, X₅, X₆, X₇, X₈, X₉, X₁₀, X₁₁, X₁₂)^T。其中，X_i~U((X_min)_i, (X_max)_i)(i=1, 2, …, 12)，X_i表示X第i维度取值，(X_min)_i表示X_min第i维度取值，(X_max)_i表示X_max第i维度取值。
给定样本个数，按照上述方式完成采样。给定飞行时间Δt，对每一个样本点输入向量通过优化算法得到样本输出向量，即摄动精确解Δv₁。
2.2 摄动解与二体解作差 在二体模型条件下，通过Bate法可快速求得二体解Δv_1TB。
首先，计算转移角度Δθ:

(4)

然后，计算常数A:

(5)

之后，通过牛顿迭代法求解飞行时间Δt的普适方程：

(6)

式中: ；；。
解得z之后，计算拉格朗日系数f和g：

(7)

于是，得到二体模型下所需的速度增量Δv_1TB:

(8)

式中：f_Bate(·)为将Bate法求解速度增量的过程函数。
摄动解与二体解作差得到摄动项引起的速度增量差为

(9)

2.3 坐标变换 已知卫星初始时刻位置r₁和速度v₁, 则

(10)

(11)

(12)

(13)

式中: H为动量矩矢量；i_o、j_o和k_o分别为轨道坐标系x、y和z的单位矢量。
于是，地心赤道惯性坐标系到质心轨道坐标系的坐标变换矩阵为

(14)

进而得到速度增量差在轨道坐标系下的表示形式：

(15)

至此，得到最终的样本输出速度增量摄动修正项(ΔΔv₁)_o。
2.4 神经网络设计 设计摄动修正网络N_HpopCorrect，其输入为r₁, v₁, r₂, v₂在惯性坐标系的三轴分量形式，即X=(r_1x, r_1y, r_1z, v_1x, v_1y, v_1z, r_2x, r_2y, r_2z, v_2x, v_2y, v_2z)^T，输出为速度增量差在轨道坐标系的三轴分量形式，即Y=((ΔΔv_1x)_o, (ΔΔv_1y)_o, (ΔΔv_1z)_o)^T。
对网络包含的隐层层数、层内神经元节点数、每层的激活函数进行遍历，假设隐层层数有L_n个取值，层内神经元节点数有N_n个取值，激活函数有A_n个取值，则获得L_n×N_n×A_n个网络，对比获得性能最优的网络。网络性能指标采用均方误差(Mean Squared Error, MSE)。

(16)

式中：M为测试样本个数；D_OUT为样本输出维度，此处为3；y_{j, k}为第j个测试样本第k个维度的输出真值；为第j个测试样本第k个维度的网络预测输出值。
2.5 应用网络进行轨道机动决策 如图 1所示，轨道机动决策流程包括了如下步骤：
步骤1??，。
步骤2??。
步骤3??。
步骤4??。

图 1 轨道机动决策流程 Fig. 1 Decision-making process for orbital maneuvers

图选项

3 仿真验证 选定样本中心点为X₀=(5 846 591.286 6, -1 761 003.240 6, 3 479 801.993 5, 3 176.387 9, -2 216.637 1, -6 458.570 7, 5 432 239.472 1, -2 636 976.590 5, -3 602 112.081 7, -3 952.612 7, 604.677 3, -6 351.835 6)^T，选定飞行时间Δt为1 000 s。
选定样本点到样本中心点各维度的最远距离如表 1所示。
表 1 样本点到样本中心点最远距离 Table 1 The farthest distance from sample point to sample center point

参数	数值
δr1x/m	1 000
δr1y/m	1 000
δr1z/m	1 000
δv1x/(m·s-1)	1
δv1y/(m·s-1)	1
δv1z/(m·s-1)	1
δr2x/m	1 000
δr2y/m	1 000
δr2z/m	1 000
δv2x/(m·s-1)	1
δv2y/(m·s-1)	1
δv2z/(m·s-1)	1

表选项






于是，样本取值范围为X_min~X_max。


在上述范围X_min~X_max内按均匀分布(X_i~U((X_min)_i, (X_max)_i)i=1, 2, …, 12)随机采取6 000个样本点，每个样本点经优化获得摄动下的精确解Δv₁，作为样本的初始输出。摄动精确解与二体解作差，再经坐标变换得到最终输出。
为比较不同输出网络的性能，设计3种神经网络，分别为速度增量网络N_Dir、摄动偏差网络N_DeTB、摄动修正网络N_HpopCorrect。其中，速度增量网络以速度增量Δv₁为网络输出，摄动偏差网络以速度增量差ΔΔv₁为网络输出，摄动修正网络以速度增量摄动修正项(ΔΔv₁)_o作为网络输出。
对3种网络的结构参数分别进行优化设计。隐层层数和层内神经元节点数取值如表 2所示。
表 2 网络结构参数取值 Table 2 Parameter value of network structure

神经网络参数	数值
网络隐层数	1，2，3，4
隐层节点数	8，16，32，64，128

表选项






隐层激活函数选取最常见的几种，分别是logsig[0, 1]，softmax[0, 1]，poslin[0, 1]，purelin[-∞, ∞]，tansig[-1, 1]。输出层激活函数取purelin[-∞, ∞]。网络训练最大迭代次数取5 000次，梯度下降下限取1.0×10^-5，误差最大连续不下降次数取50次。
在6 000个样本数据中，取前5 400个样本进行网络训练，后600个样本进行测试。
1) 直接对速度增量进行训练。统计不同参数网络的MSE记录在表 3中。由表 3可以看出，直接对速度增量进行训练，最优的速度增量网络包含3层隐层，隐层有16个神经元节点，隐层激活函数为purelin。将600个测试样本输入到最优速度增量网络，得到对应的速度增量Δv_1NDir。
表 3 不同层数、节点数和激活函数的速度增量网络的性能 Table 3 Performance of speed increment networks with different layers, units and activation functions

隐层层数/节点数	MSE
	logsig	softmax	poslin	purelin	tansig
1/8	9.95×10-6	8.38×10-5	8.23×10-9	2.06×10-9	5.76×10-6
1/16	3.05×10-6	1.81×10-5	5.42×10-6	2.02×10-9	1.45×10-6
1/32	6.40×10-6	2.20×10-6	2.17×10-9	2.02×10-9	3.83×10-6
1/64	6.41×10-6	7.47×10-6	2.73×10-7	2.01×10-9	4.27×10-6
1/128	2.46×10-4	2.62×10-6	7.91×10-6	2.02×10-9	2.97×10-6
2/8	2.26×10-5	5.98×10-5	2.86×10-5	2.02×10-9	3.10×10-5
2/16	5.66×10-6	2.69×10-5	3.86×10-6	2.02×10-9	3.70×10-6
2/32	2.60×10-6	4.37×10-6	6.06×10-5	2.02×10-9	8.53×10-6
2/64	3.33×10-4	1.54×10-6	5.09×10-4	2.02×10-9	6.34×10-5
2/128	1.89×10-3	2.91×10-6	5.29×10-3	2.02×10-9	9.91×10-4
3/8	2.94×10-5	1.18×10-4	2.96×10-1	2.02×10-9	8.34×10-6
3/16	3.23×10-6	1.60×10-5	4.88×10-5	2.01×10-9	4.08×10-6
3/32	5.15×10-6	1.64×10-5	2.43×10-3	2.03×10-9	3.28×10-6
3/64	4.24×10-4	9.36×10-6	5.23×10-3	2.02×10-9	1.25×10-4
3/128	2.90×10-3	2.93×10-6	9.91×10-2	2.02×10-9	3.13×10-3
4/8	4.02×10-5	2.98×10-1	2.96×10-1	2.02×10-9	2.57×10-5
4/16	2.45×10-6	3.03×10-1	1.05×10-2	2.02×10-9	8.95×10-6
4/32	5.82×10-6	1.31×10-5	7.84×10-3	2.03×10-9	9.21×10-6
4/64	4.66×10-4	1.13×10-5	1.35×10-1	2.02×10-9	2.91×10-4
4/128	2.65×10-3	2.98×10-1	1.50	2.02×10-9	2.94×10-3

表选项






2) 对速度增量的摄动因素进行训练。统计不同参数网络的MSE记录在表 4中。由表 4可以看出，对速度增量的摄动因素进行训练，最优摄动偏差网络包含4层隐层，每层隐层有128个神经元节点，隐层激活函数为purelin。将600个测试样本输入到最优摄动偏差网络，得到对应的ΔΔv_1NDeTB, 加上对应的二体解Δv_1TB，得到对应的速度增量Δv_1NDeTB。
表 4 不同层数、节点数和激活函数的摄动偏差网络的性能 Table 4 Performance of perturbation deviation networks with different layers, units and activation functions

隐层层数/节点数	MSE
	logsig	softmax	poslin	purelin	tansig
1/8	2.80×10-6	1.04×10-6	1.91×10-6	6.11×10-6	5.50×10-6
1/16	1.62×10-6	1.04×10-6	2.24×10-6	5.87×10-6	4.11×10-6
1/32	3.60×10-6	1.34×10-6	8.77×10-6	1.87×10-6	4.86×10-6
1/64	3.27×10-6	1.19×10-6	7.61×10-6	2.58×10-6	3.04×10-6
1/128	8.65×10-6	2.26×10-6	8.98×10-6	7.81×10-7	7.07×10-6
2/8	2.08×10-6	1.30×10-6	3.33×10-6	3.96×10-6	4.73×10-6
2/16	3.13×10-6	2.07×10-6	4.83×10-6	8.44×10-7	4.59×10-6
2/32	3.20×10-6	1.78×10-6	4.35×10-6	3.59×10-7	3.19×10-6
2/64	1.86×10-6	8.88×10-7	5.48×10-6	1.09×10-7	2.16×10-6
2/128	5.23×10-6	5.71×10-7	3.04×10-6	1.52×10-8	2.02×10-6
3/8	6.12×10-6	5.85×10-7	1.61×10-6	9.80×10-7	2.06×10-6
3/16	5.45×10-6	8.64×10-7	2.22×10-6	4.71×10-7	2.32×10-6
3/32	2.44×10-6	7.58×10-7	3.87×10-6	8.09×10-8	2.74×10-6
3/64	2.49×10-6	2.02×10-6	3.55×10-6	2.98×10-9	2.52×10-6
3/128	2.43×10-6	7.33×10-7	9.31×10-7	4.11×10-10	2.98×10-6
4/8	7.63×10-7	1.52×10-6	9.29×10-7	1.03×10-6	9.74×10-7
4/16	9.66×10-7	8.19×10-7	3.84×10-6	4.16×10-7	2.02×10-6
4/32	1.99×10-6	1.24×10-6	3.53×10-6	6.47×10-9	1.12×10-6
4/64	1.84×10-6	1.51×10-6	3.57×10-7	1.06×10-10	1.57×10-6
4/128	1.87×10-6	3.99×10-7	6.99×10-7	8.84×10-12	1.18×10-6

表选项






3) 对速度增量摄动因素坐标变换后进行训练。统计不同参数网络的MSE记录在表 5中。由表 5知，对速度增量摄动因素坐标变换后进行训练，最优的摄动修正网络包含4层隐层，每层隐层有128个神经元节点，隐层激活函数为purelin。将600个测试样本输入到最优摄动修正网络，得到每个样本的对应输出(ΔΔv₁)_o，经坐标变换得到对应的ΔΔv_{1NHpspCorrect}, 加上对应的二体解Δv_1TB，得到对应的速度增量Δv_{1NHpspCorrect}。
表 5 不同层数, 节点数和激活函数的摄动修正网络的性能 Table 5 Performance of perturbation correction networks with different layers, units and activation functions

隐层层数/节点数	MSE
	logsig	softmax	poslin	purelin	tansig
1/8	9.44×10-7	3.87×10-7	3.16×10-6	2.70×10-6	5.14×10-6
1/16	2.29×10-6	9.71×10-7	6.09×10-6	3.01×10-6	1.83×10-6
1/32	1.64×10-6	1.36×10-6	7.63×10-6	3.61×10-6	4.02×10-6
1/64	2.87×10-6	7.55×10-7	4.58×10-6	3.67×10-6	9.22×10-6
1/128	7.44×10-6	2.81×10-7	6.37×10-6	5.56×10-7	7.38×10-6
2/8	1.24×10-6	2.06×10-6	2.71×10-6	2.29×10-6	1.21×10-6
2/16	9.09×10-7	6.05×10-7	4.78×10-6	2.10×10-6	1.84×10-6
2/32	2.73×10-6	1.26×10-6	1.99×10-6	1.02×10-6	4.42×10-6
2/64	1.89×10-6	1.08×10-6	5.50×10-6	6.50×10-8	2.95×10-6
2/128	4.18×10-6	9.30×10-7	5.67×10-6	2.20×10-8	3.87×10-6
3/8	3.00×10-7	2.51×10-6	1.52×10-6	1.51×10-6	1.16×10-6
3/16	2.18×10-6	1.81×10-6	4.38×10-6	1.00×10-6	1.05×10-6
3/32	1.65×10-6	1.10×10-6	5.48×10-6	1.01×10-7	2.51×10-6
3/64	2.15×10-6	1.40×10-6	4.15×10-6	4.33×10-9	1.95×10-6
3/128	2.95×10-6	9.17×10-7	2.87×10-6	5.72×10-10	2.33×10-6
4/8	4.20×10-6	7.73×10-7	1.46×10-6	1.55×10-6	1.26×10-6
4/16	7.32×10-7	1.25×10-6	2.00×10-6	3.41×10-7	1.65×10-6
4/32	8.72×10-7	1.00×10-6	1.85×10-6	5.47×10-9	2.16×10-6
4/64	2.11×10-6	1.54×10-6	1.93×10-6	1.14×10-10	1.54×10-6
4/128	2.63×10-6	1.00×10-6	3.29×10-7	7.28×10-12	2.03×10-6

表选项






将3种最优网络得到的速度增量施加到初始时刻进行仿真，计算Δt时刻后卫星的终端位置偏差d_t。统计600个样本点对应d_t的最大值、最小值、平均值，分别记为d_max，d_min，d_mean，结果记录在表 6中。
表 6 终端位置偏差统计 Table 6 Statistics of terminal position deviation

方法	dmax/m	dmin/m	dmean/m
优化算法(样本)	3.864 0×10-3	2.340 0×10-4	2.008 0×10-3
最优速度增量网络	2.892 5×10-1	3.076 0×10-3	7.396 2×10-2
最优摄动偏差网络	1.192 3×10-2	6.240 0×10-4	4.751 9×10-3
最优摄动修正网络	1.117 8×10-2	3.150 0×10-4	4.175 8×10-3

表选项






分析表 6数据，可以得出：
1) 对比最优速度增量网络和最优摄动偏差网络的结果，不管是终端位置偏差的最大值、最小值还是平均值，最优摄动偏差网络的效果都要明显优于最优速度增量网络的效果(1个数量级)，即对速度增量摄动项因素进行训练的效果要明显优于直接对速度增量训练的效果。
2) 对比最优摄动偏差网络和最优摄动修正网络结果，不管是终端位置偏差的最大值、最小值还是平均值，最优摄动修正网络的效果都要略微优于最优摄动偏差网络的效果，即将速度增量的摄动因素进行坐标变换，投影到轨道坐标系下，训练效果更好一点。
3) 最优摄动修正网络与优化算法的结果很接近，即对速度增量摄动因素坐标变换后进行训练出的神经网络在精度方面几乎与优化算法精度一致。
仿真所用计算机性能为：内存8 GB；处理器Inter Core(TM) i7-2600@2.70 GHz；GPU NVIDIA GeForce GTX1050。优化算法给出一次决策的时间大约为200 s，综合应用最优摄动修正网络和二体解给出一次决策的时间为0.02 s左右，相比优化算法，后者耗时大大缩短，实现了精度与耗时的统一，适于星上在轨实时决策。
4 结束语 对于固定两点间固定时间的轨道机动问题，在初始位置速度和终端位置速度一定误差范围内，本文得到的最优摄动修正网络决策精度与优化算法的精度一致，但其决策耗时仅为优化算法决策耗时的0.01%，实时性高。
本文为在轨机动、行星际转移等复杂问题的快速高精度求解提供了一种思路，将精确解分解为简化模型解与残差的和，对残差进行训练得到神经网络，再与简化模型解组合获得原问题的高精度解。

参考文献

[1]	夏红伟, 李莉, 曲耀斌, 等. 卫星编队构型设计与轨道机动算法优化[J]. 中国惯性技术学报, 2013, 21(2): 186-191. XIA H W, LI L, QU Y B, et al. Satellites formation configuration design and orbit maneuver algorithm optimization[J]. Journal of Chinese Inertial Technology, 2013, 21(2): 186-191. DOI:10.3969/j.issn.1005-6734.2013.02.012 (in Chinese)
[2]	于瀚. 航天器轨道机动可达区域研究[D]. 哈尔滨: 哈尔滨工业大学, 2017. YU H. The study of reachable domain for spacecraft maneuver[D]. Harbin: Harbin Institute of Technology, 2017(in Chinese).
[3]	谭丽芬. 赤道椭圆交会轨道规划与制导方法[D]. 长沙: 国防科技大学, 2011. TAN L F. Rendezvous trajectory planning and guidance approach for equatorial elliptical orbit[D]. Changsha: National University of Defense Technology, 2011(in Chinese).
[4]	张守玉, 姜振东. 基于STK的卫星轨道机动模型设计与仿真[J]. 计算机仿真, 2004, 21(10): 25-27. ZHANG S Y, JIANG Z D. Design and simulation of satellite orbital maneuver model on STK[J]. Computer Integrated Manufacturing Systems, 2004, 21(10): 25-27. (in Chinese)
[5]	BATE R R. 航天动力学基础[M]. 吴鹤鸣, 李肇杰, 译. 北京: 北京航空航天大学出版社, 1990. BATE R R. Fundamentals of astrodynamics[M]. WU H M, LI Z J, translated. Beijing: Beihang University Press, 1990(in Chinese).
[6]	李栋林, 黄福铭. 基于Lambert问题的精确拦截与交会策略研究[J]. 飞行力学, 2008, 26(2): 57-59. LI D L, HUANG F M. Research into accurate interception and rendezvous scheme based on Lambert problem[J]. Flight Dynamics, 2008, 26(2): 57-59. (in Chinese)
[7]	桑艳, 周进. 基于Lambert算法的脉冲精确变轨策略[J]. 国防科技大学学报, 2009, 31(3): 29-32. SANG Y, ZHOU J. An approach of accurate impulse transfer based on Lambert algorithm[J]. Journal of National University of Defense Technology, 2009, 31(3): 29-32. (in Chinese)
[8]	JEZEWSKI D J. Optimal rendezvous trajectories subject to arbitrary perturbations and constraints[C]//AIAA/AAS Astrodynamics Specialist Confernce. Reston: AIAA Press, 1992: 4507.
[9]	CHANG Y, ZHOU J. Orbital correction method for two-impulse rendezvous between non-coplanner elliptic orbits considering the J2 perturbation[J]. Journal of Astronautics, 2008, 29(4): 1172-1176.
[10]	周须峰, 唐硕. 固定时间拦截变轨段制导的摄动修正方法[J]. 飞行力学, 2006, 24(4): 46-49. ZHOU X F, TANG S. Disturbed modify method of fixed-time interception's guidance in orbit-change stage[J]. Flight Dynamics, 2006, 24(4): 46-49. (in Chinese)
[11]	裴忠才, 尹丽, 王占林. 基于神经网络的仿真转台控制系统[J]. 北京航空航天大学学报, 2005, 31(9): 1045-1048. PEI Z C, YIN L, WANG Z L. Simulating turntable control system with neural network[J]. Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(9): 1045-1048. DOI:10.3969/j.issn.1001-5965.2005.09.024 (in Chinese)
[12]	CHENG L, WANG Z, JIANG F, et al. Real-time optimal control for spacecraft orbit transfer via multiscale deep neural networks[J]. IEEE Transactions on Aerospace and Electronic Systems, 2018, 55(5): 2436-2450.
[13]	ZHONG R, XU S. Neural-network-based terminal sliding-mode control for thrust regulation of a tethered space-tug[J]. Astrodynamics, 2018, 2(10): 175-185. DOI:10.1007/s42064-017-0019-0
[14]	SáNCHEZ-SáNCHEZ C, IZZO D. Real-time optimal control via deep neural networks: Study on landing problems[J]. Journal of Guidance, Control, and Dynamics, 2018, 41(5): 1122-1135.
[15]	IZZO D, SPRAGUE C I, TAILOR D V. Machine learning and evolutionary techniques in interplanetary trajectory design[M]//FASANO G, PINIER J. Modeling and optimization in space engineering. Berlin: Springer, 2019: 191-210.
[16]	FURFARO R, BLOISE I, ORLANDELLI M, et al. Deep learning for autonomous lunar landing[C]//AIAA/AAS Astrodynamics Specialist Conference. Reston: AIAA, 2018: 18-363.
[17]	董云峰, 陈士明, 苏建敏, 等. 卫星姿态控制动态模拟技术[M]. 北京: 科学出版社, 2010: 300-301. DONG Y F, CHEN S M, SU J M, et al. Dynamic simulation technology of satellite attitude control[M]. Beijing: Science Press, 2010: 300-301. (in Chinese)