回归轨道的轨迹周期性重复特性是由于轨道运动与中心天体旋转发生共振而形成。因而在轨道设计中需要考虑非球形引力摄动，而在以往多数关于地球回归轨道的研究中假设轨道主要的摄动因素来自于地球非球形引力分布。随着地球非中心引力场理论研究的进步和航天工程实践的深入发展，不少****相继尝试采用高阶引力摄动模型进行轨道设计研究。Lara[1]提出了一种基于“修正-预测”的自动求解数值方法，可得到高阶带谐项摄动下的回归轨道初值。随后，在考虑二阶带谐项和田谐项摄动下，从动力学系统的观点，Lara[5]将回归轨道考虑为分叉于地球固连坐标系下的赤道面内的周期轨道，并研究了稳定性。更进一步，通过将J2摄动下的回归轨道作为初值，Lara和Russell[6]采用微分修正算法实现了在完全地球引力摄动模型下的回归轨道快速设计。此外，Aorpimai和Palmer[2]研究了考虑带谐项摄动系数为J4和J22时，利用周转圆轨道根数求解满足回归条件的轨道。杨盛庆等[7]提出了适用于高精度轨道动力学模型的迭代修正方法可获得严格回归轨道。以上研究实现了在高阶甚至完全地球摄动引力下进行轨道设计以保证足够精度要求的目标，但是很少有研究在设计阶段直接考虑非地球引力摄动(如大气阻力、太阳辐射光压和日月引力)的作用，而它们将在轨道控制阶段加以考虑，因为非保守力尤其是大气阻力所引起的星下点轨迹漂移需要额外施加控制力进行抑制。
因此，若缺少必要的轨道维持或轨道设计初值存在一定的误差，长期运行的回归轨道在实际动力学环境中则会失去回归特性。针对该问题，一些****在轨道设计的基础上面向具体的轨道控制目标提出了一些方法。受限于星上设备体制和地面站处理能力，Aorpimai和Palmer[2]提出的多脉冲自主控制策略可将卫星由初始条件配置到回归轨道条件。基于半解析方法，Sengupta等[8]研究了在J2摄动和大气阻力作用下对地覆盖小偏心率回归轨道的控制问题。温生林等[9]同样考虑在J2摄动和大气阻力摄动下，基于Lyapunov理论设计了回归轨道卫星星下点轨迹保持的相对平均轨道根数反馈控制律。张冲难等[10]针对轨道控制时间、燃料消耗、偏心率等约束条件，给出了多脉冲轨道控制策略的具体实现方法进行回归轨道维持。针对回归轨道对地连续覆盖维持问题，Fu等[4]基于纬度幅角分析了整个星下点轨迹漂移量，并提出了维持轨迹漂移不超过给定阈值的控制策略。
本文研究了一种基于微分代数运算和考虑完全摄动因素的引力模型下高精度回归轨道设计与控制的半解析方法。通过对Poincaré映射进行高阶Taylor展开以获得轨道在一个或者多个回归周期内的状态量。该方法一方面通过研究施加于映射交点(即赤道升交点)处的速度增量对回归模式的作用以实现高精度回归轨道设计与轨道控制量优化求解；另一方面通过采用多项式运算代替传统数值积分以避免在完全引力摄动下进行长周期回归轨道递推造成的计算复杂性，从而提高回归轨道设计与轨道控制量生成效率，对星上的自主实施具有重要意义。
1 问题描述与建模 1.1 坐标系定义 图 1给出地心惯性坐标系和地心固连坐标系的定义。取赤道面为地心惯性坐标系的基本平面，
 |
