基于三维点云建模的矿山边坡稳定性分析

岳西蒙，伍法权，沙 鹏，陈宇坤

(1.绍兴文理学院土木工程学院，浙江 绍兴 312000；2.浙江省岩石力学与地质灾害重点实验室，浙江 绍兴 312000；3.中国建筑第二工程局有限公司，河南 郑州 450000)

0 引 言

随着计算机技术的飞速发展，数值模拟在边坡稳定性分析中应用越来越广泛，解决了许多实际工程问题。目前公认较好的方法有有限元法、有限差分法、离散元法。在岩土工程方面，FLAC3D数值模拟软件表现出强大的价值和实用性，但是其前处理功能并不能够满足复杂的地质条件。为了解决这一难题，崔芳鹏等[1]基于Surfer平台，提取了地表的三维地质信息，经过FLAC3D内置的fish语言对模型的单元和节点进行转换，最终导入到FLAC3D可以直接读取的模型数据文件；郭晓强等[2]利用HyperMesh软件构建不规则的三维岩体模型，并将格式转换为FLAC3D可识别的文件，建立了不规则块体的三维模型；贾曙光等[3]基于无人机倾斜摄影技术总结出了在高陡边坡地区地质信息获取的应用方法，并在此基础上提取了岩体结构面的产状信息；李术才等[4]使用数码图像量化表征了掌子面岩体结构。然而目前的摄影测量技术在工程中的应用仅限于数字表面模型重建方面，实现复杂岩体实体化进行数值模拟依然是一个难题。随着无人机摄影技术的发展，通过无人机可以快速、轻便地获得地质体数字信息，再通过内业处理可以将地质体转换为三维点云模型，为后续的工程地质分析指导提供了强有力的数据支持，为快速复杂地质建模提供了新方法。另外，传统数值模拟方法产出的位移场、应力场等岩体质量指标并不能反应出岩体工程参数。鉴于此，本文将提出一种基于无人机摄影来解决FLAC3D应对复杂地质建模的不足，并通过C++编程实现多项岩体质量量化指标为工程师提供参考。

1 三维点云模型重建

多视角图像的三维重建技术一直是计算机视觉领域的前言话题，通常是利用多视角序列二维图像来恢复目标三维信息，相对于三维激光点云建模有着低成本、高自动化、高效作业等优点。基于无人机摄影测量技术，合理航线规划，来获取高精度多视角航空二维图像，利用多视角运动恢复结构SFM(structure from motion)技术[5]来实现三维点云模型重构。目前比较流行的重建思路是利用多张序列图片进行特征点匹配，进而计算基础矩阵(基础矩阵取决于相机内参K矩阵和外参R矩阵、T矩阵)来映射到三维空间点上，并对像素点稀疏重构，利用MVS(multi view system)算法进行密集重构成密集点云模型。

图像特征点的检测与匹配是解算相机位姿和三维重建的前提条件，其检测与匹配的精度将直接影响后期建模的质量，目前特征检测的算法有很多，其中被大家公认最适合的算法是DAVID[6]总结的尺度不变特征转换算法(SIFT算法)，该算法有着很强的鲁棒性，在目标旋转、缩放、平移、光照影响和视角变化、杂物、噪声等方面表现出很强的稳定性。SIFT算法具体步骤如下所述。①多角度二维图像特征点检测，即从多张图像上提取特征点且不会因为平移、缩放、旋转、光照等噪声的影响而变化，如角点、高曲率点等。首先进行图像尺度空间极值点检测并搜索图像位置，利用高斯微分函数来识别不受噪声影响的兴趣点，然后对每个兴趣点进行评分，分数越高则越稳定，最后结合稳定点的位置关系来确定位置和尺度。②确定方位，根据图像局部梯度方向，给兴趣点分配位置和方向，以此类推后面的操作都相对于兴趣点的位置、尺度、方向进行映射，从而提供对于这些变换的不变性，将提取到的兴趣点生成兴趣点特征向量存储到计算机内存中。③采用欧式距离算法为特征向量匹配，将该算法作为图像兴趣点相似性判定标准，以实现特征点的匹配与融合。通过C++语言实现上述步骤，如图1所示，根据三角测量原理计算出三维点云的坐标信息，即得到稀疏三维点云模型。

