基于无人机调查的露天矿滑坡运动特征数值研究

何旭乾，卢才武，李 萌，陈彦钰，闫雪颂

(西安建筑科技大学资源工程学院，陕西 西安 710055)

露天矿山在开采或自然条件等因素影响下，一些边坡会发生失稳，从而引发滑坡，对矿山开采和人员生命安全造成严重威胁。因此，对滑坡风险定量分析进行研究，对于滑坡灾害的预防治理具有重要意义。随着计算机性能的不断进步，各种数值模拟软件由于自身低成本、兼容强等优点，在滑坡灾害模拟与致因分析过程中应用广泛。颗粒流(particle flow code，PFC)程序是基于CUNDALL等[1]提出的离散元理论，通过Itasca开发的二维模拟程序或三维模拟程序。TANG等[2-3]应用PCD2D分别对1941年和1999年的滑坡的运动特性、裂纹扩展、断裂特性进行分析；LO等[4]使用PCD3D对导致滑坡和破坏的运动学过程进行模拟，评估了滑坡机理、泥石流运动、沉积程度。由此可知，PFC软件能够模拟刚性粒子的运动和相互作用，允许其旋转和移动，适用于滑坡形变破坏、滑移运动、堆积特征的分析。为了反映真实的滑坡运动，曹文等[5]指出了wall-ball法具备颗粒少、节省运算时间、结果准确度高的优点，但是由于PFC软件在滑坡建模这一前期处理方面不具备优势，因而受到一定的局限。

当今无人机技术日趋成熟，其具备的高精度、高效率、低成本、低风险等优点能够很好地解决PFC软件的不足。MA等[6]使用无人机划分滑坡源区、重叠区和堆积区，结合DEM估算滑坡体积并分析DEM运动过程，定量描述了滑坡特征。因此，通过无人机建立高精度数字高程模型(digital elevation model，DEM)，能为PFC软件的滑坡运动模拟提供准确的三维离散元模型。

综上所述，本文对于滑坡运动过程离散元数值模拟主要流程为：首先通过无人机航摄影像建立滑坡区DEM并复现原始地形，构建出滑坡地形(wall)；然后以三轴压缩试验与BP神经网络实现细观参数的选取，对颗粒(ball)进行参数赋值；最后构建完整的ball-wall离散元三维滑坡模型，并从整体和局部的角度对滑坡区域的运动过程、堆积特征进行反演与分析，为滑坡的定量评估提供一套快捷、有效、完整的解决思路。

1 研究区概述

1.1 露天矿区基本工程概况

本文选取位于陕西省境内某露天煤矿作为研究对象(图1)，该露天煤矿地处黄土高原，大部分地区为风沙滩地地貌，总体地势平坦，海拔为1 060～1 230 m。当前采场非工作帮整体边坡角为24°～32°，南端帮为15°，北端帮为17°。坡体岩性由上至下分别为松散沙层、粉质黏土、砂岩，粉质黏土颜色呈浅红色或棕红色，各种中砂、细砂、粉砂为浅黄色或褐黄色。虽然整体边坡处于稳定状态，但是根据现场勘查和对地质资料分析，部分位置在风化剥蚀作用下仍出现了局部滑坡现象。

图1 露天矿现场情况与岩性分布Fig.1 Open-pit mine site situation and lithology distribution

1.2 滑坡现场地质情况

根据对采场现场的调查与分析，造成滑坡的主要影响因素包括节理裂隙、边坡角设计、边坡暴露时长等。该区域内黄土层和红土层垂直节理发育，自然风化会引起局部浅层滑坡，边坡岩体强度随着时间增长呈下降趋势，易产生滑坡。

该滑坡位于采场西北方向非工作帮的地表台阶上，台阶高度为20 m，整体地形主要由粉质黏土构成，采场区内土层垂直节理发育，经过长期暴露，边坡岩体随时间增长出现应力松弛、强度下降的现象，部分坡体向下滑落产生滑坡并形成堆积，根据滑坡后堆积DEM可知(图2)，主要滑动几乎是朝向东南方，滑移区域高约19 m，平均宽度约17 m，滑坡由顶端开始滑落，最终堆积层位于底端。

图2 滑坡区域Fig.2 Landslide area

2 数值模型建立及参数确定

2.1 露天矿边坡三维数值模型构建

