CN109697754A - 基于主方向估算的3d岩体点云特征面提取方法 - Google Patents

Abstract

本发明涉及一种基于主方向估算的3D岩体点云特征面提取方法。本发明主要包含3个阶段：(1)基于体素的聚类；(2)基于高斯核的主方向估算；(3)特征面提取。首先，利用空间二级网格将点云快速体素化并聚类成三种体素：共面体素、不共面体素和稀疏体素；然后对共面体素，使用二维高斯核自动估算点云特征面的主方向；最后根据主方向寻找种子体素，使用基于体素的区域增长，实现特征面的自动提取，得到多组特征面，并且每组特征面内的子特征面共面或平行。本发明适用于大场景岩石激光扫描点云的特征平面提取，可以高效、准确的提取岩体点云中的特征面。

Description

基于主方向估算的3D岩体点云特征面提取方法

技术领域

本发明属于计算机视觉和三维点云数据处理的交叉领域，涉及点云特征提取、岩体三维重建等技术，特别涉及基于主方向估算的3D岩体点云特征面提取方法。

背景技术

点云中的特征面提取是许多计算机图形学、图像处理和计算机视觉中的重要内容，包括三维重建、物体识别、虚拟现实等。LiDAR(Light Detection And Ranging)扫描器能够以高分辨率、高准确度，从很远的距离(达到6000m)，以超高速(每秒超过222000次)获取陆地表面的XYZ坐标等信息，LiDAR的高速发展与普及，进一步促进了对于不规则点云的特征面检测、提取和建模的研究。

三维点云特征面提取最常用的方法是霍夫变换(Hough Transform,HT)、随机采样一致性(Random Sample Consensus，RANSAC)和区域生长(Region growing,RG)。霍夫变换是2D图像中检测直线、圆和其他参数化形状的方法，也可以被用于检测三维点云中的3D物体(特征面、圆柱体等)，但传统的霍夫变换计算复杂、效率低、特征面可能不连续；RANSAC通过迭代的方式随机选取三个点，计算三点定义的特征面得分，当得分达到稳定时算法结束，但RANSAC可能检测到错误的特征面、特征面可能不连续；区域生长是对局部进行搜索，识别和扩充具有相同特征的区域，但可能导致空洞和过分割等问题。

目前，点云中的特征面提取多是对城市建筑物点云的研究，这类3D激光点云大多包含规则的几何特征面(如建筑物墙壁面、地面和桌面等)，特征面较光滑，而岩体点云表面不平坦、粗糙，形状不规则，多是以高陡边坡的形式呈现，与常见城市点云差异较大。在堤坝、地下隧道、地基和边坡等岩体工程中，岩体的稳定性分析是一项极其重要的内容，目前，岩体稳定性分析的多种方法(如块体理论、非连续变性分析、数值流形法等)都是作用于精准的岩体数值模型之上，对岩体点云的特征面提取及三维建模成为岩体工程的重要部分。

发明内容

本发明针对岩体点云表面粗糙、形状不规则等特性，为了克服现有技术在岩体点云特征面提取中效率低、准确率低、存在过分割和欠分割等问题，提出了一种基于主方向估算的3D岩体点云特征面提取方法，能搞实现高效提取特征面，且具有较高的准确率。

本发明提供一种基于主方向估算的3D岩体点云特征面提取方法，主要包涵以下步骤：

步骤1、基于体素的聚类：利用空间二级网格将点云快速体素化，并根据体素内最小样本点数量阈值和共面判定准则，将体素分类为共面体素、非共面体素和稀疏体素3类；

步骤2、基于高斯核的主方向估算：对共面体素，使用二维高斯核自动估算点云特征面的主方向；

步骤3、特征面提取：根据主方向寻找种子体素，使用基于体素的区域增长来提取特征面。

本发明所述的一种基于主方向估算的3D岩体点云特征面提取方法，作为优选方式，步骤1进一步包含以下步骤：

步骤1.1、根据给定体素长度阈值，将点云快速的体素化；

步骤1.2、根据给定的体素内最小样本点数量阈值，将体素划分为稀疏体素和非稀疏体素，对非稀疏体素内的样本点进行主成分分析(PCA)，得到样本点的特征值λ₁，λ₂，λ₃和均方误差(Mean Squared Error,MSE)，其中λ₁≤λ₂≤λ₃，根据共面判定准则λ₂＞S_αλ₁&&MSE＜ε，将非稀疏体素划分为共面体素和非共面体素；

