CN115830587A - 一种基于高精度点云数据的结构面快速自动识别方法 - Google Patents

Abstract

本发明公开了一种基于高精度点云数据的结构面快速自动识别方法，就结构面分割时，现有区域生长速度过慢、受生长准则的阈值影响较大的问题，引入模糊界限的理念，以减少阈值对区域生长的影响；同时预分配了搜索空间加快了区域生长的分割速度，形成了结构面的快速准确的分割方法。该方法在保证准确识别结构面的前提下，快速提取结构面地质信息，并以自动算法的方式提高获取结构面地质信息的效率。

Description

一种基于高精度点云数据的结构面快速自动识别方法

技术领域

本发明属于测绘工程领域，具体涉及一种基于高精度点云数据的结构面快速自动识别方法。

背景技术

三维激光扫描技术是非接触的测量技术之一，具有获取数据速度快、精度高、全自动化、全天候作业等特点，已成为地理空间数据获取的重要技术手段。为解译海量三维空间数据中蕴含的地理信息，将数据经聚类和分割后，划分为可处理的子类是最主要的解译方法。

聚类分析是无监督分类方法，在不需要先验认知或少量先验认知的情况下完成数据集的分组。而地理信息，尤其是结构面信息的界限大部分是基于工程经验的，模糊的；而模糊聚类通过概率的形式描述元素的相似度，因此更符合实际工程和经验认知。同时，模糊聚类也使用比较灵活，受噪声点影响较小因而得到了广泛运用。

区域生长是一种分割算法，具有原理简单、算法易于实现、分割细致的优点。但仍存在以下问题：①传统的区域生长法的可能生长对象是数据集的所有元素，需要逐一判断元素的生长可能性。其目的在于保证算法的准确性，但当数据集过大且大部分元素均不是可能的生长元素时，算法花费大量时间处理了不必要的数据，这会导致生长速度过慢，同时空耗大量的运行时间和存储空间。②此外，传统的区域生长准则的阈值为确定值，生长准则需要将各区域绝对区分开才能保证结果的准确性，这样的生长准则难以寻找，同时确定的生长准则导致各区域界限分明，不适合界限相对模糊的结构面分割。

在结构面识别方面已经做出了许多成果：Gigli和Casagli对点云进行空间体素分割，设计了半自动的结构面识别方法，总结了结构面地质信息的提取方法；(Gigli andCasagli,Semi-automatic extraction of rock mass structural data from highresolution LIDAR point clouds[J].International Journal of Rock Mechanics andMining Sciences.2011(48),187-198.)而Riquelme使用k-means(k均值聚类)的方法，实现了结构面自动识别方法。(Riquelme等,Identification of Rock Slope DiscontinuitySets from Laser Scanner and Photogrammetric Point Clouds:A ComparativeAnalysis[J].Procedia Engineering.2017(191),838-845.)。此外，董秀军通过点云进行HSV颜色空间着色，在进行聚类也实现了结构面的自动识别。(一种基于三维激光扫描的岩体结构面识别与信息提取方法[P].董秀军,等.中国专利:CN112529844A,2021-03-19.)

但是点云数据日益庞大复杂，而许多工程又需要及时获取结构面信息，因此对结构面地质信息的快速和准确获取提出了更高要求。国内外对结构面地质信息获取进行了初步的研究，但是对结构面地质信息的准确性和提取速度需要进一步提高。

发明内容

本发明所为了解决背景技术中存在的技术问题，目的在于提供了一种基于高精度点云数据的结构面快速自动识别方法，基于高精度点云的结构面识别的新方法，利用模糊聚类的结构面三维模型结果，优化了区域生长实现结构面的快速分割。在保证分割准确性的前提下，有效的提高了结构面识别的准确性，并且该方法能够在数分钟对百万数据的点云完成结构面识别，提高了地质测量的效率。

为了解决技术问题，本发明的技术方案是：

一种基于高精度点云数据的结构面快速自动识别方法，所述方法包括：

对预处理后的点云数据设置相匹配的三角网模型的三角形尺寸后，利用三角剖分算法建立三角网模型；

基于三角网模型计算三角网的法向量，并将三角网的法向量归一化至单位球面空间内，得到三角网的法向量的归一化结果；

