1. 引言
早期危岩体识别主要依赖人工现场作业获取结构面倾向、倾角、迹长等几何参数结合点荷载试验测定力学参数进行判定。该方法实时性强,但对人员经验依赖度高,存在信息不全、主观偏差等问题,且在高陡边坡等危险区域适用性差。
20世纪70年代,遥感技术引入工程领域,数字正射影像(Digital Orthophoto Map, DOM)、数字表面模型(Digital Surface Model, DSM)、数字高程模型(Digital Elevation Model, DEM)等遥感影像被用于室内解译危岩体。李娟[1]通过典型区遥感解译数据研究地震诱发地质灾害发育特征及分布规律;田尤等[2]基于DOM对甘肃天水麦积区幅内各类地质灾害进行目视解译,提出黄土区滑坡、崩塌、泥石流的影像解译特征;Derron等[3]利用DEM识别地貌特征与断裂集评估危岩体;Rossi等[4]结合DEM与岩性数据训练崩塌概率模型;Loye等[5]通过DEM坡角分布确定落石源区。然而,早期遥感影像分辨率不足[6],导致误判、漏判频发,且难以获取岩体结构参数。
进入21世纪后,数字摄影测量技术在我国实现了跨越式发展,数字摄影测量工作站在中国的摄影测量生产中获得了普遍的应用与推广[7]基于关联学科的进步,近景数字影像处理技术进入实用化阶段。汪磊[8]完成了一套实用的数字近景摄影测量系统设计并给出实验结果;刘子侠[9]结合控制架与VirtuoZo软件快速获取破碎岩体信息。韩东亮[10]运用数字近景摄影测技术解译大量岩体结构面几何信息,为岩体稳定性评价提供了大量可靠的基础数据。但是该方法受限于摄影距离与操作复杂度,应用范围有限。
随着计算机技术的发展,三维激光扫描作为获取空间数据的有效手段,能够快速地获取反映客观事物实时、动态变化、真实形态特性的信息[11]。Slob等[12] [13]首次用激光扫描分析岩体不连续性,可为工程提供不连续岩体内部结构特征;Feng等[14]使用全覆盖3D激光扫描技术对岩石表面进行原位测绘和记录;董秀军等[15]结合工程实例,阐述了应用三维激光扫描解决高陡边坡调查中,关于边坡快速编录和岩体结构面参数测量的原理与方法;刘昌军等[16]利用三维激光测量技术获取高陡边坡的高密度点云坐标数据,建立高边坡表面模型,提出危岩块体的识别方法。激光扫描精度高、速度快,但存在盲区与设备成本问题。
近年来无人机与三维建模技术突性发展[17],无人机凭借机动性强、成本低、精度高等优势,结合影像后处理技术,成为危岩体研究的主流工具。黄海宁等[18]以锦屏二级水电站出线场边坡落石灾害所在区域为例,将无人机摄影测量技术应用于高陡边坡危岩体调查中;马显东等[19]利用Fast-Marching算法对优势结构面进行三维重建,实现对岩体结构面几何信息和空间分布状态的特征提取;陈娜[20]对Cloud Compare的形状提取插件RSD (Ransac Shape Detection)做了二次开发,将之改编成可以适用于岩体结构面提取的插件RDD (Ransac Discontinuity Detection);Albarelli等[21]利用无人机采集的点云数据进行局部尺度的落石易发性评估。
传统方法在高位危岩识别中存在三方面不足:一是难以获取陡崖顶部等隐蔽区域的几何与裂隙细节;二是多源数据与算法协同不足,依赖人工调整导致冗余;三是传统点云算法抗噪性差,易致边界误判。本研究通过三维点云局部协方差与曲率分析增强地形细节捕捉,结合自适应聚类算法实现数据–算法深度耦合,并利用曲率差异区分危岩体、协方差滤波消除非结构性噪声,最终突破地形限制,达成高位危岩体的非接触式高精度自动化识别。
2. 研究区域概况
本研究区域位于中国广西壮族自治区桂林市临桂区会仙镇矮山村,地处扬子准地台与华南准地台构造单元交界区域。该区域历经多期次构造运动改造,形成以北东向和南北向为主导的复杂构造格局。气候特征表现为典型的中亚热带季风性湿润气候,年平均气温20.8℃,年降水量2215毫米,降雨呈现明显的季节集中性与暴雨高频特征。水系发育受珠江水系控制,主要包含桂柳运河、采空区积水及上笑河等水体单元。地貌类型以岩溶峰丛谷地为典型代表,广泛发育孤峰、残丘、溶蚀洼地等喀斯特地貌形态。斜坡坡度普遍处于45˚至80˚区间,采矿活动形成多级高陡人工边坡及临空结构面。岩体破碎特征显著,节理裂隙发育程度高,存在危岩体崩塌等潜在的地质灾害风险。地层构成以泥盆系碳酸盐岩为主,包含灰岩和白云岩两类典型岩性,岩石处于中风化状态,力学强度偏低,采矿工程扰动加剧了岩体结构劣化进程。水文地质条件具有显著复杂性,地下水系统包含裂隙溶洞水和覆盖型岩溶水两种主要类型,区域径流模数介于3~6 L/s·km2之间。地震活动性评估显示该区属于相对稳定区域,地震基本烈度评定为VI度。当前研究区内高陡岩质边坡及不稳定斜坡的潜在失稳问题,已对矿山安全生产及邻近居民区构成实质性威胁,亟待开展系统的稳定性评价与风险防控研究。研究区三维地图概览见图1。
Figure 1. Overview of the three-dimensional topography of the study area.
图1. 研究区三维地形概览
3. 技术方法
本研究针对复杂地形影像覆盖难题,构建了仿地倾斜摄影与多角度贴近摄影相融合的自适应测量方法。具体实施过程中,采用大疆Phantom 4 RTK无人机搭载高精度航测相机作为数据采集平台,通过设计分层飞行策略解决地形高差显著区域的影像覆盖问题。在贴近摄影阶段,运用坡度坡向分析方法与RANSAC平面拟合算法生成多角度航摄方案,实现目标区域的多视角数据获取。数据处理环节应用iTwin Capture Modeler Center软件完成三维重建,具体技术路线包括空三解算、密集点云生成及纹理映射等关键步骤。为提高模型精度,研究采用各向异性滤波方法与光束法平差技术进行联合优化。针对点云数据特性,实施统计滤波与半径滤波差异化处理方法以去除噪声干扰,通过八叉树采样技术优化数据密度分布。在危岩体识别方面,运用边缘跟踪技术提取裂隙特征,结合曲率阈值筛选方法定位目标点云区域,最终采用DBSCAN聚类算法实现危岩体的自动分类与识别。技术路线图见图2。
Figure 2. Technology roadmap
图2. 技术路线图
3.1. 点云数据采集
本研究采用仿地倾斜摄影测量技术获取地形基础数据,选用大疆精灵Phantom 4 RTK无人机(图3)作为数据采集设备。该矿区面积0.165平方公里,海拔区间为+99米至+304米,地形高差达205米。由于矿区面积较小且地形高差显著,采用分层飞行方案平衡影像分辨率与飞行安全。基于iflier软件的三维航线规划系统,设定120米、220米、285米三个飞行高度层(图4)。其中120米低空层重点覆盖谷底区域,获取厘米级高分辨率影像;220米中层兼顾覆盖效率与影像质量;285米高空层规避陡坡碰撞风险,其分辨率损失通过后续贴近摄影进行补偿。航向重叠度设为85%,旁向重叠度设为70%,该参数组合可确保陡坡区域影像特征点匹配成功率超过85%。各航线分辨率及拍摄效果见表1。
根据斜坡地形和结构面大致情况,规划多角度贴近航线,设计影像重叠度为航向85%,旁向70%,根据各斜坡作业面积和实际情况,设计合适的分辨率等相关参数。航线规划完成并导入无人机地面控制器后,无人机即可按照计划参数自动执行拍摄任务并采集影像。贴近摄影前后三维点云比对见图5。
Figure 3. DJI Phantom 4 RTK UAV
图3. 大疆精灵Phantom 4 RTK无人机
Table 1. Resolution and imaging performance of flight paths
表1. 各航线分辨率及拍摄效果
航线 |
分辨率(cm/像素) |
拍摄效果 |
120 m |
3.2 |
影像高分辨率 |
220 m |
6 |
兼具效率与覆盖 |
285 m |
7.8 |
规避陡坡碰撞风险,同时覆盖全地形 |
Figure 4. Flight path planning for UAV multi-layer oblique photography
图4. 无人机分层倾斜摄影航线规划图
(a)贴近摄影前 (b)贴近摄影后
Figure 5. Comparison of 3D point clouds before (a) and after (b) close-range photography (top view of Unstable Rock 1)
图5. 贴近摄影前(a)、后(b)三维点云比对(危岩1俯视)
3.2. 三维实景建模
本研究采用iTwin Capture Modeler Center (CC)软件完成边坡三维重建,生成包含三维点云数据、实景三维模型、数字正射影像(DOM)及数字表面模型(DSM)的四维地理信息成果。具体技术流程包含六个核心环节:首先执行数据导入操作,将无人机采集的影像数据载入系统;其次开展自动空三解算,通过特征点匹配算法计算相机内外方位参数并生成稀疏点云;第三阶段实施密集匹配运算,基于空三解算成果构建高密度三维点云;第四环节运用网格生成算法构建三维模型基础框架;第五步骤完成纹理映射处理,实现影像信息与三维模型的几何配准;最终进行模型优化操作,包括噪声滤除与几何空洞填补,并输出OBJ、FBX等标准格式模型文件。研究区及临近区域点云实景模型见图6。
(a)模型整体效果
(b)模型局部效果 (c)模型细节效果
Figure6. (a) Real-scene point cloud model of the study area and adjacent regions; (b) Partial view of the point cloud model; (c) Detailed view of the point cloud model
图6. (a)研究区及临近区域点云实景模型;(b)点云模型局部;(c)点云模型细节
针对地面高程数据采集的特殊性,研究区域存在植被覆盖导致的高程测量误差问题。常规处理方法需采用插值算法消除地表覆盖物的干扰,但本研究的石灰岩矿山危岩体分布于裸露开采面,不存在植被遮挡等影响因素。经实地勘验确认,研究区域数字表面模型可直接替代数字高程模型使用,该数据应用方案有效规避了插值运算引入的不确定性误差,保证了高程数据的原始精度。据此建立的研究数据集已通过精度验证,满足后续分析的可靠性要求。
3.3. 点云数据预处理
三维点云数据包含采集点的三维空间坐标与对应的RGB色彩信息。在实际数据采集过程中,由于边坡表面植被覆盖、灰尘堆积及仪器固有误差等客观因素,点云数据中普遍存在噪声点与局部数据空洞现象。这些异常数据将直接影响后续点云法向量计算的可靠性、结构面识别的准确性以及几何特征提取的精度,因此必须对原始数据进行去噪与空洞修复的预处理操作。
针对边坡区域点云数据中广泛分布的破碎小面积结构面,需在预处理阶段进行选择性筛选与剔除处理。通过开源软件Cloud Compare集成的统计滤波与半径滤波组合算法,能够有效区分并移除离散噪声点:统计滤波基于邻域距离分布特性,设定邻域点数k为50与标准差倍数σ为1.8,剔除偏离均值1.8倍标准差外的离散点;半径滤波进一步清理孤立噪声,设置搜索半径r为0.15 m与最小邻域点数n为4,确保保留连续曲面特征。对于因传感器视角限制形成的局部数据缺失区域,采用泊松曲面重建技术或径向基函数(RBF)插值方法进行修复。预处理后的点云数据量显著减少,不仅提升岩体结构面识别效率,同时通过消除噪声干扰有效提高特征识别精度,为边坡稳定性分析建立可靠的数据基础。
点云密度参数直接决定三维空间数据的解析能力,其数值大小与目标地物特征的表征精度呈现正相关关系。较高密度的点云数据能够精确还原微观地形结构与复杂几何形态,但会导致数据量呈指数级增长,显著增加计算资源消耗。研究过程中需通过密度参数的动态优化,在保证三维建模精度的基本要求下,实现计算效率与硬件资源消耗的平衡控制。
点云重采样技术包含上采样、下采样和均匀采样三种处理方式。针对高密度边坡点云数据集,采用下采样方法从原始点云中有选择性地提取部分数据点,在保留主要几何特征的前提下降低整体数据密度。该过程不仅能实现点云数据的稀疏化处理,还能同步过滤具有较大高程起伏的异常离散点。Cloud Compare软件提供三种下采样模式:随机采样模式(Random)、固定空间间隔采样模式(Space)以及基于空间划分的八叉树采样模式(Octree) (表2)。
Table 2. Comparison of three downsampling methods in CloudCompare
表2. CloudCompare三种下采样方式对比
维度 |
Random |
Space |
Octree |
计算速度 |
极快(O (n)) |
中等(O (n log n)) |
较慢(依赖树构建) |
细节保留 |
差(全局随机丢失) |
中等(局部特征保持) |
优(多尺度自适应) |
密度控制 |
仅总量控制 |
严格空间约束 |
动态分级控制 |
内存消耗 |
低(无需结构存储) |
中(需体素网格) |
高(树节点存储) |
典型应用 |
数据预览、降维 |
均匀化处理、数据融合 |
复杂结构处理、多尺度分析 |
本研究使用的初始点云数据由iTwin Capture Modeler Center软件生成,总体数据规模达到6.1 GB,包含1.8亿个三维数据点。鉴于研究涉及复杂结构面处理且需保持高细节特征,选择八叉树采样模式进行降采样处理。
Figure 7. Octree (subdivision level: 10)
图7. 八叉树(细分等级:10)
Figure 8. Point cloud data after subsampling
图8. 二次采样后点云数据
通过设定八叉树细分等级为十级,最终获得包含135.9万个数据点的优化数据集(图7、图8)。选取研究区内结构面密集交错且具有显著几何突变的区域(图9)作为重点分析对象,该区域被判定为岩体失稳高风险区,其精细化处理的点云数据为危岩体特征参数提取提供高精度基础。
Figure 9. Unstable rock identification zone
图9. 危岩识别区
3.4. 危岩体点云跟踪识别
本研究通过分析危岩体表面点云的曲率特征实现边缘裂隙识别。具体而言,当表面形态呈现不连续特征时,邻近点法线向量将呈现显著分散状态。基于此现象,采用曲率值(Surface Variation, SV)作为判别指标,通过识别高曲率区域可有效定位危岩边缘裂隙的空间分布。
曲率值的大小直接反映采样点的空间偏离程度:当SV值增大时,表明采样点偏离切线平面的距离增加,对应局部曲面弯曲度显著提升;反之,低SV值则意味着采样点趋近于切线平面,反映局部曲面平滑特征。据此建立空间判别准则,利用Cloud Compare软件的按值筛选功能,保留曲率值高于阈值的点云数据,这些点云对应岩体结构的脊谷特征点(图10~15)。
Figure 10. Fracture tracking results in unstable rock identification zone
图10. 危岩识别区岩体裂隙跟踪结果
Figure 11. Calculated SV values in unstable rock identification zone
图11. 危岩识别区SV值计算结果
Figure 12. Real-scene point cloud of Unstable Rock 1
图12. 危岩1实景点云
Figure 13. Edge tracking of Unstable Rock 1
图13. 危岩1边缘跟踪
Figure 14. Real-scene point cloud of Unstable Rock 2
图14. 危岩2实景点云
Figure 15. Edge tracking of Unstable Rock 2
图15. 危岩2边缘跟踪
在完成脊谷点云提取后,需验证其空间闭合形态。通过形态学分析判断脊谷点云是否形成闭合结构,若满足闭合条件,则将该部分点云从原始数据中剔除。剩余点云数据包含两类主要成分:由闭合脊谷点围限形成的危岩块体点云,以及构成边坡基岩的背景点云。
针对分割后的三维点云数据集,采用基于密度聚类的DBSCAN算法进行分类处理。Martin Ester提出的DBSCAN算法通过密度可达性分析实现空间聚类,其核心流程包含三个阶段:首先设定邻域半径eps和最小点数阈值minPoints,从任意数据点出发计算其ε邻域内的点数密度;其次根据密度判据划分聚类类别,若邻域内点数达到阈值则建立新聚类,否则标记为噪声点;最后迭代处理全部数据点,直至所有点云完成分类标记。
在DBSCAN算法参数优化与实验验证过程中,通过融合局部密度特征分析和全局几何约束推导,结合系统性实验框架实现三维点云分类。首先基于点云局部密度特性,计算无显著几何特征区域的平均邻域距离,获得基准值0.05以反映基础分布规律;其次结合曲率空间分布特征,分析裂缝等显著结构边界的最大延伸范围,推导出0.31的全局分割阈值。为平衡局部密度敏感性与全局结构完整性,将邻域半径确定为中间值0.25,该参数既可规避过小半径导致的噪声碎片化问题,例如eps = 0.1时背景噪声被误识别为独立聚类,又能抑制过大半径引发的分类混淆现象,例如eps = 0.3时不同曲面点云异常合并,同时将最小邻域点数固定为4以符合三维空间邻域连通性的最小拓扑需求。实验验证采用三阶段递进流程:基于0.05至0.31的参数探索范围,通过步长0.05的迭代测试发现,当eps = 0.25时目标点云形成完整闭合结构且与背景有效分离,而较小参数导致噪声碎片化,较大参数引发空间重叠错误;最终通过点云边缘跟踪技术验证聚类边界清晰度,结合闭合结构删除后的数量统计确认分离效果,跨数据集测试表明该方法在保持曲面连续性的同时实现了目标结构与背景噪声的精准分割(图16)。
Figure 16. Point cloud data after edge segmentation
图16. 分割边缘点后点云数据
在分类结果中,由曲率值分割的边缘围限的点云簇即为目标危岩体。基于Python的DBSCAN聚类库编写“点云聚类与分割”代码,将这些危岩体和非危岩体分为不同的类别,并且给每一个类别都附上显眼的颜色(图17~19)。
Figure 17. 3D visualization of DBSCAN clustering results (Python environment)
图17. DBSCAN聚类结果三维可视化(Python环境)
Figure 18. 3D point cloud classification via density-based clustering
图18. 基于密度聚类的三维点云分类结果
(a)危岩1实景点云 (b)危岩1几何点云
(c)危岩2实景点云 (d)危岩2几何点云
Figure 19. Comparison between classified point clouds and real-scene point clouds: (a) (b) Unstable Rock 1; (c) (d) Unstable Rock 2
图19. 分类后的点云与实景点云比对:(a) (b)危岩1、(c) (d)危岩2
经实测验证,危岩1体积为389.03 m3,危岩2体积为293.23 m3,其顶端高差分别达到98 m和119 m。对照《危岩防治工程技术规范》(DB45/T 1696-2018)标准,判定为高位小C型危岩。对比分析实景点云与分类结果(图19),验证了该算法在复杂岩体结构中的适用性,分类精度满足工程勘察需求。
3.5. 危岩块体失稳运动特征分析
基于识别出的危岩体(危岩1体积389.03 m3,危岩2体积293.23 m3),结合Rocscience Dips软件分析其失稳破坏模式及概率。结果表明:
平面滑动模式在极射赤平投影当中受边坡面投影、摩擦圆、潜在滑动面包络以及边界限制线所控制,破坏区域为潜在滑动面包络、摩擦圆与边界限制线所围限的区域。危岩1和危岩2所在边坡的坡向与坡角分别为168˚∠73˚和235˚∠65˚;结合附近工程所得岩体参数,内摩擦角取25˚;平面滑动主要发生在结构面倾向和边坡面倾向成一定角度的范围内,通常认为此角度为20˚~30˚,边界限制线设定为25˚,平面滑动的概率计算结果如图20、图21。
Figure 20. Planar sliding probability of Unstable Rock 1
图20. 危岩1平面滑动概率
Figure 21. Planar sliding probability of Unstable Rock 2
图21. 危岩2平面滑动概率
楔形体滑动破坏模式主要作用区域在于摩擦圆内部与边坡面外侧形成的新月形区域,次要作用区域由摩擦圆、边坡面圆和辅助平面圆之间的区域组成。楔形体滑动的概率计算结果如图22、图23。
Figure 22. Wedge sliding probability of Unstable Rock 1
图22. 危岩1楔形体滑动概率
Figure 23. Wedge sliding probability of Unstable Rock 2
图23. 危岩2楔形体滑动概率
倾倒破坏区是由边界限制线、滑动限制面以及投影网界限所限定的范围,倾倒破坏发生概率如下图24、图25。
Figure 24. Toppling failure probability of Unstable Rock 1
图24. 危岩1倾倒破坏概率
Figure 25. Toppling failure probability of Unstable Rock 2
图25. 危岩2倾倒破坏概率
综合上述三种破坏模式,得到的结果如下表3。
Table 3. Probability of failure modes for unstable rock masses
表3. 危岩体破坏模式发生概率
破坏几率/(%) |
平面滑动 |
楔形体滑动 |
倾倒破坏 |
危岩1 |
26.92 |
41.86 |
0 |
危岩2 |
5.77 |
50.96 |
0 |
危岩1平面滑动的发生概率为26.92%,若发生滑动,沿结构面J3滑动的概率为91%,危岩2发生平面滑动的概率较低;而对于楔形体滑动,两者的概率都比较高,危岩1更倾向于沿J2与J3结构面组合方向发生滑动,危岩2则更倾向于沿J1与J2结构面组合方向发生滑动,均为后缘有陡倾裂隙的滑动破坏模式;对于倾倒破坏,由于未满足倾倒破坏的条件,两者的概率均为0。最不利结构面见图26~27。
Figure 26. Most unfavorable structural plane of Unstable Rock 1
图26. 危岩1最不利结构面
Figure 27. Most unfavorable structural plane of Unstable Rock 2
图27. 危岩2最不利结构面
3.6. 适用性讨论
本文所提出的危岩识别方法有效地识别出了该高陡边坡的裸露孤立危岩体,但在本方法的基础上,仍有许多值得深入研究的问题。
1) 本文基于曲率与聚类算法的识别方法在实际复杂地形中易受噪声干扰,可能导致结果偏差。后续研究可结合多源遥感数据(如纹理、光谱等),开发自适应的邻域半径与最小点数动态优化算法,以提升算法的环境适应性与识别精度。
2) 本文结构面提取算法虽能处理局部点云数据,但在大规模场景下计算效率与内存消耗仍存在瓶颈。未来可探索分布式计算框架与并行优化技术,降低算法资源占用,同时提升复杂点云数据的工程化处理能力。
4. 结果与讨论
本研究以桂林市临桂区会仙镇大神山石灰岩矿山边坡为研究对象,通过无人机倾斜摄影与多角度贴近摄影技术获取高分辨率影像数据,结合三维重建与点云预处理方法,提出了基于DBSCAN点云聚类算法的危岩体识别技术体系。研究结果表明:
1) 采用分层飞行策略与自适应航摄方案,有效克服了高陡边坡地形复杂、高差显著导致的影像覆盖难题,实现了研究区高精度三维点云数据的采集。
2) 通过点云预处理(去噪、重采样)与曲率分析,结合DBSCAN密度聚类算法,成功将危岩体点云与基岩点云分离,识别出危岩1 (体积389.03 m3)与危岩2 (体积293.23 m3),其分类精度满足工程调查需求,验证了非接触式自动化识别方法的可行性。
3) 测量出危岩1体积为389.03m3,危岩2体积为293.23m3,危岩1顶端距离坡脚高差为98 m,危岩2顶端距离坡脚高差为119 m,均评级为高位小C型危岩。对危岩体的准确评价,能够为后续危岩治理提供基础。
NOTES
*通讯作者。