考虑岩体劣化的三峡库区顺层岩质边坡破坏分析

秦 盼 盼，黄 波 林，秦 臻，闫 国 强，白 林 丰

(1.三峡大学 防灾减灾湖北省重点实验室，湖北 宜昌 443002； 2.三峡大学 湖北长江三峡滑坡国家野外科学观测研究站，湖北 宜昌 443002； 3.中国地质大学(武汉) 工程学院，湖北 武汉 430074)

0 引 言

三峡库区广泛发育顺层基岩滑坡，对沿岸居民和长江航道带来巨大的威胁。例如，2003年7月13日发生的秭归县千将坪滑坡，体积约为2 400万m3，产生了约39 m高的涌浪，造成24人死亡或失踪[1-2]。同时，顺层岩质滑坡特殊的地质结构，导致水库蓄水对顺层滑坡的稳定性影响强烈。

Yin等认为三峡水库诱发滑坡可分为3个阶段，不同阶段其主要诱发因素不同，分别是水位变化叠加降雨和岩体劣化[3]。从顺层滑坡典型案例来看，千将坪滑坡是2003年三峡库区首次蓄水达到135.00 m诱发形成的，杉树槽滑坡是在2014年9月1日水位大幅上涨和暴雨期间发生的[4]，卡门子湾滑坡是2019年12月10日发生，由阻滑段岩体强度劣化所致[3]。

在诱发滑坡的前两个阶段，对三峡库区顺层岩质岸坡研究较多，取得了大量认识。例如，在水位变动下顺层滑坡研究方面，刘才华等研究指出库水位上升会使岩土体水致弱化，滑带抗剪强度急剧降低，坡体抗滑力降低，影响岸坡稳定性[5]；邓华锋等研究指出在长期水-岩循环作用下的力学参数弱化，会使稳定滑坡向不稳定方向转变[6]；Jian等研究发现千将坪滑坡失稳的主要诱发因素是库水位上升和持续降雨[7]；代贞伟的研究表明坡体软硬相间顺层单斜结构是藕塘滑坡形成的内在因素，降雨及水库蓄水是导致藕塘滑坡复活的主要诱因[8]。但关于第3个阶段库水位长期波动条件对顺层岸坡变形破坏的影响研究较少。Tang等通过对黄土坡滑坡进行长期现场原位监测，发现坡体位移曲线中的“跳跃点”与暴雨和库区水位快速下降的时间相吻合[9]。黄波林等研究发现岩溶岸坡岩体劣化早期，以破坏浅表、岸坡后退为主，引起的主要破坏模式为崩塌、滑移和倾倒3种类型[10]。刘广宁等也发现受库水位波动影响，坡体前缘受江水侵蚀，坡脚部位的软弱夹层向内凹陷，岩体会发生滑移、错动、坍塌等破坏模式[11]。

2008年以来，三峡库区已经经历了12个周期的145.00～175.00 m水位变动。三峡库区巫峡段内青石-抱龙一带顺层消落带岩体强烈劣化。本文以青石6号坡为研究对象，利用连续-非连续数值模拟，分析斜坡裂缝扩展、变形破坏机理和破坏运动过程。然后结合室内饱和-风干循环试验结果，讨论了岩体劣化对岸坡变形破坏过程所产生的影响。研究结果可为三峡库区顺层岩质岸坡的防灾减灾提供技术支撑。

1 青石6号坡工程地质概况

青石6号坡位于三峡库区巫峡段右岸，距上游巫山县城16 km。坡体在高程170 m处，地形起伏较大，将坡体分为上下两部分。坡体上部呈马鞍状，凸起端长度达218 m，坡体下部坡面角度逐渐转缓，形成陡缓平台。坡体两侧冲沟为上下游边界，前缘入江，受库水位波动影响(见图1)。水位为145.00 m时，斜坡纵长约410 m，宽约为170 m，总体坡向345°。坡体表层出露三叠系下统嘉陵江组3段(T1j3)白云岩夹泥质灰岩，岩层产状为340°～350°∠40°～50°(见图2)。坡体由多个岩层组成，且岩层倾角略大于坡角，为典型的顺层岩质岸坡。研究区位于楠木园背斜北翼、官渡口向斜南翼，地层稳定，无活动断裂通过。

