一种适用于半干旱半湿润地区的洪水预报方法与流程

本发明涉及半干旱半湿润地区的洪水预报研究领域，具体涉及一种基于地貌瞬时单位线改进双超模型的用于半干旱半湿润地区的洪水预报方法。

背景技术：

对于半干旱半湿润无资料地区的洪水分析，往往基于经验性公式推算设计洪水，而对洪水预报的研究相对匮乏。主要原因在于针对半干旱半湿润地区产流机制的模型较少，且水文资料的缺失致使流域汇流过程难以合理描述，因而半干旱半湿润无资料地区的洪水预报精度常常无法满足防洪需要。

针对半干旱半湿润地区的产流机制，王印杰提出了双超产流模式。他认为半干旱半湿润地区“超渗”“超持”两种产流现象同时存在，若降雨强度超过入渗能力，则由于“超渗”形成地表径流，下渗的水分补充土壤含水量，土壤水分达到饱和后由于“超持”形成壤中流和地下径流。基于双超产流模式建立的双超模型，包括双超产流模块及汇流模块的计算，模型先基于双超产流模块进行地表径流计算、壤中流和地下径流及蒸散发计算，再通过汇流模块对产流过程进行汇流计算，最终获取流域出口断面径流过程。

已有双超模型的汇流模块选用传统单位线法或滞后演算法对产流过程进行汇流计算。这两类汇流计算方法，或受到资料限制，或算法参数需依赖于降雨观测资料率定得到，使其在半干旱半湿润地区，尤其是无资料地区的应用较为困难。而随着地理信息系统(GIS)、计算机的发展与数字高程模型(DEM)的出现，以地貌特征反映流域汇流特性的地貌瞬时单位线(RGIUH)在描述流域汇流过程时展现出较强优势。RGIUH是一种基于物理基础的汇流模型，较黑箱方法推求的单位线更加能够反映流域实际情况，为无资料地区汇流机制研究提供一种新方法。

技术实现要素：

本发明旨在提供一种用于半干旱半湿润地区的基于改进双超模型的洪水预报方法，基于针对半干旱半湿润地区产流特点构建的双超模型，以具有物理基础的地貌瞬时单位线(RGIUH)改进双超模型中的汇流模块，建立适用于半干旱半湿润地区的洪水预报模型。RGIUH通过基于DEM提取的地貌特征反映流域汇流特性，为半干旱半湿润地区尤其是无资料地区的洪水预报研究提供新的方法。首先基于DEM获取流域地貌特征参数，建立研究区地貌瞬时单位线RGIUH，并以RGIUH改进双超模型中的汇流模块，构建双超—RGIUH模型，用于进行半干旱半湿润无资料地区的洪水预报研究，为提高洪水预报精度、完善半干旱半湿润地区洪灾防御体系具有重要的理论价值和现实意义。

本发明提供了一种用于半干旱半湿润地区的基于改进双超模型的洪水预报方法，包括：基于GIS平台对研究区DEM数据进行前期处理，根据斯特拉勒Strahler河流分级法获取子流域地貌特征参数；并基于地貌特征参数建立子流域R-V地貌瞬时单位线RGIUH；然后以地貌瞬时单位线法替代双超模型中的汇流计算方法构建双超-RGIUH模型，仍以双超产流模块进行产流计算，而以RGIUH汇流模块替代原汇流模块进行流域汇流计算，实现流域洪水过程模拟；洪水预报模型在参数率定后，用于进行半干旱半湿润地区的洪水预报。

上述方法中，基于GIS平台对研究区DEM数据进行前期处理的具体步骤如下：

(1)为保证水流方向判断正确、保证水流流出栅格，首先基于ArcGIS软件对原始DEM数据进行填洼处理，生成无洼地DEM；采用D8单流向算法计算水流方向，识别洼地区域，通过统计洼地区域的最低高程(指高程最小值)进而计算洼地深度，并分析洼地形成原因是否为数据误差；对于数据误差造成的洼地，需设置合理的阈值，对小于该阈值的洼地进行填充处理；该步骤需要进行多次；

(2)根据水流方向，将流域内每一个流向栅格汇聚的水量相加，获取整个流域的汇流累积量；

