CN113255677B - 一种岩体结构面及产状信息快速提取方法、设备及介质 - Google Patents

Abstract

本发明公开了一种岩体结构面及产状信息快速提取方法、设备及介质，对所有点云数据的空间范围进行子节点的分割，解决了点云数据因数据量大导致无法处理或处理速度慢的问题，提高了岩体结构面的提取速度；通过峰值聚类算法和DBSCAN算法自动提取边坡点云数据中潜在的结构面多边形，无需人工花费大量时间从海量点云中寻找结构面信息，避免了漏分、误分以及判读标准不一致等问题，极大地提升了结构面提取速度和提取精度。

Description

一种岩体结构面及产状信息快速提取方法、设备及介质

技术领域

本发明属于岩体结构面识别技术领域，尤其涉及一种基于三维激光点云数据的岩体结构面及产状信息快速提取方法、设备及介质，突破以往以人工判读为主、耗时耗力且精度不高的技术难题，主要适用于水电、水利、公路、铁路等行业的工程地质领域。

背景技术

岩体结构面是岩体内部形成的具有不同规模、形态复杂且分布错综的面状地质界面，其产状特征在一定程度上影响着岩体的强度、变形、渗流和稳定性。统计岩体结构面产状信息是地质调查工作中的重要部分，在岩体工程设计、施工等过程中均起到关键性作用。

传统的岩体结构面信息采集手段包括测线法、统计窗法等，通过作业人员在野外现场逐一量测并采集结构面信息，具有工作量大、危险程度高、精度较低及数据完整性较差等缺点。

随着数字摄影测量技术的发展，许多学者以普通数码相机为传感器获取、以全站仪等进行三维坐标量测，实现了非接触式、数字化的岩体结构面信息提取，但存在影像数据获取效率低、精度不足等问题。近期，三维激光扫描技术被广泛应用于地质信息采集中，不仅实现了非接触式、数字化技术，而且能够得到高精度、高效率、高覆盖率点云数据。已有相关研究成果实现了在点云数据上，基于人工判读并手动框选结构面，以及自动计算结构面产状信息的功能。然而，这种方法依旧需要作业人员花费大量的时间从海量点云中寻找结构面信息，存在漏分、误分、判读标准不一等现象。

发明内容

本发明的目的在于提供一种岩体结构面及产状信息快速提取方法、设备及介质，以克服传统方法需要人工从海量点云中寻找结构面信息导致漏分、误分、判读标准不一等问题。

第一方面，本发明提供一种岩体结构面及产状信息快速提取方法，包括以下步骤：

步骤1：获取描述岩体地质几何特征的点云数据；

步骤2：以所有点云数据的空间范围为根节点，将所述根节点分割为多个不同的子节点；对所有子节点进行命名，并保存为点云文件，每个所述子节点对应一个所述点云文件；

步骤3：针对单个点云文件，计算所述点云文件中所有点的法向量，再计算所有点的法向量之间的距离，得到所述点云文件对应的距离矩阵；

步骤4：根据所述距离矩阵计算所有点的法向量的局部密度，并计算每个点的法向量与比该点的法向量的局部密度大的点的法向量之间的空间距离；

步骤5：根据所述局部密度和所述空间距离对当前所述点云文件中所有点进行聚类，得到多个法向量朝向不同的点集，每个所述点集中点的法向量朝向相同或相近；

步骤6：对每个所述点集进行空间划分，得到空间上连续且朝向一致的潜在岩体结构面；

步骤7：对所有潜在岩体结构面进行空间平面拟合，得到多个潜在岩体结构面多边形；

步骤8：对所有点云文件重复步骤3～7，并挑选出具有地质意义的潜在岩体结构面多边形，得到最终的岩体结构面；

步骤9：对最终的所述岩体结构面进行产状信息的计算。

进一步地，所述步骤2中，采用八叉树空间索引将所述根节点分割为多个不同的子节点，具体实现过程为：

将所述根节点分割成八个不同的子节点；

判断每个所述子节点所包含的点云数量是否小于设定阈值，如果是，则停止该子节点的分割；否则继续将该子节点分割为八个不同的子节点，直到所有子节点所包含的点云数量小于设定阈值。

进一步地，所述步骤3中，每个点云文件对应的距离矩阵的获取过程为：

