基于数字表面模型的岩体结构面产状获取

宣程强 ，章杨松 ，许文涛

（南京理工大学土木工程系，江苏 南京 210094）

岩体结构面参数主导岩体结构稳定性并控制岩体裂隙中流体的渗流路径，对岩土工程研究有着至关重要的作用。岩体结构面参数通常由专业的工程师使用皮尺、罗盘等工具现场测量获取。但该方法工作量较大，在一些地势较为复杂的地段如峭壁、隧道难以使用，同时需要量测人员长时间置身野外环境，受环境因素影响较大，也带来了较大的安全隐患[1-2]。

目前三维激光扫描技术（light detection and ranging）[3-5]和数字摄影测量技术（digital photogrammetry）[6-8]被广泛应用于岩土工程领域。三维激光扫描技术基于激光测距原理，记录被测物体表面大量密集点的三维坐标、反射率和纹理等相关信息，可以精确重构岩体三维点云模型。然而其缺点也十分明显，性能优异的三维激光扫描仪及其配套设备软件价格十分昂贵，不利于其技术普及[9]。随着数字摄影测量技术和计算机技术的不断发展，基于运动法（structure from motion，SFM）的计算机视觉三维重建算法使获取目标物三维模型的方法变得更加简单，只需利用高清数码相机从不同角度拍摄目标物体就可以还原出其三维模型。该方法与传统摄影测量方法相比有较大的便利性，不需要事先确定拍摄相机的空间位置和姿态，通过对拍摄得到的一组有一定重叠度的目标物图像迭代计算就能自动求解出相机内外参数，十分适合有场地空间限制的复杂工程环境[10]。因此本文结合数字摄影测量技术与SFM算法，使用高清数码相机拍摄岩体影像，实现岩体三维模型的重建。

近些年来，许多研究者已经对从岩体三维模型上手动、半自动或者自动提取结构面信息展开了大量研究。杨文治等[11]通过手动选点的方式从岩体三维激光点云模型表面提取结构面。Cacciari等[12]使用激光雷达扫描仪对隧道围岩进行扫描建模，采用手动的方式提取围岩模型表面的迹线并投影到二维平面进行分析。显然，对于专业技术人员来说，采用人机交互的方式从岩体三维模型表面提取结构面参数会有更高的准确度。但是这种方式自动化率较小，无法快速获取大规模的结构面参数。目前比较常见的自动或半自动识别岩体结构面方法主要基于点云的直接分类[13]，包括另一种形式的三角网格分类[14]，都是通过对岩体点云进行平面子集的划分，然后基于平面子集的法向量方向采用不同的聚类方法（如模糊聚类方法，K均值聚类方法）实现结构面的分割。该类方法在点云模型规模不大的情况下具有较高的准确度，但当点的数目达到百万、千万级别时，每个点的计算和分类需要花费较长时间。

综上所述，本文基于数字摄影测量技术与SFM算法进行三维重构，建立精度较高岩石数字表面模型并基于此开发结构面自动识别算法。

1 基于SFM的岩体数字表面模型重建

1.1 岩体影像资料采集

传统摄影测量方法基于双目立体视觉原理通过双摄像机或者单摄像机从不同角度获取目标物体的两幅数字影像进行三维重建。然而，其前提条件是必须事先获得摄像机内外参数。相比之下，运动法（SFM）通过特征匹配从两幅或者多幅图像中提取相同的特征点并使每个特征点之间相互关联，自动求解出摄像机的内外参数，然后基于模型的一些特征点与被测目标物体之间的比例关系就可以计算出所有特征点的空间坐标[10]。针对巷道内部复杂、狭小的空间环境，为构建真实比例大小的高精度巷道三维表面模型，采用如下拍摄流程：

(1)在巷道围岩测量区域表面均匀布置一定数量的控制点和检查点，并利用全站仪获取其真实世界坐标。

(2)手持高清数码相机从不同角度拍摄巷道围岩影像。拍摄过程中尽量保证相邻影像覆盖范围的重叠度大于30%，所采集到的影像能够覆盖整个测量区域。

(3)经过一轮拍摄后，再对巷道围岩控制点所在区域以及一些关键结构面位置进行二次补拍。