(3)生成河网：根据流域实际情况设定一个适当的栅格单元阈值，当栅格的汇流累积量大于该阈值则生成栅格河网，小于该阈值则生成汇流栅格，将所有的栅格河网汇聚连接就生成了整个流域的河网水系。

上述方法中，所述根据斯特拉勒Strahler河流分级法获取子流域的地貌特征参数的具体步骤如下：

(1)基于生成的河网水系，根据Strahler河流分级法定义河流级别，从河源出发的河流为1级，同级别河流交汇所形成的河流增加1级，不同级别河流交汇以其中级别较高者作为河流级别；

(2)获取各级别河流属性数据，并计算河数率、河长率、面积率：

河数率：RB＝Ni-1Ni,i＝2,3,…,ε；

河长率：i＝2,3,…,ε；

面积率：i＝2,3,…,ε；

式中，Ni为i级河流的数目，为i级河流的平均长度，为i级河流的流域平均面积，ε为水系最高级河流的级别。

上述方法中，所述基于地貌瞬时单位线理论建立子流域的地貌瞬时单位线RGIUH的具体步骤如下：

(1)确定初始概率：Ai为排入i级河流的流域面积；Aε为全流域面积；

(2)确定河流状态转移概率：

式中，δi+1,j为指示函数，当i+1＝j时，δi+1,j取1，否则取0；Ni、Nj分别为第i、j级河流的数目；j＝2,3,…,ε，k＝2,3,…,ε；Na、Na-1分别为第a、a-1级河流的数目；

(3)路径概率：其中为水质点处于初始状态的概率，为水质点从状态xk-1转移到状态xk的概率，即转移概率；

(4)将水质点到达出口断面时所花费的时间记为为水质点在各状态的持留时间；所有路径集合下总TS的概率密度函数即为地貌瞬时单位线RGIUH的表达式：

其中，s＝<x1,x2,…,xk>，S＝{s1,s2,…}为所有路径的集合，星号*为卷积运算符号，为的概率密度函数，采用(i＝1,2,…,k)为拟合线型，其中τi为水质点在状态xi的平均持留时间的倒数，即水质点在i级河流上的汇流时间的倒数，其中V为河流平均流速，为i级河流平均长度。

上述方法中，所述以RGIUH改进双超模型中的汇流模块，构建双超-RGIUH模型，仍以双超产流模块进行产流计算，而以RGIUH汇流模块替代原汇流模块进行流域汇流计算，实现流域洪水过程模拟的步骤包括：先通过双超产流模块计算地表径流、壤中流和地下径流以及蒸散发，再通过RGIUH汇流模块以RGIUH法进行地表径流的汇流计算、线性水库法进行壤中流及地下径流的汇流计算，得到流域出口断面径流过程，具体步骤如下：

(1)双超产流模块——地表径流计算：步骤如下：

1)根据雨前土壤含水率θ0及土壤孔隙率n计算雨前充水度B0：B0＝θ0/n；

2)设定虚构微元的土壤物理参数初始值：风干土壤的吸收率Sr(mm/h^1/2)、土壤的饱和导水率Ks、土壤孔径级配参数c，取值一般为3～6，入渗能力归一化分配曲线参数：b、α0；

3)计算虚构微元第1时段可能入渗量△Fm，假定降雨开始时刻t＝t0＝0，计算式为：

式中，为流域充水度，Δt为计算时段长。

4)根据时段有效降水量△P'计算有效供水度x＝△P/△Fm。

5)计算流域时段入渗量和地表径流深△F0和△RS：

△F0＝η(x)·△Fm；△RS＝△P'-△F0

式中η(x)为反映入渗量大小的供渗函数，计算式如下：

α0≤x<1，

其中，b、α0为入渗能力归一化分配曲线α0≤α≤1的参数，α0为反映流域产流临界雨强的参数，b为曲线指数。

6)计算虚构微元累积入渗量。计算虚拟微元自0到第i时段的累积入渗量：

将Fm,i代入式Fm(t)计算式，求解下一时段计算△Fm的起始时刻t0,i：

式中，

7)令t＝t0，计算下一时段入渗量△Fm；

8)重复步骤4)～7)至降雨结束，得到地表径流过程；

