朱承金, 王述红, 任艺鹏, 孙子涵
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2019-03-26
基金项目：国家自然科学基金资助项目(51474050，U1602232)；中央高校基本科研业务费专项资金资助项目(N170108029)；辽宁省自然科学基金资助项目(20170540304，20170520341)；中建股份科技研发课题(CSCEC-2016-Z-20-8)。
作者简介：朱承金(1995-)，男，山东曲阜人，东北大学博士研究生;
王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。

摘要：为提高边坡结构面信息获取过程的安全性及测量结果的准确性, 采用无人机多角度摄影测量获取结构面信息, 得边坡结构面坐标以及出露迹线长.利用动态控制的相容Delaunay三角剖分对模型进行三角网格化, 结合块体生成规则对边坡临空面出露块体进行稳定性分析, 识别坡面出露关键块体, 并利用GeoSMA-3D系统验证方法的可行性及精度.研究结果表明:点云切割最大关键块体为38.574 8 m³, GeoSMA-3D切割最大关键块体为38.561 7 m³, 位置均在坡体中央位置, 拟合效果良好.对于不规则碎裂岩质边坡, 三维点云切割可以更好地表征边坡的块体位置, 对于指导关键块体治理具有一定优势.
关键词：无人机结构面块体岩质边坡Delaunay三角剖分
Slope Structural Plane Measurement and Block Stability Analysis Based on UAV Photogrammetry
ZHU Cheng-jin, WANG Shu-hong, REN Yi-peng, SUN Zi-han
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: SUN Zi-han, E-mail:wangshuhong@mail.neu.edu.cn
Abstract: In order to improve the safety of the information acquisition process of the slope structural plane and the accuracy of the measurement results, the multi-angle photogrammetry of the unmanned aerial vehicle(UAV) is used to obtain the structural plane information, and the coordinates of the slope structural plane and the length of the exposure are obtained. A dynamic Delaunay triangulation is proposed to triangulate the model. The stability analysis of the exposed surface of the slope is carried out in combination with the block generation rule. The key block is identified on the slope and the GeoSMA-3D is used to verify the feasibility and accuracy. The results show that the maximum key block of point cloud cutting is 38.574 8 m³, and the maximum key block of GeoSMA-3D cutting is 38.561 7 m³. Meanwhile, the blocks are both in the center of the slope. For irregular fractured rock slopes, 3D point cloud cutting can better characterize the block position of 3D slopes, and it has certain advantages for guiding key block management.
Key words: UAV(unmanned aerial vehicle)structural planeblockrock slopeDelaunay triangulation
岩质边坡的破坏不同于一般的均质、顺层土质边坡, 其稳定性主要受制于结构面的强弱, 且破坏具有一定的突发性.边坡结构面信息的精确获取是岩质边坡稳定性分析的前提.由于实际工程中某些边坡位置条件的特殊性, 使得传统近距离测量存在较大危险, 因此急需一种新型的测量方法来满足测量的需要.
三维激光和近景摄影测量技术近年来发展迅速, Li等^[1]基于三维激光技术对边坡进行扫描, 获取边坡的三维点云模型及结构面信息, 并对边坡结构面进行了聚类分组.Campbell等^[2]利用数字近景摄影测量以及投影变换等理论对结构面测量数据进行处理, 获取了结构面产状的计算模型, 完成了对迹线长的获取.王述红等^[3-4]基于非接触测量方法, 对岩体结构面信息进行获取, 并利用PFC颗粒流、非连续变形分析软件对模型进行受力分析, 得到了岩体损伤模型.
无人机航拍测量逐渐进入科学研究的领域, Niethammer等^[5]利用无人机对矿山边坡进行数据采集, 生成边坡三维点云, 完成了结构面的三维建模.Wang等^[6]基于多旋翼无人机得到岩质边坡的数字模型, 并对边坡模型进行辅助测绘, 得到了较好的效果.贾曙光等^[7]基于无人机多角度拍摄技术对边坡进行航测, 得到了边坡的三维点云模型, 并结合赤平投影对节理稳定性进行了判别.
本文采用无人机多角度摄影获取结构面信息, 用相应算法对边坡点云模型进行处理, 得到边坡结构面坐标以及出露迹长, 并按照产状因素进行分组, 采用动态控制的相容Delaunay三角剖分对模型进行三角网格化, 将结构面加入模型中, 并结合块体生成规则对边坡临空面出露块体进行稳定性分析, 对坡面出露关键块体进行识别, 并将分组结构面信息导入团队自行研发的GeoSMA-3D系统中, 导出关键块体信息, 将两种方法生成的块体信息进行对比, 验证方法的可行性及精度.
1 结构面信息获取1.1 现场多角度拍摄试验对象选取辽宁省本溪满族自治县本桓公路某岩质边坡, 该公路边坡坡体较为破碎, 节理裂隙发育, 且常有碎石掉落.传统近距离测量存在较大危险性, 故选用无人机摄影测量对其进行研究.图 1为航道平视拍摄获取的原始数据.
图 1(Fig. 1)

