CN106228131A - 行星着陆器自适应障碍检测方法 - Google Patents

Abstract

本发明提供一种行星着陆器自适应障碍检测方法，涉及图像处理领域。为了解决现有技术中当障碍地形数量大于所有地形数量的50％时无法进行平面拟合的问题，本发明通过获取高程数据，根据其中的坐标点计算初始平面，计算其他点到平面的残差值，并将残差值输入至混合高斯模型进行聚类分析，两种残差组，并将残差组进行尺度因子算法计算排除障碍点，得到非障碍点集合，根据非障碍点集合拟合平面，并将高程数据中的点于拟合平面做差，再与预设的安全阈值比较，从而判断哪些是最终的非障碍点，哪些是最终的障碍点。本发明适用于火星着陆器。

Description

行星着陆器自适应障碍检测方法

技术领域

本发明涉及自适应障碍检测方法，属于图像处理领域。

背景技术

行星着陆器动力下降段障碍检测的主要目的是确保探测器能够安全地着陆到相对平坦的区域，以火星着陆为例，由于火星表面布满陨石坑，岩石，陡坡等地形，因此如何检测出不适合着陆的障碍点并排除这些障碍点是最终着陆过程的关键所在。

现有的方法是通过激光雷达获取火星地貌高程信息，根据探测器自身尺寸以及着陆误差的总尺寸，将获取到的高程地图划分成各个小块地形，对每一块小的地形进行平面匹配，再将得到的平面与原始的高程信息相减就会得到每一小块地形的粗糙度，同时通过计算该平面的法向量与火星地形竖直方向的法向量的夹角就可以得到每一小块地形的坡度大小。探测器安全着陆要求障碍物不得高于某一数值(粗糙度阈值)坡度也不能超过某一个角度(坡度阈值)，这样通过寻找满足阈值条件的地形为着陆区，实现探测器的安全着陆。

因此，如何较为精确地对每一小块地形(局部地形)进行平面拟合是关键，传统的方法是采用最小平方中位数法来进行平面的拟合，但是这种方法不具有鲁棒性，当障碍物的丰富度大于50％的时候就会出现崩溃现象，即无法准确拟合，这种崩溃对于探测器着陆来说就有可能是致命的，然而，就一般情况而言，障碍物越丰富的地区，所含有的科学信息越丰富，这也是未来火星探测的主要降落方向，因此这又是必须克服的问题。同时，虽然已经有很多其他的回归分析方法在平面拟合上声称具有高鲁棒性，但是这些方法大多数需要根据所需要探测的地形人为的设置很多的参数，这些对于未知地形来说不契合实际。因此如何在障碍物丰富的区域进行平面拟合，以及保证拟合的鲁棒性成为着陆器障碍检测领域的重要问题。

发明内容

本发明提供一种行星着陆器自适应障碍检测方法，为了解决现有技术中当障碍地形数量大于所有地形数量的50％时无法进行平面拟合的问题。

行星着陆器自适应障碍检测方法包括：

步骤一：获取高程数据，所述高程数据是在着陆点水平坐标系下的预选着陆区范围内的三维坐标集合，所述三维坐标集合包括坐标点；

步骤二：根据随机选择的三个坐标点确定初始平面，并计算高程数据中的其他坐标点到所述初始平面的残差值，所述残差值与所述坐标点存在对应关系；

步骤三：将所述其他坐标点的残差值输入至混合高斯模型中，得到第一残差组和第二残差组，并且所述第一残差组的标准差值小于所述第二残差组的标准差值，所述第一残差组以及第二残差组均由残差值组成；

步骤四：将所述第一残差组中的残差值通过中位尺度因子算法进行迭代计算，得到候选残差值集合，所述候选残差值集合对应的坐标组成候选坐标集合；

步骤五：重复步骤二至步骤四，得到优选的候选坐标集合，重复的次数由预想的着陆区障碍物分布程度确定；

步骤六：将所述优选的候选坐标集合通过最小残差平方和准则计算得到最优坐标集合；

步骤七：将所述最优坐标组合通过平面拟合算法得到拟合平面；

步骤八：将所述高程数据中的每个坐标点与所述拟合平面做差，分别得到与所述每个坐标点对应的相对高度数据以及坡度数据，将所述每个坐标点对应的相对高度数据以及坡度数据与安全阈值进行比较，若小于或等于所述安全阈值，则将当前点判定为非障碍点，若大于所述安全阈值，则将当前点判定为障碍点。

