摘 要 ：为了适应低升阻比飞行器再入返回的大航程要求，本文针对大气跳跃再入飞行环境复杂难以直接获得解析解的特点，基于匹配渐进展开法设计了一种跳跃式再入解析预测-校正制导方法。首先分析了低升阻比飞行器大气跳跃再入轨迹的飞行剖面和制导分段方法，然后分别推导了其运动方程大气层外外解和大气层内内解的渐进展开形式，并通过匹配获得了统一的封闭解析表达式，接着基于此解析解实时预测飞行器的剩余航程，并通过不断迭代升阻比垂向分量以满足最后的落点精度，最后针对跳跃再入飞行的不同阶段设计了不同的制导策略以获得最终的倾侧角指令，相应的仿真结果证明了此方法的有效性。

关键词 ：低升阻比；大气跳跃再入；匹配渐进展开法；解析预测-校正；制导

由于目前我国还不能全球布站，测控能力有限，且着陆场的约束较强，其着陆地往往位于内陆地区（阿波罗为海上着陆），因此要求探月飞船必须具有较大的航程能力。对于低升阻比返回舱（神舟和嫦娥都属于此类），要想执行长航程再入返回飞行任务，最好的方法就是采用跳跃式再入，即探测器在进入大气层以后依靠气动力再跃出大气层飞行一段时间后，第二次再入大气层并最终着陆[1] 。

2014年10月，探月工程三期再入返回飞行试验返回器就是采用了半弹道跳跃式再入返回技术，成功返回着陆在内蒙古四子王旗回收区。由于再入返回过程速度大、航程长、高动态、姿态与弹道耦合强烈，这就要求飞行制导系统必须能够提供快速精确的跳跃再入制导指令。大气跳跃再入轨迹最初是为阿波罗登月任务设计的[2] ，现在再次考虑将其用于类似的低升阻比飞行器[3] 。一个跳跃再入轨迹起始于再入点，其具有较大的再入速度（近似第二宇宙速度）和较小的飞行路径角。进入大气层后，气动升力连续调整飞行速度矢量直至飞行器跳出稠密大气。跳出大气以后，飞行器进行开普勒飞行直至它再次进入大气。

这里主要是通过规划大气跳出点的速度和飞行路径角来控制飞行器的总航程。再入过程中倾侧角是唯一控制变量，其通过绕速度矢量旋转飞行器以改变气动升力的垂直分量。阿波罗再入制导系统采用闭环轨迹解预测跳出点速度和路径角，使用的方程为非线性运动方程的简化近似。总航程通过拼接初始大气跳跃段、开普勒段和再次进入大气段的估计航程获得，当预测航程等于剩余航程时再调整倾侧角。尽管阿波罗再入制导方法证明很成功，但其航程预测能力有限（最大4630km，对应41°航程角）[4] 。针对新型低升阻比飞行器的设计目标（包括扩展、中止和其他意外情况下的航程能力），Bairstow[5] and Brunner and Lu [6] 开发了数值积分预测制导方法，此方法通过数值积分整个非线性运动方程组向前推进跳跃轨迹，倾侧角通过迭代方式进行调整直至预测航程等于剩余航程。参考文献[7]中的预测校正制导律分别成功获得了航程10006.2千米（90°航程角）和79679千米（72°航程角），但其计算量相对较大，对机载计算机提出了较高的要求。本文基于匹配渐进展开法[8,,9,10,11] ，提出了一种基于解析再入轨迹的制导方法。首先推导了跳跃再入段大气层外外解和大气层内内解的渐进展开形式，通过匹配获得统一的封闭解析表达式，基于此解析解对跳跃段航程进行预测。然后对各飞行段的制导策略进行设计，通过不断迭代调整升阻比垂直分量，直至得到期望的航程，最终获得大气跳跃轨迹。

1 跳跃再入运动模型

1.1 运动方程质点飞行器在地心惯性系下的平面运动方程为      2 22 2 2dsindd/ 2 / sindd/ 2 / / cosdDLrVtVSC V m rtV SC V m r V rt        (1)式中， r 为到地球中心的半径； t 为飞行时间；V 飞行速度；  为飞行路径角（沿当地水平面向上为正）；  大气密度；LC 和DC 为升力和阻力系数； S 为特征面积； m 为飞行器质量；  为中心天体的引力系数。而且，大气密度假设为指数形式，即，      exp /s sr r r r H        (2)其中，密度模型的匀质大气高度值 7.9km为参考半径，取高度为61.8 km（跳跃机动的拉起高度附近）的地心矢径值。参考半径处的大气密度  sr  取 2.467e-4kg/m^3 。当大气层外飞行时大气密度为零，则方程(1)表示开普勒运动。为了获得近似解析解，方程(1)被转变成以高度为自变量的两个方程。并且引入以下的无量纲化过程，sssref///LDDsVr r Chrr CHCH rbm S  则运动方程转变为  2 2222d 2ed sin1dcos 1 1cos ed 1 21hhv bvhhbh hh v             (3)