本文使用大疆精灵4 Pro无人机对露天矿区进行航测，基于无人机影像数据，使用Agisoft Photoscan生成分辨率矿区的DEM数据。通过Global Mapper对进行影像裁剪处理，选出不稳定边坡区域，并导入Rhinoceros软件中，构建STL网格滑坡模型。在PFC3D颗粒流程序中将模型作为Geometry进行导入，Wall边界设置为模型边界。

在对滑坡区域进行反演的过程中需要对原始地形进行复现。根据DUMAN[7]的研究表明，通过对滑坡体积的估算来建立原始地形的方法具有可行性。在Rhinoceros软件中导入滑坡模型后，生成滑动后地形等高线。通过参考ALOS PALSAR 12.5 m的DEM数据来生成滑坡前地形等高线，在Rhinoceros软件根据等高线划分滑移区域。由图3可知，滑坡区域地形发生了明显改变，滑坡区域呈现下窄上宽，上部宽约23.3 m，下部宽约7.0 m，相对高差19.3 m，滑坡面积约360 m2，体积约为207 m3，1 150 m高程以下的等高线发生明显的弯曲，根据实际调查情况可以看出其下方形成了堆积。

图3 滑坡区域滑移前后的高程线Fig.3 Elevation lines before and after slippage in the landslide area

由于滑移区和堆积区存在重叠，因此在完成地形重建后，可以看出滑坡的顶端较陡，滑移区呈现倒梯形，边坡滑动面较浅，滑移区与堆积区有部分重叠，滑坡顶端与堆积体高差为14 m(图4)。

图4 滑坡区域三维地形Fig.4 3D topography of landslide area

2.2 岩土体的细观力学参数确定

在PFC3D颗粒流软件建立滑坡模型过程中，需要对模型赋予相应的参数，由于材料的细观参数与实际岩土体宏观力学参数没有直接的量化关系式，为了使得模拟结果接近实际情况，需要通过三轴压缩试验，确定合适的细观参数来准确反映实际的岩石力学特征。一般来说，这些细观参数包括颗粒弹性模量Ec、颗粒法向刚度kn、颗粒切向刚度ks、摩擦系数μ、法向黏结强度σ、切向黏结强度c，这些参数对弹性模量E、泊松比ϑ、黏聚力等宏观力学参数具有较大影响。

三轴压缩试验的目的是对岩土体的宏细观参数进行标定，以此获得实验数据。在数值模拟实验中，首先根据以往对岩土体宏细观参数的研究经验，设置初始的细观参数，然后在PFC3D软件内建立三轴试验模型，设置符合滑坡模型的颗粒大小和接触模型，并对初始细观参数进行模拟。通过对应变曲线的分析，可以获取到对应数值的宏观力学参数，将与实际试验结果相近的参数值作为输出样本。

该滑坡体主要为粉质黏土，破坏诱发因素为自然风化与长期暴露。首先，按照岩样尺寸生成边界墙体，依照墙体生成紧密颗粒，指定最大颗粒与最小颗粒的半径比，完成设定后PFC3D系统会按照预设值，将最大值和最小值随机均匀附加在模型的每个颗粒上，最终压缩岩样如图5所示。

图5 三轴试验数值模拟模型Fig.5 Numerical simulation model of triaxial test

为保证生成模型的颗粒尺寸和分布满足实验要求，使得颗粒之间具备较高的接触精度，需要对模型进行均匀化处理。均匀化处理的目的是要保证生成的各相邻颗粒之间的重叠度足够小，让颗粒与边界接触紧密，达到完全耦合的效果。根据图6的均匀化处理，利用自编程序不断调整模型边界，使得颗粒之间接触处于足够理想的状态。

图6 三轴试验均匀化处理Fig.6 Homogenization of triaxial test

接触模型选择为平行黏结模型，该模型可以同时传递力和力矩，受到破坏后刚度下降，与滑坡破坏特性相近。完成接触模型设定后，将一组细观参数值进行赋值，进行轴向加载，完成数据测定后绘制应力-应变曲线，根据测定数据计算宏观力学参数。

为了确定细观力学参数对岩土体宏观力学参数与力学特性的影响，需要对细观力学参数做敏感性分析，以此定性分析宏细观参数间的关系。根据以往学者[9-11]的研究，主要包括黏结强度、摩擦系数、刚度比、颗粒形状的影响。黏结强度决定了细观层面上颗粒间胶结强弱度，从而影响宏观的峰值强度和应力大小，若增大法向平行黏结强度与切向平行黏结强度，会引起弹性模量和峰值强度的增大。当摩擦系数增大时，弹性模量和峰值强度也在增长但涨幅较小。而颗粒的法向刚度与切向刚度之比将对泊松比产生影响，随着刚度比的增大，泊松比也会增大。模拟实验中颗粒形状采用球形，不规则颗粒形状虽然能够更好地模拟现实中岩石的几何形态，但考虑到该类块体形态难以考虑到土石间凹凸不平引起的咬合作用，使得材料模拟结构与实际差异较大引起实验结论不准确，因此后续模拟将不考虑这一因素。