读取所述点云文件，得到对应的点云的数据集S＝{x₁,x₂,…,x_i,…,x_N}，其中x_i为所述点云文件中的第i个点，N为所述点云文件中点的数量；

遍历所述点云文件中的所有点，基于KDTree索引搜索到目标点x_i的邻域点，对目标点x_i的邻域点集合进行拟合得到拟合平面，以所述拟合平面的法向量作为目标点x_i的法向量，得到所有点的法向量集合V＝{v₁,v₂,…,v_i,…,v_N}；其中v_i为点x_i的法向量；

计算两个点的法向量之间的距离，由所述距离得到距离矩阵M＝{d_ij,其中i≠j,x_i∈S，x_j∈S}，所述距离的具体计算公式为：

其中，d_ij为点x_i的法向量与点x_j的法向量之间的距离，(v_ix,v_iy,v_iz)为点x_i的法向量值，(v_jx,v_jy,v_jz)为点x_j的法向量值。

进一步地，所述步骤4中，点的法向量的局部密度的计算公式为：

其中，ρ_i为点x_i的法向量的局部密度，d_ij为点x_i的法向量与点x_j的法向量之间的距离，d_c为截断距离，N_p为点x_i在一定距离范围内邻域点的数量。

优选地，所述步骤4中，空间距离的具体计算过程为：

构建比点x_i的法向量的局部密度ρ_i大的点的集合L，如果集合L不为空集

则选取点x_j的法向量与点x_i的法向量之间距离的最小值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，其中x_j∈L；如果集合L为空集

则以距离矩阵中的最大值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，具体表达式为：

其中，δ_i为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离。

进一步地，所述步骤5中，采用密度峰值聚类算法对所述点云文件中所有点进行聚类，具体实现过程为：

对每个点的法向量的局部密度和空间距离进行归一化处理，将每个点的法向量的局部密度和空间距离映射到区间[0,1]；

由归一化处理后的局部密度和空间距离构建每个点的指数γ_i，得到所述点云文件中所有点的指数集合R＝{γ_i,i＝1,2,3,…,N}，其中，γ_i为点x_i对应的指数，γ_i＝ρ_i×δ_i，ρ_i为归一化处理后点x_i的法向量的局部密度，δ_i为归一化处理后点x_i的法向量的空间距离；

按照指数γ_i的下标对指数集合R进行倒序排列，得到倒序排列后的指数集合R'＝{γ₁,γ₂,…,γ_i,…,γ_N}；

在指数集合R'中寻找变化最大的指数γ_m，其中γ_m∈R'，指数γ_m左侧的每个指数所对应的点均为聚类中心点，指数γ_m右侧的每个指数所对应的点均为非聚类中心点，将每个所述非聚类中心点归属到与该非聚类中心点距离最近的聚类中心点，得到所述点云文件中所有点的类簇划分结果C＝{C₁,C₂,…,C_k,…,C_n}，其中C_k为第k类聚类后的点集，n为所述点云文件的类簇划分结果的数量。

进一步地，所述步骤6中，采用DBSCAN算法进行空间划分，以得到潜在岩体结构面，具体实现过程为：

步骤6.1：分别对每个点集C_k构建样本集S'＝{x₁,x₂,…,x_ii,…,x_T}，其中，T为第k个点集C_k中点的数量，x_ii为对应点集C_k中的第ii个点；

步骤6.2：搜索点集C_k下的每个点x_ii在半径参数ε内覆盖的邻域点集，如果所述邻域点集内的点数大于设定点数，则将点x_ii划分成核心对象集O，否则标记为噪声点；

步骤6.3：如果所述核心对象集O为空，则结束当前点集C_k的分类，潜在岩体结构面的数量为零；

如果所述核心对象集O不为空，则构建队列Q和新簇D，在所述核心对象集O中任选一核心对象O_t，将核心对象O_t归入到新簇D中，获取核心对象O_t在半径参数ε内的邻域点，并将该邻域点加入到队列Q中，将核心对象O_t从核心对象集O中去掉；

步骤6.4：在所述队列Q中任选一点作为目标点，将所述目标点加入到新簇D中，获取所述目标点在半径参数ε内的邻域点；

如果所述邻域点在所述核心对象集O中，则将所述邻域点加入到队列Q并从所述核心对象集O中去掉；