步骤1.3、将步骤1.2得到的非共面体素继续8等分，对每等分重复步骤1.2。

本发明所述的一种基于主方向估算的3D岩体点云特征面提取方法，作为优选方式，步骤2进一步包含以下步骤：

步骤2.1、投票所使用的累加器为半球面型累加器，半球面上每个累加器单元的面积相等，单位法向量落在每个累加器单元上的概率相等，由于法向量的方向可正可反，对于落在下半球面的法向量将之方向取反；

步骤2.2、对共面体素K计算协方差矩阵∑(x,y,z)、雅可比式J和协方差矩阵∑(φ,θ)，具体而言，

其中θ∈[0°，360°)，φ∈[0°，180°]

对协方差矩阵∑(φ,θ)进行特征分解，最小的特征值为λ_{φθ_min}，λ_{φθ_min}对应的特征向量为V_{φθ_min}，标准差std_dev＝sqrt(λ_{φθ_min})，计算投票阈值g_min，具体而言，

g_min＝Gaussian(2×std_dev×V_{φθ_min})

其中Gaussian()为二维高斯分布概率密度方程，其具体形式如下：

其中，其中μ表示(φ_k,θ_k)，表示距离中心点μ的距离，∑表示关于(φ,θ)的协方差矩阵；

步骤2.3、对共面体素K，考虑其体素尺寸和包含样本点的数量，计算其投票权重其中voxel_size表示共面体素K的体素长度，mesh_size表示为输入点云包围盒的空间尺寸，|K|表示共面体素K包含的样本数，|P|表示输入点云总的样本数，w_a,w_d是权重系数且满足w_a+w_d＝1.0。对共面体素K，在其对应的累加器单元(φ_k,θ_k)及其周围邻域计算p值，若p＞g_min，则进行投票，投票值Vote＝p×w_k；

步骤2.4、对每个共面体素重复上述步骤2.2至步骤2.3，最后使用滑动窗技术选择累加器的局部最大值，即主方向。

本发明所述的一种基于主方向估算的3D岩体点云特征面提取方法，作为优选方式，步骤3进一步包含以下步骤：

步骤3.1、利用主方向，从共面体素中选取种子体素并存入队列Q中；

步骤3.2、当Q不为空时，取出队列Q的元素q，并搜索q的邻域体素，若q的邻域体素为共面体素，则进行基于体素的区域生长，若q的邻域体素为稀疏体素或非共面体素时，对体素内的样本点进行基于点的区域生长。

本发明针对岩体点云特性的分析，利用体素化和基于二维高斯核的主方向估算，大大提高了特征面提取的效率；同时，综合利用基于体素的区域生长和基于点的区域生长进行特征面提取，提高了特征面提取的准确性。

附图说明

图1是基于主方向估算的3D岩体点云特征面提取方法流程图。

图2是共面体素示意图。

图3是使半球面型累加器投票结果示意图。

图4是特征面提取结果示意图。

具体实施方式

本发明主要基于计算机视觉和点云处理技术，提出一种基于主方向估算的3D岩体点云特征面提取方法。本发明综合利用基于二维高斯核主方向估算和基于体素的区域生长进行特征面提取，提高了特征面提取的效率和准确性，扩展了现有特征面提取方法的适用范围。

为了使本发明的目的、技术方案及优点更加清楚，以下结合附图和实施例，对本发明作进一步详细说明。本发明提供的方法能够用计算机软件技术实现流程，整体技术流程图参见图1，包括以下步骤：

步骤1、为实现将点云快速体素化并聚类成三种体素：共面体素、不共面体素和稀疏体素，首先设定体素的空间长度和体素内最少样本点的数量。实施例的具体实施过程说明如下：

步骤1.1、根据给定体素长度阈值voxel_size(本实施例中取值1.73m)，将点云体素化；

步骤1.2、根据给定的体素内最小样本点数量阈值(本实施例中取值300)，将体素划分为稀疏体素和非稀疏体素；

步骤1.3、对非稀疏体素内的样本点进行主成分分析(PCA)，得到样本点的特征值λ₁，λ₂，λ₃和均方误差(Mean Squared Error,MSE)，其中λ₁≤λ₂≤λ₃；

步骤1.4、根据共面判定准则λ₂＞S_αλ₁&&MSE＜ε(本实施例中取值S_α取值30，ε取值0.05)，将非稀疏体素划分为共面体素和非共面体素；

