岩质边坡节理信息数字识别及三维模型重构<br/>——以石家庄市南障城镇边坡为例

岩质边坡节理信息数字识别及三维模型重构
——以石家庄市南障城镇边坡为例

2023-03-15刘宁袁维高岭白雪山翟星李建朋王安礼徐江

科学技术与工程 2023年3期

刘宁,袁维*,高岭,白雪山,翟星,李建朋,王安礼,徐江

(1.石家庄铁道大学土木工程学院，石家庄 050043；2.河北省交通规划设计院有限公司，石家庄 050300；3.河北省地质环境监测院，石家庄 050022；4.河北省地质资源环境监测与保护重点实验室，石家庄 050022；5.贵州省喀斯特地区水资源开发利用工程技术研究中心，贵阳 550002)

在大型工程建设中，无论是公路铁路还是水利水电都难免会跨越或者下穿岩质边坡，所以需要保证岩质边坡具有良好的稳定性。岩质边坡稳定性分析一般需要首先进行野外勘察收集数据，然后进行相关试验分析获取力学参数，最后运用数值模拟等方法对岩质边坡进行计算评估。野外勘察收集数据通常采用人工测量和原地取试样的方法[1]，该方法会受到复杂的地形地貌条件制约，使正常的勘察调查工作无法开展，并会威胁测量人员的安全，而试验环节获取参数的过程漫长，且存在人为操作的误差，影响试验结果，因此，需要找到一种快捷、高效的获取地质数据并得到相关参数的方法。

国内外学者对上述问题进行了广泛的研究，无论是地形地貌的测量方面，还是根据地质数据进行相关参数的提取方面都取得了不少成果。梁玉飞等[2]把激光三维扫描技术应用到地质工程中，通过激光三维扫描获取研究区域滑坡的点云数据，以此为基础进行节理信息的提取，为地质灾害分析提供了数据支撑。嵇美伟等[3]利用近景摄影测量技术提取了岩石节理信息，结合Python编写的两种算法实现了用非测量相机完成高精度信息的提取。激光三维扫描技术和近景摄影测量技术近些年被应用在各种工程，同时其缺点也暴露了出来，如有障碍物遮挡或者在险峻的地形下，这两种方法实现的过程就变得复杂。随着小型无人机的发展，其以快捷、高效地获取数据等优点被广泛应用于地质调查中。Goncalves等[4]采用无人机摄影测量方法，通过设置地面控制点，进行绘制和监测了沿海地区沙丘和海滩的情况，再结合三维重建软件对沙丘和海滩的地形进行模型重构，结果表明无人机可获取高精度的图像，构建的三维模型也符合实际。周成凯等[5]利用小型无人机对高位危岩进行了数据采集，通过采用赤平投影和图像识别相结合的方法对危岩稳定性进行评价，为地质信息获取提供了一种新方法。连会青等[6]利用无人机倾斜摄影测量技术获取地质数据，然后利用Context Capture软件进行三维重构，再把三维模型导入软件Acute 3D viewer中实现了对地质信息的识别提取，结果表明得到的数据和野外地质调查的数据一致性较好，具有实用意义。戴均豪等[7]利用小型无人机对甘肃北山地区进行数据采集，再利用Photoscan软件进行三维模型的合成，对三维模型进行解译标志和地质解译，得到了精细的地质图。

