0.   引言

2008年“5·12”汶川地震产生大量的崩塌滑坡体，为泥石流活动提供丰富的固体松散物(Dai et al., 2012).震后有关研究发现泥石流暴发所需的激发雨量大为降低，增加了震区泥石流活动的危险性(Chang et al., 2014).目前国内外已有学者展开震区泥石流危险性评价的研究，Christen et al.(2010)采用泥石流发生频率方法开展泥石流危险度评价研究.Nocentini et al.(2015)结合DAN-W与FLO-2D模型，运用动态代码的方法快速评估研究地区潜在泥石流的危险性. Quan et al.(2011)利用FLO-2D模型模拟获取泥石流危险性评价结果，以此建立易损性函数曲线.唐川等通过野外实测数据确定了典型泥石流危险性评价因子(唐川等, 2011).詹钱登等(2015)充分利用FLO-2D模型对台湾地区泥石流的危险性展开评价，并通过模拟结果与现场大量实验对比获得了不同降雨模式下模拟所需的相关参数.唐川等通过对震后滑坡物源分布特征及降雨阈值强度的研究，判断不同降雨频率下泥石流的危险强度(唐川等, 2010; Tang et al., 2011; Zhou and Tang, 2014).杜军等充分结合GIS与信息量模型开展汶川县城次生地质灾害危险性评价(杜军等, 2010).崔鹏等构建了一套基于山洪泥石流动力过程的危险评价与风险管理理论与方法体系(崔鹏和邹强, 2016).

2010年8月13日强震区都江堰龙池场镇突发强降雨，导致区内48条泥石流同时暴发，大量固体松散物质伴随泥石流冲入龙溪河，造成沟口大量民房不同程度受损，道路、河堤及农田被冲毁或者掩埋，部分泥石流出现了堵江现象，其中威胁最大的为八一沟泥石流，流域内11座拦砂坝被冲毁，造成2人失踪1人受伤，约1 km路面及100余间房屋被冲毁或掩埋，直接经济损失约1 500万元.为了研究八一沟泥石流形成机制及运动模式，很多学者展开过相应的探索.马煜等通过灾后实地调查取样计算了八一沟泥石流的静力学及动力学特征，得到引发“8·13”八一沟泥石流的激发雨量相当于20年一遇，即激发八一沟泥石流的为20年一遇的降雨量(马煜等, 2011).张自光等通过现场试验法分析了八一沟泥石流的物理特性，其容重为1.77 t/m³，属于粘性泥石流.通过泥石流易发程度评分表得到八一沟泥石流得分为102分，属易发泥石流，整体处于发展期，当遇到持续降雨或短时强降雨条件，很可能再次暴发特大型泥石流(张自光等, 2010).周伟等初步分析了导致八一沟泥石流暴发的地形、物源及降雨三大条件，通过动力学参数并结合单沟泥石流危险度评价方法，开展了危险度评价.结果表明八一沟泥石流处于高危险度等级，流域内堆积大量的松散固体物质，在强降水条件下暴发大规模泥石流可能性仍然存在(周伟等, 2011).刘希林等以八一沟泥石流为例重点探讨了建立泥石流风险评估模型，根据泥石流事件发生概率与其造成的负面结果的综合为其风险性的原理，建立了泥石流综合风险评估关系式模型(刘希林和尚志海, 2012).覃怡等采用野外调查与遥感解译相结合的方法简要论述地形地貌、物源与水源三个因素相互结合共同导致了八一沟泥石流形成(覃怡和郑洪春, 2013).

经过现场调查及前人研究表明，八一沟泥石流流域内物源丰富、地形陡峭，在极端强降雨天气条件下，具备再次暴发的可能性.八一沟泥石流一旦再次暴发，将对云华村113户460名村民的生命财产安全、道路、农田及河堤带来严重威胁.为了减轻泥石流灾害带来各项损失，及时有效的提供防灾减灾对策，结合前人的工作，在前人已经获取的八一沟泥石流各项参数的基础上，采用非牛顿流体模式及中央有限差分的数值方法进行数值模拟，充分结合泥石流自身强度及暴发频率，建立八一沟危险性评价模型，得到在不同降雨强度下的八一沟泥石流危险性分布图，及时为八一沟泥石流灾后重建及重点防范区域提供可靠建议.

