基于数值模拟的耳阳河流域泥石流灾害危险性评价

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

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

泥石流是指由强降雨诱发，在地形陡峭、岩性松散的沟谷或者斜坡上形成的特殊洪流，表现为固体、液体、气体三相混合流体，一般具有历时短、强度高、破坏性大等特点[1]。山区泥石流灾害防控越来越受到重视，而泥石流危险性评价是泥石流灾害防控的关键环节[2]。泥石流危险性评价方法诸多，传统泥石流危险性评价方法主要有层次分析法、模糊综合评价、逻辑回归模型、物源可拓模型等。层次分析法确定指标权重具有很大的主观性；其它评价模型在泥石流危险性区划中以流域为基本单元，无法确定流域中精确地点的危险程度，也无法对泥石流的运动进行精确刻画；而数值方法评价的重点在于对泥石流流动过程的刻画，可识别流动速度、堆积厚度的空间分布及其随时间的变化，非常适合小流域中人口聚集区泥石流的危险性评价。宋兵等[3]运用RAMMS 模型模拟了白沙沟泥石流运动特征；Hungr[4]使用DAN3D 模型对泥石流的流动过程进行了模拟；Magirl 等[5]使用Laharz模型对圣卡塔利纳山脉7 条泥石流的运动距离和堆积范围进行了模拟，并与实际情况相验证；Horton 等[6]将二维动态泥石流模型Massflow 运用于汶川地震震中地区流域尺度的灾害图编制；Lin 等[7]采用二维洪水与土石流数值模拟软件FLO-2D 对台风诱发的松河地区泥石流进行了分析；常鸣等[8]运用FLO-2D 对都江堰八一沟泥石流进行了多种降雨频率条件下的危险性评价。上述数值分析方法在泥石流过程模拟方面均取得了较好的效果，尤其是FLO-2D 模型在泥石流危险性模拟方面潜力巨大。

泥石流灾害的致灾因子主要为强降雨，在相对脆弱的自然环境条件下，短时强降雨极易引发群发性泥石流[9]。2012年5月10日，甘肃岷县近一个小时大范围冰雹及强降雨诱发耳阳河流域突发泥石流，冲击沟口及河道两侧村落、公路、输电设备，造成18 人死亡，直接损失4 800 余万元[10]。经过与岷县多个气象站点历史降水记录相对比，确定诱发本次泥石流的降水强度约为百年一遇。

耳阳河流域泥石流沟众多，主沟平均纵比降为235‰，主沟和支沟塌岸严重。区内大量种植黄芪、党参、玉米、胡豆等作物。每年5—6月为农耕时节，农业生产后土壤松散，抗侵蚀能力差，一旦出现强降雨，容易发生群发性泥石流灾害。泥石流的威胁巨大，识别危险地段、预测致灾过程至关重要。因此，本文运用FLO-2D 模型，选取耳阳河流域6 条具有代表性的泥石流，模拟了2012年5月10日实际降雨（1%降雨概率）情况下的泥石流过程，并与地质灾害现场调查成果对比，以验证模拟的可靠性。进而，运用相同的评价模型，模拟上述泥石流在2.0%和0.2%降雨频率条件下的堆积范围、深度、流速，并进行了泥石流危险性分析，为泥石流灾害防治提供依据。

1 研究区概况

岷县位于甘肃南部，属于西秦岭北支中段，紧邻长江流域和黄河流域的分界线，洮河穿流而过，为典型的中高山侵蚀地貌，是甘肃省泥石流活跃区之一[11]。

研究区位于岷县洮河右岸支流耳阳河流域（图1），地处西秦岭山地和黄土高原的交汇区，新构造运动活跃，地貌上呈现为年轻的构造侵蚀中山地貌类型。耳阳河流域面积约63 km2，主沟长度约15.8 km，流域内最高海拔为东山顶的3 130 m，最低海拔2 310 m，位于耳阳河汇入纳纳河的沟口。耳阳河的主沟和支沟多呈“V”型谷，沟壁陡峭，临空面发育，滑坡塌岸多见；沟内耕地众多，土地松散，物源丰富。区内植被覆盖总体上差，耳阳河两岸岩土裸露，仅有农作物和零星灌木，上游局部可见人工林。

图1 耳阳河流域地形图Fig.1 Topographic map of the Eryang river watershed

