CN105761312B - 一种微地形表面重建方法 - Google Patents

Abstract

本发明公开了一种基于激光雷达的微地形重建方法，步骤如下：(1)激光扫描获得海量、离散激光雷达点云数据后，进行不同站点间的点云数据拼接，并分离地面点和非地面点；(2)进行点云数据抽稀、点云孔洞区或点云密度极低区的边界点检测；(3)计算微地形的分形特征，并确定无标度区；(4)以带法线控制的随机中点位移法对点云孔洞区和点云密度极低区进行插值；(5)利用贪婪投影三角化方法，基于地面点进行微地形重建。

Description

一种微地形表面重建方法

技术领域

本发明涉及遥感技术与测绘学技术领域，更具体涉及一种微地形表面重建方法。

背景技术

微地形是指由自然力或人类活动等作用形成的大小不等、性状各异的局部地形，具有表面细节结构复杂、地形起伏随机等特点。对于植被来说，微地形造成了水热养分等条件的差异变化，这种变化影响着地表植被的生长；对于园林景观来说，依照天然地貌或人为造出的微地形仿自然界中的起伏变化地势，决定园林中的水体、园林构筑物、树木、花草等的分布。因此，进行微地形的精细重建对于微地形环境研究、微地形改造利用具有非常重要的意义。

目前，微地形重建方法一般采用常规地形建模方法，其模型精度与计算机硬件系统、建模算法、模型三维显示算法以及地形原始数据相关。现阶段，计算机硬件以及算法已经不再是重建的障碍，但是原始地形数据的信息容量和精确度仍是制约精确地形模型建立的瓶颈。传统人工勘测方法得到的数据难以详尽描述地形微小细节变化，且人工采集数据较为费时、费力，难以采集到人难以到达区域的地形数据，数据精度也有赖于勘测人员的业务熟练度。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是如何基于大量数据，准确、自动化测量微地形数据。

(二)技术方案

为了解决上述技术问题，本发明提供了一种微地形表面重建方法，所述方法包括以下步骤：

S1、根据激光雷达获取点云数据，并对多站点的所述点云数据进行配准、拼接；将所述点云数据分类为地面点云和非地面点云；

S2、设置抽稀密度，对所述地面点云中点云密度过高的区域进行点云抽稀；利用规则格网掩膜检测算法，检测点云孔洞区和点云密度极低区；

S3、采用双对数线性拟合方法求取研究区域地形的无标度区，并基于该无标度区尺度求算地形分形维数D、表示地形表面粗糙度参数H和地形表面起伏度参数σ等地形分形参数；

S4、利用带法线控制的随机中点位移法以及所述地形分形维数D、表示地形表面粗糙度参数H、地形表面起伏度参数σ，对所述点云孔洞区和点云密度极低区进行插值；

S5、利用贪婪投影三角化方法，对步骤S4得到的地面点云进行微地形表面重建，并通过三角格网曲面模型表达微地形。

优选地，所述步骤S1中，对所述点云数据分类利用渐进不规则三角网加密法进行；

所述点云数据拼接具体为：将不同站点的点云坐标转换为同一坐标系下的坐标。

优选地，所述步骤S2中，检测点云孔洞区和点云密度极低区的具体步骤S21中，规则格网边长λ为：

其中，P为所述步骤S1中得到的地面点云的点集密度；

步骤S22中，每一所述点云对应的规则格网坐标的计算公式为：

其中，λ为所述规则格网边长，x_i,y_i为所述步骤S1中得到的地面点云的点集中点在xy平面投影坐标，x_min为x_i,取值中最小值，y_min为y_i取值中的最小值，i的取值范围为[1，n]，n表示所述地面点云的总数。

步骤S23中，检测所述点云孔洞区或点云密度极低区的边界点具体为：

遍历所述规则格网点集，当前检测点设为(x_ij,y_ij)，分别检测其临近四个点(x_ij+r,x_ij)，(x_ij-r,y_ij)，(x_ij,y_ij+r)，(x_ij,y_ij-r)的r邻域内的规则格网点云数量，r为搜索半径；

设点(x_ij-r,y_ij)的r邻域内规则格网点云数量为n₁，(x_ij+r,x_ij)的r邻域内规则格网点云数量为n₂，(x_ij,y_ij-r)的r邻域内规则格网点云数量为n₃，(x_ij,y_ij+r)的r邻域内规则格网点云数量为n₄，阈值设定为d_1,d₂，如果满足下面两情况之一，则将点(x_ij,y_ij)判定为边界点：①n₁＜d₁或者n₂＜d₁，且|n₁-n₂|＜d₂；②n₃＜d₁或者n₄＜d₁，且|n₃-n₄|＜d₂。