如果所述邻域点是噪声点，则将所述邻域点加入到新簇D中，并将所述目标点从所述队列Q中去掉；

步骤6.5：遍历所述队列Q，重复步骤6.4，直到所述队列Q为空，即完成一个潜在岩体结构面的聚类过程；

步骤6.6：遍历所述核心对象集O，重复步骤6.3～6.5，直到所述核心对象集O和队列N均为空，则结束当前点集C_k的聚类过程。

进一步地，所述步骤7中，空间平面拟合的具体实现过程为：

采用最小二乘法对所述潜在岩体结构面类簇进行空间拟合，得到空间平面方程，并根据空间平面所示点云的空间范围得到潜在岩体结构面多边形，所述空间平面方程的具体表达式为：

ex+fy+gz+h＝0

其中，e、f、g为方程系数，h为常数项，e、f、g不同时为0，且由{e,f,g}构成了该空间平面的法向量。

进一步地，所述步骤9中，产状信息包括走向、倾向和倾角；所述倾向α的具体表达式为：

所述倾角β的具体表达式为：

其中，e、f、g为空间平面方程系数。

第二方面，本发明还提供一种设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如第一方面所述的岩体结构面及产状信息快速提取方法。

第三方面，本发明还提供一种介质，其上存储有计算机程序，该程序被处理器执行时实现如第一方面所述的岩体结构面及产状信息快速提取方法。

有益效果

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

对所有点云数据的空间范围进行子节点的分割，解决了海量点云数据导致无法处理或处理速度慢的问题，提高了岩体结构面的提取速度；八叉树空间索引分割是一种快速的分割方法，进一步提高了岩体结构面的提取速度；

通过峰值聚类算法和DBSCAN算法自动获取在空间上朝向一致且空间上分离的点集，可以在不需要确定各种聚类数量的前提下得到处理任何形状的数据集，寻找出边坡点云数据中的潜在结构面多边形，极大地提升了结构面提取速度和提取精度，无需人工花费大量时间从海量点云中寻找结构面信息，避免了漏分、误分以及判读标准不一致等问题；

自动计算产状信息，进行地质信息的统计分析以及地质编录草图制图，能够为作业人员进行地质信息统计与编录节省大量时间；

本发明具有更高的精度、效率和安全性，在采集结构面信息时完整存储了现场信息。

附图说明

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

图1是本发明实施例中一种岩体结构面及产状信息快速提取方法的流程图；

图2是本发明实施例中基于八叉树空间索引进行点云数据的空间分割示意图；

图3是本发明实施例中某一分割后的点云数据示例，其包含27216个点；

图4是本发明实施例中基于DPC算法聚类得到法向量相近的点集；

图5是本发明实施例中基于DBSCAN算法提取空间分离的潜在岩体结构面；

图6是本发明实施例中对整个边坡自动化提取后得到的潜在岩体结构面多边形；

图7是本发明实施例中人工挑选出的结构面组；

图8是本发明实施例中基于结构面空间位置及产状参数自动生成的地质编录草图。

具体实施方式

下面结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

传统的岩体结构面的信息采集通过作业人员在野外现场逐一量测并采集结构面信息，具有工作量大、危险程度高、精度较低及数据完整性较差等缺点。虽然也有将三维激光扫描技术应用于地质信息采集中的相关研究，但是依旧需要作业人员花费大量的时间从海量点云中寻找结构面信息，存在漏分、误分、判读标准不一等现象。

基于上述技术问题，本发明提供一种岩体结构面及产状信息快速提取方法，通过峰值聚类和DBSCAN算法自动获取在空间上朝向一致且空间上分离的点集，无需人工聚类，大大提高了点集的获取速度和获取精度；通过对点集进行自动空间划分，实现了快速提取边坡点云数据中潜在的结构面多边形，极大地提升了结构面提取速度和提取精度，无需人工花费大量时间从海量点云中寻找结构面信息，避免了漏分、误分以及判读标准不一致等问题。

下面以具体地实施例对本申请的技术方案进行详细说明。下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。

如图1所示，本实施例所提供的一种岩体结构面及产状信息快速提取方法，包括以下步骤：

1、获取描述岩体地质几何特征的点云数据。

在选定岩体地质调查区域后，采用三维激光扫描技术与多幅定焦摄影测量技术对岩体结构周围的边坡、隧道掌子面及其边墙进行扫描和测量，获取描述岩体地质调查区表面几何形态的点云数据。

