基于Flow-R和FLO-2D的沟谷型泥石流危险性评价

林文 周伟 李婧 杨莉 唐泽

摘要：泥石流是中国西南山区频发的地质灾害之一，具有突发性、流速快以及破坏力强等特点，往往造成巨大的人员伤亡和经济损失。地质灾害危险性评估是防灾减灾管理和防治环节中的有效措施之一，用于合理量化泥石流运动堆积线路沿程泥石流灾害危险性的空间分布特征。以四川省都江堰市龙池镇麻柳沟为研究对象，选取坡度、平面曲率和上坡汇水面积作为沟谷型泥石流启动区识别因子，基于Flow-R模型识别出的泥石流启动点，耦合FLO-2D模型模拟麻柳沟泥石流运动堆积过程，将泥石流堆积深度、运动速度和暴发频率相结合分析泥石流的危险性，并采用混淆矩阵对模拟结果进行了精度评估。结果表明：基于Flow-R和FLO-2D耦合模型的麻柳沟泥石流的运动堆积模拟结果与实际情况吻合度较高。模拟泥石流堆积区面积2.19万m2，堆积体积4.92万m3，泥石流危险性由沟道中心向两侧逐渐降低，高危险区占比42.3%，中危险区占比23.1%，低危险区占比34.6%。麻柳沟形成区、流通区沟道中散落堆积物体积大约为30万m3，麻柳沟泥石流仍处于高度活跃阶段，危险性较高，今后在降雨和水动力充分的条件下大量松散物质可能将被重新激活，一旦发生泥石流将有可能造成重大地质灾害。

关 键 词：泥石流; 危险性评价; 泥石流启动点; Flow-R; FLO-2D; 都江堰市

中图法分类号： P642

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2022.05.023

0 引 言

泥石流是中国西南山区频发的地质灾害之一，具有突发性、流速快以及破坏力强等特点，往往会造成巨大的人员伤亡和经济损失。泥石流危险性评价是灾害学重点研究的课题之一。泥石流活动规模是泥石流危险性评价的一个重要指标[1]，它取决于控制流量触发的运动特征、泥石流沿程侵蚀和泥石流峰值流量等因素。FLO-2D是一种被广泛用于模拟洪水的软件，近年来也被用来模拟沟谷型泥石流的运动和堆积，便于对沟谷型泥石流进行危险性评价。阮德修等[2]将FLO-2D与3DMine进行耦合，模拟了尾矿库溃坝时的泥石流流动过程，并预测了尾矿库溃坝和半溃坝情况下的致灾程度。Chang等[3]通過实地调查，分析了汶川地震灾区泥石流的形成条件，运用FLO-2D模型，模拟了泥石流运动及其危害范围。张建石[4]通过统计分析，得出震后多年间的物源演变规律及特征，并利用FLO-2D模型模拟在不同降水频率下肖家沟泥石流的冲出特征。Nocentini[5]将FLO-2D模型与DAN-W模型进行耦合，预测了研究区潜在泥石流的危险性。黄靖玲等[6]将HEC-HMS模型与FLO-2D模型耦合，对马贵河流域进行了危险性评估。如果要获得较好的FLO-2D模型模拟效果，关键在于能找出相对准确的泥石流启动点。这就需要结合实地情况和专业经验，还应考虑到一定的客观因素。为了解决泥石流启动点识别的问题，本文引入了一种识别沟谷型泥石流启动点的软件Flow-R。该软件是一种基于GIS的重力灾害评价模型，可高效快速识别沟谷型泥石流启动点[7]。

四川省龙池镇龙溪河流域自“5·12”地震后，流域内松散物源量增加，地质灾害活动频繁，大量研究表明，泥石流灾害已成为汶川地震重灾区龙池镇最大的地质安全隐患[8]。2010年8月13日，流域在强降水条件下暴发泥石流灾害，造成巨大损失。本文以麻柳沟为研究对象，将FLO-2D模型与Flow-R模型耦合使用，来识别泥石流大多启动点，模拟泥石流淹没区范围，计算泥石流运动的最大动量和最大堆积深度[9];同时，结合流速、泥深与降水暴发频率来探究麻柳沟的危险性。研究成果可为麻柳沟的泥石流灾害的危险分析、预警及发展规划提供科学的参考依据。

1 研究区概况

麻柳沟位于四川省都江堰市龙池镇，沟口距离龙池镇约0.65 km。沟口地理坐标为103°33′20″E，31°03′55.11″N（见图1）。麻柳沟流域面积约为0.93 km2，流域相对高差为913 m，主沟长度为2.3 km，平均纵坡降401‰。巨大的相对高差对流域内的松散堆积物而言具有足够的势能条件，使得该沟易形成危害极大的高动能泥石流。

