岩溶区多软弱层复合边坡变形破坏机制分析

杨鹏程， 李思达， 韩伟歌， 涂新斌， 王彦兵， 金永军， 周鲲

(1.国网湖南省电力有限公司经济技术研究院， 长沙 410004； 2. 能源互联网供需运营湖南省重点实验室， 长沙 410004； 3. 湖南经研电力设计有限公司， 长沙 410007； 4.中国科学院地质与地球物理研究所 中国科学院页岩气与地质工程重点实验室， 北京 100029； 5. 国网经济技术研究院有限公司， 北京 102209)

针对边坡稳定性分析目前常用的方法有极限平衡法以及有限元数值模拟方法[1-2]。其中，传统的极限平衡法需要提前确定边坡滑移面的形状及位置，这对于含岩溶等不良地质体的复杂边坡而言实现难度较大。而有限元数值模拟方法则以强度折减法为主。Zienkiewicz等[3]和Matsui等[4]提出了强度折减法的概念及其实现原理，并将其与极限平衡法进行对比发现两者获取的安全系数值较为接近。此后，由于强度折减法实现方法简单，并且可获取可靠的安全系数以及边坡潜在破裂面，在边坡稳定性分析方面得到了快速发展[5-7]。强度折减法同样适用于复杂多软弱层边坡失稳分析，已有研究表明层理界面是岩石破裂失稳的主要控制因素[8]，而对于多软弱层复合边坡失稳过程中岩层界面的角色定位需要进一步研究。同样，坡体内岩溶体的存在也改变了坡体的结构及完整性，导致边坡受力更加复杂，研究难度更大。针对含岩溶边坡问题，国内外学者也从理论上开展了研究工作，吴峰等[9]对含岩溶边坡受力进行了分析；江学良等[10]探讨了地下洞室对边坡的影响规律，并建立了边坡下伏洞室计算模型。岩溶导致的破裂主要为地面沉降，因此，白日升等[11]研究了不同覆盖层厚度下的岩溶区地面塌陷规律。靳红华等[12]采用强度折减法分析了岩溶盖层稳定性。而不同溶洞形态分布及治理方法也被大量研究[13-15]。熊超等[16]针对岩溶区基坑开挖导致的层状边坡破裂机制进行了分析研究，探讨了层状边坡变形规律，但其并未考虑岩溶对边坡的影响。目前，对于含岩溶边坡失稳机制研究较少，尤其是边坡在多软弱岩层界面与岩溶共同作用下的变形失稳机制鲜有报道。

现针对某岩溶区变电站场区附近多软弱层复合边坡，采用有限元强度折减法开展含岩溶复合边坡变形破坏机制研究，为更直观地分析岩溶对边坡失稳破坏的影响，分别研究多软弱层复合边坡在不含岩溶时、岩溶无充填时以及岩溶被软弱黏土充填时的变形特征，揭示复合岩层界面以及岩溶对边坡变形的影响机制，提出含岩溶复合边坡治理建议，研究结果对含岩溶边坡支护设计具有指导意义。

1 强度折减法基本原理

强度折减法基于岩土体黏聚力和内摩擦角参数，采用摩尔-库伦(Mohr-Coulomb)准则进行单元弹塑性状态计算，当岩土未到达屈服时，通过折减公式对强度参数进行折减，然后采用折减后的黏聚力和内摩擦角重新计算岩土体弹塑性状态，从而得到新的弹塑性区分布。通过对强度参数的反复折减，直到岩土体发生失稳破坏，此时的强度折减系数则为边坡稳定系数。强度折减公式为

(1)

(2)

式中：c、φ和cn、φn分别为折减前后的黏聚力和内摩擦角；F为强度折减系数。

(2)加固整治取得了重大防洪减灾效益。近年来除险加固的堤防、水库在抗洪减灾中发挥了重要作用，减灾效益十分明显。据统计，全省防洪减灾效益达222亿元，避免36座县级以上城镇受淹。其中，水库减灾效益75.7亿元，全省各类水库拦蓄洪量52.2亿m3，减免农田受灾面积54万公顷，减免受灾人口604.7万人。

