基于数值模拟的群发性泥石流危险性评价——以甘肃岷县麻路河流域为例

曹 鹏，侯圣山，陈 亮，冯 振，王立朝，李 昂，刘军友

（1.中国地质大学（北京）水资源与环境学院，北京 100083；2.中国地质环境监测院(自然资源部地质灾害技术指导中心)，北京 100081）

0 引言

泥石流为固、液、气三相混合流，一般在山区沟谷或坡面，由强降雨、洪水等自然或人为因素引发，具有历时短，速度快、冲击大、破坏强的特点[1]。2000——2020年，岷县由于泥石流灾害造成的死亡人数共93 人。2012年5月10日，岷县大范围遭受强降雨夹冰雹，局部遭受百年一遇强降雨，引发大量泥石流灾害，造成18 人死亡[2]。

据2012年11月“岷县地质灾害详细调查”成果，麻路沟流域2012年5月10日遭受五十年一遇降雨，引发群发性泥石流，主河泥位1.2 m，灾情较严重，损失牲畜33 头，毁坏房屋5 间，农田4.13×104m2，道路3 km，直接损失约37 万元，无人员伤亡[3]。经野外现场调查，麻路河流域“5·10”泥石流堆积区特征明显，“5·10”后又暴发过多次泥石流和洪水灾害。麻路沟是迭蔵河的支流，位于岷县县城上游3 km，麻路河流域暴发泥石流，不但对沟内居民和财产造成威胁，还会威胁到沿迭蔵河分布的岷县县城区。因此，迫切需要针对麻路沟流域开展危险性评价，识别高危险性泥石流沟，划定泥石流中高危险区域，为制定准确有效的泥石流防治措施提供依据。

基于数值模拟的泥石流危险性评价方法可以具体刻画泥石流的运动特征，每个泥石流流域单元危险度是非唯一的，从而可以区分流域内不同区域的危险等级，如胡凯衡等[4]以流体力学模型为基础，用数值模拟的方法较精确地描述了泥石流流体运动和物理机制，段学良等[5]运用Massflow 对的西藏仁布杰仲沟泥石流运动特征进行分析，侯圣山等[6]运用FLO-2D 对甘肃岷县耳阳河流域内6 条泥石流沟进行危险性评价，杨涛等[7]运用FLO-2D 对岷江小流域内6 条泥石流沟进行危险性评价，颜恒明[8]对都江堰市虹口乡红色村干沟进行了数值模拟。但是，如上文所例的入汇河流的单沟泥石流或群发泥石流，并未考虑强降雨过程中河流对泥石流堆积物的冲刷携带作用，故而对准确预测泥石流堆积扇范围、淤积厚度具有不利影响。

本文选取岷县麻路河流域中6 条泥石流沟和麻路河作为研究对象，首先建立主河水流模型，并基于当地平均雨强和差异系数，运用川西地区水文模型建立1%、2%、5%降雨频率条件下泥石流流量过程曲线和流域上游支沟的清水流量过程曲线（2012年麻路河“5·10”泥石流降雨频率为2%），在主河水流模型的基础上对全流域施加不同频率降雨进行泥石流模拟，即可考虑主河对泥石流堆积物的冲刷携带作用，将2%降雨频率条件下各泥石流堆积扇均深和范围与野外调查结果对比，即可验证模型可靠性。根据模拟结果得出全流域（6 条泥石流沟，2 条清水支沟，1 条主河）危险性分区图，针对杨家沟、叶家沟、张家沟、大湾沟、峡里沟、拉路沟等6 条泥石流沟进行了单独的危险性评价，并圈定流域内受中高危险威胁的居民区，为灾害的精确防治和预警提供依据。

1 研究区概况

岷县地处甘肃南部，为陇南山地与黄土高原接壤带，属西秦岭北支中段，洮河流经县内西北部，境内泥石流灾害活动非常活跃[3]。

研究区麻路河流域位于岷县西南部，流域东南边界为黄河-长江流域分水岭，区内新构造运动强烈，形成典型的中山构造-侵蚀地貌。麻路河流域面积约33.29 km3，主河长6.5 km，平均坡降120‰，流域平均宽度4.2 km，最高海拔3 180 m，位于东山，最低海拔2 400 m，位于麻路河汇入迭藏河沟口。

