基于HEC-HMS模型的不同雨型泥石流流量变化特征

熊 江，唐 川，龚凌枫

(成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059)

5·12汶川特大地震引发近20万处滑坡[1]，造成大量松散物体堆积在坡面和沟道中，在强降雨条件下易暴发群发性泥石流。震后崩塌、滑坡地质灾害呈逐渐减弱趋势，而泥石流却进入高暴发期，持续时期长达5～10 a，影响时期长达30～40 a[2-4]。泥石流暴发主要条件为陡峭地形、丰富物源、充沛水源，而强降雨成为震区泥石流暴发最活跃影响因素。降雨频率、降雨强度、降雨历时、降雨类型对泥石流暴发频率和规模产生重要影响。

目前关于降雨对泥石流影响的研究工作主要有探讨不同雨型临界雨量对泥石流暴发影响[5]；基于临界雨量基础上，采用临界雨量线方法进行泥石流预测[5-6]；通过构建泥石流危险性评价模型对单沟泥石流不同降雨频率下危险性进行评价[7]；通过对空间雨量站重采样和分位数统计方法研究前期降雨对泥石流暴发影响[8]。另外还有学者通过软件模拟泥石流产流过程，探讨泥石流堆积深度和范围[9-10]。但是对于不同雨型条件下泥石流流量变化特征研究工作目前还没有开展。因此本文在搜集汶川震区几次典型泥石流降雨数据基础上，采用Matlab对震区雨型进行概化，并利用HEC-HMS 软件模拟高家沟流域不同雨型条件下泥石流产流过程，由此分析和探讨不同雨型条件下泥石流流量变化特征。

1 研究区概况

高家沟地处四川省汶川县银信乡北侧，海拔高度为1 040～2 845 m，流域面积3.53 km2。高家沟位于川西台陷和茂汶—丹巴地背斜交接部位，北东—南西向茂汶断裂带，地质构造复杂。其松散物源丰富。高家沟属于深切侵蚀沟谷地貌，形态呈“V”型，地形西高东低，相对高差约1 800 m；高家沟流域范围内出露地层主要为晋宁期第四期侵入岩(γo2(4))和晚更新统冰水堆积物(Q3)、第四系冲洪积层(Q4al+pl)、残坡积层(Q4el+dl)、崩坡积层(Q4dl+col)、泥石流堆积层(Q4sef)等。位于亚热带湿润季风气候区，属川西多雨中心区，暴雨频发，多年平均降水量932.6 mm，最大日降水量269.8 mm，降雨量主要集中于5—9月，占全年降雨量的78%[11]。丰富的松散物源、充沛的降雨、陡峭的地形为泥石流暴发提供充足条件，地震后高家沟流域已暴发数次泥石流，其中2011年“7.3”泥石流影响最大，且降雨数据较完整，因此选用此沟作为研究对象研究。

2 HEC-HMS模型

2.1 流域布置

HEC-HMS水文模型是由美国陆军工程兵团水文工程中心结合GIS研发的半分布式降雨—汇流模型，被广泛运用于洪水监测和预报工作中[12]。另外有少数学者将其运用到泥石流研究中，并取得较好的效果[13]。该模型主要由流域模块、气象模块、控制模块三个部分组成[9，13]。其中流域模块主要模拟降雨后植物截留、地表下渗、蒸发、产流、流域汇流等过程，每个阶段提供具体的计算方法和所需参数；气象模块模拟流域不同气候条件，包括太阳辐射、降雨量、降雨过程等；控制模块主要控制模拟时间和数据采集时间间隔[14]。本文利用ArcGIS软件基于高精度DEM进行水文分析，提取出研究区内各子流域，在此基础上，通过布置高家沟流域内子流域(subbasin)、汇流点(junction)、径流线(reach)、泻流点(sink)等要素完成流域模型布置(图1)，并对其所需参数进行测定和调查。由此对不同雨型条件下泥石流产流过程进行模拟，并探讨其流量变化特性。

图1 高家沟流域特征图Fig.1 Characteristics of the Gaojia gully

2.2 流域模块

2.2.1植被模型

