CN106228521B - 一种基于薄板样条插值的障碍物特征提取方法 - Google Patents

Abstract

一种基于薄板样条插值的障碍物特征提取方法，本发明涉及障碍物特征提取方法。本发明是要解决现有的基于平面拟合的障碍检测算法在高海拔时不能有效地提取障碍物的问题，而提出的一种基于薄板样条插值的障碍物特征提取方法。该方法是通过首先滤掉噪声点，自下而上建立控制点金字塔。通过金字塔顶层控制点建立初始薄板样条插值曲面。求该曲面与金子塔下一层控制点的残差，通过阈值处理剔除不满足条件的控制点，更新控制点矩阵，通过剩余的控制点建立新的插值曲面，重复上述过程，最终生成原始数据点的插值曲面等步骤实现的。本发明应用于障碍物特征提取领域。

Description

一种基于薄板样条插值的障碍物特征提取方法

技术领域

本发明涉及障碍特征检测方法，特别涉及一种基于薄板样条插值的障碍物特征提取方法。

背景技术

行星着陆器动力下降段障碍检测的主要目的是确保探测器能够安全地着陆到相对平坦的区域，以火星着陆为例，由于火星表面布满陨石坑，岩石，陡坡等地形，因此如何寻找出一块满足着陆要求的地形是最终着陆过程的关键所在。现有的方法都是在动力下降段末期(100米到500米)通过平面拟合求残差来探测障碍物。但是在较低的海拔范围内，一旦发现着陆区不具备安全着陆的可能，考虑到燃料的限制，探测器将很难做出大范围的机动。而且，海拔高度较低的情况下，探测器传感器的视角已经变得很小，不能够大范围的勘测地形，寻找合适的着陆点。因此需要探测器尽早地实施障碍检测操作，在动力下降段初期完成大范围的障碍检测，选择一块合适的着陆区域。待探测器在动力下降末期，在初期选择的着陆区域内选择合适的着陆点进行安全着陆。这样能够避免直接末期选择可能面临的上述问题。

然而在高海拔条件下，激光雷达获取的地形精度有限，激光脚点之间距离较大。传统的平面拟合方法并不适用。平面拟合的方法需要选择拟合窗口，窗口过大会过平滑障碍物；过小的窗口由于激光脚点距离较大，在窗口范围内采样点个数较少，受噪声影响较大。不能够有效地提取障碍物信息。

因此，如何在高海拔低分辨率，大采样距离的情况下，有效地检测障碍物是本发明重点解决的问题。

发明内容

本发明的目的是为了解决现有的基于平面拟合的障碍检测算法在高海拔时不能有效地提取障碍物的问题，而提出的一种基于薄板样条插值的障碍物特征提取方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、根据行星着陆器在动力下降段的位置、姿态及扫描激光雷达参数获取Flash雷达视场内地形的三维高程数据；其中，Flash雷达视场内地形的三维高程数据指的是地形三维坐标点，包含两个水平方向的位置和高度信息：三维高程数据为地形的三维坐标x_i,y_i,z_i；i＝1,2,3,...N；

步骤二、对步骤一中获取的三维高程数据进行数字形态学处理，去除低值噪声点得到处理过的激光点数据；

步骤三、将经过步骤二处理过的激光点数据进行降采样处理；将三维高程数据进行隔点采样，即根据2倍的相邻采样点距离，以自下而上的方式建立数据金字塔，金字塔的层数为M，确定金字塔中任意m层的数据点为m+1层2*2邻域内的最小高程数据点；将最小高程数据点称之为控制点；m∈M；

步骤四、利用步骤三得到的每一层的控制点建立薄板样条插值曲面S；

步骤五、利用步骤四得到数据金字塔第m层薄板样条曲面S与数据金字塔第m+1层的控制点做差得到残差，将残差通过阈值处理，剔除大于阈值的控制点，更新控制点数据，其中，阈值为金字塔每一层数据的均值；

步骤六、从金字塔顶层，自上而下开始利用更新后的控制点数据重复步骤四和步骤五；最终得到数据金字塔最底层的更新后的控制点，利用最底层的更新后的控制点生成最终的薄板样条插值曲面；

步骤七：利用步骤六生成的最终的薄板样条插值曲面与原高程数据r做差得到残差数据，利用残差数据建立粗糙度图；

步骤八：求步骤六生成的最终的薄板样条插值曲面梯度，利用曲面梯度求地形坡度图。

发明效果

