一种探测器自主安全着陆区确定及精确避障的策略

金 涛 章登义 蔡 波

1(武汉大学计算机学院 武汉 430072)2(武汉大学国家网络安全学院 武汉 430072)(792042768@qq.com)

发射星际探测器并实现在目标星球表面的软着陆是我国探索深空、发展航天技术的重点需求.自主安全着陆是未来着陆器进行安全着陆的重要能力，这种能力主要体现在着陆器能够对传感器采集的信息进行快速处理、判断并自动形成新的着陆指令.为了达到预定的科学目标，地面预先选择的软着陆区并不是完全平坦的，它存在着一定分布的石块或陨石坑.这些石块和陨石坑的存在会给探测器安全着陆带来危险，因此，着陆器需要具有发现和识别障碍并进行机动避障的能力.为了实现下落过程中对障碍物的躲避,需要使用探测器上携带的电荷耦合器件(charge coupled device,CCD)相机进行拍照，根据拍摄的月面图像来判断月表的地形.

至今为止,各国实施深空探测任务的着陆器还不具有对着陆区的潜在危险进行检测和规避的闭环控制能力.为了使深空探测器具有自主软着陆及避障的能力，喷气推进实验室(Jet Propulsion Laboratory,JPL)以未来火星和小天体着陆任务为背景，开发了基于光学敏感器的精确导航和障碍规避算法[1-3].欧洲航天局(European Space Agency,ESA)也开发了具有高精度着陆和自主避障能力的新一代着陆系统，以探索科学价值高且危险的区域[4-5].从当前的研究现状可以看出，基于灰度图像的避障技术和基于激光雷达(light detection and ranging,LIDAR)的避障技术已经被考虑用于探测器自主障碍中，并且得到了进一步发展[6].文献[7]综合使用多因子分析评估了月球探测的着陆区域和巡航路线.

1 相关工作

目标天体的地表障碍主要是岩石、陨石坑等.对于岩石的检测，Ma等人[8]提出了边缘光流算法，该算法是基于纹理的图像分割，光照影响较大；Matthies等人[9]通过立体视觉三维重构来进行障碍检测，该方法较为准确可靠，可以检测的障碍类型较多，但是涉及的数据计算量大，导致检测速度较慢；Huertas等人[10]提出基于模版匹配的自主障碍检测方法，这种方法的有效性取决于模板建模的精确程度，模板建模越精确，但匹配所需要的时间越长.对于最大安全区域的计算及选择并没有相应的算法；而文献[11]也只是对安全着陆圆进行了粗略的估计.为了达到预定的科学目标，初避障通过图像匹配识别所选定的初选着陆区域会存在一定的误差，导致探测器安全着陆存在一定的风险，因此，探测器需要有更高精度的自主障碍识别与规避能力.为了实现探测器的高精度避障，需要对初选着陆区域进行高精度三维地形重构，然后根据重构后的模型，进行螺旋搜索，根据探测器安全着陆阈值来选定目标着陆点.

在精避障阶段，高精度的障碍识别算法还未能成熟的应用于航天探测及仿真项目中.南加州大学的Leroy等人[12-14]都对图像匹配与识别的方法进行了研究与论述.在图像匹配与识别方法之后，激光测距仪的思想被学者们提出，即通过仪器获得探测器到目标点之间的距离，然后对采样值进行处理以获得目标星球表面信息，从而提取出星表特征进行障碍检测.这一思想最初是由JPL实验室的Andrew等人[15]提出的，之后JPL实验室又对该思想进行了进一步的研究[16]，设定了探测器障碍检测的相关参数和技术指标.之后美国 JPL实验室、MIT Lincoln实验室以及美国、英国和加拿大的多所大学开始对多传感器协同工作避障的方法展开研究.

