岩质边坡结构面交互式识别技术与系统研发

2024-03-02代永新卢敬标李子健张思远杨月明

金属矿山 2024年1期

代永新徐全卢敬标李子健张思远杨月明

(1.中钢集团马鞍山矿山研究总院股份有限公司,安徽马鞍山 243000;2.山东大学土建与水利学院,山东济南 250061;3.非煤露天矿山灾害防控国家矿山安全监察局重点实验室,安徽马鞍山 243000;4.金属矿山安全与健康国家重点实验室,安徽马鞍山 243000)

矿山岩土工程安全是矿产资源安全高效开采的基础保障。随着科技的迅猛发展,矿产资源开采已向更高、更深迈进,产生了大型露天矿山边坡,同时也产生了一系列复杂的岩土工程问题。其中,大型露天矿山边坡以岩质边坡为主,结构面的空间展布和组成特征是决定岩质边坡稳定性的重要因素。因此,准确获取边坡岩体的结构面对于边坡的稳定性分析至关重要[1-5]。

目前,结构面的获取方式大致有接触式测量和非接触式测量两种。接触式测量主要是在现场设置测窗或测线,使用罗盘和皮尺对结构面进行人工测量,该方法受限于地形条件、工作效率及工作人员作业水平,难以实现大规模、高效准确的测量[6],并且由于边坡高陡特征以及工程现场环境的复杂性往往存在较大的安全隐患。非接触式测量方法主要包括:① 钻探技术,通过钻孔方式获取岩芯确定产状或利用光学设备探入钻孔内获取岩体结构信息[7-11];② 三维激光扫描法,利用三维激光扫描仪获取岩体结构面高精度点云数据,再根据点云的几何特征识别岩体结构面[11-14];③ 数字摄影测量法,基于数字图像与摄影测量的基本原理,应用计算机三维成像技术、影像匹配、图像插值多学科理论与方法,构建拍摄对象的三维实景模型,并根据实景模型获取岩体结构面信息[15-19]。

数字摄影测量法是一种非常优秀的岩体结构面非接触式测量法,该方法的迅速发展主要得益于近年来无人机技术的快速发展和广泛使用。业内学者借助该技术主要从两个方面进行结构面识别:① 基于该技术获取岩质边坡点云或网格的法向量特征识别岩体结构面[20-23];② 生成三维实景模型,通过交互方式获取岩体结构面,国外在该方面已有相关成熟的产品,如奥地利3GSM、澳大利亚SIROVISION。目前,国内主要研究基于点云的岩体结构面获取方法,但该方法存在显著不足[20],无法有效满足工程应用需求。具体来说,无法有效识别出以迹线方式出露的结构面,易将碎石覆盖形成的平面或人工开挖形成的断口误判为结构面,对于小尺度的结构面易被忽略。对于交互识别系统,业内学者主要使用奥地利3GSM、澳大利亚SIROVISION 成熟产品或开源工具性软件进行相关研究。如毛权生等[24]使用3GSM 公司的ShapeMetriX3D 数字摄影测量系统对南泥湖钼矿采场高陡边坡进行了结构面采集和分析;姚富潭等[25]基于倾斜摄影技术获取危岩体的三维实景模型、三维空间点云等数据成果,并通过DSE 程序法[26]实现了岩体结构面产状信息的半自动提取。

基于三维实景模型的交互方式识别岩体结构面,其清晰度和精度均优于基于点云或网格的识别方式。因此,本研究以图形库技术为基础,编写能够满足三角网格、贴图信息、点线面显示等元素渲染的函数集,将倾斜摄影获取的岩质边坡三维实景模型进行可视化分析,结合Arcball 交互技术并基于实景模型以点、线、面的方式对结构面进行交互识别。同时,将结构面优势组划分、间距统计、迹长统计、边坡破坏模式分析等功能进行组合,开发了边坡岩体结构面交互识别系统。该系统能够对边坡岩体结构面进行高效准确的识别,有助于提高边坡稳定性分析的可靠性、及时性和准确性。

1 方法与原理

1.1 边坡高精度三维实景模型构建

在岩体边坡结构面识别过程中,建立高精度的三维实景模型是关键环节,实景模型包含了空间坐标和纹理等多种信息,为进一步进行岩体结构面交互识别与分析奠定了坚实基础。可以采用多种方法构建高质量的三维实景模型,如摄影测量与立体视觉、激光扫描、立体摄影与重建、全景相机以及深度学习技术等。具体的方法选择取决于场景复杂性、经费预算、精度要求和可用技术设备。典型的高精度边坡三维实景建模流程如图1 所示,包括信息采集、三维重构、三维实景模型生成与输出3 个典型环节。