图 1 坡体航道规划Fig.1 Slope channel planning

1.2 影像识别重构及点云模型生成利用sfm运动算法进行计算, sfm^[8]是一种运动恢复结构算法, 可对每张图像中特征点进行提取, 分析图像序列的几何约束关系并重建目标的三维信息, 无人机在拍摄过程中会自动记录云点坐标的数据信息对图像进行匹配和重构, 提取每张图片中的特征点, 利用sfm算法进行整个坡体的稀疏点云生成.无人机摄影测量过程中, 根据稀疏点云坐标自带位置信息对稀疏点云坐标进行渲染, 见图 2.
图 2(Fig. 2)

图 2 稀疏点云渲染图Fig.2 Sparse point cloud rendering

采用多视觉立体视觉算法PMVS^[9](patch based multi-view stereopesis)对稀疏点云进行特征点提取和匹配, 通过膨胀扩散和插值, 得到物体表面的稠密点云集, 最后通过表面点云集的过滤和重构完成坡体密集三维点云模型, 对密集点云进行渲染并去噪.
1.3 结构面信息提取倾向倾角计算.由倾向倾角的定义可知, 岩体结构面产状的测定只需知道该产状所在平面的法向量即可, 假设该结构面平面方程为

(1)

式中:A, B, C为平面参数, 可得该平面法向量n为(-A, -B, 1).
对该结构面进行点的拾取任意不共线的n个点(n>3), 得到方程:

(2)

利用最小二乘法解算为

(3)

由转换公式可得

(4)

(5)

对结构面点云参数进行选点, 按照式(1)~(3)对点云信息平面方程进行拟合, 并按照式(4)、式(5)计算结构面产状, 导出结构面产状信息, 部分产状计算结果见表 1, 结构面的无人机摄影测量与人工测量结果对比见表 2.结果表明无人机摄影测量与人工测量结果误差均在0.85%以内, 满足计算要求, 进一步证明了该方法的可行性与精度.
表 1(Table 1)

表 1 部分产状计算结果Table 1 Part yield calculation results 结构面 结构面方程参数 产状信息 A B C 倾向/(°) 倾角/(°) 迹线长/m 01 0.368 0.245 0.621 33.7 35.5 6.55 02 0.440 -1.270 0.265 109.0 78.9 6.70 03 0.337 2.824 0.324 263.2 83.5 8.53 04 -1.380 3.452 0.443 291.8 83.2 4.26 05 0.455 2.088 0.266 282.3 82.9 2.88 06 -0.100 1.480 0.185 265.9 82.9 1.20 07 -0.640 5.600 0.247 96.6 73.2 5.88 08 0.316 0.882 0.201 250.3 77.9 1.47 09 0.447 1.703 0.247 284.7 89.0 2.06 10 57.150 -28.400 0.433 26.5 89.6 1.95

