利用轨道瞬时根数预报的离轨制动控制方法

龚宇莲，何英姿，李毛毛，李克行

(1. 北京控制工程研究所，北京 100190；2. 空间智能控制技术重点实验室，北京 100190)

0 引 言

航天器离轨制动，属于变轨问题的研究范畴，它特指在轨航天器通过施加推力进行制动，进而返回或进入行星大气层的轨道改变过程[1]。离轨制动需要选择合适的轨道位置，施加恰当的推力矢量，使航天器最终满足再入/进入点约束条件。

早期的离轨制动研究主要以冲量假设和二体轨道假设为前提[2]。文献[2]研究了针对再入角需求选择不同的制动速度及制动角策略。文献[3]提出了一种求解离轨制动脉冲解析算法，并给出了确定最优制动点位置的一种迭代算法。近年来，离轨研究开始考虑固定推力、变推力情况下的变轨问题，文献[4]以制动过程攻角和侧滑角的变化率为控制量，设计了燃料最优的离轨制动方法，文献[5]分别分析了固定推力和可调推力情况下的离轨优化问题，文献[4-5]都以伪谱法为基础求解优化问题。文献[6]设计了一种组合制导方法，采用了燃料最省方案计算初始制动参数进行开环制动，而后采用相对运动控制的方式，跟踪标称轨迹实现高精度的末制导。文献[7]以大升阻比飞行器为研究对象，将再入接口条件取为地心距、速度和速度倾角，采用有限推力离轨制动方式，研究了单次“推-滑”和多次“推-滑”的离轨制动问题。文献[8-9]根据再入轨迹分析和落点经纬度，估计再入轨道最优升交点地理经度。

需要确定的制动参数，通常包含制动点火位置、制动速度大小及方向，这些参数的确定一般是在地面规划完成。为提高飞行器再入点精度，文献[10]通过在地面大量仿真训练BP神经网络建立起控制变量和返回特征参数的映射关系，改进了在轨快速预报算法。随着飞控计算机性能的发展，以及对飞行器快速机动及实时性要求的提高，在轨自主制动方法的研究也引起了学者的广泛兴趣[11-14]。此外制动过程中的制导策略对制动精度的影响同样十分显著，阿波罗飞船采用的PEG方法，在行星探测器着陆问题中得到广泛研究[15-16]。针对再入飞行器，文献[17]中采用状态空间摄动法将J2项摄动的影响引入制导律来保证再入角和再入点的精度。

以往的制动策略规划研究中，往往忽略地球引力摄动对轨道的影响。一般基于二体假设和冲量假设，规划制动脉冲、制动角及制动起始点,使得估算的制动参数难以满足实际工程应用的高精度需求。为了保证再入点精度，往往需要引入末端修正等辅助手段。但实际工程应用中，制动发动机一般对于连续开关的时间间隔有约束，且轨道再入飞行器的轨控发动机一般只能提供单方向的推力，要满足高精度修正还需通过姿态机动来确保推力矢量方向。

基于以上问题和已有方法的不足，本文考虑地球非球形引力以及有限推力的实际情况，提出一种具备在轨自主决策、高精度落点状态控制的离轨制动方案。本方法的核心思路如下：

(1)首先将制动问题转化为轨道控制问题，将再入点纬度、速度以及倾角等需求转化为对离轨过渡轨道的半长轴、偏心率以及近地点幅角需求，从而便于利用轨道机动方程实现初值求解。

(2)在飞行器运行于离轨待命圆轨时利用机载计算机空闲时间，根据自主测定轨结果，通过数值递推轨迹，迭代求解制动参数，使得轨控完成后目标轨道根数误差满足精度需求。

(3)制动开始后，根据瞬平根转换的近似解析公式，预估再入点飞行状态。当估计的再入点满足约束条件时制动结束。

1 离轨制动问题

完整的轨道根数包含了对再入点轨道面内与轨道面外的所有约束，轨道面外的约束往往需要通过回归轨道的设计[19]或轨道相位的调整来使得再入点与离轨待命圆轨共面。轨道面外的调整问题需要根据燃料与时间的约束进行规划，离轨阶段一般仅考虑轨道面内的制动问题。因此，对制动点的约束可以简化为对再入点纬度、再入点速度及速度倾角的约束。再入点地心纬度可以表达为如下：

φE=arcsin[sin(iE)sin(uE)]