步骤1.5、将步骤1.4得到的非共面体素8等分，对每等分重复步骤1.2至步骤1.4。

利用空间二级网格将点云快速体素化，并根据体素内最小样本点数量阈值和共面判定准则，将体素分类为共面体素(参见图2)、非共面体素和稀疏体素3类并用于步骤2和步骤3中。

步骤2、根据步骤1得到的共面体素(参见图2)，使用二维高斯核自动估算点云特征面的主方向，实施例的具体实施过程说明如下：

步骤2.1、使用半球面型累加器，对累加器进行离散化；

步骤2.2、对共面体素K计算协方差矩阵∑(x,y,z)、雅可比式J和协方差矩阵∑(φ,θ)，具体而言，

其中θ∈[0°，360°)，φ∈[0°，180°]

步骤2.3、对协方差矩阵∑(φ,θ)进行特征分解，将最小的特征值赋值给λ_{φθ_min}，将最小特征值对应的特征向量赋值给V_{φθ_min}，将sqrt(λ_{φθ_min})赋值给标准差std_dev，计算投票阈值g_min，具体而言，

g_min＝Gaussian(2×std_dev×V_{φθ_min})\*MERGEFORMAT (3)

其中Gaussian()为二维高斯分布概率密度方程，其具体形式如下：

其中，其中μ表示(φ_k,θ_k)，表示距离中心点μ的距离，∑表示关于(φ,θ)的协方差矩阵；

步骤2.4、对共面体素K，计算其投票权重w_k，其具体形式如下：

其中，voxel_size表示共面体素K的体素长度，mesh_size表示为输入点云包围盒的空间尺寸，|K|表示共面体素K包含的样本数，|P|表示输入点云总的样本数，w_a,w_d是权重系数且满足w_a+w_d＝1.0(本实施例中w_a取值0.75，w_d取值0.25)；

步骤2.5、对共面体素K，在其对应的累加器单元(φ_k,θ_k)及其周围邻域，使用公式(4)计算p值，若p＞g_min，则进行投票，投票值Vote＝p×w_k，参见图3；

步骤2.6、对每个共面体素重复上述步骤2.2至步骤2.5；

步骤2.7、使用N×N(本实施例中取值8)的滑动窗选择累加器的局部最大值(即主方向)，并按投票值降序存入队列S中。

利用二维高斯核和共面体素估算点云的主方向，极大地提高了计算效率，并提取出点云的主方向组，用于步骤3。

步骤3、利用步骤2得到的主方向，从共面体素中选取种子体素，并结合基于体素的区域生长和基于点的区域生长进行特征面提取，实施例的具体实施过程说明如下：

步骤3.1、初始化所有体素为“未遍历”，初始化队列Q为空；

步骤3.2、利用步骤2得到的主方向组，弹出队列S的队首元素s，计算其法向量

步骤3.3、在共面体素中搜索，选取体素法向量和的夹角小于阈值t₁(本实施例中取值5°)的所有体素，从中选取夹角最小的体素作为种子体素k，将k存入列Q并标记为“已遍历”；

步骤3.4、当Q不为空时，弹出队列Q的队首元素q，搜索q的“未遍历”邻域体素n；

步骤3.5、若体素n是共面体素，如果同时满足条件：①邻域体素n的法向量和种子区域的法向量夹角小于角度阈值t₂(实施例中取值15°)②邻域体素n的重心到种子体素区域的垂直距离小于阈值Max_dis(本实施列中取值0.3m)，则将邻域体素n与种子体素合并为同一区域，更新种子区域的法向量并把体素n压入队列Q的队尾，标记体素n标记为“已遍历”。

步骤3.6、若体素n是非共面体素或稀疏体素，则对体素内的每个样本点进行判定，如果同时满足条件：①样本点m的法向量和种子区域的法向量夹角小于角度阈值t₂(实施例中取值15°)②样本点m到种子体素区域的垂直距离小于阈值Max_dis(本实施列中取值0.3m)，则将样本点m与种子体素合并为同一区域，若体素n内超过50％的样本点与种子体素合并为同一区域，则更新种子区域的法向量把体素n压入队列Q的队尾，标记体素n标记为“已遍历”。

步骤3.7、如果队列Q不为空，则返回步骤3.4，如果队列Q为空，一个特征面提取结束，返回步骤3.3；

步骤3.8、重复步骤3.2至步骤3.7，直到队列S为空，得到特征面提取结果，参见图4.