岸坡消落带范围内岩体发育有两组产状分别为34°∠70°和270°∠80°的优势结构面，平均间距约为2.8 m，仅在坡体表层发育。岩体结构面与层面相互切割坡面岩体，形成碎块状板状六面体岩体。岩层面间以钙质胶结物填充，少量岩层面处于张开状态。在库水长期淘蚀作用下，岩层节理的张开程度逐渐扩大。在上覆岩体的重力挤压作用下，消落带范围岩体的破碎程度不断增加，岩体强度不断下降。

经长期变形演化后，坡体上出现明显的变形破坏迹象。坡体高程165 m处，间续发育有多处波状弯曲现象，从坡体上游边界间断延伸至坡体下游边界。在坡体两侧冲沟内观察，发现波状弯曲由表层向内部延伸有一定深度(见图3)。波状弯曲挠曲段由于弯折角度过大，新生大量纵向裂缝。在对青石岸坡范围内的21个斜坡进行调查后发现，在互不相邻的9个斜坡上，均发育有弯曲程度不同的波状弯曲。青石岸坡上所出现空间分布不连续、弯曲程度大小不一的波状弯曲现象说明，在青石岸坡上出现的波状现象并不是由地质营力所致，而是由坡体上部重力作用所致。

2 青石6号坡变形破坏全过程分析

受自重和库水位波动的影响，青石6号坡坡面岩体较为破碎，质量较差，具有明显变形迹象。不同破坏模式下，坡体运动方式不同，产生的灾害大小也有所差异。本文选用Munjiza提出的FEM/DEM耦合分析方法，对青石6号坡进行数值计算，预测其变形破坏全过程。

2.1 FEM/DEM原理

本文采用文献[12]中的方法，该方法结合了FEM(有限元)和DEM(离散元)优点，能够模拟出层状岩体中裂缝扩展和破碎过程，得到岩质边坡渐进式损伤和破坏过程[12-13]。

FEM/DEM耦合分析的基本思路是将整个坡体划分为三角形有限元网格，网格之间由4个节点无厚度的节理单元相连[12](见图4)，形成完整的坡体结构。通过在有限元网格中插入离散元网格表示滑体区域，允许离散元网格发生转动与位移。模型中的破坏准则包括抗剪强度、抗拉强度和断裂能。

对坡体模拟过程中，受重力作用影响，连接离散元网格的节理单元发生错动时滑体仅沿着滑动面产生蠕动变形；进一步演化时，节理单元发生分离，离散元网格相互张开，滑体开始出现破裂。通过对坡体连续-非连续的分析[13-14]，实现了变形破坏全过程模拟。上述模型的主要特点为既可对坡体初始阶段所发生的小变形进行模拟，也可对坡体在破坏过程中所发生的大规模破坏、运动岩体相互分离进行模拟，从而反映出边坡岩体变形破裂这一复杂岩体力学行为。

在岩质边坡分析方面，FEM/DEM耦合分析方法可以模拟坡体中裂缝扩展和岩体破碎过程，完整展示坡体的整个运动过程，对分析岩质边坡具有独特的优势[15-16]。国内外学者在岩体裂纹扩展[17-18]、洞室开挖[19]、岩质边坡[20]、危岩体分析[21]等方面均取得了较好成果，验证了FEM/DEM方法的可行性和普适性。但由于FEM/DEM后处理程序缺少对网格位移量的追踪，不能得到坡体运动速度和位移量，无法对坡体变形情况定量分析。

2.2 模型建立

青石6号坡数值计算模型长377 m，高297 m，由基岩岩体与滑体构成(见图5)。整个计算模型由2 000个均匀大小的计算网格组成，网格之间由节点相连。利用有限元网格单元模拟基岩岩体，离散元网格单元模拟滑体区域，节理单元模拟层面。从坡体结构来看，青石6号坡坡体中部出现了明显的凸起段，在上部岩体的重力作用下，岩层顺着层面滑动。由于坡体前缘坡面角度逐渐转缓，岩层滑移受阻，使岩层发生波状弯曲。故将坡体凸起段划分为滑体区域，其余区域在模拟过程中不发生变化，将剩余区域划分为基岩岩体。

计算模型中所需的岩体强度参数如表1所列，各部分取值先结合现场实测和经验进行综合取值，再通过对比坡体实际变形破坏情况和数值模拟结果进行不断调整优化。表中基岩岩体和滑体区域主要由灰岩岩体组成，强度较高。滑动面位于岩层层面之间，强度主要由填充其中的泥质填充物维持。经长期库水掏蚀，泥质填充物之间的胶结能力越来越差，松散程度不断增大，致使滑动面的岩土体强度大大降低。所取得滑动面内摩擦角远低于滑动面倾角，使得上覆岩体具备沿滑移面向下滑移的条件。