(1)

式中：uE为再入点纬度幅角，iE为再入点瞬时轨道倾角，φE为再入点地心纬度。

可见，在忽略轨道面外运动(即iE为固定值)的情况下，再入点地心纬度与再入点纬度幅角是一一对应的关系，对再入点纬度的需求可转化为对再入点纬度幅角的需求。

再入大气飞行阶段一般以确定的海拔高度作为转段点，地心距rE为确定值，可得再入点纬度幅角表达式：

(2)

再入点速度与速度倾角一般定义为地速及地速的倾角，而空间轨道是反映飞行器惯性系下的特性。由于地速与惯性速度存在确定的转换关系，因此，本文将速度及倾角的需求转化为惯性系下速度及倾角需求。轨道系下速度与航迹倾角分别为：

(3)

综上可以看出，忽略轨道面外运动，若飞行器所在的轨道满足再入点对aE，eE，ωE的需求，即可满足再入点的约束条件。由此离轨制动问题即可等效为轨道面内半长轴、偏心率及近地点幅角调整的变轨问题。

根据椭圆轨道运动理论，再入点与制动点满足能量守恒和动量矩守恒定理，可求得制动点的目标速度和目标倾角，给出脉冲假设下离轨制动的脉冲解,但并非任意情况下均有解。文献[7]中给出了一个确定地心距上界的公式，并提出当地心距大于临界地心距时，需要两次以上的冲量才能满足再入接口条件。从轨道控制的角度，可以给出几何意义直观的描述方法，即离轨待命轨道与离轨过渡轨道必须存在交点，才存在满足所有再入点约束的单脉冲解。当离轨待命轨道为理想圆轨时，可以简化为待命圆轨的地心距不能大于离轨过渡轨道的远地点，即：

(4)

但以上条件均是基于理想二体假设和脉冲假设而得到。实际工程中制动是一个有限推力持续作用的过程，不易求得严格的解析解。此外由于地球非球形引力使得轨道的瞬时根数在一个轨道周期内周期性波动，对于近地圆轨道卫星，半长轴波动可达到10 km以上，文献[13]中明确提出制动问题不能忽略J2项的影响。因此，本文在考虑地球非球形引力以及有限推力的前提下通过在线数值迭代求解制动参数。

2 制动参数在线确定方法

固定推力制动需要确定三个制动参数:制动角θ、制动速度增量ΔVset以及制动起始点的纬度幅角uopen。其中θ为推力矢量与飞行器速度矢量负向的夹角。可以看出，制动问题是一个多参数寻优问题。在地面规划中，多参数寻优的算法众多，本文不再一一赘述。但对于在线求解快速决策而言，解的收敛性、算法的可靠性以及算法复杂度的问题，严重制约了传统多参数寻优算法的应用。因此，本文将制动参数求解问题分解为多次单参数搜索问题，并且从工程实用的角度，在再入点速度不严格限制的情况下，将算法简化为单参数搜索问题。

2.1 基于平根数的制动参数初值求解

大量文献研究中以再入点速度、倾角以及位置为约束条件，基于能量守恒和动量矩守恒，求解了制动增量的解析解。本文基于轨道根数进行研究，因此首先给出一个基于平均轨道根数的制动初值求解算法。上文已经说明，对于制动问题，可转化为轨道面内半长轴、偏心率及近地点幅角调整的变轨问题。令待命轨道的平均根数为：a0,e0,ω0；过渡轨道的平均根数为：af,ef,ωf。本文对平均轨道根数的求解方式采用文献[9]的方法，在本文2.3节中给出了短周期项的一阶近似公式。

单脉冲有解的条件是两条轨道存在交点。交点处，两条共面轨道的纬度幅角u和地心距r一致，据此列写方程：

(5)

展开式(5),并令：

efx=efcos(ωf)，efy=efsin(ωf)，e0x=e0cos(ω0)，e0y=e0sin(ω0)。

式(5)整理为：

k1cos(uopen)+k2sin(uopen)=k3

(6)

式中:

根据式(6)可求得制动点纬度幅角。在式(6)有解的条件下，可以解出平均轨道情况下两条轨道理论的交点在轨道面内所在的位置uopen以及该位置处对应的地心距r。式(6)一般存在两个解，即过渡轨道与待命轨道存在两个交点。考虑到本公式只是提供一个可行的初值解，可选择离再入点近的交点作为理想制动点初值。