1.   八一沟流域概况

八一沟泥石流位于汶川强震区四川省都江堰市龙池镇龙溪河流域，整体流域面积约为8.6 km²，主沟全长约4.5 km，流域内最高海拔为2 456 m，最低海拔851 m，高差可达1 605 m，沟床平均纵坡降约为376‰，整体呈“树叉状”分布，主沟纵坡降较小，横断面呈“U”字型，以淤积为主，支沟纵坡降较陡，横断面呈“V”字型，以冲刷为主(图 1).八一沟流域内出露的地层岩性主要包含第四系崩坡堆积物（Q₄^e+dl）、第四系洪积物（Q₄^pl）及第四系崩滑堆积物（Q₄^del）.元古代澄江-晋宁期的中粒斜长花岗岩(γ_o2⁽⁴⁾)与震旦系下统火山岩组（Za）的灰绿色安山岩主要分布在形成区，三叠统须家河组（T₃x）的砂岩、泥岩主要分布在流通区，以碎块石、凝灰岩及安山玄武岩为主的第四系地层主要分布在堆积区.八一沟泥石流形成区内沟道冲刷明显，两侧坡体表层剥蚀严重，底部基岩揭底出露，流通区内沟道内堆积物表面颗粒粒径较大，揭底现象十分明显，大量树木折断反映其冲击力较大，堆积区淤积现象明显，大体积块石较多，致使龙溪河半堵.灌县断裂和虹口-映秀断裂穿过流域，地质构造极为复杂，岩体支离破碎，提供了大量的松散固体物质.其所在区域地震活动较为频繁，根据国家地震局2008年颁布的《四川、甘肃、陕西部分地区地震动参数区划图》得到八一沟泥石流流域基本地震烈度为Ⅷ度，地震动峰值加速度值为0.20 g，其反应谱特征周期为0.40 s.

图  1  八一沟流域平面示意

Fig.  1.  Drainage basins of the Bayi Gully

下载: 全尺寸图片 幻灯片

“5·12”汶川地震前，八一沟泥石流已经是一条老泥石流，曾经发生过3次规模较大的泥石流(沈军辉等, 2008).2008年地震后，其在当年5月14日、5月19日、2009年7月17日、2010年8月13日、8月18日受强降雨影响再次暴发泥石流灾害，其中2010年8月13日这一次泥石流规模最大，危害程度最严重，见图 2，本文以八一沟“8·13”泥石流为例，采用FLO-2D进行数值模拟.通过2011年4月26日0.5 m全色波段WorldView-2影像遥感解译，准确统计出八一沟泥石流内崩滑体的面积，按面积大小分为3种类型：①小型滑坡体：滑坡体面积 < 10 000 m²，主要以浅层滑坡为主，厚度一般不超过5 m，总面积为90.856万m²，体积约为540万m³，主要分布在泥石流沟道两岸；②中型滑坡体：10 000 m² < 滑坡体面积 < 50 000 m², 主要以中等规模的滑坡为主，其厚度为5~50 m，总面积为115.803万m²，体积约为1 031万m³，一般发育于海拔较高，坡度较为陡峭的地方；③大型滑坡体：滑坡体面积 > 50 000 m², 主要以巨型滑坡为主，其厚度可超过50 m，总面积为44.656万m²，体积约为408万m³，发育于泥石流源地，地形陡峻的地方，能够为泥石流提供充足的固体物质.根据影像解译得到的滑坡体面积及分布特征，确定数值模拟的启动及边界条件，在滑坡体相对集中且沟道陡急的地方设置集水点位置，根据流域范围设置泥石流模拟运动时的路径边界条件，通过野外获取岩土体样品实验，计算泥石流体积浓度等参数，在主沟及支沟具有明显泥痕的地方建立典型断面，量测出泥石流各沟道的深度.