可见，小型无人机在地质调查中扮演的角色越来越重要了，同时也暴露出现阶段的研究存在着不足：①从航测照片到三维模型的重建，再到地质信息的识别与提取，其过程要经过多个软件的导入与导出，十分麻烦；②上述的三维模型重构软件多为商用软件，精度达不到岩土工程的要求；③地质信息识别与提取的软件较少，且不一定是符合工程需求。这些不足之处也是目前重要的研究方向，学者们做了很多尝试。杨帅等[8]引入多视图立体视觉算法，通过对无人机位姿数据的研究提出了一种无人机影像匹配的新算法，结果表明新的方法无论在查准率还是查全率上都有明显的替身。王楠[9]引入了区域生长法用于特征点的匹配，并对用于三维模型重构的随机抽样一致性算法进行了改进，减小了三维重构过程中的误差，使得点云集的匹配结果更加准确。Riquelme等[10]提出了基于MATLAB的结构面识别算法，采用了领域寻点的方法，搜索相同特性的点分为一组并对其进行去噪处理，可以高效找到岩体结构面并计算其产状。余宏明等[11]利用无人机摄影测量的方法，获取了测区高精度的航拍影像，通过空三计算得到了测区地形地貌的三维模型，以此为基础结合GeoStudio数值模拟软件，对其渗流场和稳定性进行分析，为计算降雨条件下滑坡的稳定性提供了一种快捷的方法。宋珺敏[12]以图像数据为基础，结合点云空间数据，得到以平面形态出露的结构面，再利用点云聚类的方法得到了岩体结构面节理信息，该方法的到的节理信息较符合实际。Liu等[13]利用无人机摄影技术采集边坡数据，并以此为基础建立三维模型，然后利用三维非连续变形分析软件模拟危岩的移动，最终得到边坡稳定性结果。金爱兵等[14]利用小型无人机获取边坡地形数据并建立三维模型，将三维模型与有限差分法软件FLAC3D相结合，得到边坡稳定性结果。

综上分析，现今学者的相关研究有两方面不足，首先大多数的研究主要是利用无人机摄影技术建立三维模上的研究和基于点云数据对岩体进行节理面信息提取算法的研究，而从航拍影像到节理面信息提取的整个过程的研究相对较少，即使有也大多利用商业软件辅助完成。其次，当测区存在大面积的植被时，三维建模的难度大大提升，有关采用无人机摄影测量技术对野外露头岩质边坡进行三维重构的研究相对较少。现针对以上问题运用无人机采用倾斜摄影测量技术对野外露头岩质边坡进行数据采集，基于增量式运动结构恢复算法和深度图融合的多视点立体视觉算法(multi-view stereo，MVS)对岩质边坡进行三维模型重构，然后再此基础上进行结构面的识别，提取节理信息。根据上述过程，开发设计基于倾斜摄影的裂隙岩体边坡三维模型重构系统，为无人机在边坡工程中的应用提供了完整方法。

1 倾斜摄影技术

1.1 倾斜摄影技术要点

摄影测量是通过采集、处理、提取和表达调查区域的影像信息来获取现场信息的方法。在传统摄影测量中，一般要通过测区的分布总图布置拍摄点，然后手工拍摄，这样不仅耗时还十分烦琐，并且传统的摄影测量主要采用的是正直影像，仅对地形地貌上部区域有较好识别的信息，而无法获取沿高度方向的地貌信息。倾斜摄影的出现，弥补了正直影像的不足，通过同一架无人机实现垂直、倾斜等不同角度的同时拍摄，进而获取更加准确的测区信息，其生成模型的精度是传统摄影测量精度的2～3倍[15]。在倾斜摄影时一般需要对下视、前视、后视、左视和右视5个方向进行拍摄，现采用单镜头小型无人机进行倾斜拍摄，无人机通过规划不同的航线进行采集数据，得到测区的完整场景信息。为了确保得到的场景信息能满足三维建模的要求，需要拍摄的影像之间具有较高的重叠度和一致的地面分辨率。其中重叠度可以通过规划无人机航线的旁向和航向的交叉大小来控制，因此地面分辨率是能否准确反应出地表真实信息的重要参数之一。正直影像地面分辨率的计算公式为

(1)

式(1)中：GSD为正直影像时的地面分辨率；δ为像元尺寸；h为无人机拍摄高度；f为相机焦距。再结合无人机相机主光轴旋转角度，可得到倾斜摄影时近点、中心的和远点的大地分辨率GSDj、GSDz和GSDy，公式为

(2)

(3)

(4)

式中：a和b分别为倾斜角度和视场角的一半，其中b的计算公式为

(5)

式(5)中：γ表示相机内部图像传感器的大小在一个方向的计算。