| 图 1 地心惯性坐标系和地心固连坐标系 Fig. 1 Earth-centered inertial and Earth-centered, Earth-fixed coordinate frames |
| 图选项 |
继承文献[1]的变量定义，卫星在惯性空间中的位置由圆柱坐标(r, z, ?)确定，r为卫星到?轴的距离，z为卫星距离赤道面的高度，?为卫星子午面的瞬时经度。卫星在地心固连坐标系中的位置与速度表示为轨道状态量X=[x, y, z, vx, vy, vz]T，则其星下点轨迹的纬度φ和经度λ分别满足sin φ=z/ρ和tan λ=y/x(ρ表示卫星与地心之间的距离)。由于赤道处的星下点轨迹漂移最大，以下研究中仅需要考虑卫星向上穿越赤道面时的状态量。
回归轨道实际上为中心天体固连坐标系中的周期轨道，可通过一些数值方法求解，如文献[5-6]提出的“预测-修正”算法。该方法针对保守的中心天体引力场是有效的，但当加入非保守力的影响时，将几乎无法生成周期轨道。根据求解周期轨道的思路，回归轨道的初始状态X0须与经过特定回归圈次之后的终止状态Xf充分接近。
由于近地轨道动力学环境中存在各类摄动力的影响，特别是非保守力，即便是满足回归条件的轨道也会因摄动作用，相对于参考轨道发生漂移。在缺少必要的轨道维持下，实际轨道将逐渐偏离回归条件，并最终导致任务的失败。需要指出的是，根据不同的任务要求，轨道维持并不必过于频繁和严格，而仅需实际轨道与参考轨道的偏差不超过预设的阈值即可。
1.2 回归轨道条件 当轨道满足共振条件时，即卫星的平均角速度与地球的自转角速度可约，此时轨道为回归轨道，则具有如下关系：
 | (1) |
式中：ΔΩd为在一个交点周期Td内升交点赤经的漂移量；nM为回归周期；nN为在一个回归周期内的轨道圈数。当回归轨道具有严格的nM∶nN回归模式时，星下点轨迹在赤道位置上的经度(仅考虑升交点处)表示为
 | (2) |
式中：λ0为回归轨道起点处的经度；λj为第j圈轨道在升交点处的经度。式(2)可作为标称轨道的基准，以评估实际轨道偏离标称设计的程度。
1.3 严格和宽松回归条件 从实际任务实现的角度来说，回归轨道的任务要求可分为2类：严格和宽松。相应地，回归轨道条件可定义为一个回归周期内的精确回归解和多个回归周期内的有界解。前者表示回归轨道在经过一个回归周期后返回初始位置，后者表示回归轨道可在多个回归周期内返回初始位置。
对于精确回归解的要求，回归轨道在一个回归周期内在地心固连坐标系下的初始状态X0等于终止状态Xf。对于回归轨道有界解，回归轨道在m个回归周期后的终止状态Xf等于初始状态X0。在该情况中，起始于初始有界解，轨道在到达m个回归周期前将会出现偏离，但通过对轨道在第m个回归周期时的状态施加约束条件X0=Xf，轨道将会返回至初始状态X0附近并与之保持一定偏差，故称为有界。需要说明的是，当回归周期数m=1时，回归轨道有界解即约化为精确解。在实际轨道设计问题中，可根据期望的精度和轨道控制频率来确定采用何种解。若任务具有严格的精度要求，可根据精确解进行轨道设计并在每个回归周期内进行轨道维持；而对于宽松精度要求，可选择有界解进行设计，并在多个回归周期内进行一次轨道维持。
2 基于微分代数的高阶Poincaré映射 2.1 微分代数方法 微分代数方法起源于人们尝试利用代数手段求解解析问题，其主要思想是在计算机环境中用类似于用浮点数近似实数的方式来处理函数及其运算。
利用微分代数方法可直接在计算机环境中获得n维函数的任意k阶Taylor展开式，并计算相应函数在某点的值，而所需计算量是固定的。在微分代数框架下进行所有的计算可实现一般的常微分方程关于初始条件的直至任意阶数的相流Taylor展开。
不失一般性，考虑下述初值问题
 | (3) |
及其对应的相流ψ(t; x0)。若采用Runge-Kutta等传统数值积分方法，按照计算机的浮点运算规则对一组初始点x0仅可以获得单条轨道。而若将x0初始化为微分代数表示形式，并在微分代数框架下进行所有数值积分中涉及的运算，可以得到相流关于x0的任意阶数的Taylor展开从t0到tf的积分，其多项式表达形式为ψ(tf; x0+δx0)，记为
 |
