樟木堆积体斜坡动力稳定性与极限承载力评价

王海芝，王颂，周剑，张路青，*，黄福有

(1.北京市地质研究所，北京 100120；2.中国科学院地质与地球物理研究所中国科学院页岩气与地质工程重点实验室，北京 100029；3.中国科学院大学地球与行星科学学院，北京 100049；4.中国科学院地球科学研究院，北京 100029；5.北京工业大学 城市与工程安全减灾教育部重点实验室，北京 100124)

樟木口岸位于西藏聂拉木县的樟木镇，是中国通向南亚次大陆最大的国际贸易口岸，担负着中国和尼泊尔之间的政治、经济和文化交流。樟木镇地处喜马拉雅造山带核部的高山峡谷，是西藏地质灾害最严重的区域之一，历史上曾多次遭受滑坡、崩塌、泥石流和冰湖溃决灾害(李云贵等，1995)。由于受地形条件的限制和中尼通商的需要，樟木镇和口岸的大多数房屋建筑均修建于樟木滑坡堆积体上。樟木滑坡为一个特大复合型滑坡，坡体上不仅发育有福利院和帮村东古滑坡，而且在古滑坡堆积层上发育有现代滑坡，整个滑坡体积约1亿m3。运用自然历史分析法进行推测，樟木滑坡堆积体的形成过程大致可以分为3个发展阶段：岩质古滑坡的形成及剥蚀阶段、崩坡积物及冲洪积物在岩质古滑坡上的堆积阶段、堆积层上新老滑坡的发育阶段(易顺民等，1996；毛成文，2008)。

近年来,随着城镇化建设的推进，地质灾害频发，樟木滑坡局部出现复活，造成众多建筑物变形毁坏，道路裂缝发育，地表沉陷，严重威胁城镇的进一步发展(曾建，2006)。1986～2004年,原公安局、福利院等多处地面出现地裂缝(祝建等，2008)。2005年雨季期间，福利院滑坡变形加剧，导致管委会宿舍楼地面下沉超过30 cm、烈士陵园挡墙整体发生前移和纪念碑底座严重扭曲变形。2010年1月2日～2011年12月14日，西藏自治区地质环境监测总站樟木滑坡监测点的监测数据表明，樟木滑坡5个监测点的累计位移介于15～112 mm，局部变形强烈(尚文涛，2009)。2015年4月25日尼泊尔特大地震进一步加剧了樟木滑坡的变形(陈剑等，2016)。滑坡灾害不仅威胁到当地人民群众的生命和财产安全，也制约着口岸经济社会的发展。

祝建等(2010)对樟木滑坡地表裂缝、地表位移和深部位移进行了长期监测，结果表明,福利院古滑坡中部发生复活。胡瑞林等(2014)综合工程地质结构和岩土类型，初步研究了樟木滑坡堆积体的变形破坏规律和稳定性，认为斜坡整体稳定，但受人类工程活动等因素的影响局部可能会发生滑移。Ma等(2017)基于樟木滑坡的100个钻孔记录、电子自旋共振测年法和年轮宽度测年法，确定了不同部位的复活时间，并结合现场监测结果，发现滑坡活化与降雨密切相关。上述研究指明了樟木滑坡的影响因素及其与斜坡的相互作用关系，但还缺乏对滑坡稳定性的定量评价。

数值模拟是研究滑坡的重要手段，有利于揭示滑坡的形成机理和变形破坏过程。例如，有限元方法可以处理各向异性的复杂边坡工程问题，计算边坡的塑性区及获取位移场、应力场和应变场，并能够进行承载力分析(Lollino et al.，2010；吴应祥等，2013；Gian et al.，2016；陈豫津等，2019)；有限差分法可以考虑边坡的非线性应力和应变关系，并能够基于强度折减法对边坡进行大变形分析(Wang et al.，2020；Zhang et al.，2020；刘畅等，2020)；离散元法以牛顿第二定律和单元接触关系为基础，能够再现边坡破坏后的运动学特征(Wei et al.，2019；Zhou et al.，2019；周礼等，2019)。尚文涛(2009)在FLAC 3D软件中利用强度折减法研究了樟木滑坡的变形和破坏，认为坡体以“压-剪”破坏模式为主；左辉(2009)使用GeoStudio软件分析了帮村东滑坡在天然和饱水状态下坡体的稳定性，结果表明，在饱水时坡体处于欠稳定状态。张俊文等(2016)使用有限元强度折减法计算了福利院滑坡在自然、降雨和地震工况下的安全系数，发现强降雨和地震使得坡体的安全系数显著降低。

