CN106228521B - 一种基于薄板样条插值的障碍物特征提取方法 - Google Patents
一种基于薄板样条插值的障碍物特征提取方法 Download PDFInfo
- Publication number
- CN106228521B CN106228521B CN201610590103.7A CN201610590103A CN106228521B CN 106228521 B CN106228521 B CN 106228521B CN 201610590103 A CN201610590103 A CN 201610590103A CN 106228521 B CN106228521 B CN 106228521B
- Authority
- CN
- China
- Prior art keywords
- data
- thin
- point
- control point
- curved surface
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 43
- 230000004888 barrier function Effects 0.000 title claims abstract description 29
- 230000008569 process Effects 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims abstract description 3
- 239000011435 rock Substances 0.000 claims description 14
- 230000000007 visual effect Effects 0.000 claims description 7
- 230000003628 erosive effect Effects 0.000 claims description 5
- 238000005070 sampling Methods 0.000 claims description 4
- 230000010339 dilation Effects 0.000 claims description 3
- 238000005260 corrosion Methods 0.000 claims description 2
- 230000007797 corrosion Effects 0.000 claims description 2
- 238000004422 calculation algorithm Methods 0.000 abstract description 10
- 238000001914 filtration Methods 0.000 abstract description 9
- 238000001514 detection method Methods 0.000 abstract description 8
- 238000000605 extraction Methods 0.000 abstract description 3
- 230000000694 effects Effects 0.000 description 8
- 238000002474 experimental method Methods 0.000 description 8
- 238000000342 Monte Carlo simulation Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 238000005452 bending Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 230000008929 regeneration Effects 0.000 description 1
- 238000011069 regeneration method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/30—Polynomial surface description
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20036—Morphological image processing
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- Remote Sensing (AREA)
- Processing Or Creating Images (AREA)
Abstract
一种基于薄板样条插值的障碍物特征提取方法,本发明涉及障碍物特征提取方法。本发明是要解决现有的基于平面拟合的障碍检测算法在高海拔时不能有效地提取障碍物的问题,而提出的一种基于薄板样条插值的障碍物特征提取方法。该方法是通过首先滤掉噪声点,自下而上建立控制点金字塔。通过金字塔顶层控制点建立初始薄板样条插值曲面。求该曲面与金子塔下一层控制点的残差,通过阈值处理剔除不满足条件的控制点,更新控制点矩阵,通过剩余的控制点建立新的插值曲面,重复上述过程,最终生成原始数据点的插值曲面等步骤实现的。本发明应用于障碍物特征提取领域。
Description
技术领域
本发明涉及障碍特征检测方法,特别涉及一种基于薄板样条插值的障碍物特征提取方法。
背景技术
行星着陆器动力下降段障碍检测的主要目的是确保探测器能够安全地着陆到相对平坦的区域,以火星着陆为例,由于火星表面布满陨石坑,岩石,陡坡等地形,因此如何寻找出一块满足着陆要求的地形是最终着陆过程的关键所在。现有的方法都是在动力下降段末期(100米到500米)通过平面拟合求残差来探测障碍物。但是在较低的海拔范围内,一旦发现着陆区不具备安全着陆的可能,考虑到燃料的限制,探测器将很难做出大范围的机动。而且,海拔高度较低的情况下,探测器传感器的视角已经变得很小,不能够大范围的勘测地形,寻找合适的着陆点。因此需要探测器尽早地实施障碍检测操作,在动力下降段初期完成大范围的障碍检测,选择一块合适的着陆区域。待探测器在动力下降末期,在初期选择的着陆区域内选择合适的着陆点进行安全着陆。这样能够避免直接末期选择可能面临的上述问题。
然而在高海拔条件下,激光雷达获取的地形精度有限,激光脚点之间距离较大。传统的平面拟合方法并不适用。平面拟合的方法需要选择拟合窗口,窗口过大会过平滑障碍物;过小的窗口由于激光脚点距离较大,在窗口范围内采样点个数较少,受噪声影响较大。不能够有效地提取障碍物信息。
因此,如何在高海拔低分辨率,大采样距离的情况下,有效地检测障碍物是本发明重点解决的问题。
发明内容
本发明的目的是为了解决现有的基于平面拟合的障碍检测算法在高海拔时不能有效地提取障碍物的问题,而提出的一种基于薄板样条插值的障碍物特征提取方法。
上述的发明目的是通过以下技术方案实现的:
步骤一、根据行星着陆器在动力下降段的位置、姿态及扫描激光雷达参数获取Flash雷达视场内地形的三维高程数据;其中,Flash雷达视场内地形的三维高程数据指的是地形三维坐标点,包含两个水平方向的位置和高度信息:三维高程数据为地形的三维坐标xi,yi,zi;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所示。以拟合优度R2判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验,结果如图4所示。可以看到,基于样条插值的障碍物提取方法具有很高的鲁棒性,当障碍物分布程度达到60%时,R2依然可以达到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雷达视场内地形的三维高程数据指的是地形三维坐标点,包含两个水平方向的位置和高度信息:三维高程数据为地形的三维坐标xi,yi,zi;xi,yi为说水平方向的两个坐标,zi为高度方向上的坐标;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所示。以拟合优度R2判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验,结果如图4所示。可以看到,基于样条插值的障碍物提取方法具有很高的鲁棒性,当障碍物分布程度达到60%时,R2依然可以达到0.8以上。为了更好地展示算法的效果,又对真实火星地形(图5(a))建立了插值曲面,结果显示,本实施方式的方法可以很好的检测出岩石、坡度等障碍。图6(a)和(b)展示的是由残差值与梯度值经过步骤七、步骤八生成的障碍物图与坡度图,绝大多数的障碍物都被提取出来。
具体实施方式二:本实施方式与具体实施方式一不同的是:步骤二中对步骤一中获取的三维高程数据进行数字形态学处理,去除低值噪声点得到处理过的激光点数据的具体操作步骤如下:
步骤二一、对原始高程数据r进行开运算,开运算是对数据进行腐蚀运算后在进行膨胀预算;窗口尺寸wo为探测器最大尺寸的2倍;
ro=(r⊙wo)⊕wo (1)
其中,ro为经过开运算处理后的高程数据;⊙代表腐蚀运算,具体定义为:
其中,rer为经过腐蚀处理后的高程数据,zi为高程数据点在窗口wo下的高度值;⊕代表膨胀运算,具体定义为:
步骤二二、对ro再进行闭运算,闭运算是对数据先进行膨胀运算再进行腐蚀运算,窗口尺寸wc为探测器最大尺寸:
rc=(ro⊕wc)⊙wc (4)
rc为经过闭运算滤波后的点云数据;
激光雷达点高程据进行形态学开运算处理过后,将会剔除掉高值噪声点;
对经过开运算处理过的高程数据进行闭运算处理后,将会剔除点高程数据中的低值噪声点;
经过上述的开运算、闭运算处理后的数据点可以认为是已经不含有噪声点的数据了。其它步骤及参数与具体实施方式一相同。
具体实施方式三:本实施方式与具体实施方式一或二不同的是:步骤四中所述的建立薄板样条曲面的具体操作步骤如下:
步骤四一、给定控制点,薄板样条曲面为通过所有的控制点的弯曲最小的光滑曲面;光滑曲面弯曲程度由下述能量函数Es定义为:
式中,λ为正则参数,控制曲面的光滑度,S为光滑的薄板样条插值曲面;S={S(xi,yi)|i=1,2,3,..,N};
步骤四二、通过最小化Es后,给定滤波后的高程数据rc、正则参数λ。ω和a由下列线性方程组求得:
式中,K,C和O为子矩阵,O为3×3零矩阵,定义如下:
K={Ki,j}=R(||(xi,yi)-(xj,yj)||), (6)
其中,N为控制点的个数,R为核函数:Ki,j为i行j列的子矩阵;
公式(5)等号右侧的o为3×1的零向量。Z为一维列向量,元素为高程数据rc的高度值,定义如下:
Z=[z1 z2 .. zN]T (8)
薄板样条曲面参数向量ω和a定义如下:
ω=[ω1,ω2,...,ωN]T (9)
a=[a0,a1,a2]T (10)
式中,ω为控制点权重参数,a为多项式系数;ω和a统称为薄板样条曲面参数向量;
步骤四三、求解出ω和a后,通过式(12)即得到插值曲面S定义为:
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四:本实施方式与具体实施方式一至三之一不同的是:rc=[xi,yi,zi]T。其它步骤及参数与具体实施方式一至三之一相同。
具体实施方式五:本实施方式与具体实施方式一至四之一不同的是:步骤七中所述的利用残差数据建立粗糙度图的具体操作步骤如下:
利用步骤六得到的最终插值曲面与原始高程数据r做差得到残差数据集合ε={εi|i=1,2,3,...,N};根据探测器岩石最大许可高度Rc,将粗糙度图Rscore定义为如下正规化函数:
其中MAD为:
Rscore图中任意一点的值越大就代表该处的粗糙度越大。其它步骤及参数与具体实施方式一至四之一相同。
具体实施方式六:本实施方式与具体实施方式一至五之一不同的是:步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下:
求出步骤六确定的插值曲面的梯度数据集合β={βi|i=1,2,3,..,N};根据坡度许可阈值Sc,将坡度图Sscore定义为如下正规化函数:
式中norm(βi)为:
Sscore中任意一点的值越大,代表该处坡度越陡。为了测试本实验的具体效果,利用matlab软件模拟火星表面局部地形(图2),其中,图2中的真值曲面已知,如图3所示。以拟合优度R2判定算法的鲁棒性。共进行了2500次的蒙特卡洛模拟实验,结果如图4所示。可以看到,基于样条插值的障碍物提取方法具有很高的鲁棒性,当障碍物分布程度达到60%时,R2依然可以达到0.8以上。为了更好地展示算法的效果,又对真实火星地形(图5(a))建立了插值曲面,结果显示,本发明的方法可以很好的检测出岩石、坡度等障碍。图6(a)和(b)展示的是由残差值与梯度值经过步骤七、步骤八生成的障碍物图与坡度图,绝大多数的障碍物都被提取出来。其它步骤及参数与具体实施方式一至五之一相同。
Claims (4)
1.一种基于薄板样条插值的障碍物特征提取方法,其特征在,该方法具体是按照以下步骤进行的:
步骤一、根据行星着陆器在动力下降段的位置、姿态及扫描激光雷达参数获取Flash雷达视场内地形的三维高程数据;其中,Flash雷达视场内地形的三维高程数据指的是地形三维坐标点,包含两个水平方向的位置和高度信息:三维高程数据为地形的三维坐标xi,yi,zi;i=1,2,3,...N;
步骤二、对步骤一中获取的三维高程数据进行数字形态学处理,去除低值噪声点得到处理过的激光点数据;
步骤三、将经过步骤二处理过的激光点数据进行降采样处理;将三维高程数据进行隔点采样,即根据2倍的相邻采样点距离,以自下而上的方式建立数据金字塔,金字塔的层数为M,确定金字塔中任意m层的数据点为m+1层2*2邻域内的最小高程数据点;将最小高程数据点称之为控制点;m∈M;
步骤四、利用步骤三得到的每一层的控制点建立薄板样条插值曲面S;
步骤五、利用步骤四得到数据金字塔第m层薄板样条曲面S与数据金字塔第m+1层的控制点做差得到残差,将残差通过阈值处理,剔除大于阈值的控制点,更新控制点数据,其中,阈值为金字塔每一层数据的均值;
步骤六、从金字塔顶层,自上而下开始利用更新后的控制点数据重复步骤四和步骤五;最终得到数据金字塔最底层的更新后的控制点,利用最底层的更新后的控制点生成最终的薄板样条插值曲面;
步骤七:利用步骤六生成的最终的薄板样条插值曲面与原高程数据r做差得到残差数据,利用残差数据建立粗糙度图;
步骤八:求步骤六生成的最终的薄板样条插值曲面梯度,利用曲面梯度求地形坡度图;
步骤二中对步骤一中获取的三维高程数据进行数字形态学处理,去除低值噪声点得到处理过的激光点数据的具体操作步骤如下:
步骤二一、对原始高程数据r进行开运算,开运算是对数据进行腐蚀运算后在进行膨胀预算;窗口尺寸wo为探测器最大尺寸的2倍;
其中,ro为经过开运算处理后的高程数据;e代表腐蚀运算,具体定义为:
其中,rer为经过腐蚀处理后的高程数据,zi为高程数据点在窗口wo下的高度值;代表膨胀运算,具体定义为:
步骤二二、对ro再进行闭运算,窗口尺寸wc为探测器最大尺寸:
rc为经过闭运算滤波后的点云数据;
步骤四中所述的建立薄板样条曲面的具体操作步骤如下:
步骤四一、给定滤波后的高程数据rc、正则参数λ、计算ω和a由下列线性方程组求得:
式中,K,C和O为子矩阵,O为3×3零矩阵,定义如下:
K={Ki,j}=R(||(xi,yi)-(xj,yj)||), (6)
其中,N为控制点的个数,R为核函数:Ki,j为i行j列的子矩阵;
公式(5)等号右侧的o为3×1的零向量;Z为一维列向量,元素为高程数据rc的高度值,定义如下:
Z=[z1 z2 .. zN]T (8)
薄板样条曲面参数向量ω和a定义如下:
ω=[ω1,ω2,...,ωN]T (9)
a=[a0,a1,a2]T (10)
式中,ω为控制点权重参数,a为多项式系数;
步骤四二、求解出ω和a后,通过式(12)即得到插值曲面S定义为:
步骤七中所述的利用残差数据建立粗糙度图的具体操作步骤如下:
利用步骤六得到的最终插值曲面与原始高程数据r做差得到残差数据集合ε={εi|i=1,2,3,...,N};根据探测器岩石最大许可高度Rc,将粗糙度图Rscore定义为如下正规化函数:
其中MAD为:
2.根据权利要求1所述一种基于薄板样条插值的障碍物特征提取方法,其特征在于:rc=[xi,yi,zi]T。
3.根据权利要求1所述一种基于薄板样条插值的障碍物特征提取方法,其特征在于:步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下:
求出步骤六确定的插值曲面的梯度数据集合β={βi|i=1,2,3,..,N};根据坡度许可阈值Sc,将坡度图Sscore定义为如下正规化函数:
式中norm(βi)为:
4.根据权利要求2所述一种基于薄板样条插值的障碍物特征提取方法,其特征在于:步骤八中所述的利用曲面梯度求地形坡度图的具体操作步骤如下:
求出步骤六确定的插值曲面的梯度数据集合β={βi|i=1,2,3,..,N};根据坡度许可阈值Sc,将坡度图Sscore定义为如下正规化函数:
式中norm(βi)为:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610590103.7A CN106228521B (zh) | 2016-07-25 | 2016-07-25 | 一种基于薄板样条插值的障碍物特征提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610590103.7A CN106228521B (zh) | 2016-07-25 | 2016-07-25 | 一种基于薄板样条插值的障碍物特征提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106228521A CN106228521A (zh) | 2016-12-14 |
CN106228521B true CN106228521B (zh) | 2018-12-11 |
Family
ID=57531692
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610590103.7A Active CN106228521B (zh) | 2016-07-25 | 2016-07-25 | 一种基于薄板样条插值的障碍物特征提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106228521B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109903352B (zh) * | 2018-12-24 | 2021-03-30 | 中国科学院遥感与数字地球研究所 | 一种卫星遥感影像大区域无缝正射影像制作方法 |
CN110443770A (zh) * | 2019-08-12 | 2019-11-12 | 重庆市地理信息和遥感应用中心(重庆市测绘产品质量检验测试中心) | 基于离散粗糙度估计的机载激光点云数据噪声检测方法 |
CN112764005A (zh) * | 2021-01-05 | 2021-05-07 | 哈尔滨工业大学 | 一种结合形态学滤波的Gm-APD激光雷达低信噪比回波数据重构方法 |
CN114332631B (zh) * | 2022-01-12 | 2023-04-18 | 中铁二院工程集团有限责任公司 | 一种适用于山区危岩落石的LiDAR点云数据提取方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101908235A (zh) * | 2010-06-23 | 2010-12-08 | 苏州科技学院 | B样条曲面建模新方法 |
CN102324106A (zh) * | 2011-06-02 | 2012-01-18 | 武汉大学 | 一种顾及地表光谱信息的sfs三维重建加密稀疏dem方法 |
WO2012055752A1 (fr) * | 2010-10-27 | 2012-05-03 | Societe De Technologie Michelin | Methode de pre-traitement d'une image tridimensionnelle de la surface d'un pneumatique a l'aide de deformations b-spline successives |
CN102663402A (zh) * | 2012-04-21 | 2012-09-12 | 西北工业大学 | 基于修正扩展形态学算子的高光谱遥感图像端元提取方法 |
CN103823981A (zh) * | 2014-02-28 | 2014-05-28 | 武汉大学 | 一种数字高程模型辅助的卫星影像区域网平差方法 |
-
2016
- 2016-07-25 CN CN201610590103.7A patent/CN106228521B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101908235A (zh) * | 2010-06-23 | 2010-12-08 | 苏州科技学院 | B样条曲面建模新方法 |
WO2012055752A1 (fr) * | 2010-10-27 | 2012-05-03 | Societe De Technologie Michelin | Methode de pre-traitement d'une image tridimensionnelle de la surface d'un pneumatique a l'aide de deformations b-spline successives |
CN102324106A (zh) * | 2011-06-02 | 2012-01-18 | 武汉大学 | 一种顾及地表光谱信息的sfs三维重建加密稀疏dem方法 |
CN102663402A (zh) * | 2012-04-21 | 2012-09-12 | 西北工业大学 | 基于修正扩展形态学算子的高光谱遥感图像端元提取方法 |
CN103823981A (zh) * | 2014-02-28 | 2014-05-28 | 武汉大学 | 一种数字高程模型辅助的卫星影像区域网平差方法 |
Non-Patent Citations (1)
Title |
---|
"A Robust Filtering Algorithm of LiDAR for DTM Extraction";明洋 等;《Advanced Materials Research》;20130904;第765-767卷;第639页第一行至640页最后一行 * |
Also Published As
Publication number | Publication date |
---|---|
CN106228521A (zh) | 2016-12-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106228521B (zh) | 一种基于薄板样条插值的障碍物特征提取方法 | |
CN104867126B (zh) | 基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法 | |
Li et al. | Rainfall and earthquake-induced landslide susceptibility assessment using GIS and Artificial Neural Network | |
CN109633762B (zh) | 基于正弦函数的相关性约束条件联合反演重磁数据的方法 | |
CN109598243B (zh) | 一种月球表面安全着陆区选择方法及*** | |
Wieczorek et al. | Automatic relief classification versus expert and field based landform classification for the medium-altitude mountain range, the Sudetes, SW Poland | |
CN110110675A (zh) | 一种融合边缘信息的小波域分形红外卷云检测方法 | |
CN110673138B (zh) | 一种基于奇异值分解和模糊c均值法的探地雷达图像处理方法 | |
Hui et al. | An active learning method for DEM extraction from airborne LiDAR point clouds | |
CN107393004A (zh) | 一种获取输电线走廊中建筑物拆迁量的方法及装置 | |
CN116187168A (zh) | 基于神经网络-重力信息小波分解提高水深反演精度方法 | |
CN108460422B (zh) | 基于深度分布特征的海底地貌类型识别方法 | |
Wang et al. | A progressive morphological filter for point cloud extracted from UAV images | |
Troglio et al. | Crater detection based on marked point processes | |
CN105205450B (zh) | 一种基于非规则标识点过程的sar图像目标提取方法 | |
Pelletier et al. | Tectonic and structural control of fluvial channel morphology in metamorphic core complexes: The example of the Catalina-Rincon core complex, Arizona | |
CN111737882B (zh) | 复杂月面接近段实现自主障碍规避的着陆区选取方法 | |
Novák et al. | The Potential and Implications of Automated Pre-Processing of LiDAR-Based Digital Elevation Models for Large-Scale Archaeological Landscape Analysis | |
CN114692441A (zh) | 一种黄土滑坡稳定性预测方法、电子设备以及存储介质 | |
Di Salvo et al. | Full-waveform terrestrial laser scanning for extracting a high-resolution 3D topographic model: A case study on an area of archaeological significance | |
Moradi et al. | Individual tree of urban forest extraction from very high density LiDAR data | |
Székely et al. | Automated recognition of quasi‐planar ignimbrite sheets as paleosurfaces via robust segmentation of digital elevation models: an example from the Central Andes | |
Peyk-Herfeh et al. | Evaluation of adaptive boosting and neural network in earthquake damage levels detection | |
Anderson | Mapping Relict Charcoal Hearths in the Northeast US Using Deep Learning Convolutional Neural Networks and LIDAR Data | |
Xu et al. | Point cloud segmentation of gully based on characteristic difference using airborne lidar data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |