基于Voellmy模型的哀牢山地区滑坡泥石流灾害运动过程数值模拟

陈亮+宗士昌+李祥龙+陈红旗

摘要：我国西南山区滑坡泥石流灾害常发，滑坡泥石流灾害体运移堆积的范围是开展区域地质灾害危险性评价的重要指标之一。以云南省哀牢山地区水塘镇芭蕉树滑坡泥石流灾害为研究对象，利用Voellmy流体运动模型进行灾害体运移堆积全过程数值模拟反分析，获取数值模拟模型参数。利用反分析得到的模型参数建立与芭蕉树滑坡泥石流所处地质环境条件相似的戛洒镇南恩小学滑坡数值模型，进行滑坡泥石流灾害运移堆积过程模拟分析，模拟结果与野外地质灾害危险性评估结论符合程度较高。

关键词：滑坡；泥石流；Voellmy模型；数值模拟

中图分类号：P642 文献标识码：A 文章编号：1672-1683（2017）03-0113-07

Abstract：Landslide and debris flow are two common geological hazard types in the South-western mountainous area in China.One of the key factors for hazard assessment of landslides and debris flows is their post-initiation mobility.First，the runout process of Bajiaoshu landslide/debris flow was simulated with the Voellmy frictional model，and parameters were obtained.Then the simulation parameters were used to simulate the runout process of Nan′en School landslide for runout prediction and hazard assessment.The results proved to be consistent with field investigation results and earlier hazard assessment conclusions.

Key words：landslide；debris flow；Voellmy model；numerical simulation

滑坡泥石流災害是我国西南山区常见的地质灾害类型，其高速远距离运移堆积能力使其具有较强的破坏力和很大的影响范围。对其危险性的评价，除了需要计算评估滑坡体稳定性和泥石流发生频率之外，还需要评估滑坡泥石流发生后灾害体的运动速度和堆积范围。

21世纪以来，国内外对滑坡-碎屑流运动过程的研究主要体现在经验统计[1-4]、力学模型计算、物理模型试验及数值模拟[5-7]4个方面。A.Heim[8]首先提出滑坡架空等效摩擦因数ef 概念；A.E.Scheidegger等[9-10]在此基础上提出一系列基于体积-等效摩擦因数的滑距预测公式；王念秦等[11-12]也在我国大型滑坡研究的基础上提出了滑距预测经验公式。力学计算模型中，滑坡块体模型应用最为广泛，包括单一质点块体模型[13-14]、独立单元模型[15]以及可变块体运动二维模型[16]等；物理模型试验主要通过室内模型模拟碎屑流物质的滑动和扩散过程，分析碎屑流铲刮滑面、裹挟堆积体等机制，并为滑距计算提供试验参数[17-18]。由于滑坡在运动过程中属于非牛顿流体，在流动时不同高度的流体速度不同，因此层间会有摩擦所产生的剪应力[19-20]。浆体的流变特性是指其剪切变形时的剪应力与剪切率的关系。在描述这种非牛顿流变特性时，通常选取摩擦模型、Voellmy模型和Bingham模型[19-20]。

本文首先以云南省哀牢山地区水塘镇芭蕉树滑坡泥石流灾害为研究对象，利用Voellmy模型进行灾害体运移堆积全过程数值模拟反分析，测试该模型模拟泥石流运移堆积过程的适用性，并通过试错方法获取数值模拟模型参数；利用反分析得到的模型参数建立与芭蕉树滑坡泥石流所处地质环境条件相似的戛洒镇南恩小学滑坡数值模型，进行滑坡泥石流灾害运移堆积过程模拟分析，并与现场调查结果和危险性评估结论进行对比，进一步确定了Voellmy模型的适用范围和模拟参数选取方法。

1 Voellmy流体运动学模型原理

在连续介质模型中，泥石流流体被假设为一种非稳定以及非均质的流体，可以被2个主要流体参数来刻画，即流体高度H（x，y，t）（m）；平均流速U（x，y，t）（m/s）

2 水塘镇芭蕉树滑坡泥石流运动过程模拟反分析

2.1 滑坡基本特征