图  2  八一沟泥石流堆积扇对比

a. “5·12”地震后八一沟堆积扇全貌; b. “8·13”泥石流后八一沟堆积扇全貌

Fig.  2.  Comparison chart of Bayi gully fan

下载: 全尺寸图片 幻灯片

研究区属四川盆地亚热带湿润季风区，四季较为分明，温度偏低，主要以阴雨天气为主.每年降水量各个月份分配不均，5—9月约占全年降水量的80%，8月份为降雨最多的月份，降水量可达289.9 mm，多年平均降水量为1 134.8 mm.2010年8月13日14:00八一沟流域降雨幅度开始增加，在16:00-17:00间小时降雨量已经达到20年一遇的降雨量75 mm，12:00—23:00时间段内累计雨量达到160 mm（图 3）.

图  3  “8·13”泥石流降雨数据

Fig.  3.  Rainfall data of meteorological station on Aug. 13. 2010

下载: 全尺寸图片 幻灯片

2.   八一沟泥石流数值模拟

2.1   FLO-2D模型

八一沟泥石流数值模拟采用基于非牛顿流体及中央有限差分的FLO-2D软件进行，通过连续方程和运动方程FLO-2D控制泥石流在二维坐标系X方向和Y方向上流速和堆积深度(Mamodu and Dinand, 2013).FLO-2D模型用动力波模式及扩散波模式进行泥石流模拟，通过连续方程控制泥石流质量守恒，见式(1)，运动方程控制泥石流运动时的动量平衡，见式（2）、式（3）：

式中：I为一次降雨强度；t为时间；u为泥石流在二维模型中X方向的平均流动速度，h为泥石流的深度；v为泥石流在Y方向的平均流动速度；g为重力加速度，S_fx和S_fy为运动过程中的摩擦力，S_ox和S_oy为沟道纵坡降.基于拉格朗日描述的控制方程(质量和动量守恒方程)，可以将泥石流流体动态过程转化为一组偏微分方程(Huang et al., 2015).总应力方程包括以下两个部分：

式中：P代表了各向同性压力；τ^αβ代表了粘性应力张量；δ^αβ代表了总的应力张量；α、β代表了不同的坐标方向.FLO-2D模型具备下列假设条件：①保证模拟网格点具备唯一的高程值和粗糙度；②模拟流体运动过程属于稳态流；③流体运动属于浅水波模式及稳定流体的阻滞方程；④流体内压力的分布属于静水压力.FLO-2D数值模拟时，通过在极小的时间步长内快速计算保证数据的稳定性，时间步长主要是结合深度百分比变化参数（DPCP）和数据稳定性（CFL/FWPS）准则，当最小时间步长不满足3次连续时间间隔时，可以降低时间步长尽最大努力用合理的时间开展泥石流的模拟.

2.2   参数选取

本文主要采用5 m分辨率的DEM数据，将其转换成ASCⅡ格式的文件，用来模拟八一沟泥石流地形条件.根据现场实地调查确定有5条主要支沟，因此5个集水点的位置被确定，见图 1，集水点主要是指开始模拟泥石流的地方，同时也是实际调查过程中发现泥石流启动的地方.可以运用不同期高精度遥感影像及三维激光扫描仪对启动区物源展开动态观测，统计分析得到物源的动态演化规律，以此确定集水点所对应形成区内物源启动的难易程度.暴雨过程形成的地表径流导致坡体表面和前缘松散物质向下输移，汇入沟道，演变成危害性较大的泥石流，泥石流的深度以及支沟径流深度在本文中主要是通过野外现场调查获得的，根据主沟及支沟具有明显泥痕的地方建立典型断面，可以量测出泥石流主沟及各支沟的泥深，进一步确定数值模拟时各支沟中的参数.