倾斜摄影的几何关系如图1所示，其中T′为影像中的一个点，T为T′在的实际场景中地面的对应点，PP′为影像中的主光轴与像平面的交点，β为T和倾斜摄影中心O的连线和ON竖线的夹角。从图1中可以看出，影响大地分辨率的因素主要包括无人机拍摄的高度、焦距、像素值、倾角等，其中焦距和像素值是由无人机携带的相机确定的；飞行高度越低，大地分辨率的值越小，则精度越高；结合式(4)和图1可知，当倾角越大时，远点的分辨率的数值越小且远点和近点分辨率差值越大。

1.2 数据采集

航测区域位于石家庄市井陉县南障城镇，周围群山环绕，地形复杂。测区长度为437.0 m，宽度为232.9 m，边坡高140.0 m。收集数据选用大疆御2专业版小型无人机，无人机及地面遥控器如图2所示，无人机飞行器配备先进的全方位视觉系统及红外传感系统，对周围障碍物有感应识别的功能，这样就使得无人机可以安全完成飞行任务。无人机上搭载高精度相机，并且配合高稳定性云台，足以满足工程的拍摄要求。地面遥控器采用高清图传技术，可以控制无人机在最大8 km通讯距离内完成飞行任务，无人机技术参数见表1。结合无人机操作端自带的 DJI Pilot软件，可以实现对飞行航线的规划，得到垂直影像和倾斜影像。本文中在测区总共设置五条航线，第一条航线镜头竖直向下，其他航线镜头倾角60°，飞行航线如图3所示，每种颜色代表一条完整的航线，倾斜摄影具体飞行轨迹是无人机从航线1起点出发，飞到指定高度完成航线1的任务，然后从航线1的最后位置飞向航线2的起点继续拍摄任务，直到所有航线任务完成回到航线1起点降落，由于边坡最高140 m，故飞行高度设置为200 m。在倾斜摄影测量过程中，为了能得到详尽的数据，航向重叠度设置为80%，旁向重叠度设置为70%，等距离拍照，具体航线规划设置和航线参数信息如表2所示。

图3 航线规划图Fig.3 Route plan

表1 无人机技术参数Table 1 Technical parameters of UAV

表2 航线参数Table 2 Route parameters

2 三维模型重建及节理信息识别

2.1 基本流程

由于小型无人机搭载相机相幅较小，且对旁向重叠度和航向重叠度要求很高，得到的航片照片数量较多，为了满足能同时识别较多的照片，本文中选择基于增量式运动结构恢复算法(structure-from-motion，SFM)[16]和深度图融合的MVS算法对岩质边坡进行三维模型重构，在此基础上进行节理信息数字识别。具体技术流程见图4。

图4 技术流程图Fig.4 Technology flow chart

在开源软件COLMAP的框架的基础上，对其进行二次开发，利用了COLMAP软件开源库里的增量式SFM算法、深度图融合的MVS算法和PCL软件开源库里的随机抽样一致算法(random sample consensus，RANSAC)、区域生长法，将其重新编译为既能进行三维模型重构又能进行岩体结构面识别提取的基于倾斜摄影的裂隙岩体边坡三维模型重构系统，实现了由倾斜摄影测量技术得到的数据生成三维模型和结构面自动识别的功能，软件界面如图5所示。

2.2 稀疏重构

采用增量运动恢复结构(增量式的SFM算法)进行稀疏重构[17]，该算法通过相机的运动同时恢复相机的参数和场景结构，其主要流程包括：特征点的提取、特征点的匹配和捆绑调整优化相机姿态和场景结构。

(1)特征点提取：小型无人机进行拍摄任务时会因不可控因素使得拍摄的照片产生畸变，畸变的影响会影响到三维重构的精度。选用SIFT(scale invariant feature transform)算法对畸变的照片进行处理，能较好地完成特征点的提取任务[18]。SIFT算法是通过在尺度空间中计算极值点来检测稳定的点的特征。其中尺度空间定义为

L(x,y,σ)=G(x,y,σ)I(x,y)

(6)

(7)

式中:(x,y)为点的坐标;σ为尺度大小；G(x,y,σ)为尺度σ的可变高斯函数；I(x,y)为原图像大小。SIFT点特征检测在高斯差分尺度空间中进行，目标函数为

D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]I(x,y)

