一种时变分布式动态汇流计算方法

摘要

本发明公开了一种时变分布式动态汇流计算方法，属于径流预报领域。本发明引入地形指数刻画流域栅格蓄水容量的空间不均匀分布，可定量表征每一个栅格蓄水容量，结合产流计算结果，进而可以从栅格尺度划分流域蓄满栅格和为蓄满栅格，可为栅格流速场计算提供更高精度的输入场；本发明提出了考虑时变降雨强度和土壤含水量的汇流速度计算公式，更加精确的刻画了流速场的空间分布规律；本发明基于蓄满区域产流、不蓄满不产流的理论基础，推求了流域不同蓄满比例下的动态汇流路径，完善了单位线推求的理论基础，提高了分布式单位线汇流时间计算的准确性，解决了传统分布式单位线预报洪峰偏大的问题，提高了洪水预测的准确性。

权利要求

1.一种时变分布式动态汇流计算方法，其特征在于，包括：S1.基于SCS模型，计算流域产流量；

S2.将流域进行栅格划分，并基于地形指数计算流域内所有栅格单元的蓄水容量；

基于栅格单元蓄水容量、当前时刻降雨量和当前时刻净雨量，得到流域蓄满区域和未蓄满区域对应栅格；在栅格尺度上计算当前时刻蓄满区域面积和未蓄满区域土壤含水量；

并结合产流计算结果和流域面积，将当前时刻净雨量分配至蓄满区域；

S3.根据栅格单元土壤蓄水状态，构建时变坡面流速计算公式；其中，栅格单元土壤蓄水状态指栅格是否蓄满，以及未蓄满栅格的土壤含水量；

S4.基于不同时刻下流域蓄满区域面积，得到不同蓄满比例下流域的动态汇流路径；

S5.基于时变坡面流速计算公式和动态汇流路径，得到时变分布式动态单位线；

S6.采用所述时变分布式动态单位线和流域产流量计算结果，计算流域洪水过程。

2.根据权利要求1所述的一种时变分布式动态汇流计算方法，其特征在于，第j栅格单元的蓄水容量WMj计算公式为：式中，WMmin为全流域所有栅格中的最小蓄水容量；WMmax为全流域所有栅格中的最大蓄水容量；TImin＝min{TIj,j＝1,2…,N}；TImax＝max{TIj,j＝1,2…,N}；TIj为第j个栅格的地形指数；n为待率定参数，反映地形指数分布与蓄水容量分布的系数；N为流域栅格单元总数。

3.根据权利要求2所述的一种时变分布式动态汇流计算方法，其特征在于，t时段流域蓄满区域面积为：式中， 为流域蓄满区域的面积，αt为t时刻流域蓄满区域比例，Wj,0为j栅格初始时刻蓄水量；Pi为i时段降雨量，Ri为i时段净雨量，Aj为栅格单元面积。

4.根据权利要求3所述的一种时变分布式动态汇流计算方法，其特征在于，t时段净雨量Rt重新分配至蓄满区域面积，表达式如下：F为流域面积；Rt为t时段净雨量，即产流量，Rt'为蓄满区域净雨深，

5.根据权利要求4所述的一种时变分布式动态汇流计算方法，其特征在于，时变坡面流速计算公式为：k为流速系数，βj为j栅格单元坡度，Ic为流域平均净雨强度，Wj,t为蓄满区域对应的j栅格单元t时段土壤含水量。

6.根据权利要求2或5所述的一种时变分布式动态汇流计算方法，其特征在于，未蓄满区域对应的j栅格单元t时段土壤含水量为：Wj,t＝Wj,t‑1+Pt‑Rt

Wj,t‑1为j栅格t‑1时段蓄水量，Pt为t时段降雨量，Rt为t时段净雨量，即产流量。

7.根据权利要求5所述的一种时变分布式动态汇流计算方法，其特征在于，不同蓄满比例下的动态汇流路径为：表示蓄满栅格的汇流路径集合， Rj表示栅格单元至流域出口的汇流路径集合。

8.根据权利要求7所述的一种时变分布式动态汇流计算方法，其特征在于，步骤S5具体包括：

01.计算t时段流域蓄满比例为αt时栅格水滴在栅格中的等待时间 当水流沿栅格方向流动时 当水流沿栅格对角方向流动时 Lj为栅格边长；

沿汇流路径累加得到每个栅格至流域出口的汇流时间

