基于三维激光点云的隧道开挖岩体结构面识别与信息提取

Identification and Information Extraction of Rock Structure Surface of Tunnel Excavation Based on 3D Laser Point Cloud

摘要: 岩体结构面几何参数是评价岩体稳定性与渗流特征的重要研究基础。在隧道施工中，识别开挖面岩体结构并分析为后续施工提供了重要参考价值。采用三维激光扫描技术，获取隧道开挖面岩体点云数据，采用球面投影，对投影数据进行三角剖分，得到隧道施工开挖面三角剖分模型，利用结构面产状特征值和距离特征值，采用动态聚类算法实现了岩体结构面分类与识别，并进行结构面产状和面积信息计算。以本方法对某隧道三维激光扫描点云数据进行测试，结合人工测量数据进行比较，结果表明，采用本方法能够快速对岩体结构面进行分组识别和信息提取。

Abstract: The geometric parameters of the rock structure face are an important research basis for evaluating the stability and seepage characteristics of rock masses. In tunnel construction, the identification and analysis of the rock structure of the excavation face provides an important reference value for subsequent construction. A 3D laser scanning technique is used to obtain the point cloud data of the rock face of the tunnel excavation, a spherical projection is used to triangulate the projection data, a triangulation model of the tunnel excavation is obtained, and the structural face is classified and identified using a dynamic clustering algorithm using the structural face eigenvalues and distance eigenvalues, and the structural face yield and area information are calculated. The method was tested on 3D laser scanned point cloud data of a tunnel and compared with manual measurement data. The results show that the method can be used to quickly identify and extract information from the structural faces of rock masses.

1. 引言

结构面作为岩体结构的主体，是岩体的重要组成部分，也是控制岩体强弱和稳定的主因素，其信息是岩体结构划分以及岩体稳定性分析的基础，在很大程度上决定着隧道围岩的定性，地质工程中对施工地带的岩体结构面统计和评价极其重要。传统人工分析开挖面岩体，虽然能实现结构面的识别，但是效率低、工作环境恶劣且勘测人员随时有潜在生命危险判断识别出结构面。本文针对隧道岩体进行结构面识别，虽然隧道掌子面遭到开挖机械破坏，但结构面特征也能部分保留下来。利用表面重建得到三角网格模型，通过统计分析每个三角面片的特征参数，找出参数相同或近似的三角面片，在一定意义上，他们能够构成结构面，从而使别出结构面。

国内外学者对利用点云数据进行结构面识别进行了大量研究。Slob [1] 采用Matlab获取点云数据后，对数据三角网络重组来计算每个三角产状信息，最后采用模糊K均值聚类法实现自动识别；董秀军 [2] 通过人工指定结构面来拟合计算出该面的产状等信息，实现了半自动识别结构面；何秉顺 [3] 使用模糊K均值聚类方法来实现岩体自动识别；葛云峰 [4] 等对区域生长法进行改进，运用解析几何理论，控制阈值，对岩体结构面信息进行自动提取和识别。主要分为区域生长法和结合数码相片基于颜色的识别方法；朱云福 [5] 使用区域生长法，首先需对离散的点云数据进行三角网格表面重建，根据同一结构面的三角形法向量相近，结构面边缘部位三角形法向量突变的特点使用区域生长法实现结构面识别；刘昌军等 [6] 基于柱面投影Delaunay算法对点云数据进行三角网重构，采用模糊C平均群聚方法、FKM模糊群聚法，实现结构面的自动识别；杨文治，赵鹏 [7] 对于平面状结构面，根据结构面的方向各异性和不连续的特点，采用区域生长法识别结构面；对于迹线状结构面，采用手工选点的方式，运用最小二乘拟合对结构面进行识别。

在岩体结构面信息提取方面，描述岩体结构面几何信息的常用参数有产状、间距、尺寸、连通率、粗糙度、张开度等。P.H.S.W. Kulatilake等 [8] [9] 利用三维激光扫描仪进行了岩体结构面粗糙度参数获取；王喆等 [10] 通过三维激光扫描仪所获取的点云数据，利用RANSAC算法、向量代数法，最终得到了岩体的相关地质参数；荆洪迪 [11] 利用Z + Flasercontrol点云数据处理程序和Geomagic Studio逆向工程软件对围岩点云信息进行预处理，基于Java语言开发了由点云数据获取结构面信息的点云数据处理系统，该程序可以获取结构面倾向、倾角等主要产状要素。

岩体结构面识别是信息提取的基础，信息提取是结构面识别的目的。本文在已有研究基础上，基于三维激光扫描技术获取的点云数据，尝试利用岩体结构面的产状特征值和距离特征值，采用动态聚类算法，实现岩体结构面的识别并进行了信息提取，并给出了工作原理，便于工程实践与应用推广。