变轨前径向速度vor及法向速度vof分别为：

(7)

脉冲假设下，变轨后的径向速度vfr及法向速度vff分别为：

(8)

则制动脉冲大小为：

(9)

制动角大小为：

(10)

综上即可得到脉冲假设下制动角θ、制动点纬度幅角uopen、制动增量大小ΔVset。在制动弧段较小或精度要求低的轨道控制任务中，可直接利用脉冲假设的制动参数，将uopen设置在制动弧段的中心点。本节的计算仅提供制动参数的初值。

2.2 数值递推方程

确定飞行器的位置、速度初值和制动参数初值后，根据如下的微分方程进行制动轨道的数值递推：

(11)

为简化参数搜索的复杂度，本文在数值递推轨道过程中根据轨道根数实时预报再入点状态，通过关机时机确保每条弹道均满足再入角需求，从而避免对制动参数ΔVset的搜索。在Kepler轨道假设下，椭圆轨道未来任何时刻的位置、速度都可以通过简单的几何转换获得。但J2项的影响使得理想椭球近似下的轨道预报精度受到影响。文献[12]中给出了考虑J2项与理想Kepler轨道在再入点航迹倾角上的差异，可以看出不考虑J2项将严重影响再入点倾角精度。以下给出在考虑J2项影响下的再入点预报算法。

2.3 基于瞬平根转换的再入点状态预报

考虑J2项摄动的情况下，难以得到精确解析解，但在轨道控制理论中，轨道根数被分解为平均轨道根数和摄动项。摄动项分解为长周期项和短周期项，轨道根数预报公式可以按式(12)给出：

(12)

根据瞬平根转换的基本公式，近似解析的t时刻的瞬时轨道根数预报公式为：

(13)

式中：σ(t)为t时刻的瞬时轨道根数；σS(t)为t时刻的短周期摄动项。

本文主要参考文献[9]给出的一种无奇异点的平根数迭代法。采用的根数为σ=[a,i,Ω,ξ,η,λ]T。其中：a为轨道半长轴、i为轨道倾角、Ω为升交点赤经、ξ=ecosω，η=-esinω为轨道偏心率e在轨道面内的两个分量，λ=ω+M为近地点幅角ω和平近点角M之和。

表达式中采用了时刻t作为自变量，但对于短周期项，计算中使用的是纬度幅角，与绝对时间并无直接关系。在任一时刻，已知飞行器瞬时轨道根数时，可以通过有限次迭代计算出其短周期项，从而估计其平均轨道根数。文献[9]中给出了考虑偏心率高阶小量的短周期项计算公式，复杂度较高。对于近地轨道返回式飞行器，其离轨过渡轨道的偏心率一般较小，因此从工程实际的角度，本文给出仅考虑偏心率一阶小量的短周期项计算公式,见式(14)～式(19)。

(14)

(15)

(16)

(17)

(18)

(19)

式中：

(20)

轨道根数的长周期摄动，是随时间积分的量，平根数的一阶长周期项为[9]：

(21)

(22)

长周期项需要知道待飞时间，数值递推过程按下式估算待飞时间，从而补偿轨道根数的长周期摄动项。

(23)

根据以上平根数及长周期项，可以计算出再入点的平均轨道根数：

(24)

将式(24)得到的平均轨道根数，代入式(14)～式(19)可得预期再入点位置所对应的轨道根数短周期项σS(t)=[aS(t),iS(t),ΩS(t),ξS(t),ηS(t),λS(t)]T;进而得到预期再入点位置所对应的轨道瞬时根数σ(t)=[a(t),i(t),Ω(t),ξ(t),η(t),λ(t)]T。

根据无奇异点的瞬时根数，解算出预报的偏心率和近地点幅角:

(25)

2.4 参数迭代方案

数值递推过程中，制动开始后，首先判断当前轨道的近地点高度rp是否低于理想的再入点地心距rE。如果rp≤rE，则开始按照上一节给出的公式预报再入点的瞬时轨道根数。本文以再入角满足约束为停止制动的条件。根据轨道瞬根预报值，按下式估算再入点的真近点角：

(26)

图1 uopen迭代流程图

将真近点角f(t)=ω(t)+u(t)，代入式(3)得到再入点速度倾角γ(t)。当再入角满足γ(t)≤γE时，结束制动，并记录下累计制动速度增量ΔVset。

