考虑动态汇流路径的时变分布式单位线

易 彬,陈 璐

(1.华中科技大学土木与水利工程学院,湖北 武汉 430074；2.华中科技大学数字流域科学与技术湖北省重点实验室,湖北 武汉 430074)

流域汇流是地面径流、壤中流、地下径流汇集到流域出口断面的径流过程。净雨滴通过不同介质、不同路径、不同粗糙面向流域出口断面移动,实际汇流过程十分复杂。过去100 a间,大量学者从水动力、水文、水滴运动以及系统视角观察和分析流域汇流过程[1]。1932年,Sherman[2]在研究如何用降水量推求径流量的过程中,提出了单位线的概念,目前该方法已成为使用最为广泛的汇流计算方法。

单位线将汇流过程概化为线性时不变系统,实际上,汇流特性不仅受控于流域地形地貌和下垫面特征,也与净雨时空分布等因素有关[3]。基于流域地貌参数,Rodríguez-Iturbe等[4]提出了地貌单位线方法,但该类方法无法考虑流域地形及下垫面特征的空间分布；为了考虑下垫面空间异质性对汇流的影响,Maident等[5]提出了基于数字高程模型的分布式汇流单位线。分布式单位线计算方法将水质点在流域中的滞留时间等价为流域瞬时单位线,其核心在于速度场的计算,常用的坡面流速计算公式有SCS流速公式、曼宁公式、均匀流方程、达西公式等[6],上述流速公式考虑了植被、坡度、土地利用等因素的影响[7-8],但忽略了由时变雨强、上游入流贡献等因素引起的径流汇流非线性问题。

有学者通过建立坡面流速与时段净雨强度、上游入流贡献等因素的定量映射关系,提出了时变分布式单位线法[9-10]。孔凡哲等[11]基于改进SCS流速计算公式提出了考虑降水强度的时变分布式单位线计算方法；Bunster等[12]考虑动态的上游入流贡献,开发了时变动态上游入流贡献的分布式单位线,提高了洪峰预报精度。然而,上述方法均假定全流域水质点存在静态汇流路径,即流域上均匀分布的净雨滴通过唯一的汇流路径向出口流动,事实上,以蓄满产流为主的南方湿润地区,在流域未达到蓄满状态之前,只有蓄满区域产生净雨,其汇流路径应随流域蓄水状态而动态变化,当流域蓄满后,动态汇流路径转变为静态汇流路径。

为此,本文考虑下垫面空间分布异质性和降水强度的时变特性,根据各个时段净雨量级大小、下垫面土壤含水量状态推求流域时变汇流路径,尝试根据流域产流情况将流域静态汇流路径转换为动态汇流路径,推求同时考虑降水强度和下垫面土壤含水量时空分布特征的坡面流速计算公式,提出一种考虑动态汇流路径的时变分布式单位线。

1 研究方法

现有时变分布式单位线存在2个假定：全流域蓄水量在降水结束前可以达到饱和状态以及全流域的坡面、河网汇流路径唯一。本文从流速计算和汇流路径两方面提出解决思路,一方面,根据土壤蓄水状态提取流域动态汇流路径,另一方面,通过考虑时变土壤含水量改进现有坡面流速计算公式,在此基础上,结合改进时变坡面流速计算公式和动态汇流路径推求考虑动态汇流路径的时变分布式单位线。

1.1 产流计算

SCS模型由美国农业部水土保持局于1954年开发,特点是结构简单、参数较少,被广泛用于小型集水区径流预报[13],本文采用该模型进行时段净雨(Rt)计算。

1.2 动态汇流路径

本文将动态汇流路径定义为蓄满栅格流向流域出口过程中所经过的路径,由于降水过程的持续,流域蓄水状态实时变化导致汇流路径动态变化,获取任一时刻流域汇流路径的关键在于确定流域当前蓄水状态分布情况,新安江模型常用帕累托分布表征土壤蓄水能力的空间分布不均匀性,其缺陷在于无法获得流域蓄水能力在空间上的具体分布情况[14]。研究表明,地形指数与流域蓄水能力密切相关[15],本文采用地形指数(TI)表征流域土壤蓄水能力,根据地形指数计算栅格蓄水能力[16],表达式如下：

(1)

