煤矿工业广场开挖边坡崩塌特征分析与运动模拟

常一伦，谢婉丽，2，姬文龙，杨 惠，单 帅

（1.西北大学 地质学系，陕西 西安 710069；2.西北大学 大陆动力学国家重点实验室，陕西 西安 710069；3.陕西陕煤陕北矿业有限公司 生态环保部，陕西 榆林 719000）

0 引 言

随着生态文明被写入宪法，绿色矿山理念已上升为国家战略[1-2]。大规模矿山活动在给我国经济社会注入新鲜活力的同时[3]，不合理的矿山开发也破坏了矿区的生态与资源环境[4]，诱发崩塌、滑坡、泥石流、地面塌陷、地裂缝等系列地质灾害[5]。作为矿山灾害多发国之一，我国矿山地质灾害有范围广、种类多、损失重、影响大等特点[6]。崩塌是指较陡斜坡上的岩土体在重力作用下以倾倒、滑动、坠落等突然的方式与母岩分离，并以滑动、滚动、跳跃等方式沿斜坡向下运动，最终在坡脚堆积的地质现象[7-8]。其发生具有灾害频发性、时间不确定性、发生快速、造成危害大等特点。目前，关于崩塌的研究主要集中在形成与破坏机理、勘察设计治理等方面[9-12]。在陕北开展煤矿开挖边坡崩塌特征分析，预测崩塌的运动特征，对已建工程和拟建工程及区内车辆、设备、人员的安全防护具有重要意义。

崩塌的研究基础是现场调查，传统现场调查方式多受地形、技术[13-14]、人员、成本[15-16]等条件制约。近几年迅速发展的UAV技术具有高精度、低成本、易操作、数据易处理等优点，很好解决了这一痛点问题[17]。叶震等[18]利用多旋翼无人机航测获取了边坡三维模型，并对模型进行球形k均值聚类分析，成功提取边坡岩体结构面的倾向、倾角参数。王颂等[19]现场调查结合无人机航拍，成功提取了西藏某崩塌三维数字模型和崩塌危岩区的岩体结构面信息，进行了崩塌模拟，为沿线工程的崩塌防护提供了试验依据。

对崩塌运动特征的研究主要利用数值模拟，常用的二维崩塌数值模拟软件有ROCFALL和Rockfall Alalyst。ROCFALL可模拟崩塌运动路径，提取高度、能量、速度和到达概率等，其应用较为广泛。SAGHIR[20]对Chalk cliffs进行了基于ROCFALL和现场调查的崩塌几何尺寸、崩落点和崩落范围分析。KISTER等[21]选取崩塌发育地区的路堤，进行了ROCFALL和Rockyfor3D的崩塌模拟，并进行了比较。Rockfall Alalyst是一款基于ArcObjects和C#的GIS扩展程序[22]，可在ArcGIS下通过参数设置，利用无人机航测的高精度DEM影像底图进行崩塌的2D模拟，并在ArcScene里3D可视化。仉义星等[23]利用Rockfall Alalyst模拟了崩塌区落石的路径和危险性分析，所获得的危险性物理模型结合易发性统计模型，完成了区域崩塌滑坡的易发性和危险性等级综合分析。然而，ROCFALL是一款二维软件，其剖面选取的人为性使得结果具有很大的偶然性。Rockfall Alalyst虽可以三维可视化，但落石的距离、角度间隔等，也都是人工设置后固定不变的。以上2种基于确定性的软件都无法实现更接近真实的落石三维模拟运动。RocPro3D是一款专业三维崩塌模拟软件，使用密集点云生成的三角网格作为地形底图，可模拟落石在三维空间的运动，其功能涵盖：地面模型生成、地面介质属性定义、落石参数定义与运动模拟、数值分析与成图等，广泛应用于崩塌、泥石流等模拟与危险性评估中[24]。

以张家峁煤矿工业广场开挖边坡崩塌为研究区，调查历史崩塌的分布和危岩区危岩体特征，利用无人机进行航摄测量，获取DEM、和点云数据，进而获取岩体结构面信息。结合野外统计数据和无人机航拍，采用RocPro3D对危岩区进行三维数值模拟。旨在反演研究区崩塌模拟的各项参数，预测落石的运动轨迹、影响范围、弹跳高度、能量分布等，以期为崩塌治理、防护设计提供参考依据。

1 工程地质条件

张家峁煤矿位于神木市北30 km，坐标N 39°1′15.17″，E 110°24′24.13″。场地建设过程中对植被的破坏和原地形开挖，加上差异风化、卸荷作用等，崩塌灾害尤为严重。崩塌所在边坡坡底标高+978 m，坡顶标高+1 070 m，走向355°，宽度约500 m，长度约140 m，坡度约35°，岩体破碎度高，节理裂隙发育，植被盖度低。现有防护措施包括局部主动防护网和被动防护网、挡墙、排水渠处理（图1），为防止再次发生崩塌，进行了现场调查和无人机航拍，将东边坡分为A～F 6处典型危岩区（图2）。分别距坡脚80、80、67、55、27、31 m。坡体和坡脚均有不同数量的历史崩塌堆积。

