CN111143498B - 小河流洪水预报方法 - Google Patents

Abstract

本发明公开了一种小河流洪水预报方法，利用开源GIS代码进行二次开发，能满足洪水预报的GIS需求，同时通过模型计算引擎的分类，实现了模型的灵活配置，对流域的自动提取和节点拓扑关系的匹配关联，实现了方案的自动灵活构建，在提高中小河流预报精度的同时可很大程度节省方案作业人员的工作量。

Description

小河流洪水预报方法

技术领域

本发明涉及水利水电工程技术领域，特别是一种小河流洪水预报方法。

背景技术

我国中西部地区多为山区丘陵地带，山洪灾害频发，受灾地区域较为分散，具有流域面积较小、汇流时间较短、资料匮乏，且由于中小河流预报断面众多、上下游拓扑关系复杂，人工构建方案效率低、参数调整不便等特点。

传统中小河流预报方法的流程主要为：选择新安江三水源模型或其他水文模型(包括集总式、分布式水文模型)等，利用已有流域图层或使用ARCGIS等商业GIS工具导入DEM数据提取流域图层等相关要素；再通过经纬度录入雨量站信息，判断流域内相关代表站，通过ARCGIS等工具或者人工划分泰森多边形，求取相关权重；如果遇到上下游有断面关联，则人工建立上下游断面拓扑关系；最后选择多年洪水资料，率定模型参数，形成最终洪水预报方案。

上述预报方案建立流程，在大江大河流域面积较大、资料丰富、代表雨量站较多的情况下，可以基本满足预报要求。但在面对中小河流预报时，存在以下问题无法满足要求：1)中小河流预报流域面积较小无代表雨量站；2)中小河流预报断面众多，提取子流域、建立拓扑关系、选择模型等具有较大工作量；3)资料匮乏，汇流参数无法率定，预报精度较低；4)方案较多时，调整修改方案任务繁重。

发明内容

本发明所要解决的技术问题是，针对现有技术不足，提供一种小河流洪水预报方法，提高中小河流预报精度，节省方案作业工作人员的工作量。

为解决上述技术问题，本发明所采用的技术方案是：一种小河流洪水预报方法，包括以下步骤：

1)提取选定的试验流域的坡度、流向，并以所述坡度、流向为输入，从所述试验流域出口倒推，寻找最长汇流路径；

2)根据所述最长汇流路径和根基栅格点坡度，提取地貌单位线；

3)根据DEM数据，利用SharpMap提取相应流域，根据流域图层跟雨量站点图层交互情况，关联流域雨量站编码；为每个断面关联相应的上、下游断面雨量站编码；根据所述雨量站编码和最长汇流路径判断其他预报断面是否在汇流路径上，如果是则作为上游输入断面，并关联相关编号的上游站；如果不是，则确定没有上游输入断面；通过所述上游输入断面，将与其对应的下游断面构成上下游断面拓扑关系，确定整个流域关联的雨量站编码，构建计算下游流量的方案。

步骤1)中，利用下述方法提取试验流域的流向：利用SharpMap读取DEM每一个栅格点的数据，选择百分比法：两点的高程差与其水平距离的百分比，其计算公式如下：坡度＝(高程差/水平距离)x100％，使用百分比表示时，即：i＝h/l×100％，其中h为高程差，l为水平距离；每个栅格点均有8个流向，根据DEM每个栅格数据及邻近栅格数据，判断流向为流向最低高程的栅格点，最终形成每个栅格的D8流向。该方法简单容易实现。

所述最长汇流路径的计算过程包括：根据提取的流向，从流域出口断面处的栅格点，根据其D8方向的栅格点，倒推流向，根据每个栅格周围8个流向中流向本栅格的个数，采用倒推的方式，寻找所有汇流路径，根据累计长度最终找到最长汇流路径。

本发明中，为了便于实施，试验流域斜对角的流向距离记录为√2，上下左右记录为1。

步骤2)中，提取地貌单位线的具体实现过程包括：计算流速和坡度，根据所述流速和坡度计算径流在每一个网格内的持留时间；根据所述持留时间计算出每一个网格单元中径流的汇流时间；将所有流域内网格单元的汇流时间看做随机变量，统计分析该随机变量，得到汇流时间的概率密度分布，进而得到瞬时单位线；将所述瞬时单位线转换为无因次单位线；通过所述无因次单位线计算所述地貌单位线。利用开源控件提取，摆脱了商业软件束缚，且可以更快的提取需要的所有信息。

步骤3)中，若所述试验流域内有相应的雨量站，则选择相应的雨量站作为代表站；若没有站点，则选择临近流域雨量站作为代表站，存入数据库PSTATIONS列；对每个预报断面都进行编码，根据所提取的汇流路径，判断其他预报断面是否在汇流路径上，如果有则作为上游输入断面，关联对应编码的上游站。弥补流域缺少雨量站点无数据的缺点，且通过编码可以更快的匹配相应的雨量数据。

步骤3)之后，还包括：

4)构建分类模型库，对计算引擎进行分类存储。