对三角网模型内所有三角形排序进行排序，得到排序结果；

基于三角网的法向量的归一化结果，估算初始聚类中心V，迭代完成模糊聚类，得到模糊聚类结果；

基于排序结果和模糊聚类结果，对预判的生长对象进行预筛选，并更换区域生长种子，直到完成结构面分割，得到所有结构面的分离结果。

进一步，所述对预处理后的点云数据处理之前，所述方法还包括：

通过布置多个站点对测量区域进行多角度的高精度扫描，获取多个站点的高精度点云数据。

进一步，对多个站点的高精度点云数据进行点云拼接和降噪处理，得到预处理后的点云数据。

进一步，基于结构面识别精度的预设值，设置相应的三角网模型的三角形尺寸，再利用Delaunay三角剖分算法建立三角网模型，所述模型的N个三角面片以三角形顶点和顶点索引的形式存储；三角形顶点即点云的坐标(x，y，z)，与点云的点数量相同；顶点索引即(a_i，b_i，c_i)三个顶点以顺时针排列。

进一步，所述计算三角网的法向量，具体包括：

根据三角形顶点和顶点索引，遍历计算顶点索引为i(a_i，b_i，c_i)的三角面片的外法向向量F_i，并将其归一化至单位球面空间内；计算方法如下：

进一步，对三角网模型内所有三角形排序进行排序，具体包括：

基于三角形中心到原点的距离R对三角网内所有三角形排序，若多个三角形位于同一球面内，按高程对这些三角形排序；所述排序完成后相邻的三角形在一定厚度的球壳内，即：

sort(R)→R₁≤R₂…≤R_n。

进一步，所述估算初始聚类中心V，迭代完成模糊聚类，具体包括：

若设置了初始聚类中心和聚类数K，采用余弦距离计算球面空间数据的相异指标d，则第i个三角面片相对于第j个中心的相异度d_ij为：

d_ij ²＝1-(F_i·V_j)²，F_i●V_j为向量F_i与V_j的矢量积；

计算模糊聚类的隶属函数矩阵U，对于任意三角面片i，它对第j个聚类中心的隶属度u_ij为：

其中聚类数没有设定时将从2开始遍历；m为模糊度m＝2时能满足大部分聚类要求，为程序默认值，可以根据实际情况更改m；

若没有设置初始聚类中心和聚类数K，将采用随机算法生成隶属函数矩阵U，再计算出聚类中心V；第j个聚类中心的计算为：

然后计算此聚类中心下的隶属函数矩阵U；

通过上述计算的隶属函数矩阵U，计算新的聚类中心V*，通过新的聚类中心V*计算新的隶属函数矩阵U*，直到模糊聚类的目标函数J_m(U,V)函数为最小值时，聚类完成；模糊聚类的目标函数J_m(U,V)函数的计算为：

当两次迭代的隶属函数矩阵目标函数J_m，J_m ^*小于设置的容许误差ε时，判定算法收敛，算法终止迭代的判别式为：

J_m-J_m ^*＜ε，ε＝10^-5。

进一步，所述对预判的生长对象进行预筛选，具体包括：

当生长种子位于一小球面内(R_d1)，则预判的生长对象为大球面内(R_d2)和小球面外的球壳区域；球壳厚度由三角网最大三角形尺寸d_max决定，为保证所有预判的生长对象位于此区域内，取5倍最大三角形尺寸作为预筛选的球壳厚度R₀；

即：

R₀＝5*max(d)；

R_d2＝R_d1+R₀。

进一步，所述更换区域生长种子，直到完成结构面分割，具体包括：

I：选取一个三角面片作为生长种子；

Ii：搜索此球壳内与种子三角形共边的三角形，该三角形为生长可能需要的三角形；

Iii：定义生长原则，只有共边的三角形与种子三角形具有模糊聚类结果时才进行区域生长；

Iv：将生长的三角形作为种子三角形重复上述i-iii步骤直到生长停止(球壳内没有搜索到具有模糊聚类结果的共边的三角形)，完成一个结构面的分离；

V：更换新的生长种子重复上述i-iv步骤，直到所有三角形都完成分离，此时所有结构面完成分离。

