基于BIM技术的岩溶发育区岩质高陡边坡稳定性分析

刘均利，廖恒彬，张炳辉,2

(1. 桂林理工大学土木与建筑工程学院，广西桂林 541004； 2. 桂林理工大学广西岩土力学与工程重点实验室，广西桂林 541004)

0引 言

公路边坡失稳时有发生，严重影响公路的安全运营，尤其是岩溶发育区岩质高陡边坡，岩体溶洞、裂隙强烈发育，地下水丰富，使得岩质边坡的稳定性更为复杂[1-2]。

传统二维地质图纸对特殊地质状况的表达不够清晰，数据的遗漏、丢失、存储问题有待解决[3-4]。目前地质建模方法大多是由二维平面拓展而成，通过相邻剖面之间对应的轮廓线连接网来模拟地质形态，若是遇到复杂地质情况，相邻剖面之间轮廓线的结合往往会出现剖面“打架”的现象[5]。随着工程建筑的大规模发展，边坡工程的空间形态和地质条件越来越复杂，地质勘察的精度和广度受到实际条件的限制，边坡体内的应力状态也存在明显的空间差异性，因此传统的边坡数值模型的稳定性分析结果不足以反映真实情况[6-8]。同时地质工程师用软件模拟分析出来的结果又无法以数字化技术存储和传递，极大地降低了地质工程领域的工作效率，不便于实际工程的应用。这些“痛点”问题长期困扰着诸多学者，成为他们迫切攻克的目标[9-11]。

通过引入BIM(Building Information Modeling)技术，结合地质学、几何学、拓扑学与计算机图形学等领域的相关理论与技术，形成三维地质建模方法，可更直观、全面地展示地质特征与地质现象，推测未知地质信息。通过岩溶发育区岩质高陡边坡BIM模型，既能从宏观上判断工程地质体的特性，又能从微观上揭示地质体的内部构造、空间复杂变化等[12-13]。苏学斌等[14]采用地层/岩性混合地质建模方法，对含铀砂层非连续发育的砂岩型铀矿床地质结构进行了表达，并与岩性建模方法进行对比，验证了混合地质建模方法能够显著提高建模的精度。该方法是通过算法来实现未知地质情况的表达，算法的精确度还有待校正，而且该模型还难以灵活地修改剖面，导入分析软件步骤繁琐。

BIM技术与边坡工程的结合主要体现在构建三维数字地质模型。将地质勘察等各类信息整合到三维模型中，可以更加直观、准确、完整地展示地质情况，同时能对该地区地质情况进行全方位剖析，通过对比找到最危险二维地质剖面，实现对地质勘察资料的数据采集、复杂地质的三维可视化建模、分析与设计一体化应用[15-18]。陶志刚等[19]基于太平矿业长山壕露天金矿，构建复杂地质体高精度耦合的BIM模型，并将模型导入FLAC3D软件，分析西南采场边坡失稳破坏的原因，并提出以恒阻大变形锚索为核心的边坡加固方案，与原支护方案进行对比研究，验证了新方案具有良好的加固效果。

本文提出一种岩溶发育区岩质高陡边坡三维地质建模与有限元数值模拟相结合的方法。首先建立岩溶发育区岩质高陡边坡的三维数字地质模型(BIM模型)，精确反映复杂地质构造，然后采用布尔运算对BIM模型进行多次剖切，根据赤平投影找到最不利剖切面并建立二维有限元模型，最后采用有限元方法对边坡多种工况下的稳定性进行分析。本研究能够更好地反映地质构造，更准确地进行地质表述，并通过数值模拟分析边坡各工况下的稳定性，为公路工程的安全运营提供保障。

1工程背景

桂林市兴安县至阳朔县公路延长线工程K5+420～K5+740段开挖后形成高65 m的边坡，坡度为63°～73°，如图1所示。通过4个工程勘探孔ZK1～ZK4勘察发现，地层从上到下可以分为黏性土、破碎土层、砾岩、灰岩。其中，地表以下5～10 m岩体较破碎，以破碎土层和黏性土为主，属较软～较硬岩；地表以下10～20 m岩体裂隙中等发育，岩体较完整，以灰岩和砾岩为主，为坚硬～较硬岩。各岩层力学参数如表1所示。