2. 三维激光扫描技术

本文针对隧道施工开挖面数据进行分析，需要使用三维激光扫描仪采集完整的隧道施工开挖面数据。

三维激光扫描仪是使用激光测距技术，通过发射激光来扫描获取被测物体表面的空间三维坐标和反射光强度的仪器。三维激光扫描仪原理是：三维激光扫描仪发射器发出一个激光脉冲信号，通过反射镜旋转，发射出激光脉冲扫描目标，经被测物体表面漫反射后，沿几乎相同的路径反向传回到接收器。可以得到扫描点P到扫描仪器距离值S，根据两个均匀转动的反射镜可以获取角度值，即每个激光脉冲水平角度α和垂直角度β。三个数据用来计算被扫描目标点的空间三维坐标值。获得点P的坐标原理如公式(1)

$\begin{array}{l}X=S\cos \alpha \cos \beta \hfill \\ Y=S\cos \alpha \sin \beta \hfill \\ Z=S\sin \alpha \hfill \end{array}\}$ (1)

通过三维激光扫描仪获取隧道点云数据后，由于点云数据密度大且离散，对点云数据进行预处理，通过先球面投影后平面投影的方法将三维点云投影到平面，采用Delaunay三角网格剖分，最后还原得到便于后续分析的隧道开挖面三维模型。选取结构面产状特征值和距离特征值，采用动态聚类算法进行结构面的分组识别，并利用数学原理进行结构面的信息计算与提取。

3. 点云数据处理

3.1. 预处理

三维激光扫描仪扫描隧道得到完整的点云数据一般包含了隧道支护结构和隧道开挖面开挖区域，本文主要研究对隧道开挖面处的数据分析，故只保留开挖面区域出露的岩体部分以作后续分析，如图1所示。

图1. 三维激光扫描隧道点云数据

截取开挖面的点云数据后，对点云数据进行降噪处理。剔除掉开挖面附近与识别结构面无关的噪点。由于点云数据量庞大，对点云数据还需进行采样，本文采用模糊熵迭代的方法对原始点云均匀采样，尽可能减少点云数量的同时也保留了比较明显的特征点。

3.2. 点云三角剖分

隧道开挖面岩体不完全在同一平面上，故不能直接采用平面投影点进行三角剖分，将开挖面三维曲面重建。根据三维激光扫描仪由扫描中心向四周发射激光的特性，采用先球面投影，后三角剖分。根据映射关系，即可建立空间开挖面曲面三角网格模型。图2为开挖区域用本文方法后的表面模型展示。

Figure 2. Excavation model for excavation face

图2. 开挖面开挖模型

4. 结构面识别和信息提取

利用表面重建得到三角网格模型，通过统计分析每个三角面片的特征参数.找出参数相同或近似的三角面片，在一定意义上，他们能够构成结构面，从而识别出结构面。在同一结构面上的三角面片，其在空间中的倾角、倾向即产状近似相等，同时，近似在同一平面上。本研究对于结构面分析采用两种特征值，一个是三角形产状，另外一个是三角形的距离(到扫描中心点)。

4.1. 结构面识别

1) 产状特征

在工程地质学上，用走向、倾向和倾角三个要素来确定和描述岩层的产状；而走向和倾向在数值上有关系，走向和倾向相差90˚或者270˚，倾向一旦确定，则走向也相应确定。故本文主要用倾向和倾角来表示平面的同向特征。

倾向用于表示岩体结构面在三维空间中的倾斜方向，角度f值的范围为0~2π (0˚~360˚)。倾角用于表示岩体结构面在三维空间中的倾斜程度，大小为岩体结构面与水平面之间的锐夹角，角度θ范围为0~π/2 (0˚~90˚)。

假设(x1, y1, z1)、(x2, y2, z2)和(x3, y3, z3)是一个三角面片在大地坐标系下的三个顶点的坐标，三角面片所在平面的方程可以表示为：

$Ax+By+Cz+D=0$ (2)

其中(A, B, C)是这个平面的法向量。A、B、C、D值如下：

$\begin{array}{l}A=y1z2+z1y3+y3z3-z2y3-z3y1-z1y2\\ B=z2x3+z3x1+z1x2-x1z2-z1x3-x2z3\\ C=x1y2+y1x3+x2y3-y2x3-y3x1-y1x2\\ D=z1y2x3+z2y3x1+z3y1x2-x1y2z3-y1z2x3-z1x2y3\end{array}$ (3)

单位法向量为：

$\left(x,y,z\right)=\frac{\left(A,B,C\right)}{\sqrt{\left(A^2+B^2+C^2\right)}}$ (4)