研究区地处龙门山断裂带。2008年5月12日汶川县映秀镇发生8.0级大地震，其主发震断裂映秀-北川断裂横贯研究区，导致研究区内岩层破碎严重，成为典型的构造不稳定区域。震后流域内发生大量崩塌滑坡，为暴雨泥石流提供了充足的松散物质来源。

研究区属于中亚热带湿润季风气候区，日照时间较少，降水充沛，阴雨天气频繁。多年平均降水量为1 134.8 mm（1987～2008年），5～9月降水量占全年降水量的80%。常年充分的降水使坡面松散堆积物直接运移到沟道堆积，形成泥石流沟道物源。2010年8月13日龙池镇发生强降水事件，最大1 h降水量达75.0 mm，连续3 h降水量达到150.0 mm，导致泥石流暴发。

2 研究方法

Flow-R模型将坡度、平面曲率、上坡汇水面积作为沟谷型泥石流启动点的识别因子，然后将其转化为Flow-R模型能识别的ASCII文件，通过对不同的因子设置阈值，识别出研究区的启动点。

基于DEM数据转化为FLO-2D能识别的ASCII文件，结合现场实地调查进行网格划分。调查发现研究区泥石流沟道流动宽度大于10 m，所以将本次数值模拟网格划分为10 m的网格满足本次研究的需要，对每个网格赋值唯一的高程值、粗糙系数等基础地形因子，将Flow-R识别的启动点导入网格，输入启动点泥石流流量及流变参数，模拟麻柳沟“8·13”特大泥石流运动堆积过程，构建混淆矩阵对运动堆积范围进行精度验证。

危险性评价建模是将耦合模型结果获取的淹没区流速、泥深和降水频率相结合，将泥石流沟划分为高危险区、中危险区和低危险区。工作流程如图2所示。

2.1 泥石流启动点识别

坡度（地形条件）、松散物源可用性（物源条件）、输水量（水源条件）是影响泥石流形成的关键因素，当这些关键因素满足一定条件时即可诱发泥石流[10]。Flow-R模型考虑了泥石流形成的三大影响因素，通过设置各个因子泥石流启动识别阈值（详见3.1节），将同时满足阈值的栅格单元划分为“启动区”，将至少有一次不满足阈值的栅格单元划分为“不是启动区”。

2.2 泥石流运动流通堆积模拟

FLO-2D模型用于模拟泥石流的运动过程和跳动特征[11]。模型以网格为单元计算泥石流流动，在软件划分出大小相同的规则网格后，将高程、粗糙系数等数值赋予给各网格。二维流动是通过运动方程的数值积分和流体体积守恒来实现的，即网格中固体物质和水的比值不断变化，而总体积变化，且在模拟过程中，当低浓度流体进入网格遇到已停淤的高浓度流体时，汇流处流体发生混合并继续流动，控制方程如下。

（1） 连续方程。

ht+（uh）x+（vh）y=l（1）

式中：l表示水力坡降;h表示泥深;v表示泥石流在y方向的速度;u表示泥石流在x方向的速度。

（2） 运动方程。

Sfx=Scx-hx-VxgVxx-VygVyy-1gVxt（2）

Sfy=Scy-hy-VygVyy-VygVyx-1gVyt（3）

式中：Sfx和Sfy分别表示x、y方向的摩擦坡降;Scx和Scy分别表示x、y方向的底床坡降。

（3） 流变方程。

τ=τc+τmc+τt+τv+τd（4）

式中：τc表示黏性屈服应力;τv表示黏性剪应力;τmc表示摩尔-库伦剪应力;τt表示紊流剪应力;τd表示扩散剪应力。

2.3 泥石流危险性评价

对单沟泥石流危险性评价，根据唐川等的研究，结合泥石流运动速度和泥石流深度对泥石流强度进行划分[12]（见表1）。在此基础上，Chang等[3]以泥石流强度与暴发频率相结合的方式进行泥石流危险性评价，将研究区危险性分为高、中、低3个等级（见图3）。本文将耦合模型模拟结果结合泥石流运动深度、速度以及降水频率，对研究区的泥石流危险性进行评价。

3 結果分析

3.1 泥石流启动点识别

（1） 坡度。