(2)双超产流模块——壤中流和地下径流计算：双超产流模块根据包气带内浅层土壤和深层土壤水分含量差异性较大的特点，把包气带从上到下细分为剧变、渐变、相对稳定、稳定四层，并用每层均含有底、侧开孔的串联填土容器水箱来模拟土壤中壤中流和地下径流的产生过程；底、侧孔时段总排水量△Si的计算公式如下：

①当H0,i+ΔFi-1-ΔEi≤Hc,i时，△Si＝0；

②当H0,i+ΔFi-1-ΔEi>Hc,i且H0,i≤Hc,i时，

③当H0,i+ΔFi-1-ΔEi>Hc,i且H0,i>Hc,i时，

则侧、底孔时段排水量分别为：

△RI,i＝δi△Si i＝1,2,3,4，△Fi＝(1-δi)△Si i＝1,2,3,4

因此，壤中流之和为：△RI＝∑δi△Si i＝1,2,3,4

地下径流为：△RG＝△F4

式中，△Si为时段总排水量，i为土层号，H0,i为时段初土壤含水量，Hc,i为i土层的田间持水量，mm；△Ei为i土层的时段蒸发量，mm；αi为i土层的排水系数，它是反映层自由水排水量的大小指数，△Fi-1为i-1土层向i层的下渗水量(mm)，δi为侧向排泄份额系数；

(3)双超产流模块——蒸散发计算：双超模型蒸散发计算模块只考虑最上面两层的土壤蒸发；土壤蒸发能力ΔE'm计算式为：

式中：△Em表示时段蒸发潜能，而非实际土壤蒸发量；实际土壤蒸发量是△Ei、ΔE1、ΔE2，实际土壤蒸发量还跟实际土壤含水率等有关；△P'为有效降雨强度；

按照先上层后下层的蒸散发顺序，每层先进行重力自由水的蒸发计算，然后再进行张力水的蒸发计算，上下层的土壤蒸散发ΔE1、ΔE2的计算公式是：

式中，Gi＝θi/θc,i；Gr＝θr/θc；θi、θc,i、θr分别为第i层容器土壤含水率、饱和含水率和凋萎含水率；ΔE1,ΔE2分别为第1，2层土壤蒸发量；Gi≤1表示第i层容器中土壤水分为张力水，否则容器内为自由水；

(4)RGIUH汇流模块进行汇流计算：以RGIUH汇流模块替代原汇流模块进行流域汇流计算，即地表径流汇流采用RGIUH法，根据壤中流及地下径流退水较慢的特点，二者汇流计算采用线性水库法，汇流计算公式为：

QS(t)＝RS(t)*RGIUH

QI(t)＝CI×QI(t-1)+(1-CI)×RI(t)×U

QG(t)＝CG×QG(t-1)+(1-CG)×RG(t)×U

Q(t)＝QS(t)+QI(t)+QG(t)

式中，QS(t)、QI(t)、QG(t)分别为地表径流、壤中流、地下径流的总入流，RS(t)、RI(t)、RG(t)分别为产流模块中输出的地表径流、壤中流、地下径流产流量，U为单位折算系数，U＝F(km²)/3.6Δt(h)，CI及CG分别为壤中流及地下径流消退系数，RGIUH为推求的地貌瞬时单位线RGIUH，Q(t)为流域总出流；

(5)双超-RGIUH模型参数率定与验证，包括产流参数土壤最大吸渗率Sr，土壤饱和导水率Ks，土壤孔径级配参数c，流域归一化曲线指数b，临界雨强因子α0，土层排水系数αi，侧向排泄份额系数δi，汇流参数地下径流消退系数CG、壤中流消退系数CI、河流的平均流速V；基于实测资料率定出一套合适的参数取值，完成洪水预报模型的本地化构建，并以此进行流域洪水模拟预报。

本发明的有益效果：

在由于资料缺乏难以获取传统单位线进行汇流计算的半干旱半湿润地区，可以通过基于DEM提取地貌参数进而构建地貌瞬时单位线RGIUH来进行汇流计算，并以其改进符合半干旱半湿润地区产流特点的双超模型，构建双超—RGIUH模型，对于提高无资料地区洪水预报精度、完善半干旱半湿润地区洪灾防御体系具有重要的理论价值和现实意义。

附图说明

图1为本发明的流程图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述，但不局限于以下实施例。

实施例1：

图1显示了本发明的一种基于改进双超模型的适用于半干旱半湿润地区的洪水预报方法的流程图。包括以下步骤：