上述有关樟木滑坡稳定性的研究大多都是针对暴雨工况和采用强度折减简化的地震工况，忽略了地震的动态作用过程和人类活动对坡体的影响。笔者定量分析了地震动和车辆震动等动力变化对堆积层斜坡稳定性的影响，并计算了随动载荷作用而变化的动态安全系数。在斜坡稳定性分析的基础上，还采用普遍条分法进一步评价了樟木镇堆积层斜坡的极限承载力，以便为樟木镇建设用地的优化布局和可持续发展提供参考。

1 樟木镇滑坡概况及特征分析

1.1 工程地质条件

樟木镇隶属西藏日喀则地区聂拉木县，地处喜马拉雅山南侧的高山峡谷区。滑坡区内分布有自东向西流动的波曲河，河岸两侧地势陡峭，相对高差超过2 000 m(图1)。樟木滑坡主要发育于波曲河左岸，坡体平均坡度超过35°，其中海拔3 500 m以上区域的坡度甚至大于50°。勘探资料显示，樟木滑坡区地层岩性简单，基岩为前寒武纪曲乡岩组片岩等深变质岩，表层为崩坡积物、冲洪积物和部分人工填土。沿滑坡体不同高程处还分布有厚度小于2 m的滑带土。整个滑坡区域位于喜马拉雅地槽褶皱区内，活动断裂极其发育，构造运动异常强烈。受此影响，斜坡岩体十分破碎，风化剥蚀速度加快，仅新构造运动以来沟道下切速率超过10 mm/a。

图1 樟木滑坡堆积体位置示意图(据胡瑞林等，2014)Fig.1 Location map of Zhangmu landslide (Modified from Hu et al.,2014)

1.2 滑坡堆积体特征分析

樟木镇滑坡由福利院堆积层古滑坡与帮村东堆积层古滑坡组成(图2)，受降雨、地震和人类工程活动的影响，在古滑坡上又发育有多个具备一定活动性的现代滑坡。例如，福利院古滑体上的消防队变形区、烈士陵园变形区、中心小学变形体和帮村东现代滑坡。

图2 樟木镇福利院堆积层古滑坡与帮村东堆积层古滑坡图(据张俊文等，2016)Fig.2 Location of Fuliyuan and Bangcundong landslides (Modified from Zhang et al.,2016)

福利院滑坡位于樟木镇西北部，318国道以“Z”字形穿越滑坡体中下部，坡体陡缓相间，呈明显的“圈椅状”形态。整个地形坡度范围为25°～35°，滑坡后缘坡度可达45°，滑坡体中上部还发育有一个坡度小于15°的缓平台。现场调查发现，福利院滑坡前缘、中部和后缘均出现明显的变形迹象，整个滑坡区内大量房屋倾倒变形且地面裂缝密布。

帮村东滑坡位于樟木镇东南部，覆盖面积约占樟木镇面积的50%。滑坡体形态与福利院滑坡的类似，坡体亦是陡缓相间。地形上东南高，西北低，整体坡度范围为30°～35°，滑坡后缘坡度超过40°。该处滑坡活动以地表变形为主，滑坡区内地面和建筑物均出现裂缝，斜坡上遍布“马刀树”。

2 堆积体稳定性分析

樟木滑坡是在降雨、地震和人类工程活动等复杂因素作用下形成的复合型滑坡。考虑到堆积体稳定性的影响因素复杂，本节着重分析在地震动力和车辆动载荷作用下斜坡的稳定情况。

2.1 地震动力稳定性分析

常用的斜坡地震稳定性分析方法有拟静力法和Newmark滑块分析法。拟静力法是将地震荷载简化为静荷载进行计算，斜坡不同高程滑动块体的地震峰值加速度均为设计拟定值，每一滑块的稳定性计算均考虑水平指向坡外惯性力(由滑块质量乘峰值加速度拟定值来确定)。但拟静力法只能计算获得安全系数，而不考虑滑动面的变形，Newmark滑块分析法则可以克服上述缺陷，因此在斜坡地震安全性评价中得到了较为广泛的应用。在Newmark滑块分析法中，首先需要假定滑动面并根据极限平衡原理确定坡体的屈服加速度值，然后通过动力分析判别坡体是否产生滑移，并估算永久滑移位移。笔者假定当永久滑移位移超过5 cm时，斜坡失稳。采用这种永久位移标准能够避免重复的数值计算。

Newmark滑块分析法计算所需的地震动参数根据《中国地震动参数区划图》(中华人民共和国国家质量监督检验检疫总局等，2015)确定，樟木的地震动峰值加速度为0.20 g，地震动反应谱特征周期为0.45 s，区域地震基本烈度为Ⅷ度，由此合成的地震动反应谱和地震加速度曲线见图3和图4。计算过程中所需的岩土体参数见表1。

表1 滑坡体各岩土层参数表(据张俊文等，2016)Tab.1 Material parameters of rock and soil mass (After Zhang et al.,2016)

