四川茂县新磨村高速滑坡数值模拟研究

Numerical Simulation of Xinmo Village Expressway Landslide in Mao County, Sichuan Province

摘要: 2017年6月24日上午6时左右，四川省茂县叠溪镇新磨村发生特大高速远程滑坡。这次山体滑坡掩埋了64座房屋，造成10人死亡，另有73人失踪。这是近年来中国发生的最严重的灾难性山体滑坡事件之一，引起了国内外公众和学术界的广泛关注。高速远程滑坡是一种典型的滑坡，因为运动速度快，距离长而得名，其高速启动机理和运动特征一直是地质灾害领域研究的热点问题，而四川茂县新磨村滑坡就属于典型的高速远程滑坡–碎屑流。本文以四川省茂县叠溪镇新磨村滑坡为研究对象，通过资料收集、详细的现场调查和数值模拟，对滑坡的特征和成因机制进行了简单的分析和研究。其次，建立了新磨村滑坡的地质模型，借助DAN-3D数值软件，进行滑坡运动过程模拟。将模拟后的运动特征与堆积形态与真实情况对比，发现两者基本一致，由此得出结论，DAN-3D数值模拟可以较为真实地还原滑坡运动过程，从而可以对其它滑坡进行早期模拟与预测，以期为滑坡的识别和防治提供技术支撑。

Abstract: At around 6 am on June 24, 2017, a huge high-speed remote landslide occurred in Xinmo village, Diexi town, Mao County, Sichuan province. The landslide buried 64 homes, causing 10 deaths and leaving 73 others missing. It was one of the worst catastrophic landslides to occur in China in recent years, and has attracted wide attention from the public and academia at home and abroad. High speed and long distance landslide is a typical type of landslide, which is named for its fast moving speed and long distance. Its high speed starting mechanism and moving characteristics have been a hot topic in the field of geological disasters. The Xinmo village landslide in Mao County, Sichuan Province is a typical high-speed and long-distance landslide-debris flow. Taking Xinmo village landslide in Diexi Town, Mao County, Sichuan Province as the research object, this paper makes a simple analysis and study on the characteristics and genetic mechanism of the landslide through data collection, detailed field investigation and numerical simulation. Secondly, the geological model of Xinmo village landslide is established, and the landslide movement process is simulated with the help of DAN-3D numerical software. Comparing the simulated motion characteristics and accumulation patterns with the real situation, it is found that the two are basically consistent, and the conclusion is drawn. DAN-3D numerical simulation can restore the real landslide movement process, so that other landslides can be simulated and predicted in the early stage, in order to provide technical support for landslide identification and prevention.

1. 引言

1.1. 选题的研究意义

由于我国西南部地区高原、山区地形起伏较大，加上地震发生频繁、降雨量多且集中等因素影响，滑坡灾害频发，尤其是云贵川地区。高速远程滑坡属于一种典型的滑坡形式，这种滑坡具有超高的运动速度和超远的滑动距离，往往会给人类财产和生命安全造成巨大的损失，所以其启动机理和运动特征一直是国内外地质学术界研究的热点和难点问题。本文以四川发生的新磨村高速滑坡为研究对象，通过DAN-3D数值模拟软件来反演滑坡发生过程，以期为这类滑坡的研究提供技术支撑。

1.2. 国内外研究现状

1.2.1. 国外研究现状

国外对高速远程滑坡的研究起步较早，在经过早期的一些学者的描述和推断后，主要形成了4种模型与理论，分别是空气润滑模型[1]、颗粒流模型[2]、能量传递模型[3]、底部超孔隙水压力模型[4]。

随着计算机技术的不断发展，对高速远程滑坡的研究方法和技术也有了很大的进步，除了现场调查、统计对比、模型试验等方法，最主要的就是数值模拟方法。

关于数值模拟，最早是由Lang建立了模拟碎屑流运动过程的连续流动模型，这个模型可以模拟滑坡的初始形态、变形特征和堆积特征等[5]；

Dent写了BVSMAC的程序并反演了滑坡–碎屑流的运动过程，并且创建了双线本构模型[6]；

Hungr在拉格朗日差分法的基础上创建了DAN动力模拟方法，可以使用不同的本构方程来模拟滑坡运动[7]；

Stephen G. Evans用DAN模拟了Mount Cayley滑坡过程，其中Voellmy模型模拟结果与真实情况最为接近[8]；

总的来说，高速远程滑坡的研究发展过程缓慢，迄今取得了一些成就，但由于其本身的复杂性，其机理还不能明确。

1.2.2. 国内研究现状