| 图 2 x0邻域内初始点(x0+δx0)在tf处的近似解 Fig. 2 Approximation of initial value (x0+δx0) in the neighborhood of x0 at tf |
| 图选项 |
同时，仅需要将数值积分运算操作替换为对应的微分代数运算，即可将标准的显式积分模式转化为对应的微分代数形式。而涉及到的步长控制和误差估计只需要由Taylor展开多项式的常值部分来确定，即作为Taylor展开点的参考轨道。因此，所得到的积分结果是以数值方法得到结果的Taylor展开多项式，即相流的数值近似。
微分代数的主要优势在于：无需对变分方程进行推导和积分处理即可得到流的高阶展开式；仅需要将浮点数的代数运算替换为微分代数运算即可实现，因为微分代数方法同常微分方程是独立的。另外，利用COSY INFINITY软件可实现微分代数运算，在有限的计算时间内即可得到高阶展开式[11-12]。
2.2 高阶Poincaré映射 为了确定回归轨道的状态量，基于微分代数方法求解Poincaré映射，该映射可将截面(赤道面)上的任意点在经过一个回归周期后同样投影至该截面[13]。
假设轨道的回归模式为nM∶nN，本算法以满足回归和太阳同步约束的不动点轨道为参考点[13]。不动点状态量在经过从地心惯性坐标系到地心固连坐标系转换后的状态量为X*=[x*, y*, z*, vx*, vy*, vz*]T，回归周期取为T*=nNTd*，并令z*=0，即考虑回归轨道的起点总是在赤道面上。将状态量x、y、vx、vy、vz和回归周期T初始化为微分代数变量，并在完全引力摄动模型下进行轨道递推(时间从t=0到t=T)。由此可得到Poincaré映射的高阶Taylor展开式为
 | (4) |
通过求解下面的参数化隐式方程可消除T自由度：
 | (5) |
即令满足多项式映射(4)的分部z为0。基于微分代数方法可得
 | (6) |
将其回代至映射(4)，即可求得
 | (7) |
需要说明的是，由于轨道递推是在地心惯性坐标系下进行的，而映射(4)~(6)及参数化隐式方程(5)中涉及的轨道状态量均表示在地心固连坐标系中，故在进行微分代数运算时，需要计算从地心惯性坐标系到地心固连坐标系的转换矩阵M，坐标转换需考虑地球章动和极移效应，故矩阵M为时变的，可通过下面的一阶近似：
 | (8) |
将转换矩阵表示为微分代数形式。式中：T0为Taylor多项式(6)中的常数项，而M′(T0)为M在T0时刻的近似变化率，可通过M在T0附近的线性变化近似得到，时间变化δT取为微分代数变量。
通过Poincaré映射(7)，可将在赤道面上参考点附近的任意初始点X0经过一个回归周期后投影至赤道面，且映射(6)为所需时间(回归周期)。求解Poincaré映射需要进行关于6个变量的微分代数积分，因此相对于普通的浮点数积分需要更多的计算时间。但一旦获得了该映射，便可通过简单的多项式代入运算精确近似轨道递推，极大地减少计算量。
3 回归轨道设计 作为实现回归轨道控制的基准[14]，回归轨道设计是一个优化求解满足目标条件初值的过程，定义该多目标函数如下：
 | (9) |
式中：Δv为偏离初始猜测的速度偏差量，可直接取自微分代数变量(Δvx, Δvy, Δvz)；而x0 (v0)和xf (vf)分别为卫星的初始位置(速度)和在经过一个或几个回归周期后升交点处的终止位置(速度)；
由于在回归轨道设计中不可避免地需要进行轨道递推积分，通常不考虑更复杂的摄动因素或在复杂摄动下仅考虑到较短时间周期的轨道计算，以避免计算复杂度的增加，而本文利用微分代数运算的Poincaré映射多项式来高精度近似常用的轨道积分过程，在计算效率上更具有优势。
采用实际太阳同步回归轨道任务的回归模式作为仿真算例以阐释本节所提出的方法，其中TerraSAR-X的回归模式为11∶167，Landsat-8的回归模式为16∶233，IRS-P6的回归模式为24∶341，SPOT-7的回归模式为26∶379。对于这4种任务，一个回归周期内的实际轨道在每圈升交点处的经度与其标称值的对比如图 3所示，其中经度标称值由式(2)计算得到，而实际轨道则通过对求解优化问题(9)所得到的初值进行轨道递推得到。
 |