所述分类模型库包括产汇流模型、产流模型、河道汇流模型、坡面汇流模型、经验模型；利用所述分类模型库输出预报流量和水位。

小河流通常是指二、三级支流以下以分水岭和下游河道出口断面为界集水面积在100km2以下的相对独立和封闭的自然汇水区域。

与现有技术相比，本发明所具有的有益效果为：本发明利用了开源GIS代码进行二次开发，能满足洪水预报的GIS需求，同时通过模型计算引擎的分类，实现了模型的灵活配置，对流域的自动提取和节点拓扑关系的匹配关联，实现了方案的自动灵活构建，在提高中小河流预报精度的同时可很大程度节省方案作业人员的工作量。

附图说明

图1为方法计算流程图。

图2为方案创建流程图。

图3为方案数据库存储表结构。

图4为模型管理数据流图。

图5为模型存储界面形式。

具体实施方式

以下结合附图说明本发明方法的具体实现过程。

1、选定试验流域。

金洞流域进行验证，该站位于湘江一级支流白水中游，为湘中南山区代表站，流域面积795km²，自1959年1月开始水文测验，测验项目有水位、降水、蒸发、流量。

经雨量资料分析，金洞水文站以上区间内以及周边地区有金洞、荒塘、白果市、庙背冲、磨刀岭、内下、三湾脚、西江、下阳、晒北滩、白沙源11个站雨量资料，其中，内下、三湾脚、西江、下阳、晒北滩、白沙源6个雨量站近年来数据不完整，不具备编制洪水预报方案的要求，其它5个站具有1小时雨量摘录资料，具备编制洪水预报方案的要求。

2、选择新安江产流模型耦合地貌单位线算法。

耦合地貌单位线作为汇流计算的算法：

设想流域上有一瞬时注入的、分布均匀的净雨量P，可看作由n各水质点所组成。若每一个水质点的体积为u，则有P＝nu。假设各水质点间呈弱相互作用，即在统计上可认为是相互独立的，而且各水质点以相等的可能性降落在流域上。

令

表示第i个水质点从其所在地到达流域出口所花费的时间，i＝1,2,…n。或者说，如果水质点i在0时刻注入流域，则它将在流域中必须持留

时间才能到达流域出口断面。因此，那些在流域中持留的时间超过时间t的水质点将不能在t时刻到达流域出口断面，即对t时刻的出口断面流量无贡献，但却是t时刻流域蓄水量w(t)的一部分[1,2]，故有：

式中

指示函数，当

时，其值取1，否则取0。

由于n为无穷大，且

有独立、等可能性，故根据大数定律，必有

式中E[.]为求数学期望；

P{T_b≥t}为事件{T_b≥t}发生的概率；

综合上式可得到：

w(t)＝I₀P{T_b≥t}＝I₀[1-P{T_b≤t}]＝I₀[1-F_B(t)]

另一方面，在流域汇流阶段,其水量平衡方程为

对于仅在0时刻有净雨量瞬时注入流域的情况，上式可简化为

根据上式可得：

式中f_b(t)一持留时间T_b的概率密度函数，f_B(t)＝dF_b(t)/dt。

根据流域瞬时单位线的定义，由上式可得：

上式揭示了组成净雨量的各水质点呈弱相互作用，并且有等可能性时，流域瞬时单位线与水质点持留时间的概率密度函数等价。

3、通过利用SharpMap开源控件提供的方法，通过下载30m DEM数据作为输入，开发可自动提取坡度、流向等地貌参数的程序。

利用SharpMap读取DEM(投影坐标)每一个栅格点的数据，选择百分比法：两点的高程差与其水平距离的百分比，其计算公式如下：坡度＝(高程差/水平距离)x100％，使用百分比表示时，即：i＝h/l×100％。每个栅格点均有8个流向，根据DEM每个栅格数据及邻近栅格数据，判断流向为流向最低高程的栅格点，最终形成每个栅格的D8流向。

4、以提取的坡度、流向为输入，从流域出口倒推，寻找最长汇流路径。

根据提取的流向，从流域出口断面处的栅格点，根据其D8方向的栅格点，倒推流向，根据每个栅格周围8个流向中流向本栅格的个数，记录后继续以本方法倒推，寻找所有汇流路径，斜对角的流向距离记录为√2，上下左右记录为1，根据累计长度最终找到最长汇流路径(如图1)。

5、根据计算的最长汇流路径，栅格点坡度，提取地貌单位线，作为汇流计算方法。

通过流速公式V＝aS^b计算流速。其中坡度可根据DEM数据计算，坡度取决于表面从中心像元开始在水平(dz/dx)方向和垂直(dz/dy)方向上的变化率(增量)。

径流在每一个网格内的持留时间可由下式计算：

或

式中L为网格的边长。

可由下式计算出每一个网格单元中径流的汇流时间。

式中m为径流路径上网格的数量。

按上述的方法，可以得到流域内所有网格单元的汇流时间，将汇流时间看作是随机变量，进行统计分析，可得到汇流时间的概率密度分布。(横坐标为汇流时间，纵坐标为与某汇流时间对应的面积Ac与流域面积A的比值)。

