程雨柯，李亚虎，夏金梧，侯 赠，陈 娜
（1.长江勘测规划设计研究有限责任公司，湖北 武汉 430011；
2.水利部长江勘测技术研究所，湖北 武汉 430011；
3.湖北工业大学土木建筑与环境学院，湖北 武汉 430068）
高陡边坡上的危岩区域是水利水电工程建设中常见工程问题之一，其稳定状况事关水利枢纽工程建设和运行安全[1-2]。对边坡岩体进行地质勘察，获取边坡岩体的地质参数，是危岩体研究的首要工作。传统的地质勘察主要采用的是人工接触式方法，例如用皮尺和地质罗盘测量[3-5]。随着建设工程规模的不断增大，高陡边坡地质环境变得更加险峻和复杂，由于调查人员的可及范围有限，地质勘察工作的展开不仅效率低下，而且获取的地质参数不完整。传统的接触式调查已经难以满足实际工程的高效化和数字化要求。
在过去的20 年中，许多研究人员采用了非接触式数据采集方法来解决这一问题[6]。Bitenc 等[7]将优化后的激光雷达测量技术（light laser detection and ranging，LiDAR）应用于边坡地形测量数据采集，提高了边坡数据采集的精度；
Salvini 等[8]结合直升机摄影测量技术和地面激光扫描技术（terrestrial laser scanning，TLS）进行了边坡岩石块体和节理系统的详细地理信息绘制；
Kumhálová等[9]，贾虎军等[10]通过机载激光扫描和倾斜摄影技术，获取精准的地质灾害隐患点、精细的地形地貌勘测数据，解决了设备难以进入的陡峭地形模型数据采集的问题。然而，激光扫描设备即便借助直升机等多地形工具，在进行高差过大，时限较短的数据采集任务时仍然显得笨重，难以满足需求。
近年来，为了满足快速获取复杂边坡表面特征的需求，成本更低、灵活方便、不受地形限制的无人机得到了一定的应用，它可以解决复杂边坡死角模型难以采集的问题。Niethammer 等[11]，Riquelme 等[12]，Zeybek 等[13]采用无线电控制小型四旋翼无人机拍摄完整边坡滑坡的图像，进一步建立高分辨率的数字地形模型，实现了更快速获取坡面点云数据；
赵婷婷等[14]，梁京涛等[15]，黄发明等[16]通过无人机倾斜摄影技术（unmanned aerial vehicle oblique photography technology，UAVOP）与多图像重建技术（structure from motion，SfM）采集复杂边坡表面点云模型，实现了对紧急，复杂的边坡数据采集任务的快速反应。基于无人机的数字摄影测量技术已经足以表现复杂边坡表面的详细形态，但这些研究对精确模型的利用率仍然较低，尚未充分地通过精确点云对其所反映的边坡隐蔽、细微的安全隐患（如危岩体）进行有效的勘察。
采用无人机开展高陡边坡的安全勘察研究时，由于数据采集方案单一，设备效率、模型精度较低或高精度模型的利用率较低[17-19]，尤其在高陡边坡中，由于畸变问题，会使得模型结果与实际情况产生一些偏差。为解决这一问题，本研究提出了一种基于无人机点云的超500 m 高位边坡孤立单平面滑移式危岩体半自动提取的新方法，实现了以下三项关键的目标:
（1）基于免像控的无人机贴近摄影测量技术，设计了一套高陡边坡表面点云模型的测量方法，并构建了三维坐标精度为±10 cm、影像分辨率为1 cm 的试验区精细边坡三维点云模型。
（2）基于高陡边坡点云模型提取了异型滑移式块体。利用异型滑移式块体突出于边坡表面的细微特征，手动分割确定了异型滑移式块体的位置和其具体形态。
（3）根据异型滑移式块体模型边界点云的密度特征拟合了后壁平面。此外，自动计算了异型滑移式块体模型的各类特征值，并进一步得出稳定性系数（K）以确定块体中的危岩体。最后，结合危岩体的方量和最大高差，对危险岩体进行了危害性评价。
玉龙喀什水利枢纽工程（图1）是新疆和田玉龙喀什河山区河段的控制性水利枢纽工程，水库设计最大坝高230.5 m，总库容5.36×108m3，工程主要有最大坝高233.5 m 的面板堆石坝、“表、中、深”泄洪隧洞以及引水发电系统等组成，属大（2）型Ⅱ等工程。
图1 试验区三维地形概览Fig.1 Three-dimensional terrain overview of the test area
试验区为工程左岸坝址区的高位边坡，最大高差超过500 m。边坡地形陡峭，坡面植被不发育，山体的表面覆盖有大量灰尘，边坡表面并未出露明显的结构面。考虑施工场地狭窄、施工周期长，坝基边坡开挖对上部高位边坡扰动大，而坝址区两侧自然高位边坡存在大量危岩体和不良地质块体，在施工期构成重大安全隐患。亟须对坝址区超高位边坡进行危岩调查和稳定性分析。
经过人工初步调查，勘察到可能的危岩体多为平面滑移式危岩体（见图2）。岩体突出于边坡表面，且与边坡斜面之间存在有明显的连接边缘。本研究对试验区表面此类孤立的单平面滑移式危岩体进行进一步的调查。
图2 人工对左岸边坡的典型危岩体的初步调查Fig.2 Preliminary investigation of typical dangerous rock masses on the left bank slope by manual methods
针对试验区高边坡广泛分布平面滑移式危岩体的情况，本文设计了一套从边坡点云模型数据采集到危岩体半自动识别及评价的边坡危岩体完整的地质勘察方法（图3）。程序主要由Python 语言编写。该识别方法包括以下四个步骤：
图3 研究技术方法Fig.3 Research technical scheme
（1）数据采集阶段：采用无人机贴近摄影测量结合精细的仿地飞行航路设计、智能建模等多种技术手段，获取高边坡精细地表LAS 点云。
（2）孤立异型滑移式块体点云模型提取阶段：依据边坡突出块体边界特征将孤立的异型滑移式块体点云模型分别手动分割出来，并对模型进行点云下采样。
（3）滑移破坏后壁平面提取阶段：应用点云的密度特征自动提取突出块体与边坡连接的边界，此边界可近似为该块体的滑移破坏的平面的边界。
（4）孤立危岩体定性与危害性评估阶段：应用最小二乘法将边界点云拟合到后壁平面上，计算后壁倾角[20]。使用原块体点云和对后壁投影点之间的欧式距离计算最大高度差。利用Stokes 公式和Alpha Shape 算法确定投影点的平面边界，计算后壁平面的面积和块体点云体积。以上特征参数将用于定性孤立的危岩体（依据稳定性系数）。对于已经定性的危岩体，依据其体积与最大高差依次对所有危岩体的危害性进行进一步的评价。
2.1 数据采集及异型滑移式块体识别
本文使用的数据是新疆玉龙喀什水利工程左岸斜坡的有序点云数据，考虑到试验区山、谷、沟、梁交错的特殊地形，本文采用无人机贴近摄影测量[21-23]和全球导航卫星系统（global navigation satellite system,GNSS）实时动态差分定位技术[24]相结合的方法，综合利用免像控的无人机摄影测量技术、精细化地面模拟航路设计、智能建模等多种技术手段，获取500 m 以上边坡的精细化地表信息。
本研究使用大疆M300 RTK 无人机设备搭载Zenmuse P1 全画幅云台相机（见图4）。无人机集成了双GNSS导航定位设备，通过集成的网络RTK 模块，获取附近CORS 站提供的实时动态差分信息，可为无人机提供实时厘米级定位数据，满足了免像控标准。
图4 无人机设备Fig.4 UAV equipment
在对试验区模型采集的应用中，本研究设计的基于精细仿地飞行路线的贴近摄影测量技术方法主要采用三种综合航摄方案。
（1）初始飞行和航线规划（见图5）：采用航天飞机雷达地形测绘（shuttle radar topography mission,SRTM）采集30 m 分辨率数字高程模型（digital elevation model,DEM），初始飞行航路较高。通过空中三角测量，可以获得调查区域的图像，并生成初始数字地表模型（digital surface model,DSM）[25]。
图5 第一张正射影照片路线规划（分辨率2.5 cm）Fig.5 Route planning of the first orthophoto photograph (resolution 2.5 cm)
（2）二次精细仿地飞行规划与飞行（见图6）：基于初始模型DSM 高程数据，规划精细地面模拟飞行路线，并为该路线建立多组使用精细仿地飞行的镜头拍摄角度。为了获得高坡度倾斜摄影测量图像，对近似直立的坡面使用了贴近摄影测量技术。通过无人机三角测量和建模，生成测量区倾斜摄影测量OSGB 模型[26]。考虑到试验区边坡沟梁相间的地形，精细的仿地飞行路线沿地形梯度呈上下方向设置。
图6 精细的仿地飞行路线规划Fig.6 Detailed terrain-following flight path planning
（3）三次手控飞行：检查OSGB 模型存在的空洞和纹理拉伸的区域。如果出现这种情况，必须对关键区域进行额外的手动飞行拍摄。补拍完成后，再次进行空中三角测量和建模，生成测量区域完整的真实场景，即为最终的三维模型。针对研究区倒悬岩体等特殊场景，采用-10°等特殊视角和手控飞行的拍摄方案。
DJI Maps 被用于无人机三维模型建模，采用基于计算统一设备架构（compute unified device architecture,CUDA）的智能重构算法[27]，并将场景图文件格式模型（open scene graph binary,OSGB）转换为二进制LAS 点云模型（见图7）。
图7 无人机三维模型建模并生成点云模型Fig.7 UAV three-dimensional modeling and point cloud model generation
异型滑移式块体在边坡表面表现为一块完整的突出块体。对于边坡的点云模型，其与边坡接触处通常为一圈曲率较大的闭合点云。依据此特征在软件Cloud-Compare 内使用点云切分工具将孤立的异型滑移式块体点云模型进行提取。此外，为了减小后续的计算量，使用滤波工具进一步进行点云下采样。
2.2 块体失稳平面自动提取
2.2.1 后壁平面边界点云
本研究认为，异型滑移式块体点云的边界近似处于滑移的后壁平面上，依据边界点云的密度与密度变化率特征进行后壁平面边界的提取。对于表面均匀的点云模型来说，处于模型边缘的点云一般密度（在某范围内的点云个数）相对较小[28-30]。此外，该处的点云在某个范围内密度会明显的下降。
对于单独的异型滑移式块体点云，构造点云的K邻近（K-nearest neighbor,KNN）搜索结构，对每个点查找其搜索范围内的点云密度。对于点P1，当该点点云密度较小，在R1 范围内的点数为8（该点的密度）或者在一定的范围R2 内，存在密度远超于它的点P2（密度为19）时，本研究认为其为一枚后壁平面的近似边界点（见图8）。
图8 边界点云密度与密度变化率Fig.8 Boundary point cloud density and density change rate
2.2.2 后壁平面拟合
异型滑移式块体边界点云可近似拟合于空间中的一块平面。考虑到试验区表面被大量的灰尘覆盖，边坡模型并未出露明显的结构面，则该平面被认为是当前异型滑移式块体可能的滑动结构面（后壁平面）。图9（a）中蓝色的点云为通过边界密度分析得到的后壁平面边界点，灰色点为该块体的其他点云。利用最小二乘法将边界点拟合为图9（b）中的黑色网格平面，该平面将为异型滑移式块体点云模型的各种特征计算提供依据。
图9 异型滑移式块体的后壁平面提取Fig.9 Extraction of the rear wall plane of profiled block
2.3 孤立危岩体定性分析
2.3.1 三维特征估算
将异型滑移式块体点云整体向后壁平面进行投影，并对其进一步进行Alpha Shape 边界搜索[31]，得到平面的边界点索引。该索引方向对应的边界点逆时针排列。对边界点利用Stokes 的三维有向平面面积积分公式进行积分[32]，得到后壁平面的面积。后壁平面多边形A有N个顶点P1,P2,···,Pi,···,Pn呈逆时针方向排列，形成一个环L顶点坐标表示为Pi=(Xi,Yi,Zi)，如图10（a）所示。可求得其面积为SA=。
图10 异型滑移式块体的后壁面积与块体体积计算Fig.10 Calculation for the rear wall area and block volume of profiled block
面积由Stokes 公式求得：
式中：L——后壁平面多边形边界形成的环；
P、Q、R——代表空间中某一点依赖于x、y、z坐标的向量场三个分量函数；
dx、dy、dz——表示线积分中的微分元素；
A——后壁平面多边形；
(n,x)、(n,y)、(n,z)——多边形的单位法向量n与坐标轴x、y、z的夹角；
ds——平面上的微小面积元素。
对完整的异型滑移式块体点云应用Alpha Shape 算法得到异型滑移式块体点云的体积，如图10（b）所示。后壁平面与水平面之间的夹角为后壁倾角[33]，后壁倾角由公式（2）计算，可记为 θ，如图11（a）所示。后壁面与孤立岩体所有点之间的最大垂直距离为最大高差，见图11（b）。孤立岩体的最大高差可依据公式（3）进行计算。
图11 异型滑移式块体的后壁倾角与最大高差计算Fig.11 Calculation of the rear wall inclination angle and maximum height difference of profiled blocks
式中：θ——后壁面倾角；
式中：∆H——孤立岩体最大高差；
Ai——第i个对应的投影点云；
Bi——孤立岩体的第i个原点云。
2.3.2 危岩体定性
得到了岩体点云模型的各种空间特征后，要判断提取的异型滑移式块体是否为危险岩体[34-35]，则应结合稳定系数（K）进行定性判断。稳定系数是定性判断沿后壁平面滑移破坏的危岩体的常用依据[36]。其公式已广泛应用于单滑移面危岩体稳定性试验领域。利用《中华人民共和国新能源行业标准》中《水电工程危险岩体工程地质调查与治理规范》（以下简称《规范》）对K进行计算和评价，如表1。
表1 单平面滑动岩体稳定性评价Table 1 Stability evaluation of single-plane sliding rock mass
对于每一个可能的危险岩体对象，自动依次计算其对应的后壁面积、后壁倾角和块体方量，并对结果进行综合的稳定性分析。稳定性系数可依据式（4）计算。
式中：Q——异型滑移式块体的重量/N；
α——后壁的倾角/（°）；
φ——后壁滑动面的内摩擦角/（°）；
c——黏聚力/MPa；
S——异型滑移式块体的滑动面积/m2。
2.4 危岩体危害性评估
应用安全系数判断异型滑移式块体是否为危岩体后，还需要综合评价危岩体方量和最大高差的特征[37]，以确定这些危岩体的危险等级和危害性，其中《规范》提供的危岩体规模分类依据如表2 所示。
表2 水电工程危险岩体规模分级Table 2 Scale classification of dangerous rock mass in hydropower projects
3.1 数据采集与异型滑移式块体识别
通过精细化地面模拟航路设计、无人机贴近摄影测量、多图像三维智能建模等步骤重建获得玉龙喀什水利工程左岸500 m 以上边坡的高精度三维点云模型信息（图12）。考虑到玉龙喀什水利枢纽孤立岩体多为异型滑移式块体，研究以左岸1 块危岩体密集的区域进行危岩体半自动提取试验。
图12 左岸整体点云，左岸样本区区域Fig.12 Overall point cloud on the left bank and sample area on the left bank
在软件CloudCompare 内依据异型滑移式块体与边坡面的接触边界特征将该样本区可能的异型滑移式块体点云模型切分出来，如图13（a）。
图13 异型滑移式块体点云模型提取Fig.13 Point cloud model extraction of profiled blocks
共提取出5 块异型滑移式块体（L1，L2，L3，L4，L5）见图13（b），所有的块体均具备明显的边界特征。此外，为了减少计算成本，对点云进行下采样。图14（b）为块体L1 下采样后的点云（采样前点云个数:10 486 631，采样后点云个数：1 171 753）。
图14 异型滑移式块体下采样Fig.14 Downsampling of profiled block
3.2 精确特征提取与危岩体定性
以块体L1 为例，对于切分出的异型滑移式块体点云模型，使用点云密度与密度变化特征提取后壁平面边缘点云见图15（b），这类点可近似处地于后壁平面上。
图15 依据点云密度提取的异型滑移式块体后壁平面边缘点云Fig.15 Edge point cloud of the rear wall plane of profiled blocks extracted based on point cloud density
图15（b）中蓝色点为通过点云密度分析所提取的后壁平面边缘点云，其表现为一圈闭合的边界。而灰色点为块体点云的其他部分。五块异型滑移式块体点云边界准确模型见图15（a）。利用最小二乘法将点云拟合为平面方程，此平面为岩体破坏时可能的滑动面（见图16）。
图16 最小二乘法的后壁平面拟合Fig.16 Rear wall plane fitting using the least square method
根据《规范》，应确定异形滑块的后壁倾斜度、后壁面积与方量三项特征，以判断和评价危岩体。以异型滑移式块体L1 为例，对L1 点云边界应用最小二乘法得到平面方程，由式（4）计算后壁面与水平面的法向量夹角，即后壁面倾角。将L1 点云垂直投影到后壁上，计算原始点与投影点之间的欧氏距离。最大距离对应于L1 的最大高差（见图17）。
图17 异型滑移式块体L1 的后壁平面和点云投影Fig.17 Rear wall plane and point cloud projection of profiled block L1
L1 的面积和体积可以使用Alpha Shape 算法和Stokes 积分计算。图18（a）中的红色点云表示通过Alpha Shape 算法提取的异形滑动块的后壁有向的精细边界。该边界将用于计算后壁面积。此外，图18（b）描绘了异形滑块外表面的三角形包络网络，将用于计算块体的体积。
图18 Alpha Shape 算法分别计算异型滑移式块体的后壁面积与体积Fig.18 Alpha Shape algorithm for separate calculation of the rear wall area and volume of profiled block
如表3 所示，采用提出的方法确定各异形滑块的后壁倾角、后壁面积、岩石体积和最大高差。
表3 异形滑移式块体特征数据统计Table 3 Statistical analysis of characteristic data for profiled blocks
通过工作人员的现场调查，可获得现场大部分岩质边坡危岩体以及后壁平面的物理力学参数（如表4）。进一步，对所有识别出的异形滑块进行稳定性K计算，从而判断其是否为危险岩体（如表5），可以看出，除了L5 以外均为危岩体。
表4 异形滑移式块体物理力学参数Table 4 Physical and mechanical parameters of profiled block
表5 危险岩体稳定系数的计算与定性Table 5 Calculation and qualitative characterization of stability coefficient of dangerous rock mass
3.3 危岩体综合评价
在本研究中，在提取了试验区的危岩体后，还需要对危险岩体的危险性等级进行确定和评价。孤立岩体的稳定性和危险等级应结合稳定性系数、岩体体积和最大高差来确定。为了定量评价危险岩体的危险性，结合危险岩体的体积和高差特征对《规范》进行了评价。
如表6 所示，本研究在试验区发现了一系列小山梁状危岩体，体积均在2 000 m3左右，属大型的单平面滑动危岩体，危岩体最大高差在7～15 m，其中L4 高差为14.53 m。经过现场调查，这些危险岩体已被爆破拆除。
表6 危险岩体与岩体体积评价Table 6 dangerous rock mass and rock mass volume evaluation
3.4 适用性讨论
（1）方案没有考虑详细的地质参数。由于本研究的重点是快速识别危险岩体并提供初步评价，因此没有将岩体的具体地质特征作为依据。在未来的研究中，将自动识别过程与详细的地质参数相结合，实现对危险岩体更精细的评价。
（2）本研究只考虑了危岩体的单平面滑动的破坏模式，没有对其他滑动模型进行研究。这主要是由于研究区内普遍存在沿单一平面滑动的岩体。目前提出的半自动识别模型仅能提取沿后壁平面滑动的岩体。在未来的研究中，该模型将进一步扩展，以包含识别具有各种失稳和破坏模式的危岩体。
本研究基于无人机贴近摄影测量技术与点云模型数据处理算法，对玉龙喀什水利枢纽工程左岸边坡一块岩体密集区进行了初步的快速地质勘察，实现了试验区边坡点云模型重建、异型滑移式块体提取、稳定性分析等目标。该方法在快速获取高精度边坡数据、智能提取三维异型滑移式块体、准确计算三维特征、评估危险岩体等方面相对于传统方法具有一定的优势。本文提供了从数字边坡模型获取到危险岩体识别和初步评价的解决方法。为边坡工程数字化奠定一定的研究基础。研究得出了如下的主要结论:
（1）结合无人机贴近摄影测量、GNSS 实时动态差分定位技术，以及自主设计的三次精细综合无人机飞行路线规划和手控校验，可为航空三角测量和建模提供分辨率优于1 cm 的图像，解决了由于高差较大而产生模型畸变的问题。
（2）利用最小二乘法和Alpha Shape 计算异形滑块的后壁倾角、后壁面积、方量等特征数据，确定异形滑块的稳定系数。从计算结果可以看出，试验区大部分块体为危险岩体（K<1.05），只有1 个块体为稳定的异型滑移式块体（L4:K=1.42），解决了人工调查危岩精度较低的问题，与现场稳定性评价结果一致。
（3）对识别出的异形滑块进行特征计算，可快速初步确定危险岩体的危险等级。采用高差和体积相结合的方法对危险岩体进行评价。其中危害性最大的是一块大型危岩体（L2），高差为10.60 m，体积为2 398.60 m3，对危岩体的提取与准确的评价为边坡工程的数字化奠定了一定的基础。
猜你喜欢 后壁边界点高差 超高速撞击球形压力容器后壁损伤破坏试验研究振动与冲击(2021年23期)2021-12-20道路空间特征与测量距离相结合的LiDAR道路边界点提取算法测绘学报(2021年11期)2021-12-09高差影响下的城镇燃气管道水力计算简化公式煤气与热力(2021年10期)2021-12-02层次化点云边界快速精确提取方法研究激光技术(2021年5期)2021-08-17深刺翳风联合点刺咽后壁为主治疗脑卒中后吞咽障碍的临床观察中国民间疗法(2020年22期)2021-01-07框架结构梁板面钢筋叠合产生的高差问题探讨装饰装修天地(2020年11期)2020-07-04同时对向间接高差精密测量技术应用研究铁道勘察(2018年5期)2018-10-22地形高差较大的别墅区排水设计要点分析房地产导刊(2015年6期)2015-10-21开放式手术耳后多种自体中胚层组织填充乳突修复耳道后壁的临床效果观察疑难病杂志(2014年12期)2014-04-16一种去除挂网图像锯齿的方法及装置电脑与电信(2014年6期)2014-03-22