02.计算相邻单位线时段对应时间间隔内流域蓄满面积：表示蓄满栅格汇流时间在[(m‑1)△t,m△t]区间内的面积总和，满足△t为单位线时间间隔，单位为小时；m为单位线时段，m＝1,2,…,M，M为单位线总时段数，

03.将蓄满栅格面积转换为流域出口流量过程，得到任意时段的分布式单位线：UH(m)为第m个时段的单位线流量，单位为m3/s；H为单位净雨，单位为mm。

说明书

一种时变分布式动态汇流计算方法

技术领域

[0001] 本发明属于径流预报领域，更具体地，涉及一种时变分布式动态汇流计算方法。

背景技术

[0002] 高精度水文预报可为水库调度、防洪减灾及水资源优化配置等提供决策依据，在水资源管理、水资源开发利用和国民经济建设等领域发挥着重要作用。降雨径流建模主要包括产流和汇流两部分，单位线是应用最为广泛的汇流计算方法，常用的单位线方法包括：时段单位线、瞬时单位线、综合单位线、分布式单位线等。

[0003] 传统汇流方法的应用具有一定局限性，主要问题在于将整个流域汇流过程概化为线性系统，无法准确刻画降雨时空分布特征和下垫面的空间异质性，其次，依赖详细的水文气象数据，无法获得精度较高的径流预报结果。在此背景下，考虑了流域汇流非线性效应的时变分布式单位线方法因所需资料较少、推求过程简单、应用精度较高而备受推崇，然而，现有分布式汇流理论存在两个缺陷：其一，仅考虑降雨时空分布不均匀性和下垫面空间分布异质性对流域汇流速度的影响，其假定流域在降雨脉冲结束之前就可以达到饱和状态，这与实际情况不符，忽略了时变土壤含水量造成下垫面状态的改变；其次，现有分布式单位线假定全流域产流，其采用静态的汇流路径推求单位线，理论上，以蓄满产流为主的南方湿润地区，在流域没有达到全流域蓄满之前，只有蓄满区域产流，因此，流域汇流路径应随着蓄满区域的改变而动态变化，当全流域蓄满后，此时动态汇流路径转换为静态汇流路径。

[0004] 因此，现有分布式单位线汇流方法应用时存在高估流域坡面汇流速度，假定流域存在唯一静态汇流路径，造成洪水预报精度偏低的技术问题。

发明内容

[0005] 针对现有技术的以上缺陷或改进需求，本发明提供了一种时变分布式动态汇流计算方法，其目的在于解决现有时变分布式单位线洪水预报时存在高估汇流速度、洪水预报误差较大的技术问题。

[0006] 为实现上述目的，本发明提供了一种时变分布式动态汇流计算方法，包括：

[0007] S1.基于SCS模型，计算流域产流量；

[0008] S2.将流域进行栅格划分，并基于地形指数计算流域内所有栅格单元的蓄水容量；

[0009] 基于栅格单元蓄水容量、当前时刻降雨量和当前时刻净雨量，得到流域蓄满区域和未蓄满区域对应栅格；在栅格尺度上计算当前时刻蓄满区域面积和未蓄满区域土壤含水量；并结合产流计算结果和流域面积，将当前时刻净雨量分配至蓄满区域；

[0010] S3.根据栅格单元土壤蓄水状态，构建时变坡面流速计算公式；其中，栅格单元土壤蓄水状态指栅格是否蓄满，以及未蓄满栅格的土壤含水量；

[0011] S4.基于不同时刻下流域蓄满区域面积，得到不同蓄满比例下流域的动态汇流路径；

[0012] S5.基于时变坡面流速计算公式和动态汇流路径，得到时变分布式动态单位线；

[0013] S6.采用所述时变分布式动态单位线和流域产流量计算结果，计算流域洪水过程。

[0014] 进一步地，第j栅格单元的蓄水容量WMj计算公式为：

[0015]

[0016] 式中，WMmin为全流域所有栅格中的最小蓄水容量；WMmax为全流域所有栅格中的最大蓄水容量；TImin＝min{TIj,j＝1,2…,N}；TImax＝max{TIj,j＝1,2…,N}；TIj为第j个栅格的地形指数；n为待率定参数，反映地形指数分布与蓄水容量分布的系数；N为流域栅格单元总数。

[0017] 进一步地，t时段流域蓄满区域面积为：

[0018]

[0019] 式中， 为流域蓄满区域的面积，αt为t时刻流域蓄满区域比例，Wj,0为j栅格初始时刻蓄水量；Pi为i时段降雨量，Ri为i时段净雨量，Aj为栅格单元面积。