2、以所有点云数据的空间范围为根节点，将根节点分割为多个不同的子节点，并保存为单独的点云文件。

八叉树空间索引是一种用于描述三维空间的树状数据结构，以根节点表示整个三维空间区域，除根节点以外的每个节点表示一个立方体元素，每个节点又可以分割成八个子节点，子节点的体积之和等于父节点(父节点即为被分割的节点)的体积。

为了解决海量的点云数据计算机无法直接处理或处理速度慢的问题，需要对海量的点云数据进行分割，本实施例中，采用八叉树空间索引将根节点分割为多个不同的子节点，具体实现过程为：

2.1获取点云数据的外轮廓最小外接长方体。

遍历点云数据中所有点的三维坐标值，获取在X、Y、Z轴上的最大坐标值和最小坐标值，将X、Y、Z轴上的最大坐标值和最小坐标值进行组合得到边坡的最小外接长方体，并以此最小外接长方体为根节点。

点云数据中每个点均有三维空间坐标值，坐标值对应的坐标系为大地坐标系。

2.2获取根节点在X、Y、Z轴上的值域区间，分别沿X、Y、Z轴将对应的值域区间进行二等分，得到八个子节点，此时根节点为父节点，八个子节点的体积之和等于根节点的体积。子节点即为子节点空间，是指分割后的点云数据的空间范围。

2.3判断步骤2.2中每个子节点所包含的点云数量是否小于设定阈值，如果是，则停止该子节点的分割；否则继续将该子节点分割为八个不同的子节点(此时，该子节点为父节点，分割得到的八个子节点的体积之和等于该子节点的体积)，直到所有子节点所包含的点云数量小于设定阈值。

如图2所示的某一部分的边坡点云数据进行分割后的示意图，每个矩形框代表一个子节点。如图3所示的某一子节点的点云数据示意图，其包含27216个点。本实施例中，设定阈值为30000个，即每个子节点所包含的点云数量小于30000，则停止该子节点的继续分割。如果子节点所包含的点云数量等于0，则舍弃该子节点。

2.3、对所有子节点进行命名，并保存为点云文件，每个子节点对应一个点云文件。

当分割完成之后，根据子节点的深度和空间位置进行编码。例如，设A表示根节点的外接长方体，在第一次进行分割后，得到A1～A8 8个子节点，如果子节点所包含的点云数量大于设定阈值，则继续进行分割(比如A7)，则继续将A7进行分割并得到A71～A78 8个子节点，直到所有子节点分割完成，并得到每个子节点的编码。

将所有子节点保存为单独的点云文件，并且根据编码规则对点云文件进行命名，每个子节点对应一个点云文件，例如子节点A78对应的点云文件名为A78。

3、对于每个点云文件，计算该点云文件中所有点的法向量，再计算所有点的法向量之间的距离，得到该点云文件对应的距离矩阵。

每个点云文件对应的距离矩阵的获取过程为：

3.1读取该点云文件，得到对应的点云的数据集S＝{x₁,x₂,…,x_i,…,x_N}，其中x_i为该点云文件中的第i个点，N为该点云文件中点的数量。

3.2创建点云的KDTree对象，遍历对应点云文件中的所有点，基于KDTree索引搜索到目标点x_i的邻域点，对目标点x_i的邻域点集合进行拟合得到拟合平面，以拟合平面的法向量作为目标点x_i的法向量，得到所有点的法向量集合V＝{v₁,v₂,…,v_i,…,v_N}；其中v_i为点x_i的法向量；

KDTree是一种分割高维数据空间的索引树形结构，由二分搜索树演变而来，主要应用在高维数据查找领域中的最近邻查找和范围查找。本实施例中，目标点x_i的邻域点数量为30。

3.3基于所有点的法向量集合V＝{v₁,v₂,…,v_i,…,v_N}，计算两个点法向量之间的欧式距离，得到1/2*N*(N-1)个距离值，由这些距离值构成三角距离矩阵M＝{d_ij,其中i≠j,x_i∈S，x_j∈S}，距离的具体计算公式为：

其中，d_ij为点x_i的法向量与点x_j的法向量之间的距离，(v_ix,v_iy,v_iz)为点x_i的法向量值，(v_jx,v_jy,v_jz)为点x_j的法向量值。