制动结束后，按惯性滑行继续数值递推动力学方程，至地心距满足|r|≤rE时，结束本次递推，并记录下递推结束时的纬度幅角u(t)及速度大小。

递推过程中以再入角约束确定制动增量，因此再入角满足约束，而再入点纬度幅角则与目标值存在偏差。因此一次递推结束后，根据再入点纬度幅角u(t)与理想再入点纬度幅角uE的偏差，按下式修正制动起始纬度幅角:

uopen←uopen-ku(u(t)-uE)

(27)

式中:ku为修正系数，取值范围为0.1～0.9。

由此可以看出，一次轨迹递推过程中，ΔVset是根据再入角满足条件确定，无需迭代求解。需要搜索的参数只有制动起始点纬度幅角uopen。综上，一次制动开始点的纬度幅角搜索流程如图1所示。

参考航天飞机返回任务，给出一个具体算例。假定飞行器位于345 km×345 km的圆轨道，离轨质量约95 t，配备两台6000 N推力器，设置目标再入点为海拔120 km,惯性速度7882.8 m/s，惯性倾角-1.1°纬度幅角-14.3°。参数迭代结果如图2～图4。

图2 制动开机点纬度幅角迭代曲线

图3 制动速度增量迭代曲线

图4 再入点纬度误差迭代曲线

仿真算例可以看出，5次迭代后再入点纬度幅角误差收敛至0.01度范围内，对应的纵向位置误差在1 km以内。由于迭代过程是以轨道根数预估的惯性倾角作为制动关机点，数值递推得到的再入点惯性倾角误差仅来源于预报误差。

图5 再入点惯性速度误差迭代曲线

图6 制动角与再入速度误差的关系

图7 制动角与制动速度增量的关系

2.5 落点速度与制动角的影响

制动问题概括为通过调整ΔVset，uopen以及θ使得再入点位置、速度以及航迹倾角满足uE，VE，rE的约束条件。上一节的参数迭代方案，通过满足uE,rE约束条件得到了ΔVset,uopen，但并没有保证再入速度VE满足约束。图5即为上一节迭代过程中再入点惯性速度误差变化曲线，仿真算例显示再入点速度误差未收敛至0附近。

由于上一节制动角θ采用2.1节中脉冲假设下解析计算的初值，在迭代过程没有调整θ，本节分析调整制动角θ对速度约束的影响。

设置不同的制动角，重复采用2.4节的参数迭代方案，每次迭代收敛后记录下速度误差与制动角，以及此制动角对应的制动速度增量。图6为选择不同制动角重复采用参数迭代后对应的再入点速度误差，图7为选择不同制动角所需的制动速度增量。

从图7可以看出，对于2.4节的算例，当θ=±27°时，经过迭代求解制动速度和制动起点，可以满足再入点速度误差到0附近。但此时对应的制动速度增量约为ΔVset≈97 m/s。对比图6，在θ=0°时制动速度增量ΔVset≈88 m/s，比θ=±27°可省约10%的燃料，而再入速度误差仅增加1.5 m/s。

可见为了保证所有再入点状态约束，可能需要付出极大的能源代价。实际工程应用中，对于再入飞行器，特别是大升力体再入飞行器，对再入点初始速度偏差具备较强的适应能力。航天飞机的离轨制动任务中，再入点约束条件是使再入点径向与切向速度满足一定的线性关系[20]。对于再入点速度的限制，一般可以在任务规划阶段通过设计合理的待命轨道来保证再入点速度在一定范围之内。因此，本文不考虑对制动角进行搜索，θ按2.1节给出的计算公式求解。

3 制动关机策略

理想的制动轨迹应该是纬度幅角及航迹倾角同时得到满足。但实际飞行中，由于推力误差、比冲误差，飞行器质量误差，姿态控制误差等因素，会使得真实的轨迹与预测轨迹存在一定的差异。实际轨迹很难保证再入点落点、速度、倾角同时满足需求。因此，需要考虑以何种指标决定制动关机时机。

对于轨道再入飞行器而言，最关键的再入点状态量为航迹倾角和落点位置。本文2.3节的公式给出了基于预报再入点瞬时根数求解再入点速度倾角的算法。制动过程，同样可以根据轨道预报的瞬时根数a(t),e(t),ω(t)，代入式(2)得到再入点对应的纬度幅角u(t)。可见，制动开机过程随时可以根据当前的飞行状态，计算出此时结束制动所对应的再入点纬度幅角及再入点速度倾角。