1.2 跳跃再入轨迹类似阿波罗再入制导，大气跳跃段轨迹规划是制导律的基础。上图是纵程约为8338.5千米的跳跃再入轨迹各个制导阶段的示意图，下降段从再入点（Entry Interface，EI）处到跳跃再入最小高度处（飞行路径角  为0°），上升段从最小高度处到跃出大气处，在第二次再入大气前为开普勒段，末段从二次大气再入到开伞着陆。注意到开普勒段的高度边界并非必需在EI上，而是在某个阻力加速度充分小的高度上（阿波罗制导中定义“跃出”在阻力加速度小于0.2g时的点上）。因为匹配渐进展开法采用高度为自变量，它是开普勒段开始的“触发器”，开普勒段的开始高度对于纵程小于4632km的情况可以定在65km，对于纵程大于4632km的跳跃轨迹可以定在85km。通过采用较大的跃出高度值能提高大航程跳跃轨迹预测精度，因为大气阻力的影响变小了。

2 基于匹配渐进展开的解析跳跃轨迹

跳跃轨迹起始于重力为主的区域（外部区域），然后进入空气动力为主的区域（内部区域），最后再到达重力为主的区域。参考方程(3)，右边第一项重力起主要作用，第二项空气动力起主要作用。获得方程(3)解析解的困难之处在于没有两个区域都适合的易于获得的解析解的近似形式。解决此困难的一种方式为在每个区域独立获得的解（外部区域对应的解称为外解，内部区域对应的解称为内解），然后再匹配融合成统一的形式（统一解）。匹配渐进展开法一种奇异摄动技术[12] ，可实现此目标。