本发明涉及激光雷达数据滤波、数字图像形态学处理技术领域，具体涉及一种基于薄板样条插值的障碍特征检测方法。本发明要解决现有障碍物提取算法在探测器动力下降段初期无法有效的提取障碍物信息、不能区分高程变化来自于岩石还是地形坡度等问题。本发明涉及激光雷达数据滤波、数字图像形态学处理技术领域。步骤如下：激光雷达获取距离数据通过探测器姿态信息，重构行星表面三维高程数据；利用形态学处理去除低值噪声点。建立多层控制点金字塔模型，控制点为局部窗口内的高程最小值；利用最顶层的控制点建立薄板样条插值曲面；计算该层插值曲面与下一层控制点之间的残差，通过自适应阈值剔除不满足阈值条件的控制点；利用剩余的控制点建立该层的插值曲面；重复以上建立曲面，剔除控制点的过程，最终建立底层数据点(原始高程数据)的薄板样条曲面，该曲面即为所求；求原始高程数据与最后得到的曲面的残差，建立粗糙度图，求插值曲面的梯度建立坡度图。

本发明将形态学运算应用到激光点云的噪声滤波上，利用开运算去掉高值噪声点，利用闭运算去掉低值噪声点。经过降噪后的点云数据通过局部窗口最小建立起控制点金字塔，每一层生成薄板样条曲面。通过阈值处理，不断迭代到底层，最终得到原始高程数据的插值曲面。控制点金字塔结构的建立能够提高算法的鲁棒性，图4中的2500次蒙特卡洛实验就是鲁棒性的量化数据，在障碍物丰富度达到60％的时候依然能够保持0.72以上的拟合优度。插值曲面将岩石与地形坡度带来的高程变化有效地区分开来，坡度完全来自于插值曲面，而岩石的高度计算则来自于高程数据和插值曲面的残差。避免了传统的激光点云滤波技术无法确定引起高度变化的是地形本身坡度起伏还是岩石。

为了测试本实验的具体效果，利用matlab软件模拟火星表面局部地形(图2)，其中，图2中的真值曲面已知，如图3所示。以拟合优度R²判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验，结果如图4所示。可以看到，基于样条插值的障碍物提取方法具有很高的鲁棒性，当障碍物分布程度达到60％时，R²依然可以达到0.8以上。为了更好地展示算法的效果，又对真实火星地形(图5(a))建立了插值曲面，结果显示，本发明的方法可以很好的检测出岩石、坡度等障碍。图6(a)和(b)展示的是由残差值与梯度值经过步骤七、步骤八生成的障碍物图与坡度图，绝大多数的障碍物都被提取出来。

附图说明

图1为具体实施方式一提出的一种基于薄板样条插值的障碍物特征提取方法流程图；

图2为具体实施方式一提出的模拟火星地形图；

图3为具体实施方式一提出的模拟地形真值曲面图；

图4为具体实施方式一提出的2500次蒙特卡洛实验，1为平均拟合优度；

图5(a)为具体实施方式一提出的真实火星DEM示意图；

图5(b)为具体实施方式一提出的插值曲面图；

图5(c)为具体实施方式一提出的梯度图；

图5(d)为具体实施方式一提出的残差图；

图6(a)为具体实施方式一提出的粗糙度图；

图6(b)为具体实施方式一提出的坡度图；

图7为具体实施方式一提出的数据金子塔示意图；

图8为具体实施方式一提出的最终的差值曲面图。

具体实施方式

具体实施方式一：结合图1本实施方式的一种基于薄板样条插值的障碍物特征提取方法，具体是按照以下步骤制备的：

步骤一、根据行星着陆器在动力下降段(动力下降段是指探测器抛伞之后，反推火箭开始工作的阶段)的位置、姿态及扫描激光雷达参数获取Flash雷达视场内地形的三维高程数据；其中，Flash雷达视场内地形的三维高程数据指的是地形三维坐标点，包含两个水平方向的位置和高度信息：三维高程数据为地形的三维坐标x_i,y_i,z_i；x_i,y_i为说水平方向的两个坐标，z_i为高度方向上的坐标；i＝1,2,3,...N；

步骤二、对步骤一中获取的三维高程数据进行数字形态学处理，去除低值噪声点得到处理过的激光点数据；

步骤三、将经过步骤二处理过的激光点数据进行降采样处理；将三维高程数据进行隔点采样，即根据2倍的相邻采样点距离，以自下而上的方式建立数据金字塔，金字塔的层数为M＝4，确定金字塔中任意m层的数据点为m+1层2*2邻域内的最小高程数据点；将最小高程数据点称之为控制点；m∈M；如图7；

