基于激光点云的隧道开挖面岩体结构识别

杨 辉

(中铁第四勘察设计院集团有限公司， 武汉 430063)

在隧道施工中，由于对隧道地质情况了解不足往往会导致地质灾害的发生。而传统的人工地质素描难以快速客观地反映隧道掌子面的地质状况，需要地质测量人员深入隧道根据自己的经验来判断，而且难以分析掌子面上部的地质状况。随着科技的发展，近景摄影测量、三维激光扫描等非接触式测量技术逐渐受到关注，但在隧道光线不足、扬尘较多等恶劣的环境下，近景摄影测量拍摄的图片会出现图像不够清晰以及图像上的同名点难以匹配等情况[1]，影响后续处理；而三维激光扫描技术不依靠反射棱镜，对隧道内的光线要求不高，可直接采集被测物体表面的三维坐标及灰度值等信息，适合在隧道里推广应用。目前隧道等地下巷道已有应用，主要应用于结构场景还原[2]、结构变形监测[3]、边坡检测[4]等。

近几年，国内外学者利用三维点云数据自动提取室外岩体的结构面和危岩取得了一定的成果。Slob等[5]在获取点云数据后，利用MATLAB重建点云的三角化网络，并计算生成三角面片的产状等信息，使用模糊K聚类算法实现了自动识别；Gigli[6]等基于MATLAB软件开发平台开发了DiAna工具，主要用于2D和3D的结构面半自动识别；董秀军[7]进行了初步尝试，通过人工指定结构面来拟合计算出该面的产状等信息，实现了半自动识别结构面。刘昌军等[8-9]基于柱面投影Delaunay算法对点云数据三角重构，采用改进的模糊 C 平均群聚方法对岩体结构面实现自动识别，并编写了基于IDL语言的结构面群聚统计软件FSS。另外葛云峰等[10]、张溥文等[11]对区域生长法进行改进，控制阈值，对岩体结构面信息进行自动提取和识别。目前国内外研究较多的是针对天然的边坡岩体。天然的边坡岩体具有结构面明显，规模较大的特点。针对隧道内的研究还比较少，由于隧道施工工法和人工机械的因素，导致隧道掌子面区域的岩体较为破碎，原有的结构面遭到破坏，无法直接使用前人的研究。本文借鉴前人识别边坡岩体结构面的方法，根据隧道情况进行了优化，基本能实现隧道掌子面结构面的识别。

1 三维点云数据识别流程

三维激光扫描技术是对被测物体发射激光，然后被测物体表面反射，对被测物体进行快速精确扫描的一种技术，它能够根据扫描所得的信息，精确地重构被扫描物体的表面轮廓点云。其测试原理如下：三维激光扫描仪向被测物体的表面发射激光，仪器内部的模块可以控制横向扫描角度β和纵向扫描角度α，测量发射激光到返回之间时间差计算出被测物体与扫描仪的距离S，如图1所示，三维激光扫描目标点P坐标(xS、yS、zS)的计算公式为

图1 三维激光扫描仪原理图

(1)

通过三维激光扫描仪获取点云数据后，根据点的坐标和灰度值可以重现隧道内部的大部分面貌。一般对点云数据而言，结构面可以认为是由相邻、大致处于同一平面上且具有一定规模的点组成。本研究主要目标是能通过隧道三维点云数据，自动识别出掌子面的结构面并用不同颜色标记，但直接分析点难度较大，可将点云转为有大量三角面片形成的表面模型(TIN)，将海量点云转化为表面模型，极大便利了后续分析。

方法的流程如图2所示，分为4个步骤：①预处理数据；②点云投影；③生成三角网格模型；④识别结构面。

图2 自动识别结构面基本流程

隧道三维激光扫描点云原始数据如图3所示。

图3 隧道三维激光扫描点云原始数据

1.1 预处理数据

三维激光仪扫描的隧道包含了已有的隧道支护结构、隧道开挖区域，已支护的隧道部位表面较为平滑，以破坏原有的岩体结构，故需要先截取开挖区域出露的岩体部分以作后续分析数据的基础。依据隧道设计的标准中轴线数据和扫描仪所处的里程、坐标，计算出掌子面的里程，并保留开挖面后方3 m范围内的点云数据。

隧道内难免存在施工机械、通风管道、施工人员和扬尘等与自动识别结构面无关的噪点，噪点的存在不利于数据处理，甚至对自动识别结构面的结果具有较大影响，对于噪点较明显的可以用距离过滤的方式剔除，而对于依附在隧道表面的噪点，距离过滤无法自动剔除则使用Scene软件人工手动剔除。