根据野外调查及勘查结果，流域内地层主要为新近系青灰色砾岩砂岩及砖红色粉砂质泥岩，三叠系板岩页岩，由于构造运动活跃，风化作用强烈，岩体结构面及风化裂隙十分发育，岩体呈现为碎散状，稳定性极差，流域内第四系沉积物主要为马兰黄土、残坡积物、滑坡堆积物、泥石流堆积物等。

麻路河流域位于温带半湿润-高寒湿润气候过渡带，年均相对湿度为 65%，年均降水量在500～600 mm之间，降水时间分布不均，春夏季降水占年降水的 60%以上。由于地形和海拔影响，降雨量总体趋势为由高山区向河谷区递减，降雨历时一般在6～12 h，降雨笼罩面积越小，其降雨强度越大。

麻路河属黄河水系，为洮河的二级支流，麻路河由东向西从白土坡汇入迭藏河后，向北流6.12 km 后汇入洮河。岷县县城沿洮河和迭蔵河展布，麻路河汇入迭蔵河的位置为岷县火车站所在地，下游3 km 即为岷县县城区（图1）。

图1 研究区地理位置图Fig.1 Location map of study area

麻路河流域范围内共有8 条支沟，其中6 条为泥石流沟，均分布在主河东岸，其余2 条清水支沟位于流域源头（图2）。6 条泥石流沟均处于活跃期，中上游沟谷截面呈“V”形，下游沟谷截面呈“U”形，泥石流沟纵比降在169‰～366‰。麻路河流域地形平均坡度在20°～30°（图3），沟谷水源汇流区坡度较大，在30°～40°，据前人研究，最易发生泥石流的地形坡度在25°～45°[9]。麻路河流域地形见图2，图中沟谷编号1～8 分别代表杨家沟、叶家沟、张家沟、大湾沟、峡里沟、拉路沟、清水支沟1 和清水支沟2。麻路河泥石流沟具体特征参数示于表1。

图3 地形坡度分级图Fig.3 Topographic gradient classification map

表1 泥石流沟特征参数表Table 1 Characteristic parameters of debris flow gullies

麻路河流域泥石流固体物质来源主要有四种类型：坡面物源（耕地）、崩塌滑坡物源、沟道物源以及人工弃渣物源。每年5——6月是岷县春耕时节，坡面土地耕后松软，在强降雨条件下，大量松散物质易随坡面径流汇入沟谷，这是本区域泥石流固体物质的最重要来源，经统计，流域内耕地面积占总流域面积的58.5%，约19.47 km2，麻路河南岸耕地开垦于河流冲洪积物和崩坡积次生黄土，平均厚度为1.8 m，北岸耕地开垦于崩坡积次生黄土，部分开垦于紫红色泥岩和板岩全风化碎屑堆积物，平均厚度为1.3 m；易于被径流冲走的松散物质主要是耕后土壤，研究区耕作深度约22 cm，据此测算麻路河流域泥石流沟坡面耕地可提供的松散物质共146.94×104m3。6 条泥石流沟岸两侧塌岸非常发育，尤其各支沟上游连续塌岸现象非常普遍，沟内滑坡大多发育于古老的泥石流堆积体及残坡积体上，经野外调查，麻路河流域塌岸、崩塌、滑坡共56 处，总体积为23.96×104m3（图2，图4）。沟道物源主要来源于近期泥石流淤积、河流再搬运及流水对沟谷两侧古老泥石流堆积物的侵蚀，各沟沟底淤积深度为0.8～1.6 m，沟道松散物质总量为17.34×104m3。杨家沟流通区有烧砖厂，弃渣堆积于沟内，堆积面积约488 m2，体积约0.19×104m3，峡里沟上游左岸有一处采石场，占地面积8 500 m2，堆积碎石屑及碎石块约1.72×104m3，在强降雨条件下，碎石屑及部分石块可被坡面径流搬运至沟内。各泥石流沟松散物质储量见表2。

图2 麻路河流域平面示意图Fig.2 Sketch map of Malu River basin

图4 泥石流沟松散物质Fig.4 Loose material in debris flow gullies