一般获取细观力学参数的方法需要假定细观力学参数，通过PFC3D软件内三轴实验获取宏观力学参数并与实际室内试验数据进行对比来获取符合参数。该方法较为耗时，因此本文采用机器学习方法中的BP神经网络，对上述三轴试验的数值模拟过程获取数据进行训练。具体流程是让随机组合的细观参数，以及模拟相对应宏观力学参数作为输出样本和输入样本，通过机器学习模型进行训练，根据训练样本构建出宏观参数与细观参数之间的非线性映射关系，对细观力学参数进行反演。将10组宏观力学参数选做测试样本来验证模型，最终将实际的宏观力学参数带入模型中对细观力学参数进行反演(表1)，分别获取到最合适细观标定参数。

应用反演精度公式[11]、平均残差公式[12]、残差均方差公式[13]来验证BP神经网络模型的性能，经过验证其精度满足建模要求。将表1获得的细观力学参数导入PFC3D软件的三轴压缩试验中，获得应力-应变曲线(图7)。

表1 模型细观力学参数Table 1 Mesomechanical parameters of model

由图7可知，依照摩尔-库伦强度理论计算滑坡黏聚力和内摩擦角[14-17]，计算同室内三轴压缩实验数据对比，由表2所得的对比值可知，其宏观力学参数值非常相近，证明该方法可行，可以用于构建露天矿滑坡的颗粒流模型。

图7 应力-应变曲线示意图Fig.7 Schematic diagram of stress-strain curve

表2 模拟值的宏观力学参数对比Table 2 Comparison of macro-mechanical parameters of simulated values

3 边坡滑移过程模拟与分析

3.1 滑移运动过程模拟与分析

建立的滑坡三维模型(图8)，滑移区域共生成11 002颗颗粒，颗粒高程差约为19 m。在对颗粒附加细观参数后，施加重力即可进行滑坡运动特征的模拟。滑坡模拟过程中分别从整体区域滑移和区域特征点两个方面分析，在滑坡区域分别选取9个特征点作为监测点(图9)，在模拟中对其速度和位移进行监测，监测点的坐标分别为A(-8，13，10)、B(-3，14，10)、C(3，14，10)、D(-7.5，10，7)、E(-2.5，11，7)、F(2.5，11，7)、G(-5，5，3.5)、H(-1.5，5.5，3.5)、I(2，5.5，3.5)。

图8 原始区域PFC滑坡模型Fig.8 Original regional PFC landslide model

图9 监测点位置与编号Fig.9 Location and number of monitoring points

在将滑坡体的速度和位移清零后，完成颗粒与颗粒、颗粒与墙体间的细观参数赋值，并删除滑移坡面墙体，最后对滑坡模型施加重力，让颗粒在重力作用下进行滑移，其整体运动过程如图10所示。由图10可知，滑坡从失稳到运动停止一共经历20 s。在边坡滑移前期中，1 s时坡面发生明显形变，坡面顶端首先出现形变，坡表面颗粒速度达到了约10 m/s；到2 s时坡体迅速下滑并开始形成堆积，整体坡面继续进行滑移，坡底端颗粒速度达到15 m/s，滑坡整体运动速度达到峰值；在3 s左右随着堆积形成，滑移速度逐渐减慢，底端颗粒在堆积后继续移动，底端局部速度为17.5 m/s；到达5 s时，顶端沉降速度开始减慢，随着滑移、摩擦对能量损耗，滑坡运动速度显著减慢，底端局部颗粒仍旧向前移动。随着大部分颗粒的堆积，在10 s时滑坡体后端逐渐趋于稳定，前中端局部颗粒保持2.0 m/s速度继续移动；在15 s时滑坡体前端整体开始稳定，直至20 s时大部分颗粒基本静止，最终形成完整的堆积体，只有小部分土体以较小速度继续运动。

图10 边坡滑移过程Fig.10 Slope slip process