降噪后点云的数据量仍然是以百万、千万计，较为庞大，对于识别掌子面的结构面来说，一是处理海量点云数据对电脑硬件、处理算法的效率要求较高且耗时较长，二是识别结构面不需要高精度，可以适当对原始数据采样。本研究对原始点云均匀采样，尽可能减少点云数量的同时也保留了比较明显的特征点。

对降噪取样后的点，根据三维激光扫描工作时的标靶或定位球之间的关系，转换点云数据坐标至大地坐标系。

为了便于后续分析，依据前述截取掌子面区域时所计算中轴线的方向，对点云数据坐标旋转、平移等操作，使得隧道掌子面开挖方向为x轴正方向，扫描中心点为坐标原点。

1.2 点云投影

在隧道中扫描获得(截取)的三维数据点云，其数据在空间中分布于隧道开挖面和衬砌表面，直接采用平面投影方法对所有空间点进行投影，则投影点在任何平面上都不能反映点云在三维空间中的位置关系，若进行Delaunay三角剖分将完全失败。已有学者研究隧道结构三维点云重构，江记洲等[12]采用拼接的柱面投影还原了巷道结构，但不适合带有掌子面的隧道点云数据。根据三维激光扫描仪由扫描中心向四周发射激光的特性，可将点云投影到以扫描中心为球心的球面然后投影到平面上，但一次球面投影只适用所扫描的点云全在一个投影球的半球上，否则一次球面投影而后平面投影会导致点云相对顺序错乱，影响后续生成三角网格模型，而在隧道空间中，由于边墙和拱部的存在，导致球面投影后的点全在一个球面上的条件较为苛刻。

本文采用改进的投影方法能较好地解决这个问题。具体步骤如下：

1)将点云数据第一次投影到以原点为球心，半径为r的球面上。

2)第二次球面投影的圆心坐标设为(-r,0,0)，半径为2r。

3)对于第一次投影后x≥0的点，以原点(0,0,0)为射线的端点、第一次投影球面上的点为射线上一点的射线与半径为2r的球面的交点即为第二次球面投影的投影点。

4)对于第一次球面后x<0的点，以第一次投影球面上的点在x轴上的投影点为射线的端点、第一次投影球面上的点为射线上一点的射线与半径为2r的球面的交点即为第二次球面投影的点。

这样保证了点云经过二次球面投影后全部处于半径为2r球的一个半球上，而后投影到yz平面不会打乱点云间的相对拓扑关系，为后续生成三角网格模型提供便利。如图4所示，原始点云相对拓扑关系为ABCDEFG，第一次球面投影后为A1B1C1D1E1F1G1，而后平面投影为B2A2C2D2E2G2F2，相对关系不是原始的相对拓扑关系。改进后投影方法的如图5所示，采用这种方法，始终保持着原始的相关关系。

图4 第一次球面投影平面示意图

图5 改进的球面投影方法平面示意图

1.3 生成三角网格模型

三角剖分技术是一种应用非常广泛的面重建技术。三角剖分将一些散乱的点云数据剖分为一系列的曲面三角形网格模型。实际中应用最多的是Delaunay三角剖分，它是一种特殊的三角剖分。因为其具备的优点，使其作为目前三角剖分理论的基础，很多三维的剖分优化准则实际上是对它的扩展。本研究生成三角网格模型采用较为成熟的二维Delaunay算法。Delaunay三角判别规则有两个基本准则：

1)每个Delaunay三角形的外接圆不包含任何其他点。

2)Delaunay三角形的最小角度最大化。

在点云投影时，已经建立点云投影前后点之间的对应关系，在生成最终投影的平面点云后，采用Delaunay算法生成三角形，根据投影时建立的对应关系，即可建立三维三角网格模型。如图6所示，是整个掌子面用本研究方法后的表面模型。

图6 生成的隧道掌子面表面模型

1.4 识别结构面

掌子面结构面主要特征是三角面片的法向量大致相同且相连。本研究中以每个生成的三角面片的产状(走向、倾角)以及与三角面片到扫描中心点的距离作为识别参数，产状大致相同，说明三角面片的法向量大致相同、相互平行，但不一定处于同一平面。根据几何知识，点到同一平面的距离相同，若产状和距离都大致相同，相互相连通且成一定规模，则视为是一个结构面。

产状是识别岩体结构面较为重要的参数，计算出每个三角面片的产状，三角面片三点坐标设为Ⅰ(x1,y1,z1)，Ⅱ(x2,y2,z2)，Ⅲ(x3,y3,z3)则此三角形的平面方程为Ax+By+Cz+D=0，根据空间几何知识可以求出A、B、C、D。

设结构面的倾向线与x轴的夹角θ为结构面的倾向；结构面外法线与z轴的夹角φ为结构面倾角。由结构面方程可知，结构面单位法向量为

(2)

(3)

即φ=arccosN，θ=arctan(M/L)，在测量中倾角φ范围为0°～90°，倾向θ范围为0°～360°。