图1福利路高陡边坡段示意图Fig.1Schematic Diagram of High and Steep Slope Section of Fuli Road

表1公路高陡边坡岩土基本参数Table 1Basic Geotechnical Parameters of Highway High and Steep Slope

2岩体结构面类型及发育特征

通过钻孔电视勘察，查明各孔内的岩体裂隙发育情况，对4个钻孔内岩体的规模、产状、岩层类型、力学性质等特征进行了详细的调查和研究，如下：

(1)断层：里程K5+620附近发育一条北西走向的挤压断层，倾向北东方向，倾角50°；断层面埋深1～6 m，厚度0.1～0.5 m，充填断层泥，其附近分布一定厚度的断层角砾岩。

(2)层间错动面：受断层影响，主断层面附近岩层发生了局部层间错动，错动距离为5～20 cm，充填有黏性土或岩土碎块，其岩层产状较陡。

(3)构造裂隙：受断层构造运动影响，地表以下5 m范围内岩体破碎；5 m以下的构造裂隙一般发育，构造裂隙倾角较陡，一般大于65°，主要构造裂隙的走向与主断层走向基本一致，为北西方向，产生于断层构造运动的挤压作用，裂隙延伸深度一般在20 m左右，局部可达30 m。ZK3附近的岩体裂隙为岩溶水通道，受溶蚀影响其张开度较大，除ZK3孔附近外，其余地段的构造裂隙张开度较小，裂隙宽度一般不超过3 cm。

(4)岩体溶蚀特征：场地裂隙闭合程度较高，岩层透水性差，岩溶水运动不畅，且地下水位较低，其水位高程一般在140 m以下，岩体溶蚀动力主要来源于季节性降雨，故场地岩层溶蚀程度总体不高，5 m以下裂隙面的溶蚀作用轻微，裂隙面铁质氧化物渲染较为明显。岩体溶蚀主要发生在岩面和岩溶水通道的裂隙面，偶有小型溶孔溶洞发育于ZK3附近，溶洞发育高度为1.0～1.7 m，均为浅层溶洞，充填软塑状黏性土。场地岩溶水通道发育于ZK1孔附近，该处裂隙溶蚀现象明显，如图2所示，其中H为钻孔深度。

图2ZK1钻孔电视勘察图Fig.2ZK1 Borehole TV Survey Map

上述地质情况极有可能发生地质灾害，存在极大的安全隐患，因此需要进一步的分析验算。

3BIM的应用及建模流程

3.1技术路线

主要技术流程为：①依据IFC(Industry Foundation Class)标准，通过构件参数约束和逻辑约束建立三维地质信息模型；②以三维模型为载体，运用布尔运算拾取最不利边坡剖面模型；③将边坡模型导入有限元分析软件，在不同工况下检验其稳定性，通过综合分析提出边坡治理方案；④将各阶段的数据上传至数据管理系统，可在不同维度上对成果进行展示，也可嵌入边坡稳定性监测与预警系统用于公路运营阶段边坡失稳预警。技术路线图如图3所示。

图3边坡BIM应用流程Fig.3Slope BIM Application Process

3.2边坡三维数字地质模型的创建

采用GEO5进行三维地质建模，如图4所示，其过程分成地形建模和地质建模。地形建模是将高程点或等高线数据以*.dxf或*.dwg格式导入，生成地形模型，但该模型不包括岩体结构面及发育特征。

图4地质建模流程Fig.4Geological Modeling Process

在地层建模之前需要对钻孔数据进行检查和预处理，筛选出主要地层进行建模，以消除自动建模逻辑错误造成的不兼容。地层建模的属性信息包括高程点、等高线、坡度、坡向等，在坐标系中将地面形态以及空间位置特征使用数字坐标描述。需要注意的是在剖面建模时地层的深度不宜过小，否则会导致后期网格划分时网格点和网格线距离过小，使分析结果与实际情况有偏差。