采用强度折减法进行边坡稳定性分析时，边坡失稳破坏的判识方法至关重要。目前研究结果均表明以边坡特征点的位移拐点作为边坡失稳判据可得到误差较小的安全系数[17]。

2 工程概况

2.1 地形地貌

拟建变电站工程站址位于张家界市永定区，站址区域属丘陵地貌，地面高程452.71～493.23 m，最大高差约40.52 m，整体呈中间高，四周低，地形起伏较大。场地碳酸盐岩石地层下伏于第四系土层、风化砂岩之下，勘察期间多个钻孔均揭露有溶洞发育，场区岩溶分布不规律，属岩溶强发育地段，通过钻孔发现其内部多被灰黄色流塑-硬塑性黏土充填，物理力学性质较差。

2.2 工程地质条件

综合区域地质资料，进行现场勘察，获取钻探数据，得到拟选站址区域地层分布情况。该区域主要有第四系土层、奥陶系砂岩、泥灰岩、石灰岩等，岩性特征分述如下。

(1)黏土：灰黄色，可塑-硬塑，无摇振反应，捻面粗糙，韧性及干强度中等，分布较广，局部地段底部呈流塑-软塑状态，层厚0.30～13.00 m。

(2)强风化砂岩：下奥陶统分乡组，紫红色，泥质结构，层状构造，胶结程度一般，岩芯呈柱状，锤击无回弹，属较软岩，节理裂隙较发育，岩体较破碎，厚度未完全揭露。

(3)强风化泥灰岩：灰绿色，主要为黏土矿物、方解石，岩芯呈柱状，锤击声哑，无反弹，属较软岩，节理裂隙较发育，岩体较破碎。场地较广泛分布，厚度未完全揭露。

(4)中风化石灰岩：下奥陶统分乡组，青灰色，主要矿物成分为方解石、黏土矿物，块状结构，锤击声较清脆，不易碎，岩体较完整。场地广泛分布，厚度未完全揭露。

拟建站址附近存在多软弱层复合边坡，并且钻孔数据显示该边坡下有两处充填溶洞。因此，基于此边坡分析在岩溶区多软弱层复合边坡变形失稳机制。

3 数值模型

通过区域调查可获取拟研究边坡长度约为85 m，整体坡高约30 m。根据钻孔数据可知该边坡下存在两个充填溶洞，分别位于坡脚和坡体中部。采用三维地质建模技术，对钻孔数据进行光滑插值处理，可获取该边坡下不同岩层分界面，并估算两处溶洞的大致尺寸，从而可建立含岩溶边坡的典型地质剖面，通过AutoCAD软件绘制边坡典型地质剖面图，并将其导出为ABAQUS模拟软件可识别的DXF格式文件，实现CAD和ABAQUS模拟软件的数据对接。通过对地质剖面分区处理体现不同地层岩性的差异化，最终得到边坡剖面数值模型(图1)。该边坡主要由黏土、强风化砂岩、强风化泥灰岩、中风化石灰岩和溶洞五部分构成。

边坡内各岩土层基本物理力学参数通过土工试验获取(表1)，对边坡内五部分区域分别赋予对应材料参数。

图1 边坡数值模型Fig.1 The slope numerical model

表1 模拟参数Table 1 Simulation parameter

假设模型为理想弹塑性，对其采用Mohr-Coulomb破坏准则，并采用非对称完全牛顿进行求解。进一步对边坡模型进行网格划分，使其含12 274个四边形平面应变单元(CPE4R)(图2)。模拟时，固定坡体两侧水平位移以及模型底部的竖向位移，分别设置两个静力通用分析步，第一个分析步用于施加重力载荷，在第二个分析步中定义场变量，并通过二次开发使其控制边坡参数的强度折减，强度折减系数由2减小到0.5。

图2 网络模型Fig.2 Network model

4 边坡变形机制分析

为了更直观地研究含岩溶边坡变形破坏机制，首先进行不含岩溶时多软弱层复合边坡变形失稳分析，从而方便与含岩溶边坡进行对比研究。

4.1 不含岩溶时复合边坡变形机制分析