图1 工业广场东边坡历史工程治理措施Fig.1 Treatment measures of historical engineering for east slope of industrial site

图2 工业广场东边坡崩塌概况Fig.2 Overview of collapse of east slope

2 崩塌点基本特征

2.1 历史崩塌的分布特征

东边坡历史上曾多次发生崩塌，坡脚和坡体上堆积了大量落石。据现场调查将坡脚和坡体上落石划分为3个主要区域：北侧、中南侧和南侧。

坡脚采用人工现场量测（图3），坡顶等难以到达的地方借助高精度DOM（分辨率4.5 cm/pix）进行基于颜色、形状、大小、纹理等的人机交互目视解译（图4），多次量测取均值方式，统计体积不小于0.1 m3的落石数量分别为：73、201、75块（图5）。按落石与边坡的位置关系和粒径区间来统计落石数量，结果如图6所示。

图3 落石粒径现场调查Fig.3 Field survey of rockfall size

图4 落石粒径目视解译Fig.4 Visual interpretation of rockfall size

图5 历史崩塌堆积分布Fig.5 Distribution of historical collapse

图6 历史崩塌落石尺寸-数量分布特征Fig.6 Size-quantity distribution characteristic of historical rockfalls

由图6a知：落石崩落粒径主要集中在0.1～5.0 m3，占落石总量的97.4%；中南侧区域落石分布最多，共201块，占比57.6%，北侧、南侧相差不大，占比20.9%和21.5%。由图6b知：坡体上落石，崩落块度在0.1～0.5 m3最为集中，占坡体上落石比64.0%，这是因为初始速度一致时，体积更小的落石崩落后动能小，易受坡体摩擦与障碍物阻挡而最先滞留坡体上部，体积中等的具有更大动能到达坡体下部，另外受崩塌物源的节理面间距、沉积岩层厚度的限制，体积大于5～10 m3的落石本身即较少，因此无论坡体上还是坡脚都数量最少。由图6c知：坡脚落石，崩落块度在0.1～1.0 m3最为集中，占比79.7%。崩塌的运动方式主要有滑动、滚动、跃动和自由崩落等，受限于节理面间距、沉积岩层厚度，体积大于5 m3的崩塌大多磨圆差，不易滚动，大概率只能以较小的初始速度滑动，难以到达坡脚。总体而言，落石崩落块度分布区间，主要受控于结构面数量、展布特征及间距；坡体上的具体概率分布，则受其粒径大小、形态、磨圆、脆性、坡表性质（障碍物特征、介质特征）、坡度和微地形等共同作用。

坡脚和坡体上落石相对位置与粒径大小的分布规律，可以为崩塌模拟中的地面介质参数反演提供一定数据参考。

2.2 岩体结构面产状

由于研究区所在的边坡高度高，坡度大，难以人工进行结构面产状测量。利用无人机航摄和内业处理得到三维点云数据，通过Coltop软件可获取危岩体结构面产状（图7）。

图7 结构面信息提取统计Fig.7 Structure information extraction statistics

调查发现危岩区主要受3组结构面控制，从6处危岩区分别均匀选取若干产状平直光滑的结构面，每组结构面提取200个点的XYZ坐标值，并通过计算得到每组结构面的200个产状，将这些产状信息作图成极射赤平投影极点图（图8）。

J1:结构面1；J2：结构面2；J3：结构面3图8 岩体结构面极射赤平投影图Fig 8 Polar stereographic projection of rock structure surface

如图8所示，每组结构面产状数据相对集中，离散程度小。计算每组结构面产状的平均值，可获得6处危岩区对应的3组主要结构面产状平均值，见表1，并在图8中作出3组结构面产状平均值对应的极射赤平投影。

表1 危岩区岩体结构面产状统计

2.3 物源区危岩体特征

受岩体结构面的切割，6处危岩区均为块状结构，结构面间距直接控制了危岩体的块度。通过无人机航摄处理获取的正射影像图，统计每组结构面间距，并计算对应各危岩区的危岩体积（表2）。

由表2知，E、F处危岩体块度最大，A、B和D处较小，C处居中。由前面小节统计数据可知，北侧、中南侧、南侧3个区域的坡体上、坡脚的岩体优势粒径分别在0.1～0.5、0.1～0.5、0.1～5和0.5～5、0.1～1、0.1～0.5 m，而A～D 4处危岩物源区——除E、F 2处点状危岩区（尚未发生崩塌）外，优势体积分别为2、2、7.5、3.75 m3。

表2 危岩区主要结构面间距与体积