应用PFC程序自编FISH函数进行滑坡整体的位移与速度变化监测(图11)。由图11可知，滑坡平均速度在快速提升后又急速减缓，逐渐趋于平缓，在1.8 s时达到峰值约为5.6 m/s；平均位移不断增加，但是增长速度逐渐减缓，最终位移距离约为13 m。

图11 边坡滑移平均速度和位移Fig.11 Average speed and displacement of slope slip

完成整体滑移运动分析后，分别对滑坡各区域速度和位移变化进行分析，监测点的颗粒运动速度变化见图12。其中，图12(a)中点C峰值速度最高，约为6.9 m/s，点A、点B的速度波动较为接近，峰值速度略低于整体平均速度，三个颗粒经过10 s都停止运动；图12(b)中三个点峰值速度都与整体滑破平均速度接近，期间点F在3.5 s时有强烈波动，同顶端点运动停止时间相近；图12(c)的底端点运动时间相比于中顶端的点，运动时间和整体滑破运动时间相同，经过20 s后点H、点I静止，而点G还在缓慢移动。

图12 滑坡各部分监测点速度变化曲线Fig.12 Velocity change curve of monitoring points in each part of the landslide

对比各点速度变化图，对图13中各部分监测点位移曲线变化进行分析。其中，图13(a)中点B位移量最多为12.0 m；图13(b)中点E位移量略高于点B，在10 s时达到12.2 m后静止；图13(c)中点H、点I经过20 s后停止移动，位移量分别为13.4 m和14.4 m，而点G经过22 s后位移量还在增加。

图13 滑坡各部分监测点位移变化曲线Fig.13 Displacement change curve of monitoring points in each part of the landslide

通过对各部分颗粒速度位移变化分析可以看出，每个位置上监测点的速度变化中，靠近滑坡右侧的点C、点F、点I的速度波动频率最明显，而点C、点F两点对比相同位置上其他点位移量都最少，说明滑坡右侧更为陡峭，而滑坡左侧地势较为平坦，所以底端点G位移量最多，并在滑坡整体基本静止后仍有微量位移。滑坡顶端和中端经过10 s静止，底端颗粒在形成堆积后继续移动，监测点速度也在快速减缓后再次升高。

完成模拟后，对滑坡整体滑移峰值平均速度进行计算[15]，通过计算得到峰值为6.2 m/s，较高于模拟结果值。考虑到该理论忽略了岩体颗粒间的摩擦、碰撞的动能消耗，因此误差可以近似忽略，可以看出滑坡的反演结果和实际滑坡堆积结果基本相似，证明了实验的准确性。

3.2 滑移堆积过程分析

根据滑坡的滑移运动过程模拟分析，滑坡前边坡最大宽度约为19 m，滑坡的滑移区域和堆积区域有部分重叠。为了更直观展示滑坡的堆积特征，按照高程将这片区域的颗粒划分为不同的组，对分组颗粒进行染色。图14对比了滑坡前后的分组颗粒，能够观察到顶端颗粒、中端颗粒排序并未因滑坡发生而产生明显变化，只有低端颗粒在运动时发生明显混合。其中，顶端的两个堆积区域和滑移区域发生重叠，左侧部分颗粒被掩埋并隆起，中后端形成堆积平面。通过剖面图也可看出堆积主要产生在后端，这与滑坡后实地调查堆积特征基本一致。

图14 分组颗粒堆积分布Fig.14 Grouped particle packing distribution

通过滑坡运动轨迹和堆积特征分析，证明了离散元能够有效地模拟矿山滑坡特征与运动趋势，其数值模拟预测结能为滑坡的防治提供依据。

4 结 论

1) 本文将无人机调查同PFC3D三维数值建模相结合，通过无人机航摄构建的滑坡DEM与复现原始地形，以Geometry Import指令在PFC3D内构建出离散元三维数值模型，为后续滑坡的数值模拟提供了数据支撑。

2) 通过在三轴压缩模拟实验中加入BP神经网络的方法，完成对宏细观参数的反演与标定。经过验证，将获取的细观参数带入PFC3D内获取数据与室内试验对比，结果能够满足要求。

3) 滑坡整体历时20 s，在1.8 s时平均速度达到峰值5.6 m/s。滑坡中顶端颗粒下落过程中呈现聚集状，底端趋向于向左侧运动，对比监测点速度位移曲线运动，滑坡体右端较为陡峭。堆积体主体主要为中顶端颗粒，滑坡与堆积体有部分重叠，与实际堆积特征基本相符，证明本研究的模拟方法可以应用于滑坡的运动特征研究。