由于SFM得到的是稀疏点云模型，在三维地质建模方面并不满足精度要求，为了更加精细地恢复三维场景，需要对稀疏的点云进一步加密，比较常用的方法是基于面片的三维多视角立体视觉算法(PMVS)重建出稠密的点云[7]。该方法利用SFM算法生成的稀疏点云构造空间曲面模型，并以这些空间面片作为种子面向四周扩散，重建出周围的曲面片，最终完成真实空间的三维点云模型，见图2。

图1 稀疏三维点云模型Fig.1 Sparse 3D point cloud

图2 裁切后的稠密三维点云重建Fig.2 Dense 3D point cloud reconstruction after cutting

2 三维地质建模

2.1 规范化拍摄流程

基于无人机摄影测量技术，利用摄影测量相机对目标岩体进行表面模型数字化，其数据获取方式对重建结果至关重要。数据获取流程如下所述。①场景预判：对天气状态、地形条件、周边环境等因素综合考虑。②设备选择：对于低矮的岩体可以采用普通数码相机，对于高耸陡立的岩体可以使用无人机携带相机对目标物进行数据采集并使用5方向拍摄模式(相机正投影90°拍摄一个方向，倾斜45°拍摄四个方向)，GPS弱的峡谷地段可以采用手动控

制云台方式拍摄，但需确保图像之间有足够的重叠度，对精度要求特别高的可采取无人机搭载激光雷达传感器直接获取目标岩体三维点云信息。③规划拍摄路线：根据目标岩体大小不同，对拍摄范围进行合理取舍，对于孤立岩体可采用环绕拍摄或者5方向拍摄见图3(a)，对于岩质边坡可采取斜坡模式拍摄见图3(b)，对于拍摄范围较大的场景可采用仿地飞行模式见图3(c)。④设置飞行参数：使用地面站软件DJIGSPro设置重叠率航向方向80%旁向70%，无人机根据重叠参数可自动拍摄，但在手动拍摄的环境下需要操作人员自行判断，重叠度太高则会造成计算机资源的浪费，重叠度太低则会造成重建失败；航高设置通常比作业区域内最高物体高出50 m；飞行速度可根据航高的不同自行调节，建议设置为8～15 m/s。⑤相机参数设置：相机参数设置正确与否将直接影响三维模型质量与精度，天气状况良好情况下可采取全自动模式拍摄，其他则根据拍照时场景的不同，其IOS、快门，光圈、白平衡等参数应适当调整，应满足以下几点：图像没有局部模糊及畸变情况，在光照不足的情况下不宜将ISO调的过大，尽量避免大景深，白平衡应尽可能接近真实场景。⑥布设像控点，提高模型精度：无人机在RTK模式下小范围内(0.5 km×0.5 km)通常不需要布设地面控制点，在GPS模式下作业应布设地面控制点，可采用作业区域内的可识别辨识物，如公路中间角点，而无可辨特征点区域，应在地面放置可辨的辅助标志见图3(d)，控制点应分布均匀。规范拍摄流程如图3所示。

图3 岩体规范化拍摄流程Fig.3 Standardized photography process of rock mass

2.2 构造FLAC3D实体模型

通过以上步骤即可构造出三维地质点云模型，但是在操作过程中由于人为因素、设备误差、特征点识别和匹配过程的影响，现场获得点云数据会产生一些噪声点，如点云数据密度不均匀、遮挡等问题造成离群点、大量数据需要下采样，此类数据一般不能直接用于建模。对点云数据预处理首先要进行滤波去除周围噪声点、离群点等噪声，噪声点的去除对建模的精度影响很大，因此点云的滤波处理和三维地质模型数字化有着相辅相成的关系，也影响着后续操作。 将三维点云模型在点云处理软件CloudCompare中打开，步骤包括：①选出需要的部分对数据进行裁剪；②对目标点云进行梯度滤波处理，如图4(b)所示，包括平滑、去噪、下采样；③将三维地质模型生成网格并平滑网格(laplacian smoothing算法)，以*.stl文件格式导出，如图4(c)所示。