由图9可见历史崩塌体积在边坡的整体分布主要集中在0.1～5.0 m3，与物源区危岩体结构面控制下的平均体积具有很强关联性，说明历史崩塌最主要源自A～D危岩区。

研究历史崩塌的体积和空间分布特征，可以在崩塌模拟的调参过程中为参数范围的确定提供参考，而物源区危岩体的体积分布与展布特征，可以为数值模拟中物源区几何参数的选取提供更可靠的演示数据。

2.4 危岩体失稳方式

图9 历史崩塌优势体积在坡体的分布概况Fig.9 Distribution of dominant volume of historical collapse on slope

危岩体发生破坏的原因主要是由于主控结构面的端部岩石破裂所致[25]。目标区危岩体主要受控3组结构面而被切割，且每组结构面方向与岩体临空面方向一致，因此可根据岩体结构面之间的切割关系，初步进行危岩体的稳定性判断。如图8所示，如选取J1结构面作为临空面，则J2和J3交于A点，交线OA倾向与J1结构面倾向一致，交线倾角大于J1结构面倾角，说明危岩区在J1临空面是相对稳定的；如选取J2组结构面作为临空面，则J1和J3交于B点，交线OB倾向与J2结构面倾向一致，交线倾角小于J2组结构面倾角，说明危岩区在J2临空面是不稳定的；如选取J3结构面作为临空面，则J1和J2交于C点，交线OC倾向与J3结构面倾向一致，交线倾角小于J3组结构面的倾角，说明危岩区在J3临空面是不稳定的。但由于沉积层差异风化作用和其他次级结构面，危岩体也存在以J1结构面为临空面失稳，发生崩塌的可能存在。

倾向J1临空面的危岩体受岩体结构面切割后呈块状、板状，危岩体近乎水平，在危岩体下部抗风化能力较弱岩层发生掏蚀后，容易发生错断式崩塌。倾向J2和J3临空面的危岩体倾角多在65°～75°，受结构面切割呈高角度倾斜块状，在重力或其他外力作用下，有沿着临空面发生向下滑移、或沿着危岩体底部某一支点倾倒的趋势。据现场调查历史崩塌分布和无人机航拍，A～D 4处带状危岩区现状欠稳定，推测主要以J2、J3结构面失稳，E、F 2处点状危岩现状欠稳定，其下部岩层差异风化严重，推测主要以J1结构面失稳产生错断式崩塌。

危岩体失稳后，按所发生崩塌体的运动过程，可分为滑动、滚动、跃动、自由落体或复合等多种形式[26]，主要与危岩体的密度、磨圆、块径、高差及坡形、坡度、坡面粗糙度等有关。同一物源区发生2次崩塌，其运动路径、能量等也绝非一致，故需要借助数值手段进行模拟。

3 运动学模拟

3.1 参数反演

采用RocPro3D进行崩塌数值模拟时主要得确定2类参数：危岩体特征参数和岩土体介质表面参数。其中危岩体参数主要有shape、d、d/h、density、v和H，shape为危岩体形状，d为柱状危岩体截面直径，d/h为柱状危岩体直径与柱状危岩体母线之比，density为危岩体密度，v为危岩体初始速度，H为危岩体距离坡体表面的初始高度。岩土体介质表面参数主要有Rn、Rt、k和β，Rn为法向恢复系数，Rt为切向恢复系数，k为摩擦系数，β为滚动与跃动间的转化角。

据不同性质对岩土体进行类型划分（图10），A～F表示危岩区，①～⑥表示原始基岩面、开挖斜坡风化基岩面、生活区、G338国道、黄土路肩、河流水面，不同岩土体区域其参数取值不同。据表2和DOM测量，给出各危岩区的危岩体特征参数取值见表3。

表3 危岩体特征参数选取

首先根据Rockfor3D软件指导[27]及前人经验[28-29]选取岩土体介质参数的取值范围，然后在取值范围内试算，获取崩塌模拟的空间分布特征，最后对比模拟结果和历史崩塌分布特征，确定危岩体特征参数和岩土体介质表面参数的合理取值。最终确定参数见表4，模拟得到危岩区崩塌的空间分布如图11所示。

表4 岩土体介质表面参数设置

图10 岩土体介质及危岩区划分Fig.10 Geotechnical media and dangerous rock area division

3.2 崩塌运动分析

根据反演结果，确定模拟次数、输入的危岩体体征参数及岩土体介质表面参数，对A～F危岩区分别进行计算，分析每处危岩区的崩塌运动特征。

获得危岩区落石到达概率如图12所示，据各危岩区崩塌落石到达概率大小，作出各危岩区崩塌落石最优运动路径（图13），并给出最优路径下，各危岩区落石与坡面间的关系（图14），从而量化危岩区发生崩塌落石对其下方各区域的影响，统计结果见表5。