图3 地震动反应谱图Fig.3 Ground motion response spectrum

图4 人工合成的地震动加速度曲线图Fig.4 Synthetic seismic acceleration curve

以樟木镇中心小学强变形体为例，分别计算了天然条件和暴雨条件下斜坡在不同时刻随地震动力作用而变化的安全系数，即动态安全系数(表2)。

表2 中心小学变形体在天然条件和暴雨条件下的计算模型及动态安全系数表Tab.2 Calculation models and dynamic safety coefficients of deformation body in Central Primary School under natural and rainstorm conditions

在地震动力作用下，中心小学强变形体安全系数也呈现动态变化，天然条件下该动态系数基本大于1，地震位移约为0，但在暴雨条件下，动态安全系数大幅减小，地震位移随时间不断增大。

所有计算结果见表3。整个樟木滑坡中，福利院堆积层古滑坡在天然条件下滑移距离为6 cm，说明在地震作用下会发生失稳；中心小学变形体在暴雨条件下的地震滑移距离达1.6 m，也将发生整体失稳。

表3 樟木滑坡动态安全系数表Tab.3 Dynamic safety coefficients of Zhangmu landslide

2.2 车辆动荷载作用下的斜坡稳定性分析

樟木镇是中国与尼泊尔的重要通商口岸，公路运输是该地区贸易的唯一运输方式。中国境内的318国道建在樟木斜坡体上，行进中的车辆对斜坡稳定性必然会有一定影响，尤其是在雨季期间，斜坡岩土体的力学性质将会明显弱化，载满物品的车辆完全有可能导致路基坡体产生滑塌。为了定量认识车辆震动对斜坡稳定的不利影响，笔者将对相应的动力学问题进行初步分析。

不同载重的车辆产生的动载荷大小不同，载重量为100 kN、150 kN、200 kN、300 kN和550 kN的相关指标见表4和图5。

表4 汽车的主要技术指标(据高玉峰等，2005)Tab.4 Main technical indexes of vehicles (After Gao et al.,2005)

图5 不同载重车辆的平面尺寸图(据高玉峰等，2005)Fig.5 Plane dimensions of different vehicles (After Gao et al.,2005)

比较不同载重车辆车轮与路面的接触力，可以看出550 kN车轮的单个后轮与地面的接触压力为最大(70 kN)，因此其对高速公路软基沉降的影响更为显著。将车轮处的荷载处理为作用在公路横截面上的均布荷载，由此得到550 kN车辆后轮作用在公路横截面上的等效压应力为583.33 kPa。550 kN车辆缓慢经过公路某一截面的动荷载可以简化为半正弦动荷载(图6)。

图6 简化为半正弦的车辆动荷载图Fig.6 Diagram of vehicle dynamic load

如图7所示，318国道从福利院现代滑坡中消防队强变形区前缘经过，以该剖面为例分析重型车辆对斜坡稳定性的影响。

图7 车辆荷载作用下坡体稳定性的计算模型图Fig.7 Calculation model of slope stability under vehicle dynamic load

对图7中公路所在位置施加了图6所示的动态荷载，潜在滑体的动态安全系数、整体滑移速度及滑移距离分别见图8、图9a和图9b。计算结果表明，潜在滑体初始安全系数为1.044，550 kN的车辆缓慢行驶在国道上时，安全系数在0.974～1.065间变化。动态安全系数小于1.0时，潜在滑体产生下滑，下滑速度最大可达0.023 m/s，而最终滑动位移约为6 mm。尽管并未造成潜在滑体的失稳，但上述计算结果足以说明重型车辆在潜在滑体上部行驶时对坡体的稳定十分不利，长期作用下导致欠稳定斜坡最终失稳的可能性较大。

图8 车辆荷载作用下的斜坡动态安全系数图Fig.8 Dynamic stability coefficients of slope under vehicle dynamic load

图9 车辆荷载作用下潜在滑体的整体滑移速度和滑移距离图Fig.9 Slip velocities and distances of the sliding mass under vehicle dynamic load

3 堆积层斜坡的极限承载力分析

樟木镇整体处于地质灾害危险区内，而土地利用优化的目的是在其中找到相对安全的区块(段)，以实现口岸的持续发展。福利院和帮村东堆积层古滑坡体在暴雨条件下仍有一定的加载空间，应当分析其上部所能承受的极限荷载(斜坡极限承载力)。斜坡极限承载力是指斜坡保持稳定所能承受的极限外部荷载，可以通过最危险滑面(此时的滑体稳定性最低)下滑体的极限平衡条件来确定。这里的外部荷载是指建筑荷载、车辆荷载、挖填堆载等产生斜坡表观的静荷载。在外部荷载作用下，斜坡内部的应力状态发生改变，当某一点的剪应力达到抗剪强度时，这一点就处于极限平衡状态。随着外部荷载的继续增大，坡体内处于极限平衡状态的范围也加大，并最终发展成为贯通的滑面。