表 1 部分产状计算结果 Table 1 Part yield calculation results

表 2(Table 2)

表 2 结构面的无人机摄影测量与人工测量结果Table 2 Results of UAV photogrammetry and manual measurement on structural plane 结构面 无人机摄影测量产状信息 人工测量产状信息 对应误差/% 倾向/(°) 倾角/(°) 迹线长/m 倾向/(°) 倾角/(°) 迹线长/m 倾向 倾角 迹线长 01 33.7 35.5 6.55 33.9 35.8 6.55 0.59 0.84 0.00 02 109.0 78.9 6.70 108.7 79.1 6.68 0.28 0.25 0.30 03 263.2 83.5 8.53 263.2 83.8 8.62 0.00 0.36 0.10 04 291.8 83.2 4.26 291.5 83.4 4.25 0.10 0.24 0.26 05 282.3 82.9 2.88 282.4 82.7 2.89 0.04 0.24 0.35 06 265.9 82.9 1.20 266.0 82.9 1.20 0.04 0.00 0.00 07 96.6 73.2 5.88 96.9 73.1 5.90 0.31 0.14 0.34 08 250.3 77.9 1.47 251.0 78.0 1.46 0.28 0.13 0.68 09 284.7 89.0 2.06 285.0 89.2 2.07 0.11 0.22 0.48 10 26.5 89.6 1.95 26.5 89.2 1.95 0.00 0.45 0.00

表 2 结构面的无人机摄影测量与人工测量结果 Table 2 Results of UAV photogrammetry and manual measurement on structural plane

2 坡体出露结构面块体识别结构面与坡体表面相互贯通和切割, 可以得到坡体表面的且存在失稳可能的块体.采用Delaunay三角剖分^[10], 对三维点云模型进行三角网格划分, 并与拟合结构面进行接触判别和切割.
2.1 Delaunay剖分原理Delaunay三角剖分中常采用的是相容Delaunay方法, 利用添加新节点的方式来提高网格的质量, 但同时也会产生节点数量过多的副作用, 导致数据量臃肿, 运算效率降低.故提出动态控制的相容Delaunay三角剖分, 在增加节点的过程中通过设定阈值来动态控制节点数, 既要求网格具有一定质量, 同时又兼顾运算效率.
按动态控制的相容Delaunay三角剖分对坡体进行三角剖分, 采用在丢失线段及刻面上增加新节点的手段, 使得已丢失的线段及刻面在最后的剖分网格中显现出来.并在增加新节点的过程中动态控制其数量, 避免影响运算效率, 该方法在一定程度上克服了线段及刻面丢失的弊端.
2.2 块体切割与识别采用Wang等^[11]所阐述的方法对边坡进行切割.将拟合结构面信息导入点云模型中, 面与体交错, 生成大量的凸形及凹形块体, 需要注意的是凹形块体在程序运算过程中极难控制与识别, 故有必要将所有的凹形块体用若干个凸形块体的组合来表达.
首先将导入的结构面定义为无限大, 即结构面完全切割边坡模型, 此时边坡模型中均为凸形块体, 不存在凹形块体, 满足计算高效化的要求.
边坡体被完全切割, 但实际结构面是有固定大小的, 故需要将为方便程序计算而设置的虚拟裂隙剔除.被剔除的虚拟裂隙需要满足以下原则:该裂隙为其两侧块体的共有部分, 且两侧块体存在交集, 并且该交集不属于上述裂隙.两个块体间的虚拟裂隙被剔除后则形成一个完整块体.
通过对边坡模型完全切割及虚拟裂隙剔除, 原模型可被切割为大量的简单凸形及复杂凹形块体, 即实际情况下边坡模型中存在的独立块体.至此, 块体切割结束.
按照上述方法对边坡坡面进行切割, 并按照极限平衡法则对块体稳定性进行判别, 其中同时满足“存在出露面、几何可移动、安全系数小于规定值”三个条件的块体称为关键块体.
1) 存在出露面判定.结构面与边坡坡面相切割形成的块体均满足存在出露面要求.
2) 几何可移动判定.判定式为