进一步，在得到所有结构面的分离结果后，所述方法还包括：

计算各个结构面的地质信息：

I：产状信息：结构面法向向量n(x，y，z)为其三角面片法向量的平均值，倾向

和倾角θ的对应关系应根据法向向量n所处卦限不同分类处理，具体如下：

ii：结构面出露面积：结构面出露面积为属于它的所有三角面片面积之和。

Iii：结构面迹长：以结构面沿倾向或走向方向的最大值长度作为迹长，即将结构面出露的三角形沿倾向或走向方向投影，去投影长度的大值作为迹长。

与现有技术相比，本发明的优点在于：

本发明就结构面分割时，现有区域生长速度过慢、受生长准则的阈值影响较大的问题，引入模糊界限的理念，以减少阈值对区域生长的影响；同时预分配了搜索空间加快了区域生长的分割速度，提出了结构面的快速准确的分割方法。该方法在保证准确识别结构面的前提下，快速提取结构面地质信息，并采用自动算法的方式提高获取结构面地质信息的效率。

附图说明

图1、模糊聚类示意图；

图2：区域生长示意图；

图3：结构面测量示意图；

图4：选择测量区域图；

图5：三角网外法向量计算图；

图6：结构面识别结果图。

具体实施方式

下面结合实施例描述本发明具体实施方式：

需要说明的是，本说明书所示意的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容能涵盖的范围内。

同时，本说明书中所引用的如“上”、“下”、“左”、“右”、“中间”及“一”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

实施例1

(1)高精度点云获取。布置多个站点对测量区域进行多角度的高精度扫描，合理进行加密扫描以满足精度要求。

(2)点云拼接和降噪。对多个站点的数据进行拼接，同时删除由于外界干扰如仪器自身震颤，扬尘等作用产生噪点，以及不需要的区域。

(3)生成三角网模型。导入点云数据，根据结构面识别的精度要求，设置相应的三角网模型的三角形尺寸，再利用Delaunay三角剖分算法建立三角网模型，模型的N个三角面片以三角形顶点(点云的坐标，(x，y，z)，与点云的点数量相同)和顶点索引((a_i，b_i，c_i)三个顶点以顺时针排列)的形式存储。

(4)计算三角网的法向量。依据三角形顶点(点云的坐标)和顶点索引，计算三角面片i的外法向向量F_i，并将其归一化至单位球面空间内。计算方法如下：

(5)对三角网模型排序。按三角形中心到原点的距离R对三角网内所有三角形排序，若多个三角形位于同一球面内，按高程对这些三角形排序；排序完成后相邻的三角形总是在一定厚度的球壳内，具体见附图1。即：

sort(R)→R₁≤R₂…≤R_n

(6)设置初始聚类中心V和聚类数K。如果能根据相关资料或者前期的地质调查得到测量区域的结构面组数和优势方向，通过先验知识设置初始聚类中心和聚类数能提高结构面识别的准确性，如果没有相关数据则进行下一步。

(7)估算初始聚类中心V，迭代完成模糊聚类(如图1所示)。其过程如下：

i如果设置了初始聚类中心，采用余弦距离计算球面空间数据的相异指标d，则第i个三角面片相对于第j个中心的相异度d_ij为：

d_ij ²＝1-(F_i·V_j)²，F_i●V_j为向量F_i与V_j的矢量积；

然后计算模糊聚类的隶属函数矩阵U，对于任意三角面片i，它对第j个聚类中心的隶属度u_ij为：

其中聚类数没有设定时将从2开始遍历；m为模糊度m＝2时能满足大部分聚类要求，为程序默认值。

ii如果没有设置初始聚类中心，将采用随机算法生成隶属函数矩阵U，再计算出聚类中心V。第j个聚类中心的计算为：

然后计算此聚类中心下的隶属函数矩阵U。

iii通过前一次计算的隶属函数矩阵U，计算新的聚类中心V*，通过新的聚类中心V*又能计算新的隶属函数矩阵U*。直到模糊聚类的目标函数J_m(U,V)函数为最小值时，聚类完成。模糊聚类的目标函数J_m(U,V)函数的计算为：