4、根据距离矩阵计算该点云文件中所有点的法向量的局部密度，并计算每个点的法向量与比该点的法向量的局部密度大的点的法向量之间的空间距离。

本实施例中，采用高斯指数核函数计算局部密度，可以降低处理小规模数据集时截断距离对点的局部密度计算乃至整个聚类结果的影响。计算公式如下：

其中，ρ_i为点x_i的法向量的局部密度，N_p为点x_i在一定距离范围内邻域点的数量。为节省算法运行效率，基于KDtree索引获取点i在一定欧式距离范围内的邻域点来计算ρ_i值，在本实例中，将此邻域距离阈值设置为0.087，在单位球上对应圆心角度约为5度；d_c为截断距离，一般情况下，如果d_c取太大，则会导致点的法向量密度的区分度并不高；如果d_c值取太小，又会导致差异过大而使每个点单独分成一类。根据经验，将d_c设置为1～2％，对聚类后的结构区别度较好。在本案例中，将距离矩阵M的所有元素倒序排列，取位于2％的值为d_c的值。2％是指在排序中2％位置上的值，假如距离矩阵M有100个元素，则第2个位置上的元素值为d_c的值。

遍历该点云文件中所有点，查找比自身法向量的局部密度大的点，再计算空间距离，具体实现过程为：

构建比点x_i的法向量的局部密度ρ_i大的点的集合L，L为比点x_i的法向量的局部密度ρ_i大的所有点的集合；

如果集合L不为空集

则选取点x_j的法向量与点x_i的法向量之间距离的最小值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，其中x_j∈L；如果集合L为空集

(表明点x_j的法向量的局部密度为最大局部密度)，则以距离矩阵中的最大值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，具体表达式为：

其中，δ_i为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离(简称点x_i的空间距离)。

5、根据局部密度和空间距离对该点云文件中所有点进行聚类，得到多个法向量朝向不同的点集，每个点集中点的法向量朝向相同或相近。

密度峰值聚类算法(DPC算法)是一种基于密度峰值的聚类方法。该算法可以在不需要确定各种聚类数量的前提下得到处理任何形状的数据集，算法简单且高效，故被应用到各研究领域当中。密度峰值聚类算法核心思想如下：1)聚类中心的密度高于其邻近的样本点密度，即一般聚类中心点的密度比局部区域内相邻点的密度要大；2)聚类中心与比其密度高的数据点的距离相对较远。

本实施例中，采用密度峰值聚类算法对该点云文件中所有点进行聚类，具体实现过程为：

5.1对每个点的法向量的局部密度和空间距离进行归一化处理，将每个点的法向量的局部密度和空间距离映射到区间[0,1]。

由于每个点的法向量的局部密度和空间距离可能不在同一个数量级上，因此先进行归一化处理。

5.2由归一化处理后的局部密度和空间距离构建每个点的指数γ_i，得到每个点云文件中所有点的指数集合R＝{γ_i,i＝1,2,3,…,N}，其中，γ_i为点x_i对应的指数，γ_i＝ρ_i×δ_i，ρ_i为归一化处理后点x_i的法向量的局部密度，δ_i为归一化处理后点x_i的法向量的空间距离。

由公式(3)点x_i的空间距离δ_i可知，如果点x_i的局部密度ρ_i是最大局部密度或最大全局密度，则点x_i的局部密度δ_i远大于点x_i的最近邻点的局部密度。因此，类簇中心往往同时具备ρ_i值和δ_i值均较大的特征；反之，选取ρ_i值和δ_i值均较大的点即可得到类簇中心点。故基于ρ_i和δ_i构建指数γ_i＝ρ_i×δ_i，得到R＝{γ_i,i＝1,2,3,…,N}，其中，γ_i代表聚中心权值，γ_i值越大则越可能为簇中心。

5.3按照指数γ_i的下标对指数集合R进行倒序排列，得到倒序排列后的指数集合R'＝{γ₁,γ₂,…,γ_i,…,γ_N}。

5.4在指数集合R'中寻找变化最大的指数γ_m，其中γ_m∈R'，指数γ_m左侧的每个指数所对应的点均为聚类中心点，指数γ_m右侧的每个指数所对应的点均为非聚类中心点，将每个非聚类中心点归属到与该非聚类中心点距离最近的聚类中心点，得到该点云文件中所有点的类簇划分结果C＝{C₁,C₂,…,C_k,…,C_n}，其中C_k为第k类聚类后的点集，n为该点云文件的类簇划分结果的数量。如图4所示，基于DPC算法聚类得到法向量相同或相近的点集。