设点(x0,y0,z0)到面的距离为

(4)

本文采用的K均值聚类具有运算时间短、容易实现的优点，适合处理数以万计的点云数据。

具体实现步骤如下：

1)随机确定k的聚类中心μ1,μ2,…,μk。

2)计算每个三角面片的φ、θ、d到每个聚类中心μ1,μ2,…，μk的距离。

3)选择离每个聚类中心最近的三角面片划入子集合。

4)重新计算子集合的聚类中心。

5)重复步骤2)～4)，直到聚类中心不再发生变化或者达到指定的迭代次数。

1.5 三角面片连通性分析及结构面分组

对于结构面不是特别明显的掌子面数据，计算所得的三角面片的聚类参数总体相近，而为了控制代码的运行速度，设置了最大的迭代次数，会导致某些原本应该属于结构面的三角面片未被聚类，需要对这些未聚类的三角面片填补，尽可能地联通结构面。首先找出三角网格中所有未构成结构面的三角面片，分析每一个三角面片的相邻三角是否有构成结构面，并统计生成结构面相邻三角面片个数。对于每一个未生成结构面的三角面片，设其已聚类且相邻的三角面片个数为R，本文采用的填补优先规则如下：

1)若R=3且为同一结构面，则加入该结构面。

2)若R=2且该三角面片与已聚类三角面片夹角(此处夹角均为锐角)其中一个小于阈值，则加入该结构面，同时更新受影响的三角面片相邻已聚类三角面片的个数，若更新后满足规则1，按规则1执行。重复至没有新的三角面片加入。如图7所示，编号1的三角面片加入后，编号2的R=1更新为R=2，同样的，可按照1、6、5、2、4、3填补结构面。

图7 结构面填补示意图

3)若均不满足规则1和规则2，说明该三角面片无法聚类为结构面。

隧道开挖面由于施工影响，一是使得掌子面较为破碎，为了便于分析，去除了规模较小的结构面，保留了规模较大的结构面；二是原本较为连续的结构面变得不连续，没有较为明显的规律。为此本文将相同方向的结构面列为一组，并采用相同颜色表示，不同组采用不同颜色，企图能发现规律。

2 实际效果

为了验证本文所提方法的实用性，现场获取干巴沟隧道的开挖区域点云数据，采集数据所用的三维激光扫描仪是美国FARO公司的Focus3D X 130三维激光扫描仪。三维激光扫描仪共获取6 237 702个的原始点云数据，降噪以及采样后得到35 607个用于测试的点云数据，如图8所示，共生成的三角形个数为71 014，最终生成的表面模型如图9所示。

图8 用于测试的点云数据

图9 点云生成的表面模型

如图10所示，通过K均值聚类算法聚类，用相同的颜色标注产状接近的三角形网格单元。图11是经连通处理过后的结构面，从图中可以看出，共有4组结构面，基本成功识别了主要的结构面，没有标注颜色是不成规模且已经过滤了的三角形，其中蓝色和红色标注的结构面多块成组，根据中间蓝色标准的结构面来看，很可能有一组结构面自图的左上到右下被施工破坏。连通处理后的结构面内部有一些未被聚类的三角形，可能的原因是施工导致的区域坑坑洼洼，无法归入相应的结构面。

图10 结构面识别结果图(未填补)

图11 结构面识别结果图(已填补)

通过本研究计算出结构面规格大的结构面的产状信息，对比该结构面的实际产状信息，对比结果见表1。

表1 本研究和实际结构面倾向、倾角对比

从表1统计结果可看出，通过本研究方法计算出的产状与人工实际测量所的产状基本一致，差异值较小，处于可接受的范围内，说明本研究对自动识别隧道掌子面的结构面具有一定的效果。

3 结论

本文基于三维激光扫描技术获取隧道掌子面的三维点云数据，提出并利用改进的球面投影方法，结合二维Delaunay算法重建了隧道掌子面的三角网格化表面，根据生成的隧道表面模型，使用K均值聚类算法分析了每个三角面片的法向量及距离参数，并识别和获取掌子面主要结构面及其基本信息；而后填补应处于但未被聚类的三角面片，使得结构面有较好的连通性；对产状相似的结构面分组，并用不同颜色表示。结果表明：

1)本文所提出的一种改进的球面投影方法能很好地处理带掌子面的隧道点云数据，还可应用于其他巷道结构中。

2)本研究通过对结构面的填补和分组分析，并用同一种颜色表示同一组结构面，能填补结构面，并发现潜在的结构面走向，为掌子面地质分析提供一定的帮助。

3)采用本方法基本能识别和提取主要结构面及其产状，结构面产状误差处于可接受范围内，能高效快速地为隧道地质编录提供依据，为隧道施工安全提供有力的保障。