CN114332631B - 一种适用于山区危岩落石的LiDAR点云数据提取方法 - Google Patents

Abstract

本发明公开了一种适用于山区危岩落石的LiDAR点云数据提取方法，包括：步骤1：对山区扫描获取点云数据；步骤2：设定格网尺寸，对点云数据构建虚拟格网；步骤3：对点云数据进行去噪处理，再在虚拟格网中创建索引号，确定点云数据最大行号；步骤4：根据最大行号计算结构元素，对点云数据进行形态学开运算，得到点云数据的高差变化值；对点云数据进行曲面拟合，得到点云数据的地面数字模型，进而得到地形起伏度值；步骤5：根据高差变化值及地形起伏度值判断点云数据是否为地面点；步骤6：将结构元素与预设阈值比较，若小于阈值，则结束，否则返回步骤2更新格网尺寸。本发明对点云数据进行形态学开运算和曲面拟合，能更加准确得到地面点。

Description

一种适用于山区危岩落石的LiDAR点云数据提取方法

技术领域

本发明涉及航空摄影测量技术领域，特别涉及一种适用于山区危岩落石的LiDAR点云数据提取方法。

背景技术

在进行铁路线路勘查时，由于山区的地形起伏复杂多变，且一些区域存在着极大的落差，山区地表又覆盖有茂密的植被，常规的地面勘探与调绘无法实行，遥感技术手段也受制于气候条件和植被覆盖因素的制约，使工作人员无法准确了解高植被覆盖区下的地形情况。

发明内容

本发明的目的在于克服现有技术中所存在的上述不足，提供一种适用于山区危岩落石的LiDAR点云数据提取方法，通过利用机载激光雷达对山区进行扫描获取LiDAR点云数据，对LiDAR点云数据进行有效点云数据的提取，从而根据有效点云数据获得山区真实地形，了解山区的危岩落石情况。

为了实现上述发明目的，本发明提供了以下技术方案：

一种适用于山区危岩落石的LiDAR点云数据提取方法，包括以下步骤：

步骤1：利用机载激光雷达对山区进行扫描获取LiDAR点云数据；

步骤2：设定格网尺寸，对所述LiDAR点云数据构建虚拟格网；

步骤3：对所述LiDAR点云数据进行去噪处理，再对去噪处理后的所述LiDAR点云数据在所述虚拟格网中创建索引号，确定所述LiDAR点云数据在所述虚拟格网中的最大行号；

步骤4：根据所述最大行号计算结构元素，根据所述结构元素对所述LiDAR点云数据进行形态学开运算，得到形态学开运算前后的所述LiDAR点云数据的高差变化值；同时，根据所述结构元素对所述LiDAR点云数据进行曲面拟合，得到所述LiDAR点云数据的地面数字模型，根据所述地面数字模型得到所述LiDAR点云数据的地形起伏度值；

步骤5：利用所述步骤4所得地形起伏度值根据地形坡度公式计算高差阈值，并计算所述LiDAR点云数据的所述高差变化值与所述地形起伏度值的差，得到差值，将所述差值与所述高差阈值进行比较，根据比较结果判断LiDAR点云数据是否为地面点，若是，则提取所述LiDAR点云数据；若否，则剔除所述LiDAR点云数据；

步骤6：将所述结构元素与预设的结构元素阈值进行比较，若所述结构元素小于所述结构元素阈值，则输出提取到的LiDAR点云数据集；若否，返回步骤2更新设定的格网尺寸，重复步骤2至步骤5，以更新所提取的地面点。

本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，通过对机载激光雷达扫描山区得到的LiDAR点云数据，然后通过构建虚拟格网对LiDAR点云数据进行数据组织并进行去噪处理，再通过确定点云数据的格网最大行号计算得到结构元素，基于结构元素对LiDAR点云数据进行形态学开运算得到形态学开运算前后的LiDAR点云数据的高差变化值，从而能够初步判断地物点和地面点，在对LiDAR点云数据进行形态学开运算的同时对LiDAR点云数据进行曲面拟合，引入地形起伏度，通过地形起伏度计算得到高差阈值，通过将LiDAR点云数据进行形态学开运算前后的高差变化值与进行曲面拟合得到的地形起伏度做差，将差值与高差阈值作比较，从而能够再次对地物点和地面点进行判断，减少误判，且利用形态学开运算与曲面拟合两者相结合对山区的LiDAR点云数据进行提取，能够对地面凸起等进行有效保护，更加准确的得到地面点，从而得到山区的真实地形；使得工作人员能够通过保留的点云数据解算出被覆盖在茂密植被和复杂地形环境中的危岩落石，为工程施工提供安全保障。