类簇划分的原理为聚类中心点过度至非聚类中心点会有一个跳跃点(拐点)，通过选取斜率变化最大的点为拐点。

6、对每个点集进行空间划分，得到空间上连续且朝向一致的潜在岩体结构面。

DBSCAN算法是一种经典的基于密度的聚类算法，该算法可以不需要指定类别的数量且可以识别离群点。该算法能够发现任意形状的聚类，适合大数据量的聚类分析。

由于类簇划分结果C中，只是将法向量相同或相近的点划分成一类，但不同空间位置的点平面被聚成一类，而这可能包含了多个结构面。因此，对类簇划分结果C中的每一类C_k，基于DBSCAN算法进行空间划分，以得到潜在岩体结构面，具体实现过程为：

6.1分别对每个点集C_k构建样本集S'＝{x₁,x₂,…,x_ii,…,x_T}，其中，T为第k个点集C_k中点的数量，x_ii为对应点集C_k中的第ii个点。

6.2搜索点集C_k下的每个点x_ii在半径参数ε内覆盖的邻域点集，如果邻域点集内的点数大于设定点数，则将点x_ii划分成核心对象集O，否则标记为噪声点。

注意噪声点之后的分类中可能被发现在其它点x_jj的半径参数ε领域内，如果该半径参数ε领域中点的数量大于设定点数Nminpts，则点x_ii会被加入点x_jj的聚类中。

本实施例中，考虑到存在微小结构面的现象，半径参数ε设置为0.01m，设定点数Nminpts设置为30。

6.3如果核心对象集O为空，则结束当前点集C_k的分类，潜在岩体结构面的数量为零；

如果核心对象集O不为空，则构建队列Q和新簇D，在核心对象集O中任选一核心对象O_t，将该核心对象O_t归入到新簇D中，获取该核心对象O_t在半径参数ε内的邻域点，并将该邻域点加入到队列Q中，将该核心对象O_t从核心对象集O中去掉。

6.4在队列Q中任选一点作为目标点，将该目标点加入到新簇D中，获取该目标点在半径参数ε内的邻域点；如果该邻域点在核心对象集O中，则将该邻域点加入到队列Q并从核心对象集O中去掉；如果该邻域点是噪声点，则将该邻域点加入到新簇D中，并将该目标点从队列Q中去掉。

6.5遍历队列Q，重复步骤6.4，直到队列Q为空，即完成一个潜在岩体结构面的聚类过程。

6.6遍历核心对象集O，重复步骤6.3～6.5，直到核心对象集O和队列N均为空，则结束当前点集C_k的聚类过程。

聚类得到的新簇D即为空间上连续且朝向一致的潜在岩体结构面。

图5展示了将图4中类簇(a)单独提取出来后进行DBSCAN聚类的结果，根据空间位置的不同，共得到九个不同的潜在岩体结构面，用不同的颜色表示，黑色的点表示未被分类的噪声点。

8、对所有潜在岩体结构面进行空间平面拟合，得到多个潜在岩体结构面多边形。

空间平面拟合的具体实现过程为：

采用最小二乘法对潜在岩体结构面类簇D进行空间拟合，得到空间平面方程，并根据空间平面所示点云的空间范围得到潜在岩体结构面多边形，如图6所示，共得到235个潜在岩体结构面多边形，空间平面方程的具体表达式为：

ex+fy+gz+h＝0 (4)

其中，e、f、g为方程系数，d为常数项，e、f、g不同时为0，且由{e,f,g}构成了该空间平面的法向量。

9、挑选出具有地质意义的潜在岩体结构面多边形，得到岩体结构面。

拟合出来的潜在岩体结构面多边形不仅包括了岩体结构面，也包括了其它类型的开挖面，需要人工选出具有实际地质意义的岩体结构面，如图7所示，整个边坡可以将其划分成4组不同的结构面，分别为西北倾向的J1、J2组，和西南倾向的J3、J4组。

从潜在岩体结构面多边形中挑选出具有实际地质意义的岩体结构面，只需要人工简单的判读即可。

