考虑大孔隙下渗的HBV模型及其在山洪预报中的应用

崔逸凡，吕海深,2,3，朱永华,2，苏建宾

(1.河海大学 水文水资源学院，江苏 南京 210098；2.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；3 河海大学 全球变化与水循环国际合作联合实验室，江苏 南京 210098)

1 研究背景

近年来，由于全球气候变化导致极端天气频发，突发性暴雨频率快速增加，洪水灾害已成为造成人民生命财产损失的重大灾害，严重制约着广大地区社会经济的发展[1]。据统计，中国是一个洪水灾害频发的国家，洪水灾害造成的死亡人数占自然灾害死亡人数的比例呈上升趋势，对人民生命安全造成了严重危害[2]，因此洪水的防治已迫在眉睫。洪水防治主要包括工程措施和非工程措施，其中，洪水预报作为洪水防治最重要的非工程措施已上升到国家战略层面。目前，对于大江大河的洪水预报技术已经比较成熟，而受限于对山洪形成机制认识不全面、山区小流域水文资料匮乏和降水预报精度不高等因素，小流域山洪预报的精度往往并不能满足实际需要。水文模型模拟是进行洪水预报的重要方法，但是由于小流域山洪具有历时短强度大、陡涨陡落等特点，许多模型对于小流域山洪的模拟精度并不高，容易出现洪峰模拟较差以及洪水过程线模拟不准确等现象。

HBV模型是由瑞典国家水文气象局(SMHI)在20世纪70年代开发的水文预报模型[3]，目前已被广泛应用于全球40多个国家的数百个流域的洪水预报和水资源评估中。近年来，HBV模型在国内也已经得到广泛应用，如赵彦增等[4]利用HBV模型对淮河官寨流域进行了模拟研究，径流过程模拟结果比较理想，可在国内推广应用；王有恒等[5]利用HBV模型在白龙江流域进行了模拟研究，得出HBV模型能较好地模拟洪水对降水的响应过程。然而HBV模型对土壤层的处理较为简单，仅将土壤层分为上、下两层，其物理过程并不明确。在植被发育区，由于土壤胀缩、动物活动、植被发育及植被根系死亡等，大孔隙普遍存在。大孔隙虽然在土壤总孔隙中所占比例较小，但由大孔隙引发的优先流对降雨的快速入渗和土壤水的快速补给产生了重要贡献[6]，大孔隙的存在还会改变坡面尺度地表径流、壤中流和地下径流的形成过程和比例，提早地下径流的出流时间和增加地下径流在总径流中的比重[7]，从而对径流形成过程产生重要影响。

淮河流域位于我国南北气候过渡带，气候变化复杂，降水时空分布不均匀，淮河上游干流及其南部山区的河道坡度大，汇流时间短，极易发生短历时高强度的山洪[8]。同时，土壤大孔隙和大孔隙流在淮河流域普遍存在[9]，对径流形成过程产生重要影响。因此，本文以淮河流域上游的黄泥庄小流域为研究区，引入变动渗漏面积方法为HBV模型增添大孔隙模块，探究改进后的HBV模型在该流域山洪预报中的适用性。

2 研究区域及数据来源

2.1 研究区域

选择淮河上游黄泥庄水文站以上集水区(下文简称黄泥庄流域)为研究区，该研究区介于北纬31°06′～31°38′和东经115°21′～115°43′之间，流域控制面积805 km2。该流域地处我国南北气候过渡带，多年平均气温为11～16 ℃，流域多年平均年降水量为1 077 mm[10]，但降水的年内分布不均。流域年降水量的50%～80%集中在汛期(6-9月)[10]，且多以暴雨形式出现。汛期降水历时短、强度大、暴雨中心笼罩范围小，因此汛期时由暴雨引发的山洪频发。流域内植被发育良好，森林覆盖率高达65%，存在丰富的地下水和壤中流。流域平均海拔为479 m，但地形复杂，西北和东南以山地为主，河道坡度大，双扇形河网使区域产流快速向流域出口汇集。流域概况如图1所示。

图1 黄泥庄流域概况

2.2 数据来源