本文根据中国深空探测的现状，提出了一种基于CCD着陆相机生成表面遥感影像并进行着陆区岩石检测与安全着陆区确定的方法.这种方法不仅提出了遥感影像的生成方法，同时对着陆区岩石的检测、识别进行了研究，并且提出了确定最大安全着陆区的算法.在精避障阶段本文提出了一种基于三维重构的探测器精避障螺旋扫描方法.这种方法可以通过激光三维成像仪获取的数据，构造高精度的表面初选着陆区域三维网格模型，并针对模型的分块进行高精度螺旋扫描，结合探测器安全着陆阈值选定目标着陆点.其整个流程构成了一个完整的探测器软着陆精避障过程实现的解决方法.

Fig.1 The technical roadmap of autonomous safe landing simulation algorithm图1 自主安全着陆仿真算法技术路线

2 系统模型

软着陆的最后垂直着陆段将进行障碍物的检测与规避，由于对垂直着陆段的高度有一定的约束，因而希望能够快速地生成遥感影像，有效地检测和识别岩石障碍物，相对精确地确定安全着陆区，以满足制导控制的需要.本文首先提出了基于CCD相机的遥感影像的生成；然后用障碍识别算法对着陆表面的岩石、弹坑等障碍进行提取，主要是通过基于Canny理论的自适应边缘检测方法提取岩石轮廓，并利用椭圆拟合方法确定岩石的大小位置，简单有效；最后提出算法确定最大无障碍着陆圆，对着陆点的选择提供参考.在确定岩石轮廓的基础上，搜索最大的安全着陆区.初步确定了着陆区域后，在探测器精确避障阶段.首先利用Delaunay三角剖分算法以及分块插值策略，对初选着陆区域进行高精度三维建模，形成精细的星表DEM模型；然后针对高精度的模型，根据探测器当前的星下点位置，进行分块螺旋搜索，根据设定的探测器安全着陆阈值来选定安全的目标着陆点.螺旋搜索获得的目标着陆点是距离探测器当前位置最近的符合着陆阈值的着陆点.其技术方案如图1所示.

2.1 遥感影像生成

探测器在软着陆过程中，在下落到10～2 000 m之间时，其携带的CCD相机会对星球表面进行拍照.本节根据探测器所在位置，确定所携带的CCD相机拍摄的区域，基于光线跟踪算法结合高程数据对区域内的遥感影像图片进行仿真模拟.光线跟踪主要是判断光线与球面的相交情况，如图2所示：

Fig.2 The illustration of ray tracing图2 光线跟踪示意图

根据代数学的方法，设置光线方程：

P=P0+tV，

(1)

其中，t为待定参数，求解P的方程：

|P-O|2-r2=0.

(2)

消去O，可以得到：

|P0+tV-O|2-r2=0.

(3)

根据式(3)得到关于t的一元二次方程：

at2+bt+c=0.

(4)

求解方程组：

(5)

得到t，即可求得点P.

卫星相机拍照的效果主要受光照、气候以及相机特性影响.这里只考虑光照和相机的影响.通过计算光照模型以及相机拍照范围的计算，算法流程如图3所示：

Fig.3 Flow chart of remote sensing image generation图3 遥感影像生成流程图

2.2 障碍检测

本文提出的是基于边缘检测用椭圆拟合探测器目标表面的岩石及陨石坑等障碍进行检测.边缘是图像最基本的特征.所谓边缘是指图像像素灰度有阶跃变化或屋顶状变化的像素的集合.本文采用的是基于Canny算子的自适应边缘检测方法[17].这种方法将整幅图像分割为若干子图像，再根据非极大值抑制后的结果自适应地设定各子图像的高、低阈值.自适应Canny算法即二维高斯函数：

(6)

对图像进行平滑，计算像素的梯度幅值和方向：

(7)

其中，Px(i,j)为x方向偏导数，Py(i,j)为y方向偏导数，P135°(i,j)为135°方向偏导数，P45°(i,j)为45°方向偏导数.

(8)

通过图像的直方图进行高低阈值的确定.获得图像的边缘信息之后，然后利用椭圆拟合的方法来估测陨石坑的大小和位置.