首先，将坡体内岩溶参数设置为中风化石灰岩参数，即此时不考虑岩溶的影响。进行数值模拟，获取不含岩溶时复合边坡变形过程。ABAQUS模拟软件中可通过观察边坡塑性区分布情况来反映边坡失稳状况。因此，提取不同时刻的塑性区分布图(图3、图4)，由塑性区分布演化可以发现，在强风化砂岩和强风化泥灰岩岩层界面处最先出现局部剪切位错滑移，随后塑性区沿着岩层界面处不断向坡体上下部扩展，最终连通坡脚与坡顶，形成圆弧形潜在滑移面(图4)。由于坡体最上部黏土层较薄，在边坡滑移破裂中并未起主导作用，而第二层强风化砂岩层较厚，成为滑坡体主要组成部分。强风化砂岩与强风化泥灰岩的岩性差异使得岩层界面处最易产生滑移变形，形成潜在滑移面。采用数值模拟方法确定边坡潜在滑移面位置可为后续边坡防护设计提供重要理论支撑。

基于ABAQUS强度折减法不仅可获取边坡潜在滑移面位置，并且可得到边坡的安全系数。选取坡脚处节点作为特征点，提取特征点的水平向位移与强度折减系数，获取强度折减系数与水平位移关系曲线(图5)，由图5可以发现模拟初期坡脚处的水平位移保持不变，直到强度折减系数达到1.77时，位移发生突变，此时边坡失稳。因此，曲线拐点处的强度折减系数即为边坡的安全系数。而该边坡在不含岩溶时的安全系数为1.77，较为稳定。

图3 初始塑性区分布Fig.3 Initial plastic region

图4 潜在滑移面Fig.4 Potential slip surface

图5 强度折减系数与水平位移关系Fig.5 The relationship between strength reduction coefficient and displacement

4.2 含岩溶复合边坡变形机制分析

针对该站址区含岩溶边坡，开展含岩溶多软弱层复合边坡失稳机制分析，并与上述不含岩溶边坡进行对比分析，以便更准确获取岩溶对边坡失稳影响机制。首先获取含岩溶边坡塑性区分布图(图6)，由图6可以发现，坡体中部岩溶的存在直接影响边坡塑性区分布，从而影响到边坡潜在滑移面。而坡脚处的岩溶由于上覆岩土层较薄，岩溶受上覆岩土层压力较小，并未产生明显塑性分布。在该算例中坡体中部溶洞内出现塑性，而由于边坡体存在多软弱层，因此在岩层界面处同样出现塑性区，塑性区集中在坡体下部，而岩溶的存在抑制了塑性区向坡体上部扩展，边坡不再形成完整的圆弧形潜在滑移面。因此，此时滑坡没有大范围潜在滑移的风险，岩溶决定了边坡的稳定性，而此时岩溶在边坡整体稳定性中一定程度上具有积极的影响。

为分析产生上述现象的原因，分别提取不含岩溶时与含岩溶边坡位移云图(图7)。由图7可知，当边坡不含岩溶时坡体位移有沿着潜在滑移面向下发展的趋势，这同样也是潜在滑移面形成的原因。而对于含岩溶边坡，由于岩溶的存在产生局部地面沉降，导致位移云图有向溶洞发展趋势，因此制约了岩溶上方坡体沿着潜在滑移面下滑变形，边坡上部的位移量不足以产生塑性滑动面，从而导致边坡没有产生完整的圆弧形潜在滑移面。

进一步提取岩溶体上表面不同位置处沉降量，获取沉降曲线(图8)。由沉降曲线可以发现，岩溶顶部产生了不均匀沉降，越靠近坡脚沉降量越小，而靠近滑坡顶部沉降则呈现先增大后减小的趋势，在坡体中部位置存在最大沉降量。结合图7位移云图可以发现，由于滑坡体有沿着潜在滑移面向左下方下滑趋势，而岩溶有向正下方沉降趋势，两者综合作用下导致了沉降的不均匀分布。因此，岩溶顶部会受到上部滑坡体的影响产生不均匀沉降。针对这一特点可进行针对性的治理。

图6 塑性区分布Fig.6 Plastic region

图7 位移云图Fig.7 Displacement nephogram

图8 岩溶顶部不同位置沉降量Fig.8 Settlement at different positions on karst top