数据金子塔通过对高程数据隔点降采样获得，最底层尺寸为原始高程数据尺寸；层数由如下公式计算：

式中，length为原始高程数据的尺寸，例如原始数据尺寸为：128×128。那么length就等于128。

为了确保顶层数据能够生成大致描述地形走势的曲面，要求其尺寸不能太小，本发明选择金字塔顶层尺寸为16×16，因此需要在层数公式中减3。

步骤四、利用步骤三得到的每一层的控制点建立薄板样条插值曲面S；

步骤五、利用步骤四得到数据金字塔第m层薄板样条曲面S与数据金字塔第m+1层的控制点做差得到残差，将残差通过阈值处理，剔除大于阈值的控制点，更新控制点数据，其中，阈值为金字塔每一层数据的均值；

步骤四和步骤五所述的顶层生成薄板样条插值曲面，用该曲面与下一层数据做差，得到残差值，对残差值进行阈值处理，剔除大于阈值的点，用剩下的点再生成薄板样条插值曲面(第二层)。用第二层曲面和第三层的数据做差，再阈值处理，用剩下的点再生成薄板样条插值曲面(第三层)。用第二层曲面和最低层的数据做差，再阈值处理，用剩下的点再生成薄板样条插值曲面，即为最终的差值曲面，如图8；

步骤六、从金字塔顶层，自上而下开始利用更新后的控制点数据重复步骤四和步骤五。最终得到数据金字塔最底层的更新后的控制点，利用最底层的更新后的控制点生成最终的薄板样条插值曲面如图5(b)；

步骤七、利用步骤六生成的最终的薄板样条插值曲面与原高程数据r做差得到残差数据，利用残差数据建立粗糙度图如图5(c)；

步骤八、求步骤六生成的最终的薄板样条插值曲面梯度，利用曲面梯度求地形坡度图如图5(d)。

本实施方式效果：

本实施方式涉及激光雷达数据滤波、数字图像形态学处理技术领域，具体涉及一种基于薄板样条插值的障碍特征检测方法。本实施方式要解决现有障碍物提取算法在探测器动力下降段初期无法有效的提取障碍物信息、不能区分高程变化来自于岩石还是地形坡度等问题。本实施方式涉及激光雷达数据滤波、数字图像形态学处理技术领域。步骤如下：激光雷达获取距离数据通过探测器姿态信息，重构行星表面三维高程数据；利用形态学处理去除低值噪声点。建立多层控制点金字塔模型，控制点为局部窗口内的高程最小值；利用最顶层的控制点建立薄板样条插值曲面；计算该层插值曲面与下一层控制点之间的残差，通过自适应阈值剔除不满足阈值条件的控制点；利用剩余的控制点建立该层的插值曲面；重复以上建立曲面，剔除控制点的过程，最终建立底层数据点(原始高程数据)的薄板样条曲面，该曲面即为所求；求原始高程数据与最后得到的曲面的残差，建立粗糙度图，求插值曲面的梯度建立坡度图。

本实施方式将形态学运算应用到激光点云的噪声滤波上，利用开运算去掉高值噪声点，利用闭运算去掉低值噪声点。经过降噪后的点云数据通过局部窗口最小建立起控制点金字塔，每一层生成薄板样条曲面。通过阈值处理，不断迭代到底层，最终得到原始高程数据的插值曲面。控制点金字塔结构的建立能够提高算法的鲁棒性图4中的2500次蒙特卡洛实验就是鲁棒性的量化数据，在障碍物丰富度达到60％的时候依然能够保持0.72以上的拟合优度。插值曲面将岩石与地形坡度带来的高程变化有效地区分开来，坡度完全来自于插值曲面，而岩石的高度计算则来自于高程数据和插值曲面的残差。避免了传统的激光点云滤波技术无法确定引起高度变化的是地形本身坡度起伏还是岩石。

为了测试本实验的具体效果，利用matlab软件模拟火星表面局部地形(图2)，其中，图2中的真值曲面已知，如图3所示。以拟合优度R²判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验，结果如图4所示。可以看到，基于样条插值的障碍物提取方法具有很高的鲁棒性，当障碍物分布程度达到60％时，R²依然可以达到0.8以上。为了更好地展示算法的效果，又对真实火星地形(图5(a))建立了插值曲面，结果显示，本实施方式的方法可以很好的检测出岩石、坡度等障碍。图6(a)和(b)展示的是由残差值与梯度值经过步骤七、步骤八生成的障碍物图与坡度图，绝大多数的障碍物都被提取出来。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中对步骤一中获取的三维高程数据进行数字形态学处理，去除低值噪声点得到处理过的激光点数据的具体操作步骤如下：