[0020] 进一步地，t时段净雨量Rt重新分配至蓄满区域面积，表达式如下：

[0021]

[0022] F为流域面积；Rt为t时段净雨量，即产流量，Rt'为蓄满区域净雨深，[0023] 进一步地，时变坡面流速计算公式为：

[0024]

[0025] k为流速系数，βj为j栅格单元坡度，Ic为流域平均净雨强度，Wj,t为蓄满区域对应的j栅格单元t时段土壤含水量。

[0026] 进一步地，未蓄满区域对应的j栅格单元t时段土壤含水量为：

[0027] Wj,t＝Wj,t‑1+Pt‑Rt

[0028] Wj,t‑1为j栅格t‑1时段蓄水量，Pt为t时段降雨量，Rt为t时段净雨量，即产流量。

[0029] 进一步地，不同蓄满比例下的动态汇流路径为：

[0030]

[0031] 表示蓄满栅格的汇流路径集合， Rj表示栅格单元至流域出口的汇流路径集合。

[0032] 进一步地，步骤S5具体包括：

[0033] 01.计算t时段流域蓄满比例为αt时栅格水滴在栅格中的等待时间 当水流沿栅格方向流动时 当水流沿栅格对角方向流动时 Lj为栅格边长；

[0034] 沿汇流路径累加得到每个栅格至流域出口的汇流时间

[0035]

[0036] 02.计算相邻单位线时段对应时间间隔内流域蓄满面积：

[0037]

[0038] 表示蓄满栅格汇流时间在[(m‑1)△t,m△t]区间内的面积总和，满足△t为单位线时间间隔，单位为小时；m为单位线时段，m＝1,2,…,M，M为单位线总时段数，

[0039] 03.将蓄满栅格面积转换为流域出口流量过程，得到任意时段的分布式单位线：

[0040]

[0041] UH(m)为第m个时段的单位线流量，单位为m3/s；H为单位净雨，单位为mm。

[0042] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果。

[0043] 本发明针对传统模型中的蓄水容量曲线只能定性描述流域蓄水容量空间不均匀分布的缺陷，本发明引入地形指数刻画流域栅格蓄水容量的空间不均匀分布，可定量表征每一个栅格蓄水容量，结合产流计算结果，进而可以从栅格尺度划分流域蓄满栅格和为蓄满栅格，可为栅格流速场计算提供更高精度的输入场；

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

[0045] 本发明建立了考虑动态汇流路径的时变分布式单位线，针对传统单位线假定全流域均匀产流的不合理假定，本发明基于蓄满区域产流、不蓄满不产流的理论基础，推求了流域不同蓄满比例下的动态汇流路径，完善了单位线推求的理论基础；

[0046] 本发明结合所提改进时变流速计算公式和动态汇流路径理论，计算了不同蓄满比例下栅格汇流时间分布场，具体而言，通过改进流速公式提高了汇流时间计算的精度，通过划分蓄满区域和未蓄满区域，将汇流方法的汇流路径从静态汇流路径发展为动态汇流路径，总体而言，提高了分布式单位线汇流时间计算的准确性，更具物理机制；

[0047] 综上，本发明提出了一种动态路径汇流理论，推求了同时考虑降雨强度和下垫面土壤含水量时空分布特征的坡面流速计算公式，从栅格尺度上更加精细地刻画了流域流速分布场，在此基础上，推求了时变分布式动态单位线，解决了传统分布式单位线预报洪峰偏大的问题，从物理成因层面解析了流域汇流速度的影响机理，提高了洪水预测的准确性。

附图说明

[0048] 图1本发明的时变分布式动态汇流方法流程图；

[0049] 图2本发明实施例提供的研究区域图；

[0050] 图3研究区域植被分布图；

[0051] 图4研究区域地形坡度分布图；

[0052] 图5研究区域地形指数分布曲线；

[0053] 图6中(a)‑(c)不同蓄满比例下下垫面饱和面积分布示意图；

[0054] 图7所提基于动态汇流路径的时变分布式单位线；

[0055] 图8传统时变分布式单位线；

[0056] 图9所提方法与传统方法净雨对比图；

[0057] 图10场次洪水预报效果对比图。

具体实施方式

[0058] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

[0059] 参考图1，本发明提供的一种时变分布式动态汇流计算方法，包括以下步骤：

[0060] (1)基于SCS模型，进行流域产流计算；

[0061] 步骤(1)所述基于SCS模型，进行流域产流计算，具体表达式如下：

[0062]

[0063]