从落石空间分布来说，A区崩塌主要从坡体正前方通过，大概率对生活区形成威胁，最远也仅到达生活区内，不会对停车场、G338国道、黄土路肩和河流水面构成威胁；B区崩塌主要从坡体正前方通过，C、D区崩塌主要从坡体正前方和侧沟谷通过，B、C和D危岩区大概率对停车场和少部分生活区形成威胁，最远可到达靠近G338国道的边缘处，不会对G338国道、黄土路肩和河流水面构成威胁；E、F区大概率对停车场形成威胁，最远也仅到达停车场内，不会对生活区、G338国道、黄土路肩和河流水面构成威胁。

图11 危岩区崩塌空间分布Fig.11 Spatial distribution of collapse in dangerous rock area

从运动规律来说，落石从危岩区向下运动过程中，随着落石与坡面的碰撞，损失的势能转化成内能、热、声能等。在坡降大的位置，落石距离坡面的竖向距离也大，碰撞后弹跳的高度往往也较大。落石距离坡面的竖向最大距离多位于坡体的中下部，其势能逐渐减小，动能先增后减。从坡降较小到坡降较大的位置，落石与坡面间的竖向距离持续增加，其动能通常也随着下降运动，在发生下一次碰撞前的最后一刻，达到峰值。以B区落石最优运动路径为例，落石平均尺寸2 m×1 m×1 m，通过坡面和停车场时，其落石的质心距离坡面的距离6.39 m，距离停车场地面距离6.5 m，对停车场的冲击动能可达1 600～1 800 kJ，因此会对其停车场构成威胁。

据数值模拟结合现场调查的历史崩塌分布，6处危岩区物源丰富明确，节理发育，岩体风化破碎程度较高，属于高易发区，将持续发育崩塌。从模拟结果来看，A危岩区落石最优路径下，其下方的被动防护网可以起到一定的防护作用，B危岩区落石最优路径下，其下方被动防护网与落石在防护网处的高度基本持平，难以绝对起到防护作用，C、D、E、F等4个危岩区下部除有少量树木，无其他防护，有必要进行防治。可利用数值模拟的结果为其防治提供依据。据图8—图10模拟所得崩塌落石影响范围、堆积位置、最优运动路径，可明确防护设施的布设范围。据图14模拟所得最优运动路径下各危岩区的崩塌落石距离目标区的高度和动能，可明确防护的高度、角度和强度。

图13 落石最优运动路径预测Fig.13 Prediction of optimal movement path of rockfalls

图14 最优路径下各危岩区落石与坡面关系Fig.14 Relationship of rockfall and slope in each dangerous rock area under optimal path

由图14知，受6个危岩区影响最大的目标区主要是生活区和停车场，G338国道、黄土路肩、河流水面均不受影响，从图14的A危岩区可以看出，在落石最优运动路径下，被动防护网能较为有效地拦截崩塌落石；由图14的B危岩区可以看出，在落石最优运动路径下，被动防护网拦截崩塌落石的几率大幅降低。由于被动防护网仅在生活区东侧坡体中下部有所设置，因此具体的防护措施布设范围应选取从生活区与停车场交界开始的整个停车场区域，并根据模拟结果结合工程投资等，给出本研究区更为合理的崩塌防治方案。

4 结 论

1）野外调查结合无人机航测技术、点云结构面参数获取技术，可以实现对陕煤张家峁煤矿工业广场东边坡崩塌落石的历史崩塌分布、岩体结构面产状和物源区危岩体特征的快速提取，有效克服了传统野外崩塌调查耗时长、难度大、精度低等弊端，为广泛工况下的崩塌地质灾害的识别、预测、防治提供了思路和技术方法。

2）基于无人机自动航摄的外业和内业流程化处理，提取高精度点云和DEM数据，并以现场调查所统计的数据为依据，利用参数反演确定崩塌模拟的各项参数，进行了崩塌的三维数值模拟。结果表明模拟结果与现场较为吻合，为运动学模拟提供了可靠依据。

3）利用三维数值模拟所获得的崩塌落石在单位栅格上的到达路径的概率分布，可预测落石在坡体上的最优运动路径，有效避免了2D崩塌模拟软件选取剖面时的主观干扰因素，显著减小误差，有效的获取了崩塌的影响区域。结果表明，6处危岩区主要威胁工业广场东北侧的生活区和东南侧的停车场，不会对工业广场中部的G338国道、工业广场西侧的黄土路肩和河流水面构成威胁；落石距离坡体的最大高度多集中在坡体中下部，其运动的动能先增大，后减小；位于生活区东侧，坡体中下部的被动防护网能较为有效地防护A危岩区的落石，但对B危岩区的落石防护作用不显著，C、D、E、F等4个危岩区主要威胁停车场无防护措施。以上，可为崩塌的防治设计提供一定的指导依据。