(6)

式中:w是块体所受驱动力; n_i是块体表面的方向矢量(指向块体内部); i是多边形在块体内的编号(出露面不计).
3) 安全系数小于规定值判定.根据Mohr-Coulomb准则, 块体在其自重力作用下的安全系数计算式为

(7)

式中:G为块体自重; α_i为滑动面的倾角; S_i为滑动面i的面积; c_i, φ_i分别为滑动面上的黏聚力和内摩擦角.关键块体分布位置如图 3所示.
图 3(Fig. 3)

图 3 块体识别结果Fig.3 Block recognition results

图中具有滑落危险的块体, 由深色标识标出, 共有11块, 分布在坡体出露面.块体安全系数等参数见表 3.
表 3(Table 3)

表 3 点云关键块体计算结果Table 3 Calculation results of point cloud key block 块体编号 滑动面数目 块体体积/m3 安全系数 01 2 15.778 4 0.967 02 3 38.574 8 0.882 03 1 5.556 7 0.969 04 2 7.413 4 0.858 05 2 9.503 7 0.962 06 1 23.809 4 1.136 07 2 16.144 8 0.945 08 2 12.510 0 0.923 09 1 27.237 5 0.984 10 2 6.571 5 0.953 11 2 25.169 1 0.991

表 3 点云关键块体计算结果 Table 3 Calculation results of point cloud key block

将结构面产状, 出露迹线长等参数输入到自行研发软件GeoSMA-3D中, 采用结构面Beacher圆盘对坡体进行切割, 按照该坡体的情况对三维边坡进行建模.切割块体、结构面迹线切割计算结果见图 4，图 5.
图 4(Fig. 4)

图 4 切割块体结果Fig.4 Cutting block result

图 5(Fig. 5)

图 5 结构面迹线切割结果Fig.5 Cutting results of structural plane trace

表 4为GeoSMA-3D关键块体计算结果.对比两种计算方法, 点云切割得到11组块体, 其中滑动面组成为三滑动面的块体有1组, 双滑动面块体7组, 单滑动面块体3组.具有滑动危险的块体为10组.GeoSMA-3D系统进行切割, 得到10组块体, 单滑动面块体6组, 双滑动面块体4组, 具有滑动危险的块体9组.
表 4(Table 4)

表 4 GeoSMA关键块体计算结果Table 4 GeoSMA key block calculation results 块体编号 滑动面数目 块体体积/m3 安全系数 01 2 5.421 6 0.913 02 1 38.561 7 0.856 03 1 7.146 9 0.932 04 2 8.167 4 0.978 05 1 22.871 6 1.017 06 1 24.598 3 0.953 07 1 8.669 1 0.959 08 2 25.237 5 0.906 09 1 4.871 6 0.987 10 2 4.183 7 0.975

表 4 GeoSMA关键块体计算结果 Table 4 GeoSMA key block calculation results

对识别块体体积进行统计, 结果见图 6.由表 2，表 3和图 6可以看出, 两种方法切割块体体积在0到40 m³之间, 点云切割最大块体为38.574 8 m³, 安全系数为0.882, GeoSMA-3D切割最大块体为38.561 7 m³, 安全系数为0.856, 二者均分布在坡体中央位置.
图 6(Fig. 6)

图 6 块体体积统计结果Fig.6 Statistical results of block volume