图1 边坡高精度三维实景建模流程Fig.1 High precision 3D real scene modeling process of slope

图1 中,“三维重构”是核心流程,主要采用SFM算法进行图像三维重建和测量[27-32]。 SFM 算法通过分析多幅图像中的视觉特征和相机运动,可以估计出场景中的三维点云和相机姿态,从而实现对物体或地面的三维建模。该算法具体实现流程包括特征提取与匹配、对齐与稀疏重构、优化与稠密重建、网格生成、纹理映射与模型生成。然而,SFM 算法也存在一些不足:① 在获取精度较高且范围较广的模型时,三维重建流程非常耗时,且对计算机的硬件条件有一定要求;② 模型坐标以初次拍摄照片时相机的局部坐标系为基准,导致模型在空间方位上与实际测量中的东北天(ENU)坐标系不符合;③ 模型尺度与真实尺度存在缩放关系,无法直接量取真实距离。

因此,在“信息采集”流程中,除了需要采集图像信息以外,还需要获取方位、尺度信息来解决SFM 法存在的尺度与方位不确定性问题。在实际应用中,通常针对不同尺度的边坡场景,选择不同的手段构建高精度三维实景模型。对于小尺度边坡场景,可直接使用相机或手机拍摄多角度图像(图1),并结合位姿校正设备获取建模信息,如文献[28]使用的“L”形标尺或3GSM 中使用的圆盘。对于大尺度边坡场景,可借助无人机进行倾斜摄影获取图像,并通过空三解算构建大场景边坡三维实景模型,其中无人机自带的RTK 或人工布设的像控点可用于位姿控制。

1.2 岩体结构面识别与解译

结构面在岩体露头中主要以迹线、面的方式呈现,如图2 所示。基于高精度三维实景模型的结构面识别思路,主要是通过观察高精度岩质边坡三维实景模型,并选取结构面上不共线的3 个或多个控制点(特征点),再对结构面的产状、间距、迹长等信息进行解译。

图2 岩体中结构面出露方式Fig.2 Appearance of structural plane in rock mass

1.2.1 结构面产状求解

岩体结构面产状一般用倾向α(0°～360°)和倾角β(0°～90°)表示,岩体结构面倾向和倾角的空间几何关系如图3 所示。图3 所示三维空间坐标系中,X轴正方向对应地理正东方向,Y轴正方向对应地理正北方向,Z轴正方向为垂直向上方向。高精度岩质边坡三维实景模型中的坐标系与图3 坐标系相同,可在三维实景模型中根据所选结构面的控制点,采用最小二乘法拟合最优平面,求得最优平面的法向量N(a,b,c) ,并据此求解结构面产状(α,β) 。具体过程如下:

图3 结构面法向量与产状的空间关系Fig.3 Spatial relationship between normal vector and occurrence of structural plane

设在岩质边坡实景模型中选取了n个特征点,记为Pi(Xi,Yi,Zi),其中i∈0 ～n,n个特征点的最优拟合平面F的表达式为

式中,A、B、C、D为平面方程的待定系数。

记第i个特征点Pi(Xi,Yi,Zi)到平面F的距离为di,即:

式(3)的矩阵形式为

由式(4)可求得平面F的待定系数A、B、C、D。平面F的单位法向量N(a,b,c) 可通过下式求解:

由于在三维空间中,平面F法向量N的方向是可以指向上半空间或下半空间,为满足岩体结构面产状表达的需要,当法向量N指向下半空间(c＜0)时,则将法向量N取反,即令a、b、c分别取负。

结构面法向量与产状的空间关系如图3 所示。根据图3 所示空间关系,结构面产状(α,β) 可通过下式进行求解:

1.2.2 结构面迹长求解

根据选取的n个特征点集Pi(Xi,Yi,Zi),除了可以计算出结构面的产状外,还可以根据最小包围圆算法[29]计算出结构面的迹长,如图4 所示。具体步骤为:① 将n个特征点集Pi(Xi,Yi,Zi)映射到平面F,得到点集P′i;② 在点集P′i中任意取3 个点P1、P2、P3;③ 以选取的3 个点构造最小包围圆O;④ 在点集P′i中查询距离O的圆心最远的点pv;⑤ 判断pv是否在圆O内,若在,则终止流程,最小包围圆O的直径DO即为迹长;否则,继续执行步骤⑥;⑥ 在{P1、P2、P3、Pv}中选取3 个点,构造包含这4 个点的最小包围圆O′,转到步骤③,其中选取的3 个点尽可能是边界上的点。

图4 结构面迹长求解示意Fig.4 Schematic of structural plane trace length solution

1.2.3 结构面间距求解