该技术发明从运算效率和特征面提取准确率方面比较传统特征面提取方法有极大提高，表1给出了对多组点云进行特征面提取的时间消耗，在特征面提取精度方面，经过定量分析，3个评估度量Precision,Recall和F1得分分别达到了91.92％、91.67％和91.8％。

表1特征面提取时间消耗

本文以上所述仅的具体实施例仅仅是对本发明精神做举例说明，应当指出：本领域技术人员在不脱离发明原理的前提下可以对所描述的具体实施例做出任何修改、变化或等效，但这些修改、变化和等效都将落入本发明的保护范围之内。

Claims (4)

1.一种基于主方向估算的3D岩体点云特征面提取方法，其特征在于主要包涵以下步骤：

步骤1、基于体素的聚类：利用空间二级网格将点云快速体素化，并根据体素内最小样本点数量阈值和共面判定准则，将体素分类为共面体素、非共面体素和稀疏体素3类；

步骤2、基于高斯核的主方向估算：对共面体素，使用二维高斯核自动估算点云特征面的主方向；

步骤3、特征面提取：根据主方向寻找种子体素，使用基于体素的区域增长来提取特征面。

2.根据权利要求1所述的一种基于主方向估算的3D岩体点云特征面提取方法，其特征在于：所述步骤1进一步包括以下步骤：

步骤1.1、根据给定体素长度阈值，将点云快速的体素化；

步骤1.2、根据给定的体素内最小样本点数量阈值，将体素划分为稀疏体素和非稀疏体素，对非稀疏体素内的样本点进行主成分分析(PCA)，得到样本点的特征值λ₁，λ₂，λ₃和均方误差(Mean Squared Error,MSE)，其中λ₁≤λ₂≤λ₃，根据共面判定准则λ₂＞S_αλ₁&&MSE＜ε，将非稀疏体素划分为共面体素和非共面体素；

步骤1.3、将步骤1.2得到的非共面体素继续8等分，对每等分重复步骤1.2。

3.根据权利要求1所述的一种基于主方向估算的3D岩体点云特征面提取方法，其特征在于：所述步骤2进一步包括以下步骤：

步骤2.1、投票所使用的累加器为半球面型累加器，半球面上每个累加器单元的面积相等，单位法向量落在每个累加器单元上的概率相等，由于法向量的方向可正可反，对于落在下半球面的法向量将之方向取反；

步骤2.2、对共面体素K计算协方差矩阵∑(x,y,z)、雅可比式J和协方差矩阵∑(φ,θ)，具体而言，

其中θ∈[0°，360°)，φ∈[0°，180°]

对协方差矩阵∑(φ,θ)进行特征分解，最小的特征值为λ_{φθ_min}，λ_{φθ_min}对应的特征向量为V_{φθ_min}，标准差std_dev＝sqrt(λ_{φθ_min})，计算投票阈值g_min，具体而言，

g_min＝Gaussian(2×std_dev×V_{φθ_min})

其中Gaussian()为二维高斯分布概率密度方程，其具体形式如下：

其中，其中μ表示(φ_k,θ_k)，表示距离中心点μ的距离，∑表示关于(φ,θ)的协方差矩阵；

步骤2.3、对共面体素K，考虑其体素尺寸和包含样本点的数量，计算其投票权重其中voxel_size表示共面体素K的体素长度，mesh_size表示为输入点云包围盒的空间尺寸，|K|表示共面体素K包含的样本数，|P|表示输入点云总的样本数，w_a,w_d是权重系数且满足w_a+w_d＝1.0，对共面体素K，在其对应的累加器单元(φ_k,θ_k)及其周围邻域计算p值，若p＞g_min，则进行投票，投票值Vote＝p×w_k；

步骤2.4、对每个共面体素重复上述步骤2.2至步骤2.3，最后使用滑动窗技术选择累加器的局部最大值，即主方向。

4.根据权利要求1所述的一种基于主方向估算的3D岩体点云特征面提取方法，其特征在于：所述步骤3进一步包括以下步骤：

步骤3.1、利用主方向，从共面体素中选取种子体素并存入队列Q中；

步骤3.2、当Q不为空时，取出队列Q的元素q，并搜索q的邻域体素，若q的邻域体素为共面体素，则进行基于体素的区域生长，若q的邻域体素为稀疏体素或非共面体素时，对体素内的样本点进行基于点的区域生长。