摘要
针对单追踪星在轨服务多目标星的高精度交会轨迹优化问题，提出一种基于多轮迭代修正的三步串行优化方法。采用差分进化算法获得基于J2摄动解析模型的交会序列与机动脉冲的近似解，通过序列二次规划算法进行局部优化得到第一次交会机动脉冲的高精度解，构建多轮迭代修正模型更新后续目标星状态及机动脉冲，以减小不同模型间误差的扩散影响。仿真结果表明，提出的方法能够以较少的优化迭代次数，较稳定得到机动成本近似最优的高精度多星交会机动参数。
Abstract
For the high-precision rendezvous trajectory optimization problem of a single chase satellite serving multiple target satellites on orbit, a three-step serial optimization method based on the multi-round iterative correction was proposed. The rendezvous sequence and approximate impulse solutions based on the analytical J2 perturbation model were obtained by using the differential evolution algorithm, and the high-precision solution of the first rendezvous maneuvering impulses was obtained by using the local optimization of the sequential quadratic programming algorithm. Moreover, a multi-round iterative correction model was constructed to update the subsequent target satellites′ states and maneuvering impulses to reduce the influence of the error diffusion between different models. Simulation results show that the proposed method can stably obtain the high-precision maneuvering parameters with approximate least cost in a few iterations.
航天技术飞速发展，在轨航天器数量不断增多，空间任务规模和复杂度也日渐增加[1-2]。在一次飞行任务中访问多个空间目标，可以提高任务收益，降低任务平均成本，这一思路在小行星探索中尤为明显且已经成功运用[3]。针对近地空间，在轨加注、碎片清除等在轨服务任务同样对一次任务访问多个目标提出了需求[4-6]。多星交会轨迹优化是实现上述任务场景的关键技术[7]。
多星交会轨迹优化问题是航天动力学领域的一个经典问题，诸多学者针对该问题开展了研究，大多可以从动力学模型和算法框架两个方面加以总结[8-9]。作为大规模空间目标访问的基础，对单转移段速度增量的高效度量能够提高全局优化过程中目标分配和序列优化的效率。Alfriend等[10]揭示了小倾角地球同步轨道（geostationary orbit，GEO）多星交会问题的转移速度增量与赤道面内角动量矢量投影点相对距离的正比关系，并将其转化为平面点集最短路径求解问题。Huang等[11]、Chen等[12]分别采用半解析和解析的方法基于轨道参数差值对J2摄动下两目标之间交会所需最优速度增量进行估算。Zhang等[13]对J2摄动条件下的轨道动力学模型进行线性化近似用于求解多脉冲交会方案。Yang等[14]提出了一种基于Lambert算法做同伦打靶的轨道两点边值问题求解方法，可用于求解摄动条件下航天器长时间转移交会问题。
对于多星交会轨迹优化问题所涉及的混合整数非线性优化问题，学者们通过对优化算法的适应性改进及应用得到了不错的优化结果。Bakhtiari等[15]针对二体轨道下的多星交会问题，采用改进的粒子群（particle swarm optimization，PSO）算法对转移脉冲进行优化。Dutta等[16]针对多星在轨服务问题，提出了等价分布式的交会策略，将其刻画为最小成本网络流问题，通过局部搜索算法进行最优求解。Yu等[17]针对GEO非共面在轨服务场景，同时以燃料成本与转移时间为指标，采用多目标粒子群算法进行优化。Zhu等[18]建立了基于服务星-燃料站协同在轨服务网络的多星交会模型，考虑了J2摄动对燃料成本的影响，通过一种基于密度的混合空间聚类算法进行多星访问优化。此外，Huang等[19]基于解析速度增量估算模型推导出目标函数相对于优化变量的偏导数以及约束表达式，并以此引导非线性规划问题的求解。由于多星交会问题的复杂性和层次性，许多学者采用了分层优化的方法进行模型的建立和求解。Zhang等[20]、Zhao等[21]以及Li等[22]在处理J2摄动下的多星交会问题时，都采用了分层优化策略，将序列问题和脉冲问题进行有效解耦，改善了全局优化效果。对于末端速度约束更加自由的多星飞越规划问题，彭晨远等[23]采取了两层优化方法得到不同精度模型下的脉冲解。
已有研究大多关注二体或J2解析摄动模型下的多星交会，还未看到基于高精度摄动轨道积分模型的多星交会研究。高精度轨道模型的计算量远大于二体及解析摄动模型，采用进化计算等优化算法直接基于高精度轨道模型进行优化往往计算效率过低。对单次交会常采用微分修正将二体及解析摄动模型改进为高精度结果，但存在迭代过程对初值敏感而难以收敛的问题。
本文提出一种高精度多星交会优化方法：基于J2摄动解析模型得到多星交会序列以及脉冲近似解，对近似解进行微分修正得到第一次交会的高精度解；采用多轮迭代修正的方法对后续交会进行逐段状态更新并重新规划，以得到满足交会精度要求及速度增量近最优的多星交会轨迹方案。
1 多星交会任务场景
单追踪星对多目标星的交会服务场景如图1所示，多颗目标星处于相邻但不相同的轨道面以及轨道位置（为方便展示，图1中刻画为同一圆轨道），各目标星的序号标记为p∈{1，2，3，···，N}，N为目标星总数量。
图1多星交会场景示意图
Fig.1Sketch map of multi-satellite rendezvous scenario
追踪星由运载火箭发送至交会序列第一颗目标星相同空间位置速度状态，与第一颗星交会过程的机动消耗在多星交会任务中不做考虑。任务开始时，追踪星从第1颗目标星所在位置出发，经不同转移轨道逐个与目标星交会，总共存在N-1个转移弧段。该过程分为三个主要步骤：首先，追踪星施加脉冲机动经转移轨道与目标星交会；其次，追踪星与目标星保持相同位置速度进行自由飞行，并在规定的时段完成在轨服务任务；最后，追踪星施加脉冲机动离开该目标星进入转移轨道前往下一个目标星。所有转移弧段按照上述步骤循环进行，直到完成N颗目标星的交会任务。为满足工程实际需求，追踪星与目标星的轨道计算均基于高精度摄动轨道积分模型。
对于第n∈{1， 2， 3， ···，N-1}次交会任务，将每个转移弧段的目标星序号、交会转移时间以及在轨服务时间分别设置为pn、Δttrans，n和Δtserv，n。对于单次交会的轨道转移，多脉冲转移较两脉冲转移往往成本更小[24]。在交会过程中采用多脉冲形式，每个转移弧段脉冲数量相同，脉冲序号为i∈{1，2，3，···，M}，M为单次交会脉冲数量，每次机动的脉冲时刻和脉冲矢量分别为tn，i、Δvn，i。对于该多星交会任务，其初始时刻设置为t0，追踪星初始空间状态为EC（t0），对于序号为p的目标星，其初始空间状态为，交会过程中时间和轨道状态关系如下所示。
1.1 时间关系
第n次交会任务的初始时刻和结束时刻分别为：
(1)
(2)
1.2 轨道状态关系
第n次交会任务的目标星初始和结束时刻状态分别为：
(3)
(4)
式（3）~（4）中G为轨道外推传播函数，考虑多种摄动项影响时，基于式（5）所示常微分方程的数值积分实现功能。
(5)
其中，r为惯性坐标系下的位置矢量，r为空间距离，μ为地球引力常数，aJ为地球非球形摄动，as为太阳引力摄动， am为月球引力摄动，ad为大气阻力摄动。
第n次交会任务中，追踪星的交会初始轨道状态如式（6）所示：
(6)
其中，第1次任务的初始轨道状态与第1颗目标星相同，这是因为该场景假设追踪星部署于第一颗目标星附近。
2 多星交会任务规划模型与方法
2.1 任务规划模型
第1节给出的多星交会问题可以看成是一个动态旅行商优化问题，下面对其设计变量、优化目标与约束条件进行建模。
2.1.1 设计变量
高精度多星交会任务轨迹规划最终得到完整多星交会最优方案，其优化变量包括两部分：首先是交会次序变量S，如{1，2，3，···，n!}，为整数变量；然后是交会过程精确脉冲时间tn，i以及交会过程精确脉冲矢量Δvn，i，为实数变量。
2.1.2 优化目标
多星交会优化目标为整个任务过程所有脉冲机动的速度增量总和最小：
(7)
2.1.3 约束条件
约束条件主要包括脉冲时间、脉冲矢量以及交会末端位置等约束。
1）交会时间等式约束。采用固定时间交会转移模式，交会转移时间及在轨服务时间之和需与交会间隔时间相等：
(8)
其中，Δtinterv为两次交会间隔时间，为尽可能通过轨道面漂移来减小升交点赤经差ΔΩ从而节省航天器轨道机动成本，下面对间隔时间做出规定。
首先，给出交会上限时间计算表达式，其定义为目标航天器飞行X圈次所对应的时间：
(9)
其中，aT为目标星轨道半长轴。
其次，追踪星与目标星升交点赤经差值完全由轨道面漂移消除所需时间决定，即
(10)
其中，下标“C”和“T”代表相关参数分别属于追踪星和目标星。
结合式（9）和式（10）给出两次交会间总间隔时间：
(11)
其中，对于|ΔΩ|不断扩大的情况，Δtinterv取，这样既能保证交会时间较充足，防止Lambert解不收敛或者出现速度增量过大的情况，又使得|ΔΩ|增加量受到限制。
2）在轨服务时间下限约束。考虑交会后在轨服务具体操作时间需求，设置时间下限：
(12)
3）单次脉冲值上限约束。受发动机性能影响，单次脉冲大小存在上限：
(13)
4）交会末端相对距离上限约束。高精度多星交会对于每个交会段末端，均设置最大交会距离误差：
(14)
其中，、分别为第n次交会时刻目标星和追踪星的空间位置矢量，ε为最大交会距离误差。
2.2 优化方法
2.2.1 三步串行优化方法
针对高精度摄动轨道积分过程计算量大，难以直接进行全局优化的问题，本文提出三步串行优化方法，如图2所示，对于第i次转移，有以下优化流程。
图2高精度交会三步串行优化方法
Fig.2Three-step serial optimization method for high precision rendezvous
步骤1：初始化所有目标星轨道状态矢量。
步骤2：利用差分进化算法基于半解析估算方法[10]对于当前各星状态进行多星交会序列规划问题的全局优化，得到J2摄动下速度增量估计量最小的最优交会序列，该步骤优化变量为交会序列S，优化目标函数为脉冲速度增量估计量总和，如式（15）所示：
(15)
其中，为每个转移弧段速度增量估计量。
步骤3：选择步骤2最优交会序列中前两颗星作为当前转移段追踪星与目标星，利用差分进化算法，基于J2摄动Lambert同伦打靶方法[13]，对得到的最优交会序列中第一个交会转移弧段，进行J2摄动条件下的转移弧段脉冲时刻及近似脉冲矢量的全局优化。优化变量为近似脉冲时刻，以及近似脉冲矢量，优化目标函数为脉冲近似速度增量总和：
(16)
步骤4：基于高精度轨道模型，利用序列二次规划算法对近似脉冲矢量进行局部微分修正，以消除末端交会误差，得到当前序列第一个交会转移弧段脉冲精确解。优化变量为交会过程脉冲变量tn，i、Δvn，i，目标函数为交会末端相对距离：
(17)
步骤5：若相对距离指标J3超出公差范围未能满足约束，则返回步骤3；若满足约束，则输出当前交会弧段脉冲精确解。
2.2.2 多轮迭代修正优化模型
上述优化方法能得到单转移弧段高精度交会方案，但直接应用于多星交会连续多个转移弧段优化时，存在模型误差扩散，难以收敛的问题。基于三步串行方法提出多轮迭代修正的高精度优化模型，具体流程如图3所示。
步骤1：输入追踪星以及未完成交会的目标星的初始位置速度。
步骤2：通过三步串行优化方法获得单次交会目标星对象及其高精度解。
步骤3：输出当前交会弧段高精度解，若还有未交会目标星，更新高精度轨道摄动条件下目标星轨道状态，转到步骤1重复上述步骤进行剩余交会弧段求解，直至依次求解出所有交会弧段高精度解。
该优化模型能够基于上一次高精度规划的结果对后续目标星进行更新，并在此基础上对剩余目标星交会近似解进行重规划，进而消除目标星在J2摄动模型和高精度摄动模型的累计偏差。在这种修正下，J2摄动近似解并不会过于偏离高精度解，能够改善高精度解迭代过程的收敛性。
图3多轮迭代修正优化模型
Fig.3Multiple rounds of iterative correction optimization model
3 算例分析
3.1 算例配置
考虑单颗追踪星与多颗目标星之间的序列交会，追踪星的起点状态为交会序列第一颗目标星的空间状态，设置算例Ⅰ为5颗邻近轨道的低轨卫星，具体参数（半长轴a、偏心率e、轨道倾角i、升交点赤经Ω以及纬度幅角u）如表1所示。设置算例Ⅱ为7颗较大轨道差异的低轨卫星，由于算例Ⅱ中目标星偏心率不为零，所以将纬度幅角u拆分为近地点角距ω以及真近点角f。任务开始的历元时刻为2022年11月3日00：00（UTC），交会上限时间采用的是第一颗目标星的30个轨道周期。
表1算例Ⅰ目标星参数
Tab.1Target parameters of example Ⅰ
摄动轨道积分采用6×6阶的EGM96地球非球形引力场模型，同时考虑大气阻力摄动与太阳和月球的第三体引力摄动。
3.2 单交会弧段仿真结果分析
在多星交会轨迹优化问题中，单交会弧段的高精度脉冲交会求解，是整体优化方法中的重要部分。以目标星S1到S2的交会弧段为例，在高精度摄动条件下，采用三脉冲Lambert打靶迭代再进行局部修正，得到脉冲交会高精度解。
图4和图5分别为脉冲近似解全局优化和高精度解局部修正迭代过程，两类优化都能在一定迭代次数内实现最优值良好收敛，最终获得高精度解：t=6 704.813 1 s、t=9 504.456 3 s、t=83 304.419 4 s时，脉冲量为Δv=149.590 m/s、Δv=234.056 m/s、Δv=53.811 m/s。该脉冲方案在J2000.0地心惯性坐标系下的三维仿真结果如图6所示。在末端时刻两星的空间位置几乎重合，两者末端轨道参数偏差为ΔE={-8.605×10-7 km，6.230×10-13，3.797×10-14（°），-7.772×10-15（°），5.002×10-12（°）。
图4近似解全局优化迭代过程
Fig.4Approximate solution global optimization iterative process
图5精确解局部修正迭代过程
Fig.5Accurate solution local correction iterative process
图6高精度交会轨迹仿真
Fig.6High precision rendezvous trajectory simulation
3.3 多星交会仿真结果分析
3.3.1 多星交会方案及仿真结果
算例Ⅰ多星交会优化结果如表2所示，所对应脉冲方案仿真结果部分轨道参数变化情况如图7所示。
表2算例Ⅰ多星交会优化结果
Tab.2Optimization result of multi-satellite rendezvous of example Ⅰ
图7从航天器半长轴、轨道倾角、升交点赤经三个轨道参数的角度对两航天器的交会过程进行展示，追踪星在每个转移过程中都分别施加了三次脉冲，并在终端时刻追踪星与各个目标星轨道参数在收敛标准内完全吻合，实现了与多颗目标星的高精度交会。
图7算例Ⅰ多星交会过程部分轨道参数变化情况
Fig.7Change of partial orbital parameter in multi-satellites rendezvous process of example Ⅰ
此外，在满足交会末端精度要求的前提下高精度多星交会优化目标是实现总速度增量，即机动成本最小。图8展示了不同轨道模型速度增量累积量的变化情况，其变化趋势基本一致，J2摄动解析解和高精度最优解之间差别较小。因此，在J2摄动解邻域范围内迭代计算高精度解具有合理性。该算例中三种速度增量最终累积量接近，估算速度增量累积量与高精度实际速度增量累积量相对误差为6.83%，J2摄动速度增量累积量与高精度实际速度增量累积量相对误差为4.21%。该结果表明前两步全局优化保证了多星交会机动成本消耗全局近似最优。
图8算例Ⅰ速度增量累积变化过程
Fig.8Variation process in velocity increment cumulative quantity of example Ⅰ
不同的交会上限时间Δtmax对于脉冲速度增量存在显著影响，按式（9）将其分别设置为目标星飞行X=10圈、X=20圈、X=30圈的时间并进行对比，结果如图9所示。当Δtmax较大时，整体速度增量累积量更小，这是由于连续变量设计空间较大时，允许较长时间的Lambert转移，从而可能搜索到脉冲速度增量较小的单个转移弧段。
图9算例Ⅰ不同交会上限时间速度增量对比
Fig.9Comparison of velocity increments at different rendezvous time upper limits of example Ⅰ
3.3.2 特性分析及对比
为了进一步说明优化方法的适用性，进一步选取7颗目标星且初始参数差异较大的算例Ⅱ进行计算仿真，目标星参数如表3所示，最终优化结果如表4所示，该算例计算结果说明：所提方法在较大轨道面差时仍然能满足多星交会精度要求。
高精度解的机动成本是否最小直接取决于简化模型近似解的全局最优性，为进一步阐述方法特性，将本文提出分步优化方法同混合整数优化方法进行对比，各自针对J2解析摄动模型近似解进行5次优化，计算结果如表5所示。其中分步优化方法5次优化中最优速度增量比混合优化方法大2.17%，但平均速度增量比混合优化方法小1.79%，两种方法优化能力相近，但分步方法平均计算时间仅为混合优化方法的16.93%。结果说明本文分步方法能够保证简化模型近似解全局最优且具有更好的计算效率。
表3算例Ⅱ目标星参数
Tab.3Target satellites′ parameters of example Ⅱ
表4算例Ⅱ多星交会优化结果
Tab.4Optimization result of multi-satellite rendezvous of example Ⅱ
表5不同优化方法J2近似解计算结果对比
Tab.5Comparison of calculation results of J2 approximate solutions using different optimization methods
简化模型近似解与高精度解存在一定差距，对算例I中目标星按各自交会时间基于不同摄动模型进行外推预报，图10、图11分别从升交点赤经和纬度幅角的角度展示了各个目标星在交会时刻时二体模型与J2摄动模型相较于高精度模型的偏差。对于纬度幅角，二体模型与J2摄动模型都存在一定的偏差，最后一次交会时二体模型的纬度幅角偏差比J2摄动模型大42.9°，而图10显示二体模型与高精度模型存在持续增长的显著的升交点赤经偏差，该偏差造成轨道面角度差的不断扩大，在交会轨迹优化过程中当轨道面角度超出了局部修正算法的调整能力上限时，则会加剧计算不收敛的情况。
图10不同目标星交会时刻模型RAAN误差对比
Fig.10Comparison of RAAN errors of different models in rendezvous time of every target satellite
图12给出了两种简化模型近似解直接进行局部修正的计算结果，两种方法迭代结束后的末端相对距离偏差都存在持续扩散的趋势，且二体模型近似解局部修正结果因轨道面偏差值增大会出现突变扩散的情况。该结果说明，以J2摄动模型脉冲方案作为中间近似解，且进行目标星轨道状态的多轮迭代更新，以减小收敛难度是必要的。
图11不同目标星交会时刻模型纬度幅角误差对比
Fig.11Comparison of latitude angle errors of different models in rendezvous time of every target satellite
图12不同模型近似解直接局部修正结果对比
Fig.12Comparison of direct local correction results for approximate solutions based on different models
以末端相对距离为指标，将上述两种方法以及多轮迭代修正方法分别运行5次，计算结果对比如表6所示。
算例仿真中将交会末端可容许的误差设置为5 km。前两种直接修正的方法都存在交会末端位置偏差扩散的情况，且后续弧段末端位置平均偏差并未收敛至公差范围内，而基于多轮迭代修正方法进行计算，各段交会末端位置平均偏差均小于1×10-6 km，达到了高精度多星交会任务所需的精度要求。
表6不同优化方法的收敛情况对比
Tab.6Convergence situation comparison of different optimization methods
4 结论
针对单追踪星在轨服务多目标星的高精度交会轨迹优化问题，本文提出了一种基于多轮迭代修正的三步串行优化方法，仿真结果表明：
1）提出的三步优化串行方法能够在较少迭代次数下得到单个转移弧段的脉冲交会高精度解，追踪星和目标星末端交会状态吻合良好。
2）目标星不同模型轨道的偏差随时间扩散，以多星交会简化模型近似解为基础直接局部修正，难以得到满足精度要求的有效收敛。
3）提出的多轮迭代修正三步优化方法能够获得摄动条件下多星交会任务的近最优高精度轨道机动参数。
最后，本文主要针对远距离交会场景的多星交会优化方法进行研究，尚未考虑实际任务中近距离交会轨道调整策略和约束条件，下一步将针对综合考虑远近距离交会约束的多星交会任务进行整体优化研究。