优选的，所述虚拟格网是以规则格网为基础进行建立，是对所述LiDAR点云数据在XY平面进行格网剖分但不进行内插的数据组织形式。

优选的，所述步骤3中，采用经验模态分解的方法对所述LiDAR点云数据进行去噪处理。

优选的，所述步骤4中，对所述LiDAR点云数据进行形态学开运算，具体包括：基于所述结构元素对所述LiDAR点云数据先进行腐蚀运算，再进行膨胀运算。

优选的，所述步骤4中，采用克里金插值的方法对所述LiDAR点云数据进行曲面拟合。

优选的，所述步骤5中，利用所述步骤4所得地形起伏度值根据地形坡度公式计算高差阈值，具体包括：利用所述地形起伏度值根据地形坡度公式值进行计算，得到地形坡度，再利用所述地形坡度进行计算，得到高差阈值。

优选的，将所述差值与所述高差阈值进行比较，根据比较结果判断LiDAR点云数据是否为地面点，具体包括：当所述差值小于或者等于所述高差阈值，判断所述LiDAR点云数据为地面点，当所述差值大于所述高差阈值，判断所述LiDAR点云数据为地物点。

与现有技术相比，本发明的有益效果：

本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，通过对机载激光雷达扫描山区得到的LiDAR点云数据，然后通过构建虚拟格网对LiDAR点云数据进行数据组织并进行去噪处理，再通过确定点云数据的格网最大行号计算得到结构元素，基于结构元素对LiDAR点云数据进行形态学开运算得到形态学开运算前后的LiDAR点云数据的高差变化值，从而能够初步判断地物点和地面点，在对LiDAR点云数据进行形态学开运算的同时对LiDAR点云数据进行曲面拟合，引入地形起伏度，通过地形起伏度计算得到高差阈值，通过将LiDAR点云数据进行形态学开运算前后的高差变化值与进行曲面拟合得到的地形起伏度做差，将差值与高差阈值作比较，从而能够再次对地物点和地面点进行判断，减少误判，且利用形态学开运算与曲面拟合两者相结合对山区的LiDAR点云数据进行提取，能够对地面凸起等进行有效保护，更加准确的得到地面点，从而得到山区的真实地形；使得工作人员能够通过保留的点云数据解算出被覆盖在茂密植被和复杂地形环境中的危岩落石，为工程施工提供安全保障。

附图说明：

图1是本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法的流程示意图。

图2是本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法中对LiDAR点云数据进行形态学开运算的初始模拟地形图。

图3是本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法中对LiDAR点云数据进行形态学开运算的腐蚀运算结果模拟图。

图4是本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法中对LiDAR点云数据进行形态学开运算的膨胀运算结果模拟图。

图5是本发明所述的一种适用于山区危岩落石的LiDAR点云数据提取方法中对LiDAR点云数据进行形态学开运算前后的高差变化值与设定的高差阈值进行比较的示意图。

具体实施方式

下面结合试验例及具体实施方式对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。

实施例1

如图1所示，一种适用于山区危岩落石的LiDAR点云数据提取方法，包括以下步骤：

步骤1：利用机载激光雷达对山区进行扫描获取LiDAR点云数据；

步骤2：设定格网尺寸，对所述LiDAR点云数据构建虚拟格网；

步骤3：对所述LiDAR点云数据进行去噪处理，再对去噪处理后的所述LiDAR点云数据在所述虚拟格网中创建索引号，确定所述LiDAR点云数据在所述虚拟格网中的最大行号；

步骤4：根据所述最大行号计算结构元素，根据所述结构元素对所述LiDAR点云数据进行形态学开运算，得到形态学开运算前后的所述LiDAR点云数据的高差变化值；同时，根据所述结构元素对所述LiDAR点云数据进行曲面拟合，得到所述LiDAR点云数据的地面数字模型，根据所述地面数字模型得到所述LiDAR点云数据的地形起伏度值；

步骤5：利用所述步骤4所得地形起伏度值根据地形坡度公式计算高差阈值，并计算所述LiDAR点云数据的所述高差变化值与所述地形起伏度值的差，得到差值，将所述差值与所述高差阈值进行比较，根据比较结果判断LiDAR点云数据是否为地面点，若是，则提取所述LiDAR点云数据；若否，则剔除所述LiDAR点云数据；