地质建模采用钻孔数据和剖面数据混合建模。其中，钻孔数据建模需要将地质勘察数据导入其中，然后设置一个尽可能容纳所有地层的层序控制孔，以便建模时钻孔能够更好地兼容。此外还可添加特征钻孔对未探测位置的地形进行局部精细化设计，然后对其他参与建模的所有钻孔添加虚拟层，使所有钻孔处于激活状态，参与三维建模。剖面建模需要在二维工作界面完成，允许工程师依据工程经验在工作界面自定义地层线和断面线对特殊地质进行表达，使三维地质模型更加精确。决定复杂地质建模精确度的几个关键环节包括：

(1)断层及层间错动：通过软件中地质剖面模块进行表达。将勘察所得具有特殊地质特性的ZK1添加至剖面上，由勘察信息确定断层及层间错动的方向、角度、深度、厚度及填充材料，结合地质工程师对实际地质情况的了解，逐一添加地层多段线，并在添加地层面中将多段线定义为断层或层间错动，对未识别多段线的区域，还可以针对区域指定新的地层。

(2)构造裂隙及岩体溶洞溶蚀：通过生成地质模型模块对已知钻孔的连接顺序进行调整，达到兼容状态参与建模。根据钻孔取样试验对裂隙中的填充物进行岩土基本参数的调整，达到能够参与后期有限元分析的要求。特殊地质表达如图5所示，最终生成的边坡模型如图6所示。

图5断层及层间错动Fig.5Faults and Interlayer Dislocation

图6边坡地质模型Fig.6Slope Geological Model

4有限元分析

根据赤平投影结果对三维数字地质模型研究区域内的剖面进行剖切，观察剖面的地层分布、岩土属性和连续情况，根据地质工程师经验，选择最不利剖面建立二维边坡模型，如图7所示。分别采用OptumG2(以下简称G2)有限元极限分析方法[20-22]和摩根斯坦-普拉斯法(Morgenstern-Price法)[23-25]进行对比分析。

图7最不利剖面模型Fig.7Most Unfavorable Section Model

4.1G2有限元分析方法

有限元极限分析方法采用理想塑性模型，研究材料达到极限状态时的力学关系，通过不断降低材料强度(按照一定的比例降低岩土黏聚力、内摩擦因数)或增大荷载，求出材料的极限荷载或稳定安全系数，同时也能够求解材料中的破坏形态。该方法不需要预先知道潜在滑动面，也不需要知道滑动面上的滑动力与抗滑力，就可以计算出潜在滑动面位置和安全系数。对于实际工程的岩土强度与稳定问题而言，目的就是要求出潜在滑动面位置和安全系数，传统极限分析就能够满足实际工程的计算精度要求，没有必要引入本构关系，从而化繁为简[26-27]。

有限元强度折减法的基本原理是将坡体强度参数黏聚力c和内摩擦角φ的正反切值同时除以一个折减系数Ftrial，得到折减后的c′、φ′值[式(1)],然后作为新的资料参数输入，再进行试算，当计算不收敛时，对应的Ftrial被称为坡体的最小稳定安全系数，此时坡体达到极限状态，发生剪切破坏，同时可得坡体的破坏滑动面位置，以及位移、塑性应力和应变值，能够针对性的分析边坡失稳破坏的原因[28-31]。

(1)

G2有限元分析方法是将有限元极限分析方法和有限元强度折减法相结合。先进行有限元极限分析，得到材料的极限荷载和安全系数，然后对其进行有限元强度折减迭代计算，当边坡进入失稳状态时停止计算，得到边坡失稳破坏的位移、应变、应力云图。2种方法的结合既能减少计算时间，又能得到边坡破坏的位移、塑性应变和应力值，同时满足实际工程与科学研究的需要。

4.1.1 数值分析设置

将最不利边坡模型导入G2后，材料类型采用Mohr-Coulomb模型，输入各岩层基本参数，再对模型各地层进行岩层信息的赋予。根据郑颖人边坡与滑坡工程治理[24]中对有限元法计算模型的边界界定，取坡脚到左端边界的距离为坡高的1.5倍，坡顶到右端边界的距离为坡高的2.5倍，且上下边界总高不低于2倍坡高，建立有限元模型。根据模型和实际地质条件确定边界条件为：模型左右边界施加水平约束，底部边界施加完全约束，模型顶端及边坡部分为自由边界。采用网格自适应的方法建立三角形单元，并对网格进行平滑处理，如图7所示，单元数量为11 423，节点数量为3 861。