但为了防止过度迭代，当两次迭代的隶属函数矩阵目标函数J_m，J_m ^*小于设置的容许误差ε时，认为算法收敛。算法终止迭代的判别式为：

J_m-J_m ^*＜ε，ε＝10^-5。

(9)预筛选可能的生长对象。当生长种子位于某一小球面内(上)R_d1，则可能的生长对象为大球面内(上)R_d2和小球面外的球壳区域见附图2。球壳厚度由三角网最大三角形尺寸d_max决定，为保证所有可能的生长对象位于此区域内，取5倍最大三角形尺寸作为预筛选的球壳厚度R₀。即：

R₀＝5*max(d)；

R_d2＝R_d1+R₀；

(10)更换区域生长种子，直到完成结构面分割(如图2所示)。完整的步骤如下：

i get{R_d1}，{growth}+＝{R_d1}；

选取一个三角面片作为生长种子；

ii R_d2＝Rd1+R0，{growth0}＝(R_d1：R_d2)∩R_d1；

搜索此球壳内与种子三角形共边的三角形，该三角形为生长可能需要的三角形；

iii{growth}+＝(Id{growth0}＝＝Id(R_d1))；

定义生长原则，只有共边的三角形与种子三角形具有模糊聚类结果时才进行区域生长；

iv{R_d1}＝{growth}，go loop(i：iii)，

end，cluster＝{growth}，delete R{growth}；

将生长的三角形作为种子三角形重复上述步骤直到生长停止(球壳内没有搜索到具有模糊聚类结果的共边的三角形)，完成一个结构面的分离；

v go loop(i：iv)，

end；

更换新的生长种子，直到所有三角形都完成分离，此时所有结构面完成分离。

(11)计算个结构面的产状。结构面法向向量n(x，y，z)为其三角面片法向量的平均值，倾向

和倾角θ的对应关系应根据法向向量n所处卦限不同分类讨论。具体如下：

(12)结构面出露面积和迹长(如图3所示)。结构面出露面积为属于它的所有三角面片面积之和，以结构面沿倾向或走向方向的最大值长度作为迹长，即将结构面出露的三角形沿倾向或走向方向投影，去投影长度的大值作为迹长。

(13)输出结构面信息。

实施例2：

下面结合具体实例，说明本方法的具体实施方式。

(1)现场勘查，进行相关外业工作，做好相关记录；获取测量区域的点云；

(2)完成点云的拼接和降噪工作，选择要测量的区域(见附图4)；

(3)根据现场调查成果和测量精度要求设置三角网模型的三角形尺寸，生成并优化三角网模型(见附图5)；

(4)计算三角网外法向向量，(见附图5)；

(5)算法自动进行模糊聚类，区域生长，完成结构面识别。聚类结果见附图6；

(6)测量结构面信息；

(7)输出结构面测量数据。

本实施例中算法在i5-10400F CPU，8GB RAM的计算机上处理了约78万的三角网数据，总耗时5分钟。

附表:部分结构面信息表

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

上面对本发明优选实施方式作了详细说明，但是本发明不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

不脱离本发明的构思和范围可以做出许多其他改变和改型。应当理解，本发明不限于特定的实施方式，本发明的范围由所附权利要求限定。

Claims (10)

1.一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述方法包括：

对预处理后的点云数据设置相匹配的三角网模型的三角形尺寸后，利用三角剖分算法建立三角网模型；

基于三角网模型计算三角网的法向量，并将三角网的法向量归一化至单位球面空间内，得到三角网的法向量的归一化结果；

对三角网模型内所有三角形排序进行排序，得到排序结果；

基于三角网的法向量的归一化结果，估算初始聚类中心V，迭代完成模糊聚类，得到模糊聚类结果；

基于排序结果和模糊聚类结果，对预判的生长对象进行预筛选，并更换区域生长种子，直到完成结构面分割，得到所有结构面的分离结果。

2.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述对预处理后的点云数据处理之前，所述方法还包括：

通过布置多个站点对测量区域进行多角度的高精度扫描，获取多个站点的高精度点云数据。

3.根据权利要求2所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，对多个站点的高精度点云数据进行点云拼接和降噪处理，得到预处理后的点云数据。

