自动化实现小流域山洪洪水流量和淹没预报的方法

1.本发明涉及洪水预报技术领域，具体涉及一种自动化实现小流域山洪洪水流量和淹没预报的方法。


背景技术：

2.洪涝灾害是比较严重的自然灾害之一，小流域山洪作为一种常见的洪涝灾害，具有发生频率高、破坏性突发性强和预报预测难度大等特点，对预报预警的准确性和模型的计算效率有较高的要求。
3.小流域山洪预报预警方式有：临界雨量和水文水动力模型。临界雨量采用经验性公式，以概化参数反映降雨的不均匀性，预报预警的精度和准确性易受到降雨的空间分布以及地形地貌的影响。集总式水文模型对于指定出口断面的流量预报效果较好，缺乏对河道其它断面的流量预报。分布式(栅格)水文模型能对多个河道断面进行流量预报，但编写复杂，且无法获得洪水淹没范围。


技术实现要素：

4.本发明为了克服以上技术的不足，提供了一种自动化实现小流域山洪洪水流量和淹没预报的方法，实现从降雨到流量预报到淹没范围的自动化实时预报。
5.术语解释：
6.1、dem：digital elevation model，数字高程模型。
7.2、scs-cn模型：the soil conservation servicerunoffcurvenumber method。
8.3、k值：流速系数。
9.4、cn值：curve number值。
10.本发明克服其技术问题所采用的技术方案是：
11.一种自动化实现小流域山洪洪水流量和淹没预报的方法，包括步骤：
12.步骤s1、获取预设流域的参数，包括栅格水流流向、栅格坡度、流域河网、河网宽度；
13.步骤s2、获取预设流域内每个栅格的cn值和k值；
14.步骤s3、计算预设流域内的每个栅格产流量和汇流情况；
15.步骤s4、构建基于开源模型lisflood-fp的山洪淹没模型；
16.步骤s5、通过山洪淹没模型自动计算并输出流域流量预报值和洪水淹没范围。
17.进一步地，步骤s1中，通过dem获取栅格水流流向、栅格坡度、流域河网、河网宽度，具体包括：
18.步骤s1.1、通过公共网站下载预设流域的dem数据；
19.步骤s1.2、采用arcgis对下载的dem进行填洼处理；
20.步骤s1.3、采用arcgis计算栅格水流流向；
21.步骤s1.4、采用arcgis计算相邻栅格的高程差，求得每个栅格的坡度；
22.步骤s1.5、采用arcgis，通过统计栅格水流流向，计算流入各个栅格上游的栅格数，获得栅格的集水面积，然后设置集水面积阈值，提取预设流域的河网；
23.步骤s1.6、采用arcgis，指定预设流域的一点，根据栅格水流流向寻找流域的出水口，综合分析河网、流域出水口和栅格水流流向，获得预设流域的河网；
24.步骤s1.7、获取河网宽度。
25.进一步地，步骤s2具体包括：
26.步骤s2.1、通过公共网站获取预设流域的土壤利用情况和土壤类型，将预设流域的土壤按照下渗率分为四类，根据每一类土壤的土壤利用情况查找cn值关系表，调查前期土壤含水情况，根据不同湿润情况和cn值关系表获得栅格的cn值；
27.步骤s2.2、根据土壤利用情况，查找坡度流速系数表，读取栅格的k值。
28.进一步地，步骤s3具体包括：
29.步骤s3.1、获取预设流域的雨量站的数据，通过插值计算得到面降雨数据；
30.步骤s3.2、通过scs-cn模型计算每个栅格的产流量；
31.步骤s3.3、基于步骤s1获取的栅格水流流向，寻找流域的出水口，标记为0，流入标记为0的栅格为1，流入标记为1的栅格的栅格为2，依此类推，直至遍历整个流域，读取栅格的编号，获得栅格汇流的先后顺序；
32.步骤s3.4、计算水流在栅格的流速，流速乘以时间步长，得到水流流过的距离，水流距离用以判断在时间步长内水流到达的栅格；基于步骤s3.2获得的每个栅格的产流量，根据栅格汇流的先后顺序，对流域进行汇流计算，栅格的汇流量等于当前栅格的产流量加上从其它栅格流入该栅格的水量；
33.步骤s3.5、采用马斯京根法，计算河道栅格流量。
34.进一步地，步骤s3.2中，计算每个栅格的产流量具体如下：
[0035][0036]
上式中，runoff为径流量，单位是mm；a为初损系数；p为时段降雨量；s为流域的最大可能滞留量，计算公式如下：
[0037][0038]
上式中，cn值由土壤类型、土地利用情况和土壤前期含水情况决定。
[0039]
进一步地，步骤s3.4中，计算水流在栅格的流速具体为：
[0040]
v＝ksb(3)
[0041]
上式中，v为流速；k为流速系数，至少由土地利用特征决定；b为经验系数；s为坡度。
[0042]
进一步地，步骤s4中，山洪淹没模型采用开源模型lisflood-fp进行计算，其中每个河道栅格点的河道水深采用浅水方程计算：
[0043]
[0044][0045]
上式中，q
x
和u
x
为x方向的流量和速度；a为河道断面面积；g为重力加速度；n为曼宁糙率系数；r为水力半径；h为水深；z为河床高程。
[0046]
进一步地，步骤s4具体包括：
[0047]
步骤s4.1、根据栅格水流流向，统计相邻河道栅格中流入该河道栅格的栅格数，流入起点河道栅格的河道栅格数为0，流入支流汇入栅格点的河道栅格数大于等于2，流域的出口点不汇入流域的其它河道栅格；
[0048]
步骤s4.2、从步骤s3计算得到的河道栅格流量过程中，提取河道的起点、出口以及支流汇入点的流量过程；
[0049]
步骤s4.3、根据河道特征和洪泛区特征，通过查找表格，获得河道曼宁糙率系数和洪泛区曼宁糙率系数；
[0050]
步骤s4.4、将预设流域的dem以及步骤s4.2得到的河道栅格流量过程、河道宽度、河道曼宁糙率系数、洪泛区的曼宁糙率系数、输入输出文件的地址、山洪洪水预报总时长、初始时间步长以及河道水流模拟是否采用扩散波，按照lisflood-fp输入要求准备数据。
[0051]
进一步地，步骤s5具体包括：
[0052]
步骤s5.1、运行lisflood-fp模型，生成exe文件；
[0053]
步骤s5.2、将步骤s4.4准备的数据作为输入文件，调用步骤s5.1生成的exe文件，生成淹没数据；
[0054]
步骤s5.3、读取淹没数据，将淹没数据以geotiff格式存储。
[0055]
进一步地，步骤s5.2中，lisflood-fp模型输出的淹没数据的格式默认为wd文件，通过预设算法将其改为输出csv、txt或excel格式文件。
[0056]
本发明的有益效果是：
[0057]
1、本发明实现了从降雨到流量预报到淹没范围的自动化实时预报。
[0058]
2、本发明的山洪淹没模型，以土壤类型、土地利用情况及前期水文条件等信息作为模型输入，使得本方法具有非常强的物理基础。
[0059]
3、本发明的山洪淹没模型采用栅格数据，能充分考虑空间分布的不均匀性。
[0060]
4、本发明中，降雨过程为面降雨过程，能较全面的考虑降雨的空间分布情况。
附图说明
[0061]
图1为本发明实施例所述的自动化实现小流域山洪洪水流量和淹没预报的方法的流程图。
[0062]
图2为本发明实施例所述的预设流域的edm的示意图。
[0063]
图3为本发明实施例所述的栅格水流流向的示意图。
[0064]
图4为本发明实施例所述的预设流域的cn值的示意图。
[0065]
图5为本发明实施例所述的预设流域的k值的示意图。
[0066]
图6为本发明实施例所述的预设流域的坡度的示意图。
[0067]
图7为本发明实施例所述的预设流域的河网的示意图。
[0068]
图8为本发明实施例所述的预设流域河网的宽度的示意图。
[0069]
图9为本发明实施例方法计算的某一时刻河网栅格的流量的示意图。
[0070]
图10为本发明实施例方法计算的某一时刻研究区域的淹没水深的示意图。
具体实施方式
[0071]
为了便于本领域人员更好的理解本发明，下面结合附图和具体实施例对本发明做进一步详细说明，下述仅是示例性的不限定本发明的保护范围。
[0072]
如图1所示，本实施例所述的一种自动化实现小流域山洪洪水流量和淹没预报的方法，包括步骤：
[0073]
步骤s1、获取预设流域的参数，包括栅格水流流向、栅格坡度、流域河网、河网宽度；
[0074]
步骤s2、获取预设流域内每个栅格的cn值和k值；
[0075]
步骤s3、计算预设流域内的每个栅格产流量和汇流情况；
[0076]
步骤s4、构建基于开源模型lisflood-fp的山洪淹没模型；
[0077]
步骤s5、通过山洪淹没模型自动计算并输出流域流量预报值和洪水淹没范围。
[0078]
基于栅格水流流向和流域河网，计算栅格汇流的拓扑关系；计算每个栅格水流流速，根据栅格拓扑关系，计算每个栅格的水流流量；对于河道栅格，采用马斯京根法，计算河道流量。
[0079]
对河网进行分析，获得河网的起始点、支流交汇点以及流域出口，提取这些栅格点的流量过程。对这些点的河宽进行测量，插值获得整个河网栅格的河宽。将栅格点的流量过程、河道栅格河宽、河网及流域其它栅格点的糙率系数和dem输入山洪淹没模型。
[0080]
具体地，下面以某区县为例，对本实施例所述的一种自动化实现小流域山洪洪水流量和淹没预报的方法进行示例性说明，但并不局限于本实施例。本实施例中，输入到山洪淹没模型的数据如图2-图8所示，其他输入数据包括预设流域(本实施例所述的预设流域是指本实施例所研究的流域)的面降雨以及预设流域任意一点。流域内给定的点用以寻找研究流域。将这些数据输入山洪淹没模型，自动完成计算，输出结果如图9和图10所示，即流域的洪水淹没范围。这里需要说明的是，图2-图10中，横坐标和纵坐标的单位均为栅格数。
[0081]
步骤s1、获取预设流域的参数，包括栅格水流流向、栅格坡度、流域河网、河网宽度。具体如下：
[0082]
步骤s1.1、通过公共网站下载预设流域的dem数据，可以从下面3个网站获取免费的dem数据：aster全球数字高程模型、航天飞机雷达地形任务(srtm)、全球数字地表模型(alos world 3d-30m)。
[0083]
步骤s1.2、为了能提取完整的河网，采用arcgis对下载的dem进行填洼处理，经本实施例填洼处理后的预设流域的dem如图2所示，图2右侧的图例数值表示高程，单位是米。
[0084]
步骤s1.3、采用arcgis计算栅格水流流向，具体是调用arcgis工具箱中的d8算法计算栅格水流流向，如图3所示，图3右侧的图例数值表示水流流向，数值和水流流向的对应关系为1：东、2：东南、4南、8：西南、16：西、32：西北、64：北、128：东北。
[0085]
步骤s1.4、采用arcgis计算相邻栅格的高程差，求得每个栅格的坡度，如图6所示，图6右侧的图例数值表示坡度角的正切值。
[0086]
步骤s1.5、采用arcgis，通过统计栅格水流流向，计算流入各个栅格上游的栅格
数，获得栅格的集水面积，然后设置集水面积阈值，提取预设流域的河网，如图7所示，白色的线条为河网，图7中包括多个子流域的河网。
[0087]
步骤s1.6、采用arcgis，指定预设流域的一点，根据栅格水流流向寻找流域的出水口，综合分析河网、流域出水口和栅格水流流向，获得预设流域的河网。
[0088]
步骤s1.7、获取河网宽度，可以通过实测获取河网宽度或通过卫星图像获得河网宽度，所得河网宽度数据如图8所示，图8右侧的图例数值表示河宽，单位是米。
[0089]
上述步骤s1.1-步骤s1.6除了采用上述所述的方法外，也可以自己编写代码完成。
[0090]
步骤s2、获取预设流域内每个栅格的cn值和k值。具体包括：
[0091]
步骤s2.1、通过公共网站(比如，清华大学10m/30m数据http://data.ess.tsinghua.edu.cn/等)获取预设流域的土壤利用情况和土壤类型，将预设流域的土壤按照下渗率分为a、b、c、d四类，a类土壤的下渗率大于7.5毫米/小时，b类土壤的下渗率3.8～7.5毫米/小时，c类土壤的下渗率1.3～3.8毫米/小时，d类土壤的下渗率小于1.3毫米/小时，根据每一类土壤的土壤利用情况查找cn值关系表，调查前期土壤含水情况，根据不同湿润情况和cn值关系表获得栅格的cn值，如图4所示，图4右侧的图例数值为cn值，通过图得到的数据量比较大且数据更直观。
[0092]
步骤s2.2、根据土壤利用情况，查找坡度流速系数表，读取栅格的k值，如图5所示，图5右侧的图例数值为k值，即流速系数。
[0093]
步骤s3、计算预设流域内的每个栅格产流量和汇流情况。具体包括：
[0094]
步骤s3.1、获取预设流域的雨量站的数据，该数据为点数据，通过插值计算得到面降雨数据，插值方法可以选取克里金插值、双线性插值、最近邻插值或反距离插值，根据插值效果，选取合适的插值方法。
[0095]
步骤s3.2、通过scs-cn模型计算每个栅格的产流量，计算每个栅格的产流量具体如下：
[0096][0097]
上式中，runoff为径流量，单位是mm；a为初损系数；p为时段降雨量；s为流域的最大可能滞留量，计算公式如下：
[0098][0099]
上式中，cn值由土壤类型、土地利用情况和土壤前期含水情况决定。
[0100]
步骤s3.3、基于步骤s1获取的栅格水流流向，寻找流域的出水口，标记为0，流入标记为0的栅格为1，流入标记为1的栅格的栅格为2，依此类推，直至遍历整个流域，读取栅格的编号，获得栅格汇流的先后顺序。
[0101]
步骤s3.4、计算水流在栅格的流速，流速乘以时间步长，得到水流流过的距离，水流距离用以判断在时间步长内水流到达的栅格；基于步骤s3.2获得的每个栅格的产流量，根据栅格汇流的先后顺序，对流域进行汇流计算，栅格的汇流量等于当前栅格的产流量加上从其它栅格流入该栅格的水量。
[0102]
具体地，计算水流在栅格的流速具体为：
[0103]
v＝ksb(3)
[0104]
上式中，v为流速；k为流速系数，至少由土地利用特征决定；b为经验系数；s为坡度。
[0105]
步骤s3.5、采用马斯京根法，计算河道栅格流量，如图9所示，图9右侧的图例数值为流量，单位是立方米每秒，从图9可以直观的看到流量的具体数值。
[0106]
步骤s4、构建基于开源模型lisflood-fp的山洪淹没模型，即，山洪淹没模型采用开源模型lisflood-fp进行计算，其中每个河道栅格点的河道水深采用浅水方程计算：
[0107][0108][0109]
上式中，q
x
和u
x
为x方向的流量和速度；a为河道断面面积；g为重力加速度；n为曼宁糙率系数；r为水力半径；h为水深；z为河床高程。
[0110]
步骤s4具体包括如下步骤：
[0111]
步骤s4.1、根据栅格水流流向，统计相邻河道栅格中流入该河道栅格的栅格数，流入起点河道栅格的河道栅格数为0，流入支流汇入栅格点的河道栅格数大于等于2，流域的出口点不汇入流域的其它河道栅格，通过该方法寻找河道的起点；
[0112]
步骤s4.2、从步骤s3计算得到的河道栅格流量过程中，提取河道的起点、出口以及支流汇入点的流量过程；
[0113]
步骤s4.3、根据河道特征和洪泛区特征，通过查找表格(比如一些公共网站或专业技术手册中都有曼宁糙率系数的表格)，获得河道曼宁糙率系数和洪泛区曼宁糙率系数；
[0114]
步骤s4.4、将预设流域的dem以及步骤s4.2得到的河道栅格流量过程、河道宽度、河道曼宁糙率系数、洪泛区的曼宁糙率系数、输入输出文件的地址、山洪洪水预报总时长、初始时间步长以及河道水流模拟是否采用扩散波，按照lisflood-fp输入要求准备数据。
[0115]
步骤s5、通过山洪淹没模型自动计算并输出流域流量预报值和洪水淹没范围。具体包括：
[0116]
步骤s5.1、运行lisflood-fp模型，生成exe文件。
[0117]
步骤s5.2、将步骤s4.4准备的数据作为输入文件，调用步骤s5.1生成的exe文件，生成淹没数据，lisflood-fp模型输出的淹没数据的格式默认为wd文件，本实施例可以通过预设算法将其改为输出csv、txt或excel格式文件。
[0118]
步骤s5.3、读取淹没数据，将淹没数据以geotiff格式存储，如图10所示，图10右侧的图例数值为淹没水深，单位为米，从图10可以直观的看到淹没水深的具体数值。
[0119]
采用本实施例所述的方法可以实现从降雨到流量预报到淹没范围的自动化实时预报。
[0120]
以上仅描述了本发明的基本原理和优选实施方式，本领域人员可以根据上述描述做出许多变化和改进，这些变化和改进应该属于本发明的保护范围。

技术特征：
1.一种自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，包括步骤：步骤s1、获取预设流域的参数，包括栅格水流流向、栅格坡度、流域河网、河网宽度；步骤s2、获取预设流域内每个栅格的cn值和k值；步骤s3、计算预设流域内的每个栅格产流量和汇流情况；步骤s4、构建基于开源模型lisflood-fp的山洪淹没模型；步骤s5、通过山洪淹没模型自动计算并输出流域流量预报值和洪水淹没范围。2.根据权利要求1所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s1中，通过dem获取栅格水流流向、栅格坡度、流域河网、河网宽度，具体包括：步骤s1.1、通过公共网站下载预设流域的dem数据；步骤s1.2、采用arcgis对下载的dem进行填洼处理；步骤s1.3、采用arcgis计算栅格水流流向；步骤s1.4、采用arcgis计算相邻栅格的高程差，求得每个栅格的坡度；步骤s1.5、采用arcgis，通过统计栅格水流流向，计算流入各个栅格上游的栅格数，获得栅格的集水面积，然后设置集水面积阈值，提取预设流域的河网；步骤s1.6、采用arcgis，指定预设流域的一点，根据栅格水流流向寻找流域的出水口，综合分析河网、流域出水口和栅格水流流向，获得预设流域的河网；步骤s1.7、获取河网宽度。3.根据权利要求1所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s2具体包括：步骤s2.1、通过公共网站获取预设流域的土壤利用情况和土壤类型，将预设流域的土壤按照下渗率分为四类，根据每一类土壤的土壤利用情况查找cn值关系表，调查前期土壤含水情况，根据不同湿润情况和cn值关系表获得栅格的cn值；步骤s2.2、根据土壤利用情况，查找坡度流速系数表，读取栅格的k值。4.根据权利要求1所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s3具体包括：步骤s3.1、获取预设流域的雨量站的数据，通过插值计算得到面降雨数据；步骤s3.2、通过scs-cn模型计算每个栅格的产流量；步骤s3.3、基于步骤s1获取的栅格水流流向，寻找流域的出水口，标记为0，流入标记为0的栅格为1，流入标记为1的栅格的栅格为2，依此类推，直至遍历整个流域，读取栅格的编号，获得栅格汇流的先后顺序；步骤s3.4、计算水流在栅格的流速，流速乘以时间步长，得到水流流过的距离，水流距离用以判断在时间步长内水流到达的栅格；基于步骤s3.2获得的每个栅格的产流量，根据栅格汇流的先后顺序，对流域进行汇流计算，栅格的汇流量等于当前栅格的产流量加上从其它栅格流入该栅格的水量；步骤s3.5、采用马斯京根法，计算河道栅格流量。5.根据权利要求4所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s3.2中，计算每个栅格的产流量具体如下：
上式中，runoff为径流量，单位是mm；a为初损系数；p为时段降雨量；s为流域的最大可能滞留量，计算公式如下：上式中，cn值由土壤类型、土地利用情况和土壤前期含水情况决定。6.根据权利要求4所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s3.4中，计算水流在栅格的流速具体为：v＝ks
b
(3)上式中，ν为流速；k为流速系数，至少由土地利用特征决定；b为经验系数；s为坡度。7.根据权利要求1所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s4中，山洪淹没模型采用开源模型lisflood-fp进行计算，其中每个河道栅格点的河道水深采用浅水方程计算：的河道水深采用浅水方程计算：上式中，q
x
和u
x
为x方向的流量和速度；a为河道断面面积；g为重力加速度；n为曼宁糙率系数；r为水力半径；h为水深；z为河床高程。8.根据权利要求1所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s4具体包括：步骤s4.1、根据栅格水流流向，统计相邻河道栅格中流入该河道栅格的栅格数，流入起点河道栅格的河道栅格数为0，流入支流汇入栅格点的河道栅格数大于等于2，流域的出口点不汇入流域的其它河道栅格；步骤s4.2、从步骤s3计算得到的河道栅格流量过程中，提取河道的起点、出口以及支流汇入点的流量过程；步骤s4.3、根据河道特征和洪泛区特征，通过查找表格，获得河道曼宁糙率系数和洪泛区曼宁糙率系数；步骤s4.4、将预设流域的dem以及步骤s4.2得到的河道栅格流量过程、河道宽度、河道曼宁糙率系数、洪泛区的曼宁糙率系数、输入输出文件的地址、山洪洪水预报总时长、初始时间步长以及河道水流模拟是否采用扩散波，按照lisflood-fp输入要求准备数据。9.根据权利要求8所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s5具体包括：步骤s5.1、运行lisflood-fp模型，生成exe文件；步骤s5.2、将步骤s4.4准备的数据作为输入文件，调用步骤s5.1生成的exe文件，生成淹没数据；步骤s5.3、读取淹没数据，将淹没数据以geotiff格式存储。
10.根据权利要求9所述的自动化实现小流域山洪洪水流量和淹没预报的方法，其特征在于，步骤s5.2中，lisflood-fp模型输出的淹没数据的格式默认为wd文件，通过预设算法将其改为输出csv、txt或excel格式文件。

技术总结
本发明公开了一种自动化实现小流域山洪洪水流量和淹没预报的方法，包括步骤：步骤S1、获取预设流域的参数，包括栅格水流流向、栅格坡度、流域河网、河网宽度；步骤S2、获取预设流域内每个栅格的CN值和K值；步骤S3、计算预设流域内的每个栅格产流量和汇流情况；步骤S4、构建基于开源模型LISFLOOD-FP的山洪淹没模型；步骤S5、通过山洪淹没模型自动计算并输出流域流量预报值和洪水淹没范围。本发明可以实现从降雨到流量预报到淹没范围的自动化实时预报。降雨到流量预报到淹没范围的自动化实时预报。降雨到流量预报到淹没范围的自动化实时预报。


技术研发人员：王乃玉 汪英俊 鄢尚
受保护的技术使用者：浙江大学
技术研发日：2023.01.06
技术公布日：2023/8/9