国内关于高速远程滑坡的研究相对较晚，但国内的学者在此方面还是做了大量的研究。胡广韬首次分析了高速滑坡–碎屑流的运动特征，提出了“碎屑流的持续效应”、“临床弹性冲动加速效应”等[9]；周迎庆研究了1988年西宁滑坡的形成演变及高速滑动机理，将其破坏过程分为塑流–拉裂变形、蠕滑–拉裂变形、岩体剪断，滑移面贯通等阶段[10]；殷跃平对2000年西藏林芝地区波密县易贡藏布河发生的巨型高速滑坡进行研究，提出了一些减灾的工程措施和方案[11]。张龙对重庆市鸡尾山高速远程滑坡的运动过程进行了PFC3D模拟，发现在一定数值的摩擦系数与粘结强度时，模拟结果较为合适[12]。费建波等人引入颗粒底部摩擦模型，并采用圣维南方程模拟高速远程滑坡，最终发现用颗粒流模型的计算结果能更好地反映实测的滑坡体最终堆积情况[13]。

由于高速远程滑坡的体积大、具有“尺寸效应”、运动机制的影响因素的复杂性，高速远程滑坡的运动机理和运动规律还不能明确，所以高速远程滑坡的研究还有很长的一段路要走。

2. 工程地质概况

2.1. 地形地貌

新磨滑坡位于青藏高原与四川盆地的过渡地带，主要发育岷山隆起，除东南边界为龙门山的尾段外，大部分地区属于邛崃和岷山。主要发育的河流为松坪沟河，松坪沟河是岷江的支流之一，自西北向东南流经松坪沟乡。黑水河、赤不苏河和松坪沟河分别在大两河口、小两河口和叠溪镇流入岷江。松坪沟河两岸广泛分布斜坡和湖相沉积。滑坡沉积物主要分布在海拔2260~2400 m的松坪沟河道内。

2.2. 地层岩性

研究区露头地层主要为三叠系中统杂古脑组(T2z)，主要出露基岩为灰色厚层块状石英砂岩、千枚岩和板岩，岩层倾向与坡向一致，说明松坪沟左岸为顺层斜坡。由于构造挤压作用，地层呈波状变形，千枚岩、板岩等中间断面形成易滑动的软弱夹层，造成斜坡不稳定，易产生滑坡灾害。受密集的节理和裂缝的影响，沿这些节理和裂缝，风化作用非常活跃。

2.3. 地质构造

研究区位于较场弧形构造带的前弧西翼，主要活动断裂构造为岷江隆起的东界–岷江断裂、西界–虎牙断裂和位于弧形西翼的松坪沟断裂。滑坡区构造上处于NWW向松潘–甘孜褶皱带、近EW向西秦岭构造带与NE向龙门山断裂带的三角交汇部位，发育一系列紧密线状弧形同斜倒转褶皱[14]。研究区构造活动强烈，褶皱发育。

3. 新磨滑坡数值模拟分析

3.1. DAN-3D介绍

DAN-3D (Dynamic Analysis of Landslide in Three Dimensions)最早是由加拿大学者Hungr提出并制作的等效流体动力分析软件，DAN-3D用于滑坡、泥石流等动态分析。DAN-3D基于连续介质力学原理将运动中的滑体等效为流变性质的流动体，通过设定其运动路径，流变模型及参数来模拟运动过程的三维数值模拟软件。

DAN-3D包括五种流变模型，分别是Frictional model (摩擦模型)、Plastic model (塑性模型)、Newtonian model (牛顿模型)、Bingham model (宾汉体模型)及Voellmy model，Stephen G. Evans、R. Sosio等多个学者通过不同模型实验模拟，表明Frictional model、Voellmy model能更好的模拟滑坡–碎屑流的运动和堆积过程[8] [15]。

3.2. 高程模型

新磨滑坡发生之后，许多地质学家都进行了现场调查和分析，本文采用滑坡前和滑坡后的地形图对滑坡运动过程进行分析，地形图是通过遥感影像数据获得，比例尺为1:10000。通过对地形图进行处理，得到了滑坡前和滑坡后的3D数字高程模型。

3.3. 参数选取及模型建立

本文选用Frictional model模拟滑坡滑源区，Voellmy model模拟流通区和堆积区。模拟过程中，上部边界是自由边界，下部边界是正应力和剪应力：

Frictional(摩擦)模型假设抗剪强度仅与正应力的关系：

$\tau =\sigma \left(1-r_u\right)\tan \phi$ (1)

