葛丹桐,崔平远,高 艾
(1. 北京理工大学宇航学院,北京100081;2. 深空自主导航与控制工信部重点实验室,北京100081;3. 飞行器动力学与控制教育部重点实验室,北京100081)
火星安全着陆轨迹快速生成的能控集法
葛丹桐1,2,3,崔平远1,2,3,高 艾1,2,3
(1. 北京理工大学宇航学院,北京100081;2. 深空自主导航与控制工信部重点实验室,北京100081;3. 飞行器动力学与控制教育部重点实验室,北京100081)
针对火星着陆器在动力下降过程可能存在的着陆轨迹实时修正和着陆点再选择等问题,本文通过离线计算着陆器在动力下降段开始时的能控集范围,并在实际着陆过程中判断着陆器当前状态与能控集间的关系,实时确定最终着陆点并快速搜索相应着陆轨迹。若着陆器无法到达预定着陆点,则在视野范围内根据安全因子重新选择着陆点规划运动轨迹。仿真表明,基于能控集的快速轨迹规划法可根据着陆器的实际初始状态快速确定最终着陆点并获得相应着陆轨迹,以有限的燃耗实现火星安全软着陆的目标。
火星动力下降;安全着陆;能控集;轨迹生成;安全因子
火星着陆探测是我国未来深空探测的主要目标,尤其是实现在具有科学价值地区的精确软着陆,然而这样的地区大多环境复杂、地形崎岖[1-2],同时在整个着陆过程中还存在着系统累积误差与环境扰动[3],所以对GNC系统提出了很高的要求。动力下降段作为着陆器着陆火星的最后阶段,除了将下降速度减至零外,还需尽可能修正开伞点误差以及伞漂所造成的与预定轨迹的位置偏差[4],以减小最终的着陆误差,这一过程主要通过着陆器的横向机动来实现[5]。由于实际过程中着陆器相对于标称轨迹发生的偏移量难以提前准确估计,在转移至目标着陆点时消耗的燃料有超出工程约束的风险[6]。因此在动力下降段开始时,需要结合着陆器当前状态及能控范围,确定最终着陆点[7]并快速搜索出相应着陆轨迹。当预定着陆点超出着陆器能控范围时,则需在可见范围内重新选取着陆点并在线规划轨迹,以提高任务的成功率。
以上方法涉及到对着陆器能控集的离线计算与分析,能控集指能够到达给定末端状态的初始状态集合,其概念最早产生于数学系统理论中,近些年,一些学者结合实际任务特点对其进行了改进并将之用于行星着陆探测研究中[8-9]。考虑到计算的复杂性,传统的能控集分析多用于前期任务设计,并未考虑其在实际任务执行过程中的作用。同时,受制于星载机的计算与存储能力,目前能够在线生成的轨迹大多由形式简单的解析制导律产生,如多项式制导[10]、能量最优制导[11]等,这些方法出于实时计算需求,只能满足基本的约束如初始末端状态约束,而在提前预判是否超出燃料上限以及提高安全着陆系统的灵活性等方面还有待进一步提高。
本文在分析速度增量的基础上,得到火星动力下降段位置及速度能控集的快速生成方法,并将结果进一步用于下降过程中着陆轨迹的快速生成,为安全着陆的实现提供参考,其具体应用示意图如图1所示。
通过离线分析能控集,星载计算机可在实际下降过程中根据着陆器当前状态与标称值间的偏差大小,判断其是否在提前存储好的可控范围内,实现以一定精度到达预定着陆点的目标,从而确定任务的最终着陆点,并进行轨迹的在线快速搜索。在极端情况下,假设目标着陆点及备用点均无法满足需求,超出着陆器的控制范围,则将重新在线规划着陆点与着陆轨迹使任务得以继续进行,从而提高着陆的安全性,实现行星表面的软着陆。
假设采用四条着陆腿着陆器,不能转向的发动机固连在着陆器机体外侧,指令推力矢量的方向通过改变着陆器姿态来实现,在理想情况下,发动机推力矢量方向能够迅速机动到与制导加速度指令一致的方向上去,此外,由于目前对火星着陆器尚无滚转要求,因此忽略推力矢量在滚转方向的动力学。整个动力下降段动力学方程如下所示[12]
(1)
式中:(x,y,z)T为位置矢量,(u,v,w)T为速度矢量,Γ为比冲的大小,姿态角θ和ψ如图2所示,表示了推力矢量在惯性系下的夹角。
采用好奇号任务在动力下降段使用的多项式制导方法,假设每个方向上的加速度都是时间的二次多项式[10]
a(t)=C0+C1t+C2t2
(2)
通过积分可得到速度和位置关于时间的函数。再结合初始和末端状态约束即可解得相应系数。在下降过程中通过不断测量当前状态,并将测量结果输入给制导模块作为初始状态,即可得到闭环的多项式制导律,其系数为
(3)
为了进一步得到下降时间,假设竖直方向加速度为关于时间的线性函数,即令C2z为0,同时假设rzf=vzf=azf=0,得到剩余时间表达式
(4)
多项式制导形式简单、计算方便,在需要大量存储离线轨迹时,只需要存储多项式的系数即可,免去了记录每个节点状态的巨大内存需求,为轨迹的在线快速搜索与生成提供了便利。参考“好奇号”任务,仿真所用到的着陆器相关参数如表1所示。
表1 仿真参数Table 1 Simulation parameters
速度增量ΔV的大小表征了下降过程中着陆器为了减速所消耗的燃料多少。本节针对多项式制导方法,对速度增量表达式进行分析,得到着陆在不同位置的燃耗分布规律,从而避免了逐点积分计算着陆轨迹对应的燃耗,仅通过判断着陆点所在区域即可快速得到相应轨迹的速度增量大小,为下降段能控集分析打下基础。
给定目标着陆点,下降过程的速度增量可由下式计算得到[13]
(5)
式中:
(6)
(7)
即可得到两个水平方向速度极值点对应的时刻,除了剩余时间tgo外,另一个解为
(8)
在每个水平方向,如果0 (9) (10) 由于tgo>0,故 (11) (12) (13) 联立以上结果,关于时间的判断条件转化为关于空间的表达式 当vx0>0时, 当vx0≤0时, (14) 当vy0>0时, 当vy0≤0时, (15) 区域1和其他区域最大的区别在于该区域内的速度增量仅仅只取决于初始状态,而在区域2~4中,速度增量是初始状态和末端着陆点位置共同作用的结果。同时,随着着陆点与中心点距离的增大,速度增量呈现出非线性变化特点。区域1的中心点为其边界的平均值,即(rx0+1/3vx0tgo,ry0+1/3vy0tgo)。如图4所示,在x-z平面存在以下几何关系 (16) 代入式(4),上式可进一步转化为 (17) 类似结果也可在y-z平面内得到。可以看出,区域1的中心点正是初始速度矢量与地表平面的交点,将中心点设为最终着陆点,着陆器将沿着初始位置与末端位置的连线飞行。 一般来说,任务在设计时会留有一定余量,使得着陆器在紧急情况如避障时仍有剩余燃料来实现机动[14],因此按照预定轨迹运动的着陆器所消耗的燃料便可以控制在合理范围内。然而,由于伞降段无控制执行机构,开伞点误差会经传递而不断累积[15],当动力下降段开始时,着陆器的实际状态可能会与预定值出现极大的偏差,为了实现精确软着陆,着陆器需要在减速的同时横向转移很大的距离,以尽可能小的误差着陆在目标着陆点附近,这一过程往往需要消耗大量的燃料,在极端情况下甚至有超出约束范围的可能。因此分析动力下降段开始时的能控集,明确着陆器在有限燃耗下能够到达给定目标着陆点的位置速度范围,以进行在线着陆点确定和轨迹搜索或规划就显得尤为必要。 对于软着陆任务,着陆器要实现的首要目标是以零速度安全降落在火星表面,因此末端着陆速度均设为0,根据表1中着陆器参数,得到所允许的最大速度增量[16]: (18) 着陆器若要实现在给定的燃料质量下以一定精度到达目标着陆点,动力下降段的初始位置不能距离着陆点太远,速度也不能过大,为了定量描述能够到达目标点的允许状态范围,本节分别针对动力下降段考察着陆器的位置能控集和速度能控集,并对于每个可行的状态讨论其对应轨迹的存储方式。 表2 不同区域速度极值及速度增量表达式Table 2 Velocity extremum and velocity increment of different areas 3.1 能控集计算 在分析位置能控集时,仅考虑初始位置变化对速度增量的影响,固定初始速度与末端状态为 v0=[20,20,-80]Tm/s 当初始位置改变时,对应轨迹的速度增量也会发生改变。根据上文的推导结果,针对每一个初始位置判断目标着陆点位于哪个区域,并采用表2中相应区域的速度增量表达式对该轨迹对应的ΔV进行求解。将高度变化范围设置为20m到2000m,对每一个高度分别进行计算,保留最大速度增量350m/s以内的区域,得到的初始位置能控集如图5所示,图中随着高度的降低,位置能控集呈现出先增大后减小的趋势,在大约1500m处达到最大。然而动力下降段的主要目的在于减速,其初始高度并不能无限降低到地表,这主要受发动机能够提供最大推力的限制,过低的初始位置将严重影响有效减速高度,导致着陆器即使在燃耗范围内始终以最大推力飞行也无法实现零速度着陆。考虑到多项式制导在竖直方向的加速度是关于时间的一次函数且末端值为0,因此其最大值出现在下降段的最开始,即 azmax=az0=az(t=0)=C0z= (19) (20) 因此最小的可行高度可由下式计算得到 (21) 类似方法同样可以用于速度空间,得到相应的动力下降段初始速度可行域与沿标称轨迹下降过程中的速度允许偏离范围。具体来说,为了得到初始速度可行域,首先固定初始位置与末端状态 r0=[-500,-800,1500]Tm 接着针对每一个初始速度判断目标着陆点位于哪个区域,并用相应区域的速度增量表达式对ΔV进行求解。将下降速度变化范围设置为-50m/s到-100m/s,对每个下降速度分别进行计算,保留最大速度增量350m/s以内的区域,得到的初始速度能控集如图6所示,图中随着下降速度的减小,速度能控集不断缩小。 3.2 轨迹簇离线存储 类似地,当有多个备选着陆点时,可依次得到每个着陆点对应的能控集范围,这些能控集的并集表明了在动力下降段开始时由系统误差及外界干扰造成的着陆器相对于标称状态的允许偏离范围。能控并集中每一个状态量对应一条满足工程约束的着陆轨迹,通过离线计算并存储这些轨迹,可实现着陆器星上轨迹实时快速搜索,从而减少了在线计算量与生成轨迹所需要的时间。同时,由于多项式制导形式简单,位置与速度均为关于时间的多项式函数,因此只需要存储多项式的系数即可得到相应的着陆轨迹,对内存的需求也大大降低,对于非多项式形式的制导方法,可采取存储轨迹节点的方式以降低存储需求。 不失一般性,以下假设初始速度已知,仅考虑位置能控并集及其对应的轨迹簇。假设提前选取的三个预定着陆点X,Y,Z,其坐标分别为rX=[0,0,0]Tm、rY=[1000,500,0]Tm、rZ=[-1000,-1000,0]Tm,考虑科学价值、工程约束等因素,其按优先级排序为X→Y→Z,即X为目标着陆点,Y、Z为备选点。根据第3.1节分别计算着陆点X、Y、Z在高度为1500m处对应的横向初始位置能控集,得到的结果如图7所示。 注意到三个着陆点在同一高度的初始位置能控集出现了交叉重叠的情况,当着陆器实际位置落在重叠区域时,根据优先级顺序选取着陆点。将能控并集每个点的位置信息与对应轨迹的多项式系数以及相应着陆点优先级顺序按照以下方式存储在矩阵M中 在线判断着陆器当前状态与提前存储的能控集间的关系存在两种结果:若着陆器状态位于能控集范围内,意味着着陆器能够利用有限燃耗到达预定的目标着陆点,此时仅需要搜索出符合当前状态的轨迹参数,即可得到下降轨迹;若着陆器无法到达预定着陆点,即着陆器状态由于系统误差或外界扰动影响超出能控集时,系统需要根据当前情况实时在线规划出新的着陆点与着陆轨迹,在不违背工程约束的前提下实现着陆器的安全着陆。 4.1 在能控集范围内 Ck0=[1.6593,7.3481,2.8444]T 则相应的制导律与轨迹为 (22) 相应的加速度及位置变化曲线如图8~9所示,着陆器通过在线控制推力器推力大小,使其沿生成的轨迹运动,即可保证在有限燃耗的情况下到达目标着陆点。 4.2 超出能控集范围 当着陆系统根据敏感器信息判断出着陆器位于能控集范围外,意味着即使消耗尽所有的燃料也无法将着陆器送至预定目标点时,系统需要对可行的着陆点进行快速重新评估,并根据当前状态与新选取出的着陆点在线生成着陆轨迹。 假设下降段初始位置r0=[3000,-800,1500]Tm,由图7可知,此时无论到达哪个目标着陆点消耗的燃料都将超出着陆器的最大转移能力,因此需要在线对视野范围内的地形进行快速评估,使得着陆器利用有限燃耗降落在安全的地方[17]。考虑在以目标着陆点为中心的6000m×6000m地形范围内,结合着陆器初始状态,分别以每个像素点作为着陆点计算各点速度增量,剔除高于最大速度增量350 m/s的区域,对剩余的地形区域进行测量筛选,利用敏感器信息重建地形,得到行星表面DEM信息,如图10所示,其中深色部分为视野范围内不可达区域,浅色部分为提取的可行区域。 根据得到的DEM信息计算可行区域内各点的坡度φ和表面粗糙度dl[18],结合着陆器能够容忍的障碍尺寸如岩石高度、坡度大小等约束,利用文献[19]提出的安全因子概念以及基于安全因子的着陆点评估方法,通过以下简化的地形安全评分标准对该区域各点进行安全性评估,从而得到最适宜着陆的地区 (23) 假设ηφ=ηd=0.5,φsafe=15°,dsafe=0.5m,不失一般性,不可达区域的安全评分此处取S(xi,yi)=20,将计算结果按安全评分高低进行排序,评分越低意味着对应区域越安全,得到的新着陆点rf=[1730,-1870,0]Tm(即全局最小值所在位置)如图11所示。此后着陆器再根据当前状态与新的末端状态利用多项式制导重新生成轨迹,得到多项式系数 C0=[-6.9499,-6.1914,2.8444]T 相应的加速度及位置变化曲线如图12~13所示。同理,当着陆器沿该轨迹运动时,可实现有限燃耗下的安全着陆。 另外,尽管传统的多项式制导可在下降过程中实现实时燃耗估计并更新着陆点,但这一过程往往涉及到对着陆器的机动范围进行在线分析与评估,对星载机的计算速率提出了很高的要求;而本文提出的方法其主要优势体现在可提前离线完成着陆器到达不同备选着陆点的燃耗评估预测,以尽可能减少在线运算量为目标,给出了一套快速搜索生成轨迹的方法,从而在不增加星载机计算量的前提下提高了任务的安全性。 火星动力下降过程中着陆器为了消除横向偏差,往往需要进行大范围的横向转移,这对其所携带的有限燃料提出了挑战。为提高火星着陆的安全性,本文在离线分析着陆器到达不同着陆点的速度增量分布规律的基础上,结合工程参数,提前计算出动力下降段的初始状态能控集范围,并将其与对应的轨迹参数进行存储。在实际飞行过程中,着陆器通过判断当前状态与存储的能控集间的关系以确定最终着陆点及相应下降轨迹:当着陆器状态位于能控集中时,通过在提前存储的矩阵中进行一维线性搜索得到最终着陆点坐标及相应轨迹多项式系数;当着陆器状态位于能控集外时,则需要在线根据地形安全评价标准——安全因子重新选取着陆点,快速进行轨迹重规划,以保证着陆器在有限燃耗下着陆在安全的地形表面。仿真结果表明,基于能控集的快速轨迹规划法能够在下降过程中快速确定最终着陆点并生成下降轨迹,可实现着陆器在初始状态无法提前准确预知情况下的安全着陆,从而提高了着陆器的自主性与任务的成功率。 本文以多项式制导为例,给出了轨迹快速生成能控集法的具体应用。相比于更为复杂的难以实现实时在线计算的制导方法,本文方法对在线确定着陆点并快速生成轨迹将更具有优越性。未来可围绕更为先进的制导方法(如凸优化等)研究不同制导律下的轨迹快速生成策略,为安全着陆系统轨迹设计提供参考。 [1] 崔平远, 胡海静, 朱圣英. 火星精确着陆制导问题分析与展望[J]. 宇航学报, 2014, 35(3):245-253. [Cui Ping-yuan, Hu Hai-jing, Zhu Sheng-ying. Analysis and prospect of guidance aspects for Mars precision landing[J]. Journal of Astronautics, 2014, 35(3):245-253.] [2] Epp C D, Robertson E A, Carson J M. Real-time hazard detection and avoidance demonstration for a planetary lander[C]. AIAA SPACE 2014 Conference and Exposition, San Diego, USA, August 4-7, 2014. [3] 于正湜, 崔平远. 行星着陆自主导航与制导控制研究现状与趋势[J]. 深空探测学报, 2016, 3(4):345-355. [Yu Zheng-shi, Cui Ping-yuan. Research status and developing trend of the autonomous navigation, guidance, and control for planetary landing[J]. Journal of Deep Space Exploration, 2016, 3(4): 345-355.] [4] 任高峰, 高艾, 崔平远,等. 一种燃料最省的火星精确着陆动力下降段快速轨迹优化方法[J]. 宇航学报, 2014, 35(12):1350-1358. [Ren Gao-feng, Gao Ai, Cui Ping-yuan,et al. A rapid power descent phase trajectory optimization method with minimum fuel consumption for Mars pinpoint landing[J]. Journal of Astronautics, 2014, 35(12):1350-1358.] [5] Harris M W, Aç1kmee B. Maximum divert for planetary landing using convex optimization[J]. Journal of Optimization Theory and Applications, 2014, 162(3):975-995. [6] Dueri D, Aç1kmee B, Scharf D P, et al. Customized real-time interior-point methods for onboard powered-descent guidance[J]. Journal of Guidance Control & Dynamics, 2017, 40(2):197-212. [7] 田阳, 崔平远, 崔祜涛. 基于图像的着陆点评估及着陆器运动估计方法[J]. 宇航学报, 2010, 31(1):98-103. [Tian Yang, Cui Ping-yuan, Cui Hu-tao. Landing site assessment and probe motion estimation based on image[J]. Journal of Astronautics, 2010, 31(1):98-103.] [8] Eren U, Dueri D, Aç1kmee B. Constrained reachability and controllability sets for planetary precision landing via convex optimization[J]. Journal of Guidance Control & Dynamics, 2015, 38(11):1-17. [9] Long J, Gao A, Cui P. Controllable set analysis for planetary landing under model uncertainties[J]. Advances in Space Research, 2015, 33(2):281-292. [10] Steinfeldt B A, Grant M J, Matz D A, et al. Guidance, navigation, and control system performance trades for Mars pinpoint landing[J]. Journal of Spacecraft & Rockets, 2012, 47(1):188-198. [11] D’Souza C. An optimal guidance law for planetary landing[C]. AIAA Guidance, Navigation, and Control Conference, 1997. [12] Topcu U, Casoliva J, Mease K D. Minimum-fuel powered descent for Mars pinpoint landing[J]. Journal of Spacecraft & Rockets, 2007, 44(2):324-331. [13] Wong E C, Singh G, Masciarelli J P. Guidance and control design for hazard avoidance and safe landing on Mars[J]. Journal of Spacecraft & Rockets, 2015, 43(2):378-384. [14] 王大轶, 李骥, 黄翔宇,等. 月球软着陆过程高精度自主导航避障方法[J]. 深空探测学报, 2014, 1(1):44-51. [Wang Da-Ji, Li Ji, Huang Xiang-yu, et al. A pinpoint autonomous navigation and hazard avoidance method for lunar soft landing[J]. Journal of Deep Space Exploration, 2014, 1(1):44-51.] [15] Quadrelli M B, Wood L J, Riedel J E, et al. Guidance navigation and control technology assessment for future planetary science missions[J]. Journal of Guidance Control & Dynamics, 2015, 38(7):1165-1186. [16] Ploen S R, Seraji H, Kinney C E. Determination of Spacecraft landing footprint for safe planetary landing[J]. Aerospace & Electronic Systems IEEE Transactions on, 2009, 45(1):3-16. [17] 崔平远, 葛丹桐. 一种行星安全着陆点综合评估方法[J]. 深空探测学报, 2016, 3(4):363-369. [Cui Ping-yuan, Ge Dan-tong. An integrated evaluation of planetary safe landing site[J]. Journal of Deep Space Exploration, 2016, 3(4):363-369.] [18] 吴伟仁, 王大轶, 黄翔宇,等. 月球软着陆自主障碍识别与避障制导方法[J]. 中国科学:信息科学, 2015, 45(8):1046-1059. [Wu Wei-ren, Wang Da-yi, Huang Xiang-yu, et al. Autonomous hazard detection and avoidance guidance method for soft lunar landing[J]. Scientia Sinica Informationis, 2015, 45(8):1046-1059.] [19] Cui P, Ge D, Gao A. Optimal landing site selection based on safety index during planetary descent[J]. Acta Astronautica, 2017, 132: 326-336. 通信地址:北京市海淀区中关村南大街5号北京理工大学宇航学院(100081) 电话:(010)68918910 E-mail: gedt@bit.edu.cn 崔平远(1961-),男,教授,博士生导师,主要从事飞行器自主导航与控制、深空着陆器自主技术与轨道设计。本文通信作者。 通信地址:北京市海淀区中关村南大街5号北京理工大学宇航学院(100081) 电话:(010)68918611 E-mail: cuipy@bit.edu.cn (编辑:牛苗苗) Rapid Generation of Mars Safe Landing Trajectory Based on Reachability Set GE Dan-tong1,2,3, CUI Ping-yuan1,2,3, GAO Ai1,2,3 (1. School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China; 2. Key Laboratory of Autonomous Navigation and Control for Deep Space Exploration, Ministry of Industry and Information Technology, Beijing 100081, China; 3. Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, Beijing 100081, China) To solve the problems of possible real-time trajectory correcting and landing site re-designating during Mars powered descent, the paper computes the reachability set at the beginning of the powered descent phase offline and determines the final landing site online by comparing the current state of the lander with the obtained reachability set. The corresponding descent trajectory is then rapidly searched. If the lander fails to reach any of the given targets, a new landing site in the field of view will be selected through safety index and a trajectory will be planned. The simulation result shows that the proposed reachability set-based rapid trajectory planning method manages to determine the final landing site and obtain the landing trajectory in real-time according to the practical initial state of the lander, fulfilling the goal of Mars soft landing with limited fuel consumption. Mars powered descent; Safe landing; Reachability set; Trajectory generation; Safety index 2017-03-06; 2017-03-31 国家自然科学基金(61374216, 61304226, 61304248, 61603039);中国博士后科学基金(2016M591087) V448.2 A 1000-1328(2017)05-0497-09 10.3873/j.issn.1000-1328.2017.05.008 葛丹桐(1992-),女,博士生,主要从事行星探测制导与控制、安全着陆与障碍规避。3 能控集分析及离线轨迹存储
rf=[0,0,0]Tm
vf=[0,0,0]Tm/s
rf=[0,0,0]Tm
vf=[0,0,0]Tm/s4 轨迹在线快速生成方法
Ck1=[-0.1559,-0.5605,-0.0506]T
Ck2=[0.0022,0.0076,0]T
C1=[0.4563,0.4024,-0.0506]T
C2=[-0.0059,-0.0052,0]T5 结 论