坡度是泥石流形成的主要因素。大多数地区发生泥石流的地形坡度大于15°[6]，Ortigao等[13]认为，泥石流启动时的地形坡度介于20°～25°之间。麻柳沟坡度主要分布在20°～50°，占整个流域的86.41%。结合麻柳沟地形坡度特征，将泥石流启动的地形坡度阈值定为20°。

（2） 平面曲率。

麻柳沟泥石流是典型的沟谷型泥石流，泥石流的物源一般分布在曲率为凹形的区域。平面曲率表示地形的凹凸变化，可用于识别泥石流沟道等负地形，间接体现泥石流的物源。现有研究表明，泥石流启动时的平面曲率阈值大多在-0.5/100 m-1到-2/100 m-1之间[14]。Delmonaco等[15]在平面曲率为-1/100 m-1的凹地形中识别出泥石流物源。麻柳沟负地形平面曲率统计结果表明，该沟平面曲率大于-3/100 m-1的比例达到84%。本文将研究区平面曲率阈值定为-3/100 m-1。

（3） 汇水面积。

上坡汇水面积，即汇流累积量，表示DEM中流经每个栅格的流水累积量，旨在识别研究区内任何活跃的河流、沟渠或暗流。Park等通过观测实际发生的泥石流物源，构建了泥石流源区坡度与汇水面积的阈值预测模型[16]。该模型如下：

θ=15e32.34A-160.85（5）

式中：A表示汇水面积，θ表示坡度。参考公式，通过设定100，500 m2和1 000 m2汇水面积阈值提取水系，结果表明1 000 m2阈值最符合实际情况。对麻柳沟流域进行了汇水面积统计，结果显示，该沟78%的汇水面积在0～1 000 m2范围内，因此，将研究区汇水面积率阈值定为1 000 m2。

由图5启动点模拟结果可知：麻柳沟泥石流启动区主要分布在沟道及沟道两侧，启动点坡度分布在20°～45°之间，其中，超过一半的启动点分布于30°～40°之间，启动点高程主要分布在1 100～1 300 m。沟顶细小汇水沟道为泥石流提供了丰富的水源，促使下部沟道沉积物转化为泥石流。

3.2 泥石流运动流通堆积模拟

参考水文地质手册，结合实地考察，泥石流启动点主要分布在杨家沟和矶子沟。采用雨洪法分别计算了杨家沟和矶子沟处启动点的洪峰流量。泥石流在运动过程中挟带有大量的松散物质，参考FLO-2D使用手册[11]，泥石流中水和沉积物总混合物的体积可以通过将水体积乘以放大系数（BF）来确定。放大系数 BF 与泥石流的体积浓度密切相关[10]，采用野外配浆实验获取体积浓度为0.58，计算出麻柳沟启动点流量。对洪峰流量采用简单概化的五边形方法来求取泥石流峰值流量过程线[17]。

泥石流的屈服强度τy及黏滞参数η与堆积范围相关[18]。τy、n和η可由下式计算：

n=0.33C-0.15VeCv-0.15lnh（6）

η=α1eβ1Cv（7）

τy=α2eβ2Cv（8）

式中：α1、α2、β1和β2为经验系数;Cv为泥石流体积浓度;h为泥深;n为曼宁系数。曼宁系数常被用来确定水和泥石流固体物质之间的接触关系。

对麻柳沟进行现场调查发现，沟道内植被稀疏，而且固体堆积物平均厚度在10 cm以上。曼宁系数可以间接表现流体的流变特征，根据公式（6）计算出的麻柳沟曼宁系数为0.12。层流阻滞系数K、α和β均为经验系数，一般通过查阅 FLO-2D 使用手册或进行专项试验获取。本次研究，根据前人对经验系数的取值[19]以及FLO-2D使用手册[11]中的参数值列表，获得了麻柳沟FLO-2D模拟的主要参数（见表3）。

麻柳沟数值模拟结果如图6所示。由图6（a）可知：沟口的堆积面积为2.19×104 m2，占总流域面积的2.3%，最大堆积厚度为4.5 m，平均堆积深度为2.25 m，冲出总量为4.92×104 m3。研究区最大堆积深度达5.90 m，最小堆积深度为0.04 m，沟道中心泥深最大，向沟道两侧泥深减少，沟道中和堆积扇前段泥深较大，速度较快，沟道两侧和堆积扇翼端堆积深度较小，运动速度较慢。

利用混淆矩阵对模拟结果进行精度评价，并采用准确率（A）和敏感度（S）指标对模拟结果进行评价，准确率和敏感度越高，表明模拟结果越好[20]。

A=TP+TNN（9）

S=TPTP+FN（10）