| 图 3 升交点处经度实际值与标称值的对比 Fig. 3 Comparison of actual and nominal longitude values at ascending nodes |
| 图选项 |
由图 3所示，对于不同回归模式的回归轨道，经度的实际值处于标称值所表示的“□”内，说明由本节所提出的设计方法获得的轨道初值精度可得到保证。进一步地，将图 3中每一圈升交点处经度的实际值与标称值之间的误差表示出来，如图 4所示(图中nN为无量纲值)。对比发现，实际经度偏离标称值的误差大小不超过0.008°，对应在赤道上的漂移距离为0.89 km。同时可以发现，TerraSAR-X模式的轨道具有最大的经度漂移，这是因为其对应轨道高度(半长轴为6 883.513 km)相对于其他3种模式最低，受到的大气阻力摄动作用最强。需要指出的是，本节求解得到的初值和轨道递推均在完全引力摄动模型下进行且没有消除短周期项，因此经度误差的变化出现了短周期的振荡。
 |
| 图 4 一个回归周期内每圈轨道升交点处经度实际值与其标称值对比的漂移量 Fig. 4 Drift value of actual longitude compared with nominal one at ascending nodes of each cycle of orbit during one repeat cycle |
| 图选项 |
4 脉冲轨道控制策略 4.1 控制策略设计 对于任何采用回归轨道的空间任务，要面临的一个主要问题是：当卫星偏离参考轨道一定范围时，需要施加周期性控制以恢复至回归轨道条件，否则任由偏差增大将导致任务失败。根据第3节所确定的回归轨道初值，卫星在经过一个或若干个回归周期后终止状态将会偏离初始状态。因此，为消除该偏差，本节设计一种脉冲轨道控制策略以进行轨道维持，步骤如下：
步骤1 ??根据在轨道设计阶段所提出的设计方法，可得到满足多目标函数(9)的回归轨道初值。
步骤2 ??对所得到初值进行轨道积分，得到与设计阶段采用的时间长度(一个或者多个回归周期)相同的轨道状态，并利用此时的轨道状态重新构造微分代数映射。
步骤3 ??为得到下一个(或多个)回归周期内的初值，即脉冲控制的目标值，利用在步骤2中重新构造得到的微分代数映射，通过优化方法对控制问题，即式(10)进行求解。
步骤4 ??根据步骤2得到的第一个(或者多个)回归周期结束时的末状态量与步骤3得到的目标值的速度差值，即可获得轨道控制所需要的速度脉冲。
在由上述轨道控制方法得到的速度脉冲作用下，所有的星下点轨迹将会维持在标称轨迹附近期望的偏差阈值内。需要说明的是，控制目标值是根据控制精度要求而决定的，为得到该目标值可将控制问题表示为
 | (10) |