式中：WM,j为j栅格单元的最大蓄水容量；WM,min为全流域所有栅格中的最小蓄水容量；WM,max为全流域所有栅格中的最大蓄水容量；TI,min=min{TI,j,j=1,2,…,N}，N为流域栅格单元总数；TI,max=max{TI,j,j=1,2,…,N}；n为待率定参数,反映蓄水容量不均匀分布的系数,该值越大,流域蓄水容量的空间分布越不均匀。

(2)

式中：Fαt为流域蓄满部分的面积；Pi为i时段降水量；Ri为i时段净雨量；Wj,1为j栅格初始时段蓄水量；Aj为栅格单元面积。

采用动态汇流路径进行流域汇流计算时,由于仅蓄满栅格发生产流,而SCS模型认为全流域均匀产流,因此,t时段净雨量(Rt)需重新分配至蓄满面积,根据水量平衡原理：

(3)

(4)

1.3 基于动态汇流路径的时变分布式单位线

分布式单位线方法已较为成熟,被广泛用于中小流域汇流计算,具体流程详见文献[11],本文所提方法步骤如下：

(1) 根据流域DEMs将流域栅格化,确定栅格边长,在此基础上,采用D8算法确定栅格水流方向,并根据栅格水流方向提取该栅格至流域出口的汇流路径集合Rroad,j(第j个栅格至流域出口路径上的栅格单元集合),则全流域所有栅格汇流路径集合可表示为

Rroad={Rroad,j|j=1,2,…,N}

(5)

由于该路径集合是静态的,然而实际汇流情况为只有蓄满区域有汇流路径,为此,结合分布式蓄水容量,推求不同蓄满比例(αt)下蓄满栅格的汇流路径集合(Rroad,αt),满足：

Rroad,αt={Rroad,j|Aj∈Fαt}

(6)

式中：Rroad,αt为流域蓄满区域的汇流路径集合,Rroad,αt⊂Rroad。随着降水过程的持续,αt不断变化,因此,Rroad,αt构成流域不同蓄水状态下的动态汇流路径集合。

(2) 计算每个栅格逐时段的汇流速度(Vj,αt),获得不同流域土壤含水量状态下的流速分布场。现有分布式单位线仅考虑降水时空分布不均匀性和下垫面空间分布异质性对流域汇流速度的影响,其假定流域在降水脉冲结束之前可达到蓄满状态,这与实际情况不符,忽略了时变土壤含水量造成下垫面状态的改变,已有研究表明,坡面流速与植被覆盖、雨强、坡度、土壤含水状态有关[7-8,17]。为此,针对蓄满单元,研究采用考虑雨强的SCS流速计算公式[11],针对处于汇流路径的未蓄满栅格单元,提出考虑降水强度和土壤含水量的时变流速计算公式：

(7)

式中：Vj,αt为当流域蓄满比例为αt时j栅格单元的汇流速度；k为流速系数,依据下垫面植被、土壤类型进行确定；βj为j栅格坡度；Ic为流域参考雨强；Wj,t为j栅格t时段蓄水量。

由式(7)可知,当某一栅格属于河道单元时,流速固定为2 m/s[6]；当该栅格处于未蓄满区域且未处于蓄满栅格汇流路径上时,栅格流速为0,不参与汇流；当该栅格处于蓄满区域时,流速跟时变雨强有关；当该栅格处于未蓄满区域,但上游存在蓄满区域时,需要同时考虑雨强和土壤含水量对栅格流速的影响。

其中,未蓄满栅格的土壤含水量表达式如下：

(8)

式中：Wj,t-1为j栅格t-1时段蓄水量；Qj,t为j栅格t时段的上游入流贡献量;Pt为t时段降水量。

(9)

(10)

式中：Lj为栅格边长。

(11)

在此基础上,将蓄满栅格面积转换为流域出口流量过程,得到任意时段的分布式单位线：

(12)

式中：Qm为第m个时段的单位线流量,m3/s；H为单位净雨,mm。

1.4 参数率定与目标函数

研究全面考虑洪水过程低流量误差、高流量误差、洪量误差等多目标因素,其中,高流量误差强调洪峰预报误差[18]。选择Kling-Gupta系数(EKG)[19]、纳什效率系数(ENS)[20]、纳什效率系数对数形式及均方根误差与标准差的比值(RSR)4个目标参与优化率定,优化算法采用SCE-UA算法[21]。其中,EKG和ENS包含的平方误差项强调了高流量的预报精度；RSR系数通过均方根误差与标准差的比值量化了洪量误差；ENS通过取对数压缩变量尺度,强调了低流量预报精度。为全面考虑4个目标,建立了涵盖上述指标的综合目标函数,依据Brunner等[22]研究成果,采用多目标加权的方式对水文模型进行优化率定,表达式为