结构面间距表示相邻且同组结构面间的法向距离,如图5 所示。在模型中测量时,可根据测线法进行求解;还可以将每组结构面映射到最优产状对应的平面,再求解相邻平面的法向距离D。

图5 结构面间距求解示意Fig.5 Schematic of solving structural plane spacing

本研究岩质边坡结构面交互式识别系统将采用后者方法进行结构面间距求解,过程如下:

(1)选取结构面1,设结构面1 所属优势组产状为(α,β) ,对应的平面法向量为N(a,b,c) ,则结构面1 的平面表达式为A·X+B·Y+C·Z+D1=0。

(2)寻找与结构面1 同组且相邻的结构面2, 则结构面2 的平面表达式为A·X+B·Y+C·Z+D2=0。

(3)求解结构面1 和2 的间距D,公式为

式中,D1和D2分别为结构面1 和结构2 平面表达式中的系数。

1.3 岩质边坡运动学分析

岩质边坡的稳定性往往由影响岩体的不连续面控制,在这种条件下,结构面产状与边坡方向在空间上的几何关系决定了边坡的稳定性和潜在的破坏模式。岩质边坡运动学分析是指在不考虑岩体中所存在的地应力情况下,仅对岩石块体可能产生的几何运动进行评估。岩质边坡运动学分析具有简单和快速的特点,无需进行复杂建模且计算耗时较少,通常用来评估平面破坏、楔形破坏、弯曲破坏、倾倒破坏4 种破坏模式[30]。

系统设计中,主要根据裸露在岩质边坡表面所统计的结构面以及边坡方向进行运动学分析,具体思路如下:

(1)平面破坏。不连续结构面倾向αd与岩质边坡倾向αs平行或接近平行(通常±20°以内),不连续结构面倾角βd小于边坡倾角βs,不连续结构面倾角βd大于结构面摩擦角φ。

(2)楔形体破坏。 2 个不连续结构面相交,且相交线倾向αi与边坡方向相同,交线的倾角βi大于结构面的摩擦角φ,交线的倾角βi小于边坡的倾角βs。

(3)弯曲破坏。不连续结构面反倾向边坡,不连续面倾角小于(90°-斜坡倾角+摩擦角φ)。

(4)倾倒破坏。 2 个不连续结构面相交,且相交线反倾向边坡,相交线倾角小于(90°-边坡倾角βs)。

2 岩质边坡结构面交互式识别系统设计与开发

2.1 系统设计

本研究系统设计融合了机器视觉、计算机图形学、统计学等相关领域理论,从底层技术开发,编写能够满足三角网格、贴图信息、点线面显示等元素渲染的函数集,开发三维实景模型的可视化模块,通过输入边坡高精度三维实景模型,并结合计算机视觉、数学几何以及交互式可视化技术,实现边坡岩体结构面拾取、聚类、赤平投影分析以及可视化等功能,系统构成如图6 所示。该系统设计了多个模块和功能丰富的工具箱,提供了多样化的数据处理与应用,包括结构面识别、优势组求解、间距统计、距离测量等功能。整个系统采用C#语言开发,具有直观、友好的界面和高度交互性。

图6 边坡岩体结构面识别系统构成Fig.6 Structure plane identification system of slope rock mass

该系统为用户提供了一个直观且功能丰富的工作平台,主界面如图7 所示。主要界面包括可视化窗口、菜单栏、工具栏、边坡对象管理窗口、边坡属性窗口以及结构面对象管理窗口等核心组件。

图7 结构面识别系统界面Fig.7 Structural plane identification system interface

2.2 功能模块

系统主要功能模块包括输入模块、结构面识别模块、识别结果解译模块和工具箱模块。

(1)输入模块。输入模块(图8(a))的主要输入信息包括边坡高精度三维实景模型和位姿文件。其中,边坡高精度三维实景模型基于前文1.1 节中三维重构算法生成,采用标准的OBJ 文件格式保存模型信息,包括模型本体(. OBJ)、材质文件(. mtl)、贴图文件(.jpg,.png)3 种文件。位姿文件包含边坡高精度三维实景模型的方位与尺度信息。

图8 系统典型功能窗口Fig.8 Typical function windows of the system

(2)结构面识别模块。结构面识别模块是结构面识别系统的核心模块,主要根据前文1.2 节中的相关原理采用点、线、面方式对结构面进行交互式识别。

(3)识别结果解译模块。识别结果解译模块是对结构面结果信息进行计算解译的主要模块,包括结构面优势组划分(图8(b))、间距统计(图8(c))、迹长统计、迹线素描图生成、破坏模式分析(图8(d))。