10、对岩体结构面进行产状信息的计算。

对挑选出的岩体结构面提取需要的产状信息，产状信息包括走向、倾向和倾角。

走向为岩体结构面与水平面的交线，在空间上是一条无限延伸的线，没有方向性，故以北偏东(NE)或者北偏西(NW)进行表示，具体为：

其中，N、E、W分别表示北、东、西方位。

倾向是垂直走向线沿倾斜层面向下方所引的倾斜线，倾斜线的水平投影线所指的层面倾斜方向就是岩层的倾向，倾向α的具体表达式为：

结构面的倾斜线与水平投影线之间的夹角即为结构面的倾角，倾角β的具体表达式为：

其中，e、f、g为空间平面方程系数。

以平行边坡走向的平面为投影面，将岩体结构面多边形的点投影到投影面上，并标注各组结构面的产状信息，得到地质编录草图(如图8所示)，帮助作业人员进行地质编录图件时提供矢量图支撑。

以上所揭露的仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或变型，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种岩体结构面及产状信息快速提取方法，其特征在于，包括以下步骤：

步骤1：获取描述岩体地质几何特征的点云数据；

步骤2：以所有点云数据的空间范围为根节点，将所述根节点分割为多个不同的子节点；对所有子节点进行命名，并保存为点云文件，每个所述子节点对应一个所述点云文件；

步骤3：针对单个点云文件，计算所述点云文件中所有点的法向量，再计算所有点的法向量之间的距离，得到所述点云文件对应的距离矩阵；

步骤4：根据所述距离矩阵计算所有点的法向量的局部密度，并计算每个点的法向量与比该点的法向量的局部密度大的点的法向量之间的空间距离；

步骤5：根据所述局部密度和所述空间距离对当前所述点云文件中所有点进行聚类，得到多个法向量朝向不同的点集，每个所述点集中点的法向量朝向相同或相近；

步骤6：对每个所述点集进行空间划分，得到空间上连续且朝向一致的潜在岩体结构面；

步骤7：对所有潜在岩体结构面进行空间平面拟合，得到多个潜在岩体结构面多边形；

步骤8：对所有点云文件重复步骤3～7，并挑选出具有地质意义的潜在岩体结构面多边形，得到最终的岩体结构面；

步骤9：对最终的所述岩体结构面进行产状信息的计算。

2.如权利要求1所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤2中，采用八叉树空间索引将所述根节点分割为多个不同的子节点，具体实现过程为：

将所述根节点分割成八个不同的子节点；

判断每个所述子节点所包含的点云数量是否小于设定阈值，如果是，则停止该子节点的分割；否则继续将该子节点分割为八个不同的子节点，直到所有子节点所包含的点云数量小于设定阈值。

3.如权利要求1所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤3中，每个点云文件对应的距离矩阵的获取过程为：

读取所述点云文件，得到对应的点云的数据集S＝{x₁,x₂,…,x_i,…,x_N}，其中x_i为所述点云文件中的第i个点，N为所述点云文件中点的数量；

遍历所述点云文件中的所有点，基于KDTree索引搜索到目标点x_i的邻域点，对目标点x_i的邻域点集合进行拟合得到拟合平面，以所述拟合平面的法向量作为目标点x_i的法向量，得到所有点的法向量集合V＝{v₁,v₂,…,v_i,…,v_N}；其中v_i为点x_i的法向量；

计算两个点的法向量之间的距离，由所述距离得到距离矩阵

M＝{d_ij,其中i≠j,x_i∈S，x_j∈S}，所述距离的具体计算公式为：

其中，d_ij为点x_i的法向量与点x_j的法向量之间的距离，(v_ix,v_iy,v_iz)为点x_i的法向量值，(v_jx,v_jy,v_jz)为点x_j的法向量值。

4.如权利要求1所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤4中，点的法向量的局部密度的计算公式为：

其中，ρ_i为点x_i的法向量的局部密度，d_ij为点x_i的法向量与点x_j的法向量之间的距离，d_c为截断距离，N_p为点x_i在一定距离范围内邻域点的数量；

优选地，所述步骤4中，空间距离的具体计算过程为：

构建比点x_i的法向量的局部密度ρ_i大的点的集合L，如果集合L不为空集

则选取点x_j的法向量与点x_i的法向量之间距离的最小值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，其中x_j∈L；如果集合L为空集