以上得到的三维地质模型仅为数字表面模型，并不是实体模型，需要对数字表面模型实体化来构造三维实体模型，以实现后续的数值模拟。将得到的三维网格文件(*.stl)导入Rhinoceros软件中，步骤包括：①三维网格进行CAD曲面构建(图4(d))；②曲面进行实体拉伸及布尔运算，并进行犀牛内网格划分(图4(e))；③利用Itasca公司开发的Rhinoceros高级网格生成插件Griddle对三维实体网格模型进行再啮合网格划分，并保存成FLAC3D软件可以识别的*.f3grid网格文件(图4(f))，以浙江省某孤立岩体为例，具体流程如图4所示。

通过分析研究可知，相对于采用有限元软件网格划分的方法，此方法可以快速处理体积庞大、结构复杂的三维地质体，且不需要网格单元、节点转换程序。此方法不仅适用于无人机倾斜摄影测量，对于三维激光扫描同样适用，参照3.2章节建立的实体模型即可；在输出文件格式方面，可导入FLAC3D、3DEC及有限元软件ANSYS、ABAQUS、NASTRAN。

图4 基于无人机摄影的三维地质建模流程Fig.4 3D geological modeling process based on digital photography technology

3 应用实例

3.1 现场数据获取

绍兴市越城区某一大型露天采石场分为东西两个矿坑，废弃矿山几乎包含整个山体(图5)，分布标高32.0～182.5 m，最大高差124.5 m，平面总投影面积281 630 m2，其中宕底面积180 670 m2，人工边坡面积100 960 m2。人工边坡形态较复杂，东侧矿坑和西侧矿于中部贯通相连，形成了平面上呈“镜框”型的东西两个宕口。 两宕口最高均已接近自然山顶，其中东侧宕面最高点高程180 m，边坡坡度一般50°～75°，局部直立状，宕底高程一般为55～67 m，边坡最大高差近120 m；西部宕面最高点高程182 m，边坡地形坡度一般45°～70°，局部呈直立或反倾状，宕底高程46～50 m，边坡最大高差130 m。宕底大多为基岩裸露，局部堆放的浮土(弃渣)厚度一般为0.5～2.0 m。该矿区内节理发育，各组节理统计见表1。根据相关规范，综合地区经验，区内岩石物理力学参数见表2。

通过现场踏勘，采用Phantom 4RTK无人机利用斜坡模式对整个矿区进行摄影测量，为了保证三维模型重建的精度，现场布设了10处控制点，采用UniStrongG970II GNSS RTK对地面点的坐标进行量测，以方便后面数据矫正，其相对误差控制在厘米以内，符合《低空数字航空摄影规范》(CH/Z 3005—2010)的要求[8]。该矿山共拍摄了267张航片，通过SFM-PMPS原理利用ContextCapture软件对航片进行三维重建，共获得67 820 945个三维点云，地质表面点云模型见图6，数字高程见图7。

图5 调查区航拍全貌Fig.5 Aerial photos of survey area

表1 节理产状野外调查统计表Table 1 Statistical table of field survey on joint occurrence

表2 岩石物理力学参数指标Table 2 Physical and mechanical parameters of rock

3.2 地质体建模

三维重建后的模型，受数据量庞大、噪点的影响，一般不能直接进行三维地质建模，要将上述点云曲面模型导入CloudCompare软件进行点云去噪，下采样(本文使用八树叉算法)等处理并生成三角网格，网格的大小决定模型的精细成度，继续将网格做平滑处理(本文使用拉普拉斯算法)，修正后的三维网格模型共包含261 851个小三角网格，如图8所示，将网格以*.stl文件格式保存。

图6 三维表面模型Fig.6 Three dimensional surface model

图7 高程图Fig.7 Digital elevation map(DEM)

图8 修正后部分三维网格模型Fig.8 Modified 3D mesh model

图9 曲面片模型Fig.9 Surface model

图10 有限单元网格模型Fig.10 FEM model

将上述的网格文件*.stl导入到Rhinoceros软件，利用犀牛RhinoResurf命令将网格模型拟合成曲面模型，拟合结果如图9所示。拟合后的曲面片可以直接拉伸成实体，由于边坡并不是平面的，故此时拉伸后的模型底部并不是平面，考虑到数值模拟边界条件的设定，需要对底部及其他部位进行裁平处理，将其与其他实体模型布尔运算实现地质模型实体化。经过多次尝试对新建的地质体模型进行有限元网格剖分，发现ANSYS、MIDAS等有限元分析软件对于大体积、高曲率模型不仅读取效率低下，而且在布尔运算、网格划分时经常报错，严重阻碍了使用人员的工作效率及积极性。Itasca公司开发的Griddle程序是基于Rhinoceros的三维CAD软件通用交互式网格生成插件，结合Rhino的强大复杂结构几何体处理能力及Griddle再啮合Rhino表面网格功能，可以构造出高质量的四面体或者六面体网格。利用Griddle对实体模型进行面划分和体划分，如图10所示。该模型有57 878个节点和296 808个单元体组成，最后以FLAC3D可识别的*.f3grid网格文件格式导出。在FLAC3D软件中通过impgrid xxx.f3grid命令将*.f3grid网格文件导入，即可生成由无人机摄影测量技术得到的精细复杂地质体数值计算模型(图11)。