式中：优化变量Δv=[Δvx, Δvy, Δvz]T；xt、vt和
需要说明的是，采用微分代数方法所得到的Taylor多项式映射通常对多个回归周期(一般为2~3个)均有效，可用来近似真实的轨道递推结果。该映射可以通过地面离线计算得到，并在卫星入境可见时上注至星载计算机，并在随后的1~2个回归周期内(直至多项式精度发散)由星上进行多项式运算即可。由于多项式的计算仅涉及乘法和加法运算，且无需每个回归周期内均上注轨道控制指令，在线计算并不需要过多消耗星上有限的CPU计算资源和占用地面上行资源，便于星上进行轨道递推计算，该优点对星上自主轨道控制具有重要的作用。
4.2 严格精度情形 本节将说明具有严格精度要求的控制策略。在本情形下的算例中，取位置漂移、速度漂移和升交点赤经漂移的阈值分别为10-6 km、10-3 km/s和10-7(°)/s以确定具有严格精度要求的轨道目标状态量。为保持一致性，本节同样以在设计阶段所采用的4个实际回归轨道任务模式为例。
实施脉冲轨道控制前后的2个回归周期的状态量如表 1~表 4所示，各表的第2列和第4列中的位置、速度分别为开始第1个和第2个回归周期的初始条件，而第3列和第5列中的位置、速度分别为第1个和第2个回归周期结束时的终止状态。在第1个回归周期结束时，通过施加脉冲控制使卫星到达目标状态，即第2个回归周期的初始状态，并随后开始第2个回归周期。所需的速度增量Δv只需通过对比第3列和第4列的速度分量即可，各回归模式算例(TerraSAR-X、Landsat-8、IRS-P6和SPOT-7)所需大小分别为6.817 8 cm/s、6.607 0 cm/s、9.728 1 cm/s和13.847 6 cm/s。通过以上脉冲机动，即可将轨迹偏差维持在给定阈值内，并满足严格的精度要求。
表 1 TerraSAR-X回归模式算例轨道控制结果 Table 1 Orbital control results of TerraSAR-X repeat pattern
| 状态量 | 第1个回归周期 | 第2个回归周期 | |||
| 初始 | 终止 | 初始 | 终止 | ||
| x/km | 6 784.509 118 | 6 784.508 908 | 6 784.508 908 | 6 784.508 491 | |
| y/km | -1 198.825 274 | -1 198.825 404 | -1 198.825 404 | -1 198.825 332 | |
| z/km | 0 | 0 | 0 | 0 | |
| vx/(km·s-1) | -0.268 267 730 | -0.268 211 934 | -0.268 272 691 | -0.268 244 003 | |
| vy/(km·s-1) | -1.463 279 438 | -1.463 277 386 | -1.463 278 583 | -1.463 459 538 | |
| vz/(km·s-1) | 7.544 055 446 | 7.544 025 069 | 7.544 055 978 | 7.544 001 930 | |
| T/s | 950 399.38 | 950 399.64 | |||
表选项
表 2 Landsat-8回归模式算例轨道控制结果 Table 2 Orbital control results of Landsat-8 repeat pattern
| 状态量 | 第1个回归周期 | 第2个回归周期 | |||
| 初始 | 终止 | 初始 | 终止 | ||
| x/km | 6 975.552 101 | 6 975.551 558 | 6 975.551 558 | 6 975.551 381 | |
| y/km | -1 232.582 639 | -1 232.582 799 | -1 232.582 799 | -1 232.582 778 | |
| z/km | 0 | 0 | 0 | 0 | |
| vx/(km·s-1) | -0.285 046 636 | -0.285 007 355 | -0.285 073 308 | -0.285 030 544 | |
| vy/(km·s-1) | -1.561 515 791 | -1.561 706 363 | -1.561 705 676 | -1.561 855 639 | |
| vz/(km·s-1) | 7.426 314 342 | 7.426 284 548 | 7.426 288 429 | 7.426 264 784 | |
| T/s | 1 382 399 | 1 382 399.6 | |||
表选项
表 3 IRS-P6回归模式算例轨道控制结果 Table 3 Orbital control results of IRS-P6 repeat pattern
| 状态量 | 第1个回归周期 | 第2个回归周期 | |||
| 初始 | 终止 | 初始 | 终止 | ||
| x/km | 7 091.011 029 | 7 091.010 546 | 7 091.010 546 | 7 091.010 247 | |
| y/km | -1 252.984 275 | -1 252.984 999 | -1 252.984 999 | -1 252.984 943 | |
| z/km | 0 | 0 | 0 | 0 | |
| vx/(km·s-1) | -0.295 622 153 | -0.295 550 878 | -0.295 647 647 | -0.295 609 033 | |
| vy/(km·s-1) | -1.623 266 169 | -1.623 412 825 | -1.623 403 153 | -1.623 752 336 | |
| vz/(km·s-1) | 7.356 222 412 | 7.356 200 355 | 7.356 202 783 | 7.356 149 740 | |
| T/s | 2 073 599.4 | 2 073 600.5 | |||
表选项
表 4 SPOT-7回归模式算例轨道控制结果 Table 4 Orbital control results of SPOT-7 repeat pattern
| 状态量 | 第1个回归周期 | 第2个回归周期 | |||
| 初始 | 终止 | 初始 | 终止 | ||
| x/km | 7 096.348 844 | 7 096.348 569 | 7 096.348 569 | 7 096.348 374 | |
| y/km | -1 253.927 469 | -1 253.928 444 | -1 253.928 444 | -1 253.928 404 | |
| z/km | 0 | 0 | 0 | 0 | |
| vx/(km·s-1) | -0.296 089 274 | -0.296 009 0 | -0.296 111 439 | -0.296 062 414 | |
| vy/(km·s-1) | -1.625 981 449 | -1.626 113 0 | -1.626 098 844 | -1.626 408 561 | |
| vz/(km·s-1) | 7.353 025 818 | 7.353 006 575 | 7.353 009 099 | 7.352 962 120 | |
| T/s | 2 246 398.4 | 2 246 399.8 | |||
表选项
不同于采用调节半长轴平根的方式进行轨道维持，本文方法直接通过速度修正进行轨道控制，改变瞬时轨道根数，如图 5~图 7所示。图中各轨道圈次升交点处半长轴a、轨道倾角i和升交点赤经漂移率
 |