式中：$\tau$ 为滑体的抗剪强度(Pa)；$\sigma$ 为正应力(Pa)；$r_u$ 为孔隙水压力系数；$\phi$ 为动摩擦角(˚)。

(1)中，孔隙水压力系数和动摩擦角之间存在如下关系：

${\phi }_b=\arctan \left(1-r_u\right)\tan \phi$ (2)

式中：${\phi }_b$ 为有效摩擦角(˚)。

Voellmy模型假设滑体受到的阻力为滑体受到的摩擦力和湍流流动造成的阻力之和：

$\tau =\sigma f+\rho g{v}^2/\xi$ (3)

式中：$f$ 是摩擦系数；$\xi$ 是湍流系数；$g$ 是重力加速度(m/s²)；$v$ 是滑体的平均速度(m/s)；$\rho$ 是滑体的密度(kg/m³)。

根据以上参数特点和新磨滑坡的真实调查情况，新磨滑坡选用Frictional(摩擦)模型和Voellmy模型复合模型进行模拟，经过多次模拟，反复调整，最终得到模型及参数如下(表1)：

Table 1. Models and parameters

表1. 模型及参数

	滑源区Ⅰ	流通区Ⅱ	堆积区Ⅲ
模型	Frictional模型	Voellmy模型	Voellmy模型
容重/kN∙m^−^3	28	28	28
摩擦角/(˚)	30	-	-
摩擦系数	-	0.10	0.15
湍流系数/m∙s^−^2	-	600	400
内摩擦角/(˚)	-	30	35
孔隙水压力系数$r_u$	0.10	-	-

3.4. 模拟结果及分析

新磨滑坡是典型的高速远程滑坡，根据掌握的滑坡资料，可将滑坡分为三个区，分别是滑源区、流通区和堆积区。通过DAN-3D模拟得到的不同时间新磨滑坡堆积体形态如图1所示：

Figure 1. Simulation results diagram of Xinmo Village

图1. 新磨村模拟结果图

3.4.1. 滑坡发生过程

一、模拟结果

根据上述滑坡–碎屑流堆积形态随时间的变化规律，将滑坡全过程分为三个阶段：

启动阶段：0~10秒，由图，可明显看出10秒时，滑体部分相对最初已发生明显位移变化。大约3.97 × 10⁶ m³的岩块从海拔高度为3400 m的位置高速剪出，沿NE-SW方向运动，滑体前缘高程为3000 m，此时滑体前缘和后缘相对高差为400 m。

刮铲流通阶段：10~40秒，这段时间内滑体高速运动，刮铲带走了流通区内大约5.48 × 10⁶ m³的松散堆积物，滑坡体积倍增，40秒时滑坡前缘已经到达大概高程2300 m的位置，而滑坡后缘才刚离开滑源区到达流通区，此时后缘和前缘位置相对高差达700 m。

流动堆积阶段：40~120秒，这个阶段又可分为两段，第一段是40~80秒，这段时间内滑坡堆积体形态变化差异大，从40秒时开始，滑坡呈扇形向前方扩散运动，50秒时滑坡前缘位置到达新磨村，60秒时滑坡完全覆盖新磨村，80秒时滑坡已经到达最前缘部分；第二段时间是80~120秒，此时滑坡基本已经停止运动，形态变化不大。在这个过程内滑体冲过了松坪沟河，然后在陡峭狭窄的河谷两侧停下来堆积，并形成了堰塞湖。

二、实际情况

新磨滑坡的发生大致时间是从5时39分到5时40分59秒，从5时39分开始，滑坡先发生了7秒钟的小规模崩塌，5时39分07秒，发生一次较大规模的滑动，这次滑动持续了5秒钟，此后到5时39分52秒，大规模的滑动引起了39秒的小规模崩塌。从5时39分53秒到5时40分22秒，这段时间记录了滑坡的主要刮铲运动过程，这个过程仅用时29秒，堆积过程持续了37秒。本文模拟滑坡发生总时间为120 s，主滑坡发生时间为80秒，实际滑坡发生总时间为119秒，主滑坡时间为66秒，模拟情况与实际发生基本一致，模拟情况良好。

3.4.2. 滑坡速度特征

一、模拟结果

DAN-3D模拟得到的各个时间段速度如图2所示：

Figure 2. Velocity diagram of the Xinmo landslide at different times

图2. 新磨滑坡不同时刻速度图