表1 计算模型中岩体强度参数Tab.1 Rock strength parameters of calculation model

在计算过程中，使用张拉破坏断裂能和剪切破坏断裂能对岩质边坡破坏前和破坏过程中所出现的岩体破裂、损伤和软化现象进行模拟[22]，通过不断计算节点位置得到岩质边坡的运动后状态。对模型底部和左侧边界节点/单元施加固定约束，固定坡体位置。对坡体内部节点/单元施加重力，模拟坡体在重力作用下变形破坏过程。计算模型中的最大不平衡力与初始最大不平衡力比值低于10-5或离散网格的物理位移趋近于一个常数不再变化便视为模型已达到应力平衡，单元网格不再发生变化[23]。

2.3 数值模拟结果分析

本文中的工况仅考虑坡体在重力作用下的变形演化过程，并未涉及其他因素。基于数值模拟结果，可以将其变形演化过程划分为裂缝新生、扩展、弯曲-隆起、裂缝贯通、大规模滑移等5个阶段。初始阶段，坡体结构完整，无网格单元张开(见图6)，以连续分析为主。基岩岩体中节理单元强度较大，在接下来的演化过程中处于稳定状态，计算全过程为连续分析。

裂缝新生阶段：坡体中部岩体(主动传力区)在重力作用下内部应力逐步调整。坡体中下部岩体(被动挤压区)发生蠕滑，滑体发生剪张破坏，剪力向临空面方向作用，形成张拉裂缝(见图7)。此阶段，仅坡体中下部岩体(被动挤压区)的节理单元断裂，少量网格相互分离，滑体以连续分析为主，占坡体变形破坏总时步的33%。

裂缝扩展阶段：在坡体中部岩体重力的持续作用下，坡体中应力不断调整，坡体中上部岩体(主动传力区)继续蠕变，中下部岩体(被动挤压区)中的新生裂缝由内部向两端不断张开扩展，向上延伸长度增加(见图8)。此阶段，被动挤压区中节理单元断裂数量不断增加，使得三角形单元张开度逐渐增加。滑体下部开始转为非连续分析，上部仍保持连续介质。这一阶段占坡体变形破坏总时步的23%。

弯曲-隆起阶段：随着作用时间的增加，被动挤压区垂直于层面的压应力远小于顺层压应力，由于两者压力差过大，致使坡面轻微弯曲隆起。此阶段中，在坡体内部形成明显的应力集中区[24]和“X”形节理。弯曲-隆起区域位于坡面弯曲起伏转折段，使得滑移段岩体下滑受阻，有利于此现象出现。随着岩层弯曲隆起程度不断增大，当弯折角度过大时，弯曲部位开始出现纵向裂缝，岩体逐渐碎裂化发生松动，在弯曲部位前缘出现岩体崩落(见图9)。此阶段，被动挤压区中节理单元的断裂数量持续增加，滑体与基岩岩体发生错动进而发生节理单元变形，滑体以非连续分析为主，占坡体变形破坏总时步的33%。

裂缝贯通阶段：随着坡体中部岩体(主动传力区)蠕滑位移的持续增加，对坡体中下部岩体(被动挤压区)的挤压程度不断增大，弯曲-隆起程度增加。隆起部位内的岩体裂缝相互贯通形成完整剪出面，坡体向临空面已产生一定位移(见图10)。此阶段，滑体中三角形网格的张开度不断扩大，使得相邻张开区域相连，滑体以非连续分析为主，占变形破坏总时步的11%。

滑移破坏阶段：坡体中下部岩体(被动挤压区)阻抗能力达到极限，滑移面贯通，形成完整剪出口。坡体顺着坡面方向发生“滑移-剪断”引起大规模滑移，最终堆积在坡脚(见图11)。此阶段，滑体中的三角形单元相互分离逐渐碎片化，以非连续分析为主。