2010年8月13日下午16时随着降雨幅度的增加，八一沟泥石流暴发，其中16:00至17:00一小时降雨量可达到75 mm，相当于研究区20年一遇的降雨量，泥石流数值模拟危险性评价不仅要考虑暴发强度，还要考虑到泥石流的暴发周期.结合往年的降雨数据及《四川省中小流域暴雨洪水计算手册》记录的降雨资料，重点研究八一沟泥石流流域中降雨强度为5年一遇、20年一遇、50年一遇、100年一遇、200年一遇5种情况下的降雨强度.通过积分的数学思想，由降雨首先推导出清水的总流量，即结合降雨数据和流域面积确定不同降雨强度下的洪水历时及八一沟流域内每一个集水区域内总的清水流量，具体见式(5)，通过《泥石流灾害防治工程勘查规范》以及《四川省中小流域暴雨洪水计算手册》提供的式(6)计算不同降雨强度下每个集水区域的泥石流最大洪峰流：

式中：W_P代表了集水流域内洪水总量（10⁴ m³）；H_Tp代表总时间为T的降雨总量（mm）；F代表了一个集水点所在的流域面积（m²）；H代表了支沟径流深度（cm）；Q_P代表频率为P的最大洪峰流量，（m³/s）；ψ代表了径流系数；i代表了平均最大暴雨强度；S代表了雨力系数；N代表了一次降雨指数；τ代表了降雨汇流时间（h）.

根据W_P，Q_P再结合《四川省中小流域暴雨洪水计算手册》中洪水流量过程线计算式(7)、式(8)，从而得到八一沟泥石流一个集水点流域范围内的清水流量过程线.

式中：T_P代表了清水流量过程线洪峰历时；T代表了清水流量过程线中不同时刻值（h）；Q代表了清水流量过程线中各个时刻相对应的流量值（m³/s）；x、y代表了四川省小流域洪水计算手册中清水流量过程线系数值.八一沟流域内共设置5个集水点，对应有5条流量过程线，而不同的降雨强度导致流量过程线不同，模拟开始后将同时进行.由于各支沟流域内的地形地势都有区别，因此同时进行数值模拟时，每个集水点的模拟参数及流量过程线都不同，通过野外测量、室内试验及上述公式求得，以③号集水点为例计算了它在降雨强度为5年一遇、20年一遇、50年一遇、100年一遇、200年一遇下的清水流量过程线，具体见图 4.

图  4  八一沟③号集水点不同降雨强度下清水流量过程线

Fig.  4.  The water flow process of ③ set point in Bayi Basin under different rainfall frequency

下载: 全尺寸图片 幻灯片

泥石流是由水沙混合的固液两相流，这些固体颗粒在流动过程中会呈现出不均匀分布的状态，体积浓度决定了泥石流的破坏程度，也决定数值模拟时泥石流流量过程线的基本参数，一般固体物质体积与泥石流总体积比值为体积浓度，具体参见式（9)：

式中：C_V代表了泥石流的体积浓度；V_W代表了固体物质的体积含量；V_D代表了泥石流总的体积含量.一般来说，泥石流的体积浓度与沉积物的含沙量、沉积物密度、泥石流混合物的密度以及大块颗粒含量密切相关.由于泥石流暴发时难以直接取样测试，根据马煜等对八一沟泥石流开展的物理实验研究，计算得到在堆积区上游容重为2.01 g/cm³，堆积区下游容重为1.88 g/cm³，结合数值模拟集水点位置位于堆积区的上游，推导其在20年一遇降雨下的体积浓度为60% (马煜等, 2011).在FLO-2D数值模拟中需要得到泥石流流量过程线，因此需要在洪水流量过程线的基础上乘以流量放大系数BF，由式(10)求得流量放大系数.此外，泥石流的屈服强度及粘滞参数也与模拟最终的堆积范围相关，具体参见式（11)：

式中：τ_y代表了屈服强度；η代表了粘滞参数；α₁，α₂，β₁，β₂代表了求得泥石流屈服强度及粘滞参数的系数值，主要通过实验或者参考FLO-2D手册获取，八一沟数值模拟参数主要根据詹钱登等进行一系列实验后得到的不同体积浓度下泥石流对应各种流变系数进行取值(詹钱登等, 2015)，α₁=0.811，α₂=0.004 62，β₁=13.7，β₂=11.24.