4.3 无充填岩溶对边坡变形影响分析

拟建站址场区处于强岩溶发育区段，岩溶在场区内分布不均匀，其中大部分岩溶均被黏土充填，但也存在少量无充填岩溶。因此，研究无充填岩溶对边坡稳定性的影响至关重要。基于上述模型，将坡体中部岩溶内填充物去除，设置成无填充岩溶，进行含岩溶空洞的边坡失稳分析，得到位移云图及塑性区分布图(图9)。

由图9可以发现当边坡内含有无充填岩溶时，主要破裂模式为坍塌破坏，岩溶区域出现强烈的地面沉降，导致上部坡体直接垮塌。而在岩溶左侧边界处产生塑性区，即潜在破裂面。上部坡体沿着该破裂面沉降垮塌，而在岩层界面处则存在较小的顺结构面滑移区，此时坡体以垮塌破坏为主，伴有少许滑移破裂，呈现复合破裂形式，无充填岩溶对边坡失稳起主导作用。

图9 无充填岩溶边坡位移云图及塑性区分布Fig.9 Displacement nephogram and plastic region distribution of karst slope without filling

5 含岩溶复合边坡治理建议

由上述研究结果可以发现：边坡与岩溶两者处于相互制约相互影响的关系中，边坡的存在影响了岩溶顶部的不均匀沉降，而岩溶的沉降又反过来影响坡体的潜在滑移面。当岩溶无充填或充填物强度较弱时，将产生较大沉降量，此时上部坡体会产生大范围垮塌，此时边坡呈现崩塌和滑移复合形式破坏。而岩溶内充填物强度较强时则对边坡有较小的影响。因此，边坡与溶洞的关系中存在一个平衡点，在某种特殊情况下，坡体中后方的岩溶对滑坡体的稳定将有积极性的影响。即岩溶的沉降量在允许范围内，并且允许范围内的沉降量可有效减少中后方坡体向坡脚滑移的位移量，而向坡脚位移量的减少又有利于边坡的稳定。当然，这种情况受到岩溶大小、位置、内部充填物强度的影响。因此，在处置含岩溶边坡时，可首先结合三维地质建模技术与数值模拟方法对边坡进行稳定性分析，确定拟研究边坡与溶洞是否可达到平衡关系，如果可以则可采用注浆方式治理岩溶，通过向岩溶内注浆控制沉降量，使其达到该平衡点，此时可大幅度节省工程成本。而当岩溶注浆填充强度过高时，则岩溶对边坡的影响可忽略不计，此时主要考虑不同岩层界面对滑坡体稳定性影响即可。而治理多软弱层复合边坡最有效的方法则是采用抗滑桩[18]，通过数值模拟方法准确获取潜在滑移面后，可针对性地进行抗滑桩设计，获取最优布桩方式，节约施工成本。

6 结论

采用ABAQUS数值模拟软件，基于有限元强度折减法研究了含岩溶多软弱层复合边坡失稳机制，得到如下结论。

(1)多软弱层复合边坡首先在岩层界面处产生剪切滑移破裂，最终顺着岩层界面滑移贯通，产生圆弧形潜在破裂面。

(2)边坡与岩溶两者相互制约相互影响，边坡导致岩溶顶部不均匀沉降，而岩溶的沉降又反过来影响坡体稳定性。岩溶无充填或充填物强度较弱时，边坡呈现垮塌为主和岩层界面滑移为辅的破裂形式，岩溶充填物较强时对边坡稳定性无较大影响。

(3)理想状态下坡体内岩溶在某些特殊情况下可对边坡稳定性产生积极影响。岩溶引起的地面沉降可减小坡体下滑趋势，若该沉降在允许范围内，则该岩溶的存在对边坡稳定性是有益的。而特定的条件则与岩溶位置、大小、充填物有关，具体影响规律还有待进一步研究。

(4)得益于三维地质建模技术的发展，针对含岩溶多软弱层复合边坡，可有效获取坡体内部岩层分布、岩溶状况，并可获取较为准确的三维地质模型，通过三维地质模型与数值模拟软件的数据对接，可采用数值模拟手段对拟研究问题进行超前探索，从而有效指导后续施工设计，实现经济效益的最大化。