岩石检测方法则是基于阴影检测，该检测技术主要分为2个过程，即阴影检测和岩石区域确定.阴影检测过程就是图像分割过程，图像阈值分割是一种广泛使用的图像分割技术，它利用了图像中要提取的目标物与其背景在灰度特性上的差异，把图像视为具有不同灰度级的2类区域的组合，选取一个合适的阈值，以确定图像中每一个像素点应该属于目标还是背景区域，从而产生相应的位置图像.这里我们用椭圆拟合的方式利用岩石阴影和太阳光角度来估计岩石区域的范围.

岩石的大体形状可以通过岩石阴影估计出来，岩石和阴影模型如图4所示.最合适的拟合椭圆能够有效的表示岩石的形状和大小[18].

Fig.4 The rock and shadow model图4 岩石和阴影模型

椭圆拟合技术是指找到一个参数集合使得给定的数据点到椭圆的距离最短，求解过程:

FK(X)=KX=ax2+bxy+cy2+dx+ey+f=0,

(9)

其中

(10)

这里定义FK(X)为点(x,y)到曲线FK(X)=0的代数距离.可以计算多个点到曲线代数距离平方和的最小值来拟合接近一般的圆锥曲线：

(11)

为了避免零解或者不同K值表示相同的圆锥曲线，对参数向量K应该有一定的限制.为了在拟合椭圆的同时保证解决线性最小平方问题的效率，这里对向量K有一定的限制从而使得一般二次曲线表示为一个椭圆.最适合的限制众所周知，即b2-4ac≥0.尽管在一般情况下，这种不等式的约束条件实行是比较困难的，在这种情况下，我们可以自由任意地缩放参数从而简单地引入一定的等式约束4ac-b2=1.这个二次的等式限制也可以表示为矩阵形式KTCK=1，即：

(12)

根据式(9)的定义，我们构造矩阵D：

(13)

可以重写为

SK=λCK,

(14)

KTCK=1，

(15)

(16)

2.3 着陆区的确定

将整个图像区域分为2部分：障碍区域和非障碍区域，目标是要求在非障碍区域寻找最大圆.关键一点是要确定最大区域的质心.由于整个图像仅分成2部分，可以用黑白颜色值来表示，白色区域表示非障碍区域.而图像形态学经常在通过阈值化获得的二进制图像中完成，因而本节将介绍利用形态学操作在白色区域中确定最大区域的质心.

基于数学形态学的图像处理是用具有一定形态的结构元素去度量和提取图像中的对应形状以达到对图像分析和识别的目的[20].在形态学中，结构元素的形态决定了运算所提取的图像中形态信息.这里利用数学形态学的膨胀和腐蚀基本变化操作，分割出其中的白色区域，并确定白色区域中连续的最大区域.假定结构元素为g，图像为f，则灰度形态膨胀[21]：

(17)

灰度形态腐蚀:

(18)

膨胀是指将一些图像与核进行卷积.核可以是任何的形状或大小，核可以视为模版或掩码，膨胀是求局部最大值的操作.腐蚀是膨胀的反操作.通过膨胀黑色区域或是腐蚀白色区域，不断逼近缩小白色区域从而确定最大的白色区域；然后利用椭圆拟合来确定不规则区域的质心作为最大圆的中心；最后从该中心扩散，寻找能够仅仅覆盖白色区域的最大圆.

算法1.基于数学形态学的着陆区确定算法.

① 设置搜索半径radus的初值，一般设为图片大小较短边值的一半；

② 以radus为搜索半径，找到质心点center(x,y)中心，在区域(x±radus,y±radus)内开始搜索；

③ 若碰到障碍区域的点Pointi(xi,yi)，设置radus=distance(pointi,center)-1，转步骤②，否则继续搜索；

④ 设radus_found=radus，radus_found≤0表示算法搜索失败(但由于是在最大连续非障碍区域中寻找，不会出现radus_found≤0)；

⑤ 确认以center(x,y)为中心、radus_found为半径的圆形区域均为非障碍区.