在(6)中确定的f的值将始终产生介于0和π/2 (0˚和90˚)或介于0和−π/2 (0˚和−90˚)之间的值，而地理坐标系需要顺时针方向的介于0和2π (0˚和360˚)之间的值。对应关系如下：

若A ≥ 0且B ≥ 0，Q = 0˚；

若A ≥ 0且B < 0，Q = 360˚；

若A < 0且B < 0，或者A < 0且B ≥ 0，Q = 180˚；

2) 距离特征

产状只表示了三角面片的方位，无法表示位置，产状大致相同的三角面片由于互相之间存在距离，不一定属于同一个结构面，根据空间几何知识，点到相互平行但不重合的平面距离不相同，故结构面到扫描中心的距离也可作为特征值。

依据三角形平面方程 $Ax+By+Cz+D=0$，点(x₀, y₀, z₀)到平面的距离d为：

$d=\frac{\left|A{x}_0+B{y}_0+C{z}_0+D\right|}{\sqrt{A^2+B^2+C^2}}$ (7)

3) 聚类分组

聚类是根据所设特征集将具有相似性的元素划为一个类别。因此，聚类分析可用于整合数据点。在本文中，通过设置产状特征值和距离特征值，采用基于K均值聚类方法改进的ISODATA动态聚类算法进行聚类分析，得到岩体结构面的分组情况。ISODATA算法的基本程序如下：

① 设定初始值及初始聚类中心数，但也可在迭代运算过程中人为修改。

② 计算各类中样本的距离函数等指标。

③ 按给定的要求，将所获得的聚类进行分裂和合并处理，以获得新的聚类。

④ 反复重复，重新计算各项指标，判别聚类结果是否符合要求，经过多次迭代运算以后，结果收敛则运算结束。

4) 拟合

对于已分组结构面，由于其方向一致，可依据三角形的距离聚类，提取出每个结构面，采用最小二乘法对平面进行拟合。

设拟合平面方程为： $Ax+By+Cz+D=0$ ( $C\ne 0$)

则 $z=-\frac{A}{C}x-\frac{B}{C}y-\frac{D}{C}$

令 $a_0=-\frac{A}{C}$， $a_1=-\frac{B}{C}$， $a_2=-\frac{D}{C}$

则： $z=a_{0}x+a_{1}y+a_2$

已知结构面上n个点(n ≥ 3)坐标： $\left(x_i,y_i,z_i\right)$， $i=0,1,\cdots ,n-1$

设 $S={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}{\left(a_0{x}_i+a_1{y}_i+a_2{z}_i-1\right)}^2}$，令 $\frac{\partial S}{\partial a_k}=0$， $k=0,\text{1},\text{2}$

$\{\begin{array}{l}{\displaystyle \sum 2\left(a_0{x}_i+a_1{y}_i+a_2-z_i\right)x_i=0}\hfill \\ {\displaystyle \sum 2\left(a_0{x}_i+a_1{y}_i+a_2-z_i\right)y_i=0}\hfill \\ {\displaystyle \sum 2\left(a_0{x}_i+a_1{y}_i+a_2-z_i\right)z_i=0}\hfill \end{array}$

$\{\begin{array}{l}{a}_0{\displaystyle \sum x_i^2+}{a}_1{\displaystyle \sum x_i{y}_i+}{a}_2{\displaystyle \sum x_i={\displaystyle \sum x_i{z}_i}}\hfill \\ a_0{\displaystyle \sum x_i{y}_i+}{a}_1{\displaystyle \sum y_i^2+}{a}_2{\displaystyle \sum y_i={\displaystyle \sum y_i{z}_i}}\hfill \\ a_0{\displaystyle \sum x_i+}{a}_1{\displaystyle \sum y_i+}{a}_{2}n={\displaystyle \sum z_i }\hfill \end{array}$

得到：

$\left(\begin{array}{ccc}{\displaystyle \sum x_i^2}& {\displaystyle \sum x_i{y}_i}& {\displaystyle \sum x_i}\\ {\displaystyle \sum x_i{y}_i}& {\displaystyle \sum y_i^2}& {\displaystyle \sum y_i}\\ {\displaystyle \sum x_i}& {\displaystyle \sum y_i}& n\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right)=\left(\begin{array}{c}{\displaystyle \sum x_i{z}_i}\\ {\displaystyle \sum y_i{z}_i}\\ {\displaystyle \sum z_i}\end{array}\right)$

解该方程组，得a₀、a₁、a₂即得结构面拟合方程组。

4.2. 结构面信息提取

1) 产状

依据结构面拟合的平面方程 $z=a_{0}x+a_{\text{1}}y+a_{\text{2}}$，则法向量为(a₀, a₁, −1)，则该结构面倾角θ和倾向φ为：

$\theta =\arctan \left(\frac{-1}{\sqrt{a_0^2+b_0^2}}\right)$ (8)

$\phi =\arctan \left(\frac{a_1}{a_0}\right)+Q$ (9)

其中，若a₀ ≥ 0且a₁ ≥ 0，Q = 0˚；

若a₀ ≥ 0且a₁ < 0，Q = 360˚；

若a₀ < 0且a₁ < 0，或者a₀ < 0且a₁ ≥ 0，Q = 180˚。

2) 面积