流域内出露的地层主要是泥盆系、侏罗系的板岩、页岩和砂岩等，经过构造作用和风化作用，岩体破碎，稳定性较差；区内第四纪沉积物主要为黄土、冲洪积物、泥石流堆积物、残坡积物，耳阳河下游东岸（右岸）剖面图见图2。

图2 耳阳河右岸地质剖面图Fig.2 Geological section of the right bank of the Eryang River

研究区为高寒湿润气候带，多年平均降水量为560.8 mm，6—8月降水占总降水量的 60%以上。受地形及海拔高度影响，降水量空间分布差异较大，呈现出由河谷到山顶降水量逐渐增大的趋势。

耳阳河属于黄河水系，为黄河的三级支流：耳阳河在茶埠镇沟门村汇入纳纳河，之后向西流约500 m，在茶埠镇阳坡汇入洮河。经过野外调查，耳阳河流域内共发育23 条泥石流沟，其中左岸9 条，右岸14 条，泥石流沟线密度为1.45 条/km。

本研究选取的6 条泥石流沟对称分布于耳阳河两岸（表1），均处于活跃期，呈“V 形谷”，主沟纵比降为189.7‰～316.9‰，平均坡度35°～50°。松散物质类型主要为坡面的残坡积物以及耕作土，其次为沟道两岸的小型崩塌、滑坡及沟岸坍塌，其余为沟底再搬运形成的松散物质。6 条沟谷中都有大面积耕地，占沟谷总面积的51%，义子沟、接哈沟、扎龙沟上游农田已经退耕，牌嘴堡北沟（牌嘴沟）和接哈沟上游有人工林。6 条沟沟底情况相似，可见泥石流相砂砾，碎石以及泥球，碎石粒径多在2～10 cm，最大可见1.3 m，母岩主要为砂岩、砂质板岩、灰岩、灰质板岩。6 条泥石流沟的补给段长度占各自沟道总长度的61.8%～69.2%。接哈沟堆积扇挤压耳阳河河道，使耳阳河出现一处牛轭状河曲（图3）。

表1 泥石流沟特征参数Table 1 Property list of the key debris flow gullies

图3 本研究选取的六条泥石流平面分布图Fig.3 Distribution of the six gullies in this study

2012年“5·10”特大暴洪泥石流发生之后，在牌嘴沟中修建了10 道拦挡坝。2019年9月现场调查发现，前6 道拦挡坝已被泥石流堆积物、两侧崩滑堆积物淤满、掩埋，后4 道拦挡坝仍具有拦挡作用。在扎龙沟沟口修建了约400 m 长的排导渠，宽10～15 m，深3 m。

2 泥石流危险性评价

FLO-2D 模型由O’Brien[12-13]提出，采用有限差分计算垂向深度和流动速度，预测泥石流流动和堆积范围，可用于泥石流、洪水、溃坝、城市淹没等过程的模拟，并可用于灾害危险性评价。

FLO-2D 软件在模拟泥石流运动过程中，将地形分为若干等大网格，每一个网格中，其曼宁系数和高程值都是唯一的，通过运动方程以及连续方程，可以计算出每个网格中流体深度、流量，进而得知流体的运动范围，通过动量方程计算出相邻网格间流体的速度变化。FLO-2D 模型方程如下：

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

x—X轴方向距离/m；

y—Y轴方向距离/m；

h—流体深度/m；

I—有效降雨强度/（mm·h-1）；

u—X轴方向流速/（m·s-1）；

v—Y方向流速/（m·s-1）。

连续方程控制了泥石流运动时每个网格内质量守恒。

式中：g—重力加速度/（m·s-2）；

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

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

X、Y方向的运动方程保持了运动平衡。

式中：Sf—摩擦坡降；

Sy—屈服坡降；

Sv—黏性坡降；

Std—分散坡降；

η—黏滞系数；

k—层流阻滞系数；

γm—泥石流容重/（kN·m-3）；

ntd—等效曼宁系数。

流变方程考虑了泥石流运动时颗粒之间的碰撞对泥石流流动阻力的影响。

2.1 数据处理与导入

2.1.1 地形数据处理

地形数据来自1∶5 万数字高程模型（DEM）。将高程数据转化为ASCII 格式，导入FLO-2D，将之剖分为10 m×10 m 的评价单元，通过插值计算，确定每个评价单元的高程，完成地形数据处理。

用ArcGIS 将研究区流域边界、出水点、拦挡坝等主要地形要素矢量化，导入FLO-2D。

2.1.2 确定泥石流参数