在确定初始的着陆区后，采用图像金字塔继续对选中的区域进行处理.图像金字塔是以多分辨率来解释图像的一种有效的结构.一幅图像的金字塔是一系列金字塔形状排列的分辨率逐步降低的图像集合.

3 三维重构算法

在着陆区域基本确定之后的精避障阶段，探测器在精避障初始时刻，开启激光三维成像仪，对初选目标着陆区域进行激光三维成像，测量目标范围内散乱点的高程数据，然后建立初选目标着陆区域的TIN模型，再根据建立的TIN模型，利用插值算法，计算出规则网格模型的坐标及高程值，形成初选着陆区域的DEM规则网格模型，最后通过三维引擎对规则网格模型进行渲染，绘制出初选着陆区域的三维地形网格.其主要流程如图5所示：

Fig.5 Three-dimensional reconstruction algorithm flow图5 三维重构算法流程

3.1 三维TIN模型构建

算法2.TIN模型构建算法(如图6所示).

Fig.6 The TIN model building algorithm图6 TIN模型构建算法

① 初始化三角网格，寻找散点集最外围的4个点构建矩形外包围框，连接任意一条对角线，对划分出的2个三角形及其边编号.

② 判断是否有未处理的散乱点，如果所有散乱点均处理完毕，转到步骤⑥.

③ 插入一个新的散乱点，得到所有外接圆包含该点的三角形，记录三角形的编号.

④ 遍历步骤③ 中记录的三角形，删除公共边，形成多边形空腔.

⑤ 连接点P与多边形空腔的各个顶点，形成新的TIN三角网，转到步骤②.

⑥ 遍历所有边，删除包含矩形外包围框中点的边，获得最终的三维TIN模型.

3.2 DEM规则网格构建

形成不规则网格TIN模型后，需要根据实际仿真的需求，确立所需DEM规则网格模型的精度，从而确定待插值DEM规则网格点的坐标.然后，开始进行DEM规则网格的构建.主要分2个步骤：

1)计算待插点在不规则三角网格模型中的三角形编号；

2)根据记录编号的三角形各个顶点的高程值，利用插值算法计算出待插点的高程值，形成DEM规则网格模型.

3.2.1 判断点在TIN模型三角形中算法

通过插值算法计算DEM规则网格点的高程值，首先需要判断规则网格点在TIN模型的哪个三角形中，针对N个散乱点所构建的TIN模型一共有2N-6个三角形，且DEM规则网格模型的精度要求较高，若每个待插值点都扫描一遍TIN模型中的所有三角形，那么重复冗余的计算量是相当大的.

为了减少搜索的次数，降低程序的冗余度，采用TIN模型分块搜索的策略，即将TIN模型按照经纬度分为大小相等的正方形块，然后只扫描遍历一次各个三角形的顶点，将每块中包含的三角形的编号记录下来(只要三角形有一个顶点在该块中则认为该块包含该三角形)，之后判断规则网格点时，先根据经纬度确定规则网格点在TIN模型的哪个分块中，然后将点P所在的块以及其邻接块所包含的三角形作为搜索的范围.针对点所在的块的位置：内部块，需要搜索9个分块；边界块，需要搜索6个分块；顶角块，需要搜索4个分块.由此，可以大大减少TIN模型中三角形搜索量，提高DEM规则网格模型的构建效率，算法示意图如图7所示：

Fig.7 The judgment point algorithms in TIN model triangle图7 判断点在TIN模型三角形中算法

3.2.2 DEM内插建网算法

算法3.DEM内插建网算法(如图8所示).

① 当AB与AC不平行于x轴时，过待插值规则网格点P作1条平行于x轴的直线与AB和AC(或其延长线)分别相交于点L和点R.当AB和AC其中有1条边平行于x轴时，则过点P作平行于x轴的直线与另外2条边分别相交于点L和点R；

② 计算点L和点R的坐标，然后计算点L和点R的高程值:

ZL=ZA+|ZB-ZA|×|XL-XA|/|XB-XA|,

(19)