不同植被类型其截留率、持水能力不同[15]。HEC-HMS水文模型提供了网格简单植被、动态植被、简单植被、无植被四种模型。其中网格简单植被模型充分考虑了不同植被类型截留能力差异；而动态植被模型考虑了植被在不同生长阶段截留能力和作物系数变化；简单植被模型代表简单植被林冠层。根据野外调查发现，高家沟流域地震后植被恢复主要是混合型杂木，没有明显的植被类型分带，且模拟时间较短，林冠层截留能力和作物系数随时间变化不大。所以本文选择简单植被模型用于高家沟植被截留模拟。所需参数为植被初始截留量、最大截留量、作物系数。本文模拟降雨历时为24 h，降雨总量为336.9 mm，初始截留量设置为0；根据植被截留量与降雨历时、降雨量之间关系(图2)[16]，最大截留量设置为12 mm；作物系数(Kc)是反应植物生产不同时期需水量与蒸发量比值，作物生长旺盛时期最大值约为1.0[17]，震后高家沟流域植被恢复正处于旺盛阶段，综合实地情况本次参数设置为0.9。

图2 植物截留总量与降雨历时、降雨量关系Fig.2 Relationship between the total plant interceptionand rainfall duration and rainfall

2.2.2表面模型

HEC-HMS水文模型充分考虑了松散堆积体性质对下渗影响，主要涉及参数为初始含水率和饱和时所需最大降雨量。通过野外取样测验，地表松散堆积体初始含水率为6%，达到饱和时所需降雨量为14 mm。

2.3 产流模型

产流模型采用SCS 曲线法。该模型具有结构简单，所需参数少的特点，且充分考虑土壤质地、土地利用方式、前期土壤含水量等对产流影响[13]。该模型只需要一个反应流域综合特性参数CN值。该模型建立在水量平衡方程和两个基本假定基础上：假定流域实际地表径流量(Q)与潜在最大径流量比值等于实际下渗量(F)与潜在蓄水量(S)比值；初损量(Ia)属于潜在蓄水量的范畴。原理如下：

P=Ia+F+Q

(1)

(2)

Ia=λ·S

(3)

(4)

Q是由区域性地理因子和气候因子决定的参数，其取值范围为0.1～0.3[18]，在HEC-HMS半分布式水文模型中该值通常设置为0.2，且能准确预报产流流量[19]。因此式(4)转化为：

(5)

(6)

CN值充分考虑土壤类型和植被覆盖情况，其取值范围为0～100。结合CN值查算表和相关学者研究成果[19-20]，在充分对比和考虑流域实际情况基础上，本次CN值取85。

2.4 汇流模型

汇流模型选择SCS单位水位线法，该方法主要通过净雨量、流域面积、洪峰时间汇流时间计算单位线洪峰流量[21-22]：

(7)

(8)

(9)

(10)

式中：QP——净雨/mm；

F——流域面积/km2；

R′——单位净雨量/mm；

tp——降雨峰值时间/h；

tc——汇流时间/h；

L——滞后时间/h；

l——河道长度/m；

y——平均坡降/%。

通过ArcGIS软件提取流域坡降、径流线长度等结合式(10)可以计算出降雨形成径流后滞留时间。

2.5 河道演化模型

河道演化模型选择LAG模型。该模型适用于径流路径短，演变时间不受降雨深度变化影响。该模型结构简单，只需要一个演变时间参数[14]。高家沟属于急陡型沟道，流域内各子流域径流路径短，纵比降大，河道演变时间随降雨深度变化影响不大，且子流域地表粗糙程度基本一致，所以文章河道演化模型选择LAG模型是可行的。滞时演变模型方程为式(11)，结合滞时演变法与扩散模拟法之间参数关系式(12)、(13)可以计算出演变时间参数[23]：

(11)

(12)

(13)

式中：I(t)——入流量/(m3·s-1)；

Q(t+τ)——出流量/(m3·s-1)；

K′——线性水库滞时/s；

τ——线性渠道滞时/s；

t——时间/s；

D——扩散系数。

2.6 雨型模型