步骤6：将所述结构元素与预设的结构元素阈值进行比较，若所述结构元素小于所述结构元素阈值，则输出提取到的LiDAR点云数据集；若否，返回步骤2更新设定的格网尺寸，重复步骤2至步骤5，以更新所提取的地面点。

试验中，虚拟格网是以规则格网为基础而建立的，而规则格网是最简单且易于实现的点云组织方法，其实质是采用一定的插值算法将离散点云重采样为规则排列的格网点云。这种数据组织方法具有三个优势：首先，这种方式能够简化数据结构、提高计算效率。其次，规则格网能够有效的降低点云数据的密度，从而减小参与计算的数据量，降低计算机占用内存。最后，这种组织方式还可以利用现有的图像运算方法，丰富点云数据的处理手段。但规则格网存在以下问题：第一，对数据缺失区域通常采用内插的方式进行数据填充，容易引入内插误差；第二，对点云进行重采样会破坏原有的数据结构，不可避免的会带来精度损失，影响后续处理精度。

虚拟格网是对点云在XY平面进行格网剖分但不进行内插的数据组织形式。任何一个点

在格网中的位置可通过以下公式计算获得：

式中，X_k、Y_k是某点云P的X、Y坐标，X_min是点云X方向的最小值，Y_min是点云Y方向的最小值，cellsize是划分格网的边长，即格网尺寸，X_i和Y_i表示点P在格网中的位置。

试验中，通过经验模态分解(Empirical Mode Decomposition，EMD)的方法对LiDAR点云数据进行去噪处理。该方法最初是为了精确描述信号频率随时间的变化而提出的直观且自适应好的瞬时频率分析方法。其基于信号本身，将信号分解为一系列数据序列，而这些数据序列通常具有不同的尺度特征。每个数据序列称之为本征模态函数(IntrinsicMode Function，简称IMF)。分解的结果由若干个本征模态分量imf_i(t)以及1个残余分量r_n(t)组成，如以下公式所示：

每个本征模态分量都是平稳的窄带信号，低阶分量与高频信号相对应，高阶分量与低频信号相对应。每个本征模态分量需符合两点：首先，对于全局数据，相较于过零点的个数，极值点的个数与其相等或者最多相差一个；其次，对于数据的每一点，上下包络的算数平均值为零。

首先将点云(x，y，z)变换为二维格网数据p(x，y)，并保留映射关系R；接着对p(x，y)采用EMD分解，每次分解时采用形态学膨胀运算和腐蚀运算获取数据的上包络u(x，y)和下包络l(x，y)，分解得到n个本征模态分量和1个残余分量；然后采用最大类间方差法的思想确定噪声主导模态的分界点k，并对前k个低阶本征模态分量采用软阈值约束，获取去噪后的各分量imf₁′，...，imf_k′；用所有的本征模态分量和残余分量进行重构，获得数据p’(x，y)，并根据上述的映射关系R，获得重构后的点云数据(x，y，z’)，最后对点云数据中的各个点按以下公式判断：S＝{p|abs(z(p)-z′(p))＞th1&&Num(abs(z(p)-z(p_i))＜th2)＜th3}

式中，S是噪声点集，p是待判定点，z(p)是p点原有高程值，z’(p)是p点重构后新的高程值，th1是噪声点高差阈值，一般取值为5m，abs(z(p)-z′(p))＞th1表示待判定点重构前后高程变化量要大于噪声点高差阈值。p_i表示点云中的任意一点，th2为临近点高差阈值，一般取值为1m，Num()表示满足条件点的个数，th3为临近点个数阈值，一般取值为3，Num(abs(z(p)-z(p_i))＜th2)＜th3表示待判定点一般为孤点，其临近点个数应小于临近点个数阈值。若满足上式条件，则将该p点判定为噪声点并进行剔除。

对噪声点进行剔除后，为方便检索，每个格网只存储格网内点云的索引号，而不存储坐标数据。点云的索引号ID是根据点在原点云中的位置而定义的。例如，点云P可以看作是N×Dim的向量，N是点云的个数，Dim表示点云的维数。

若为第k个点，那么该点的ID就为k。每个格网Gird(i，j)所存储点云ID的集合cell(i，j)可按下式进行定义：

式中，ID_k为点