ZR=ZA+|ZC-ZA|×|XR-XA|/|XC-XA|.

(20)

③ 依照点P的坐标值，计算出点P的高程值:

ZP=ZL+|ZR-ZL|×|XP-XL|/|XR-XL|.

(21)

4 精避障扫描策略

由于初避障过后，探测器是位于初选着陆区域中心上方的，为了使精避障目标着陆区域搜索的时间效率最高，我们采用螺旋搜索的策略，即从探测器的星下点开始搜索计算，计算的最小单位是探测器所要求的控制精度的范围作为边长的一个小正方形分块，然后逐层向外螺旋式扩展搜索，一旦搜索到符合需求的着陆点分块，就停止搜索.这样搜索到的目标着陆点离探测器的距离最近，且搜索时间最短，可以较优的实现探测器软着陆精避障的目的.策略如图9所示.

针对螺旋搜索算法流程中判断分块是否符合软着陆需求的问题，我们需要根据探测器的安全着陆参数来采取相应的安全判定规则.一般情况下，主要通过最大最小高度差h1和h2、坡度δ和粗糙度α来衡量当前区域是否符合软着陆的需求：

(22)

在判断式(22)阈值的过程中，按照最大最小高度差、坡度和粗糙度的顺序依次判定，当某一项不符合探测器安全着陆参数的要求时，立刻停止搜索当前块，不再计算后面的参数.这样可以缩短判定安全区域的时间，一旦当前块不符合着陆要求，则立即计算下一块的参数，从而尽可能地优化精避障安全着陆区域判定的性能.

4.1 最大最小高度差

最大高度差和最小高度差为当前分块的DEM规则网格点中的最大高程值和最小高程值与当前分块的DEM规则网格点的平均高程值的差.根据上述定义，首先求得当前分块的DEM规则网格点的平均高程值.

(23)

其中，havr表示当前分块的平均高程值，n表示当前分块DEM规则网格点的行列号，hi,j表示当前分块DEM规则网格模型中第i行、第j列的点的高程值.

然后，搜索当前分块的DEM规则网格模型中所有点的高程值，查找出max(hi,j)和min(hi,j)，将其与havr做差值，求得最大最小高度差：

h1=max(hi,j)-havr，
h2=min(hi,j)-havr，
|h1|≤hmax，
|h2|≤hmax,

(24)

其中，max(hi,j)和min(hi,j)是当前分块DEM规则网格模型中点的高程值中的最大值和最小值，hmax是探测器安全降落参数所规定的最高或最低障碍物阈值.

4.2 坡度值

坡度表示地表地形的倾斜程度，如图10所示：

Fig.10 Slope value图10 坡度值

一般采用度数法来表示相应地形分块的坡度值，即月表当前分块与水平面的夹角，等同于月表分块的法向量与z轴的夹角.

首先计算其法矢量ni,j:

(25)

然后计算其坡度值，即法向量ni,j与z轴的夹角φ.

(26)

4.3 粗糙度

粗糙度是反映地表凹凸程度的指标，本文采用平面夹角法来计算分块的粗糙度，即将分块的4个顶点对角线相连，确定4个平面，计算斜边相邻的2个平面的夹角，作为粗糙度参数.

设2个平面的夹角为θ，且θ∈[0,π]，则定义粗糙度R=1-cosθ∈[0,2].当2个平面的夹角θ=0°时，粗糙度R=0；当2个平面的θ=180°时，粗糙度R=2.这里计算2组平面的夹角粗糙度，取2组结果中较大的一个作为最终结果，提高了单一方法计算粗糙度的精度.具体计算流程如下.

1)计算构成4个平面的向量

例如平面ABC的构成向量为

IAB=(Xi,j+1-Xi,j,Yi,j+1-Yi,j,Zi,j+1-Zi,j)，
IBC=(Xi+1,j+1-Xi,j+1,Yi+1,j+1-Yi,j+1,
Zi+1,j+1-Zi,j+1).

(27)

如式(2)，分别计算出平面ABC，ADC，ABD，CBD的构成矢量.