如图1所示，黄泥庄流域有相对密集的降水观测站网，流域及其周边的12个雨量站点的降水观测数据可以从淮河流域水文年鉴资料中获取。但是获取的长历时资料主要以日尺度观测为主，仅在必要时进行小时尺度加密观测。因此，该项研究仅选择了2010 -2015年间有小时尺度观测数据的特定洪水事件进行研究。考虑到12个雨量站点在黄泥庄流域分布相对均匀，在对雨量站的降水数据仔细检查后，使用泰森多边形面积权重方法获取小时尺度流域面平均雨量。另外，由于雨量站缺乏气温观测且流域范围内没有国家气象站，因此使用了高分辨率的欧洲中期天气预报中心再分析产品(ECMWF)的逐小时气温数据作为HBV模型的驱动，该产品已经被证明在中国区域有很好的应用效果[11-12]。黄泥庄水文站的径流观测资料来自淮河流域水文年鉴洪水要素摘录表。为了与模型输出结果进行匹配，研究中使用了线性插值方法对加密观测的洪水摘录资料进行插值处理，进而得到洪水事件期间逐小时的流量观测资料。

考虑到下游的防汛承受能力，研究中以500 m3/s为阈值，在2010-2015年间共选择了6次洪水事件(详细信息如表1所示)，且在选定的洪水事件期间，12个雨量站都有完善的加密降水观测资料。从表1中可以看出，2010-2015年期间每年都至少有1次较大洪水，且洪水事件主要集中在每年的6-7月间。其中，最严重的1次发生在2013年，洪峰流量高达1 628.49 m3/s。

表1 黄泥庄流域2010-2015年洪水事件概况

3 HBV模型改进

3.1 概述

HBV模型为概念性水文模型，它具有输入参数少、适用性强、模拟精度高等优点，可以适用于各种复杂气候条件下的水文模拟，目前已被广泛应用于水文预报、气候变化研究、水资源评估、地下水模拟和水量平衡研究等方面[13-14]。HBV模型以降水量、气温和潜在蒸散发为主要驱动，输出出口断面径流量。如图2所示，HBV模型主要由3个部分组成，分别是：积雪与融雪模块、土壤湿度计算模块和径流响应模块。考虑到大孔隙对HBV径流响应过程的影响，研究中使用变动渗漏面积方法对上层水库进行改进，详细描述如下。

图2 考虑大孔隙的HBV模型结构

3.2 积雪和融雪模块

HBV模型采用阈值方法进行雨雪分割(此次研究中使用0 ℃为雨雪分割阈值)，且模型假定融雪和积雪速率与温度成正比。当温度低于阈值温度时，降水主要以雪的形式累积为积雪，反之则为降雨。当温度一直保持在阈值温度以下时，则不会产生径流，但温度一旦超过阈值温度，融雪径流开始产生，其中融雪速率如下式：

Sm=DD(T-Tt)

(1)

式中：Sm为融雪速率，mm/h；DD为度日因子，mm/(℃·h)；T为小时平均气温，℃；Tt为阈值气温，℃。

3.3 土壤湿度计算模块

流域的降水通常分为两部分：一部分渗入土壤；另一部分用于形成地表径流，即为时段径流深ΔQ。一般来说，土壤含水量越高，降水对径流形成的贡献就越大，当土壤含水量接近最大土壤含水量时，下渗量减少，降水对产流量的贡献增加。时段径流深的计算如下：

(2)

式中：ΔQ为时段径流深，mm；SM为实际土壤含水量，mm；FC为最大土壤含水量，mm；β为形状系数。

由于蒸散发与气温息息相关，模型采用气温折算方法获取小时尺度潜在蒸散发。如公式(3)所示，利用小时平均气温与日平均气温之间的差值来修订日平均潜在蒸散发，进而获取小时尺度的潜在蒸散发。考虑到实际蒸散发受土壤含水量影响很大且与修正后的潜在蒸散发呈线性关系，因此当土壤含水量大于或等于土壤凋萎含水量时，实际蒸散发就等于修正后的潜在蒸散发；当土壤含水量小于土壤凋萎含水量时，受限于土壤水分不足，实际蒸散发将使用土壤含水量与凋萎含水量的比例进行折算。

PEa=(1+C(T-Tm))PEm

(3)

(4)

式中：PEa为修正后的潜在蒸散发，mm；Tm为日平均气温，℃；PEm为日平均潜在蒸散发，mm；C为模型参数；Ea为实际蒸散发，mm；PWP为土壤凋萎含水量，mm。