周伟等将汶川震区泥石流激发雨型概括为快速激发型、中速激发型、慢速激发型三类[5]。在峰值型、均匀型、递增型、递减型等基础降雨模型中峰值型最符合三种激发雨型降雨过程[13]。通过搜集汶川震区典型泥石流降雨数据，采用詹氏修正法对雨场进行分割[5]。分割结果见图3a～c，与快速激发、中速激发、慢速激发三种激发雨型具有较好对应关系。通过概化结果发现采用三次方进行雨型概化结果最接近雨场分割结果，因此将雨型概化为三次峰值早到型、三次峰值型、三次峰值晚到型(图3d～f)。为了探究不同雨型下泥石流流量变化特性，降雨总量和降雨历时需保持一致。根据《四川省中小流域暴雨洪水计算手册》暴雨量等值线图、变异系数计算出高家沟地区50年一遇降雨条件下，24 h降雨量为336.9 mm[11]。因此本文降雨总量设置为336.9 mm，降雨时长为24 h。其中三次峰值早到型降雨峰值设置在降雨开始后6 h、三次峰值型设置在12 h、三次峰值晚到型设置在18 h。各雨型小时雨量计算公式，三次峰值早到型：

qi=0.106955i3(i=1,…,6)

qi=0.003961(25-i)3(i=7,…,24)(14)

三次峰值型：

qi=0.013396i3(i=1,…,12)

qi=0.013396(25-i)3(i=13,…,24)(15)

三次峰值晚到型：

qi=0.003961i3(i=1,2,…,18)

qi=0.106955(25-i)3(i=19,…,24)(16)

3 泥石流流量计算

目前关于泥石流流量计算有配方法、经验公式法、成因法、形态调查法、雨洪修正法等[25]。其中雨洪修正法认为泥石流流量是清水流量和固体物质流量叠加之和，且充分考虑了泥石流体积浓度和沟道堵溃系数[25]。汶川地震后，震区泥石流流域内由于崩塌、滑坡等地质灾害造成大量松散堆积体不均匀堆积在沟道中，形成堵塞，泥石流暴发后沿程流量、流速、重度不断发生变化且具有放大效应。因此在计算震区泥石流流量时，堵溃系数是必须考虑的重要因素。基于此，本文采用雨洪修正法计算高家沟泥石流[10，23]：

Qc=Qp+Qs

(17)

式中：Qc——泥石流流量/ (m3·s-1)；

Qp——清水流量/ (m3·s-1)；

Qs——固体物质流量/ (m3·s-1)；

Cv——泥石流体积浓度/%；

Dc——堵溃系数。

通过现场野外调查，高家沟流域共有13处堵塞点，按《泥石流灾害防治工程勘查规范》调查方法对各堵塞点进行调查，结果见表1。本文取其平均数作为整条泥石流流域堵塞系数参与泥石流流量计算，其值大小为2.14。

表1 高家沟流域堵塞系数与重度

体积浓度是影响泥石流流速、冲出范围、冲出量的重要参数，其大小与泥石流沟道内物质颗粒级配、沟道纵比降、植被发育、物源分布等相关[10]。该值计算一般采用现场配浆法先计算泥石流重度，然后进行泥石流体积浓度计算：

(18)

式中：Cv——泥石流体积浓度/%；

γ——泥石流重度/ (kN·m-3)；

γw——水的重度，取9.8 kN/m3；

γs——泥沙重度，一般取26.5 kN/m3[10]。

在各堵塞处进行配浆实验，由泥石流暴发时现场目击者认为达到暴发时的浓度为止，现场对配浆体质量、体积进行测量，获取泥石流重度，见表1，取其平均值参与计算，其值为19.4 kN/m3。

图3 降雨模型概化图Fig.3 Rainfall model and generalized graph(a)三次峰值早到型；(b)三次峰值型；(c)三次峰值晚到型

4 结果与分析

4.1 模拟结果