式中：TP为实际淹没区与模拟淹没区均覆盖的栅格数;TN为实际淹没区与模拟淹没区均未覆盖的栅格数;FN为实际淹没区覆盖而模拟淹没区未覆盖的栅格数;N为研究区栅格总数。根据麻柳沟2011年4月26日谷歌影像提取麻柳沟泥石流淹没区范围（见图6（b）），根据公式（9）和公式（10）计算出的模拟结果准确率为91.0%，敏感度为72.3%。表明耦合模型应用于麻柳沟泥石流的运动堆积模拟的结果与实际情况吻合度较高。

3.3 泥石流危险性评价

模拟还原了麻柳沟“8·13”泥石流运动堆积过程。由麻柳沟泥石流危险性分区图（见图7）可知：高危险区主要分布在沟道中心，面积约为58 300 m2，占危险区总面积的42.3%。由沟道中心向两侧泥石流危险性逐渐降低，中危险区面积为31 200 m2，占总面积的23.1%;低危险区面积为45 400 m2，占总面积的34.6%。高危险区域分布范围较大，低危险和中危险分布范围相对较小。沟道两侧和堆积扇翼端主要为低危险区域;泥石流形成区和堆积区边缘主要为中危险区;高危险区则主要集中在沟道中和堆积扇前段。在形成区和流通区的松散堆积体总量约30万m3，在极端降水天气下极易发生次生泥石流灾害，应该加强防灾减灾工作和泥石流风险管理工作。

4 结 论

本文以麻柳沟为研究对象，基于Flow-R模型识别泥石流启动点后耦合FLO-2D模型，用以模拟泥石流运动;采用速度、泥石流堆积深度和降水暴发频率相结合的方法，对单沟泥石流危险性进行了评价，可以得出如下结论。

（1） 麻柳沟泥石流启动区绝大部分布在沟道及沟道两侧，启动点坡度分布在20°～45°之间，其中超过一半的启动点分布于30°～40°之间，高程主要分布在1 100～1 300 m范围内。

（2） 模拟泥石流堆积区面积为2.19万m2，堆积体积为4.92万m3。计算形成区、流通区沟道中散落堆积物体积大约为30万m3，耦合模型应用于麻柳沟模拟的准确率为91.0%，敏感度为72.3%，模拟精度较高，可以为土地资源规划、灾害防治设计提供参考依据。

（3） 麻柳沟高危险区分布范围较大，占危险区面积的42.3%，沟道中心危险性较高，沿沟道中心向两侧危险性降低，高危险区主要集中在沟道中和堆积扇前沿，低危险区主要分布于沟道两侧和堆积扇翼端，泥石流上沟段和堆积区边缘主要为中危险区。

综上所述，基于Flow-R和FLO-2D耦合模型将泥石流源区和泥石流运动相结合，可以得出泥石流運动堆积线路沿程灾害危险性空间分布特征。根据危险分布特征图可知，麻柳沟泥石流仍处于高度活跃状态，今后在降水和水动力充分的条件下，大量的松散物质可能被重新激活，一旦发生泥石流将可造成重大的地质灾害。本文研究成果可为泥石流灾害监测预警提供建设性考虑。本文研究未考虑物源量、岩性等其他因子对泥石流启动点识别的影响，以及未考虑泥石流在运动过程中侧蚀和下蚀的影响，使得数值模拟的结果稍显偏小，且较为保守，在今后的研究中需要对其进行补充和完善。

参考文献：

[1] LUNA Q B，REMATRE A，VAN ASCH T W J，et al.Analysis of debris flow behavior with a one dimensional run-out model incorporating entrainment[J].Engineering Geology，2012，128：63-75.

[2] 阮德修，胡建华，周科平，等.基于FLO-2D与3DMine耦合的尾矿库溃坝灾害模拟[J].中国安全科学学报，2012，22（8）：150-156.

[3] CHANG M，TANG C，VAN ASCH T W J，et al.Hazard assessment of debris flows inthe Wenchuan earthquake-stricken area，southwest China[J].Landslides，2017，14（5）：1783-1792.

[4] 张建石.汶川县肖家沟泥石流物源演变及冲出规模研究[J].人民长江，2020，51（8）：37-43.

[5] NOCENTINI M，TOFANI V，GIGLI G，et al.Modeling debris flows in volcanic terrains for hazard mapping：the case study of Ischia Island（Italy）[J].Landslides，2015，12（5）：831-846.

[6] 黄靖玲，徐粒，周华真，等.基于历史灾害情景再现的小流域山洪风险研究[J].灾害学，2020，35（3）：194-200.