步骤1，收集与处理研究区基础数据，包括降水径流数据及空间信息数据。空间信息数据包括DEM数据、土地利用分布图。

步骤2，DEM数据前期处理，包括生成无洼地DEM、获取汇流累积量及河网生成。此过程包括以下步骤：

(1)基于ArcGIS水文分析模块，采用D8单流向算法计算水流方向，识别洼地区域，通过统计洼地区域的最低高程进而计算洼地深度，并分析洼地形成原因是否为数据误差；设置合理的洼地深度阈值，对小于该阈值的洼地进行填充处理；多次进行本过程最终生成无洼地DEM。

(2)根据水流方向，将流域内每一个流向栅格汇聚的水量相加，获取整个流域的汇流累积量；根据流域实际情况设定一个适当的栅格单元阈值，当栅格的汇流累积量大于该阈值则生成栅格河网，小于该阈值则生成汇流栅格，将所有的栅格河网汇聚连接就生成了整个流域的河网水系。

步骤3，根据斯特拉勒(Strahler)河流分级法获取子流域的地貌特征参数。此过程包括以下步骤：

(1)将从河源出发的河流为1级，同级别河流交汇所形成的河流增加1级，不同级别河流交汇以其中级别较高者作为河流级别；

(2)获取各级别河流属性数据，并计算河数率、河长率、面积率等地貌特征值。

河数率：RB＝Ni-1/Ni,i＝2,3,…,ε

河长率：i＝2,3,…,ε

面积率：i＝2,3,…,ε

式中，Ni为i级河流的数目，为i级河流的平均长度，为i级河流的流域平均面积，ε为水系最高级河流的级别。各参数求均值后取常数参与后续计算。

步骤4，基于地貌特征参数建立子流域的地貌瞬时单位线RGIUH。此过程包括以下步骤：

(1)确定初始概率：Ai为排入i级河流的流域面积；Aε为全流域面积；

(2)确定河流状态转移概率：

式中，δi+1,j为指示函数，当i+1＝j时，δi+1,j取1，否则取0；Ni、Nj分别为第i、j级河流的数目；j＝2,3,…,ε，k＝2,3,…,ε；Na、Na-1分别为第a、a-1级河流的数目；

(3)路径概率：水质点通过某一路径s到达流域出口的概率表示为

其中为水质点处于初始状态的概率，为水质点从状态xk-1转移到状态xk的概率，即转移概率；

(4)识别水质点到达流域出口处的所有可能路径S＝{s1,s2,…}。以四级河流为例，共有8种可能路径，各路径的概率为：

式中参数根据(1)(2)中公式计算。

(5)计算概率密度函数即地貌瞬时单位线RGIUH的表达式：

其中，s＝<x1,x2,…,xk>，S＝{s1,s2,…}为所有路径的集合，星号*为卷积运算符号，为水质点持留时间的概率密度函数，采用(i＝1,2,…,k)为拟合线型，其中τi为水质点在状态xi的平均持留时间的倒数，即水质点在i级河流上的汇流时间的倒数，其中V为河流平均流速，为i级河流平均长度。

(6)将步骤3中所提取的地貌参数代入以上步骤中可最终获取地貌瞬时单位线RGIUH的表达式。

步骤5，以RGUIH改进双超模型中的汇流模块，构建双超-RGIUH模型，仍以双超产流模块进行产流计算，而以RGIUH汇流模块替代原汇流模块进行流域汇流计算，实现流域洪水过程模拟。通过以下步骤实现：

(1)双超产流模块——地表径流计算：逐时段计算地表径流的步骤如下：

1)根据雨前土壤含水率θ0及土壤孔隙率n计算雨前充水度B0：B0＝θ0/n；

2)设定虚构微元的土壤物理参数初始值：风干土壤的吸收率Sr(mm/h^1/2)、土壤的饱和导水率Ks、土壤孔径级配参数c，取值一般为3～6，入渗能力归一化分配曲线参数：b、α0；

3)计算虚构微元第1时段可能入渗量△Fm，假定降雨开始时刻t＝t0＝0，计算式为：

式中，为流域充水度，Δt为计算时段长。

4)根据时段有效降水量△P'计算有效供水度x＝△P/△Fm。

5)计算流域时段入渗量和地表径流深△F0和△RS：

