月球软着陆自主避障技术研究

2021-09-01陈少松牛津桥

载人航天 2021年4期

高峰，陈少松，牛津桥

(1.中国空间技术研究院，北京 100094； 2.中国人民解放军 63993 部队，北京 101149)

1 引言

随着深空探测技术的不断发展，人类对月球探测的目的发生了本质改变，由原来的飞行技术验证、地外科学实验向特定区域(月背/两极/阴影区)探测、表面大尺度巡视转变，这对着陆系统的精确性、强鲁棒、自适应性提出了更高的要求。

近几十年来，软着陆制导技术获得了长足的发展。从早期的苏联月球9 号采用的建立月球垂线法[1－2]到半开环半闭环的重力转弯制导方法[3－6]；从结合轨迹优化[7]与轨迹跟踪[8－11]的标称轨迹制导方法到能够在线应用的显示制导方法[12－14]以及最近兴起的智能学习制导方法[15－18]。从以上技术发展变化可以总结出软着陆制导技术的发展趋势呈现3 个阶段[19]:基于人为智能的给定轨道直接着陆、基于机器智能的轨迹优化以及基于人工智能的学习制导方法。

传统的给定轨道及优化思想实现的软着陆制导很大程度上依赖于先验知识，比如模型、约束、环境等，同时容易受建模偏差、不确定扰动等因素影响，且现有的障碍模型不具备较好的泛化能力。因此在月球这样多障碍、安全区域狭小和不确定性大的环境中，传统方法难以满足软着陆所需的的精度要求和安全性要求。近几年，凸优化技术逐渐成为研究的热点，基于凸优化的轨迹优化方法成为解决上述问题以及在线应用的最佳方案[20]。 Acikmese 等[21－22]提出了将非凸优化问题进行凸化转化的思想，解决了以发动机推力幅值的上下界作为约束条件的燃耗最优问题；Blackmore 等[23]的研究已经扩展到寻找最小的着陆误差轨迹的燃耗最优问题；Bai 等[24]提出了利用非凸函数对障碍物进行合并，将这些非凸障碍约束转化为凸障碍约束的线性化方法。

嫦娥工程中面对复杂月背、极区、阴影区等条件，安全着陆区域很小，本文针对精确自主避障技术开展深入研究，在凸锥模型的基础上建立半球模型和通用凸包模型来描述月面障碍物，利用凸优化技术，得到一种在多重障碍约束环境下的燃料最优的制导算法，并进行仿真实验。

2 软着陆问题描述

2.1 软着陆数学模型

采用三自由度模型描述着陆器运动，考虑以月面固联坐标系为基础进行建模[21]，ox方向垂直于月面指向外空间，oy沿月面指向月球北极，oz与其他2 个坐标轴构成空间直角坐标系，如图1所示。

图1 坐标系定义Fig.1 Definition of the coordinate

基于此坐标系，着陆器动力学可以描述为式(1):

式中，r∈IR3是着陆器的位置向量，v∈IR3是着陆器的速度向量，g∈IR3是月球的引力加速度向量，Tc∈IR3是净推力矢量，m是着陆器质量，α是燃料消耗率。

推力器认为是可变有限的，总推力受式(2)约束:

式中，T1、T2是推力的上下限。在着陆器下降过程中，需要满足式(3)的路径约束:

考虑边界条件、推力约束、着陆器倾角和表面高度，以燃料最优为优化指标，软着陆轨迹优化问题可描述为式(4)所示最优控制问题:

式中，t0为起始时刻，tf为终端时刻，m0为起始质量，r0，rf分别为起始位置和终端位置，v0，vf分别为起始速度和终端速度。

2.2 障碍物数学模型

行星表面的障碍物具有不同的形状和高度。Bai 等[24]将障碍物建模为高度可调、半锥角可调的凸锥，本文在此基础上将其扩展为半球和通用凸包模型。

在凸锥的情况下，着陆器不与障碍物碰撞的条件如式(5)所示:

式中，ne＝(1，0，0)T是单位向量，φ为半锥角，锥顶的位置向量为H＝(Hx，Hy，Hz)T，着陆器的位置向量为P＝(rx，ry，rz)T。

除了将障碍物描述为凸锥，还可以将其描述为包围障碍物的以其几何中心点S＝(Sx，Sy，Sz)为球心、RS为半径最小外接半球，如图2 所示。

图2 半球障碍模型Fig.2 Hemisphere obstacle model

在这种情况下，只要着陆器在半球的半径范围外就是安全可行的，那么着陆器不与障碍物碰撞的约束条件如式(6)所示:

将障碍物建模为凸锥或者半球，虽然处理起来非常方便，然而缺少通用性，不能反映不同类型障碍物的几何特征，故需要一种统一的建模方式，即通过传感器获得的点云数据直接建立三维凸包。本文使用随机增量式算法来构造点集T的凸包。首先定义凸包可见区域和地平线，如图3所示。

图3 凸包可见区域和地平线Fig.3 Visible area and horizon of convex hull