在模拟泥石流运动过程中，曼宁系数反应流域地表粗糙程度对水流运动的影响，根据Quan et al. (2013)分析数值模拟试验参数值选取范围，将该八一沟的曼宁系数设定为0.05.根据以上公式得到八一沟泥石流数值模拟的各项参数，以③号集水点20年一遇为例，具体见表 1.

表  1  八一沟泥石流③号集水点20年一遇参数选取

Table  Supplementary Table   The parameter of ③ set point in Bayi basin under 20 years rainfall frequency

参数	数值
流域面积F（km2）	2.08
沟道长度L（km）	2.5
径流系数ψ	0.918
洪峰流量QP（m3/s）	37.46
洪水流量WP（104m3）	5.62
体积浓度CV	0.6
径流深度H(cm)	27
汇流时间τ(h)	1.2
曼宁系数	0.05
屈服强度τy(MPA)	4 903
粘滞系数η	5 704
放大系数BF	2.5

下载: 导出CSV 

| 显示表格

2.3   数值模拟及验证

根据八一沟泥石流各个集水点的参数取值，模拟了在不同降雨强度下泥石流的运动情况，得到了不同频率下的泥石流堆积范围，见图 5，由图可以看出八一沟泥石流在5年一遇的降雨情况下松散固体物质不会冲出沟口，20年一遇及以上降雨强度可以冲出沟口掩埋农田、房屋等土地利用类型，特别是100年及200年一遇的降雨强度下可引发泥石流堵江现象，产生灾害链效应.

图  5  八一沟泥石流不同降雨强度下数值模拟结果及验证

Fig.  5.  Numerical simulation results of Bayi debris flow under different rainfall frequency

下载: 全尺寸图片 幻灯片

根据野外现场调查及覃怡等计算得到引发“8·13”八一沟泥石流的激发雨量相当于20年一遇(覃怡等, 2013)，因此通过数值模拟获取的20年一遇“8·13”暴雨诱发八一沟泥石流堆积危险范围、野外实际测量堆积危险范围及两者重合的危险范围验证八一沟泥石流数值模拟的精确程度，具体参见式（12):

式中：A为八一沟泥石流数值模拟的精确程度；S_a为野外实地测量八一沟泥石流堆积危险范围；S_o为八一沟泥石流实际堆积危险范围与数值模拟堆积危险范围的重叠面积；S_n为八一沟泥石流数值模拟后的堆积危险范围，具体参数见表 2. A的取值范围在0~1之间，当其取值越近1时，越能够表明数值模拟结果精确，经过计算八一沟泥石流模拟的精度A为77%，符合数值模拟的精度要求，见图 5, 由此得到采用FLO-2D模型进行泥石流数值模拟能够反应八一沟泥石流的运动特征.

表  2  八一沟泥石流数值模拟精度表

Table  Supplementary Table   Accuracy of numerical simulation in Bayi gully

泥石流名称	“8·13”泥石流堆积面积
	实测值(103 m2)	模拟值(103 m2)	重叠区(103 m2)	精确程度
八一沟泥石流	136.48	112.28	109.09	0.78

下载: 导出CSV 

| 显示表格

3.   八一沟泥石流危险性评价

本文结合泥石流强度及降雨强度开展八一沟泥石流危险性评价.通过FLO-2D模拟泥石流运动过程能够得到不同降雨强度下的八一沟泥石流流速及堆积深度，参照台湾学者彭思顯等人对泥石流强度的研究(Peng and Lu, 2012)，采用5年一遇、20年一遇、50年一遇、100年一遇和200年一遇降雨强度下的八一沟泥石流堆积深度及流速分布情况制定泥石流强度划分标准.研究表明泥石流强度大小主要体现在流速及堆积深度两方面(常鸣和唐川, 2014)，因为流速和堆积深度的大小直接决定了泥石流的破坏能力，因此根据泥石流的堆积深度及堆积深度与流速的乘积进行强度高低的界定，具体见表 3.