步骤二一、对原始高程数据r进行开运算，开运算是对数据进行腐蚀运算后在进行膨胀预算；窗口尺寸w_o为探测器最大尺寸的2倍；

r_o＝(r⊙w_o)⊕w_o (1)

其中，r_o为经过开运算处理后的高程数据；⊙代表腐蚀运算，具体定义为：

其中，r_er为经过腐蚀处理后的高程数据，z_i为高程数据点在窗口w_o下的高度值；⊕代表膨胀运算，具体定义为：

步骤二二、对r_o再进行闭运算，闭运算是对数据先进行膨胀运算再进行腐蚀运算，窗口尺寸w_c为探测器最大尺寸：

r_c＝(r_o⊕w_c)⊙w_c (4)

r_c为经过闭运算滤波后的点云数据；

激光雷达点高程据进行形态学开运算处理过后，将会剔除掉高值噪声点；

对经过开运算处理过的高程数据进行闭运算处理后，将会剔除点高程数据中的低值噪声点；

经过上述的开运算、闭运算处理后的数据点可以认为是已经不含有噪声点的数据了。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤四中所述的建立薄板样条曲面的具体操作步骤如下：

步骤四一、给定控制点，薄板样条曲面为通过所有的控制点的弯曲最小的光滑曲面；光滑曲面弯曲程度由下述能量函数E_s定义为：

式中，λ为正则参数，控制曲面的光滑度，S为光滑的薄板样条插值曲面；S＝{S(x_i,y_i)|i＝1,2,3,..,N}；

步骤四二、通过最小化E_s后，给定滤波后的高程数据r_c、正则参数λ。ω和a由下列线性方程组求得：

式中，K，C和O为子矩阵，O为3×3零矩阵，定义如下：

K＝{K_i,j}＝R(||(x_i,y_i)-(x_j,y_j)||), (6)

其中，N为控制点的个数，R为核函数：K_i,j为i行j列的子矩阵；

公式(5)等号右侧的o为3×1的零向量。Z为一维列向量，元素为高程数据r_c的高度值，定义如下：

Z＝[z₁ z₂ .. z_N]^T (8)

薄板样条曲面参数向量ω和a定义如下：

ω＝[ω₁,ω₂,...,ω_N]^T (9)

a＝[a₀,a₁,a₂]^T (10)

式中，ω为控制点权重参数，a为多项式系数；ω和a统称为薄板样条曲面参数向量；

步骤四三、求解出ω和a后，通过式(12)即得到插值曲面S定义为：

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

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：r_c＝[x_i,y_i,z_i]^T。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤七中所述的利用残差数据建立粗糙度图的具体操作步骤如下：

利用步骤六得到的最终插值曲面与原始高程数据r做差得到残差数据集合ε＝{ε_i|i＝1,2,3,...,N}；根据探测器岩石最大许可高度R_c，将粗糙度图Rscore定义为如下正规化函数：

其中MAD为：

Rscore图中任意一点的值越大就代表该处的粗糙度越大。其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下：

求出步骤六确定的插值曲面的梯度数据集合β＝{β_i|i＝1,2,3,..,N}；根据坡度许可阈值S_c，将坡度图Sscore定义为如下正规化函数：

式中norm(β_i)为：

Sscore中任意一点的值越大，代表该处坡度越陡。为了测试本实验的具体效果，利用matlab软件模拟火星表面局部地形(图2)，其中，图2中的真值曲面已知，如图3所示。以拟合优度R²判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验，结果如图4所示。可以看到，基于样条插值的障碍物提取方法具有很高的鲁棒性，当障碍物分布程度达到60％时，R²依然可以达到0.8以上。为了更好地展示算法的效果，又对真实火星地形(图5(a))建立了插值曲面，结果显示，本发明的方法可以很好的检测出岩石、坡度等障碍。图6(a)和(b)展示的是由残差值与梯度值经过步骤七、步骤八生成的障碍物图与坡度图，绝大多数的障碍物都被提取出来。其它步骤及参数与具体实施方式一至五之一相同。

Claims (4)

1.一种基于薄板样条插值的障碍物特征提取方法，其特征在，该方法具体是按照以下步骤进行的：

