基于无人机摄影测量的露天矿边坡非连续变形分析

李金威

(黄河交通学院继续教育学院,河南 焦作 454950)

露天矿岩质边坡是由众多节理面、裂隙面、断层面等地质构造切割岩体所形成的复杂结构体系,与土质边坡构造具有较大差异,其稳定性不仅受制于岩石本身的强度,而且在很大程度上依赖于岩体结构面的力学性质,故岩质边坡结构面信息获取的准确性在其稳定性分析中尤为重要[1-2]。

传统结构面信息测量方法一般通过罗盘等工具直接对边坡进行接触测量,该方法工作效率较低、工作量偏大,具有一定的危险性,且对于高陡边坡,测量人员难以进行接触实测[3-4]。工程领域中三维激光和近景摄影测量技术逐渐被广泛应用,LI 等[5]利用激光扫描技术完成边坡三维建模,并根据获取的结构面信息进行了深度学习分组。王述红等[6]将近景摄影测量技术应用于岩体结构面信息采集,并据此建立岩体三维模型,实现了岩体的稳定性分析与关键块体治理。

非连续性是岩体变形过程固有的属性。DDA 方法[7]具有严密的理论及高效的计算方法,在模拟大变形、大位移等不连续变形方面具有独特优势。DDA方法在本质上与有限元方法(FEM)类似,它使用有限元类型的网格,但不同的是其所有块体都是真正的相互独立块,并由预先设置的不连续性构造限制其运动。DDA 块体可以是任意凸形或凹形,当块体间发生接触时,利用摩尔-库伦准则对接触界面进行处理,并且通过加载函数及步长设定平衡方程并求解。在FEM 方法中,未知数是所有节点的自由度之和,在DDA 方法中,未知数是所有块体的自由度之和,因此从理论上讲,DDA 方法是FEM 方法的推广[8]。

基于DDA 独具的优势,近年来许多研究人员对其初始代码进行了一些修改,以提高其性能并解决DDA 中的一些问题。AMIR 等[9]将DDA 中的现有块体作为可变形盘引入系统中进行求解,提高了求解效率。JIANG 等[10]提出了一种有效方法来解决由于大尺度转动而导致的体积增加误差问题,取得了良好的效果。MORGAN 等[11]提出了一种利用非连续变形分析进行水力压裂的地质力学模型,与Griffith 分析模型能够很好地匹配。FU 等[12]利用非连续变形分析方法开发了一种用于地下洞穴挖掘的系统。WU等[13]研究了大角度转动条件下的弹性扰动问题。YU 等[14]使用高阶位移函数描述块体的运动,提高了算法精度。ZHENG 等[15]引入了离散体结构和组件的双重形式解决了由于开—闭迭代和应用硬弹簧而引起的问题。FAN 等[16]通过虚功原理提高了变形计算的准确性,降低了大角度转动的误差。在DDA 能量耗散研究中,HATZOR 等[17-18]、TSESARSKY 等[19]深入探讨了DDA 动力系数与其运算结果的关系,证明了在DDA 中添加适当阻尼会使得其模拟与实际更加吻合。KOO 等[21]通过在DDA 计算步骤中添加校正函数来减少由于大块体旋转引起的误差,并通过增加质量比例阻尼来考虑能量损失,以满足大块体平移要求。姜清辉等[21]在模拟准静态过程时,运用变分机理来考虑阻尼力,以此表征其能量耗散形式。刘永茜等[22]采用自动调整时间步长的方法,在不同速度运动阶段选用不同的阻尼运算规则,取得了较好的运算结果。

露天矿边坡稳定性分析的关键在于求得安全系数及其滑落面,而这一求解过程的关键又在于岩体结构面信息的精准获取;再者岩质边坡破坏是一个动态过程,其中的能量耗散有必要充分考虑。以往研究鲜有针对结构面信息精准获取并将其用于非连续变形分析计算中,多数默认其精度满足计算要求,但在实际情况下结构面形态各异,须重点关注信息采集过程及其精度,并考虑其动态破坏过程的能量耗散,才可保证后续非连续变形分析的准确性。为了对上述问题提供解决方案,本研究借助无人机非接触摄影测量优势,提出多层次全方位无人机摄影测量方法获取边坡结构面信息,保证其结构面信息的采集精度,同时施加适当的黏性阻尼,确保非连续变形分析过程的能量耗散,并将改进后的方法应用到某露天矿岩质边坡的非连续变形分析中。