本发明的技术效果是，在进行平面拟合之前先进行聚类分析，排除掉标准差较大的坐标集合，再对由此得到的剩余坐标集合进行平面拟合，这样在实际进行平面拟合时，已经排除掉了较多的障碍点，可以解决由于障碍点过多导致无法拟合的问题，从而提高了障碍检测系统的鲁棒性。同时本发明不需要设置过多的参数，所有步骤中仅需要一个对障碍物丰富度的粗略估计值即可达到障碍检测的目的，并且具有很强的自适应性。

附图说明

图1为本发明的行星着陆器自适应障碍检测方法的流程图；

图2为使用matlab模拟的火星局部地形图；

图3为蒙特卡罗模拟结果；

图4为火星局部真实地形图；

图5为本发明的方法与现有技术方法的效果对比图。

具体实施方式

具体实施方式一：本实施方式的行星着陆器自适应障碍检测方法，包括以下步骤：

步骤一：获取高程数据，所述高程数据是在着陆点水平坐标系下的预选着陆区范围内的三维坐标集合，所述三维坐标集合包括坐标点；

步骤二：根据随机选择的三个坐标点确定初始平面，并计算高程数据中的其他坐标点到所述初始平面的残差值，所述残差值与所述坐标点存在对应关系；

步骤三：将所述其他坐标点的残差值输入至混合高斯模型中，得到第一残差组和第二残差组，并且所述第一残差组的标准差值小于所述第二残差组的标准差值，所述第一残差组以及第二残差组均由残差值组成；

步骤四：将所述第一残差组中的残差值通过中位尺度因子算法进行迭代计算，得到候选残差值集合，所述候选残差值集合对应的坐标组成候选坐标集合；

步骤五：重复步骤二至步骤四，得到优选的候选坐标集合，重复的次数由预想的着陆区障碍物分布程度确定；

步骤六：将所述优选的候选坐标集合通过最小残差平方和准则计算得到最优坐标集合；

步骤七：将所述最优坐标组合通过平面拟合算法得到拟合平面；

步骤八：将所述高程数据中的每个坐标点与所述拟合平面做差，分别得到与所述每个坐标点对应的相对高度数据以及坡度数据，将所述每个坐标点对应的相对高度数据以及坡度数据与安全阈值进行比较，若小于或等于所述安全阈值，则将当前点判定为非障碍点，若大于所述安全阈值，则将当前点判定为障碍点。

具体而言，在步骤二中，随机从高程数据集X＝{x_i},i＝1,2.......m选三点生成平面，m为高程数据中坐标点的个数，高程数据中的坐标点通过激光雷达获取，获取到的坐标点个数m为有限值。平面方程为：

n﹒x+d＝0

式中，平面法向量n由下式求出:

n＝(x₁-x₂)×(x₃-x₂)

平面截距d由下式求出:

d＝-nx₂

x₁，x₂，x₃为随机选出的三点。

每一点的残差值r_i定义为剩余点带入该平面方程中所得值的平方，即

{r_i}＝(n﹒x_i+d)²

在步骤三中，将每一点的残差值输入至混合高斯模型，最终得到两类残差组，其中标准差较小的为第一残差组，我们认为第一残差组中的可能拥有较多的安全着陆点，因此保留这部分用于下一步骤的计算。

在步骤四中，进一步剔除障碍点，将满足如下公式的残差值对应的坐标定义为安全着陆点的候选点。

$\left\{r_i\right\}=\left\{r_i\right|r_i<2.5{\widehat{\mathrm{\&sigma;}}}_{1,i-1},i=1,\mathrm{2...}m\}$

上式中，为上一步迭代的残差尺度阈值，其具体形式在具体实施方式五中给出。

在步骤五中，重复次数的确定方法为：

随机从数据点中取三点，共取N次，每一次取点，这三点正好是非障碍点的概率P为：

P＝1-(1-(1-ε)³)^N

其中，ε为预想的着陆区障碍物的分布程度。因此可以求出需要重复的次数N

N＝log(1-P)/log(1-(1-ε)³)

ε的取值越高，表示障碍点占的比重越多，也就表示我们对障碍点的估计越保守，ε为0到1之间的任意值。

在步骤六中，通过最小残差平方和准则计算得到最优坐标集合的具体过程为：统计候选坐标集合的残差值的平方和；对于外点，则为最终尺度因子的平方与外点个数的乘积。