从以上分析可知，青石6号坡裂缝主要是由于岩层相互挤压形成。滑体对阻滑端岩体形成挤压，岩层发生弯曲-隆起，随着弯曲程度的增加，岩层破碎程度增大，逐渐形成完整的滑移面，沿着剪出口发生大规模滑移。青石6号坡数值分析过程中所出现的滑移-弯曲现象，符合现阶段坡体表面所出现的变形迹象，可以用于预测斜坡的破坏过程。对于青石6号坡这类库水型岸坡，水位变动对消落带岩体所产生的损伤劣化作用，在岸坡的变形演化过程中起到了极为重要的作用[25]。如进一步结合库水位变动对消落带岩体强度和结构所产生的不良影响，讨论不同劣化情况下岸坡变形破坏模式如何变化，所得结果会更贴近库区的实际工况，适用范围更加广泛。

3 库水位长期变动对岩体劣化的影响

3.1 库水对岸坡的损伤劣化作用

自三峡水库开始正式运行后，库水位经周期性调度后，消落带岩体长期处于“饱和-风干”条件，水岩相互作用会造成岩体强度劣化。刘新荣等通过干湿循环试验模拟砂岩在库水位涨落情况下，砂岩抗剪强度劣化规律[26-27]。邓华锋等发现经循环加载作用后，内部损伤的砂岩岩样在浸泡-风干作用下抗剪、抗压强度劣化效应明显，具有明显的时间效应和非均匀性[28]。以上说明库水位周期性升降会使岩体强度逐步劣化，但是对岸坡稳定性和变形破坏过程会产生何种影响，仍然值得探讨。

本文选取青石-抱龙段4个不同区域的灰岩(T1j3)岩样，按照试验规程要求进行“浸泡-风干”循环试验[29]。其中试样1，2位于岸坡消落带之上，试样3，4位于消落带范围内，已经历过多次水位变动，由此试样3，4的初始单轴抗压强度低于试样1，2的初始单轴抗压强度。每循环5个周期对岩样进行力学强度测定，测定经过0，5，10，15，20，25，30，35，40，45，50次循环试验后灰岩(T1j3)岩样的饱和单轴抗压强度，得到在不同循环周期时灰岩单轴抗压强度下降曲线(见图12)。经过50次循环试验后，4个试样饱和单轴抗压强度的总下降率约为18%～26%，单次循环平均下降率约为0.37%～0.51%(见图13)，总体相差不大。

受尺寸效应的影响，室内岩样的下降速率并不能表征岸坡岩体下降速率。闫国强等根据灰岩的劣化损伤特征，引入核衰变劣化控制方程[30]：

N(t)=N0e-λt

(1)

式中：N0为初始物理强度参数，λ为劣化常数，t为循环周期。

式(1)可以较好地描述灰岩岩体物理力学参数的衰减过程。本文利用核衰变劣化控制方程对上述结果进行修正，得到在不同循环次数下灰岩的饱和单轴抗压强度，如表2所列。

由表2可知，4个试样饱和单轴抗压强度的总下降率近似相等，在50次循环周期后约为29.88%。此下降率明显大于室内干湿循环作用下的下降率结果，考虑到了岩体的实际物理力学属性，更加贴近实际结果。

表2 采用劣化控制方程得到的灰岩饱和单轴抗压强度Tab.2 Saturated uniaxial compressive strength of limestone by deterioration control equation

由此将消落带岩体(坡面高程145～175 m内)单个循环周期强度下降率近似取值为采用劣化控制方程下灰岩饱和单轴抗压强度平均下降率，约为0.6%。在数值模拟计算时，对消落带范围岩体黏聚力、内摩擦角与滑动面黏聚力以0.6%的下降率折减。统计在不同循环周期下，岩体强度下降率逐渐增加时，坡体失稳所需时步(坡体内部裂缝相互贯通时坡体失稳)和变形破坏模式是否发生变化。

经计算发现，随着消落带范围岩体强度下降率增加，坡体失稳所需时步逐渐减少，变形破坏模式并不会发生变化(见表3)。消落带岩体强度降低，坡体阻滑段岩体阻抗能力下降，是造成坡体失稳速度加快的主要原因，坡体出现弯曲-隆起的时间也逐步提前。消落带岩体强度下降程度的高低，仅会影响岸坡变形演化速度，并不会改变坡体变形破坏模式，仍以滑移-弯曲为主。此劣化方程未考虑到岩体结构面对岩体物理力学强度的影响，会导致所取下降率依然偏小，坡体失稳所需时步增大。如在后续的研究中考虑在水位波动条件下，岸坡岩体结构面劣化效应对岩体强度下降率的影响，会更加贴近实际工况。