则以距离矩阵中的最大值作为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离δ_i，具体表达式为：

其中，δ_i为点x_i的法向量与比该点x_i的法向量的局部密度ρ_i大的点x_j的法向量之间的空间距离。

5.如权利要求1所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤5中，采用密度峰值聚类算法对所述点云文件中所有点进行聚类，具体实现过程为：

对每个点的法向量的局部密度和空间距离进行归一化处理，将每个点的法向量的局部密度和空间距离映射到区间[0,1]；

由归一化处理后的局部密度和空间距离构建每个点的指数γ_i，得到所述点云文件中所有点的指数集合R＝{γ_i,i＝1,2,3,…,N}，其中，γ_i为点x_i对应的指数，γ_i＝ρ_i×δ_i，ρ_i为归一化处理后点x_i的法向量的局部密度，δ_i为归一化处理后点x_i的法向量的空间距离；

按照指数γ_i的下标对指数集合R进行倒序排列，得到倒序排列后的指数集合R'＝{γ₁,γ₂,…,γ_i,…,γ_N}；

在指数集合R'中寻找变化最大的指数γ_m，其中γ_m∈R'，指数γ_m左侧的每个指数所对应的点均为聚类中心点，指数γ_m右侧的每个指数所对应的点均为非聚类中心点，将每个所述非聚类中心点归属到与该非聚类中心点距离最近的聚类中心点，得到所述点云文件中所有点的类簇划分结果C＝{C₁,C₂,…,C_k,…,C_n}，其中C_k为第k类聚类后的点集，n为所述点云文件的类簇划分结果的数量。

6.如权利要求1～5中任一项所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤6中，采用DBSCAN算法进行空间划分，以得到潜在岩体结构面，具体实现过程为：

步骤6.1：分别对每个点集C_k构建样本集S'＝{x₁,x₂,…,x_ii,…,x_T}，其中，T为第k个点集C_k中点的数量，x_ii为对应点集C_k中的第ii个点；

步骤6.2：搜索点集C_k下的每个点x_ii在半径参数ε内覆盖的邻域点集，如果所述邻域点集内的点数大于设定点数，则将点x_ii划分成核心对象集O，否则标记为噪声点；

步骤6.3：如果所述核心对象集O为空，则结束当前点集C_k的分类，潜在岩体结构面的数量为零；

如果所述核心对象集O不为空，则构建队列Q和新簇D，在所述核心对象集O中任选一核心对象O_t，将核心对象O_t归入到新簇D中，获取核心对象O_t在半径参数ε内的邻域点，并将该邻域点加入到队列Q中，将核心对象O_t从核心对象集O中去掉；

步骤6.4：在所述队列Q中任选一点作为目标点，将所述目标点加入到新簇D中，获取所述目标点在半径参数ε内的邻域点；

如果所述邻域点在所述核心对象集O中，则将所述邻域点加入到队列Q并从所述核心对象集O中去掉；

如果所述邻域点是噪声点，则将所述邻域点加入到新簇D中，并将所述目标点从所述队列Q中去掉；

步骤6.5：遍历所述队列Q，重复步骤6.4，直到所述队列Q为空，即完成一个潜在岩体结构面的聚类过程；

步骤6.6：遍历所述核心对象集O，重复步骤6.3～6.5，直到所述核心对象集O和队列N均为空，则结束当前点集C_k的聚类过程。

7.如权利要求1～5中任一项所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤7中，空间平面拟合的具体实现过程为：

采用最小二乘法对所述潜在岩体结构面类簇进行空间拟合，得到空间平面方程，并根据空间平面所示点云的空间范围得到潜在岩体结构面多边形，所述空间平面方程的具体表达式为：

ex+fy+gz+h＝0

其中，e、f、g为方程系数，h为常数项，e、f、g不同时为0，且由{e,f,g}构成了该空间平面的法向量。

8.如权利要求7所述的岩体结构面及产状信息快速提取方法，其特征在于，所述步骤9中，产状信息包括走向、倾向和倾角；所述倾向α的具体表达式为：

所述倾角β的具体表达式为：

其中，e、f、g为空间平面方程系数。

9.一种设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1～8中任一项所述的岩体结构面及产状信息快速提取方法。

10.一种介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1～8中任一项所述的岩体结构面及产状信息快速提取方法。

Images