3.4 改进后的径流响应模块

HBV模型将径流形成过程概化为一个响应函数，该函数由一个上层非线性水库和一个下层线性水库组成，即把径流成分概化为两个响应盒子。上层盒子有两个径流出口，其响应函数通过消退系数K0、K1来体现，下层盒子有一个径流出口，其响应函数通过消退系数K2来体现。当上层盒子水位超过阈值L时，上层盒子迅速形成地表径流(Q0)，并从第一个出口流出。上层盒子的第二个出口和下层盒子出口的径流响应相对较慢，形成壤中流(Q1)和基流(Q2)。流域出口径流深(Q)即为各部分径流深之和。

(5)

(6)

Q=Q0+Q1+Q2

(7)

式中：Su为上层水库水位，mm；Sl为下层水库水位，mm；L为阈值水位，mm；K0为地表径流消退系数；K1为壤中流消退系数；K2为基流消退系数；Q0为地表径流，mm；Q1为壤中流，mm；Q2为基流，mm；Q为总径流深，mm。

但是，HBV模型中并未考虑大孔隙对降雨入渗的促进作用以及大孔隙的存在对径流组成成分的影响，所以模拟出的洪水过程线可能不准确，尤其是对洪峰的模拟精度较差。考虑到这一点，研究中将研究区划分为设置模拟大孔隙下渗的变动渗漏面积和一般透水面积两部分，且变动渗漏面积与土壤含水量有关。参考肖庆元等[9]的研究结果，变动渗漏面积和一般透水面积划分如下：

(8)

式中：AA为变动渗漏面积所占比例；BB为一般透水面积所占比例；IA为AA的上限；n为经验指数。

在增加大孔隙模块后，改进的HBV模型中，将上层盒子分为两个部分：变动渗漏部分和一般透水部分。在变动渗漏部分，下渗能力很强，一般情况下水可以沿着大孔隙直接进入下层盒子而无损失；在一般透水部分，仍然依据HBV模型的渗透常数方法进行，改进后的上层土壤渗流可以表述为：

Qprec=Su·Kprec+AA·ΔQ

(9)

式中：Kperc为渗流消退系数；Qprec为渗流，mm。

3.5 模型评价标准

根据《水文情报预报规范》(SL 250-2000)，采用Nash-Sutcliffe效率系数即Nash系数、洪峰流量相对误差、洪量相对误差和峰现时间误差来衡量模型的模拟效果，它们的计算公式和最优值如表2所示。Nash系数反映的是径流模拟值与观测值之间的吻合度，洪峰流量相对误差反映的是对洪峰的模拟效果，洪量相对误差反应的是对洪水总量的模拟效果，峰现时间误差反映的是对洪峰出现时间的模拟效果。

表2 模型评价指标

4 模拟结果和分析

4.1 参数敏感性分析

在改进的HBV模型中，IA表示变动渗漏区在流域中所占的最大比例；n为经验指数，反映随着土壤湿度变化，变动渗漏区中实际大孔隙区域的变化。由公式(8)可知，当IA为定值时，随着n的增加，流域中大孔隙实际所占的面积逐渐减少；反之，当n为定值时，随着IA的增加，实际大孔隙面积比例增加。当IA=0时，变动渗漏面积所占比例为0，可认为土壤中不存在大孔隙；当IA=1，n=0时，变动渗漏面积所占比例为1，可认为土壤中全是大孔隙。为分析模型改进效果，并确定新增大孔隙模块主要控制参数的敏感性，不失一般性的选择20100705次洪水为例，采用控制变量法分析参数IA和n对洪峰流量和洪水总量的影响，即使IA在0.1～1之间以0.1的步长均匀变化，n在0～0.8之间以0.2的步长均匀变化，统计IA和n的变化对洪峰流量和洪水总量的影响，结果如图3和4所示。

由图3和4可知，洪峰流量和洪水总量对新增加大孔隙模块的两个参数都很敏感。随着IA的增加和n的减小，流域中大孔隙所占面积逐渐增加，洪峰流量和洪水总量均随之减小且洪峰流量的变率远超洪水总量。这主要是由于流域孔隙度的增加促进了土壤的快速入渗进而产生了更多的壤中流，而壤中流的汇流时间较长使得模拟洪水过程线坦化，降低了洪峰的模拟值。深入分析可以发现，IA和n之间有着复杂的相互关联。当IA较小时，变动渗漏区域在流域中的占比较小，洪峰流量和洪水总量对n的变化响应强度较弱，随着IA的增加，曲线的分布更加分散，表明此时主要的控制参数为经验参数n。总的来说，对于IA而言，其高值区间比低值区间更加敏感，而对n而言则是低值区间更加敏感。