表2 麻路河流域泥石流沟松散物质量（单位：104 m3）Table 2 Volume of loose material in debris flow gullies of Malu River basin (unit: 104 m3)

水是泥石流形成不可缺少的重要条件。麻路河流域“5·10”降雨由17 点到20 点，持续约4 h，总降雨量为43.8 mm，达到50年一遇的强度，其中18 点到19 点一小时降雨量为21.8 mm，而据前人研究，岷县可能发生泥石流的24 h 的界限雨量H24(H)=30 mm，1 h 的界限雨量H1(H)=15 mm[3]，麻路河流域“5·10”降雨远远超过界限雨量。除大湾沟外，其他5 条泥石流沟内都有季节性水流，流量在6——8月份达到最大值，约0.3～0.9 m3/s。

2 FLO-2D 模拟

O’BRIEN[10]提出并验证了FLO-2D 数值模型，FLO-2D 可以用于二维流体模拟，预测流体运动特征，在实际应用中，可以应用到洪水、泥石流的模拟、溃坝分析以及城市淹没等方面。

泥石流流体运动特征和堆积特征主要由流体流动速度、流量、流深、流动范围等四个指标表示。FLO-2D 实现四个指标解算的原理是：将评价区域划分为若干方形单元，每个单元均有各自的高程和粗糙度，连续方程（式1、2）和流体动量方程（式3）可计算出每个单元的流深、流量，通过各单元的流深数值可以确定泥石流流动范围，流体动量方程可以确定相邻单元之间流体的流速。

式中：g——重力加速度;

Sfy、Sfx——X、Y方向摩擦坡降；

Sox、Soy——X、Y方向床底坡降。

式中：t——流体运动时间/s；

h——流深/m；

i——有效雨强/(mm·h-1)；

u——X方向流速/(m·s-1)；

v——Y方向流速/(m·s-1)。

2.1 参数选取

地形数据采用高精度1∶10 000 数字高程模型，将DEM 转化为同等精度的ASCII 码，结合ArcGIS 对流域边界、水系、泥石流起动点、居民区进行矢量化，将ASCII 码和矢量数据导入到FLO-2D 后，把高程模型划分为20 m×20 m 的网格并进行插值计算，最后圈定计算区域，完成数据处理。

根据勘查成果，麻路河泥石流堆积物平均重度为2.53 g/cm3，泥石流性质为黏性。据郭富赟等[11]对岷县耳阳河泥石流流体重度测算结果17.0～17.6 kN/m3，研究区与耳阳河流域相邻，泥石流堆积物物质组成和粒度成分相近，故本文泥石流流体容重取值17.3 kN/m3，求得流体平均泥沙体积浓度为52.29%，由式（4）可求得，泥石流堆积物体积膨胀因子BF为2.10。

黏性泥石流流体体积浓度高（一般大于50%），由于屈服应力，会产生层流作用，层间摩擦作用以层流阻滞系数表示，根据WOOLHIS（1975）的方法[12]测算结果（表3），研究区层流阻滞系数取值范围为1 000～4 000，且数值随地面粗糙度增大而增大，确定研究区层流阻滞系数为2 250。

表3 层流阻滞系数表Table 3 Laminar retardation coefficient

在流变模型中，流体屈服应力τ 和黏滞系数η 与泥沙体积浓度cv具有指数关系，即流体浓度越大，其屈服应力和黏滞性越强。根据堆积物的物质组成及黏粒含量，将麻路河泥石流进一步划分为中阻黏性泥石流。王裕宜[13]根据不同阻力的黏性泥石流的泥砂比，获得统一的流变参数计算式：

式中：nc——曼宁系数；

R——泥沙比（低阻泥石流：R=0.96；

R——中阻泥石流：R =0.75，高阻泥石流： R =0.29）；

h——泥深/m。

杨家沟、叶家沟、张家沟、大湾沟、峡里沟、拉路沟沟道泥深分别为2.4，2.5，2.5，1.0，1.9，1.3 m，由式（5）可求出：上述沟谷的曼宁系数分别为0.046，0.048，0.048，0，0.033，0.014。主河道曼宁系数根据FLO-2D 手册设为0.15。已知R=0.75，联立式（7）（8）（10）（11）得：α1=0.000 16、β1=17.41、α2=0.023、β2=18.28。

