基于D-InSAR和PFC2D技术的白格滑坡稳定性分析
靳立周,1), 王盈2), 常文斌3), 田颖颖1), 袁仁茂,1),*
1)中国地震局地质研究所, 北京 100029
2)中国地震应急搜救中心, 北京 100049
3)中国地震局兰州地震研究所, 兰州 730000
摘要
2018年10月11日和11月3日, 位于金沙江上游右岸的西藏江达县波罗乡白格村先后发生了2次特大规模的滑坡堵江事件, 经人工干预后险情得以解除。运用差分干涉合成孔径雷达(D-InSAR)技术对白格滑坡后缘的潜在危岩体进行了7.5d连续不间断的监测, 结果表明: 滑坡体上侧区域目前仍然不稳定, 具有再次发生滑坡堵江的危险。基于此, 文中利用颗粒流软件PFC2D模拟了滑坡后缘潜在危岩体在自身重力、 强降雨、 地震条件下的稳定性状况。模拟结果表明: 后缘潜在危岩体在静力作用下不会产生明显的失稳滑动, 在强降雨和强地震动条件下会发生失稳破坏, 可能会再次堵塞金沙江并形成堰塞湖。根据文中的模拟结果可对滑坡稳定性做出科学评价, 同时为以后类似滑坡的防灾减灾提供参考。
引用格式
靳立周, 王盈, 常文斌, 田颖颖, 袁仁茂. 基于D-InSAR和PFC2D技术的白格滑坡稳定性分析[J]. 地震地质, 2023, 45(1): 153-171 DOI:10.3969/j.issn.0253-4967.2023.01.009
JIN Li-zhou, WANG Ying, CHANG Wen-bin, TIAN Ying-ying, YUAN Ren-mao. STABILITY ANALYSIS OF THE BAIGE LANDSLIDE USING D-INSAR AND PFC2D MODELING[J]. Seismology and Geology, 2023, 45(1): 153-171 DOI:10.3969/j.issn.0253-4967.2023.01.009
0 引言
2018年10月11日, 受大量和持续降雨产生的影响, 西藏江达县白格村发生特大规模山体滑坡, 大量滑坡体下滑堵塞了金沙江干流, 形成堰塞坝, 而后湖水上涌, 淹没了多条道路。2日后, 堰塞坝通过自然泄流形成泄流槽, 险情得以解除。2018年11月3日, 滑坡后缘的大裂缝处再次发生滑动破裂, 下滑的碎屑物质沿途刮铲所经岩体形成规模巨大的碎屑流, 碎屑流高速下滑后再次堵塞金沙江, 形成堰塞坝。2次滑坡形成的坝体总方量约为3.2×107m3(冯文凯等, 2019)。第2次滑坡所形成的堰塞坝高度比第1次高50m (郭晨等, 2020), 堰塞湖的蓄水量预测为7.9×108m3 (陈祖煜等, 2020), 对上、 下游人民的生命财产造成了极大的安全隐患。通过人工挖掘泄流通道、 修建导流槽, 坝体上、 下游得以贯通, 消除了溃坝危险。
白格滑坡属于高位-远程滑坡, 滑坡剪出口高于坡脚地面, 滑坡体从高陡斜坡上部剪出并加快跌落, 同时不断铲刮下部岩体, 使体积明显增加, 伴随着滑坡体的瓦解粉碎, 进而转换成快速远程碎屑流滑动。高位滑坡通常具有运移速度快、 能量大、 破坏力强的特点, 并逐渐发展成高位启动、 高速远程和灾害链的成灾模式(胡卸文等, 2009; 许强等, 2009; Hao et al., 2014; Yuan et al., 2014; 殷跃平等, 2017)。除滑坡带来的危害外, 滑坡堵江本身及其伴随的堰塞湖溃决问题可能会对滑坡区上游和下游地区产生严重的影响(柴贺军等, 2000)。滑坡发生后, 堰塞湖上游水位持续上升, 江达县的波罗乡、 白玉县金沙乡等位于金沙江上游的乡镇被淹。第2次泄洪后金沙江发生较大规模洪水, 使四川、 云南等中下游部分临江区域的城镇被淹, 公路桥梁被冲毁, 给人民群众的生命财产和水电站等基础设施的安全带来了极大威胁。
国内众多学者对白格滑坡的形成机制和滑坡残留危岩体的稳定性进行了相关研究。冯文凯等(2019)认为白格滑坡后缘受到分支断层控制作用的影响, 在长期重力卸荷、 降水和地下水不断浸润的综合影响下, 山体边坡造成累积时效性形变, 最终整体失稳下滑。王立朝等(2019)认为金沙江构造断裂带岩性破碎及特殊的高山峡谷地形地貌是白格滑坡产生的基础地质条件, 地震、 冻融和降雨是本次滑坡的诱发因素。邓建辉等(2019)对滑坡后缘裂缝区的危险性展开分析, 发现滑坡后缘牵引区现阶段出现较为严重变形的方量约为5.5×106m3, 存在再次发生滑坡和堵江的风险。Li等(2020)基于地基雷达变形观测和星载SAR数据像元偏移跟踪, 发现2次滑坡后残余的岩体仍然在移动。周礼等(2019)利用PFC3D对2次滑坡的产生、 运动和堆积过程进行了数值模拟, 在反演结果的基础上, 模拟和预测了白格滑坡滑源区剩余潜在不稳定部分未来的运动路径和堆积范围。本文将结合野外现场调查资料、 边坡雷达和滑坡后缘潜在危岩体的形变监测数据, 利用离散元软件PFC2D模拟白格滑坡后缘危岩体在强降雨、 地震条件下的失稳破坏情况。所得结果表明, 滑坡后缘危岩体在上述情况下会继续发生失稳破坏, 可能再次堵塞金沙江并产生堰塞湖。根据本文的模拟结果能够对滑坡的稳定性做出科学评价, 也可为灾害发生后应采取的应急措施和防灾减灾提供理论指导。
1 滑坡区地质环境条件
白格滑坡位于西藏自治区江达县波罗乡白格村金沙江右岸的山脊上, 滑坡体的中心坐标为(31°4'56″N, 98°42'18″E)(图1), 滑坡区顶部高程为3700m, 河谷高程为2880m, 河流侵蚀严重, 河谷强烈下切呈“V”形, 属构造侵蚀高山峡谷地貌。滑坡区为大陆性季风高原型气候, 多年平均气温为7.5°, 历年平均降水量为600mm。降水量时空分布不均, 大约80%的降水集中在610月之间。
图1 白格区域构造背景图(改自Ren等, 2018)
红框为白格滑坡发生区域
白格滑坡区的地表被第四纪沉积物所覆盖, 滑坡周边主要出露三叠系和元古代地层(图2), 岩性为: 三叠系上统金古组(T3jn)灰岩带, 三叠系上统松多组(T3x2)的碳酸盐岩和碎屑岩, 石炭系上统生帕群(C2sh)岩屑砂岩、 粉砂岩、 灰岩等, 元古界雄松群(ptxna)片麻岩组, 华力西期蛇纹岩带(φw4), 燕山期戈坡超单元(ηγ)二长花岗岩和燕山期则巴超单元(γδ)花岗闪长岩, 滑坡顶部和前缘出露地层的岩性大多为元古界雄松(ptxna)片麻岩组, 产状为 235°∠40°, 滑坡中部的滑床区域出露抗剪强度较差的华力西期蛇纹岩。受周边断层长期活动的影响, 滑坡区坡顶的片麻岩破碎, 岩体节理裂隙发育, 地层结构面复杂, 主要发育2组片麻岩节理面, 产状为60°~80°∠75°~85°、 100°~115°∠80°(许强等, 2018)。
图2 滑坡区地质图
滑坡区地处金沙江缝合带上, 构造活动强烈, 区域地质环境复杂, 由一系列NW向的褶皱和断层组成。断层以逆断层为主, 主要有雪青-龙岗断裂(F4)、 竹英-贡达断裂(F3)、 波罗-木协断裂(F2)和贡达-迪中断裂(F1)。其中, 波罗-木协断裂(F2)刚好从滑坡源区顶部通过, 对滑坡的影响较大。
2 滑坡概述
2.1 2次滑坡堵江事件的基本特征
白格滑坡的平面状态具有明显的圈椅状特点: 滑坡的后缘和前缘较陡, 而滑床所在的中间区域则比较平缓(图3), 总体平均坡度为40°~65°, 由斜坡顶部到底部呈现陡缓陡的特征。滑坡剪出口的高程约为3000m, 滑坡后缘的高程为3700m, 坡脚高程为2850m, 高差达850m。白格滑坡-堰塞体的地质剖面如图4所示, 巨大的高差使得滑坡体产生很大的势能, 受到地形影响, 滑坡高位下滑后在重力作用下高速运动, 滑坡体以翻滚、 碰撞、 滑跃的运动形式快速堆积在金沙江河道, 在高速撞击金沙江左岸后迅速爬高, 并向右岸回落堆积, 堵塞金沙江形成堰塞坝。第1次滑坡下滑岩土体的方量约为2.4×107m3, 形成的堰塞坝平均高度为120m, 蓄水量达2.9×108m3。滑坡发生2日后, 堰塞坝发生漫顶溢流, 逐渐发展成泄流槽, 险情得以解除(Ouyang et al., 2019a)。
图4 白格滑坡地质剖面图(II')
第2次滑坡主要发生在第1次滑坡源区的后缘陡壁处, 滑坡后缘裂缝由于局部滑动发生贯通, 边坡后缘产生崩塌, 整个滑坡体迅速下滑, 滑坡体底滑面上堆积的松散物质被刮铲而形成凹槽状, 滑坡体解体后形成碎屑流, 快速下滑的碎屑流堆积于第1次滑坡形成的堰塞体上, 形成堰塞坝, 再次阻断了金沙江。第2次滑坡下滑岩体的总方量为8.5×106m3, 蓄水量>5×108m3。后由政府部门组织修建导流槽进行人工泄洪, 消除了溃坝危险。
2.2 滑坡后缘潜在危岩体的变形特征
通过对历史遥感影像进行分析可知, 滑坡区岩体的变形过程至少经历了50a(许强等, 2018)。2011年, 斜坡的后缘和前缘发生变形, 后缘的裂缝逐渐发生贯通形成拉裂面; 20112015年, 滑源区斜坡后缘和前缘的变形区域面积有所增加, 但没有发生剧烈的位移变化; 2018年9月, 卫星影像显示, 滑坡后缘产生明显的整体位错, 滑源区中部有较大规模的滑塌; 滑坡前缘形成了明显的剪出口, 最后于2018年10月11日斜坡发生失稳破坏, 造成滑坡堵江事件。第1次滑坡发生后, 在滑坡后缘及两侧发现3处潜在的不稳定岩体(图5), 同时, 第1次滑坡又为第2次滑坡的发生提供了陡峭的临空面。在2次大规模的滑坡发生后, 根据雷达监测发现滑坡后缘的3处潜在不稳定体仍在持续变形, 具有明显的活动性, 存在再次堵塞金沙江的风险(赵程等, 2020)。
3 边坡监测结果分析
3.1 D-InSAR技术概述
差分干涉合成孔径雷达(D-InSAR)技术利用SAR图像去除干涉图中的非形变相位, 以实现地表形变的高精度测量。该技术不会受到云、 雾天气的影响, 可获取大范围、 高精度、 全天候、 全天时的形变信息。它有能力探测到厘米甚至毫米级的地表微小变形, 能够获得基于面范围的地面形变信息, 并能获取边坡面区域的连续形变图(王桂杰等, 2010; 齐麟等, 2020), 从而直观地看到边坡产生形变的全过程。目前, 其已被广泛应用于地震活动、 火山灾害、 地面沉降、 滑坡等地质灾害的形变过程监测中(Ouyang et al., 2019b)。
3.2 监测区的数据采集及分析
一般在山体滑坡发生后, 边坡体仍然处于不稳定阶段。为了进一步验证现阶段边坡的稳定性, 我们在白格滑坡发生7个月后使用差分干涉合成孔径雷达对滑坡后缘的潜在不稳定岩体进行了现场监测。根据白格滑坡监测边坡区域的分布情况, 将雷达的监测点放置于边坡正对面(图6), 全天24h无间断监测, 通过扫描滑坡区的三维地形, 对边坡的形变进行三维实时监测, 监测期从2019年6月4日上午8点开始至2019年6月11日21点结束, 监测期间平均10min获取1轨数据, 共采集轨数为1084轨。截至2019年6月11日, 雷达监测图显示滑坡面的局部区域出现明显的形变, 在滑坡体上侧区域发生小面积滑移, 从雷达监测图中可以明显分辨出滑坡区(图7)。雷达监测过程的时间序列结果可有效反映出滑坡体的渐进变化(图8)。
图8 金沙江白格滑坡地基雷达形变时间序列图
从连续监测结果(图7)来看, 目前滑坡后缘的潜在危岩体依然存在形变, 7.5d的连续监测结果显示场景的最大形变量为230mm(雷达监测结果为负值), 平均形变速率为30mm/d。从图7可以看出, 整个边坡体存在4个主要的形变区域, 均位于滑坡体的中上缘, 2号区域的形变量相对其他位置较大。
一般而言, 重点监测区域的重点点位可以代表某区域的形变特点。为此, 我们选择场景中1、 2、 3、 4区域内的点位进行历史形变分析(表1, 图6)。从4个区域点位的累积形变结果看, 这些点位的累积形变量值均约为200mm。其中, 区域2的形变量最大, 达227mm; 区域4的形变量最小, 达170mm(图9)。为了更好地描述打点(即雷达布设点位)累积形变量的趋势, 我们对形变曲线进行了滤波平滑处理, 滤波结果可较好地代表观测值的趋势, 图10是打点位置的形变速率曲线, 可以看到打点位置区域的形变速率几乎在-2~0mm/h附近波动。从几个单点的历史形变曲线和速率曲线可以看出, 点位的形变趋势基本呈现直线缓慢下降, 并没有出现形变加速增大的现象, 表明目前滑坡体没有进入临滑状态。
通过对金沙江白格滑坡进行现场雷达监测发现, 在观测期2019-06-042019-06-11的7.5d内, 滑坡体的上侧区域发生了小面积滑移, 几个点位均揭示出量值约为200mm的接近雷达方向的形变, 顶部的形变速率可达30mm/d。这表示在金沙江白格滑坡发生后7个月, 滑坡体仍然不稳定, 可能会再次发生滑坡, 堵塞金沙江。因此, 有必要对滑坡后缘潜在危岩体的运动及堆积过程进行数值模拟, 并对边坡的稳定性进行科学分析。
4 滑坡危岩体数值模拟分析
4.1 PFC数值模拟方法简介
颗粒流法(Particle Flow Code, PFC)是于20世纪70年代提出的基于通用离散单元模型(DEM)框架的一种定量评价方法。PFC程序定义颗粒为有质量和表面刚性的圆形或球体, 允许离散颗粒产生位移和旋转, 伴随着计算过程的进行可自动检索新的接触, 颗粒球能够通过特定的黏结模型组合在一起(Cundall et al., 1992; Manzella et al., 2013)。该方法的基本核心是以力-位移定律和牛顿第二定律为基本理论模拟颗粒的运动和相互作用, 通常用颗粒接触本构模型来描述颗粒之间的相互作用。
该方法不会受到变形量的限制, 可方便、 快速地解决非连续介质力学问题, 使得宏观岩土体的变形破坏特征可根据细观颗粒的运动特征得到体现, 并有效模拟介质的开裂、 分离等非连续现象, 能直接反映岩土体的破坏机理、 运动过程和结果, 对大规模运动条件下的远程滑坡模拟具有较好的应用性。
4.2 模型细观参数的确定
颗粒流法将颗粒的细观参数赋予岩土体, 包括颗粒属性、 摩擦系数、 切向刚度、 法向刚度等(周健等, 2009)。正确选取模型的细观参数值是获取精确模拟结果不可缺少的因素, 也是进行滑坡稳定性数值模拟分析的关键。
室内物理实验得到的宏观参数与PFC颗粒流模型中的细观参数并没有直接的量化关系式, 通过室内物理试验得到岩土体的宏观力学参数(内摩擦角, 黏聚力等)无法直接应用于颗粒流数值模拟中, 因此需要在PFC中采用单轴压缩试验对模型的细观参数进行标定, 以获取颗粒的细观力学参数(石崇等, 2012)。所得PFC中的单轴压缩模拟曲线如图11所示, 数值模拟细观参数和室内物理试验宏观响应的关系是非线性的, 模拟所得结果与室内测定的岩石强度参数接近(表2)。最终数值模拟模型所用颗粒的细观参数见表3。
4.3 模型细观参数的确定
通过对白格滑坡的滑动特征进行现场调查, 我们发现在滑坡后缘附近存在2次滑坡滑动后形成的滑动陡坎面及拉张裂缝, 同时滑坡滑床下部的岩体由于被刮铲而形成一个位置较好的临空面。经现场雷达检测数据分析发现, 位于滑坡体后缘滑动陡坎面附近的2号区域和3号区域的形变量相对其他位置较大, 虽然4号区域的形变量较小, 但其地形较陡、 坡度较大, 具有较好的势能条件, 在外界因素的扰动下容易发生垮塌。因此, 我们选取AA'和BB'为研究剖面进行数值模拟(图12), 发现AA'和BB'剖面的模拟结果相似。我们以AA'剖面的数值模拟为例进行讨论, 并分析边坡的稳定性。
滑坡数值模型建立在野外实地勘察的基础上。本文采用PFC中的平行黏结模型进行数值模拟, 并根据Ball-Ball法构建滑坡模型。Ball-Ball法的优点是可以模拟滑坡的启动及滑动过程, 有利于更好地理解滑坡的失稳及滑动机制(曹文等, 2017)。建立二维Ball-Ball模型具体步骤如下:
(1)根据AA'剖面在AutoCAD中绘制 1︰1 边坡几何模型, 并保存成DXF文件, 在剖面图中设置2道先前滑移造成的滑裂面。
(2)通过PFC2D中的“几何导入”命令将DXF文件导入, 生成滑坡边界墙体。
(3)在相应的范围内初步生成颗粒, 根据表2中取得的颗粒细观力学参数对颗粒进行赋值, 模型平衡后, 再通过削坡的方式删除超出范围的颗粒, 得到Ball-Ball模型(图13)。
4.4 边坡稳定性分析
在模拟过程中, 首先通过上述模型及参数清除初始边坡模型中颗粒的位移和速度, 然后对颗粒施加重力作用, 模拟边坡在静力作用下的稳定性。图14为边坡在静力作用下最终变形结果的位移云图, 从图14左侧的位移标尺可以看出, 边坡模型的主体位移量在200~600mm之间, 通过地基雷达监测结果来看, 目前现场监测到的累积形变量值大部分约为200mm, 而数值模拟监测到的200~600mm变形量是边坡变形演化的最终结果。模拟结果表明, 模型在静力作用下只是缓慢变形, 不会发生明显的失稳滑动, 边坡的初始变形量与雷达监测的位移结果基本一致, 表明模拟结果较为合理。图15为边坡模型在静力作用下的位移趋势图, 可以看出模型在滑裂面附近和右侧边坡的位移量较大, 说明该斜坡的滑体在静力作用下仍然具有滑动趋势, 并不稳定, 在强降雨、 地震等诱发因素的影响下, 当边坡强度超过临界状态时, 该边坡可能会发生失稳破坏。因此, 本文将继续考虑强降雨和地震2种工况对边坡的破坏情况, 进一步分析边坡的稳定性。
采用强度折减法进行降雨作用下的边坡稳定性模拟。强度折减法是利用数值分析方法分析边坡稳定性的一种简单有效的方法, 其基本思想主要是拆减介质强度参数(主要是黏聚力和内摩擦角), 当系统达到不稳定状态时得到的折减系数就是边坡的安全系数(雷远见等, 2006)。强降雨能够使边坡内部产生更高的孔压, 导致边坡更容易发生失稳破坏。根据相关室内试验结果, 强降雨状态可使碎块石土颗粒的黏结强度参数折减22%、 摩擦因数折减20%(梁鑫等, 2011)。基于此, 对离散元模型中颗粒的细观参数进行调整。添加该降雨折减条件后, 根据前述方法建立强降雨工况下边坡失稳破坏的预测模型, 模拟边坡在强降雨条件下的变形破坏过程。循环时步至边坡达到稳定状态终止, 模拟记录的边坡破坏过程如图16所示。
从图16中可以看出, 当计算运行到6×104时步时, 两侧边坡在重力作用下整体开始发生向下的位移变形; 运行到3.6×105时步时, 左侧边坡中上部滑裂面附近的抗剪强度较弱, 岩体较破碎, 向坡脚方向发生失稳下滑, 位移明显; 运行到9.0×105时步时, 两侧坡体的颗粒冲破坡脚的锁固作用, 颗粒向前滑出; 运行到1.68×106时步时, 右侧边坡上缘发生局部剪切破坏, 产生拉张裂缝; 运行到1.74×106时步时, 右侧边坡上缘失去支撑, 发生失稳破坏, 颗粒逐渐在坡脚处堆积; 最终, 在3.6×106时步时边坡基本停止运动, 颗粒堆积于坡脚处。模拟结果表明: 随着时步增加, 边坡的变形位移量也在不断增加, 在运行过程中, 坡体逐渐出现了拉张裂缝, 不断发展直至被剪切贯通, 发生失稳破坏。
跟据《中国地震动参数区划图》(GB18306-2015), 白格滑坡区50a内超越概率10%的地震动峰值加速度为0.20g, 对应的地震基本烈度为Ⅷ度, 地震动加速度反应谱的特征周期为0.40s。据《水电工程区域构造稳定性勘察技术规程》(NB/T35098-2017)区域构造稳定性评价标准可知, 白格滑坡的区域构造稳定性较差(赵永辉等, 2020)。基于此, 本文在模拟地震工况部分采用2013年岷县-漳县地震的地震波, 其垂直向地震动峰值加速度为0.21g, 水平向地震动峰值加速度为0.096g, 原始地震动数据如图17所示。在离散元数值模型中, 地震动以速度的形式从边坡模型底部输入, 图18为0.21g峰值加速度作用下的边坡破坏情况。
地震动力的失稳形式不同于一般条件的静力失稳, 一旦出现裂缝便迅速扩展(刘婧雯等, 2014)。在地震动作用下, 2~6s时, 右侧边坡后缘先产生裂缝并迅速拉裂, 滑坡体刮铲下部的颗粒一起整体向前加速下滑并解体, 局部出现块体抛射现象; 8~10s时, 右侧边坡顶部变形加剧, 开始出现新的拉张裂缝, 发生失稳下滑; 10~20s右侧边坡下滑的颗粒高速冲上左侧边坡, 在颗粒碰撞之后发生回弹现象, 左侧斜坡出现小范围松动滑动, 但整体基本处于稳定状态。模拟结果表明: 地震作用对边坡稳定性的影响主要体现在边坡右侧部分, 主要是由于右侧边坡的坡度较陡, 具有良好的临空面条件, 在地震作用下容易发生崩塌、 滑落现象。当模型运行时, 右侧边坡后缘先产生拉张裂缝, 随着坡体中部发生剪切破坏与后缘的拉张裂缝发生贯通, 坡体整体发生失稳破坏。左侧边坡的坡度相对较缓, 势能条件较差, 在施加地震波条件下, 坡面有部分岩体破碎滑落, 但整体上基本处于稳定状态。
5 结论
(1)通过差分干涉合成孔径雷达连续7.5d的24h不间断监测, 在滑坡体上缘识别出4个主要的形变区域, 并发现滑坡体上侧区域发生小面积滑移, 选取的4个测量点位均揭示出量值约为200mm的形变。顶部的形变速率可达30mm/d, 表明在金沙江白格滑坡发生的7个月后, 滑坡体仍不稳定。
(2)运用颗粒流离散元软件PFC2D对变形区域的边坡进行了模拟分析, 揭示了其在自身重力、 强降雨及强地震动条件下的变形破坏特征。结果表明, 滑坡后缘的潜在不稳定危岩体在强降雨、 强地震动等诱发因素下仍然会发生失稳破坏, 应当采取适当削坡等工程措施对潜在危岩体进行治理, 并继续加强对滑坡的实时监测, 设置监测预警系统, 及时消除潜在危岩体下滑并再次堵塞金沙江的隐患。
转自：“测绘学术资讯”微信公众号
如有侵权，请联系本站删除！