△F0＝η(x)·△Fm；△RS＝△P'-△F0

式中η(x)为反映入渗量大小的供渗函数，计算式如下：

α0≤x<1，

其中，b、α0为入渗能力归一化分配曲线α0≤α≤1的参数，α0为反映流域产流临界雨强的参数，b为曲线指数。

6)计算虚构微元累积入渗量。计算虚拟微元自0到第i时段的累积入渗量：

将Fm,i代入式Fm(t)计算式，求解下一时段计算△Fm的起始时刻t0,i：

式中，

7)令t＝t0，计算下一时段入渗量△Fm；

8)重复步骤4)～7)至降雨结束，得到地表径流过程；

(2)双超产流模块——壤中流和地下径流计算：双超产流模块根据包气带内浅层土壤和深层土壤水分含量差异性较大的特点，把包气带从上到下细分为剧变、渐变、相对稳定、稳定四层，并用每层均含有底、侧开孔的串联填土容器水箱来模拟土壤中壤中流和地下径流的产生过程；底、侧孔时段总排水量△Si的计算公式如下：

①当H0,i+ΔFi-1-ΔEi≤Hc,i时，△Si＝0；

②当H0,i+ΔFi-1-ΔEi>Hc,i且H0,i≤Hc,i时，

③当H0,i+ΔFi-1-ΔEi>Hc,i且H0,i>Hc,i时，

则侧、底孔时段排水量分别为：

△RI,i＝δi△Si i＝1,2,3,4；△Fi＝(1-δi)△Si i＝1,2,3,4

因此，壤中流之和为：△RI＝∑δi△Si i＝1,2,3,4

地下径流为：△RG＝△F4

式中，△Si为时段总排水量，i为土层号，H0,i为时段初土壤含水量，Hc,i为i土层的田间持水量，mm；△Ei为i土层的时段蒸发量，mm；αi为i土层的排水系数，它是反映层自由水排水量的大小指数，△Fi-1为i-1土层向i层的下渗水量(mm)，δi为侧向排泄份额系数；

(3)双超产流模块——蒸散发计算：双超模型蒸散发计算模块只考虑最上面两层的土壤蒸发；土壤蒸发能力ΔE'm计算式为：

式中：△P'为有效降雨强度；

按照先上层后下层的蒸散发顺序，每层先进行重力自由水的蒸发计算，然后再进行张力水的蒸发计算，上下层的土壤蒸散发ΔE1、ΔE2的计算公式是：

式中，Gi＝θi/θc,i；Gr＝θr/θc；θi、θc,i、θr分别为第i层容器土壤含水率、饱和含水率和凋萎含水率；ΔE1,ΔE2分别为第1，2层土壤蒸发量；Gi≤1表示第i层容器中土壤水分为张力水，否则容器内为自由水；

(4)RGIUH汇流模块进行汇流计算：以RGIUH汇流模块替代原汇流模块进行流域汇流计算，即地表径流汇流采用RGIUH法，根据壤中流及地下径流退水较慢的特点，二者汇流计算采用线性水库法，汇流计算公式为：

QS(t)＝RS(t)*RGIUH

QI(t)＝CI×QI(t-1)+(1-CI)×RI(t)×U

QG(t)＝CG×QG(t-1)+(1-CG)×RG(t)×U

Q(t)＝QS(t)+QI(t)+QG(t)

式中，QS(t)、QI(t)、QG(t)分别为地表径流、壤中流、地下径流的总入流，RS(t)、RI(t)、RG(t)分别为产流模块中输出的地表径流、壤中流、地下径流产流量，U为单位折算系数，U＝F(km²)/3.6Δt(h)，CI及CG分别为壤中流及地下径流消退系数，RGIUH为推求的地貌瞬时单位线RGIUH，Q(t)为流域总出流；

(5)双超-RGIUH模型参数率定与验证，包括产流参数土壤最大吸渗率Sr，土壤饱和导水率Ks，土壤孔径级配参数c，流域归一化曲线指数b，临界雨强因子α0，土层排水系数αi，侧向排泄份额系数δi，汇流参数地下径流消退系数CG、壤中流消退系数CI、河流的平均流速V；基于实测资料率定出一套合适的参数取值，完成洪水预报模型的本地化构建，并以此进行流域洪水模拟预报。