通过流域出水口处sink元件对各雨型产生清水流量进行监测，将各雨型模拟结果与雨洪法计算结果进行统分析，结果见表2；将监测数据绘制成流量—历时曲线(图4)。结果显示三次峰值早到型于5时开始产流，初始阶段流量增长缓慢，6时之后流量突变加快，8时左右呈近直线快速增长至峰值，于12时流量达到峰值，其值为33.5 m3/s，误差为-19.4%，峰值滞后6 h，峰值之后流量快速降低，在16时处流量降低速度变缓。三次峰值型于8.5时开始产流，经过缓慢增长后于15时突变加快，于18时达最大值，为41.5 m3/s，误差为-0.2%，峰值滞后时间6 h，峰值后流量持续快速降低至0。三次峰值晚到型于13时开始产流，经过5h缓慢增长后转为近直线快速增长，于次日凌晨达流量最大值为45.8 m3/s，误差为10%，峰值滞后6 h，峰值后流量以近直线方式持续降低至0。另外从图4看出，三次峰值晚到型峰值流量最大、其次为三次峰值型、最小为三次峰值早到型，即降雨峰值推迟，开始产流时间越靠后，峰值流量越大。

表2 沟口清水流量模拟结果

图4 沟口清水流量历时线Fig.4 Line of clean water flow in the outlet

4.2 不同雨型泥石流流量变化特性分析

4.2.1峰值流量变化分析

为了探究不同雨型下峰值流量变化，基于沟口清水流量模拟结果，结合式(17)、(18)计算泥石流流量，将计算结果整理形成泥石流流量—历时曲线(图5)。同时将泥石流峰值流量计算结果与雨洪法计算结果进行统计分析，结果见表3。结果显示三次早到型峰值流量为166.83 m3/s、误差为-25.6%；三次峰值型峰值流量为206.67 m3/s、误差为-7.8%；三次晚到型峰值流量为228.084 m3/s，误差为1.7%。该结果表明降雨峰值越靠后，泥石流峰值流量越大，越接近雨洪法计算结果。另外由图5看出，峰值流量随降雨峰值推迟不仅呈现出逐渐增加现象，并且其增长幅度逐渐变小。这是由于降雨峰值越靠前，前期累计降雨量较少，降雨开始后，峰值降雨迅速到来，使得降雨强度超过土壤下渗能力，土壤饱和带深度较浅，孔隙水压力较小，因此其暴发规模、流量都较小；相反随降雨峰值推迟，前期降雨强度越小，降水充分下渗进入松散堆积体，土壤饱和带深度越深，孔隙水压力越大，因此泥石流暴发后规模、揭底深度越大，从而导致峰值流量逐渐变大，且其增长幅度逐渐变小。

图5 泥石流流量历时曲线Fig.5 Flow duration curve of debris flow

编号类型峰值流量/(m3·s-1)模拟法雨洪法误差/%1三次峰值早到型166.83-25.62三次峰值型206.67224.27-7.83三次峰值晚到型228.0841.7

4.2.2流量演变分析

为了探究整个流域从上游至沟口流量演变特征，将高家沟流域内各子流域出水口及沟口峰值流量整理和拟合分析，结果见表4、图6。由拟合结果可知，三种雨型下泥石流流量从上游到沟口呈现出线性演变模式，其相关系数均大于0.97，说明该拟合结果接近现实情况。其次从拟合结果发现流量演变过程中同一雨型保持稳定速率演变；而不同雨型间，随降雨峰值推迟，演变速率不断增大，而增长幅度逐渐降低。分析认为这与高家沟急陡沟道特性有关，由于高家沟流域面积小、物源丰富、纵比降大、侵蚀强烈，因此在强降雨情况下，各子流域快速产流并迅速输出出水口，并无长时间滞留，因此从上游至沟口流量表现出线性增长。同时由于高家沟属于急陡型沟道，泥石流暴发物源条件和地形条件极其充分，其泥石流暴发和汇流特征主要受降雨条件影响和控制。根据前文分析，随降雨峰值推迟，前期降雨能充分下渗，松散堆积体充分饱和，孔隙水压力越大，泥石流暴发后其流量、规模也越大，因此泥石流流量随降雨峰值推迟其演变速率逐渐增加，而增长幅度逐渐降低。

表4 流量演变模式