由模拟结果图可以看出，滑坡从高处滑出，10秒时最大速度达42 m/s，向下运动时，重力势能转化为动能，速度不断增大，最大速度达到了60 m/s，沿途不断刮铲，到达流通区的老滑坡堆积物时，碎屑流体势能转化为动能，推动松散堆积物向前运动，最终摧毁了新磨村，经过新磨村位置时最大速度可达48 m/s，之后速度开始不断减小，到110秒时，平均速度只有1 m/s，120秒时，速度基本为0，滑坡前缘已经到达最远位置。

3.4.3. 滑坡体积特征

一、模拟结果

滑坡初始滑动体积为3.97 × 10⁶ m³，启动阶段滑坡体积变化不大，刮铲流通阶段，体积迅速增大，滑坡离开流通区后，体积继续快速增大，到60秒时，体积基本达到稳定值，60秒到80秒，体积缓慢增大，80秒到120秒，体积几乎无变化，最终得到堆积物总体积为1326万m³。

滑源区在0时刻平均厚度为22.5 m，0~10 s，此时岩体崩解从高处滑出，滑体的厚度迅速减小，10~40 s，滑体不断刮铲流通区内的松散堆积物，厚度又逐渐增大，40 s时，滑坡–碎屑流呈扇形扩散，滑体的厚度不断减小，直到运动至堆积区的老滑坡堆积物上，滑坡厚度又开始缓慢增大最终滑坡的堆积物的平均厚度为13.2 m。

滑坡在流通区内刮铲物质最多，且中部偏动位置刮铲深度最深，形成的堆积物最厚，经过流通区，滑坡越过陡坎后，开始呈扇形扩散，在堆积区，滑坡主要堆积于东侧，这与实际情况相符合。

二、实际情况

滑源区最大厚度超过40 m，平均厚度为25 m，滑体总体积为5.47 × 10⁶ m³，最终大约有3.97 × 10⁶ m³的岩块滑动剪出，约150 × 10⁶ m³万的留在原处堆积；流通区内老滑坡堆积体体积为8.74 × 10⁶ m³，被铲刮物质估算体积约5.48 × 10⁶ m³，最终留下的堆积物体积为3.26 × 10⁶ m³；到达堆积区后，加上原有的老滑坡堆积物，形成的堆积体体积为1400万m³。滑体铲刮最终形成了U型沟槽，且刮铲后流通区形成了西高东低的地势。

模拟结果与实际情况对比发现(图1和图2)，体积相差不大，厚度堆积形态差别不大，但是模拟结果平均厚度偏高，可见DAN-3D模拟整体可靠，但存在误差。

4. 结论与展望

4.1. 结论

2017年6月24日的新磨滑坡属于一次重大地质灾害，此次灾害的发生对人民生命和财产安全都造成了重大伤害(10人死亡、73人失踪、3人受伤)，本文通过现场调查、资料收集等各种技术手段分析了该滑坡发生的工程地质条件、成因机制和运动过程。通过DAN-3D软件模拟了滑坡的发生过程，将模拟结果与滑坡真实调查情况做了对比分析。最终得出以下结论：

1) 新磨滑坡属于岩质滑坡，形成滑坡–碎屑流–堰塞湖的灾害型链。该滑坡分为滑源区、流通区和堆积区三个部分。

2) 本文在一些掌握的资料和数据的基础上，选用Frictional模型模拟滑坡滑源区，Voellmy模型模拟流通区和堆积区，成功建立了新磨村滑坡的模型，运用DAN-3D数值模拟软件进行了滑坡运动过程的模拟，最后将模拟结果与滑坡发生的真实情况做对比，发现二者基本一致，说明DAN-3D软件适用于此类滑坡的模拟。希望可以借助此软件对其它不稳定边坡进行模拟分析，最后可以采取一些针对性的防治措施，避免一些灾害的发生。

4.2. 展望

本文基于一些资料和数据对新磨高速远程滑坡成因机制进行了简单的分析研究，并提供了一些观点和DAN-3D滑坡模拟情况。但由于一些主观因素和客观因素，本文的研究还存在着许多不足之处：

1) 本文对滑坡的分区也很粗略，尤其是堆积区，有明显的岩石分层情况，但本文没有进行深入分析。

2) 本文的模拟也存在着不足，只采取了一种方法进行了模拟，模拟结果部分存在偏差，参数选取部分也没有校准，结果存在着偶然性。

最后，传统的基于实地调查的地质灾害评价方法很难对山区潜在的不稳定岩体进行早期识别，尤其是地震灾区，岩体受到的内部损伤很难察觉。所以建议使用一些现代技术，例如卫星、无人机的遥感数据，InSAR技术等，再通过地质专家的分析研究，希望可以成功地识别出潜在灾害并进行快速评估。