郭富赟等[14]对耳阳河流域泥石流流体重度进行了测算，为17.0～17.6 kN/m3。本文采用文献[14]的泥石流流体重度测算结果，并结合堆积物平均重度，计算出泥石流泥沙体积浓度为46.67%～50.67%，取中间值48.67%。

高浓度泥石流在屈服应力的作用下出现层流特征，层流之间的摩擦作用即可用层流阻滞系数（k）表示，本文参考前人研究实例，用工程地质类比法，k取值定为2 250。

黏滞系数（η）和屈服应力（τy）的大小主要取决于泥沙体积浓度，关系式如下：

式中：α1、α2、β1、β2—试验系数；

Cv—体积浓度。

根据泥石流堆积物、流体的容重和物质组成，判断研究区泥石流为黏性，且为中阻性，泥沙比（Rns）为0.75。王裕宜等[15]统一的泥沙比-体积浓度-流变参数关系式如下：

式中：nc—曼宁系数；

h—泥深/m（本文以平均堆积厚度表示）。

将各泥石流沟泥深带入式（7）即可得到拉龙沟、扎龙沟、拉路沟、接哈沟、义子沟、牌嘴沟的曼宁系数，分别为0.019 3，0.022 4，0.004 5，0，0，0.022 4。将式（5）（6）（8）（9）联立，求解得到α1、α2、β1、β2分别为0.000 16，0.002 30，17.41，18.28。

2.2 “5·10”泥石流数值模拟

根据野外调查和访问得知，耳阳河流域泥石流于2012年5月10日17 时50 分许开始爆发，18 时左右达到洪峰，21 时逐步回落，泥石流过程历时约3 h。降雨量数据选用岷县麻子川自动站降雨观测数据，如图4所示[16]。流量过程曲线计算选用《四川省中小流域暴雨洪水计算手册》中川西地区水文模型，本次研究区位于甘肃省南部，和川西地区距离较近，地质、地貌、气候均有相似性，此计算手册中的参数在本区降雨-径流分析中具有良好的适用性。运用推理公式法概化各泥石流沟清水流流量过程曲线，并乘以放大系数BF（BF=1/（1-Cv）=2.10），得到泥石流流量过程曲线，见图5。

图4 岷县麻子川自动站降水曲线图（2012年5月10日）Fig.4 Precipitation curve of the Mazichuan monitoring station（May 10,2012）

从图5可以看出，各泥石流沟于爆发15～30 min后达到洪峰，约3 h 后流量逐步回落；这与现场调查和访问了解到的“5·10”当时的情况基本相符。

图5 清水流量过程曲线（a）和泥石流流量过程曲线（b）Fig.5 Flow process curve of（a）clear water and（b）debris flow

模拟结果显示：泥石流流体在流通区运动速度较快，最大速度可达7 m/s；在沟口附近速度迅速下降，大部分泥石流汇入耳阳河，部分区域堆积深度较大，堵塞河道（图6）。

图6 “5·10”泥石流深度（a）和流速（b）Fig.6 “5·10” debris flow（a）depth and（b）velocity

通过FLO-2D 模拟，得到了6 条泥石流沟的堆积扇面积和泥石流堆积平均厚度，如表2所示。实地调查测量得到的面积和厚度同时列于表2，并将模拟结果和实际面积进行对比。可以看出：除拉龙沟和拉路沟外，其他各沟堆积扇面积误差率在-6.2%～28.5%之间、堆积扇均厚误差率在-7.5%～23.3%之间。模拟堆积厚度大部分低于实际厚度，主要原因是自“5·10”泥石流爆发后，每年有少量松散物质冲出沟口，堆积扇厚度逐年加大，实测的厚度比“5·10”泥石流实际堆积厚度大。拉龙沟与拉路沟堆积面积误差率分别达到189.8%和2 827.5%，模拟堆积扇面积远远大于实际面积，主要原因为村民在沟口附近建房、开垦土地，对泥石流堆积扇进行了人为改造。

表2 模拟结果与调查结果对比表Table 2 Comparison of the simulation and survey results

2.3 “5·10”泥石流危险性评价

当前泥石流危险性分级指标主要有流速和泥深[17]、冲击力和堆积厚度[18]、频率和强度[19-20]等，耳阳河流域泥石流的危害性主要在于泥石流携带的固体物质淤埋农田和房屋，以及中高速黏性洪流对沿岸的冲击和侵蚀。因此，本次选取流速、泥深等指标进行危险性分级，采用如表3所示的逻辑关系，确定泥石流危险性的空间分布。

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