2.2 建立模型

本研究的研究区处于甘肃省南端，地形、地质、气象条件均与川西地区情况接近，岷山周边的泥石流特征也与陇南、川西地区特征相似，因此本研究运用川西水文模型（《四川省中小流域暴雨洪水计算手册》）计算流域内6 条支沟的基流量，基于各沟基流量通过FLO-2D 模拟主河水运动状态，建立水系模型（图5），并在野外实际观测点位置（图5中A、B 两点）提取对应单元流速、流深数据（图6）。

图5 降雨前流域水系模型Fig.5 River system model before rainfall

由图6可知，泥石流暴发3.5 h 左右附近支沟水流首先汇流到A 点，之后由于上游水流汇入主河，导致4-6 h 内流速与流深发生变化，约6.5 h 后流速与流深达到稳定状态，整个流域的汇流时间为6.5 h。结合A、B 两点流速与流深，可计算两点流量分别为5.41 m3/s和4.00 m3/s，与野外观测所得情况基本相符。

图6 观测点流速、流深Fig.6 Flow velocity and depth at observation point

通过岷县地区泥石流历史发现，本地群发性泥石流为单峰型，涨水快退水慢。基于岷县当地最大降雨强度（表4），运用推理公式法，得到六条泥石流沟在不同降雨频率条件下的清水流量过程曲线，已知泥石流体积膨胀因子BF为2.1，可求出不同降雨频率条件下泥石流流量过程曲线（图7）。由图7可知，2%降雨频率条件下麻路沟流域6 条泥石流在0.5～1 h 内先后达到洪峰，在3～4 h 后逐渐衰落，总体历时约4 h，麻路河流域 “5·10”泥石流暴发于10日18 时左右，18 时30 分左右达到洪峰，21 时先后回落。

表4 不同频率的降雨参数（单位：mm）Table 4 Rainfall parameters at different frequencies（unit:mm）

图7 不同降雨频率条件下泥石流流量过程曲线Fig.7 The flow process curve of debris flow under different rainfall frequency

2.3 模拟结果及验证

首先对2%降雨频率条件下（即2012年“5·10”泥石流降雨条件）麻路河流域各泥石流沟和主河进行模拟，得出流域中泥石流、洪流的流速及流深情况（图8）。

图8 2%降雨频率条件下泥石流模拟结果Fig.8 Simulation results of debris flow under the condition of 2% rainfall frequency

为验证模型的可靠性，将2%降雨频率条件下模拟堆积扇面积、均厚、位置及形状与野外调查结果相对比：堆积扇面积和均厚验证方法为求取模拟与调查结果的误差比率；堆积扇位置、形状验证方法为对模拟和调查堆积扇范围进行几何运算，求取数值模拟精度因子Ac（式8），Ac越接近1，表明数值模拟与实际情况越吻合。

式中：S0——野外调查堆积扇与模拟堆积扇重叠面积/m2；

Sa——野外调查堆积扇面积/m2；

Sn——模拟堆积扇面积/m2。

h——泥深/m。

由对比结果可知（表5、6）：除峡里沟和拉路沟外，堆积扇面积的误差比率为2.15%～18.03%，堆积扇均厚的误差比率为-28.31%～-11.74%，数值模拟精度因子Ac为0.71～0.80，模拟结果较准确。数值模拟泥石流堆积扇面积稍大且厚度稍小，主要原因为堆积扇边缘部分经过人为改造而难以确定，故平均厚度的取值往往趋向于堆积扇中部，且自2012年以来每年有少量泥沙冲出沟口也导致了堆积扇厚度增大。峡里沟和拉路沟由于地形限制形成共同的泥石流堆积扇，且堆积扇人为改造作用强烈，几乎除河道外皆建有房屋，沟道人为扰动强烈，可识别堆积扇范围小，导致堆积扇范围及均厚误差比率高达35.75%和-67.91%，数值模拟精度因子仅为0.63。

表5 模拟-调查误差率表Table 5 Simulation & survey comparison table

表6 数值模拟精度表Table 6 Numerical simulation accuracy table