3.1 极限承载力的计算分区及分析方法

在斜坡极限承载力的实际计算过程中，采用增量法逐级增加外荷载，直到所加荷载沿着该滑动面的整体安全系数趋近于1。对于安全系数的计算采用普遍条分法，其满足2个静力学平衡条件：水平方向与铅直方向的力平衡；所有方向的力围绕一点的力矩之和为0。

依据地形条件、地质特征及已有滑坡分布状况，可以将樟木斜坡堆积体划分为23个区块；由于8区与9区都包含着强变形区，依据变形程度分别划分出4个亚区，即亚区8-1、8-2、9-1和9-2(图10)。

图10 极限承载力计算区块划分及计算结果图Fig.10 Block division and calculation results of ultimate bearing capacity for Zhangmu landslide

其分析方法如下：

(1)针对既定的计算区块，搜索斜坡的优势滑面及最危险滑面，确定潜在的最危险滑体，计算相应的安全系数(图11a)。

(2)计算潜在滑体各条块底部倾角与内摩擦角的相对大小，确定地表等效均布荷载的作用位置(图11b)。

(3)逐渐增加等效均布荷载，直至潜在滑体的安全系数为1±0.01，对应的等效均布荷载即为斜坡极限荷载(图11c)。

图11 极限承载力分析方法图Fig.11 Diagram of analysis method of ultimate bearing capacity

3.2 堆积层参数选取

鉴于樟木镇所在区域降雨丰沛，雨季较长，相应的极限承载力分析将采用饱水状态下的岩土体强度参数(表1)。

3.3 极限承载力分区计算结果

依据计算区块划分(图10)及樟木镇的工程地质平面图(图2)，选择已有的勘查剖面分别计算1～23号区块中的极限承载力。为保证计算结果的代表性，所确定的剖面包括岩质滑坡剖面、福利院堆积层古滑坡剖面和帮村东堆积层古滑坡剖面。

福利院堆积层古滑坡上部发育的消防队强变形区(区块9-1)、烈士陵园强变形区(区块8-1)和中心小学变形体(区块11)稳定性较差，地表不宜继续加载，直接将相应的极限承载力取为0。上述区块以外的其他区块还有不同程度的承载空间，相应的极限承载力计算是建设用地适宜性评价的主要依据之一。按照前述承载力计算方法，计算各区块在等效均布荷载作用下的斜坡极限承载力，部分计算模型见图12，统计各区块结果并标注于图10中。

图12 极限承载力计算模型示例图Fig.12 Examples of ultimate bearing capacity calculation models

计算结果表明，樟木堆积层斜坡极限承载力相对较高的部位主要有3个区段：福利院堆积层古滑坡与帮村东堆积层古滑坡分界线附近的山梁(其中，区块12～15的极限承载力分别为230 kPa、160 kPa、270 kPa和620 kPa)；帮村东堆积层古滑坡的东部(其中，区块22和23的极限承载力分别为310 kPa和300 kPa)；樟木镇堆积体后缘(其中，区块3和4的极限承载力分别为330 kPa和500 kPa)。堆积体古滑坡的前缘(区块7、10和16)和后缘(亚区9-2、区块19、20)的承载力偏低(低于100 kPa)，其他区块(亚区8-2、区块1、2、5、6、21)的极限承载力处于100～200 kPa。

4 结论

综合宏观调查研究和模拟计算可知，樟木滑坡堆积体的稳定性受区域地形、地质、降水、地表径流、地震及人类工程活动等因素的控制。地表水与地下水入渗、工程活动与地震是滑坡的激发因素。相关结论如下：

(1)地震动力稳定性分析表明，仅地震作用下福利院堆积层古滑坡局部失稳；在极端条件下(地震+暴雨)，动态安全系数进一步减小，中心小学变形体将发生整体失稳。

(2)车辆动荷载作用下的斜坡稳定性分析表明，重型车辆缓慢行驶在潜在滑体上部时安全系数在0.974～1.065之间变化，潜在滑体下滑速度最大可达到0.023 m/s，最终滑动位移约为6 mm。

(3)依据地形、地貌、地层岩性及新老滑坡边界，将樟木堆积层划分为23个区块，利用极限平衡法计算了饱水条件和等效均布荷载作用下不同斜坡区块的极限承载力。现代滑坡及强变形区接近极限平衡状态，没有进一步的承载空间，应禁止增加表观荷载；堆积层古滑坡的前缘和后缘区块的承载力偏低，其他区块还具有一定极限承载力。

相关计算结果可为樟木滑坡灾害的治理和樟木镇建设用地的优化布局提供科学依据。