若不考虑制动过程中调整制动角，则制动一旦开始，控制系统只有惟一的控制量——制动关机时间。单一的控制量不能保证再入点纬度幅角和速度倾角同时满足约束。如果以再入点纬度幅角满足约束为条件结束制动，则再入点倾角会存在一定范围的散布；反之，如果以再入点倾角满足约束为条件结束制动，则再入点纬度幅角会存在一定范围的散布。飞行任务中，以何种指标决定制动关机时间，主要取决于该指标对偏差的敏感度。下一章中，针对2.4节提出的参考航天飞机再入返回轨道的算例，通过数值仿真对比两种不同的制动关机策略对应的再入点散布情况。

4 仿真分析

针对2.4节提出的算例，考虑如下的状态散布：

1)初始轨道半长轴：±3 km随机误差；

2)初始轨道平均偏心率：0.004随机误差；

3)制动过程姿态：±1°随机扰动；

4)飞行器质量：1%随机误差(-950～950 kg)；

5)制动推力：3%随机误差(-360～360 N)。

制动参数在线确定，闭环制动仿真过程分别以预报的再入点纬度幅角满足约束关机(策略1)和以预报的再入点航迹倾角满足约束关机(策略2)。分别进行200次蒙特卡洛打靶仿真，仿真结果如图8～图17。其中图8～图10为20条按策略1仿真的制动过程速度、高度和航迹倾角变化曲线，图11～图13为20条按策略2仿真的制动过程速度、高度和航迹倾角变化曲线。图14～图17为两种策略200组蒙特卡洛仿真的再入点纬度幅角和惯性倾角散布。不同的关机策略产生的误差散布范围如表1所示：

表1 不同关机策略的误差散布Table 1 Errors dispersion with different shutdown strategies

图8 策略1惯性速度变化曲线

图9 策略1高度变化曲线

图10 策略1惯性倾角变化曲线

图11 策略2惯性速度变化曲线

图12 策略2高度变化曲线

图13 策略2惯性倾角变化曲线

图14 策略1再入点纬度幅角偏差散布

图15 策略2再入点纬度幅角偏差散布

图16 策略1再入点惯性倾角偏差散布

图17 策略2再入点惯性倾角偏差散布

从仿真结果可以看出，策略1以预报纬度幅角满足约束条件为制动关机条件时，再入点纵向位置散布可以达到极高的精度(小于3.5 km)，误差源包含预报误差及主引擎关机后滑行过程中姿控推力器产生的轨迹扰动。策略2以预报再入点倾角满足约束为制动关机条件，再入点倾角误差散布极小(≤0.0015°)，但纬度幅角的散布则可达到120 km。

可见在制动角不控的工况下，通过实时预报再入点状态确定制动关机时机，至少可以保证再入点纬度幅角或再入点倾角中的一个量具有较高的精度。工程应用中，可根据实际任务对再入角和纬度幅角误差的容忍程度来确定制动关机的时机。

5 结 论

本文研究了一种在考虑地球非球形引力J2项影响下，基于瞬平根转换关系对再入点状态进行预报的自主离轨制动方法。算法首先根据再入点纬度和倾角的约束，通过在线数值迭代搜索制动参数，确定制动角、制动开机时间和制动速度增量大小。制动过程中通过近似解析的预报算法，计算再入点倾角和纬度幅角实时确定制动关机时刻。

在考虑离轨待命圆轨半长轴、偏心率存在初始散布，飞行器质量、推力、姿态均存在偏差的情况下，按照本文设计的离轨制动方案，以保证再入点航迹倾角约束和保证再入点纬度幅角约束的两种制动策略进行蒙特卡洛打靶仿真，仿真结果显示出再入点飞行状态精度较高。可以为大气层内再入飞行创造良好的初始条件。

离轨制动工程实践中还存在碰撞、敏感区域过境风险等与飞行任务顶层规划密切相关的问题，难以完全依靠在轨自主决策。实际应用中可通过实时下传自主规划的制动参数，通过地面设备对轨迹进行复核复算。当飞行轨迹确实存在风险时，可通过注入新的目标再入点状态或直接注入制动角调整轨迹，从而增强本方法在工程实践中的灵活性。