1 结构面信息采集

1.1 多层次全方位摄影测量

以某露天矿岩质边坡为研究对象,该边坡基岩裸露,节理发育,且时常发生块体滑落或局部滑塌,失稳趋势明显,对该区域造成了极大安全隐患。坡体主要由条带状砂页岩和砂岩组成,上覆0.6 m 左右岩体强风化严重,岩体总体呈东西走向,坡高约48 m,边坡角为76°～86°。采用无人机多层次全方位摄影测量采集边坡结构面信息,设置合理航道,通过三维摄影获取结构面信息,如图1所示。

图1 多层次全方位摄影测量Fig.1 Multilevel and omnidirectional photogrammetric measurement

其中多层次指在距离边坡面1、3、5 m 处分别构建无人机飞行平面;全方位指在各个方位对边坡体进行摄影测量,区别于以往地面定点摄影测量,用覆盖边坡体的百分比来表征该变量。相较于传统测量方法,多层次全方位摄影测量可以最大程度地还原带有坐标的实物模型,并运用计算机算法获取结构面信息,精度高,对于复杂边坡体结构面信息采集有很强的鲁棒性。多层次全方位摄影测量精度见表1。

表1 多层次全方位摄影测量精度Table 1 Precision of multilevel omnidirectional photogrammetry%

由表1 可知:无人机多层次全方位摄影测量精度与多层次量及全方位率成正比,当多层次量为3、全方位率为100%时,与全站仪测量的相对误差小于1%,满足计算要求。

1.2 结构面信息计算

根据结构面产状倾向、倾角概念可得,其产状由从属面的法向量控制,故任意结构面可表示为

式中,A、B、C分别指代结构面参数,且从属面的法向量n可由(-A,-B,1)表示。

根据三点成面原则,将不在同一直线的n个点(n>3)提取可得一结构面,结构面可表示为

则[A,B,C]T可表示为

进而确定其产状信息倾向、倾角分别为

边坡实际点云坐标输出结果及出露迹线提取结果如图2、图3所示。

图2 点云坐标输出结果Fig.2 Output results of point cloud coordinates

图3 出露迹线提取结果Fig.3 Extraction results of exposed trace lines

根据该计算过程所得结构面产状信息及其与人工测量所得信息对比结果见表2,对比结果显示两者误差均在0.70%以内,反映出本研究算法具有良好的精度。

表2 结构面产状信息的无人机摄影测量与人工测量结果对比Table 2 Comparison results of structural plane information between UAV photogrammetry and manual measurement

2 考虑能量耗散的DDA 法

2.1 DDA 基本原理

任意块体i中任意点(x,y)处的位移(u,v)在两个维度上与6 个位移变量相关,即:

式中,(u0,v0)为块体内特定点(x0,y0)的刚性位移;r0为块体绕特定点(x0,y0)的转动角度;εx,εy为该块体的正应变;γxy为该块体的切应变。

对于块体内任意点(x,y),其位移(u,v)可以表示为

以最小势能原理为基础,组建的块体系统整体平衡方程为

式中,Kij是6×6 子矩阵;Kii为块体的几何与物理参数;Kij(i≠j)指块体间接触参数;Di和Fi为6×1 子矩阵;Di为块体i的变形参数(d1i,d2i,d3i,d4i,d5i,d6i);Fi为块体i上分拨给6 个变形参数的应力矩阵。

2.2 能量耗散

岩体破坏过程是不可逆过程,该过程存在岩体变形损耗的能量及结构面相对滑动所损失的能量等,即过程中伴随着能量耗散。非连续变形分析过程可施加适当的黏性阻尼来完成能量耗散机制的设定,并通过经典斜坡案例测试动量耗散阈值,岩体参数取值见表3,边坡模型如图4所示。与工程实验室试验动态位移结果对比,结果(图5)显示:在设置2.45%～2.53%黏性阻尼时与室内试验结果最为匹配,故选取2.5%黏性阻尼作为能量耗散阈值,动力系数k取0.985,用以完成能量耗散设定,使得边坡动态破坏过程与实际破坏过程更加吻合。

