李俊峰,陈红旗,刘红岩
(1. 中国地质环境监测院,北京 100081;2. 中国地质大学(北京)工程技术学院,北京 100083;3.西藏大学工学院,西藏 拉萨 850000)
2011年8月26日,位于峨眉山市九里镇兴阳村九沙河右岸“王山-抓口寺”不稳定斜坡因连续降雨发生滑动破坏,滑坡体总体积约5.0×106m3[1]。近年来该斜坡一直处于不稳定状态,2015年6月14日,该滑坡再次启动,27户共52人安全受到威胁,滑体前缘冲进九沙河,形成厚度达50 m左右的堰塞坝体,汛期对下游存在重大安全隐患[2]。该滑坡前后两次避险成功得益于排查及时、预案完备、应急响应和处置及时,灾害隐患点的早期发现和跟踪监测发挥重要作用。应急现场实际经验说明,滑坡运动机制分析,以及综合利用数值模拟手段验证概念模型的准确性和灾害发生过程的反演和预测有利于科学防灾减灾行动,有助于减少人民生命财产的损失。
滑坡工程地质问题强调以地质力学的观点、自然历史的观点从灾害成因分析、发展过程推演和力学机制分析定性认识灾害体现状和掌握或预测变形体发展趋势,最终达到对不稳定斜坡的定量评价[3]。长期以来,公路边坡问题、库岸边坡问题以及工民建人工切坡问题引发崩滑流突发地质灾害接连不断,利用数值模拟技术进行稳定性分析[4-8]、参数敏感性分析[9-11]和滑动过程推演[12-15]是滑坡问题研究的重要方向。滑坡稳定性分析和滑动过程模拟依然是滑坡灾害反演和预测的重难点问题,其中依然存在大部分数值方法并不能将真实时间和数米级以上的大变形分析统一的不足。
非连续大变形分析(Discontinuous Deformation Analysis, DDA)是由美籍华人石根华博士提出的一种非连续介质数值方法,其对力学现象的数学和数值描述与块体真实运动相一致,模拟过程非常接近实际,且兼有真实时间变量和大变形问题处理的优势,是一种可靠的纯动力学方法,满足研究滑坡稳定性和滑坡全过程特征的要求[16]。本文首先通过野外现场调查、非接触式量测等手段对王山-抓口寺滑坡形成机制进行了详细的分析;然后采用DDA方法对滑体进行稳定性及滑坡滑动过程进行详细的模拟,揭示了滑坡滑动过程及机理,加深了对灾害发生原因和过程的理解,有助于防灾减灾的科学决策。
滑坡区域内地下水按含水层性质及埋藏条件,分为第四系松散地层中的孔隙潜水和基岩裂隙水两大类型:
(1)松散地层孔隙潜水:主要赋存于斜坡崩残坡积层孔隙中,为孔隙潜水,其结构松散,地层的透水性好,不利于地下水的富集,地下水贫乏。地下水主要受大气降水的补给,随季节变化幅度较大,部分入渗下部基岩,沿斜坡地形低洼处排泄,汇入冲沟。
(2)基岩裂隙水:主要赋存于二叠系上统峨眉山玄武岩组裂隙中,地下水不易赋存。根据调查玄武岩柱状裂隙发育,节理裂隙连通性好,在陡崖和斜坡山脊地带,上覆土层较薄,基岩部分段直接出露,大气降水和地表水通过近似垂直的裂隙直接入渗,沿贯通裂隙径流或储藏于裂隙中,部分以泉的形式从结构面或在低洼地段溢出,导致部分冲沟段断流。
图1 王山-抓口寺滑坡工程地质平面图[2]Fig. 1 Engineering geological drawing of Wangshan-zhuakousi Landslide
王山-抓口寺滑坡位于九沙河右岸(图1),属构造剥蚀中低山斜坡地貌,山体总呈东西向展布,地势南高北低。滑坡整体呈“箕”形(图2),具有顺向坡地质结构,坡纵向前缓后陡,坡度约26°,坡向330°,滑坡竖向高程区间在530~815.6 m,相对高差285.6 m,中部最宽达445 m,前缘最窄为245 m,滑体最大纵长710 m,面积30×104m2,滑体厚度10~40 m,平均厚度20 m,滑体体积6.0×106m3左右[17]。滑体二次滑动后,滑体前缘堆积厚度约40~50 m的堆积体,堵塞九沙河,形成堰塞湖,汛期对下游存在重大安全隐患。根据现场调查滑坡后壁光面长度及其他标志物运动情况,推测滑坡整体滑距在40 m左右。
该滑坡2011~2015年已发生两次滑动,滑坡后缘及左右侧边界变形明显。滑坡形成较典型微地貌有:滑坡壁、滑坡台阶、左右侧拉裂坎、滑坡舌(逆冲反翘),滑体表面出现“马刀树”。滑坡后缘呈“圈椅”形(图1和图2),滑坡后壁本次未发现明显变形(图3),但后缘平台加宽;左侧以冲沟为界,出现拉裂坎,本次滑动变形加剧,部分位置玄武岩基岩出露,冲沟现已被滑体超覆;右侧边界上部形成高陡坎,本次滑动变形加剧,陡坎高度增加,边界有羽状裂缝,滑体滑舌部分向下滑动遇北侧山体堆积形成反翘鼓丘(图4),至原九沙河河床高度50 m左右。
图2 滑坡滑动后全貌[2]Fig.2 Spatial morphology of landslide
图3 滑坡后壁Fig.3 Main scrap of landslide
图4 滑体前部反向堆积Fig.4 Reverse accumulation of landslide toe
1.3.1形成条件分析
滑坡形成与发展受滑坡区地层岩性及其结构、构造特征、内外营力特征、气象水文特征等多方面影响。
首先,该区出露的峨眉山玄武岩,中间含凝灰岩夹层,属软硬相间的地层。凝灰岩夹层呈薄层状产出,抗风化能力弱,遇水易软化,遭遇地下水长期浸泡极易产生夹层泥化现象,造成层间带物理力学软化严重,是控制边坡变形及破坏方式的内因。
其次,滑坡为顺向坡,层面与坡面以小角度相交(主滑方向近SN向)。区内构造简单,未见其它大的断裂裂隙;岩体节理裂隙发育,主要发育一组陡倾X节理(见上文),长度数十厘米到数十米不等,面多平直光滑。因此,层面与节理裂隙的不良组合是控制滑动的关键,提供了边坡变形失稳边界。
再者,该区受新构造活动影响较强,特别是地震作用的影响。挽近以来本区地壳强烈抬升,两岸坡体急剧抬升,岩体中地应力释放,斜坡岩体向临空方向回弹膨胀,引起应力重分布和应力集中,岩体中原生及次生结构面产生横张拉裂。岩体完整性遭到较大程度的破坏,为地表水及地下水的入渗和运移提供了良好的通道,一方面,水的楔劈力促进了结构面的扩张破坏,降低岩体强度,地下水的静水压力不仅使滑面上有效法向应力降低,从而降低了滑面抗滑力,而且切割面中的静水压力还增加了滑坡体的下滑力,边坡稳定性降低;另一方面也为凝灰岩的遇水崩解提供了可能。研究区受“5·12”汶川地震活动影响强烈,使得岩体产生松动破裂,结构遭到进一步的破坏。
此外,在自重应力场的长期作用下,中上部浅表层岩体发生沿软弱层(带)的蠕滑变形,使得裂隙(包括陡倾节理)进一步拉开拉长,节理裂隙连通性增加,风化作用加剧,岩体结构面的强度进一步低,挤压引起下部岩层向临空方向发生隆胀变形;九沙河河水的淘蚀作用使得下部岸坡的岩体结构遭到严重破坏,坡脚部位差异卸荷变形产生近水平向的剪裂隙,不仅破坏岩体的完整性,也为边坡的最终失稳滑落提供了剪出口。
最后,滑坡调查区降雨丰沛,且强降雨集中,滑坡岩体物理、化学风化作用比较强烈,对岩体的结构和强度均不利;河水位涨落变幅较大,水位反复变动带内岩体强度低、边坡稳定性差。
1.3.2滑坡发生机制分析
工程地质定性分析要求对地质体进行地质过程的机制分析,是定量评价的基础。通过现场调查,初步判断滑坡体演化过程(图5)如下:
(1)斜坡演化的平衡阶段(图5(a)):斜坡演化初期,河谷下切,经过一定的地质构造作用,此阶段斜坡保持动态平衡。
(2)滑坡体蠕滑拉裂阶段(图5(b)):斜坡陡倾裂隙与缓倾角结构面组合(图6)在坡体重力作用下发生时效变形。陡裂与缓裂的组合方式,在地形地貌上就表现为陡坎和阶地。
(3)滑坡体加速变形阶段(图5(c)):在“5·12”强震作用下,岩土体破碎,有利于地表水入渗滑坡体发生滑移拉裂,后缘拉裂缝增宽,在降雨作用下滑坡局部失稳。
(4)滑坡体快速滑动阶段(图5(d)):汛期持续降雨,大量降水经表层堆积物,沿玄武岩裂隙通道入渗至凝灰岩上层面,凝灰岩夹层相对隔水,水沿凝灰岩上界面流动,且凝灰岩夹层中含蒙脱石、伊利石、高岭土等水敏性矿物遇水软化,导致凝灰岩夹层强度急剧降低。坡脚河道流量增大,河水的浸泡、侧蚀、冲刷作用强烈,降低了斜坡体稳定性。双重因素的作用下,滑坡体发生高速滑动解体。
图5 滑坡启动至滑动解体过程演化示意图Fig.5 Schematic drawing of the process from starting to sliding of rock slope
综上,王山-抓口寺滑坡最终形成以后部推移为主、前部牵引为辅的“强降雨-入渗-岩体软化-推拉”破坏模式,属典型滑移-拉裂大型顺层岩质滑坡。
图6 滑坡体陡倾裂隙和缓倾角结构面组合Fig.6 Combination of landslide fractures and gentle dip angle plane
根据上述对王山-抓口寺滑坡周边地质条件、滑坡体结构特征和变形破坏机制的分析,滑体表层为碎石土,滑体中下部玄武岩岩体被近于正交的两组节理切割破碎,滑体实际运动过程中,离散作用明显,为模拟实际块体运动中平移、翻转和碰撞等运动过程,结合非连续变形数值方法大变形和真实时间模拟的优势,选用DDA方法对滑坡启动和滑动全过程进行应急数值模拟研究。由于DDA的计算原理在相关文献[5,7,12-13,16]中已有详细的阐述,因此这里就不再赘述。
2.1.1计算模型概化
根据现场踏勘、滑动面现场调查、基岩出露情况和已有资料综合分析,结合现场三维激光扫描非接触式量测技术[17],获得王山-抓口寺滑坡地质剖面图,并根据现场估计堆积方量和基岩出露情况推测原始坡面位置(图7)。根据滑坡结构和变形特征,滑坡基岩出露清晰,滑床上部基岩面即为滑动面,由现场地质调查可知滑坡滑动原因主要为集中降雨入渗对滑体结构面强度参数弱化,故综合上述情况将滑坡数值模型概化为滑体和滑床两部分:滑床为基岩且结构较为完整,划分为一个块体;滑体部分实际主要包括松散碎石土堆积层、被两组裂隙切割破碎玄武岩岩体和厚度在0~20 cm的凝灰岩层滑带。根据现场实际滑坡滑动和变形情况,综合研究目的对滑坡地质概化模型进行适当简化,将滑体组成的碎石土、切割破碎玄武岩和凝灰岩薄层滑带简化为岩体和软弱结构面两部分,而后在ANSYS中进行网格剖分,导出节点和单元文件,利用自编接口程序将其转为DDA程序能识别的数据结构,生成DDA数值模型(图8)。模型中共划分单元303个,节点223个,该滑坡为大型顺层岩质滑坡,为充分反映滑体力学参数和载荷条件等非结构因素对滑坡的影响,块体单元尺寸应尽量均一[7],并在滑坡体前部、中部和后部各选取监测单元,单元编号分别为40#、11#和5#,具体见图8。
图7 滑坡地质剖面图[18]Fig.7 Geological profile of landslide
图8 滑坡数值模型Fig.8 Numerical model of landslide
2.1.2参数及边界条件
滑坡物理力学参数的取值综合考虑两方面因素:(1)现场踏勘,经验取值;(2)数值参数反演分析。物理力学参数类型的选取和取值直接关系数值计算的有效性和计算的效率。首先采用正交试验设计理念[17],选取最大位移比、结构面黏聚力和结构面内摩擦角等三个因素,建立三因素四水平(L16(43))(表1)的正交试验。通过拟合滑坡现场滑距和滑体前缘堆积范围及厚度,最终确定滑坡灾情发生时物理力学参数(表2)。
表1 正交试验参数及取值表
根据现场对滑坡体的调查和判断,拟合滑坡变形和标志物运动情况,通过参数反演最终确定灾情发生时结构面强度参数取值(表2)。同时,考虑计算模型简化结果,将滑坡参数分为岩体和结构面参数两类。通过改变结构面强度参数模拟实际降雨对节理和层面等软弱面的弱化情况,具体表1中结构面1表示天然工况下节理和滑动面强度参数取值;结构面2表示天然+降雨工况下节理和滑动面强度参数取值。
表2 滑坡力学参数取值
边界范围的圈定以现场调查为主,以滑坡变形和直接影响范围为主要指标,主要考虑位移边界条件,具体设置如下:模型的左右边界和底部边界设为固定边界,上部自由边界以模拟滑体在自重作用下的滑移。
2.2.1稳定性分析
经现场踏勘和对比历史资料,天然工况下结构面强度参数如表2所示。DDA数值程序能实现对滑坡关键块体位移和应力应变的实时跟踪,因为滑坡前缘存在高10余米的临空面,即便在滑坡物理力学参数很大的情况下,依然可能出现滑体前缘块体的局部运动,不利于反映滑体稳定的真实状况;滑体后缘块体在形成拉裂面后也存在块体局部稳定的可能,亦不能真实反映滑体发生不稳定滑动的整个过程稳定性系数变化,故综合考虑滑坡现场地形和滑坡变形运动基本特征,选取滑体中部11#块体为关键块体较为合适,能比较真实反映滑体稳定状况。
天然状态下,关键块体稳定性系数如图9所示。初始状态滑体稳定性系数fs约2.4左右,在约45 s监测时间内,稳定性系数由急剧降低到迅速升高到3.0,之后呈阶梯式下降,fs值总体均大于1,说明滑坡在天然0状态整体稳定。整个过程,滑体最大速度0.6 m/s,最大滑距0.63 m(图10)。总体来看,fs的不断减小说明滑体天然状态下存在缓慢蠕动,存在滑坡隐患,在诱发因素作用下极有可能运动加剧。模拟结果与现场应急调查掌握情况一致,说明了天然状态岩体和结构面强度参数取值的正确性和合理性。
图9 天然工况11#块体速度及稳定性系数与时间关系曲线Fig.9 Relationship curve of 11# block about speed and stability coefficient with time in natural conditions
图10 天然工况11#块体滑距-时间关系曲线Fig.10 Relationship curve of 11# block about speed with time in natural conditions
降雨工况下,通过对滑体前中后部分别布置监测块体,对各监测块体的稳定性系数、速度和运动距离进行动态监测(图11),其中表征滑体稳定性的11#关键块体的具体布置位置见图8。分析结果显示滑坡运动前fs在1.5左右,滑坡启动后fs急剧降低,说明滑坡进入快速滑动阶段;在20~50 s时间范围内,fs出现振荡上升,说明此时滑体处于缓慢堆积停止动态过程。从图11中可以看出,11#关键块体的稳定性系数和速度大小的变化呈现此消彼长的关系,很好的反映了滑坡从启动、滑动到堆积停止的整个动态发展过程。
图11 降雨工况11#块体速度及稳定系数时间曲线Fig.11 Curve of speed and stability coefficient with time in rainfall condition
2.2.2结构面强度参数取值影响性分析
在降雨工况结构面2参数取值前提下,利用控制变量法分别单独研究降雨工况下结构面强度参数黏聚力c和内摩擦角φ降低对滑坡滑距和速度的影响。从图12和图13可以看出,在0~15 s时间范围内滑距随时间变化曲线基本一致,而最大滑距则与c值呈负相关;而滑体运动最大速度总体上随着c值的升高而降低,同样在初始0~10 s时间范围内速度-时间曲线趋势一致,局部速度呈现上下一定幅度的变动,很好的反映了实际块体运动过程中的平移、翻转和碰撞等运动形式。分析图14和图15结果,可以发现随着φ值的增加,最大滑距和最大速度均不断减小,且滑体运动初期滑距和速度变化情况差别明显,说明φ值变化对滑体影响很大,消耗滑体动能作用明显;此外,从图14看出φ取0°时,在滑体达到最大滑距之后出现小幅度的降低,同一时间段速度出现负值,很好的反映了在降雨导致结构面强度参数极度弱化的情况下,滑体前缘遇对面山体阻碍发生逆冲反翘,导致中后部块体发生反向运动堆积,与现场实际情况一致。总体来看,c值变化对滑坡滑距和速度影响有限,φ值对滑体运动影响显著,这一结果与葛云峰等[11]研究相同。
图12 c值变化滑距-时间曲线Fig.12 Curve of sliding distance with time by the change of cohesion
图13 c值变化速度-时间曲线Fig.13 Curve of speed with time by the change of cohesion
图14 φ值变化滑距-时间曲线Fig.14 Curve of sliding distance with time by the change of angle of internal friction
图15 φ值变化速度-时间曲线Fig.15 Curve of speed with time by the change of angle of internal friction
2.2.3滑动过程分析
数值模拟结果显示,模拟结果表明滑坡运动过程分为启动、快速滑动和堆积三个阶段(图16):启动阶段(图17),由于滑体前缘10余米的临空面,导致滑体前缘块体首先发生滑动,随后滑体后缘出现拉裂面,滑体整体向前滑动;快速滑动阶段,滑体表层块体松动,由前部逐渐向中后部扩展,部分块石飞向空中,模拟结果很好的反映了块体飞起和下落堆积过程,布置在滑体前、中、后三个部位的监测单元速度-时间曲线显示滑体最大速度约在6.6 m/s;堆积阶段(图18),滑体局部最大位移92 m左右,结合实际模拟堆积结果,滑体整体滑距约38 m,与现场调查滑体后缘滑距约40 m保持一致,前缘堆积厚度53 m,滑动持续时间70 s左右。模拟滑坡变形运动过程与工程地质调查机制分析结果基本一致,验证了滑坡概念模型的确性,有助于理解灾害发生过程及机制。
图16 滑坡滑动过程变形图[18]Fig.16 Deformation pattern of landslide sliding process
图17 降雨工况各监测块体滑距-时间曲线Fig.17 Curve of sliding distance with time in rainfall condition including all monitoring blocks
图18 降雨工况各监测块体速度时间曲线Fig.18 Curve of speed with time in rainfall condition including all monitoring blocks
(1)现场调查发现王山-抓口寺滑坡中下部玄武岩体在被发育的X型节理切割和地震共同作用下结构破碎,滑坡裂隙发育且张开度好,降雨入渗导致滑坡结构面强度参数严重弱化为滑体滑动过程的控制性影响因素,局部出露泥化页岩层面,滑床上部出露基岩面即为滑坡整体滑动面;且滑体两侧发育大量冲沟有利于雨水侧向入渗,加速滑坡变形;此外,滑坡前缘高达十余米的陡坎,提供了绝佳的临空条件,在上述几种因素的共同作用下,滑体发生以“强降雨-入渗-岩体软化-推拉”破坏模式为主的典型滑移-拉裂顺层岩质滑坡滑动。
(2)综合现场地质调查和数值模拟结果,通过对比分析滑坡滑动过程和堆积状况,可以发现对滑体包含碎石土、节理切割破碎玄武岩和薄层凝灰岩滑带简化为岩体和结构面两部分符合客观实际,达到了预期目的。
(3)结构面强度参数取值的影响性分析显示内摩擦角φ和黏聚力值c的大小与滑坡滑距和最大滑速呈负相关,且内摩擦角φ对灾害体位移和速度影响显著,消耗滑坡动能作用明显,相同条件下黏聚力则不明显。
(4)天然状态滑坡整体基本稳定,但滑体部分块体发生缓慢蠕动,存在安全隐患;滑坡在降雨诱发下发生滑动,滑动过程中最大速度6.6 m/s,滑体整体滑距38 m,滑体前缘堆积厚度53 m,持续时间70 s左右,模拟结果与现场调查基本一致。