通过u(Δt,t_k)＝S(t_k)-S(t_k-Δt)将上述计算的瞬时单位线转换为无因次单位线。

通过

计算出10mm单位线。

6、根据DEM，利用SharpMap提取相应流域，根据流域图层(面图层)跟雨量站图层(点图层)交互情况，关联流域雨量站编码；根据节点为每个断面关联相应的上下游断面编码(见图2)。

每个雨量站都有编码，如果流域内有相应的雨量站，则选择相应的雨量站作为代表站。如果没有站点则选择临近流域雨量站作为代表站，存入数据库PSTATIONS列。对每个预报断面都进行编码，如金洞断面编码ZX000001，根据上一步提取的汇流路径，判断其他预报断面是否在汇流路径上，如果有则作为上游输入断面，如ZX069201,ZX069207。则关联ZX069201,ZX0692073为ZX000001的上游站。以逗号隔开的方式存储在数据库ARCD列中，其中-1表示上游没有输入断面，具体存储形式如图3所示。

7、构建分类模型库，对计算引擎进行分类保存，包括产汇流模型、产流模型、河道汇流模型、坡面汇流模型、经验模型等。通过不同的产汇流模型组合和扩展，实现模型的灵活配置(如图4)。

设定标准的输入输出格式，形成单个的模型库，每个模型库进行分类(如图5)。

8、分析计算结果，对比分析利用地貌单位线汇流模型耦合新安江产流模型的方法和传统水文预报之间的精度误差。

传统方法选择新安江模型+滞后演算法进行对比，可以看出地貌单位线在洪峰误差等方面预报精度更高。

Claims (8)

1.一种小河流洪水预报方法，其特征在于，包括以下步骤：

1)提取选定的试验流域的坡度、流向，并以所述坡度、流向为输入，从所述试验流域出口倒推，寻找最长汇流路径；

2)根据所述最长汇流路径和根基栅格点坡度，提取地貌单位线；

3)根据DEM数据，利用SharpMap提取相应流域，根据流域图层跟雨量站点图层交互情况，关联流域雨量站编码；为每个断面关联相应的上、下游断面雨量站编码；根据所述雨量站编码和最长汇流路径判断其他预报断面是否在汇流路径上，如果是则作为上游输入断面，并关联相关编号的上游站；如果不是，则确定没有上游输入断面；通过所述上游输入断面，将与其对应的下游断面构成上下游断面拓扑关系，确定整个流域关联的雨量站编码，构建计算下游流量的方案。

2.根据权利要求1所述的小河流洪水预报方法，其特征在于，步骤1)中，利用下述方法提取试验流域的流向：利用SharpMap读取DEM每一个栅格点的数据，选择百分比法：两点的高程差与其水平距离的百分比，其计算公式如下：坡度＝(高程差/水平距离)x100％，使用百分比表示时，即：i＝h/l×100％，其中h为高程差，l为水平距离；每个栅格点均有8个流向，根据DEM每个栅格数据及邻近栅格数据，判断流向为流向最低高程的栅格点，最终形成每个栅格的D8流向。

3.根据权利要求2所述的小河流洪水预报方法，其特征在于，所述最长汇流路径的计算过程包括：根据提取的流向，从流域出口断面处的栅格点，根据其D8方向的栅格点，倒推流向，即根据每个栅格周围8个流向中流向本栅格的个数，采用倒推的方式，寻找所有汇流路径，根据累计长度最终找到最长汇流路径。

4.根据权利要求3所述的小河流洪水预报方法，其特征在于，所述试验流域斜对角的流向距离记录为

上下左右记录为1。

5.根据权利要求1所述的小河流洪水预报方法，其特征在于，步骤2)中，提取地貌单位线的具体实现过程包括：计算流速和坡度，根据所述流速和坡度计算径流在每一个网格内的持留时间；根据所述持留时间计算出每一个网格单元中径流的汇流时间；将所有流域内网格单元的汇流时间看做随机变量，统计分析该随机变量，得到汇流时间的概率密度分布，进而得到瞬时单位线；将所述瞬时单位线转换为无因次单位线；通过所述无因次单位线计算所述地貌单位线。

6.根据权利要求5所述的小河流洪水预报方法，其特征在于，步骤3)中，若所述试验流域内有相应的雨量站，则选择相应的雨量站作为代表站；若没有站点，则选择临近流域雨量站作为代表站，存入数据库PSTATIONS列；对每个预报断面都进行编码，根据所提取的汇流路径，判断其他预报断面是否在汇流路径上，如果有则作为上游输入断面，关联对应编码的上游站。

7.根据权利要求1所述的小河流洪水预报方法，其特征在于，步骤3)之后，还包括：

4)构建分类模型库，对计算引擎进行分类存储。

8.根据权利要求7所述的小河流洪水预报方法，其特征在于，所述分类模型库包括产汇流模型、产流模型、河道汇流模型、坡面汇流模型、经验模型；利用所述分类模型库输出预报流量和水位。

Images