芭蕉树滑坡泥石流灾害2002年发生于云南省哀牢山地区水塘镇乡邦迈村原芭蕉树村民小组上方坡体。滑坡位于哀牢山中部，为构造剥蚀高中山地貌。滑坡发生在哀牢山群下亚组（Ptaa）地层之中，滑体岩性为黑云斜长片麻岩，片理层较薄，最薄片理厚度仅为2 cm，片理面可见白色的拉伸线理和眼球状构造，片理产状在50°～65°∠36°～45°之间，岩层顺倾坡外。片麻岩底部为一层厚度不足5 cm的白云母片岩夹层，抗剪强度和摩擦系数均较小，该夹层成为顺岩体片理面的易滑面，成为坡体失稳的优势面。作为红河深大断裂的一支，哀牢山断裂横穿滑坡体。滑坡源区位于哀牢山断裂下盘，滑坡堆积区则覆盖于哀牢山断裂上盘[21]。

该滑坡的平面形态呈“长椭圆形”（见图1），滑坡最大纵向长度为540 m，滑坡后缘陡坎至前缘剪出口处最大长度为126 m，地形呈圈椅状，滑体分为两部分，主堆积区平均宽度120 m，其前缘左侧为平均长度50 m，宽26 m的次级堆积区，整个滑体呈长椭圆形，中部较宽，两头较窄，滑坡面积为5.36万m2。滑坡体后缘高程1 170 m，前缘高程1 002 m，整体落差168 m。滑坡平均长度157 m，平均宽度120 m，整体平均厚度6～8 m，体积约13.2万m3。通过对滑坡结构特征的分析，可将滑坡分为5个区域（见图2）：滑源区、滑体飞行区、碎屑流堆积区、弹射体堆积区以及边缘气浪影响带。

图1 芭蕉树滑坡地质平面

Fig.1 Geological plan of Bajiaoshu landslide

图2 芭蕉树滑坡地质剖面

Fig.2 Bajiaoshu Landslide profile

通过对滑坡附近居民的走访及查问，滑坡所处山体斜坡近年来一直比较稳定，并未发现贯通裂缝及局部剥、坠落等斜坡变形迹象。但2002年8月，滑坡所处区域发生了一次罕见的极端降雨过程，从8月13日开始，滑坡所处的山坡及坡脚出现洪水，洪水深度接近成年人膝盖，山坡表面的残坡积物等松散堆积层被洪水冲刷带走形成较小规模的坡面泥石流。13日上午8时30分，该滑坡开始启动，伴随着轰鸣的爆裂声，滑坡体整体从前缘剪出口冲出，向前下方飞行了50 m左右后坠地解体，滑坡主体分解成为碎屑流继续向前下方滑动364 m后停止运动；滑坡左侧一小部分滑体坠地后发生弹射，再次飞行170 m后坠地分解；滑坡主体分解成为碎屑流运动中，其前缘和侧缘产生近30秒的强大气浪，这股气浪将附近的芭蕉树小组摧毁。

2.2 运动过程数值模拟反分析

2.2.1 模型建立与模拟参数取值

数字高程模型是滑坡泥石流数值模拟最基础、最重要的步骤，这直接决定了模拟结果的精度以及正确性。通过对工作区地图进行矢量化后，利用ArcGIS软件，建立TIN（Triangulated Irregular Network）最优化不规则三角网模型，通过对TIN模型数字高程点的采集，获取模型网格点数字高程，并进行内插计算，最终获得研究区DEM。通过GIS模拟系统内置分析代码，可以得到流体的流速高度过流量等关键参数，并将计算所得结果以GIS格式数据的形式进行三维效果展示。在前述连续介质模型原理基础上，利用动力数值模拟计算软件，对该滑坡泥石流进行实例研究。

模型长度=752 m，宽度=650 m，共计48 995个节点，48 552个单元根据现场地质灾害调查结果，估算滑坡初始固体物质总量为121 469 m3，建立的三维模型见图3，为滑坡滑动前平面位置，图中彩色部分为滑坡体，灰色为地表。