2)计算相应平面的法向量

平面ABC法矢量n1和平面ADC的法矢量n2为

(28)

同理，计算出平面ABD和平面CBD的法矢量n3和n4.

3)计算法矢量的夹角余弦值

例如平面ABC的法矢量n1和平面ADC的法矢量n2的夹角的余弦值为

(29)

同理，计算出平面ABD和平面CBD的法矢量夹角的余弦值cosθ2.

4)计算2组地表粗糙度值

R1=1-cosθ1,
R2=1-cosθ2.

(30)

5)确定最终星球表面分块粗糙度

根据2组粗糙度取较大值的原则，确定最终的月表分块粗糙度.若R1

5 本文算法实验结果

选取月面不同区域的地形测量数据作为输入参数，对本文算法进行了实验，验证了本算法的正确性，结果如表1所示，在每组图中从左至右依次为：遥感影像生成、自适应边缘检测图、椭圆拟合障碍检测图、最大安全着陆区确定图.

Table 1 Algorithm Verification Results表1 算法验证结果

图11为DEM规则网格三维重构数据，左侧为散乱点经纬度及高程值，右侧为算法执行后构建的DEM规则网格点的经纬度及高程值.

Fig.11 DEM 3D reconstructed data图11 DEM三维重构数据

图12为利用三维引擎所绘制出的三维重构后的高精度初选着陆区域DEM规则三角网格模型.

Fig.12 DEM 3D mesh model图12 DEM三维网格模型

图13为螺旋搜索安全判定策略所计算的各项参数，依次为当前分块中心点经纬度、最大高度差、最小高度差、粗糙度和坡度值.

为了度量算法的准确性及算法效率，选取了月面不同区域的地形测量数据作为输入参数，对算法进行了实验，并与基于立体视觉和阴影分析的避障算法以及基于多因子着陆区分析评估方法进行了实验分析对比.实验结果如表2所示，结果表明本文算法在着陆区选择准确性上高于基于多因子着陆区分析评估方法的精度，达到了基于立体视觉和阴影分析算法的实验精度.在计算效率方面均高于基于多因子着陆区分析评估方法和基于立体视觉和阴影分析算法.

Fig.13 Spiral search security decision policy parameters图13 螺旋搜索安全判定策略参数

Table 2 Efficiency and Accuracy Comparison of Algorithms表2 算法效率及准确性对比

6 结 论

探测器软着陆障碍检测与安全着陆是未来的深空探测的关键技术之一，本文提出了一种两阶段的探测器自主安全着陆仿真算法.首先给出了一种基于遥感影像生成的着陆区岩石障碍的检测以及安全着陆区的确定方法;然后根据探测器软着陆精避障的基本要求，提出了一种基于三维重构的探测器精避障螺旋扫描方法.对初选着陆区域进行高精度月面三维重构，螺旋搜索算法提高了目标着陆点搜索的时间效率25%以上，且寻找到最近的符合安全着陆阈值的目标着陆点，节约探测器能源.

本文的算法方案通过确定CCD相机拍摄的区域，根据高程数据对区域内的遥感影像图片进行模拟，利用自适应的边缘检测方法和阴影检测以及椭圆拟合技术对着陆区的障碍进行了检测.提出了基于数学形态学的最大着陆圆算法，该算法快速有效，可用于实时障碍检测，同时有效确定相对安全着陆区，为探测器的精确避障提供了良好的依据.该算法在探测器下降过程中，能够对多帧图像进行障碍检测，从而保证着陆区选择的有效性.基于三维重构的探测器精避障螺旋扫描方法.对初选着陆区域进行高精度月面三维重构提高了软着陆的避障精度.本文算法准确性达到93%以上.经对比分析，本文算法在着陆区选择准确性上高于基于多因子着陆区分析评估方法的精度，达到了基于立体视觉和阴影分析算法的实验精度.在计算效率方面均高于基于多因子着陆区分析评估方法和基于立体视觉和阴影分析算法.