(4)工具箱模块。工具箱模块主要包含系统比较常用的一些应用性功能,如距离测量、面积测量、产状测量、已识别结构面选择、已识别结构面查找、结构面信息输出等。

3 系统应用

3.1 工程背景

某露天矿山位于河北省滦州市城区约3 km,地理坐标为东经118°45′40″,北纬39°38′20″～39°39′42″,北距京山铁路滦州市车站8 km,西距迁(迁安)曹(曹妃甸)铁路菱角山站4 km,具体位置如图9 所示。矿区属暖温带半湿润大陆性季风气候,四季分明,多年平均气温10.5 ℃,极端最低气温-23.1 ℃,极端最高气温39.9 ℃,多年平均降雨量为668 mm。该矿地表长2 920 m,宽1 500 m,采场露天底部标高为-562 m。

图9 露天矿边坡(研究区)地理位置Fig.9 Geographical location of the open-pit slope(study area)

近年来,随着采矿生产不断推进,该矿露天采坑深度逐渐增加,导致边坡高度逐步上升。目前,采矿区域北侧边帮已经达到-157 m 以上的高度,已接近边界位置。与此同时,在采矿区域北侧拟建立皮带运输系统,其中帮高-127 m 平台上,计划建设皮带运输转运站。然而,经过对现有资料的分析以及多次实地踏勘,发现该区域边坡倾向约204°,整体边坡角约45°,台阶边坡角约60°,部分岩体受到不同程度的风化剥蚀,并在开挖和爆破震动作用下,存在局部不稳定的因素(图10),易发生单台阶滑塌、楔形破坏等灾害。为此,本研究通过无人机对该类潜在风险区域开展倾斜摄影测量,以建立高精度的三维实景模型。同时,采用所开发的岩质边坡结构面交互式识别系统,对该区域的岩体结构面进行准确识别,并进一步进行边坡潜在破坏模式分析,以便更有效地预测可能出现的问题。

图10 隐患区域分布Fig.10 Distribution of the hidden danger areas

3.2 高精度三维实景模型获取

本研究采用无人机开展倾斜摄影进行研究域三维实景建模,技术流程如图11 所示。主要步骤包括航线规划、照片采集、SIFT 特征检测与匹配、对齐与稀疏重建、优化与稠密重建、场景重建。

图11 研究区三维实景建模技术流程Fig.11 Technical flow of 3D real scene modeling in the study area

本研究摄影测量采用大疆M3E 型号无人机,并搭载RTK 模块。所构建的研究区三维实景模型如图12 所示。倾斜摄影地面影像分辨率优于0.05 m,空三加密基本定向点平面位置残差不大于0.25 m,高程残差不大于0.10 m。

图12 研究区三维实景模型Fig.12 3D real scene model of the study area

3.3 岩体结构面识别与应用

将前文3.2 节所获取的边坡高精度三维实景模型导入系统,采用点、线、面交互方式对边坡中的结构面进行识别,识别结果如图13 所示,共识别出938 条结构面。

图13 结构面识别结果Fig.13 Identification resluts of structural plane

结构面识别结束后,对已识别的结构面进行优势组划分,识别结构面赤平投影及优势组划分结果如图14 所示。根据结构面分布情况,共划分2 组优势结构面:第1 组J1,图中蓝色点集,倾向204°,倾角50°;第2 组J2,图中绿色点集,倾向320°,倾角110°。

图14 优势组划分结果Fig.14 Division results of dominant group

同时,采用前文1.2 节方法对结构面进行解译,结果见表1。由表1 可知:第1 组J1,共识别了8 2 8条,占总识别数量的88.3%,平均迹长Lm=6.0 m,平均间距Dm=0.09 m;第2 组J2,共识别了110 条,占总识别数量的11.7%,平均迹长Lm=7.05 m,平均间距Dm=3.92 m。

表1 结构面识别统计结果Table 1 Statistical results of structural plane identification

根据现场踏勘可知,研究区无层状发育的结构面,因此仅采用前文1.3 节中的运动学分析方法对该区域进行平面和楔形体2 种破坏模式分析,结果如图15 所示。由图15(a)可知:优势组J1的表征点在风险区域内,其倾向与边坡倾向基本一致,倾角小于台阶边坡倾角、大于结构面摩擦角,台阶边坡存在沿J1顺层滑动的可能;优势组J2的表征点在风险区域外,其倾向与边坡倾向大角度相交,倾角大于台阶边坡倾角,基本不会发生沿J2顺层滑动的可能。分析图15(b)可知:J1与J2投影大圆的交点位于破坏区域内,交线倾角小于台阶边坡倾角,且大于结构面摩擦角,台阶边坡存在沿交线方向产生楔形滑动的可能。