表  3  八一沟泥石流强度划分表

Table  Supplementary Table   Debris flow intensity partition in Bayi gully

八一沟泥石流强度	堆积深度(m)	关系式	堆积深度与流速乘积(m2/s)
高	H≥2.5	OR	VH≥2.5
中	0.5≤H＜2.5	AND	0.5≤VH＜2.5
低	0.0≤H＜0.5	AND	VH＜0.5

下载: 导出CSV 

| 显示表格

通过模拟八一沟泥石流不同降雨强度下的堆积危险范围，结合其实际暴发情况，对泥石流暴发概率进行划分，降雨强度20年一遇以内都为高暴发概率，20年至100年一遇为中暴发概率，100年至200年一遇为低暴发概率，高于200年一遇泥石流暴发概率非常低.将泥石流强度和暴发概率标准进行综合分析，得到泥石流危险性评价标准，见图 6.

图  6  八一沟泥石流危险性评价标准

Fig.  6.  Debris flow hazard assessment standards in Bayi gully

下载: 全尺寸图片 幻灯片

根据泥石流危险性评价标准，将八一沟泥石流的危险性划分为高危险性，中危险性及低危险性3类.本文详细统计了八一沟流域内这3种类型的面积及特征，具体见图 7.

图  7  八一沟泥石流危险性评价

Fig.  7.  Debris flow Hazard assessment result in Bayi gully

下载: 全尺寸图片 幻灯片

① 高危险性，该类面积占62%，由于该危险类型主要分布在沟道两侧及堆积区平坦的地方，因此大量的居民房屋、农田、道路及河堤分布在该类型中，极易受到泥石流的破坏，需要重点防护，必要时应针对一些大型重点居民聚集区、重要交通设施等保护对象进行工程防治，及时展开实时监测预警，确保遇到泥石流灾害时及时疏散当地居民.

② 中危险性，该类面积占28%，该类型所占研究区面积不大，主要有农田、少量居民点等土地利用类型分布，它们较易受到泥石流的破坏，主要体现为淤埋现象的中等程度损毁，需要对该类型范围内的房屋道路等土地利用类型进行重点防护，也要对重点部位展开相应的监测预警.

③ 低危险性，该类面积占10%，该类型在研究区中占极少的面积，主要分布有农田及河堤等土地利用类型，在泥石流暴发时该类型范围内的土地利用类型可能受到轻微破坏，但一般不会造成太大的经济财产损失，必要时可以适当增加预警措施.

4.   结论

受“5·12”汶川地震影响，强震区产生了大量松散固体物质，一旦遇到强降雨极有可能转换成泥石流灾害，给当地居民带来重大经济财产损失.为了详细研究强震区不同降雨频率下的泥石流危险性，便于及时开展监测预警，本文以都江堰龙池场镇八一沟泥石流为例，结合5年、20年、50年、100年及200年一遇的降雨强度，采用FLO-2D数值模拟的方法建立泥石流危险性评价模型.在降雨强度为20年一遇下的八一沟泥石流灾害事件中，通过野外实测数据验证，数值模拟的精度高达78%.本文充分考虑到泥石流暴发的强度效应及时间效应，利用数值模拟得到八一沟泥石流运动过程中的堆积深度和流速，以此建立泥石流强度划分标准，再结合不同降雨强度下的泥石流暴发概率最终建立强震区泥石流危险评价标准.本文主要针对泥石流的运动过程展开数值模拟，在今后的探索中将进一步加强对形成区泥石流启动机理及危险性的研究.

根据建立的强震区泥石流危险评价标准，绘制八一沟泥石流危险性分布图.结果表明，八一沟泥石流危险范围内属于高危险性区域占62%，中危险性区域占28%，低危险区域占10%，在中、高危险范围内房屋、道路、河堤及农田等土地利用类型很容易受到泥石流的破坏，需要加强监测预警及重点防护，当泥石流来临时确保及时疏散这些地区的居民，该危险分布图能够为八一沟泥石流灾后重建及重点土地利用类型的监测预警提供有效的参考依据.