的索引号，X_max是点云X方向的最大值，Y_max是点云Y方向的最大值，W为格网的最大行号，H为格网的最大列号，ceil(n)为不小于n的最小整数。

值得注意的是，如果虚拟格网内点云的个数为零，那么该虚拟格网的ID存储集合为空，这也是与规则格网不同的地方之一。虚拟格网既避免了规则格网因内插而引起的系统误差，同时还保留着规则格网组织结构简单、易于实现的特点，同时现有的图像处理方法同样也适用于虚拟格网组织的数据。如今，基于虚拟格网组织的点云数据处理方法已经有很多，如点云滤波、点云分割等。

利用形态学开运算的方法来实现地面点云的提取，一般的，形态学开运算由膨胀和腐蚀两种基本算法组成，对于机载LiDAR点云数据来说，即是对点云数据进行滤波，膨胀运算是取滤波窗口内(滤波结构元素为B)点云高程的最大值作为该点新的高程值，公式定义如下：

[δ_B(f)(x，y)＝max{f(x+i，y+i)|i，j∈[-W，W]}；(x+i)，(y+i)∈D_f]，

腐蚀运算则是取滤波窗口内点云高程的最小值作为该点新的高程值，公式定义如下：

[δ_B(f)(x，y)＝min{f(x+i，y+i)|i，j∈[-W，W]}；(x+i)，(y+i)∈D_f]，

式中，B是结构元素，大小为(2×W+1)×(2×W+1)，D_f是f的取值范围。

试验中，对LiDAR点云数据进行形态学开运算，先进行腐蚀运算，后进行膨胀运算，如图2所示，为一初始模拟地形图，包含有起伏的地面、以及小的危岩落石。首先对该区域地形进行腐蚀运算，腐蚀运算即是取邻域内点云高程的最小值作为一点新的高程值，而邻域内的最低点通常都为地面点，所以区域内的危岩体以及危岩落石会被剔除，然而部分地面点因地形起伏也会被剔除掉，腐蚀运算结果如图3所示。接着对腐蚀运算后的结果进行膨胀运算，可以发现除了危岩体和危岩落石之外，被腐蚀掉的地面点得到了有效的还原，如图4所不。

按照公式：dH＝f-δ_B[ε_B(f)]，用原有的点云高程值与形态学开运算后的高程值作差，得到形态学开运算前后点云数据的高差变化值dH，设定最大高差变化值dH_max为高差阈值T，将高差变化值dH大于高差阈值T的点(山体)判定为地物点并进行剔除，将高差变化值dH小于高差阈值T的点(危岩落石)判定为地面点并进行保留，如图5所示。

试验中，采用克里金插值的方法对LiDAR点云数据进行曲面拟合。

克里金插值是一种利用区域已知变量数据以及变异函数的结构特点，来对待定点进行无偏、最优估计的一种方法。设待定点p四周有n个已知采样点p₁，p₂，...，p_n，对应值为f(p_i)(i＝1，2，...，n)，则普通克里金插值法对该待定点p的估值为：

式中，λ_i为权重，表示没个样本值对待定点的重要程度。

由以上公式可以看出只要知道λ_i，便可求出未知点的估值。由克里金插值的定义可知，该方法应满足无偏性以及方差最小性，即：

采用拉格朗日乘常数法可以推出用变异函数表示的普通克里金方程组：

上式也可以用矩阵形式来表达：

Kλ＝D

λ＝K^-1D

从以上公式可以看出，要求出权重系数λ，需首先求得变异函数γ(h)。

本实施例选用球状模型作为理论变差函数模型，该模型通常表示为：

式中，c₀为块金值，c为拱高，h为滞后距，a为变程。

球状模型的各个系数可通过对已知样本点进行多项式回归求得。不同结构元素对应不同层级，已知样本点即为在当前层级滤波结构元素为B条件下点云数据中高程最小值对应的点，也可称之为控制点。由于在利用形态学开运算对LiDAR点云数据进行处理时，滤波结构元素的尺寸是与格网最大行号相关联，而格网最大行号与格网尺寸相联系，当通过设定的格网尺寸使格网最大行号从大到小变化时，那么结构元素便也从大到小变化，且最大滤波结构元素大于该区域最大危岩体的尺寸，因此在当前层级滤波结构元素为B条件下获取的控制点可以确定为地面点，而基于此类控制点进行克里金插值拟合出的曲面可以确定为当前层级对应的数字地面模型(digital terrain model，简称DTM)。获取DTM后，地形起伏度

便可通过公式