[0064] 式中，R为地表总径流量，mm；P为降雨总量，mm；S为流域当前可能最大滞留量，mm；CN值无量纲，反映降雨前期流域特征；Ia为初始损失量，mm，Ia＝mS，m为待率定参数。

[0065] 进一步地，t时段净雨量可通过下式推求：

[0066]

[0067] 式中，Pi为i时段降雨量；Rt为t时段净雨量(产流量)；

[0068] (2)推求考虑流域地形指数的分布式流域蓄水容量曲线，在此基础上，统计当前时刻流域蓄满面积，计算未蓄满栅格单元当前土壤含水量，进一步结合产流计算结果将流域平均净雨重新分配至蓄满区域；

[0069] 步骤(2)所述推求考虑流域地形指数的分布式流域蓄水容量曲线，首先，基于流域DEM数据提取地形指数，具体表达式如下：

[0070]

[0071] 式中，TIj为第j个栅格的地形指数；Asj为上坡区通过单位等高线汇集到j栅格单元的面积；βj为j栅格单元坡度。

[0072] 进一步，根据流域地形指数，推求反映流域土壤水空间异质性的分布式蓄水容量曲线，表达式为：

[0073]

[0074] 式中，WMj为j栅格单元的蓄水容量；WMmin为全流域所有栅格中的最小蓄水容量；WMmax为全流域所有栅格中的最大蓄水容量；TImin＝min{TIj,j＝1,2…,N}；TImax＝max{TIj,j＝1,2…,N}；n为待率定参数，反映地形指数分布与蓄水容量分布的系数；N为流域栅格单元总数。

[0075] t时段流域蓄满面积为：

[0076]

[0077] 式中，Wj,0为j栅格初始时刻蓄水量； 为流域蓄满部分的面积；Aj为栅格单元面积。

[0078] 进一步地，未蓄满区t时段土壤含水量为：

[0079] Wj,t＝Wj,t‑1+Pt‑Rt                   (7)

[0080] 式中，Wj,t为j栅格t时段蓄水量；Wj,t‑1为j栅格t‑1时段蓄水量。

[0081] 进一步地，由于仅蓄满区发生产流，而SCS模型认为全流域均匀产流，因此，t时段净雨量Rt需重新分配至蓄满面积，表达式如下：

[0082]

[0083] 式中，F为流域面积；Rt'为蓄满区域净雨深，其余参数意义同上。则实际蓄满区净雨深可表示为：

[0084]

[0085] (3)考虑流域下垫面状态时空变异特性，提出时变坡面流速计算公式；

[0086]

[0087] 式中， 为当流域蓄满比例为α时j栅格单元的汇流速度；k为流速系数，依据下垫面植被、土壤类型进行确定；Ic为流域平均净雨强度，其余参数意义同上。

[0088] (4)根据所提时变坡面流速计算公式和流域蓄满面积，推求考虑土壤含水量的分布式动态单位线；

[0089] 01.根据D8算法确定栅格水流方向，并根据栅格水流方向提取该栅格至流域出口的汇流路径集合Rj(第j个栅格至流域出口路径上的栅格单元集合)，则全流域所有栅格汇流路径集合可表示为：

[0090] R＝{Rj|j＝1,2,…,N}                   (11)

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

[0092]

[0093] 式中， 表示蓄满栅格的汇流路径集合， 表示集合 真包含于集合集合R表示流域静态汇流路径的集合；其余参数意义同上。由于流域不同蓄满比例下不同，因此，其构成不同蓄满比例下的动态汇流路径。

[0094] 02.根据所提流速计算公式计算每个栅格逐时刻的汇流速度 获得不同流域蓄满比例αt下的流速分布场；

[0095] 03.计算t时段流域蓄满比例为αt时栅格水滴在栅格中的等待时间 并沿汇流路径累加得到每个栅格至流域出口的汇流时间

[0096]

[0097]

[0098] Lj为栅格边长，当水流沿栅格方向流动时 当水流沿栅格对角方向流动时

[0099] 04.计算相邻单位线时段对应时间间隔内流域蓄满面积：

[0100]

[0101] 表示蓄满栅格汇流时间在[(m‑1)△t,m△t]区间内的面积总和，满足△t为单位线时间间隔，单位为小时；m为单位线时段，m＝1,2,…,M，M为单位线总时段数，

[0102] 05.将蓄满栅格面积转换为流域出口流量过程，得到任意时段的分布式单位线：

[0103]

[0104] UH(m)为第m个时段的单位线流量，单位为m3/s；H为单位净雨，单位为mm；