步骤S24中，连接所述边界点成为边界线具体为：

将整个检测区域沿着竖直坐标方向划分为N层，每一层的间距设置为所述规则格网边长为λ，每层中上端最后一个点为p_i-start，下端最后一个点为p_i-end；

检测在当前层存在的所有的所述边界点：分别计算各所述边界点与p_i-start和p_i-end的距离，分别以s_i-start和s_i-end表示，如果s_i-start小于s_i-end则将当前检测的所述边界点插入p_i-start；否则，将其插入p_i-end后；

按照层数递减的顺序依次搜索完所述整个检测区域，所有的所述边界点按照其邻接关系排成一个有序的点集，完成所述点云孔洞区和点云密度极低区的提取。

优选地，所述步骤S3具体包括以下步骤：

S31、计算所述计算地形分形维数D、表示地形表面粗糙度参数H和地形表面起伏度参数σ，

对于真实地形表面X(x,y)，其任意布朗运动公式为：

对Hlog(Δx²+Δy²)^1/2取不同的值，即取不同的(Δx,Δy)，经最小二乘法计算即可求出H和σ的值，地形分形维数D＝3-H；

S32、采用双对数线性拟合方法确定表征地形满足分形特征尺度范围的无标度区；根据分形理论中对于无标度区计算的要求，以logE[|X(x+Δx,y+Δy)-X(x,y)|]为纵轴坐标，以log(Δx²+Δy²)^1/2)为横轴坐标拟合曲线，求取所述曲线中线性相关性最好的一段所对应的[log(Δx²+Δy²)^1/2 _min，log(Δx²+Δy²)^1/2 _max]作为无标度区，其中log(Δx₂+Δy₂)^1/2 _min为步骤S1中的地面点云的分辨率，log(Δx²+Δy²)^1/2 _max为最终的目标。

优选地，所述步骤S4中，插入点云的坐标为：

法线控制的基本思想通过已知点的法线方法推知新插入点的法线方向，利用法线对插值过程进行控制，使新插入的点法线与临近点法线存在一致性。这使得经过插值后生成的曲面在保持其原始表面特性的同时，也使插值点与周围点存在值的传递关系，能够有效的解决传统的中点位移法会产生折痕效应的问题。

优选地，所述步骤S5具体为：

将三维点通过法线投影到某一平面；然后，对投影得到的点云，采用Delaunay空间区域增长方法作平面内点的三角化，从而得到各点的连接关系；最后，根据投影点云的连接关系确定各原始三维点间的拓扑连接，所得三角格网即为重建的曲面模型。

(三)有益效果

本发明提供了一种微地形表面重建方法，具有以下的优点：通过点云数据抽稀、点云孔洞区或点云密度极低区的边界点检测，解决微地形建模中数据量与建模精度之间的矛盾；并以带法线控制的随机中点位移法，有效的解决传统的中点位移法所产生的“折痕问题”，实现了精确、自动化建模。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的一种微地形重建方法的流程图；

图2a、2b、2c、2d为本发明的点云分类结果剖面图；

图3a、3b、3c为本发明点云孔洞区和点云密度极低区检测结果示意图；

图3d为本发明点云孔洞区和点云密度极低区检测结果与原始点云重叠效果图；

图4为本发明的一个子点云数据集双对数坐标直线拟合结果图。

具体实施方式

以北京市门头沟区王平镇内一处煤矿复垦区微地形重建过程为例，如图1所示，本发明的具体实施方法如下，如图1所示：

步骤1、在研究区内多站架设，激光扫描获取点云数据，对多站点的点云数据进行配准、拼接，并将点云数据分类为地面点云和非地面点云。其中对点云的分类结果如图2a、2b、2c、2d所示。

根据要保证全覆盖且相邻站点之间有一定比例的重复扫描区域、尽量减少遮挡、易于选择同名点、扫描精度适中的原则，确定激光扫描方案；在相邻站点的激光雷达之间放置4-6个标靶球作为点云数据拼接同名点，所述标靶球位置选择遵循尽量全覆盖、分布于两扫描站点之间的公共区域等原则；设置合适的水平扫描分辨率和垂直扫描分辨率，依次对选取的站点进行激光扫描，直至完成实验区的全覆盖扫描。

激光扫描完成后首先进行多站点的点云数据拼接。由于激光雷达自身不带GPS定位系统，因此扫描获得的多站原始点云坐标为彼此独立的三维坐标系，需要通过点云拼接将不同站点扫描所得的点云整合到同一个坐标系下，具体做法为：以一个站点作为基准站，通过选取控制点(即标靶球所在位置)构建旋转矩阵将多个站点的坐标进行统一；然后，采用渐进不规则三角网加密算法进行激光雷达点云分类，通过不断迭代构建不规则三角格网，判断满足阈值点的类别，将进行配准、拼接后的点云分类为地面点与非地面点。