表3 边坡模型岩体参数Table 3 Mechanical parameters of the slope rock

图4 滑块沿斜坡滑动模型Fig.4 Model of sliding block along with slope

图5 不同黏性阻尼下DDA 与试验位移结果对比Fig.5 Comparison between DDA and experimental displacement under different viscous damping

通过上述计算分析,考虑能量耗散的DDA 程序运算流程如图6所示。

图6 DDA 程序运算流程Fig.6 Calculating flow of the DDA procedures

3 露天矿岩质边坡非连续变形分析

3.1 计算模型及参数

将无人机采集的结构面信息输入至DDA,边坡体由结构面分割成块体,作为DDA 计算的基本单元,DDA 具体数值计算分析模型如图7所示,其中1#、2#孔为工程现场测斜孔,内置测斜仪用以获取边坡位移,边坡模型岩体参数取值见表4。

表4 边坡模型岩体参数Table 4 Mechanical parameters of the slope rock

图7 露天矿岩质边坡DDA 数值计算分析模型Fig.7 Numerical calculation and analysis model of DDA for open pit rock slope

3.2 现场位移监测及边坡非连续变形分析

测斜孔布置参考图7,测点间隔为1 m,1#、2#测斜孔内不同深度所发生的累积位移如图8所示。

图8 测斜孔测点累积总位移值Fig.8 Cumulative displacements of monitoring points of the inclination hole

由图8 可知:1#孔附近深度为0～12 m 的岩体位移随时间变化较大,深度12 m 附近以下的岩体变形量较小,且趋于稳定状态,说明1#孔内0～12 m 深度附近岩体已发生相对滑动,随时有发生失稳的可能,且滑动面就在该孔深度为12 m 附近位置。随着时间的变化,2#孔内前8 m 深度左右测点总位移呈现规律式上升,且越靠近上坡面的岩体位移量越大,深度8 m以下岩体位移趋于稳定,说明滑动面位置在深度8 m附近,且根据位移上升的规律推断该位置岩体存在倾倒破坏可能。

为预测失稳面确切位置及相应安全系数,并采取相应措施治理该露天矿边坡,本研究运用考虑能量损耗的非连续变形分析方法进行计算,具体动态失稳过程如图9所示。

图9 公路岩质边坡非连续变形分析过程Fig.9 Discontinuous deformation analysis of highway rock slope

由图9 可知:在到达第1.67×106时步前,边坡整体处于稳定状态,当到达第1.67×106时步时,即以折减系数为1.167 参与程序运算时,边坡触发变形失稳过程,各块体发生不同程度位移,且后续变形愈演愈烈,最终发生图9 中所示的大型滑坡,该边坡安全系数为1.167,危险性较高。分析可得,本研究方法计算出的滑落面位置与根据1#、2#测斜孔所得数据预测的滑落面位置基本一致,从现场监测数据角度印证了该方法的可行性及精度。

该露天矿岩质边坡案例可归为牵引式滑动破坏模式,滑坡体基本为强风化砂页岩,存在极大失稳可能性,预失稳岩体均为坡体表层4 m 以内岩体。考虑上述因素及实际情况,建议按照本研究方法计算得出的滑落带位置,剥离上覆预失稳岩体,从源头上解决该露天矿边坡的滑坡灾害问题。

4 结论

(1)提出了一种多层次全方位无人机摄影测量获取边坡结构面信息的方法,为边坡结构面信息采集提供了一种新手段。

(2)考虑了露天矿边坡非连续变形分析过程中的能量耗散,当施加2.45%～2.53%黏性阻尼时模拟边坡动态破坏过程与实际破坏过程吻合度较好。

(3)基于多层次全方位无人机摄影测量及考虑能量耗散的非连续变形分析方法,分析了某露天矿岩质边坡的稳定性,讨论了其最危险滑落面位置,并考虑到其安全系数及预失稳岩体体积等因素。建议按照本研究方法计算得出的滑落带位置,剥离上覆预失稳岩体,从源头上解决该边坡的滑坡灾害问题。