步骤一、根据行星着陆器在动力下降段的位置、姿态及扫描激光雷达参数获取Flash雷达视场内地形的三维高程数据；其中，Flash雷达视场内地形的三维高程数据指的是地形三维坐标点，包含两个水平方向的位置和高度信息：三维高程数据为地形的三维坐标x_i,y_i,z_i；i＝1,2,3,...N；

步骤二、对步骤一中获取的三维高程数据进行数字形态学处理，去除低值噪声点得到处理过的激光点数据；

步骤三、将经过步骤二处理过的激光点数据进行降采样处理；将三维高程数据进行隔点采样，即根据2倍的相邻采样点距离，以自下而上的方式建立数据金字塔，金字塔的层数为M，确定金字塔中任意m层的数据点为m+1层2*2邻域内的最小高程数据点；将最小高程数据点称之为控制点；m∈M；

步骤四、利用步骤三得到的每一层的控制点建立薄板样条插值曲面S；

步骤五、利用步骤四得到数据金字塔第m层薄板样条曲面S与数据金字塔第m+1层的控制点做差得到残差，将残差通过阈值处理，剔除大于阈值的控制点，更新控制点数据，其中，阈值为金字塔每一层数据的均值；

步骤六、从金字塔顶层，自上而下开始利用更新后的控制点数据重复步骤四和步骤五；最终得到数据金字塔最底层的更新后的控制点，利用最底层的更新后的控制点生成最终的薄板样条插值曲面；

步骤七：利用步骤六生成的最终的薄板样条插值曲面与原高程数据r做差得到残差数据，利用残差数据建立粗糙度图；

步骤八：求步骤六生成的最终的薄板样条插值曲面梯度，利用曲面梯度求地形坡度图；

步骤二中对步骤一中获取的三维高程数据进行数字形态学处理，去除低值噪声点得到处理过的激光点数据的具体操作步骤如下：

步骤二一、对原始高程数据r进行开运算，开运算是对数据进行腐蚀运算后在进行膨胀预算；窗口尺寸w_o为探测器最大尺寸的2倍；

其中，r_o为经过开运算处理后的高程数据；e代表腐蚀运算，具体定义为：

其中，r_er为经过腐蚀处理后的高程数据，z_i为高程数据点在窗口w_o下的高度值；代表膨胀运算，具体定义为：

步骤二二、对r_o再进行闭运算，窗口尺寸w_c为探测器最大尺寸：

r_c为经过闭运算滤波后的点云数据；

步骤四中所述的建立薄板样条曲面的具体操作步骤如下：

步骤四一、给定滤波后的高程数据r_c、正则参数λ、计算ω和a由下列线性方程组求得：

式中，K，C和O为子矩阵，O为3×3零矩阵，定义如下：

K＝{K_i,j}＝R(||(x_i,y_i)-(x_j,y_j)||), (6)

其中，N为控制点的个数，R为核函数：K_i,j为i行j列的子矩阵；

公式(5)等号右侧的o为3×1的零向量；Z为一维列向量，元素为高程数据r_c的高度值，定义如下：

Z＝[z₁ z₂ .. z_N]^T (8)

薄板样条曲面参数向量ω和a定义如下：

ω＝[ω₁,ω₂,...,ω_N]^T (9)

a＝[a₀,a₁,a₂]^T (10)

式中，ω为控制点权重参数，a为多项式系数；

步骤四二、求解出ω和a后，通过式(12)即得到插值曲面S定义为：

步骤七中所述的利用残差数据建立粗糙度图的具体操作步骤如下：

利用步骤六得到的最终插值曲面与原始高程数据r做差得到残差数据集合ε＝{ε_i|i＝1,2,3,...,N}；根据探测器岩石最大许可高度R_c，将粗糙度图Rscore定义为如下正规化函数：

其中MAD为：

2.根据权利要求1所述一种基于薄板样条插值的障碍物特征提取方法，其特征在于：r_c＝[x_i,y_i,z_i]^T。

3.根据权利要求1所述一种基于薄板样条插值的障碍物特征提取方法，其特征在于：步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下：

求出步骤六确定的插值曲面的梯度数据集合β＝{β_i|i＝1,2,3,..,N}；根据坡度许可阈值S_c，将坡度图Sscore定义为如下正规化函数：

式中norm(β_i)为：

4.根据权利要求2所述一种基于薄板样条插值的障碍物特征提取方法，其特征在于：步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下：

求出步骤六确定的插值曲面的梯度数据集合β＝{β_i|i＝1,2,3,..,N}；根据坡度许可阈值S_c，将坡度图Sscore定义为如下正规化函数：

式中norm(β_i)为：