[0105] (5)采用所提时变分布式动态单位线，结合产流计算结果，计算流域洪水过程。

[0106] 以龙虎圩流域为研究对象，研究区域图如图2所示。由于本发明关注流域汇流方法计算精度，采用SCS模型进行流域产流量计算，选取历史多场典型洪水进行模型参数率定，模拟时间步长为1h，流域内CN值根据土壤类型、土地利用、降水前土壤湿润情况确定，取值参考美国SCS模型CN取值表，初损取值为0.25S，研究流域植被覆盖情况见图3。

[0107] 将流域以30m×30m为单元进行栅格划分，对DEM图层进行填洼、提取坡度、提取流向、栅格化河网处理，研究区域地表坡度分布如图4所示；

[0108] 提取流域地形指数，在此基础上，推求分布式蓄水容量曲线，龙虎圩流域地形指数分布如图5所示，地形指数分布区间为[2.66,21.98]，以新安江模型蓄水容量分布曲线为参考，率定得到流域栅格最小蓄水容量和最大蓄水容量分别5mm、103mm，参数n取值为0.95，则蓄水容量分布曲线表达式为：

[0109]

[0110] 式中，SMj为j栅格单元的蓄水容量。

[0111] 根据蓄水容量分布曲线，当流域蓄满面积比例为α时，若通过SCS模型计算得到下一时段净雨为P，则下一时段蓄满区域净雨深为 进一步可以通过地形指数分布情况计算得到流域蓄满区域分布情况，以流域蓄满比例0.25、0.5、0.75为例，对应的流域蓄满区域分布如图6中(a)‑(c)所示。

[0112] 根据下垫面植被分布情况确定流速系数，相应植被对应地坡面流速系数k取值如表1所示。

[0113] 表1 不同植被类型对应的流速系数

[0114]

[0115] 确定流域参考雨强，根据研究区域历史平均雨强计算参考雨强，目标流域参考雨强为10mm/h，当流域蓄满面积比例为α时，该时刻流域栅格流速计算公式如下：

[0116]

[0117] 由上式可知，当某一栅格属于河道单元时，流速固定为2m/s；当该栅格处于蓄满区域时，流速仅跟雨强相关；当该栅格处于未蓄满区域，但处于蓄满栅格的流动路径上时，考虑土壤含水量对栅格流速的影响；当该栅格处于未蓄满区域且无径流时，栅格流速为0，不参与汇流路径统计。

[0118] 进一步，为简化所提时变单位线应用，减少单位线条数，本发明实施例将降雨强度和流域蓄满比例离散化，净雨量和流域蓄满比例因子离散结果如表2、3所示：

[0119] 表2 每个时段的雨强 对应离散值

[0120]

[0121] 表3 每个时段的土壤含水量状态αt对应离散土壤含水状态αs

[0122]

[0123] 由于流域蓄满比例与土壤含水量可以互相转换，给定蓄满比例时，联合分布式蓄水容量曲线即可推求未蓄满区域土壤含水量。

[0124] 根据所提时变坡面流速公式和流域饱和面积，推求考虑土壤含水量的分布式动态单位线。结合步骤(2)推求的不同蓄满比例下的动态汇流路径，以及步骤(3)计算的时变栅格流速，统计获得流域时间‑面积柱状图，获得的时变分布式动态单位线如图7所示，传统时变分布式单位线如图8所示，其主要区别在于本发明所提方法采用了动态汇流路径，某一蓄满状态下流域汇流路径只与蓄满栅格有关，因此，同一雨强下不同土壤含水状态对应的单位线有显著差异，由此可见，所提时变分布式动态单位线在物理机制上更为清晰。

[0125] 采用所提时变分布式动态单位线进行场次洪水预报，并将预报结果与传统分布式单位线进行对比，以20130713场次洪水为例，参数CN为73，根据SCS模型计算得到该场降雨净雨量为92mm，该场次洪水前期土壤含水量为24mm，结合分布式蓄水容量曲线可知，流域初始蓄满比例达70％，采用所提方法和传统分布式单位线方法计算得到的逐时段净雨对比如图9所示，通过将全流域净雨分配到时变的蓄满区域，各个时段的净雨均出现不同程度的增加，随着流域接近饱和，所提方法和传统方法净雨深趋于相等；洪水预报结果如图10所示，结果表明采用所提方法在洪水全过程和起涨阶段均优于传统时变分布式单位线，总体上，所提预报方法精度优于传统时变分布式单位线。

[0126] 本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。