滑坡体的勘察作业中进行了对滑体物质的钻探取样和室内土体剪切试验，试验结果见表1，包括土样的快剪强度和残余强度。由于此模拟关注于当滑坡启动之后发生运移和堆积的过程，在模拟中将滑坡体定义为已经产生滑动破坏的黏性岩土体，因此模拟反分析的参数取值区间以残余剪切强度为依据，而滑体密度采用了土工试验值，滑体厚度采用工程地质剖面图中滑体的平均厚度10 m。地表参数表示滑坡泥石流物质的运动路径上的地表的摩阻力，由于地表的复杂和不均匀性，无法直接用试验获取该参数，因此作为变量，利用迭代反分析获得。滑坡初始速度设置为零。经过多次反分析模拟，通过对滑坡泥石流模拟运动过程和堆积形态与实际滑坡泥石流堆积形态进行比对，最终芭蕉树滑坡泥石流模拟参数取值见表1。

2.2.2 模拟反分析结果

数值模拟反分析结果见图4、图5，其中图4为芭蕉树滑坡泥石流运动堆积模拟过程。从图5中可以看出，模拟的滑坡经历了沿滑动面下滑、沿平缓地表运移和完成堆积的过程，在整个运动过程中滑坡体产生了很大的变形，不再保持滑动前的形态，岩土体流体化，形成了沿一条較窄的路径堆积的形态。

图5为滑坡沿运移堆积方向的堆积体纵剖面图，图中绿色曲线为地形剖面线，红色曲线为滑坡泥石流堆积体表面剖面线，其纵坐标读取右侧坐标轴高程值；灰色区域为滑坡堆积体剖面，其纵坐标读取左侧坐标轴堆积厚度值。将图5结果与图2实测堆积剖面做对比，可以看出滑坡物质在下滑后形成的运移堆积体堆积形态与图中实际现场调查得出的形态吻合度较高，可以看出数值模拟得到的堆积剖面基本吻合，说明模拟建模和参数取值是合理的。

图4 芭蕉树滑坡泥石流运动堆积模拟过程

Fig.4 Flow simulation process of Bajiaoshu landslide and debris flow

图5 芭蕉树滑坡堆积体纵剖面

Fig.5 Longitudinal profile of Bajiaoshu landslide

3 南恩小学滑坡泥石流运动过程模拟

3.1 滑坡基本特征

该滑坡体位于云南省哀牢山地区戛洒镇南恩河东南侧的斜坡体之上，为一由片麻岩大块石与坡积土混合组成的崩塌堆积体滑坡，见图6。滑坡后缘发育在老滑坡堆积体中，即老滑坡堆积体前缘过渡部位，呈半圆弧状，覆盖层由第四系含碎石粉质黏土和层内夹的崩塌片麻岩孤石组成，下伏基岩为片麻岩，无基岩出露，地面标高1 906～1 908 m，坡度约40°，地形较陡。滑坡主滑方向为60°，后缘高程1 770 m，前缘高程1 604 m，高差166 m，纵向长571 m，横向宽160～227 m，面积为6.5万m2，滑体厚度一般18～26 m，平均厚度约22 m，规模约143.0万m3，为大型中层牵引式土层滑坡。

滑坡体的发展分为两个阶段，第一阶段是20世纪90年代初S307省道的修建和1999年7月南恩小学的修建，开挖坡脚减小了坡体的阻滑段，改变了坡体内部的应力状态，随着坡体内的应力调整，堆积体蠕滑变形逐渐加大；第二阶段是暴雨触发阶段，2002年8月11日-14日的暴雨是这次滑坡形成的最直接的诱因，据调查数据显示，2002年8月11日-14日，境内发生强降雨，8月11日至12日，降雨量为4 mm，13日升至27 mm，到14日降雨量猛升至101 mm，在连续降雨，特别是14日的强降雨情况下，大量的雨水下渗，导致滑坡区土体饱和，滑体重量加大，基覆界面处土体的抗剪强度急剧降低，滑坡从前缘开始逐渐滑动，南恩小学操场出现张裂缝、教室开裂、挡墙拉裂及S307省道挡墙拉裂，后缘的居民房屋墙体位移明显。斜坡体前部的强烈变形导致老滑坡堆积体后缘出现下错，并在随后的自身调整中达到了自我平衡，目前处于基本稳定状态。但考虑到滑坡体在强降雨条件下可能会产生滑动并演化为泥石流，对下游南恩河沿岸多个村寨存在很大威胁，因此需要进行滑坡泥石流启动后运动堆积距离和破坏范围进行模拟预测，以进行危险性评价。