4.2 模拟结果分析

基于改进后的HBV模型，对2010-2015年间选定的6次洪水事件进行模拟，其中前4次洪水事件用来校准模型，采用试错法获取模型的最优参数[15]，后2次洪水事件用来检验模型结果。

表3汇总了6次洪水事件模拟结果的统计指标，总的来说，改进后的HBV模型在黄泥庄流域洪水模拟中模拟效果较好。率定期的4次洪水事件Nash系数分别为0.82、0.92、0.91和0.89；验证期的两次洪水事件Nash系数分别为0.92和0.76，能较好地匹配两个时期洪水过程。除明显低估了20130705次洪水的洪峰流量外(洪峰流量相对误差为-28.38 %)，其余洪峰流量和洪量相对误差均在20 %以下，能较好地捕获洪峰流量和洪量，对灾害预警和灾情评估有重要意义。此次率定期和验证期的峰现时间误差均在2 h以内，改进后的HBV模型能够准确地再现洪峰出现时间。依据《水文情报预报规范》(SL 250-2000)，6次洪水事件的模拟精度均达乙级及以上标准，满足实际工程应用精度，可用于该流域洪水模拟与预报。

图3IA-n-洪峰流量相关图 图4IA-n-洪水总量相关图

表3 黄泥庄流域洪水模拟结果

图5和6分别展示了率定期和验证期洪水过程线的观测值和改进后HBV模型的模拟结果。如图1所示，黄泥庄流域由南北两个扇形子流域构成，洪水涨落迅速，主要以单峰的尖细洪水过程为主(图5和6)。尽管改进后的HBV模型能较好地模拟洪水过程线，但依然倾向于低估洪峰。这可能是受限于集总式模型结构特征，改进后的HBV模型依然难以刻画径流过程对流域降水空间异质性的响应。同时，黄泥庄流域的雨量站点均分布在靠近河道的相对平缓的谷地，对坡面地形雨及短历时暴雨的捕获率不足，难以准确反映洪峰对应期流域降水，这也进一步加剧了HBV模型对洪峰流量的低估。另外，可能是受到前期土壤湿度的影响，黄泥庄流域水文过程线对小雨并不敏感，这也导致HBV模型高估了洪水过程的起涨点(例如20100705和20130705洪水)，后续的研究需要针对这一过程进行进一步的改进。集总式HBV因其具有模型结构相对简单、计算量小、模型参数和所需驱动数据较少等优点，已经被广泛应用于洪水模拟与预报，但对降雨径流过程概化不足也进一步增加了模型的不确定性。总体而言，添加对物理过程的精细化描述能进一步提高模型的模拟精度。

图5 模型率定期观测和模拟的洪水过程线

图6 模型验证期观测和模拟的洪水过程线

5 结 论

无资料小流域的山洪预报一直是洪水预报中的一个难点，土壤大孔隙及其引发的大孔隙流通过改变流域产汇流过程，进一步增加了这类地区洪水模拟与预报的难度。本文使用变动渗漏面积法改进HBV模型，并以淮河流域上游黄泥庄小流域为例，对2010- 2015年的6次洪水事件进行模拟，得出了以下结论：

(1)大孔隙模块的控制参数IA和n的变化会对洪峰流量和洪水总量产生较大影响。

(2)率定期的Nash系数在0.82～0.92之间，验证期的Nash系数在0.76～0.92之间，洪峰流量相对误差和洪量相对误差基本控制在20 %以内，峰现时间误差不超过2 h。3次洪水事件的模拟精度达到洪水预报甲级标准，3次达到乙级标准，说明改进后的HBV模型能够很好地模拟黄泥庄流域的洪水过程，满足该流域的实际应用需求。

(3)本研究可为淮河流域山区的防洪减灾、水资源管理调度等提供参考。同时，引入变动面积方法对HBV模型的改进和对存在大孔隙的相似小流域山洪预报研究有重要的借鉴意义。