总体上看，在天然水系上叠加降雨更加符合实际情况，考虑了主河洪水与冲出泥石流的相互作用，且为整个流域提供了洪水、泥石流的运动情况，经验证天然工况下水系模型及降雨工况下泥石流模型模拟结果与现场符合性较好，模型及参数较为可靠。

3 危险性评价与讨论

由于2%降雨条件下模拟结果较好，故运用同样方法及参数对5%和1%降雨频率条件下泥石流进行模拟，根据模拟结果对三种情况进行危险性评价。

泥石流的破坏形式主要为淤埋、冲击和摩擦三种方式[14-17]。泥石流影响强度则和泥石流深度、流速、堆积有关。本研究运用基于泥石流影响强度的危险性分区方法（表7）进行危险性分区（图9）。

表7 泥石流危险性分区指标Table 7 Risk classification of debris flow

图9 1%、2%、5%降雨频率下危险性分区图Fig.9 Hazard zoning map in 1%、2%、5% rainfall frequency

由危险性分区统计表（表8）得知，在不同降雨频率条件下危险性最高的泥石流沟为峡里沟和拉路沟，其共同堆积扇高危险区面积比例超过60%，大湾沟和杨家沟堆积扇高危险区面积和比例较小，危险性最低。

表8 1%、2%、5%降雨频率条件下泥石流堆积区危险性分区统计表Table 8 Statistical table of hazard zoning of debris flow accumulation area in 1%、2%、5% rainfall frequency

对1%降雨频率条件下危险性分区图进一步分析，圈定出6 处受高度威胁的居民区见图9（a），分别位于叶家沟、峡里沟、拉路沟泥石流堆积区，张家沟泥石流流通区和麻路河中下游河曲附近，危险区分布细节示于图10，已识别出百年一遇降雨条件下泥石流危险区范围内的具体房屋（图10）。图10 中Ⅱ、Ⅲ、Ⅵ三个区域，泥石流破坏形式主要为冲击和淤埋，距离泥石流沟口和沟道近的房屋建筑首先受到强力冲击，距离远的房屋建筑受到不同程度淤埋；图10 中Ⅰ、Ⅳ、Ⅴ三个区域，房屋建筑建于河曲一侧，暴涨的洪水冲刷携带泥石流堆积物，在河曲处由于离心力的作用冲出河道，冲击、淤埋沿岸居民建筑。在GIS 中测量，图10 中6 处高危险区受泥石流严重威胁的房屋面积共约为49 753.13 m2。

图10 百年一遇降雨情况下高危险区域细节图Fig.10 The details of high risk area details (P=1%)

4 结论

本文基于岷县当地降雨数据和川西水文模型，通过FLO-2D 建立麻路河流域水系模型和6 条群发性泥石流模型，考虑主河洪水对泥石流堆积物的冲刷携带作用，模拟2%降雨频率条件下流域中泥石流的运动和堆积特征，以泥石流模拟结果与实际“5·10”泥石流调查结果相对比，证明了泥石流堆积扇模拟范围、厚度、形状、位置的准确性，故进一步以相同模型对1%和5%降雨频率条件下泥石流进行模拟并进行危险性分区，确定麻路河流域高危险泥石流沟，选取百年一遇的强降雨的情况（P=1%），进一步在遥感图上圈定高危险居民区并分析各自潜在受灾类型，为泥石流灾害综合防治提供依据。

（1）1%、2%、5%降雨频率条件下流域内堆积扇危险区总面积分别为0.126，0.115，0.109 km2，即降雨频率越小，堆积扇危险区总面积越大，但相差数值较小，故流域内构造-侵蚀微地貌对泥石流堆积扇范围具有控制作用。

（2）根据1%、2%、5%降雨频率条件下泥石流模拟结果进行危险性分区并对堆积扇危险面积进行统计，得出在不同降雨频率条件下危险性最高的泥石流沟为峡里沟和拉路沟，所以应将峡里沟和拉路沟作为泥石流防治工作的重点。

（3）通过重点受灾分析及主河流速图，麻路河河曲处洪水流速明显加大，且多为受灾部位，应将河曲部位作为泥石流、洪流重点防治区域。

（4）针对全流域的群发性泥石流危险性评价，考虑主河洪水作用的方法是可行可靠的。