图11 三维地质体模型Fig.11 3D geological model

3.3 危岩体稳定性分析

基于FLAC3D软件提供的二次开发接口，利用C++语言将统计岩体力学本构模型[9-10]及拓展功能计算模块写入头文件和源文件并编译成FLAC3D软件可自动调用的动态链接库文件JointModel.dll，应用JointModel裂隙岩体本构模型将上述节理几何参数、力学参数及岩石力学参数赋予到地层中，采用静力的方式对边坡稳定性分析，分析结果如下所述。

在两矿山中间类似于“鱼鳍”形状的部位最大位移达到了3.15 cm，且该区域岩体在节理的切割下较破碎，在雨水的冲刷滋润和植物根劈作用下，发生地质灾害的可能较大(图12)。

图12 总位移云图Fig.12 Displacement nephogram

图13为矿山边坡的变形模量云图。图中z方向的变形模量为1 GPa的部位基本上都集中在“鱼鳍”形状附近，再次证明了该区域为地质灾害高风险区域。图14为岩体质量分布云图。图14中“鱼鳍”形状附近的岩体质量分数RMR(SMRMz)为15～30，根据邬爱清等[11]通过大量的工程实测数据分析提出了岩体质量分数BQ与RMR之间的换算关系，计算出z方向的SMRM等级为15～30所对应的BQ值为172.2～263.6，相当于Ⅴ～Ⅳ级的岩体质量，估值RMR为25，为Ⅳ级岩体，此结果与现场勘察结论一致。由此可知，SMRM分级等值线图提供了更多有关岩体质量空间变化的信息，这也表明固定的岩体质量等级并不适合实际情况。

图13 变形模量云图Fig.13 Nephogram of deformation modulus

图14 质量分布云图Fig.14 Nephogram of SMRM quality rating

图15 破坏概率云图Fig.15 Nephogram of failure probability

图16 稳定性系数云图Fig.16 Nephogram of stability cofficient

4 结 论

1) 使用无人机摄影测量技术对矿山边坡进行快速数据采集，总结规范化拍摄流程，实现最优化对矿山边坡三维点云数据的获取。阐述SFM和VMS算法原理，利用ContextCapture三维重建软件实现边坡表面岩体数字化，为后续岩体结构解算、岩体形貌建模及边坡数据建档提供很好的数据支撑。

2) 将岩体三维点云模型结合CloudCompare、Rhinoceros软件，可以快速、精确地实现三维地质模型构建，该方法操作简单、高效、精确，极大地推进了地质体数字化研究进程；应用到FLAC3D数值模拟计算：该矿山岩体受节理影响最大位移位于两矿之间“鱼鳍”形状附近，最大位移为3.15 cm，在雨水滋润冲刷及风化作用下，该部位岩体极有可能出现滑塌、倾倒破坏。

3) 基于FLAC3D数值计算平台利用C++编程对裂隙岩体本构模型开发拓展，将裂隙岩体的变形模量、质量分级、破坏概率、稳定性系数等计算模块写入其中。通过计算，该矿山在结构面切割和人工开挖的影响下最不稳定区域位于两矿连接处“鱼鳍”形状附近，该区域的稳定性系数为1.0以下；变形模量为1 GPa；岩体质量分数SMRMz(RMR)为15～30，并通过经验公式换算出质量分数BQ值为172.2～263.6，为Ⅴ～Ⅳ级的岩体质量级别，该矿山最大破坏概率位于两矿连接处“鱼鳍”形状附近，破坏概率达到43%。通过对比，计算结果与勘察结论基本一致，验证了该方法的有效性，为当地绿色矿山管理部门提供数据参考。