由于开挖面模型是采用三角化模型进行分析，所以结构面面积近似认为是三角面片的面积求和计算得到。通过海伦–秦九韶三角形面积公式即可求得三角形的面积。

设三角形的三个顶点分别为(x₁, y₁, z₁)、(x₂, y₂, z₂)和(x₃, y₃, z₃)，则三条边边长为：

$\begin{array}{l}{l}_1=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2}\\ l_2=\sqrt{{\left(x_1-x_3\right)}^2+{\left(y_1-y_3\right)}^2+{\left(z_1-z_3\right)}^2}\\ l_3=\sqrt{{\left(x_2-x_3\right)}^2+{\left(y_2-y_3\right)}^2+{\left(z_2-z_3\right)}^2}\end{array}$ (10)

根据公式(10)计算三角形周长的一半p为：

$p=\frac{l_1+l_2+l_3}{2}$ (11)

根据公式(10)和公式(11)按照海伦–秦九韶三角形面积公式求出第i个三角形面积S_i为：

$S_i=\sqrt{p\left(p-l_1\right)\left(p-l_2\right)\left(p-l_3\right)}$ (12)

设组成结构面的三角形面积为 $\left\{S_1,S_2,\cdots ,S_n\right\}$，则该结构面尺寸为：

$S={\displaystyle \underset{i=1}{\overset{n}{\sum }}{S}_i}$ (13)

4. 实例分析

为验证所提方法的实用性，对某隧道施工现场连续采集开挖面区域三维点云数据，三维激光扫描仪型号为FAROFocus3DX130。对原始点云图像截取局部开挖面区域点云数据，经降噪处理等得到隧道开挖面重建模型，通过动态聚类对开挖面的产状特征值聚类分组，共得到了5个聚类集合，如图3表示了结构面的分组结果，结构面被分为5组，每组用同一种颜色表示。

从图中可以看到开挖面岩体结构面较为明显，开挖面右半部分结构面组数较多，且小结构面较多，极有可能是由于爆破干扰产生的，这些人为造成的结构面表面往往凹凸不平且数量众多。

为验证本研究所采用的计算产状与面积方法，选取上图开挖面上主要结构面，对比本研究自动计算产状与面积和人工测量的结构面产状与面积，结构面编号如图4所示，对比结果见表1。

Figure 3. Grouping of rock structure faces at the excavation face

图3. 开挖面岩体结构面分组

图4. 主要结构在隧道开挖面中的位置

Table 1. Position of the main structure in the tunnel excavation face

表1. 自动计算与人工测量结果

从表中可以看出，自动计算的倾向、倾角与人工测量的差距不大，最大倾向最大差3.49˚，倾向最大差2.68˚，均在可接受范围内，说明本研究所采用的方法能正确计算结构面的产状；自动计算的面积与人工测量的面积有一定误差；从表中可以看出，大部分的人工计算值比自动计算大，且一般结构面的面积越大，误差越大，主要是由于本研究的结构面是以三角面片为最小单位，可能结构面内部存在较大的凹陷或者凸起，未被归为该结构面，而人工计算时，只能大致估计，并不能如此精细，如结构面S₁；但结构面内部较为平整时，误差并不大，如结构面S₂和S₉。

5. 结论

基于三维激光扫描的隧道点云数据，通过点云处理得到曲面重建模型，利用结构面的产状特征值和距离特征值，采用动态聚类算法对开挖面结构面进行分组识别，再计算出结构面的产状和面积特征信息，实现了隧道开挖面岩体结构面的快速识别和信息提取，主要结论如下：

1) 提出了利用产状特征值和距离特征值进行结构面分组，采用基于K均值聚类改进的动态聚类算法进行聚类研究。

2) 通过分组得到的结构面进行信息提取，得到岩体结构面的产状和面积信息通过工程实例进行分析对比，验证了本方法的实用性，不足之处为本文只对结构面进行了分组识别，而结构面和迹线是有一定的关联的，因此若将结构面和结构面迹线结合分析，相互印证，可提升结果的可靠性。

基金项目

中国中铁股份有限公司科技研究开发计划项目(2019-重大-08-03)。