图6 流域流量演变结果Fig.6 Results of the watershed flow evolution

4.2.3流量增长和衰退阶段分析

为了探究泥石流流量增长和衰退阶段特性，需要确定泥石流暴发时间。通过计算高家沟泥石流暴发临界流量，将其在泥石流历时曲线中对应时间作为泥石流暴发时间节点[26]：

(19)

(20)

(21)

式中:qc——临界流量/(m3·s-1)；

ρs——松散堆积体密度/(t·m-3)；

ρ——水密度/ (t·m-3)；

θ——沟床坡度/(°)。

根据高家沟流域调查和堆积体筛分试验获得颗粒级配曲线(图7)，堆积体密度ρs为1.97 t/m3。计算出高家沟流域泥石流暴发临界流量为8.36 m3/s，将该值绘制到泥石流流量历时曲线图5中。结果显示三次峰值早到型于7时左右暴发泥石流；三次峰值型于12时左右暴发泥石流；三次峰值晚到型于15时左右暴发泥石流。为了探究不同雨型下泥石流流量增长和衰退阶段特征，将泥石流暴发后流量—历时曲线划分为增长阶段和衰退阶段，采用最小二乘法进行拟合分析，结果见表4、图8。

图7 颗粒级配曲线Fig.7 Particle grading curve

结果显示流量增长阶段三种雨型均以三次函数形式增长，且相关性系数均在0.99以上。另外通过图8a可以看出，三次峰值早到型经过1 h缓慢增长后转

为快速增长，且其增长速率逐渐降低，从产流开始到峰值所需时间约5.5 h；三次峰值型经过约3 h缓慢增长后转为快速增长，到峰值所需时间6 h；三次峰值晚到型缓慢增长阶段时长约为1 h，到峰值所需时间约9 h。即三种雨型下，泥石流流量增长阶段所需时间随降雨峰值推迟而增大。

流量衰退阶段，三种雨型流量都以对数形式衰退(图8b)，即初始阶段衰退速度较快，然后逐渐变缓。其中三次峰值早到型衰退阶段所需时间约为14.5 h；三次峰值型约为8 h；三次峰值晚到型约为2.5 h。即降雨峰值越靠后，其衰退所需时间越少，衰减速度越快。出现这种现象原因主要是由于降雨峰值越靠后，前期降雨小时雨量越小、累计雨量越大，使得降水充分下渗，土壤水含水量越高，松散堆积体更加饱和，孔隙水压力越大，因此泥石流暴发规模越大，流速越快。另外由于降雨总量和时间一致，因此峰值越靠后，峰值后累计雨量越小，导致泥石流衰退阶段速度更快。

表5 泥石流流量增长和衰退阶段拟合结果

图8 不同雨型流量增长和衰退阶段Fig.8 Stages of flow growth and decline under different rain types(a)流量增长过程；(b)流量衰退过程

5 结论

(1)三次峰值早到型、三次峰值型、三次峰值晚到型峰值清水流量分别为33.5 m3/s、41.5 m3/s、45.8 m3/s，误差分别为-19.4%、-0.2%、10%；泥石流峰值流量分别为166.83 m3/s、206.67 m3/s、228.08 m3/s，误差为-25.6%、-7.8%、1.7%。

(2)三种雨型下泥石流流量从上游到沟口以线性模式不断演变，且随降雨峰值推迟，演变速率不断增大，幅度逐渐减小。

(3)计算得出高家沟泥石流暴发临界流量为8.36 m3/s，推算出三次峰值早到型、三次峰值型、三次峰值晚到型三种雨型泥石流暴发时间分别为7时、12时、15时左右，即出现在峰值降雨前后，属于降雨激发型泥石流。

(4)泥石流流量增长阶段，三次峰值早到型以三次下凹型增长，而三次峰值型和三次峰值晚到型以三次上凹型增长。且达到峰值所需时间分别为5.5 h、6 h、9 h，即随降雨峰值推迟所需时间越长。

(5)流量衰退阶段，三种雨型均以对数上凹型衰退。衰退所需时间降雨峰值推迟而逐渐减少，分别为14.5 h、8 h、2.5 h。