O=w1(1-ENS)+w2(1-EKG)+w3lg(1-ENS)+w4RSR

(13)

式中：w1、w2、w3、w4分别为4个目标权重。本文重点关注汇流精度,洪峰准确度较为重要,因此,ENS和EKG的比重相对较高,RSR和log(1-ENS)的比重相对较小。

2 例证研究

2.1 研究区域及数据

龙虎圩水和东石河是中国南方山区小河流,位于广东省韩江流域,分别汇入程江和柚树河,地处山区丘陵地带,地形复杂,属亚热带气候,气候高温湿热、暴雨频繁,但时空分布不均,具有河床坡降陡、天然落差大、洪水易涨易退等特点,降水径流存在显著非线性关系,容易引起山洪灾害。为此,研究选择以上2个小流域为研究对象,其中,龙虎圩水河长20.8 km,流域面积130.8 km2,河流平均比降2.92‰；东石河全长23.6 km,流域面积152.4 km2,河流平均比降3.56‰。流域DEMs分辨率为30 m×30 m,在此基础上提取流域栅格坡度、流向等信息;流域植被覆盖数据来自清华大学植被数据库(data.ess.tsinghua.edu.cn)；用于参数率定和检验的28场降水径流数据由广东省水文局梅州水文分局提供,其中,20场用于参数率定,8场用于结果检验。

2.2 参数率定结果

选取2个流域近年各10场洪水进行参数率定。模型率定分为产流参数率定、单位线参数率定、汇流参数率定3个阶段。需要率定的参数有初损系数(η)、降雨前期流域特征参数(CN)、Ic、表征蓄水容量空间分布均匀性的系数(n)、WM，min及WM，max,其余参数(流域坡度、流速系数)可根据流域下垫面信息确定。

(1) 根据历史场次洪水过程推算流域净雨量,以实际净雨量与SCS模型预报净雨量差值最小为依据率定SCS模型产流参数CN和η,实际净雨量通过流域出口断面流量过程反算得到。

(2) 采用试算法率定分布式单位线参数Ic。由于目标流域场次洪水的降水强度多集中在5～15 mm/h,根据水量平衡原理重新分配后,蓄满区雨强范围多为20～60 mm/h,在参数Ic率定时,先根据流域历史平均降水强度假定几组离散值,本文假定20 mm、40 mm、60 mm,然后在其中优选,具体参见文献[11]；此外,推求分布式单位线还需获取流速系数、流域坡度等信息,其中,流速系数依据植被覆盖情况进行确定,不同植被类型对应的流速系数见文献[6],流域植被类型、地形坡度通过Arcgis工具提取。

(3) 根据时段净雨及式(12)进行汇流计算,该过程需要依据流域蓄满比例选择单位线,流域蓄满比例的计算关键在于确定流域蓄水容量空间分布,因此,采用1.4节所提方法率定蓄水容量参数n、WM,max及WM,min,目标函数权重确定参考文献[22],w1—w4分别为0.5、0.25、0.15、0.1。根据率定后的流域CN值计算流域潜在蓄水能力(S),以该值作为流域平均蓄水容量,进一步试算参数n和WM,max,研究发现,在流域平均蓄水容量不变的前提下,参数WM,min不敏感,变化范围较小,本文取值5 mm,参数n、WM,max对流域蓄水容量的空间分布有显著影响,参数n取值越大,蓄水容量空间分布越不均匀,WM,max取值也越大。本文参数率定结果如表1。

表1 龙虎圩流域和东石流域参数取值方案Table 1Parameter schemes of the Longhuwei River basin and Dongshi River basin

表2 龙虎圩流域和东石流域场次洪水率定结果及方案选择Table 2Calibration accuracy and schemes selection results in the Longhuwei River basin and Dongshi River basin

由表2可知,龙虎圩流域和东石流域率定期洪峰相对误差在±20%以内,峰现时间误差均在±6 h之间,2个流域的ENS均值分别为0.86和0.82,场次洪水合格率超80%。

2.3 动态汇流路径的时变分布式单位线推求