(8)

式(8)中：k为比例系数。

利用高斯差分对数据进行计算，得到取得极大值与极小值时的三维点，这些三维点为离散的潜在特征点，然后将这些潜在特征点进行曲面拟合，拟合曲面的峰谷或者峰底即为最终特征点。通过本文研发的软件检测的特征点结果如图6所示(见图中红色点为特征点)。

图6 SIFT特征检测结果图Fig.6 SIFT feature detection result graph

(2)特征点匹配：由于采用无人机进行采集数据，数据中包含定位系统坐标(global positioning system，GPS)和飞行姿态，故选择空间匹配的方法。空间匹配的原理是首先给定一个查询点，然后利用K最近邻算法计算距离查询点距离最近的点记为候选点，最后判断候选点到查询点距离是否超过了最大距离，没超过最大距离的候选点加入特征点集，超过的最大距离的候选点剔除，完成特征匹配。

(3)增量的运动恢复结构：首先选取一对图像进行重建，再利用三角量测的方法对相机进行旋转、平移和三角化计算得到一个空间三维点，再把这个三维点重新投影到两个视角中计算重投影误差，误差较小则重建成功。其中重投影误差函数为

(9)

式(9)中：Pj为相机参数；Mi为三维点坐标；mij为第j幅影像上的第i个观测点；f(Pj,Mi)为观测点对应的投影点。

在此基础上增加新的图像重复上述过程，直至计算完所有照片，得到由空间三维点组成的稀疏场景模型，增量式SFM过程见图7。

图7 增量式SFM过程图Fig.7 Incremental SFM process diagram

完成上述所有计算即可得到稀疏三维点云模型，通过本文研发的软件得到的稀疏点云图如图8所示，图8中红色区域描述了无人机采样时的位置及相机角度。

2.3 稠密三维重建

由图8可以看出稀疏三维点云模型整体比较粗糙，不能直观地反映出岩质边坡的状态。故以稀疏三维点云模型为基础，采用基于深度图融合的方法，对岩质边坡进行稠密三维重建。

基于深度图融合的方法是通过对每幅影像重建深度图进而建立点云模型，由于不同影像拍摄的位置不同，看到的场景就会有所偏差，这些偏差即为视差，视差可以反映出物体在三维空间的深度。深度图融合的MVS算法的基本流程：首先，在航测照片集中选取航测的照片，对每一张航测的照片选取若干张照片作为邻近照片，用稀疏三维重建阶段的结果作为判断是否具有较高的相似性的标准，确保稠密三维重建的准确性。然后对初始的稀疏特征点开始深度估计并且把这些稀疏特征点作为种子点，用PatchMatch算法进行深度图估计，该算法是通过像素点q的像素坐标(qx,qy)和一个参数fq(aq,bq,cq)描述视差qd，公式为

qd=aqqx+bqqy+cp

(10)

PatchMatch算法中，每个像素初始化深度值是指用随机数为fq的三个参数设置初始值，然后把每个像素的深度信息传递到邻近像素，根据匹配代价判断是否更新邻近像素深度信息，最后不断循环迭代得到场景的深度信息。该过程中深度信息优化的数学模型为

(11)

式(11)中：IR为参考影像深度值；Ik为第k个视角的深度值；ck为第k个视角的深度尺度因子；(i,j)为位置坐标。

通过上述步骤，计算出每幅照片的深度图，如图9所示(不同颜色代表岩质边坡不同位置距离无人机的不同距离)。当计算出数据集中所有的图的深度图，把他们相互融合即可得到高程图，如图10所示，从高程图中可以看出，最高处的高程为140 m与实际测量值吻合。

图9 深度图及其原图Fig.9 Depth map and its original

图10 高程图Fig.10 Digital elevation map

最后，将深度图进行融合就能得到稠密三维模型及高程图。通过本文研发的软件得到稠密点云三维模型如图11所示，可以看出稠密三维重建的模型能较好地还原岩质边坡。

图11 稠密点云图和航拍全景图Fig.11 Dense point cloud and aerial panorama

2.4 节理信息识别