假设Ti落在凸包lh(Ti－1) 之外，在点Ti的位置朝lh(Ti－1) 看去，在lh(Ti－1) 的表面上，那些可见的小平面连在一起组成的连通区域称之为点Ti在凸包lh(Ti－1) 上的可见区域。而围成这个区域的lh(Ti－1) 上的边组成的折线称之为点Ti在凸包lh(Ti－1) 上的地平线[25]。具体算法过程如图4 所示。

图4 凸包生成算法Fig.4 Generation algorithm of convex hull

容易证明在三维空间中，任意n个点的凸包的复杂度为o(n)[25]，得到凸包以及凸包所包含的点集后，给出建立障碍约束的方法。假设最终得到的凸包由NT个凸顶点组成，定义函数f(P)＝∑‖P(t)－Ti‖/NT，其中Ti是凸顶点的位置向量。该函数反映了点与凸包之间的相对距离关系，是IRn上的凸函数，也就是说只要函数大于凸包边界点函数的最大值，则点在凸包外。为了更均匀地反映函数与点的位置关系，可将最大值用平均值替代，则着陆器不与障碍物相撞的条件如式(7)所示:

式中，NT为凸包障碍点集个数。综上，将一般性的月球软着陆避障问题描述为式(8)所示:

式中，H为障碍锥锥顶的位置向量，S为半球障碍球心的位置向量，T为凸点集的位置向量，ϑ是半锥角φ的补角。

3 基于凸优化的轨迹优化方法

3.1 凸化转化

2.2 节式(8)的非凸性会造成计算求解困难，Acikmese 等[21]给出了将非凸约束转化为凸约束的方法，令u＝Tc/m，σ＝Γ/m，z＝ln(m) ，其中Γ为推力松弛向量，则式(4)中的非凸约束将被转化为式(9):

然而线性函数中省略了f(P) 的泰勒展开式中的二阶项，令向量(Pξ－Q) 和向量(P－P0)之间的夹角为l，容易求得f(P) 和其二阶余项的相对误差如式(12)所示:

通过在节点处连续线性化障碍约束，式(13)便可以表示为式(14)的形式:

式中，s为松弛因子。根据式(12)可知f(P)的二阶余项为正，这使得约束更加严格。图5 可以看出，着陆器的实际可飞行空间区域大于理论可飞行空间区域，故计算所得解是有效的。

图5 障碍转化前后关系Fig.5 Relationship before and after barrier transformation

线性化后的问题由于缩小了解的范围有可能导致问题无解，故引入松弛因子来适当放宽约束条件，以在避障和最优解之间取得平衡。

3.2 二阶锥优化

在整个控制时域中，每个离散时间点的状态满足式(17):

4 实验仿真与结果

4.1 仿真条件

本文算法通过Matlab 和Yalmip 联合求解，初始参数选择如表1 所示。

表1 初始参数Table 1 Initial parameters

4.2 仿真验证

本文在障碍物分别为凸锥、半球和凸包的多障碍环境下对凸优化算法进行验证。

4.2.1 凸锥多障碍

选择3 个障碍物作为约束条件，并且这3 个障碍物最好能够排成1 个平面以凸显避障的效果，其参数如表2 所示。

表2 凸锥障碍参数Table 2 Parameters of convex hull obstacle

仿真结果如图6 所示，多重障碍物之间明显存在可行路径，但是安全区域较小。从图6 可以看出着陆器连续绕开了3 个具有不同半锥角的障碍物，着陆位置为[0.0697，0.1070，0.9918]m，着陆速度为[－0.5184，－0.1682，－0.1468]m/s，消耗燃料384.4954 kg，得到了全局最优解，说明面对多重障碍时凸优化算法有效。

4.2.2 半球多障碍

相对于凸锥障碍而言，多重半球障碍环境下着陆器的可行安全飞行空间更窄，通常情况下需要穿越2 个相邻半球之间的缝隙，并且经常需要采取爬升的策略进行避障。半球障碍参数如表3所示，其仿真结果如图7 所示，可见飞行轨迹穿越黄色和橙色半球之间的缝隙，绕开绿色半球从而到达目标点；着陆位置为[0.0697， 0.1071，0.9918] m，着陆速度为[－0.4482， 0.1854，－0.121]m/s，消耗燃料286.8325 kg。

图7 半球仿真结果Fig.7 Results of hemisphere simulation

4.2.3 凸包多障碍

点云图是由计算机随机生成的，以模拟传感器所扫描的障碍的点云数据，如图8(a)所示，由点云生成凸包的结果如图8(b)所示。

图8 凸包生成结果Fig.8 Results of convex hull generation

得到的凸包由2 个不规则多面体组成，相比前2 种模型而言，凸包更能反映障碍物的几何特征，可飞行区域也是最精确的。仿真结果如图9所示，着陆器为了穿越狭窄的空间，首先采取了拉低的策略穿过低矮的峡谷，而后面对斜挡的障碍采取爬升的策略飞掠障碍表面，在极小的飞行空间内找到了一条安全可行的轨迹；着陆位置为[0.0697， 0.107， 0.9918] m，着陆速度为[－0.1762，－0.4675，－0.1171] m/s，消耗燃料361.8431 kg。