3.2 運动过程数值模拟

利用南恩小学滑坡勘察资料，建立三维地形模型和滑坡体初始位置模型见图7。由于南恩小学滑坡所处位置地质环境条件与芭蕉树滑坡泥石流基本一致，因此直接应用芭蕉树滑坡泥石流数值模拟反分析得到的模拟参数（表1），进行南恩小学滑坡运动过程模拟。

假设滑坡体已启动，初始速度为零。数值模拟得到的滑坡泥石流运移堆积模拟过程见图8。从图8南恩小学滑坡运移堆积模拟过程平面图中可以看出，滑坡运移堆积产生了分叉，分别沿北侧和南侧两条山谷进行由西向东的运动。图9为勘察结果与模拟结果对比图，底图为勘察所评估的南恩小学滑坡泥石流运动路径，底图西南角覆盖层为数值模拟滑坡泥石流堆积范围，可以看出模拟结果与评估结果符合程度很高。图10和图11分别为北侧和南侧山谷内部滑坡泥石流堆积体剖面图。对比模拟结果和地质灾害调查结果，可以发现模拟中滑坡泥石流的运移路径存在分叉，分为南北两条，北部堆积体又存在二级分叉。这一模拟结果体现出了以GIS为平台运用Voellmy模型的优势，即以数字栅格为模拟单元，进行单元内质量守恒和整体动量守恒的数值计算，不受传统有限体积法模拟中滑坡模型的网格剖分限制，可以实现滑坡体的离散性模拟。

4 结论

（1）Voellmy流体运动模型遵守动量守恒定律，可以较好得模拟滑坡泥石流启动后的运移和堆积过程，从而得到滑坡泥石流灾害的破坏范围和破坏强度。

（2）在相似的地质环境条件下，对同一类型的地质灾害进行破坏运动过程数值模拟所需的计算参数可以取同一组数值。

（3）Voellmy模型参数的取值会显著影响模拟结果，由于模型参数不易通过实验获取，因此合理的方法是选取所需模拟的滑坡泥石流所在区域内已发生的滑坡泥石流灾害进行数值模拟反分析，通过迭代法得到Voellmy模型参数取值，作为模拟参数。

参考文献（References）：

[1] 〖ZK（#]Prochaska A B，Santi P M，Higgins J D，et al.Debris-flow runout predictions based on the average channel slope （ACS）.Eng Geol[J].Engineering Geology，2008，98（1）：29-40.DOI：10.1016/j.enggeo.2008.01.011

[2] Berti M，Simoni A.DFLOWZ：A free program to evaluate the area potentially inundated by a debris flow[J].Computers & Geosciences，2014，67（2）：14-23.DOI：10.1016/j.cageo.2014.02.002

[3] Gartner J E，Cannon S H，Santi P M，et al.Empirical models to predict the volumes of debris flows generated by recently burned basins in the western U.S.[J].Geomorphology，2008，96（s 3-4）：339-354.DOI：10.1016/j.geomorph.2007.02.033

[4] Dahl M P J，Mortensen L E，Jensen N H，et al.Magnitude-frequency characteristics and preparatory factors for spatial debris-slide distribution in the northern Faroe Islands[J].Geomorphology，2013，188（436）：3-11.

[5] Chen H X，Zhang S，Peng M，et al.A physically-based multi-hazard risk assessment platform for regional rainfall-induced slope failures and debris flows[J].Engineering Geology，2015.DOI：10.1016/j.enggeo.2015.12.009

[6] Calvo B，Savi F.A real-world application of Monte Carlo procedure for debris flow risk assessment[J].Computers & Geosciences，2009，35（5）：967-977.DOI：10.1016/j.cageo.2008.04.002

[7] Aronica G T，Biondi G，Brigandì G，et al.Assessment and mapping of debris-flow risk in a small catchment in eastern Sicily through integrated numerical simulations and GIS[J].Physics & Chemistry of the Earth Parts A/b/c，2012，49（211）：52-63.DOI：10.1016/j.pce.2012.04.002