计算得到。

利用地形起伏度值根据地形坡度公式计算高差阈值

根据地形起伏度值计算得到高差阈值T，计算当前层级点云数据在形态学开运算前后高差变化值dH与当前层级点云数据所在虚拟格网的地面起伏度值

的差值，

若dH与

之间的差值大于该层级高差阈值，则将该点判定为地物点并进行剔除，否则，判定为地面点并进行保留。

预设结构元素阈值的大小为3×3，当前层级滤波完成后，比较当前结构元素与预设结构元素阈值，若当前结构元素小于预设结构元素阈值，则结束，否则，返回步骤2重新设定格网尺寸，从而计算得到新的结构元素，再次对点云数据进行处理，直到结构元素小于预设结构元素阈值结束处理，这样便能准确提取出地面点，从而根据保留的点云数据得到山区的真实地形。

其中，首先确定该点所在的格网，进而获得该格网对应的地面起伏度并按照以下的地形坡度公式计算出该层级对应的地形坡度，进而计算出高差阈值T：

式中，s(i，j)表示格网(i，j)处的地形坡度，i是行号，j是列号。s_c为地形增益，通常设置为常量。z_i，j+1、z_i，j-1、z_i+1，j和z_i-1，j分别为各个格网在拟合DTM中所对应的高程值；根据地形坡度计算高差阈值：

式中，B_k为第k次迭代结构元素的尺寸大小，cellsize为格网尺寸，dH₀为最小高差变化值，s为地形坡度，dH_max为最大高差变化值。

以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但本发明不局限于上述具体实施方式，因此任何对本发明进行修改或等同替换；而一切不脱离发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，包括以下步骤：

步骤1：利用机载激光雷达对山区进行扫描获取LiDAR点云数据；

步骤2：设定格网尺寸，对所述LiDAR点云数据构建虚拟格网；

步骤3：对所述LiDAR点云数据进行去噪处理，再对去噪处理后的所述LiDAR点云数据在所述虚拟格网中创建索引号，确定所述LiDAR点云数据在所述虚拟格网中的最大行号；

步骤4：根据所述最大行号计算结构元素，根据所述结构元素对所述LiDAR点云数据进行形态学开运算，得到形态学开运算前后的所述LiDAR点云数据的高差变化值；同时，根据所述结构元素对所述LiDAR点云数据进行曲面拟合，得到所述LiDAR点云数据的地面数字模型，根据所述地面数字模型得到所述LiDAR点云数据的地形起伏度值；

步骤5：利用所述步骤4所得地形起伏度值根据地形坡度公式计算高差阈值，并计算所述LiDAR点云数据的所述高差变化值与所述地形起伏度值的差，得到差值，将所述差值与所述高差阈值进行比较，根据比较结果判断LiDAR点云数据是否为地面点，若是，则提取所述LiDAR点云数据；若否，则剔除所述LiDAR点云数据；

步骤6：将所述结构元素与预设的结构元素阈值进行比较，若所述结构元素小于所述结构元素阈值，则输出提取到的LiDAR点云数据集；若否，返回步骤2更新设定的格网尺寸，重复步骤2至步骤5，以更新所提取的地面点。

2.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，所述虚拟格网是以规则格网为基础进行建立，是对所述LiDAR点云数据在XY平面进行格网剖分但不进行内插的数据组织形式。

3.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，所述步骤3中，采用经验模态分解的方法对所述LiDAR点云数据进行去噪处理。

4.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，所述步骤4中，对所述LiDAR点云数据进行形态学开运算，具体包括：基于所述结构元素对所述LiDAR点云数据先进行腐蚀运算，再进行膨胀运算。

5.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，所述步骤4中，采用克里金插值的方法对所述LiDAR点云数据进行曲面拟合。

6.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，所述步骤5中，利用所述步骤4所得地形起伏度值根据地形坡度公式计算高差阈值，具体包括：利用所述地形起伏度值根据地形坡度公式值进行计算，得到地形坡度，再利用所述地形坡度进行计算，得到高差阈值。

7.根据权利要求1所述的一种适用于山区危岩落石的LiDAR点云数据提取方法，其特征在于，将所述差值与所述高差阈值进行比较，根据比较结果判断LiDAR点云数据是否为地面点，具体包括：当所述差值小于或者等于所述高差阈值，判断所述LiDAR点云数据为地面点，当所述差值大于所述高差阈值，判断所述LiDAR点云数据为地物点。

Images