4.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，基于结构面识别精度的预设值，设置相应的三角网模型的三角形尺寸，再利用Delaunay三角剖分算法建立三角网模型，所述模型的N个三角面片以三角形顶点和顶点索引的形式存储；三角形顶点即点云的坐标(x，y，z)，与点云的点数量相同；顶点索引即(a_i，b_i，c_i)三个顶点以顺时针排列。

5.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述计算三角网的法向量，具体包括：

根据三角形顶点和顶点索引，遍历计算顶点索引为i(a_i，b_i，c_i)的三角面片的外法向向量F_i，并将其归一化至单位球面空间内；计算方法如下：

6.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，对三角网模型内所有三角形排序进行排序，具体包括：

基于三角形中心到原点的距离R对三角网内所有三角形排序，若多个三角形位于同一球面内，按高程对这些三角形排序；所述排序完成后相邻的三角形在一定厚度的球壳内，即：

sort(R)→R₁≤R₂…≤R_n。

7.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述估算初始聚类中心V，迭代完成模糊聚类，具体包括：

若设置了初始聚类中心和聚类数K，采用余弦距离计算球面空间数据的相异指标d，则第i个三角面片相对于第j个中心的相异度d_ij为：

d_ij ²＝1-(F_i·V_j)²，F_i·V_j为向量F_i与V_j的矢量积；

计算模糊聚类的隶属函数矩阵U，对于任意三角面片i，它对第j个聚类中心的隶属度u_ij为：

其中聚类数没有设定时将从2开始遍历；m为模糊度，m＝2时能满足大部分聚类要求，为程序默认值，可以根据实际情况更改m；

若没有设置初始聚类中心和聚类数K，将采用随机算法生成隶属函数矩阵U，再计算出聚类中心V；第j个聚类中心的计算为：

然后计算此聚类中心下的隶属函数矩阵U；

通过上述计算的隶属函数矩阵U，计算新的聚类中心V*，通过新的聚类中心V*计算新的隶属函数矩阵U*，直到模糊聚类的目标函数J_m(U,V)函数为最小值时，聚类完成；模糊聚类的目标函数J_m(U,V)函数的计算为：

当两次迭代的隶属函数矩阵目标函数J_m，J_m ^*小于设置的容许误差ε时，判定算法收敛，算法终止迭代的判别式为：

J_m-J_m ^*＜ε，ε＝10^-5。

8.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述对预判的生长对象进行预筛选，具体包括：

当生长种子位于一小球面内(R_d1)，则预判的生长对象为大球面内(R_d2)和小球面外的球壳区域；球壳厚度由三角网最大三角形尺寸d_max决定，为保证所有预判的生长对象位于此区域内，取5倍最大三角形尺寸作为预筛选的球壳厚度R₀；

即：

R₀＝5*max(d)；

R_d2＝R_d1+R₀。

9.根据权利要求8所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，所述更换区域生长种子，直到完成结构面分割，具体包括：

I：选取一个三角面片作为生长种子；

Ii：搜索此球壳内与种子三角形共边的三角形，该三角形为生长可能需要的三角形；

Iii：定义生长原则，只有共边的三角形与种子三角形具有模糊聚类结果时才进行区域生长；

Iv：将生长的三角形作为种子三角形重复上述i-iii步骤直到生长停止(球壳内没有搜索到具有模糊聚类结果的共边的三角形)，完成一个结构面的分离；

V：更换新的生长种子重复上述i-iv步骤，直到所有三角形都完成分离，此时所有结构面完成分离。

10.根据权利要求1所述的一种基于高精度点云数据的结构面快速自动识别方法，其特征在于，在得到所有结构面的分离结果后，所述方法还包括：

计算各个结构面的地质信息：

I：产状信息：结构面法向向量n(x，y，z)为其三角面片法向量的平均值，倾向

和倾角θ的对应关系应根据法向向量n所处卦限不同分类处理，具体如下：

ii：结构面出露面积：结构面出露面积为属于它的所有三角面片面积之和。

Iii：结构面迹长信息：以结构面沿倾向或走向方向的最大值长度作为迹长，即将结构面出露的三角形沿倾向或走向方向投影，去投影长度的大值作为迹长。

Images