根据2.2节确定的单位线参数可以推求考虑动态汇流路径的时变分布式单位线和现有时变分布式单位线,为方便工程实际应用,结合相关研究[10-11,23],本文采用忽略上游入流贡献的式(8)进行未蓄满栅格土壤含水量计算；同时,为减少单位线条数,将降水强度和流域蓄满比例离散化,净雨量和流域蓄满比例因子离散区间见表3和表4。

表3 每个时段的雨强对应的离散值Table 3Discrete value Is of excess rainfall intensity It/40 corresponding to each time interval

表4 每个时段的土壤含水量状态对应离散土壤含水状态Table 4Discrete value αs of soil moisture content αt corresponding to each time interval

以龙虎圩流域为例,获得蓄水容量空间分布后,根据其蓄满比例获得具体蓄满区域,不同蓄满比例下(以αt取0.25、0.50、0.75为例),流域蓄满区域分布如图1所示,结合表3、表4和流域蓄满面积可以推求考虑动态汇流路径的时变分布式单位线。以推求蓄满比例为0.25时的单位线为例,即图1(a)中绿色区域所产生的单位净雨在流域出口形成的流量过程线(图2中绿色的单位线),相应地,可以推求不同蓄满比例下的汇流单位线,现有时变分布式单位线结果如图3所示。

图1 龙虎圩流域不同蓄满比例对应的产流面积空间分布Fig.1 Runoff areas corresponding to different full storage ratio in the Longhuwei River basin

由图2可知,在同一雨强下,考虑动态汇流路径的时变分布式单位线在流域蓄满比例较小时单位线峰值也较小,单位线上涨历时略有提前趋势,但不显著；在同一蓄满比例下,随着降水强度的增加,流域汇流速度增加,汇流时间减少,单位线峰现时间提前,总历时减小。本文所提考虑动态汇流路径的时变分布式单位线的物理意义为单位时段、单位净雨均匀降落在蓄满区域,在流域出口形成的流量过程线；图3为现有时变分布式单位线,其物理意义为单位时段、单位净雨均匀降落在全流域,在流域出口形成的流量过程线。两者在降水及下垫面时空分布、汇流速度及汇流路径等方面均存在显著差异。

图2 龙虎圩流域考虑动态汇流路径的时变分布式单位线Fig.2 Time-varying distributed unit hydrographs considering dynamic flow routing paths for the Longhuwei River basin

图3 龙虎圩流域现有时变分布式单位线Fig.3 Current time-varying distributed unit hydrographs for the Longhuwei River basin

2.4 与现有分布式单位线对比分析

根据SCS模型产流计算结果,分别采用本文所提单位线和现有分布式单位线进行模型检验,选取龙虎圩流域和东石流域各4场洪水进行预报,当流域前期较为干燥时选用参数方案1,否则采用参数方案2,预报结果如表5。

表5 本文方法和现有方法预报结果对比Table 5Comparisons of the results predicted by the proposed and current methods

由表5可知,采用考虑动态汇流路径的时变分布式单位线进行汇流计算后,龙虎圩流域和东石流域场次洪水平均纳什效率系数在0.85以上,洪峰相对误差均值约为8%,平均峰现时间为1.65 h,预报精度达到工程应用乙级及以上精度,且优于现有时变分布式单位线预报结果,其中龙虎圩流域预报效果优于东石流域。以龙虎圩流域20030517号、20120527号场次洪水为例展开分析,采用本文所提动态汇流单位线进行径流预报时,这2场洪水预报精度有较大提升,2场洪水预报径流过程如图4所示；龙虎圩流域20060601号、东石流域20161021号、20190609号场次洪水存在精度提升不明显或各指标没有同步提升现象,原因如下：龙虎圩流域20060601号洪水前期土壤含水量达106 mm,流域接近蓄满状态,动态流径的时变分布式单位线和现有方法推求的单位线差距较小；东石流域20161021号、20190609号洪水产流参数不具有普适性。

图4 场次洪水预报流量与实测流量对比Fig.4 Comparisons of the predicted and observed flow for flood events 20030517 and 20120528

由于采用相同的产流参数,采用2种单位线进行预报的产流量是一样的,因此,预报洪水过程误差主要由单位线的差异造成。由图4可知,采用现有时变分布式单位线进行汇流计算时,20030517号洪水和20120527号洪水预报误差均高于所提方法,其中,20030517号洪水洪峰相对误差超过14%,ENS为0.88；20120527号洪水峰现时间误差为3 h,ENS为0.87。采用所提单位线进行预报后,2场洪水过程的ENS分别提高至0.94和0.93。结果表明,考虑动态汇流路径的时变分布式单位线可以更好地表征流域的汇流过程。