$\mathrm{\&Psi;}=\underset{i=1}{\overset{n_{in}}{\&Sigma;}}{r_i}^2+\underset{i=1}{\overset{m-n_{in}}{\&Sigma;}}2.5{{\widehat{\mathrm{\&sigma;}}}_1}^2$

其中，n_in表示第i次迭代时的内点个数。本发明中所表述的内点均指非障碍点，外点均指障碍点。最小Ψ函数所对应的那一组内点数据就是最终真实的内点数据，利用平面拟合算法进行平面拟合，然后进行步骤七，确定障碍物，完成行星着陆器在下降段的障碍检测。

可以将障碍点坐标赋予属性值0，将安全着陆点坐标赋予属性值1，于是得到一个能够表示行星地表障碍分布的矩阵，后续可以根据安全着陆点最大连通区域以及着陆器自身的燃料消耗等情况具体决定应该着陆在哪片区域。

具体实施方式二：本实施方式与具体实施方式一不同的是：高程数据根据行星着陆器导航信息以及激光雷达参数确定；预想着陆点为在行星着陆器发射前确定的着陆椭圆的中心。

本发明的自适应障碍检测算法是行星着陆器在动力下降段末期开始执行的，在动力下降段末期激光雷达开始工作，其原理是雷达波的反射。通过激光雷达获取高程数据还需要考虑到姿态信息和着陆器距离预想着陆点的位置。其中姿态角包括俯仰、滚转和偏航角。预想着陆点是在着陆器发射之前，工作人员预先设定的着陆椭圆的中心点。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一至二之一不同的是：在步骤三中将残差值输入至混合高斯模型并得到第一残差组和第二残差组的具体过程为：

根据所述混合高斯模型对所述残差值的概率密度分布p(x_i)进行估计，所述混合高斯模型的子高斯模型N(x_i；μ_j,Σ_j)的个数M为2，表示将所述残差值分为两类，一类为障碍点，另一类为非障碍点，所述残差值的概率密度分布p(x_i)为子高斯模型的加权之和，即：

$p\left(x_i\right)=\underset{j=1}{\overset{M}{\&Sigma;}}{\mathrm{\&omega;}}_{j}N(x_i;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)$

$N(x_i;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)=\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^m\left|{\mathrm{\&Sigma;}}_j\right|}}\exp \&lsqb;-{(x-{\mathrm{\&mu;}}_j)}^T{\mathrm{\&Sigma;}}_j^{-1}(x-{\mathrm{\&mu;}}_j)\&rsqb;$

上式中ω_j为第j个子高斯模型的权重系数，μ_j为其对应的均值，Σ_j为其对应的方差，ω_j、μ_j以及Σ_j的求解方法为最大期望值算法。

将残差的概率密度分布G写成含有两个子高斯分布F、H的混合分布模型

G＝(1-ω)·F+ω·H

ω为混合高斯分布的权重系数，对于第i个残差数据，它由第k个model生成的概率为：