2.1 外解与内解的渐进展开高度 h 和高度比  的值决定了跳跃轨迹上的某个特定点是在内部还是外部区域。当高度 h 相对于高度比  很大的时候，方程(3)中的空气动力项是小的，并且运动近似为开普勒轨道。因此在外部区域，空气动力项的影响被看做小阶的扰动。而当高度 h 和高度比  相当时，空气动力不能再被看作是小量，并且可认为空气动力远大于重力，令轨迹不能近似为开普勒轨道。因此在内部区域内，把重力看作一个低阶的摄动项。匹配渐进展开法的一个关键特征为不同区域的近似方程来自于对同一方程组的限制处理，并且此限制过程都需要小量趋于零。因此在内区域，方程(3)表示成扩展变量的形式，hh假设外解为如下形式，     2 2 2 2 20 1 2v v h v h v h       (4)     20 1 2cos cos cos cos h h h           (5能够证明20v 和0cos  满足方程，2 201 12 1 2 1iiv vh h   (6)   0 01 cos 1 cosi i ih v h v      (7)其中下脚标 i 标识初始条件。定义复合常值为，2 221iiihv vh (8)  1 coscosi i ih vv  (9)方程(6)和(7)被写成，2 2021hv vh (10) 00coscos1vv h (11)高阶项为，2 21 20 v v     (12)1 2cos cos 0       (13)表达式(4)和(5)的首项即为精确的外部解。实际上，方程(6)和(7)是开普勒运动积分，分别表示能量守恒和角动量守恒。假设内解为如下形式，     2 2 2 2 20 1 2v v h v h v h             (14)     20 1 2cos cos cos cos h h h                (15)其首项为20v 和0cos   ，满足关系式，2 002exp v C     (16)0cos e cos2hb     (17)其中   cosL LC C    为升阻比的竖直分量， C 和    为积分常数，并且与初始条件 , ,j j jh v    有关，22expjjC v    (18)cos cos e2jhjb     (19)高阶项似乎不能以闭环形式获得。参考文献[8]以积分形式给出了一阶解，这里只考虑首项。

2.2 外解与内解的匹配当获得了外解和内解的首项以后，开始构造统一有效解的首项。合成解通过外界和内解相加然后减去它们的公共部分获得，这样得到的合成解将包含四个未知常数 ,cos , ,cos v C      。由于根据已知的初始条件     , ,i i ih v h h    只能解算两个，故还需要增加另外两个边界条件，可以通过在重叠区域对外解和内解进行匹配获得。当前，匹配意味着外解对于小的高度 h 和内解对大的扩展高度 h 之间的协调。匹配外解和内解的首项得到下面的关系式。22exp v C     (20)cos cos     (21)方程左边表示当 h 趋于0（ 0 h  ）时对外解的限制，方程右边表示当扩展高度 h 趋于无穷（ h ）时对内解的限制。得到一致有效合成解为 2 22 2exp1hh          (22)   1/222cos cos 1 2 1exp2h h hb h              (23)其中，cos cos exp2b h         (24)此解近似于方程 (3) 的统一解。

2.3 跳跃段解析解计算流程为了进一步的研究，考虑在给定初始条件    , ,i i ih v h h     （即再入大气前轨迹上的一点，并在外区域）的情况下，如何使用方程 (22) 和(23)获得 近似跳跃轨迹。因为根据假设，在初始点空气力可以忽略，故常数 v  和 cos   的值可以通过方程 (8) 和 (9) 计算，常数 cos    可以通过方程(21)计算。通过这些常数值，由方程 (22) 和 (23) 可计算 v 和  的值。但是，如果实际轨迹为跳跃的，那么初始为负的飞行路径角  随着高度 h 的减小而增加，并且在高度降至最低minh 时为零。但是，存在的问题是内解   （实际上是 cos   ）可能比合成解  先变成零，因为它们的  阶项不同。这种情况下，在 0   附近可能出现条件 cos 1    的情况。为了避免这种情况，不应该用方程 (21) 而应该用方程(25) 计算    ， minmincos 1 cos exp2h bh         (25)其中minh 为飞行路径角  为零时的高度 h 值。一旦 v  和 cos   计算得到，通过方程 (23) 进行迭代可以确定最小高度minh 。这种方式获得的    值将导致方程 (21) 左右两边的  阶项不一致。这是允许的，因为匹配只要求阶一致。在无量纲高度最小minh h  （此高度将跳跃轨迹的下降段 0   和跳出段 0   分开），需要对常数   ,cos ,cos v     进行调整，这种情况产生是由于 v 和  为关于 h 的双值函数。由于飞行路径角在跳出阶段为正，导致    (当 h 时内解的极限）也应为正。更确切地讲，内解关于最小高度minh 的对称性表明，   exitr entry      (26)再入和跳出段的常数关系必须令再入段和跳出段的合成解在高度最小处匹配。相应关系为    2 24expentryexit entryv v      (27)    min122min minmin2cos 1 exp22 11exitexith bh hhv                (28)现在整条跳跃轨道的近似已经完成了。方程(22) 和 (23) 的合成解是近似的基础，方程 (8) 、 (9)以及 (25) ~ (28) 用于计算常数。注意通过方程 (8)和 (9) 计算 v  和   需要给出的初始条件在外区域。如果给出的初始条件沿着再入段轨迹离再入点比较远（空气动力不能忽略），那么需要同时解算方程(22)、(23)和(25)以计算再入常数 v  、  和    。对于跳出段情况是一样的。考虑再入段常数 v  和   。由方程 (8) 和 (9) 可知，他们与大气层外的某点状态     ,i iv h h     有关，并且再入段飞行路径角   应为负。将方程 (8)和 (9) 重写成开普勒转移轨道的能量和角动量积分的形式。vh  (29)  1 cos 1 cosi i iv h v     (30)因此 v  和   的物理意义分别为再入段在0 h （即参考半径sr ）处的速度和飞行路径角。当然，只有忽略飞行过程中空气动力，由常数表示的状态才能够获得。由于可能存在 cos 1    的情况，这本身在数学上并不矛盾，因为 cos   是包括初始条件（方程(9)）的合成解的任意标注。但是，如果使用方程(21)确定 cos    ，那么矛盾就会产生，因为需要计算 cos    的反余弦。正如上面所建议的，通过方程(25)计算    ，不仅能避免一个难题，也能避免此数学矛盾。另一个补救措施为增加sr 的值。但是，如果sr 过大，位于sr 半径圆附近的  阶的边界层将不再是气动力为主。为了使 cos 1    ，sr 将取大气层内相对较高的一个值。那么当 h 为  阶时空气动力为主的假设将不再成立。轨道衰减最好使用规则摄动，而非奇异摄动进行近似。对于一阶近似，规则摄动理论下的轨迹仅由重力决定，而空气动力被当作小的摄动项，即只存在一个外区域。总之，对于再入之前的某点，当初始条件不严格要求满足cos 1    时，此条件的不满足表明这点在匹配渐进展开获得的近似解的适用区域附近。类 似 的 ， 对 于 跳 跃 段 的 常 数   exit v  和  cosexit  为不考虑空气动力时在sr r  或 0 h 时的速度和飞行路径角。一旦通过式 (8) 和 (9) 以及初始状态确定了常值   和   ，便可以应用条件 cos 1   ，通过式(23)确定跳跃轨迹最小高度。因式(23)是关于 h的超越方程，故采用牛顿迭代法确定minh 。庆幸的是，式(23)关于 h 的解析导数可获得，牛顿迭代在迭代4～5次后即可收敛。

3 跳跃再入轨迹制导策略

3.1 纵向航程预测本文提出的跳跃制导通过对升阻比垂向分量 进行迭代，使预测总纵程匹配到达目标落点的需要航程。在再入的初始点，  取一个标称试用值（当纵程不匹配时），通过式(22)和(23)解析的计算出一条大气跳跃轨迹（实际上迭代得到的就是一条参考轨迹）。预测总纵程表示为：total down up Kepler finalR R R R R     (31)跳跃再入上升段和下降段的航程角通过积分航程角对速度的微分来确定：d d dVdt dt dV  (32)cos d Vdt r  (33)2sinDdVadt r    (34)/ d dV  的离散值通过求在高度EIh 和minh（下降段）以及minh 和exith （上升段）之间的解析跳跃解（式(22)和(23)）来计算，采用梯度积分法积分 / d dV  。然后各个纵程可用地球半径乘上行航程角来计算（例如，down down ER r   ，up up ER r   ）。应当注意到上升段和下降段的航程角能通过纯解析式计算，这些解析计算有赖于对速度、飞行路径角以及minh 和exith 处的阻力加速度的预测。开普勒段的航程角通过二体假设运动计算： 21 exitKepler2exit1 cos2cos1 2 cosPP P     (35)其中， 2exitexit circ/ P V V  是预测跃出速度exitV 与当地圆轨速度exitcircV 的平方比，exitcirc exitV r   。跃出速度exitV 和飞行路径角exit 由匹配渐进展开解式(22)和(23)在跃出高度exith 计算得到。末段纵程finalR 通过对存储的参考轨迹线性内插来计算。这里采用了阿波罗制导的标准末段剖面，其 通过常 值升阻 比垂 向分量final0.18  （ 53  ）获得。标准参考轨迹采用速度作为自变量，因此  final exitR f V  ，因为假定在二次再入时的速度等于跃出大气时的速度。实际的需要航程goR 通过飞行器的当前位置向量和目标落点区的位置向量来计算。通过估计剩余飞行时间考虑地球自转效应对目标位置向量的影响。升阻比垂向分量  采用正割迭代法调整，直到预测纵程偏差total goR R  小于46.3km（约25海里）。迭代程序运行两次：首次是在再入点EI，第二次是在最小高度minh 处。在minh 处重计算跳跃轨迹的剩余上升段极大地提高了对大气跃出状态exitV 和exit 的预测，因为匹配渐进展开解的复合常值   、   和    能通过采用飞行器最低点当前状态重新计算。数值仿真结果表明，在所有尝试情况下，对  的迭代在3～5次迭代后收敛（注意到阿波罗制导也采用了总迭代次数限制在7次之内的正割迭代法）。

3.2 倾侧角指令生成由解析航程预测迭代得到的升阻比垂向分量 的开环值通过采用阿波罗反馈策略的闭环项进行扩展[4] 。在下降和上升段，指令升阻比垂向分量为：   ref ref ref C V RDK V V K r r          (36)其中，常值ref 为航程预测迭代得到的结果，参考速度refV 以及参考高度变化率refr 通过式(22)和(23)的当前高度函数形式的封闭解来计算。反馈增益VK 和RDK 基于阿波罗制导增益，速度增益VK （单位s/m）与当前阻力加速度Da 成比例：minexitexit0.0112D DVhD Da aKa a(37)其中，exit 21.22m/sDa  为跃出高度附近的近似阻力加速度，minhDa 为在上拉高度处的估计阻力加速度。高度变化率增益RDK （单位s/m）为：min2exitexit0.0381D DRDhD Da aKa a     (38)在开普勒段，采用常值倾侧角飞行，这里取0.2503   （ 33.44   ）。在末段（二次再入），指令升阻比垂向分量基于线性摄动理论和存储的参考轨迹：  final go pred/CL DR RR   (39)其中预测航程为：   refpred final ref D DDR RR R r r a ar a       (40)该段所有的参考值都通过线性内插已存储的以速度为自变量的参考剖面来计算。偏导数通过伴随法和一个线性化系统来计算，并且它们也存储为速度的函数。阿波罗最终段制导参考轨迹的参数如下表所示。最后得到倾侧角为1maxcos/CL D   (41)其中max/ 0.3 L D  。这里采用了阿波罗侧向制导逻辑，如果横程偏差超过了某一边界值，则进行指令倾侧翻转，其中，边界值的界限与2V 成比例，以限制翻转次数。

4 仿真分析

4.1 仿真初值以 阿 波 罗 指 令 舱 （ Apollo CommandModule， CM）为对象进行制导算法的验证，采用文献【NASA TN D-6725， MAR. 1972】中的数据，参考面积211.9m S  ，质量 5511kg m 。指令 舱 为 质 心 偏 置 对 称 体 ， 配 平 攻 角 范 围147 ~158 （ 22 ~ 23  ）。升力系数、阻力系数通过以马赫数M为自变量对配平气动系数线性插值得到。期望的落点经度为75°，相应的航程为8348.95919km，对落点纬度不进行要求。

4.2 计算结果表 4 再入制导落点参数Table 4 Landing Parameters of Entry GuidanceLanding parameters  Value  UnitLanding Height  617.29  mLanding Longitude  1.309  °Landing Latitude  23.4832  °相应的航程为8348.62090km，落点误差为0.33829km。整个再入过程的高度和倾侧角变化如下图所示。由仿真结果可知，本文所给出的载人飞船返回舱再入制导算法可以很好地实现跳跃式再入，航程达到了8348.62千米，且具有较高精度。图3中倾侧角变化曲线不连续光滑，主要是因为倾侧角指令是分段进行设计的，而段间的衔接问题没有考虑，这一问题将在后续的研究中进行解决。

5 结论

本文基于匹配渐进展开法，设计了一种满足低升阻比飞行器大航程要求的跳跃再入解析预测-校正制导方法。在不同的跳跃再入飞行阶段，设计了不同的制导策略，并给出了相应的仿真结果分析，仿真结果证明了此制导算法的适用性。

参 考 文 献[1]  Du X, Li H Y . Design of Skip Entry Guidance for Lu-nar-Return Vehicle[J]. Collection of the Eighth Aca-demic Conference of Chinese Society of AstronauticsSpace Exploration Technology Committee, 2011: 33-36.(in Chinese)杜昕, 李海阳. 探月飞船跳跃式再入制导律设计[C].中国宇航学会深空探测技术专业委员会第八届学术年会论文集, 2011: 33-36.[2]  Graves C A, Harpold J C. Apollo Experience Report—Mission Planning for Apollo Entry[R]. NASA TN D-6725, March, 1972.[3]  Lu P. Predictor–Corrector Entry Guidance for LowLifting Vehicles[J]. Journal of Guidance, Control,and Dynamics, 2008, 31(4): 1067–1075.[4]  Morth R. Reentry Guidance for Apollo[R]. Massachu-setts Inst. of Technology, Rept. R-532, V ol. 1, Cam-bridge, MA, Jan. 1966.[5]  Bairstow S H. Reentry Guidance with Extended RangeCapability for Low L/D Spacecraft[D]. M.S. Thesis,Massachusetts Inst. of Technology, Cambridge, MA,Feb. 2006.[6]  Brunner C W, Lu P. Skip Entry Trajectory Planning andGuidance[J]. AIAA Paper 2007-6777, Aug. 2007.[7]  Lu P. Predictor-Corrector Entry Guidance for Law-Lifting Vehicles[J]. Journal of Guidance, Controland Dynamics, 2008, 31(4): 1067-1075.[8]  Shi Y Y , Pottsepp L. Asymptotic Expansion of a Hyper-velocity Atmospheric Entry Problem[J]. AIAA Journal,1969, 7(2): 353–355. [9]  Shi Y Y , Pottsepp L, Eckstein M C. A Matched Asymp-totic Solution for Skipping Entry into Planetary Atmos-phere[J]. AIAA Journal, 1971, 9(4): 736–738.[10] Vinh N X, Busemann A, Culp R D. Hypersonic andPlanetary Entry Flight Mechanics[J]. Univ. of MichiganPress, Ann Arbor, MI, 1980: 254–272.[11] Calise A J, Melamed N. Optimal Guidance of Aeroas-sisted Transfer Vehicles Based on Matched AsymptoticExpansions[J]. Journal of Guidance, Control, and Dy-namics, 1995, 18(4): 709–717.[12] Liu R C. Some Methods for Singular Perturbation Prob-lems [D]. Jilin: Institute of Mathematics, Institute of Ji-lin University, 2008. (in Chinese)刘日成. 奇异摄动问题中的若干方法[D]. 吉林: 吉林大学数学研究所, 2008.