对于特殊地质而言，G2中的裂隙可以表征岩体结构面及土体内部裂隙等，但这并不是岩土体的竖向拉张裂缝。对于断层而言，在G2中一般采用剪切节理的方式来模拟，剪切节理是夹在两种材料之间无限薄的另一种材料。对于裂隙而言，在Mohr-Coulomb模型中实体模块选择定义1个或2个面，输入参数。需要注意的是，该步骤的倾角α1、α2为从水平面逆时针转向竖直面的角度，如图8所示。特殊地质面参数如表2所示。图9为模型自适应网格划分。

图8倾角示意图Fig.8Schematic Diagram of Inclination

4.1.2 分析结果

对结果提取水平向位移云图、竖向位移云图、总位移云图及塑性应变图(图10)进行分析，得到边坡的安全系数为0.79(小于规范所定的1.35)，可知边坡处于失稳状态。图10中深色区域即为位移发生区域，颜色越深，位移越大。

表2特殊地质岩土基本参数Table 2Basic Parameters of Special Geology Rock

图9模型自适应网格划分Fig.9Model Adaptive Mesh Generation

图10自然工况下的位移及应变等值云图Fig.10Displacement and Strain Contour Map Under Natural Conditions

从图10(a)可以看出，在自重工况下，破坏面形状大致呈圆弧状，最大横向位移出现在距离坡脚35 m的位置，最大位移达到1 m，位移变形集中于边坡的破碎土层堆积处。从图10(b)可以看出，在Y方向上，边坡整体均存在向下移动的趋势，其中竖向位移的最大值出现在坡顶处，为0.312 m。从图10(c)可以看出，在XY方向上，最大位移达到1.03 m，而且在破碎土层和黏性土之间出现了贯通的塑性区，有相对滑移的趋势，这是重点关注的支护区域。从图10(d)可以看出，塑性应变沿破碎土层、黏性土与砾岩层交界面分布，形成局部塑性应变集中带。结合以上分析可知，边坡在破碎土层和黏性土、角砾岩交界处的位移最大，该处边坡失稳破坏主要受到岩体溶蚀、构造裂隙、层间错动的影响较大，极有可能沿着岩层交界处发生剪切-滑出现象。在钻孔附近的断层和裂隙由于处在谷底位置，没有向外发展剪切-滑出破坏的趋势，在自重条件下没有向外发展的趋势，因此对上部边坡的稳定性影响不大。

4.2摩根斯坦-普拉斯法对比验证

为了验证G2有限元分析的准确性，采用摩根斯坦-普拉斯法对该边坡的安全系数进行稳定性验算。

该方法基本流程与G2有限元分析近似。在搜索最不利圆弧滑面时，需要在边坡模型上绘制任意曲线滑动面，这里需要注意的是绘制任意滑动面时需要贯穿整个边坡，而且要避免所画圆弧出现反倾。由图11可以看出，边坡安全系数为0.8，不满足要求，且安全系数小于规范所定的1.35，说明边坡在天然工况下处于不稳定状态，易发生边坡滑移事故。

综上，摩根斯坦-普拉斯法分析结果与G2有限元分析结果近似相等，验证了有限元极限分析方法的准确性。

图11边坡稳定性分析结果Fig.11Slope Stability Analysis Results

5边坡治理措施的分析

由上述分析可知，在自然工况下的边坡不稳定需要进行处治。根据地质地形特征，结合项目区域地质灾害特点及工程治理经验，采用削坡和锚杆相结合的方式进行处治。

5.1削坡处理

挖方边坡每10 m设置一级平台，平台宽度1.5 m，挖方路段设置碎落台，碎落台宽度1.5 m。前三级使用0.3的坡率进行放坡，后坡全部采用0.5边坡坡率进行放坡。在第一级边坡坡顶处设置被动防护网，用于防治坡面碎石落到行车道上，边坡坡顶处设置被动防护网，用于防治坡顶落石。在距被动防护网5 m处设置截水沟，以防止边坡被雨水径流冲刷，确保边坡稳定。

5.2削坡分析结果