${\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_i^k=\frac{{\mathrm{\&omega;}}_{k}N\left(x_i\right|{\mathrm{\&mu;}}_k,{\mathrm{\&sigma;}}_k)}{{\mathrm{\&Sigma;}}_{j=1}^2{\mathrm{\&omega;}}_{j}N\left(x_i\right|{\mathrm{\&mu;}}_j,{\mathrm{\&sigma;}}_j)}$

首先假设高斯混合模型参数已知，由初始值或上一步迭代确定，此步称之为E步骤；接下来，采用最大似然估计计算各个自模型的均值与方差，并更新权重系数：

${\mathrm{\&mu;}}_k=\frac{1}{m}\underset{i=1}{\overset{m}{\&Sigma;}}{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_i^k{r}_i$

${\mathrm{\&sigma;}}_k=\frac{1}{m_k}\underset{i=1}{\overset{m}{\&Sigma;}}{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_i^k(r_i-{\mathrm{\&mu;}}_k){(r_i-{\mathrm{\&mu;}}_k)}^T$

$m_k=\underset{i=1}{\overset{m}{\&Sigma;}}{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}_i^k$

此步骤称之为M步骤。重复上述E、M两步骤直到算法收敛。

其它步骤及参数与具体实施方式一至二之一相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：在步骤四中通过中位尺度因子算法进行迭代计算，得到候选残差值集合具体过程为：

根据如下公式进行计算：

${\widehat{\mathrm{\&sigma;}}}_{1,i}=\frac{\left|\right|r_K|-\mathrm{\&mu;}|}{{\mathrm{\&Phi;}}^{-1}\left(0.5\right(1+\frac{n_K}{n_i}))\&CenterDot;({(1-\mathrm{\&omega;})}^2+{\mathrm{\&omega;}}^{2}T)}$

其中，Φ^-1为逆高斯密度累计函数，|r_K|为步骤二中得到的所有残差值按照从小到大顺序排列后的中位值，μ为所述步骤二中确定的所有残差值按照从小到大顺序排列后的均值，T为所述第二残差组的标准差与所述第一残差组的标准差的比值，n_i为上一步迭代内点的个数，n_i的初始值为第一残差组中残差的个数，n_K为残差值满足如下条件的数据点个数：

$\left|{\tilde{r}}_i\right|<2.5{\widehat{\mathrm{\&sigma;}}}_{1,i-1}$

为第i次迭代时得到的残差值集合，通过迭代n_i值，直到尺度因子收敛到固定值再通过公式确定本次循环的最终非障碍点的残差值集合。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤一在行星着陆器距离火星地表100米至500米时开始执行。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤五中的预想的着陆区障碍物分布程度通过先前火星轨道器拍摄的火星地表照片进行人为估算得出，属于人为给定的参数，但这个值对本发明的障碍检测效果影响不大，一般来说这个值越高代表着陆区的障碍物数量越多，也就需要更加的保守估计。因此，在实际应用中如果需要保守估计，我们可以将这个值预设为较大的值，例如70％。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：在步骤七中，所述平面拟合算法为最小二乘法。其中，使用其他的平面拟合算法也是允许的，由于本发明已经剔除了较多的障碍点，对于剩余的点，使用较为简单的平面拟合算法即可，因此本方案选择了最小二乘法进行平面拟合。

其它步骤及参数与具体实施方式一至六之一相同。

本发明的障碍检测方法可以通过matlab软件进行仿真实验。图2是用matlab软件模拟出的火星表面局部地形图。其中，图2中的基准平面已知，以拟合平面与基准平面的夹角作为误差量，进行了500次的蒙特卡罗实验，比较了自适应障碍检测算法与传统的最小平方中位数法(LMedsq)的误差大小，如图3所示。图3中横轴为障碍物的丰富度，纵轴为拟合误差，方形节点对应的曲线为本发明的自适应障碍检测方法，圆形节点对应的曲线为现有技术的最小平方中位数法。为了更好地展示自适应障碍检测算法的效果，又以图4中的真实火星局部地形为样本，利用激光雷达获取的三维高程图进行平面拟合，同样比较了上述两种方法的效果，如图5所示。图5中处于上方的平面为本发明的自适应障碍检测方法拟合出的平面，下方的平面为现有技术的最小平方中位数法的平面，可看出现有技术的拟合方法出现了明显的偏差，特别是被椭圆圈出的区域，偏差尤其严重。通过上述的数值模拟和真实数据拟合，可以看到，自适应障碍检测方法相比现有技术的方法具有更好的表现。

Claims (8)

1.一种行星着陆器自适应障碍检测方法，其特征在于，包括以下步骤：

步骤一：获取高程数据，所述高程数据是在着陆点水平坐标系下的预选着陆区范围内的三维坐标集合，所述三维坐标集合包括坐标点；

步骤二：根据随机选择的三个坐标点确定初始平面，并计算高程数据中的其他坐标点到所述初始平面的残差值，所述残差值与所述坐标点存在对应关系；

步骤三：将所述其他坐标点的残差值输入至混合高斯模型中，得到第一残差组和第二残差组，并且所述第一残差组的标准差值小于所述第二残差组的标准差值，所述第一残差组以及第二残差组均由残差值组成；

步骤四：将所述第一残差组中的残差值通过中位尺度因子算法进行迭代计算，得到候选残差值集合，所述候选残差值集合对应的坐标组成候选坐标集合；

步骤五：重复步骤二至步骤四，得到优选的候选坐标集合，重复的次数由预想的着陆区障碍物分布程度确定；

步骤六：将所述优选的候选坐标集合通过最小残差平方和准则计算得到最优坐标集合；

步骤七：将所述最优坐标组合通过平面拟合算法得到拟合平面；

步骤八：将所述高程数据中的每个坐标点与所述拟合平面做差，分别得到与所述每个坐标点对应的相对高度数据以及坡度数据，将所述每个坐标点对应的相对高度数据以及坡度数据与安全阈值进行比较，若小于或等于所述安全阈值，则将当前点判定为非障碍点，若大于所述安全阈值，则将当前点判定为障碍点。

2.根据权利要求1所述的行星着陆器自适应障碍检测方法，其特征在于，所述高程数据根据行星着陆器距离预想着陆点的位置、姿态角以及激光雷达参数确定；所述预想着陆点为在所述行星着陆器发射前对着陆点区域的预测值。

3.根据权利要求1所述的行星着陆器自适应障碍检测方法，其特征在于，在步骤三中将残差值输入至混合高斯模型并得到第一残差组和第二残差组的具体过程为：

根据所述混合高斯模型对所述残差值的概率密度分布p(x_i)进行估计，所述混合高斯模型的子高斯模型N(x_i；μ_j,Σ_j)的个数M为2，表示将所述残差值分为两类，一类为障碍点，另一类为非障碍点，所述残差值的概率密度分布p(x_i)为子高斯模型的加权之和，即：

$p\left(x_i\right)=\underset{j=1}{\overset{M}{\&Sigma;}}{\mathrm{\&omega;}}_{j}N(x_i;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)$

$N(x_i;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)=\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^m\left|{\mathrm{\&Sigma;}}_j\right|}}\exp \&lsqb;-{(x-{\mathrm{\&mu;}}_j)}^T{\mathrm{\&Sigma;}}_j^{-1}(x-{\mathrm{\&mu;}}_j)\&rsqb;$

上式中T为所述第二残差组的标准差与所述第一残差组的标准差的比值，m为所述高程数据中坐标点的个数，ω_j为第j个子高斯模型的权重系数，μ_j为其对应的均值，Σ_j为其对应的方差，ω_j、μ_j以及Σ_j的求解方法为最大期望值算法。

4.根据权利要求2或3所述的行星着陆器自适应障碍检测方法，其特征在于，在步骤四中通过中位尺度因子算法进行迭代计算，得到候选残差值集合具体过程为：

根据如下公式进行计算：

${\widehat{\mathrm{\&sigma;}}}_{1,i}=\frac{\left|\right|r_K|-\mathrm{\&mu;}|}{{\mathrm{\&Phi;}}^{-1}\left(0.5\left(1+\frac{n_K}{n_i}\right)\right)\&CenterDot;\left({\left(1-\mathrm{\&omega;}\right)}^2+{\mathrm{\&omega;}}^{2}T\right)}$

其中，为第i次迭代的尺度因子，Φ^-1为逆高斯密度累计函数，ω为混合高斯分布的权重系数，|r_K|和μ分别为所述步骤二中确定的所有残差值按照从小到大顺序排列后的中位值和均值，n_i为第i次迭代时非障碍点的个数，n_i的初始值为第一残差组中残差的个数，n_K为残差值满足如下条件的数据点个数：

$\left|{\tilde{r}}_i\right|<2.5{\widehat{\mathrm{\&sigma;}}}_{1,i-1}$

其中为第i次迭代时得到的残差值集合；通过迭代n_i值，直到尺度因子收敛到固定值再通过公式确定本次循环的最终的非障碍点的残差值集合。

5.根据权利要求1或4所述的行星着陆器自适应障碍检测方法，其特征在于，步骤一在行星着陆器距离火星地表100米至500米时开始执行。

6.根据权利要求5所述的行星着陆器自适应障碍检测方法，其特征在于，步骤五中的预想的着陆区障碍物分布程度通过先前火星绕飞轨道器拍摄的火星地表照片或高程数据进行人为估算得出。

7.根据权利要求1或6所述的行星着陆器自适应障碍检测方法，其特征在于，在步骤七中，所述平面拟合算法为最小二乘法。

8.根据权利要求7所述的行星着陆器自适应障碍检测方法，其特征在于，所述重复的次数确定的方法为：

随机从所述坐标点中取三点，共取N次，每一次取点，这三点正好是非障碍点的概率P为：

P＝1-(1-(1-ε)³)^N

其中，ε为预想的着陆区障碍物的分布程度，重复的次数N由下式得出：

N＝log(1-P)/log(1-(1-ε)³)。