本发明中激光扫描使用的激光雷达仪器是由奥地利RIEGL公司2012年推出的、具有全波形回波技术的新型脉冲式地基激光雷达扫描仪RIEGL VZ-1000，扫描过程中激光雷达水平扫描分辨率设置为0.02°，垂直扫描分辨率设置为0.03°。

步骤2、设置抽稀密度，对所述地面点云中点云密度过高的区域进行点云抽稀；利用规则格网掩膜检测算法，检测点云孔洞区和点云密度极低区，如图3a、3b、3c所示。图3d为本发明点云孔洞区和点云密度极低区检测结果与原始点云重叠效果图。

为了兼顾距离激光扫描站点较远处的扫描精度，扫描站点周围10m范围内的激光点云密度可达每平方米上万个扫描点，微地形建模中并不需要如此大的点云密度，因此需要对站点周围一定范围内的点云进行抽稀，以使整个扫描区具有较为一致的点云密度。设置抽稀密度为p，即每平方米有p个点，对点云密度过高的区域进行点云抽稀。

设置规则格网边长λ为0.045m，检测由于遮挡等原因造成的点云孔洞区和点云密度极低的区域，边界点检测过程中，阈值d₁,d₂的判定与x轴邻域内点云差值或y轴邻域内点云数量相关，根据已知点云的密度为p，确定(r＝0.1,d₁＝15,d₂＝2)、(r＝0.3,d₁＝130,d₂＝10)、(r＝0.5,d₁＝300,d₂＝35)为比较优的参数组合，图3a、3b、3c为利用这三种参数组合检测出的点云孔洞区和点云密度极低区的边界，可以看出：当r取0.1时检测出的点云孔洞区和点云密度极低区边界点中引入了过多的噪点，许多非点云孔洞点也被误判为点云孔洞区或点云密度极低区边界点；当r取0.5时可以明显看到一些较小的点云孔洞区和点云密度极低区的边界并未被检出，且检测边界点不完整；当r取0.3时，检测出的点云孔洞区和点云密度极低区边界点较为理想，能基本提取出完整的点云孔洞区和点云密度极低区边界点。图3d为r＝0.3时点云孔洞区和点云密度极低区检测结果与扫描得到的地面点云重叠度。

步骤3、确定表征地形满足分形特征尺度范围的无标度区，并计算地形分形维数D、表示地形表面粗糙度参数H和地形表面起伏度参数σ；

地形的分形特征只在一定尺度下满足，我们将这个地形满足分形特征的尺度范围称为无标度区，本发明采用双对数线性拟合方法确定表征地形满足分形特征尺度范围的无标度区，图4为双对数坐标直线拟合结果，纵轴坐标为logE[|X(x+Δx,y+Δy)-X(x,y)|]，横轴坐标为log(Δx²+Δy²)^1/2)。从图4可以看出：当log(Δx²+Δy²)^1/2)取值为(1,2.9)时(logE[|X(x+Δx,y+Δy)-X(x,y)|],log(Δx²+Δy²)^1/2)有较好的线性相关性，因此取(1,2.9)作为无标度区，计算得到地形分形维数D为2.12、地形表面粗糙度参数H为0.82、地形表面起伏度参数σ为1.58。

步骤4、利用带法线控制的随机中点位移法，对所述点云孔洞区和点云密度极低区进行插值。

法线控制的基本思想通过已知点的法线方法推知新插入点的法线方向，利用法线对插值过程进行控制，使新插入的点法线与临近点法线存在一致性。这使得经过插值后生成的曲面在保持其原始表面特性的同时，也使插值点与周围点存在值的传递关系，能够有效的解决传统的中点位移法会产生折痕效应的问题。插入点云的坐标为：

步骤5、利用贪婪投影三角化方法，对步骤S4得到的地面点云进行微地形表面重建，并通过三角格网曲面模型表达微地形。

利用贪婪投影三角化方法是局部进行的，首先沿一个随机点云的法线将该点云投影到局部二维坐标平面内并连接其他悬空点，然后再进行下一点云，从而将越来越多的邻居点云包含进三角化网格，使网格逐渐增长，形成三角格网曲面。

贪婪投影三角化方法是通过处理一系列可以使网格“生长扩大”的点(边缘点)，延伸这些点直到所有符合几何正确性和拓扑正确性的点都被连上，之所以称其为贪婪投影三角化方法，是因为它总是将一个点与尽可能多的邻居点云连接起来。