| 图 5 施加脉冲机动前后半长轴的演化 Fig. 5 Evolution of semi-major axis before and after impulsive maneuvers |
| 图选项 |
 |
| 图 6 施加脉冲机动前后轨道倾角的演化 Fig. 6 Evolution of inclination before and after impulsive maneuvers |
| 图选项 |
 |
| 图 7 施加脉冲机动前后升交点赤经漂移率的演化 Fig. 7 Evolution of drift rate of right ascension of ascending node before and after impulsive maneuvers |
| 图选项 |
因此，基于严格回归轨道条件的轨道设计及每个回归周期施加一次脉冲机动的轨道控制方法可以用来完成具有严格精度要求的回归轨道任务。
4.3 宽松精度情形 考虑到在实际任务设计中，存在着对精度要求比较宽松的情况，如允许设计轨道同标称轨道之间存在一定的偏差，此时不必如严格精度情形那样在每个回归周期内均进行一次脉冲机动控制，从而对燃料的消耗也将随之大大减少。
本节以萨瑞卫星技术有限公司实施的UoSAT-12地球观测任务[15]为例，对宽松精度情形下的轨道控制策略进行说明，并将所得结果同已知结果[2]进行对比。
UoSAT-12卫星质量为300 kg，轨道倾角为64.75°，高度约为650 km，轨道回归模式为7天102圈。由于该轨道不是太阳同步轨道，在求解作为初始猜测的不动点时，需将太阳同步约束替换为倾角约束。基于微分代数方法求解得到的Poincaré映射(7)在多个回归周期内也具有可靠的精度，故只需要在完全引力摄动模型下进行一次轨道递推得到映射(7)的Taylor多项式，并随后通过代入多项式计算即可近似轨道递推结果(在本例中需代入求解2次)。由此，可以确定3个回归周期内的全部轨道状态。采用宽松回归条件进行轨道设计，容许实际轨迹相对标称轨迹存在一定的漂移，但在第3个回归周期时轨道状态将返回至初始轨道状态附近。
根据本例的控制要求，仅需初始位置与3个回归周期后的终止位置偏差保持一定的距离范围即可，则控制优化问题(10)转化为
 | (11) |
且在本例中将距离阈值xt设定为2 km，以此说明回归轨道宽松情形的控制策略。
由于采用代入Taylor多项式(7)的方法来近似求解本应由轨道递推得到的实际状态量，先需要检验该方法的精度。对比采用微分代数方法得到的每一个回归周期处的星下点轨迹经度与在相同时机采用轨道递推方法计算实际的经度，结果如图 8所示。图中：回归周期个数m无量纲。可以发现，在3个回归周期内采用微分代数方法得到结果能较好地吻合利用轨道递推方法所得到的结果。
 |