[7] HORTON P，JABOYEDOFF M，RUDAZ B，et al.Flow-R，a model for susceptibility mapping ofdebris flows and other gravitational hazards at a regionalscale[J].Natural Hazards and Earth System Sciences，2013，13（4）：869-885.

[8] 許强.四川省8.13特大泥石流灾害特点、成因与启示[J].工程地质学报，2010，18（5）：596-608.

[9] CHEN H，DADSON S，CHI Y.Recent rainfall-induced landslides and debris flows in northern Taiwan[J].Geomorphology，2018，77：112-125.

[10] RICKENMANN D，ZIMMERMANN M.The 1987 debris flows in Switzerland：documentation and analysis[J].Geomorphology，1993，8（2-3）：175-189.

[11] FLO-2D Software Inc.FLO-2D，2-Dimensional Flood Routine Model Manual[Z].State of Arizona：FLD-2D Software，Inc.，2004.

[12] 唐川，周钜乾，朱静，等.泥石流堆积扇危险度分区评价的数值模拟研究[J].灾害学，1994，9（4）：7-13.

[13] ORTIGAO J A R，KANJI M A.Landslide classification and risk management[M].Berlin：Springer，2004.

[14] BLAHUT J，WESTEN C，STERLACCHINI S.Analysis of landslide inventories for accurate prediction of debris-flow source areas[J].Geomorphology，2010，119（1-2）：36-51

[15] DELMONACO G，LEONI G，MARGOTTINI C，et al.Large scale debris-flow hazard assessment：a geotechnical approach and GIS modeling[J].Natural Hazards and Earth System Sciences，2003，3：443-455.

[16] PARK D W，LEE SR，VASU N N，et al.Coupled model for simulation of landslides and debris flows at local scale[J].Natural Hazards，2016，81：1653-1682.

[17] 夏添.震区泥石流危险性评价及预警减灾系统研究[D].成都：成都理工大学，2013.

[18] O’BRIEN J S，JULIEN P，FULLERTON W.Two-dimensional water flood andmudflow simulation[J].Journal of Hydraulic Engineering，1993，119：244-261.

[19] 程思.都江堰市龙溪河流域震后多沟同发泥石流危险性及易损性研究[D].成都：成都理工大学，2015.

[20] CARRANZA E，CASTRO O T.Predicting lahar-inundation zones：case study in West Mount Pinatubo，Philippines[J].Natural Hazards，2006，37：331-372.

（编辑：赵秋云）

Gully debris flow hazard assessment based on Flow-R and FLO-2D coupling models

LIN Wen，ZHOU Wei，LI Jing，YANG Li，TANG Ze

（State Key Laboratory of Geohazard Prevention and Geoenvironment Protection，Chengdu University of Technology，Chengdu 610059，China）

Abstract：

Debris flow is one of the frequent geological disasters in southwest  mountainous areas of China.It is characterized by sudden occurrence，fast velocity and strong destructive force，which often causes huge casualties and economic losses.Risk assessment of geological hazards is one of the effective measures in disaster prevention and mitigation management and prevention，which can reasonably quantify the spatial distribution characteristics of debris flow hazard along the debris flow movement and accumulation route.Taking Maliu gully in Longchi Town，Dujiangyan City，Sichuan Province as the research object，slope，plane curvature and upslope catchment area were selected as the identification factors for the initiation area of gully debris flow.Based on the starting point of debris flow identified by Flow-R model，the FLO-2D model was coupled to simulate the movement and accumulation process of debris flow in Maliu gully.The debris flow hazard can be determined by combination of accumulation depth，velocity and outbreak frequency of debris flow.In addition，the accuracy of the simulation results was assessed by confusion matrix.The results showed that the simulation results of Maliu gully debris flow by Flow-R and FLO-2D coupling models was in good agreement with the actual conditions.The simulated debris flow accumulation area is 2.19×104m2，and the accumulation volume is 4.92×104m3.The debris flow risk gradually decreases from the center of the channel to both sides.High risk area accounts for 42.3%，medium risk area accounts for 23.1%，and low risk area accounts for 34.6%.In formation area and circulation area of Maliu gully，the volume of scattered accumulation materials in the channel is about 300000 m3.The Maliu gully debris flow is still in a highly active stage with high risk.In the future，under the condition of sufficient rainfall and hydrodynamic force，a large quantity of loose materials may be reactivated，which may cause major geological disasters once it occurs.

Key words：

debris flow;hazard assessment;starting point of debris flow;Flow-R;FLO-2D;Dujiangyan City