在上述三角化过程中用到了基于Delaunay的空间区域增长算法，该方法通过选取一个点云三角片作为初始曲面，不断扩张曲面边界，最后形成一张完整的三角网格曲面，最后根据投影点云的连接关系确定各原始三维点云间，即步骤(1)中的地面点云的拓扑连接，所得三角格网即为重建得到的曲面模型。在建网的过程中，每个点云有可能被赋予这四种状态：free、fringe、boundary、completed，还有一种特殊的状态none，用来检查某些值的有效性，这四种状态的具体意义为：

free：未加入任何三角形。

fringe：已加入一些三角形，还有加入另一些三角形的可能性。

boundary：已位于三角化网格的边界，即已加入若干三角形，不能加入更多三角形，但未被这些三角形完全包围。

completed：已加入若干三角形，不能加入更多的三角形，并被这些三角形完全包围。

本发明的方法与现有技术相比，具有以下的优点：通过点云数据抽稀和点云空洞边界检测，解决微地形建模中数据量与建模精度之间的矛盾；并以所研发的带法线控制的随机中点位移法，有效的解决传统的中点位移法所产生的“折痕问题”。

以上实施方式仅用于说明本发明，而非对本发明的限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行各种组合、修改或者等同替换，都不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于激光雷达的微地形重建方法，其特征在于，所述方法包括以下步骤：

S1、激光扫描获取激光雷达点云数据，并对多站点的所述点云数据进行配准、拼接；将配准、拼接后的点云数据分类为地面点云和非地面点云；

S2、设置抽稀密度，对所述地面点云中点云密度过高的区域进行点云抽稀；利用规则格网掩膜检测算法，检测点云孔洞区和点云密度极低区；

S3、采用双对数线性拟合方法求取研究区域地形的无标度区，并基于该无标度区尺度求算地形分形维数D、表示地形表面粗糙度参数H和地形表面起伏度参数σ等地形分形参数；

S4、利用带法线控制的随机中点位移法以及计算得到的地形分形参数，对所述点云孔洞区和点云密度极低区进行插值；

S5、利用贪婪投影三角化方法，对步骤S4得到的地面点云进行微地形表面重建，并通过三角格网曲面模型表达微地形；

其中，步骤S4中，在已知两点或者多点中点处插入一个新的点，新插入点云的坐标为：

法线控制的基本思想通过已知点的法线方法推知新插入点的法线方向，利用法线对插值过程进行控制，使新插入的点法线与临近点法线存在一致性；这使得经过插值后生成的曲面在保持其原始表面特性的同时，也使插值点与周围点存在值的传递关系，能够有效的解决传统的中点位移法会产生折痕效应的问题。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1中，对所述点云数据分类利用渐进不规则三角网加密法进行。

3.根据权利要求1所述的方法，其特征在于，所述步骤S2中，检测点云孔洞区和点云密度极低区的具体步骤为：

S21、根据所述步骤S1中得到的地面点云密度确定规则格网边长；

S22、利用所述规则格网边长，计算各点云对应的规则格网坐标，构成规则格网点集；

S23、根据所述规则格网点集的拓扑关系检测所述点云孔洞区或点云密度极低区的边界点；

S24、连接所述边界点成为边界线，完成所述点云孔洞区和点云密度极低区的提取。

4.根据权利要求1所述的方法，其特征在于，所述步骤S3为：

S31、采用双对数线性拟合方法确定表征地形满足分形特征尺度范围的无标度区；根据分形理论中对于无标度区计算的要求，以logE[|X(x+Δx,y+Δy)-X(x,y)|]为纵轴坐标，以log(Δx²+Δy²)^1/2)为横轴坐标拟合曲线，求取所述曲线中线性相关性最好的一段所对应的[log(Δx²+Δy²)^1/2 _min,log(Δx²+Δy²)^1/2 _max]作为无标度区，其中log(Δx²+Δy²)^1/2 _min为所述步骤S1中的地面点云的分辨率，log(Δx²+Δy²)^1/2 _max为最终的目标；

S32、计算所述计算地形分形维数D、表示地形表面粗糙度参数H和地形表面起伏度参数σ；

对于真实地形表面X(x,y)，其任意布朗运动公式为：

对Hlog(Δx²+Δy²)^1/2取不同的值，即取不同的(Δx,Δy)，经最小二乘法计算即可求出H和σ的值，地形分形维数D＝3-H。

5.根据权利要求1所述的方法，其特征在于，所述步骤S5具体为：

将所述步骤S4中得到的三维点云通过法线投影到某一平面；然后，对投影得到的点云，采用Delaunay空间区域增长方法作平面内点的三角化，从而得到各点的连接关系；最后，根据投影点云的连接关系确定所述步骤S1中的地面点云间的拓扑连接，所得三角格网即为重建的曲面模型。