经上述计算，得到了岩质边坡三维模型，模型中包含着空间几何特征信息，因此可以对三维点云模型进行结构面的识别与提取相关信息。

结构面识别过程包括两个方面，分别为结构面的区域搜索和结构面的平面拟合。采用区域生长法[19]进行结构面区域搜索，由于输入的数据为稠密点云模型，其三维特征点的数目庞大，使得计算时间过长且效率低下，所以首先减少三维点云的数据且保持三维点云的形状特征不变，该过程就是对稠密点云模型进行下采样，下采样后把候选点中曲率最小的点作为初始生长点，搜索邻域，比较邻域中点的法向量和初始生长点法线之间的夹角与阈值的大小，小于阈值的点保留到聚类数据，再比较聚类数据中的点的曲率与曲率阈值的大小，小于阈值的点并入到初始生长点的区域，并将该点作为下一个生长点向外生长，直至到没有可以归并的点为止。

经过结构面的区域搜索后可以得到各个结构面的点集合，然后利用RANSAC(random sample consensus)进行岩体结构面拟合和转换[20]。在用RANSAC算法时，首先要选取有效的初始数据，在此基础上不断迭代增加数据集，最后选取数据最多的模型为最终模型。RANSAC平面拟合的流程如下。

(1)在三维点云数据集中选取初始点云数据，计算其对应平面方程为

Ax+By+Cz+D=0

(12)

式(12)中：A、B、C、D为平面方程参数。

(2)计算所有点至该平面的代数距离di，选取阈值d，若di

(3)重复以上步骤，根据数据集计算出迭代结束参数，根据该迭代结束参数判断迭代是否停止。

在开源软件COLMAP中加入了PCL编程库中Sample-consensus模块,实现了RANSAC平面拟合。Sample-consensus模块中模型定义为平面模型，模型设置了4个参数normal_x、normal_y、normal_z和d分别对应着平面方程参数A、B、C和D。本文中利用基于倾斜摄影的裂隙岩体边坡三维模型重构系统，得到的结构面识别图如图12所示。

图12 结构面识别图Fig.12 Structural plane identification map

由于三维模型的拓扑结构是基于无人机倾斜摄影拍摄的照片构建的，其坐标系和客观世界坐标系相同，故可以根据空间解析几何原理可以推导出结构面产状的关系[21]，推导关系如图13所示，表达式为

图13 结构面产状几何关系示意图Fig.13 Schematic diagram of the geometric relationship of structural planes

(13)

(14)

(15)

式中：θ和φ分别代表结构面的倾角和倾向;ω为同组结构面之间的间距。

从图13可以看出，该岩质边坡有多组结构面，利用研发的软件得到节理信息数据，对岩质边坡识别出125个结构面数据，给出部分展示数据如图5(b)所示。整理所有数据剔除错误结构面，得到7组显著结构面和若干个单独结构面，由于内容有限，仅提供显著结构面数据，平面方程参数见表3。为证明获得的岩体结构面倾向、倾角的可靠性，用人工现场测量的方法进行比较，如表4所示，可以看出对于7组结构面的产状，用两种方法计算的最大误差为3.19，具有一致性。根据表4中的结构面产状数据，将点云划分成有限区域，得到结构面划分图，展示了岩质边坡结构面的具体分布状况，如图14所示。

图14 结构面划分图Fig.14 Structural plane division diagram

表3 岩体结构面数据Table 3 Rock mass discontinuity data

表4 岩体结构面产状结果比对Table 4 Comparison of rock mass discontinuity occurrence results

3 边坡数值模拟计算

小型无人机对岩质边坡采集好数据，经过基于倾斜摄影的裂隙岩体边坡三维模型重构软件计算，得到地形网格模型及节理信息数据。在此基础上，利用FLAC3D有限差分法软件进行数值模拟。数值模拟中采用Hoek-Brown准则[22]，其表达式为

(16)

式(16)中：σ1、σ3分别为岩体最大主应力、最小主应力；σc为岩体抗压强度；mb、s、α为岩体的Hoek-Brown常量。

3.1 三维模型建立