两种切割方法产生的块体具有一定的差异.比较两种方法可以得到, 传统的GeoSMA-3D与点云切割模型相比, GeoSMA-3D无法表征三维复杂边坡, 而三维点可以较好地表征三维复杂边坡块体的位置, 尤其对于不规则碎裂岩质边坡, GeoSMA-3D块体预判位置具有不确定性, 对提前支护造成一定困难, 且GeoSMA-3D边坡切割块体为不规则块体, 相应支护措施的选取较为复杂, 而点云切割则避免了该问题的产生, 模型可以更好表征三维边坡的块体位置, 对于指导关键块体治理具有一定优势.
3 结论1) 采用无人机对边坡结构面信息进行采集, 并进行边坡点云模型的构建, 通过结构面点云坐标对结构面进行拟合, 得到边坡结构面信息, 相比于传统测量方法, 不受时空和位置的限制, 测量危险性大大减小.
2) 点云切割模型与GeoSMA-3D切割模型得出的块体体积及块体切割位置大致相符, 为关键块体滑落预测及边坡支护提供了理论支持.
3) 对于不规则碎裂岩质边坡, 三维点云切割可以更好表征三维边坡的块体位置, 对于指导关键块体治理具有一定优势.
参考文献

[1]	Li H B, Li X W, Li W Z, et al. Quantitative assessment for the rockfall hazard in a post-earthquake high rock slope using terrestrial laser scanning[J].Engineering Geology, 2019, 248(8): 1–13.
[2]	Stephanie C, Robert S, Jane R, et al. Using near-surface photogrammetry assessment of surface roughness(NSPAS)to assess the effectiveness of erosion control treatments applied to slope forming materials from a mine site in West Africa[J].Geomorphology, 2018, 322(1): 188–195.
[3]	王述红, 昝世明, 王存根, 等. 基于非接触测量技术裂隙岩体隧道开挖局部损伤分析[J].现代隧道技术, 2014, 51(3): 110–116. ( Wang Shu-hong, Zan Shi-ming, Wang Cun-gen, et al. Analysis of local damage induced by tunnelling in a fractured rock mass based on noncontact measurement Techniques[J].Modern Tunnelling Technology, 2014, 51(3): 110–116.)
[4]	王述红, 朱承金, 张紫杉, 等. 基于动态强度折减DDA法的边坡多滑面稳定性分析[J].煤炭学报, 2019, 44(4): 1084–1091. ( Wang Shu-hong, Zhu Cheng-jin, Zhang Zi-shan, et al. Stability analysis of multi-slip surface of slope based on dynamic strength reduction DDA method[J].Journal of China Coal Society, 2019, 44(4): 1084–1091.)
[5]	Niethammer U, James M R, Rothmund S, et al. UAV-based remote sensing of the Super-Sauze landslide:evaluation and results[J].Engineering Geology, 2012, 128(4): 1–11.
[6]	Wang S H, Zhang Z S, Wang C G, et al. Multistep rocky slope stability analysis based on unmanned aerial vehicle photogrammetry[J].Environmental Earth Sciences, 2019, 78(8): 260–275.DOI:10.1007/s12665-019-8145-z
[7]	贾曙光, 金爱兵, 赵怡晴. 无人机摄影测量在高陡边坡地质调查中的应用[J].岩土力学, 2018, 39(3): 1130–1135. ( Jia Shu-guang, Jin Ai-bing, Zhao Yi-qing. Application of UAV oblique photogrammetry in the field of geology survey at the high and steep slope[J].Rock and Soil Mechanics, 2018, 39(3): 1130–1135.)
[8]	Lindeberg T. Feature detection with automatic scale selection[J].International Journal of Computer Vision, 2008, 80(2): 189–210.DOI:10.1007/s11263-007-0107-3
[9]	Shi L M, Guo F S, Hu Z Y. An improved PMVS through scene geometric information[J].Acta Automatica Sinica, 2011, 37(5): 560–568.
[10]	Lee D T, Lin A K. Generalized Delaunay triangulation for planar graphs[J].Discrete & Computational Geometry, 1986, 1(3): 201–217.
[11]	Wang S H, Ni P P. Application of block theory modeling on spatial block topological identification to rock slope stability analysis[J].International Journal of Computational Methods, 2014, 11(1): 903–914.