| 图 8 基于微分代数方法和轨道递推方法的经度对比结果 Fig. 8 Comparison of longitude values obtained by DA-based orbital propagation and numerical integration-based orbital propagation |
| 图选项 |
如前所述，采用微分代数方法利用宽松回归条件，仅涉及1个回归周期内的轨道递推和另外2次Taylor多项式的代入求值，相比于直接采用3个回归周期内的轨道递推方法，大大减少了计算量，实现了轨道的快速和准确的设计。在3个回归周期内升交点处的实际和标称经度值的对比结果如图 9(a)所示，相应地，二者的偏差表示如图 9(b)所示，图中nN无量纲。
 |
| 图 9 三个回归周期内经度的实际值与标称值对比 Fig. 9 Comparison of actual and nominal longitude values during three repeat cycles |
| 图选项 |
根据本节的轨道控制策略和宽松精度要求，可得到轨道控制前后6个回归周期内的轨道状态量如表 5所示，回归轨道的初值如第2列所示，在前3个回归周期内起始于该状态的轨道将会首先发生漂移，如表 5中的第3列状态量所示，该过程的经度偏差已表示在图 9(b)中；在第3个回归周期时施加脉冲机动Δv，其结果由第4列与第3列的速度差值计算得到。图 10给出了在随后的3个回归周期内，实际轨迹与初始状态的经度偏差，可见不超过0.001 63°，约为1.81 km。相应地，轨道控制施加前后轨道在真赤道坐标系下的半长轴和轨道倾角的变化关系如图 11所示。控制频率为每3个回归周期(21天)进行1次轨道控制，所施加的脉冲控制量大小Δv为3.858 3 cm/s，平均每天仅需要1.837 3 mm/s的控制消耗。该结果同文献[2]给出的结果(每天1次轨道控制，平均每次轨道控制量为3.6 mm/s，轨迹漂移平均值为2.72 km)一致，但是本文提出的方法基于完全引力摄动模型，且所需轨道控制频率低、轨道控制计算效率高同时具有更好的轨迹抑制效果(偏差不超过1.81 km)。
表 5 宽松精度情形下轨道控制结果 Table 5 Orbital control results for loose-accuracy scenario
| 状态量 | 前3个回归周期 | 后3个回归周期 | |||
| 初始 | 终止 | 初始 | 终止 | ||
| x/km | -1 684.591 611 | -1 684.515 457 | -1 684.515 457 | -1 682.767 179 | |
| y/km | -6 825.069 865 | -6 825.112 359 | -6 825.112 359 | -6 825.595 939 | |
| z/km | 0 | 0 | 0 | 0 | |
| vx/(km·s-1) | 2.638 596 174 | 2.638 601 972 | 2.638 618 204 | 2.638 975 345 | |
| vy/(km·s-1) | -0.666 368 889 | -0.666 251 883 | -0.666 255 966 | -0.665 570 472 | |
| vz/(km·s-1) | 6.802 699 118 | 6.802 626 184 | 6.802 660 947 | 6.802 471 580 | |
| t/s | 1 794 301.5 | 1 794 297.6 | |||
表选项
 |
| 图 10 轨道控制后轨迹经度与其初始值的偏差 Fig. 10 Deviation of actual orbital longitude from initial longitude after orbital control |
| 图选项 |
 |