基于倾斜摄影的裂隙岩体边坡三维模型重构软件可以提取点云的X、Y、Z各坐标数据，利用劳德内算法还原地形得到边坡地形图，以“ply”格式进行保存，如图15所示。利用Rhino软件对地形网格模型进行优化，生成“.f3grid”格式文件，然后把三维模型导入FLAC3D，三维有限元模型如图16所示，三维模型长400 m、宽400 m、高140 m。

图15 地形网格图Fig.15 Topographic grid

图16 三维有限元模型Fig.16 Three dimensional finite element model

3.2 三维模型参数

在进行边坡稳定性分析时，岩土体均采用Hoek-Brown本构模型，计算时除了模型上表面以外，对各边界均处采用固定端约束。式(16)中mb、s、α三个常量通常利用地质强度指标(geological strength index，GSI)对岩体参数进行确定的方法，其表达式为

(17)

(18)

(19)

式中：GSI为地质强度指标；D为岩体扰动因子，其取值范围为0.0～1.0，现场无扰动岩体为0，而非常扰动岩体为1.0；mi为岩体霍克布朗常数，可以通过查表5获得。研究区域为中等砂岩，现场无扰动岩体，查表5可知式(19)中mi=19，D=0。

表5 由岩石类型所决定的常数miTable 5 Constant mi determined by rock type

地质强度指标GSI 与岩体破裂程度和结构面的状态有关。从上述三维重构及节理信息数字识别得出的结果出发，选用岩体体积节理数JV和结构面条件因子JC对地质强度指标GSI进行计算。根据已有的节理数据推导出适用于本研发软件的岩体体积节理数的公式。

对于岩体体积节理数JV的统计，通常计算公式[22]为

(20)

式(20)要求每组结构面的数量两个及两个以上，才能算出结构面的间距。在实际测量中存在单个结构面，忽略它们将导致测量结果不准确，假设这些单个结构面的间距取各组结构面平均间距累加的平均值，在此基础上提出了计算岩体体积节理数JV的公式，即

(21)

式(21)中：N为单独的结构面个数；A为测区面积。

根据现场勘查可以知道本测区面积A=8 000 m2，结合表3和表4中的数据，利用式(21)，求得JV=25.62。结合Sonmez等[23]得出的地质强度指标表，如图17所示，地质强度指标表中可以看出结构面条件因子JC的大小根据结构面表面的状况分为1～18，JC越大表明结构面表面条件越好。根据如图11中的航拍实景图和稠密点云三维模型，结合采集数据时的现场勘查结果可知，该岩质边坡出露面风化严重且有大小不一的坑洞，表面条件较差，由此估计出JC=6。根据JV和JC的值查图17，得到GSI的值为52。为了验证结果的可靠性，根据Palmstrom的[24]公式进行计算，求得JV=25.60，可见JV值的大小仅差0.02，两者结果非常相近，证明了本公式的可靠性。

图17 GSI量化表Fig.17 GSI quantitative table

利用式(17)～式(19)可计算出Hoek-Brown准则中的参数mb、s和α，结果见表6。

表6 Hoek-Brown准则参数Table 6 Hoek-brown criterion parameters

3.3 数值计算结果

基于强度折减法对岩质边坡的安全系数进行计算，将岩质边坡土体的抗剪切强度参数逐渐减小，直到体系不平衡力与典型内力比率小于10-5，认为体系达到了力的平衡，即完成计算，利用model factor-of-safety命令即可实现该过程。在求解过程中FLAC3D软件能自动搜索出岩质边坡最危险的位置，分析其在临界状态下的稳定性。计算结果表明：图18为临界状态下岩质边坡总体位移云图，可以看出，岩质边坡最危险的位置在靠近山顶的陡坡处，最大位移为51 cm。岩质边坡失稳通常是沿着剪应变增量最大的位置发生，因此需要找寻岩质边坡剪应变的集中带，图19为临界状态下岩质边坡最大剪应变增量云图，可以看出，只有靠近坡脚处小部分的剪应变增量较大，其他处剪应变增量都较小且变化均匀，说明该边坡整体稳定性较好，但是在剪应变增量较大的区域容易发生剪切破坏，与现场勘察结果一致。