2.5 流域汇流时间规律解析

分布式单位线的核心为栅格汇流时间的推求,为解析动态汇流路径对流域汇流时间的影响机制,下面深入分析流域汇流时间的动态变化规律。采用动态汇流路径后,有一部分蓄满栅格由于流经未蓄满区域,导致汇流时间发生变化,为比较流域不同蓄满比例下土壤含水量对栅格汇流速度的影响,离散净雨强度设定为最大值Is=2.0,通过提取流域最先达到蓄满状态的25%流域面积(F0%—F25%)的栅格汇流时间,分析该部分面积上栅格在不同流域蓄满比例下的汇流时间变化,并与采用现有时变流速计算公式的汇流时间进行比较,以龙虎圩流域为例,流域最先蓄满的25%面积上的栅格汇流时间分布如图5。

图5 龙虎圩流域25%蓄满栅格汇流时间分布Fig.5 Flow routing time distribution of the 25% full storage grid in the Longhuwei River basin

图5展示了图1(a)中绿色区域的栅格汇流时间,αt=1.00的橙色区域代表现有时变流速计算公式的结果,横坐标代表栅格,按照蓄水容量从小到大排序,以图5F18%处对应的栅格为例,其在蓄满状态为0.25、0.50、0.75、1.00时对应的汇流时间分别为9 h、7 h、6 h、5 h,可以看出随着流域蓄满比例的增加,蓄满栅格的汇流时间逐渐减少,这是因为随着降水的持续,处于蓄满路径的未蓄满栅格单元也逐渐趋于饱和,未蓄满栅格单元逐渐减少,由式(7)可知,栅格蓄满后流速达到最大值,因此,随着流域蓄满比例的增加,栅格汇流时间逐渐缩短,直至流域达到全蓄满状态,此时动态汇流路径转化为静态汇流路径,流域汇流时间趋于稳定。

图6为龙虎圩流域不同蓄满比例下考虑动态汇流路径时计算的汇流时间和静态汇流路径计算的汇流时间箱线图,结果表明,现有时变流速公式会高估坡面汇流速度,汇流时间不随蓄满比例而改变,考虑土壤含水量对流域汇流速度的影响后,随着流域逐渐达到全部蓄满状态,即αt趋于1.00时,2种方法汇流时间趋于一致,汇流时间分布更加集中；此外,随着流域蓄满面积的增加,土壤含水量对汇流时间的影响有减小趋势,所提汇流速度计算方程对前期较为干旱的情况影响较大。

图6 本文所提方法和现有方法不同蓄满比例下流域汇流时间箱线图Fig.6 Box plots of flow routing time under different storage ratio of the proposed and current methods

3 结 论

(1) 针对现有模型中的蓄水容量曲线只能定性描述流域蓄水容量空间不均匀分布的缺陷,引入地形指数刻画流域栅格蓄水容量的空间不均匀分布,定量表征每一个栅格蓄水容量,结合产流计算结果,进而从栅格尺度划分流域蓄满栅格和未蓄满栅格,可为栅格流速场计算提供更高精度输入场。

(2) 提出了考虑时变降水强度和土壤含水量的汇流速度计算公式,通过划分4种栅格状态(蓄满栅格、未蓄满但处于汇流路径栅格、未蓄满且不处于汇流路径栅格、河道栅格),更加精确地刻画了流速场的空间分布规律。

(3) 建立了考虑动态汇流路径的时变分布式单位线,针对传统单位线假定全流域均匀产流的不合理假定,本文以南方山区龙虎圩流域和东石流域为研究对象,基于蓄满区域更易产流的理论基础,推求了流域不同蓄满比例下的动态汇流路径,完善了单位线推求的理论基础。

(4) 结合本文所提改进时变流速计算公式和动态汇流路径理论,计算了不同蓄满比例下栅格汇流时间分布场,具体而言,通过改进流速公式提高了汇流时间计算的精度,通过划分蓄满区域和未蓄满区域,将汇流方法的汇流路径从静态汇流路径发展为动态汇流路径,应用结果表明,所提方法应用精度较高,提高了分布式单位线汇流时间计算的准确性。