危险性评价结果见表4、图7。“5·10”泥石流总危险区面积为57 900 m2，占总流通区面积的60.2%，其中高危险区面积占总危险区面积比例较高的是接哈沟和扎龙沟，分别为34.23%和38.60%。

表4 泥石流危险性分区统计表Table 4 Statistics of debris flow hazard zoning

由于“5·10”泥石流过后重新规划了房屋建筑，目前耳阳河沿岸的大部分建筑在1%降雨频率条件下较为安全，但扎龙沟、接哈沟沟口和义子沟沿岸的部分房屋建筑易受到泥石流威胁，共涉及约21 处房屋（见图7中的局部放大）。

图7 “5·10”泥石流危险性分区图Fig.7 Hazard zoning map of “5·10”debris flow

2.4 不同降雨频率条件数值模拟和危险性评价

通过对“5·10”过程进行模拟，并与实际情况进行对比，验证了模型和参数的可靠性，说明FLO-2D 技术能够较好地重现泥石流过程、模拟出泥石流流速、堆积深度的时空分布。使用同样的参数及方法，对P=2.0%及P=0.2%降雨频率条件下泥石流情况进行了模拟，并开展了上述两个工况的泥石流危险性评价。P=2.0%及P=0.2%降雨参数取值见表5，模拟得到的流量过程曲线见图8。

表5 不同频率的降雨参数Table 5 Rainfall amounts of different frequencies

图8 不同降雨频率条件下泥石流流量过程曲线Fig.8 Flow process curve of debris flow in different rainfall frequencies

利用表3中的逻辑运算方法，得到了P=2.0%及P=0.2%降雨概率的泥石流危险性分区图（图9、表6）。

图9 2.0%（a）和0.2%（b）降雨频率条件下泥石流危险性分区图Fig.9 Hazard zoning map of debris flow in（a）2.0% and（b）0.2% rainfall frequencies

表6 2.0%和0.2%降雨频率条件下泥石流堆积区危险性分区统计表Table 6 Statistics of hazard zoning of debris flow accumulation area in 2.0% and 0.2% rainfall frequencies

模拟结果显示：在2.0%降雨频率条件下，泥石流总危险面积为37 900 m2，牌嘴沟、拉路沟、义子沟无高危险区，加之牌嘴沟已修建道拦挡坝，危险性不高。低危险区占总危险区的72.2%。扎龙沟高危险区占比最高，为13.95%，主要分布于沟口，距离左岸居民建筑最近约30 m，威胁较小。但义子沟沿岸仍有4 户居民处于中危险区内，接哈沟沟口仍有5 户居民处于中危险区内。

在0.2%降雨频率条件下，泥石流危险区总面积为60 500 m2，与百年一遇条件下相比，危险区总面积仅增加了4.5%（2 600 m2）；但中-高危险区面积占比增大，其中接哈沟堆积区高危险面积约4 000 m2，占总面积的34.78%；其次为义子沟和扎龙沟，高危险区面积分别为3 100 m2和2 500 m2，占比分别为12.82%和35.82%。考虑到“5·10”后，在牌嘴沟已增设拦挡坝，牌嘴沟堆积区无高危险区。除拉路沟外，其他各沟沟口和沟道附近居民均受到威胁，共约38 户，其中中-高危险区居民总户数约22 户。

3 结论

（1）利用自动站降雨观测数据，运用FLO-2D 模拟重现了“5·10”泥石流运动和堆积特征，得到了泥石流堆积范围、堆积深度、发展过程。与野外调查所得的现场情况进行了对比，排除人为大规模改造作用，堆积扇均厚模拟结果误差在-7.5%～23.3%之间，堆积扇面积模拟结果误差在-6.2%～28.5%之间，数值模拟结果和实际情况基本相符，模拟效果较好。

（2）通过模拟，以泥石流堆积深度和冲击速度为标准，划分了不同降雨频率场景下的泥石流危险区；在此基础上圈定了易于受泥石流威胁的房屋、财产等潜在受灾对象，为泥石流风险防控及治理提供了参考。

（3）FLO-2D 模型能够较好地进行泥石流过程重现和危险性区划，通过不同降雨频率下的数值模拟，可以得到泥石流流体深度和速度的空间分布和变化过程，能够对泥石流过程及冲击场景进行直观表达。