表3 不同循环周期下坡体失稳所需时步变化Tab.3 Time step of slope instability under different cycle periods

3.2 岩体劣化对岸坡变形破坏模式的影响

岩体强度下降仅是消落带岩体劣化的表现形式之一，更为常见的是对岸坡结构的破坏。在岩层面发育的顺层岸坡上，经长年累月的库水掏蚀，逐步发育有大量沿层面与垂直于层面的节理裂隙切割岩体，使得坡面逐渐碎块化，并被侵蚀掏空形成凹腔。将消落带范围岩体掏蚀深度定为0.97，2.68，4.49，6.74 m，并作为数值模拟的4种工况(见图14)。得到不同岩体结构劣化结果下，对坡体变形失稳所产生的影响。

当掏蚀深度达0.97 m时，在消落带区域形成凹腔，使得坡体阻滑段岩体的阻抗能力急剧下降，失稳所需时步大幅度缩减。由于掏蚀区域较少，并未影响坡体变形失稳模式。当掏蚀深度达2.68 m时，形成的凹腔区域扩大，阻滑段岩体变薄。滑体前缘临空，顺着层面直接滑移，破坏模式由滑移-弯曲转变为滑移(见图15)。当淘蚀深度分别为4.49，6.74 m时，凹腔区域不断扩大，阻滑段岩体逐步消失，滑体前缘临空程度增大，坡体更易直接滑移，失稳所需时步进一步缩减(见表4)。

表4 不同劣化程度下坡体失稳模式变化情况Tab.4 Change of slope instability mode under different deterioration degree

对青石6号坡而言，由于其阻滑段岩体处于三峡库区消落带范围内，长期库水位波动的影响会对阻滑段岩体强度和结构产生劣化效应，影响坡体稳定性。从上述分析可知，消落带岩体强度的劣化并不会改变坡体的变形失稳模式；而坡体结构对劣化程度的响应达到临界点时，坡体变形破坏模式会发生变化。相比岸坡其他区域，消落带区域岩体所处应力环境更复杂，劣化效应极为明显，影响到未来库岸边坡变形演化方式。黄波林等也发现：“水-岩-岸”互馈作用，加快了层面贯通；同时由于趋表-后退效应，不断下切的消落带岩体，改变了消落带岩体结构，两者共同作用加快了岸坡变化速度[25]。王健等也发现坡脚岩体在不同下切深度下，岩质顺层岸坡的变形破坏模式会发生相互转换[31]，印证了本文通过降低消落带岩体强度和破坏岩体结构，模拟在长期库水波动条件下岸坡变形破坏过程的可行性。

4 结论与建议

以FEM/DEM耦合分析方法作为青石6号坡数值模拟手段，结合室内饱和-风干试验结果，对坡体变形破坏机理、库水对岸坡损伤劣化作用与变形破坏模式的影响等方面进行研究分析，得出以下结论：

(1) 以FEM/DEM耦合分析方法预测得到了青石6号坡变形破坏全过程，将其划分为裂缝萌生、扩展、弯曲-隆起、裂缝贯通、大规模滑移等5个阶段，坡体变形破坏模式为弯曲-隆起。模拟结果中出现的弯曲-隆起现象与坡体目前出现的变形迹象相对应，证明了模拟结果的准确性和可靠性。

(2) 根据室内“饱和-风干”试验，得到灰岩(T1j3)岩样饱和单轴抗压强度下降速率，以此对消落带岩体进行强度折减。计算表明，岩体强度的下降仅会加快岸坡变形演化速度，并不会影响坡体变形破坏模式。

(3) 随着坡体结构劣化程度的增加，致使阻滑段岩体逐渐变薄，阻滑能力降低，滑体前缘临空程度不断增大，使得青石6号坡变形破坏模式由滑移-剪断向顺层滑移转变。

目前青石6号坡正处于蠕变阶段，但是随着库水位对消落带岩体的长期劣化作用，会使得岸坡岩体的破碎程度逐渐增加，加快坡体的变形演化速度。整个青石-抱龙段存在类似的20个顺层岩质斜坡，应当对其进行长期的监测。针对各个斜坡目前所处阶段，进行相应的监测或治理。否则一旦发生失稳破坏，滑体入江会带来严重的涌浪灾害，对长江航道与沿岸居民带来严重的威胁。