拍摄所用的数码相机为尼康COOLPIX P310，镜头焦距为5 mm。总共采集了161张巷道围岩照片，每一张照片的像素为4 608×3 456。

1.2 基于SIFT的图像特征匹配

在利用SFM算法进行三维重建之前，需要对获取的岩体图像进行特征点的检测和匹配，主要利用尺度不变特征变换(Scale-Invariant Feature Transform，SIFT）算法实现。该算法由Lowe提出，通过尺度空间的极值检测从图像中提取特征点的尺度和旋转不变量，并利用高斯滤波器(Difference of Gaussian，DOG)计算得到特征点的位置，然后通过比较特征点的描述子（Descriptor）来匹配多幅图像中的公共点。

尺度空间即为试图在图像领域中模拟人眼观察物体的概念与方法，目的是检测图像在尺度变化时的稳定特征。为了在尺度空间中检测到稳定的关键点，通过使用不同的尺度参数对图像进行高斯模糊处理生成图像的高斯金字塔，然后对高斯模糊后的相邻图像进行差分得到图像的高斯差分金字塔，再判断差分后图像中每一个点与其26个邻域像素点的关系是否满足最大或者最小并保留满足的特征点为候选关键点。在确定了图像关键特征点以后，需要对特征点进行方向赋值。通过计算获取关键特征点邻域像素的梯度方向以及幅值，比较各个方向的幅值大小，幅值最大的方向即为该特征点的主方向[15]。图1为利用SIFT算法计算得到的围岩图像局部特征点及其主方向矢量，图中箭头的朝向即为特征点的主方向，长度代表其尺度大小。

图1 SIFT特征点向量Fig.1 SIFT feature point vectors

在获取图像关键特征点位置和特征向量后，通过比较不同图像特征点向量的欧氏距离来判断特征点是否匹配。为提高匹配的正确率，采用设立距离阈值的方法比较特征点之间欧氏距离的最近邻距离与次近邻距离，当距离小于该阈值时，则认为匹配成功[16]。当对所有图像完成关键特征点获取与匹配后，通过计算就能得到大量目标物体图像特征点的内外方元素和坐标。

1.3 表面模型构建

从目标物图像到其三维模型生成主要经历以下几个流程：①稀疏点云构建；②点云稠密化；③曲面重构。

稀疏点云的构建主要依靠SFM算法实现，首先从大量经过特征匹配的图像中选取种子图像进行初始位姿估计，然后向整个系统中不断添加新的图像进行特征匹配，同时对新添加图像中匹配到的特征点进行三角定位和位姿估计，直至所有图像添加完毕[17]。种子图像的选取很大程度上决定了后续建模的精准度，其选取原则主要有以下两点：(1)两图像间匹配的特征点个数；(2)图像间的基线宽度。通常两个图像间匹配的特征点越多越有利于相机位姿的求解，得到的结果也更加精确。巷道围岩稀疏点云模型如图2(a)所示。

点云稠密化主要基于面片的三维多视角立体视觉(Patch-Based Multi-View Stereo，PMVS)算法实现，PMVS算法的输入是经SFM处理得到的稀疏点集、图像以及相关参数，通过匹配、膨胀、过滤三个步骤完成在局部灰度一致性下的特征点密集匹配，最终输出带有颜色的密集点云[18]。如图2(b)所示为通过点云加密得到的巷道稠密点云模型。该算法几乎实现了图像中每一个像素点的匹配，重建了每一个像素点在三维空间中的位置，使获取的点的密集程度基本达到了原有图像的清晰度。

图2 巷道三维重建模型Fig.2 3D reconstruction model of roadway

巷道三维模型的曲面重构由MeshLab网格编辑软件[19]实现，利用该软件中强大的三角网格编辑功能可以轻松实现由点云到数字表面模型的转化。经统计最终得到的巷道数字表面模型表面共有1 376 877个顶点，2 709 674 个三角网格。

1.4 建模精度分析

为验证上述结合摄影测量技术与SFM算法重构巷道三维数字表面模型方法的有效性，对5个检查点的外业实测坐标与重构得到的三维坐标进行了误差比较，如表1所示。从表1中可以看出，检查点的中误差都在0.01 m以内，模型的精度可以满足实际工程需要。

表1 三维模型检查点误差统计Table 1 Error statistics of the 3D model checkpoints

2 结构面识别算法

结构面识别算法程序基于Matlab2016b平台编写，测试环境个人PC主要配置如下：64位WIN10操作系统，CPU为Intel core i5-9300H，GPU为NVIDIA GeForce GTX1660Ti，内存16 GB。算法基本流程如图3所示。

图3 结构面识别算法流程图Fig.3 Flow chart of discontinuities recognition

2.1 平滑处理

通过激光雷达扫描仪或者摄影测量方法获取的目标物三维表面模型由于一些人为操作问题和环境因素影响或多或少会产生一些噪声[20]。显然，这些噪声相对于整个模型来说占比十分渺小，对于目标物模型的三维可视化产生不了很大的问题，但对于后续算法分析组成整个三维表面模型的每一个三角网格来说，噪声所带来的高频偏差会使得其表面法向量与实际产生较大的出入，从而导致结构面搜索的中断。其中一些孤立的离群噪声点比较容易区分，可以利用点云编辑软件，如CloudCompare，通过自动或者手动的方式去除。另一些存在于模型表面的噪声由于其特征不清楚，难以通过上述方式消除，所以本文利用拉普拉斯平滑算法对表面模型进行平滑处理以降低噪声对后续算法开展的影响。

拉普拉斯平滑算法是计算机数字几何处理领域中比较常用的网格平滑算法，该算法以三角网格中的伞状结构为单元，先求出模型中所有伞状结构顶点Poriginal和 构成其相邻面的顶点的重心Pmean，然后将该重心坐标作为新的三角网格顶点Pnew替换原来顶点Poriginal的位置，从而达到平滑三角网格的作用。本文采用阻尼迭代拉普拉斯平滑算法对岩体模型三角网格进行平滑处理，该算法是由拉普拉斯平滑算法变形而来，由于基础的拉普拉斯平滑算法在平滑网格时可能会因其单次平滑处理对模型表面点的位置有较大影响而使模型表面三角网格发生重叠，所以引入一个加权系数 λ，通过控制该系数的大小改变单次平滑算法对模型顶点位置影响，降低三角网格发生重叠的概率，其定义如下：

平滑算法阶段需要用户控制的变量只有一个迭代次数I，首先建立一个临时点集用于存储平滑处理后模型中所有点的位置，当迭代次数从Ii→Ii+1时，对临时点集中的点进行同样的平滑处理得到新的顶点位置并重新存入临时点集，直至迭代次数到达用户设定的上限，平滑即算完成。随着迭代次数的增加，三角网格顶点的位置分布就更加均匀，从而得到更加平滑的模型。图4为加权系数 λ=0.3时在一小块巷道围岩数字表面模型表面分别作用0，10，30次平滑迭代处理的效果对比图。

图4 平滑岩体数字表面模型Fig.4 Smoothing rock mass DSM

2.2 平面筛选

平面筛选阶段的作用主要有两点：①分隔不同的平面区域；②平坦性检测，为后续结构面区域生长算法挑选种子点集。主要由边缘检测和平面检测两个步骤实现。

2.2.1 边缘检测

为了保证检测到的平面只代表一个潜在的结构面，必须利用平面边缘进行分割，防止在结构面区域生长阶段出现跨区域搜索现象。本文首先通过最小二乘法拟合模型伞状结构单元上所有节点组成的平面以获取每一个三角网格顶点的法向量vi，然后基于每一个顶点，设定一个较小的搜索半径(由生成模型的网格间距决定)，获取搜索半径内所有顶点法向量的平均值为顶点个数），计算每个顶点法向量与平均法向量的角距离，由式(2)可得，并将其平均值与用户设定角度阈值 βmax比较，如果大于该阈值βmax，则认为该点靠近边缘。边缘检测结果如图5(a)所示，图中绿色区域即为边缘检测算法搜索到的岩体数字表面模型上的平面边缘区域，剩余的红色区域将用于下一个阶段的平面检测。

2.2.2 平面检测

平面检测方法和边缘检测类似，只是需要设定较大的搜索半径和较小的角度阈值 βmin，如果计算得到的每一个三角网格顶点的法向量vi与搜索半径内所有顶点法向量的平均值vmean的平均角距离小于阈值 βmin，则认为该点在平面内。由于在平面检测前已经将平面边缘部分剔除，所以最终检测到的平面应当是独立分开的。图5(b)中红色区域部分即为最终检测到的潜在结构面，将用于下一步结构面区域生长算法进行最终判定。

图5 结构面识别过程Fig.5 Processes of discontinuities recognition

2.3 结构面搜索

结构面搜索基于区域生长法原理，区域生长法的基本思想是：在一个区域中，确定一个种子点作为生长起点，将种子点邻域内的候选点与种子点比较，当候选点与种子点具有相似属性时，合并候选点并将其作为下一个生长点进行向外生长。本文采用的结构面区域搜索算法主要包括三个步骤：内搜索、外搜索、面积判定。

内搜索主要目的是判断潜在结构面区域内部的每个点是否满足构成结构面的要求。首先为每个潜在结构面区域建立临时点集，然后在每个独立的潜在结构面区域内确定合适的种子点并赋予其中心法向量vcenter，vcenter通过计算各自所在平面区域范围内所有顶点法向量的平均值得到。计算种子点中心法向量vcenter与其邻域内候选点法向量vi的角距离并与设定的角度阈值 β进行比较来判断是否合并该候选点。如果其角距离小于阈值 β，则将该候选点存入临时点集并作为下一个生长点继续向外生长以不断扩大结构面区域。根据Riquelme等[13]关于平面子集共面性的结论，通过反复试算，本文最终采用的阈值 β=20°。当所有结构面区域大小均不再发生变化时，即可判定内搜索已经完成。

外搜索的目的是对结构面边缘进行补全，由于在平面筛选阶段为了分隔平面剔除了一部分靠近平面边缘的点使其没有参与后续计算，这些点中的部分有可能也满足结构面区域判定条件，可以被归入相应的结构面区域中。采用相同的判定条件对结构面边缘附近的点进行判定并将符合条件的点存入相应的结构面临时点集中。

最后计算每个结构面临时点集构成的三角网格区域大小，与设定的面积阈值Smin进行比较，当网格区域面积大于Smin时，则认为该结构面有效。图5(c)为在岩体数字表面模型搜索到的结构面，每个结构面都用不同的颜色在模型中表征，同图5(b)比较，可以看到，结构面的边缘部分得到了很好的补全，保证了后续计算结构面产状结果的准确性。

2.4 岩体结构面拟合与转换

在搜索得到构成每一个结构面的点集后，基于随机采样一致性（random sample consensus, RANSAC）拟合计算岩体结构面法向量。

RANSAC算法最大的特点就是随机性和假设性，其假设一组样本数据中包含正确数据和噪声，同时也假设存在相应的方法可以计算出符合这些样本数据的模型参数[21]。通过随机抽取组成结构面点集中的任意三个点构成一个参考平面作为样本模型参数，比较点集中其余点与该参考平面的距离。这里需要提前设置一个距离阈值，当点与参考平面的距离小于阈值时，视该点为正确点，反之，则该点为无效点。经过多次随机抽取，存在最多正确点的参考平面就是结构面的最佳拟合平面。

由于暴露在岩体表面的结构面通常比较粗糙，成波浪状，而RANSAC算法考虑了点的波动性[14]，其平面拟合结果更加符合实际情况，所以采用该算法进行岩体结构面拟合是十分合适的。图6是RANSAC算法拟合结构面示意图。

图6 RANSAC算法拟合结构面点云Fig.6 Fitting a plane to the point cloud of a discontinuity by RANSAC

获取了结构面法向量之后就要将其换算成产状，如图7所示建立三维空间坐标系，X轴指向正东，Y轴指向正北，Z轴竖直向上。

图7 法向量与产状关系Fig.7 Relationship between normal vector and orientation

用一个圆盘代表结构面，v为结构面的单位法向量，θ 和 α分别是结构面的倾角和倾向。根据空间解析几何理论，很容易得到结构面产状与法向量的计算关系：角 θ 和倾向 α的公式：

相应地，可以推导出岩体结构面法向量转换成倾

2.5 可靠性验证

为验证本文提出结构面产状提取算法的有效性，分别使用人工现场测量方法和Riquelme等[13]开发的DSE软件对同一部位岩体的6个结构面（在图5(c)中用黑色圆圈标号）产状结果进行采集。其中DSE软件的结构面分割效果如图8所示，最终的结构面产状对比结果如表2所示。从表2中可以看出，本文方法所获取的产状同DSE软件提取的结果较为接近，与人工现场测量结果相比，倾向最大误差为6.535°，平均误差为2.812°；倾角最大误差为3.789°，平均误差为2.121°。

图8 DSE分割结构面点云Fig.8 Segmenting discontinuities point cloud using DSE

表2 结构面产状测量比较Table 2 Comparison of the discontinuities orientation measurements

由于岩体结构面通常较为粗糙，呈弯曲状，而人工接触式测量往往只对结构面单点进行数据采集，其结果准确性无法得到有效保证。本文所提出的产状提取方法考虑了单个结构面上所有点的影响，其结果是可以被接受的。

3 工程应用

北山坑探（图9）位于甘肃北山旧井块段的十月井Ⅰ级断裂带附近，是研究地下岩体特性、结构面分布规律以及进行现场测试的基础试验平台[22]。北山坑探主要由斜坡道和水平巷道组成，在开挖巷道揭露特征明显的花岗岩结构面。对巷道开挖部分进行三维模型重构，在此基础上进行智能化提取裂隙节理参数，可为后续的岩体质量评价、裂隙网络建模提供重要的基础参数，同时也为本文提出结构面智能化识别算法的适用性提供一个典型的案例研究。

图9 北山坑探现场测量工作照片Fig.9 Field measurement work photo in the Beishan exploration tunnel

本文中结合数字摄影测量技术与SFM算法重构了一段长约8 m的巷道数字表面模型，如图2(c)所示。通过对模型旋转、放大，可以更清楚地观察巷道内部各个位置岩体的节理特征。

对整个巷道数字表面模型进行平滑、边缘检测、平面检测和结构面搜索，共计识别得到120 个具有明显特征的面状结构面。每个结构面都用不同的颜色加以区分表征在模型上，如图10所示。

图10 巷道围岩结构面识别结果Fig.10 Roadway surrounding rock discontinuities recognition result

图11 结构面产状聚类分组极点图Fig.11 Pole distribution of discontinuities orientation cluster grouping

表3 结构面产状分组结果Table 3 Results of discontinuities orientation grouping

将巷道DSM上的结构面按照极点图分组结果使用不同颜色进行表征，如图12所示，可以让工程人员方便、快捷地判断巷道内部结构面的分布情况。

图12 巷道围岩结构面分组表征Fig.12 Roadway surrounding rock discontinuities mapped by groups

4 结论

(1)精度较高、特征点识别丰富的岩体DSM模型对岩体结构面识别工作十分重要，在对目标岩体进行摄影测量时应尽量保证每张照片之间有一定的重叠度，拍摄的仰角不要过大，最好使镜头与岩体结构面保持平行或近似平行。

(2)结构面产状的提取流程主要为：①平滑岩体数字表面模型；②边缘检测；③平面检测；④基于区域生长原理进行结构面搜索；⑤基于RANSAC算法进行结构面拟合。通过与人工实地测量方法以及现有的结构面识别软件对比，该方法具有较高的准确度，能够满足实际工程需要。然而，本文提出的方法并不是要代替地质工程师的知识和经验，旨在为他们提供一个高效、准确的测量岩体结构面参数的工具，减少由于人为因素而产生的错误。最后将结构面按照聚类结果表征在岩体三维表面模型上也提供给工程人员一种随时随地查看现场结构面分布的途径。

(3)考虑到算法的局限性，本文所提出的结构面识别方法还只适用于识别形状接近于平面的面状结构面，还无法有效识别闭合的线状结构面，这也是作者未来需要突破的重点。