[8] HEIM A.Landslide and human lives [M].Vancouver：Bitech Press，1989：57-79.

[9] SCHEIDEGGER A E.On the prediction of the reach and velocity of catastrophic landslides[J].Rock Mechanics，1973，5（4）：231-236.DOI：10.1016/0148-9062（74）91709-4

[10] COROMINAS J.The angle of reach as a obility index for small and large landslide[J].Canadian Geotechnical Journal，1996，33（2）：260-271.

[11] 王念秦，张倬元.一种典型黄土滑坡的滑距预测方法[J].西北大学学报，2003，33（1）：111-114.（WANG Nian-qin ，ZHANG Zhuo-yuan.A forecasting method of sliding distance on typical loess landslides[J].Journal of Northw est University，2003，33（1） ：111-114.（in Chinese）） DOI：10.3321/j.issn：1000-274X.2003.01.030

[12] 王思敬，王效宁.大型高速滑坡的能量分析及其灾害预测[C].滑坡论文集.成都：四川科技出版社，1989：117-124.（WANG Si-jing，WANG Xiao-ning.Energy analysis and disaster prediction of large - scale high - speed landslide[C].Landslide Proceedings.Chendu：Sichuan Science and Technology Press，1989：117-124.（in Chinese））

[13] FANNIN R J，WISE M P.A method for calculation of debris flow travel distance[C]// Proceedings of the 48th Canadian Geotechnical Conference.[S.l.]：[s.n.]，1996：643-650.

[14] MCLELLAN P J，KAISER P K.Application of two-parameter model to rock avalanches of the Mackenzie mountains[C]// Proceedings of the 4th International Symposium on Landslide.[S.l.]：[s.n.]，1984：135-140.

[15] CUNDALL P A，STRACK O D L.A discrete numerical model for granular assemblies[J].Geotechnique，1979，29（1）：47-65.

[16] 劉忠玉，马崇武，苗天德，等.高速滑坡远程预测的块体运动模型[J].岩石力学与工程学报，2000，19（6）：742-746.（LIU Zhong-yu，MA Chong-wu，MIAO Tian-de，et al.Kinematic block model of long run-out prediction for high-speed landslides[J].Chinese Journal of Rock Mechanics and Engineering，2000，19（6）：742-746.（in Chinese）） DOI：10.3321/j.issn：1000-6915.2000.06.012

[17] EGASHIRA S，HONDAB N，ITOHC T.Experimental study on the entrainment of bed material into debris flow[C]// Physics and Chemistry of the Earth，Part C：Solar，Terrestrial and Planetary Science.[S.l.]：[s.n.]，2001：645-650.DOI：10.1016/S1464-1917（01）00062-9

[18] MANZELLA I，LABIOUSE V.Qualitative analysis of rock avalanches propagation by means of physical modelling of non-constrained gravel flows[J].Rock Mechanics and Rock Engineering，2008，41（1）：133-151.DOI：10.1007/s00603-007-0134-y

[19] 马宗源，廖红建，张骏.Bingham型黏性泥石流流体的三维数值模拟[J].西安交通大学学报，2008，42（9）：1146-1150.（MA Zong-yuan，LIAO Hong-jian，ZHANG Jun.Three dimensional numerical simulation of bingham viscous debris flow fluid[J].Journol of Xi′an jiaotong university，2008，42（9）：1146-1150.（ in Chinese）） DOI：10.3321/j.issn：0253-987X.2008.09.018

[20] 齐超，邢爱国，殷跃平，等.东河口高速远程滑坡-碎屑流全程动力特性模拟[J].工程地质学报，2012，20（3）：334-339.（QI Chao，XING Ai-guo，YIN Yue-ping，et al.Numerical simulation of dynamic behavior of Donghekou rockslide-debris avalanche[J].Journal of Engineering Geology，2012，20（3）：334-339.（ in Chinese）） DOI：10.3969/ j.issn.1004-9665.2012.03.005

[21] 中国地质环境监测院.云南新平县地质灾害详细调查成果报告[R].2007，106-111.（ China Institute of Geo-Environment Monitoring.Report on the findings of geological hazard in Xinping county，Yunnan province[R].2007，106-111.（in Chinese））