| 图 11 六个回归周期内半长轴和轨道倾角的演化 Fig. 11 Evolution of semi-major axis and inclination within six repeat cycles |
| 图选项 |
5 结论 1) 面向严格和宽松任务约束的回归轨道任务，提出了基于微分代数的高阶Poincaré映射方法，解决了在完全高精度引力模型下回归轨道设计与控制问题。
2) 针对严格任务约束，在每个回归周期内均要求初始与终止轨道状态量相等，且在每个回归周期内施加一次脉冲机动。针对宽松任务要求，通过Poincaré映射代替轨道积分，实现了在多个回归周期内的一次轨道设计和在此时间内的一次脉冲控制。
3) 本文方法可极大地减少计算量，并具有对计算资源依赖较小的优势，便于星上自主进行轨道递推计算，可用于星上自主轨道维持。
参考文献
| [1] | LARA M. Searching for repeating ground track orbits: A systematic approach[J]. Journal of the Astronautical Sciences, 1999, 47(3): 177-188. DOI:10.1007/BF03546198?utm_medium=affiliate&utm_content=meta&utm_campaign=DDCN_1_GL01_metadata |
| [2] | AORPIMAI M, PALMER P L. Repeat-ground track orbit acquisition and maintenance for earth-observation satellites[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(3): 654-659. DOI:10.2514/1.23413 |
| [3] | NADOUSHAN M J, ASSADIAN N. Repeat ground track orbit design with desired revisit time and optimal tilt[J]. Aerospace Science and Technology, 2015, 40: 200-208. DOI:10.1016/j.ast.2014.11.007 |
| [4] | FU X, WU M, TANG Y. Design and maintenance of low-earth repeat-ground-track successive-coverage orbits[J]. Journal of Guidance, Control, and Dynamics, 2012, 35(2): 686-691. DOI:10.2514/1.54780 |
| [5] | LARA M. Repeat ground track orbits of the Earth tesseral problem as bifurcations of the equatorial family of periodic orbits[J]. Celestial Mechanics and Dynamical Astronomy, 2003, 86(2): 143-162. DOI:10.1023/A:1024195900757 |
| [6] | LARA M, RUSSELL R P. Fast design of repeat ground track orbits in high-fidelity geopotentials[J]. Journal of the Astronautical Sciences, 2008, 56(3): 311-324. DOI:10.1007/BF03256555 |
| [7] | 杨盛庆, 杜耀珂, 陈筠力. 基于迭代修正方法的严格回归轨道设计[J]. 宇航学报, 2016, 37(4): 420-426. YANG S Q, DU Y K, CHEN J L. Design of strictly-regressive orbit based on iterative adjustment method[J]. Journal of Astronautics, 2016, 37(4): 420-426. DOI:10.3873/j.issn.1000-1328.2016.04.007 (in Chinese) |
| [8] | SENGUPTA P, VADALI S R, ALFRIEND K T. Satellite orbit design and maintenance for terrestrial coverage[J]. Journal of Spacecraft and Rockets, 2010, 47(1): 177-187. DOI:10.2514/1.44120 |
| [9] | 温生林, 闫野, 张华. 低轨回归轨道卫星轨迹漂移特性分析与控制[J]. 系统工程与电子技术, 2015, 37(3): 613-619. WEN S L, YAN Y, ZHANG H. Analysis and control of groundtrack drift for recursive low Earth orbit satellites[J]. Systems Engineering and Electronic, 2015, 37(3): 613-619. (in Chinese) |
| [10] | 张冲难, 卞燕山, 王西京, 等. 严格回归轨道自动生成算法及实现[J]. 推进技术, 2018, 39(7): 1472-1478. ZHANG C N, BIAN Y S, WANG X J, et al. Algorithm and implementation of automatic generation of rigorous recursive orbit[J]. Journal of Propulsion Technology, 2018, 39(7): 1472-1478. (in Chinese) |
| [11] | BERZ M. Modern map methods in particle beam physics[M]. London: Academic Press, 1999. |
| [12] | BERZ M, MAKINO K. COSY INFINITY version 8.1 user's guide and reference manual[Z]. East Lansing: Michigan State University, 2002. |
| [13] | HE Y, ARMELLIN R, XU M. Bounded relative orbits in the zonal problem via high-order Poincaré maps[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(1): 91-108. DOI:10.2514/1.G003612 |
| [14] | 杨盛庆, 杜耀珂, 王文妍, 等. 严格回归轨道的管道导航方法研究[J]. 中国空间科学技术, 2017, 37(6): 10-16. YANG S Q, DU Y K, WANG W Y, et al. Study on the tube navigation of strictly regressive orbit[J]. Chinese Space Science and Technology, 2017, 37(6): 10-16. (in Chinese) |
| [15] | FOUQUET M, SWEETING M. UoSAT-12 minisatellite for high performance Earth observation at low cost[J]. Acta Astronautica, 1997, 41(3): 173-182. DOI:10.1016/S0094-5765(97)00181-1 |