削坡后采用有限元极限分析方法得到边坡安全系数为0.536，仍处于失稳状态，且产生了从第5级边坡到坡顶的塑性贯通区域(图12)，有可能发生剪切-滑出破坏，且最大位移出现在第5级边坡处，因此对削坡后的边坡采用锚杆支护进行加固。

图12削坡稳定性分析图Fig.12Stability Analysis Diagram of Cutting Slope

5.3边坡支护方案

根据削坡后的有限元分析结果，采用如下边坡处治方案：①在第1级边坡坡顶处设置被动防护网，防止坡面碎石落到行车道上；②在距离被动防护网5 m外处设置截水沟，以防止边坡被雨水径流冲刷；③第5～8级边坡范围采用锚杆挂网喷混凝土防护措施，锚杆布置为2.5 m(行距)×2.5 m(排距)，呈梅花状布设，单根锚杆长L=3～18 m(具体视开挖后节理裂隙发育特征确定)，设计锚孔孔径φ90，锚杆入射角(俯角)为20°，锚筋采用φ25钢筋，锚杆孔内灌注M30水泥砂浆或纯水泥浆，钢筋网采用φ6.5钢筋，网度200 mm×200 mm，喷射细混凝土厚度为8 cm，混凝土强度等级为C20，挂网边缘用φ18钢筋，锚固长度为1 m，将钢筋网固定；④局部坡面岩石较完整地段采用C20素混凝土进行5 cm的喷锚防护。

5.4锚杆支护分析结果

支护后的边坡稳定分析结果如图13所示，边坡安全系数达到1.405，满足规范要求。经过支护后边坡的最大X方向位移为0.796 6 m，Y方向位移为0.1 m，XY方向总位移为0.797 m，但在位移云图中可以看出，边坡位移云图整体颜色较浅，实际边坡整体X方向位移只有0.16 mm，Y方向位移为0.38 mm，XY方向总位移几乎为0；同时借助该软件对边坡失稳破坏的模拟，发现边坡在支护工况下没有发生剪切-滑出破坏，变形极小，将位移云图的显示范围调至最大发现，最大位移发生在坡面处，对整体边坡稳定性影响不大。

图13支护后网格划分及位移云图Fig.13Mesh Generation and Displacement Nephogram After Support

可见，锚杆支护极大地约束了边坡的侧向位移，有效平衡边坡的滑动力，限制了有效塑性区域的发展，提高了边坡的整体稳定性。

6结语

(1)结合BIM技术，通过剖面建模与钻孔建模混合方法建立的三维数字地质模型具有以下优势：①该方法是一种新的数据经验融合建模的方法，发挥了有限地质勘察数据的最大价值，弥补了只靠二维剖面或钻孔信息建模的单一性，为地质建模提供了新的思路；同时也解决了传统建模相邻轮廓线之间混乱识别、建模精度等难题；②混合建模方法可以准确表达特殊地质情况，实现对模型任意位置的剖切观察，与传统二维图纸相比，该方法对地质结构的描述及可视化方面都有着巨大的优势，同时也能为三维空间内部信息统一储存、汇总提供了渠道。

(2)三维地质数字模型与有限元数值分析方法结合，使得三维数字地质模型在表达特殊地质情况的同时还能够应用于有限元极限分析，模拟边坡的失稳过程，通过摩根斯坦-普拉斯法验证有限元极限分析的准确性，提高有限元极限分析方法的可靠程度和计算效率，对工程应用具有重要的指导意义。

(3)在未知潜在滑动面的情况下，可以通过G2有限元分析方法模拟边坡失稳破坏情况，得到潜在滑动面位置、安全系数、位移及应力应变，分析边坡失稳破坏的原因，提出针对性治理方案。该项目岩溶发育区岩质高陡边坡在自然工况下安全系数为0.79，削坡工况下安全系数为0.536，经过支护加固后边坡的安全系数达到1.405，此时边坡处于稳定状态，位移也由1.03 m减小至0.16 mm，可见此支护方案对边坡稳定性起到显著效果。通过该方法对特殊地质情况的边坡进行分析，提出合理、准确、直观的评价，为边坡设计施工提供了依据，具有广泛的应用前景。