一种基于曼宁公式的河道枯水流量计算方法

本发明涉及河道水动力，尤其涉及一种基于曼宁公式的河道枯水流量计算方法。

背景技术：

1、随着人类活动和全球气候变化的影响，枯水期补水调度成为了流域生态管理最为重要的工作之一。最大可能节约用水、科学有效调度有限的水资源以精准满足流域生态用水需求具有现实意义，而这项工作离不开河道枯水流量的计算。因此，高效、准确模拟出河道枯水流量的演进过程是流域实施水量调度和生态保护等工作的前提。

2、河道水流模拟计算目前主要有经验方法、水文学方法和水动力学方法三大类。其中，经验方法依据统计相关分析，构建以线性系统为主的经验性模型；水文学方法从水量平衡原理和河道槽蓄角度出发，建立河道水流演进的半经验性模型；水动力学方法则主要通过求解圣维南方程组实现非恒定水流数值计算。其中经验方法类似于黑箱模型，对资料要求低，计算精度也相对较低，一定程度上可解释为运动波的理论求解，方法主要运用于水文模型中的产汇流计算，在流域枯水、补水调度等水流计算中很少应用；水文学方法是水动力学方法的简化，方法没有对地形资料的要求，是运动波方程的差分形式，计算远没有动力波求解复杂，但模型结果在精度上不占优势，在河道洪水演算中应用较为广泛；而水动力学方法通过求解完整的圣维南方程组计算得到流量水位时空分布信息，能综合反映各水力要素之间的动态联系，演算结果更能满足枯水期水量调度的精度需求。尽管如此，枯水期水流小、水位低、甚至有干河床情形，采用水动力学方法常出现数值不稳定、迭代不收敛等问题，即便采用一些激波捕捉技术或者窄缝方法进行处理也无法解决频繁的干湿交替问题。

技术实现思路

1、本发明的目的在于提供一种基于曼宁公式的河道枯水流量计算方法，从而解决现有技术中存在的前述问题。

2、为了实现上述目的，本发明采用的技术方案如下：

3、一种基于曼宁公式的河道枯水流量计算方法，包括如下步骤，

4、s1、数据收集及格栅化处理：

5、对收集或获取的河道地形数据、河床下渗数据和模型参数进行格栅化处理，以将各河道划分为多个单元格，并为各单元格进行数据输入和参数赋值；

6、s2、单元格水流流向计算：

7、基于d8算法，通过计算当前单元格与其周围单元格之间的坡降，将单元格水流方向与相邻八个单元格之间坡降最大值代表的方向作为该单元格的水流流向；

8、s3、单元格水流流速计算：

9、基于单元格的水流流向及对应坡降，利用曼宁公式计算单元格的水流流速，并对流速大的急流情况进行修正；

10、s4、单元格水量交换量和单元格水位计算：

11、分别计算单元格水量交换量所包含的单元格流入流出量和单元格渗漏水量，基于获取的单元格水量交换量，利用水量平衡原理，计算各单元格水位；

12、s5、单元格水量修正：

13、所有单元格水位计算完毕后，利用各单元格流入扣除流出和下渗水量得到的单位面积不平衡水量对所有单元格逐个进行水量平衡修正；

14、s6、单元格水位信息更新：

15、更新所有单元格水位信息，进入下一水流计算时段。

16、优选的，步骤s1具体包括如下内容，

17、s11、收集河道地形数据和河床下渗数据，并根据河道河床情况确定模型参数；

18、s12、将河道地形数据、河床下渗数据和模型参数进行格栅化处理，按一定分辨率将河道划分为多个正方形的单元格，并对所划分的每一个单元格进行数据输入和参数赋值。

19、优选的，步骤s12中，单元格的最大行数由河道沿程最大宽度决定，单元格的最大列数由河道长度决定；

20、m＝bmax/ds    (1)

21、n＝len/ds    (2)

22、其中，m为河网划分单元格的最大行数；n为河网划分单元格的最大列数；bmax为沿着水流方向最大河道宽度；len为河道长度；ds为河道网格划分分辨率，即所划分的单元格的边长。

23、优选的，步骤s2具体包括如下内容，

24、s21、以中心单元格为计算对象，依次将其周围的单元格进行编码；

25、s22、根据上一时段计算的每个单元格水量结果，计算当前第i行、第j列的单元格(i,j)与周围单元格之间的坡降；

26、slope(i,j)＝dz/ds    (3)

27、

28、其中，z为单元格水位，dz为中心单元格和相邻八个单元格的水位差；单元格(i,j)与相邻单元格之间的距离与方向有关，对于位于单元格(i,j)正上、正左、正右、正下方向上的相邻单元格，单元格(i,j)与相邻单元格之间的距离为ds，坡降计算公式采用方程式(3)；对于位于单元格(i,j)左上、右上、左下、右下方向上的相邻单元格，单元格(i,j)与相邻单元格之间的距离为坡降计算公式采用方程式(4)；

29、s23、根据单元格(i,j)坡降计算结果，假设每个计算时间步长内水流方向唯一，单元格(i,j)水流方向为与该单元格与相邻八个单元格间坡度最大值代表的方向；

30、aspect(i,j)＝f(max slope(i,j))    (5)

31、其中，aspect(i,j)为计算单元格与周围单元格的水流交换方向代号；f()为最大坡度对应单元格指示的水流交换方向代号。

32、优选的，步骤s23中最大坡度对应单元格指示的水流交换方向代号f()采用如下规则确定，

33、若maxslope(i,j)≤0，则赋以aspect(i,j)为当前单元格代号值，说明当前单元格水流不流出；

34、若maxslope(i,j)>0,且只有一个相邻单元格对应坡降最大，则赋以aspect(i,j)为对应相邻单元格的编码值；

35、若maxslope(i,j)>0,且有两个及以上相邻单元格对应坡降最大值，则按照顺时针方向赋以aspect(i,j)编码号。

36、优选的，步骤s3具体包括如下内容，

37、s31、根据单元格水流流向及对应坡降，由曼宁公式计算单元格水流流速，

38、

39、其中，v为单元格水流流速；ns为曼宁糙率系数；r为水力半径，近似为单元水深；

40、s32、根据弗洛德数判定流态的方法，对速度限定上限阀值，以对流速大的急流情况进行修正；

41、

42、其中，g为重力加速度；h(i,j)为单元格(i,j)的水深。

43、优选的，步骤s4具体包括如下内容，

44、s41、基于单元格水流流速，计算单元格流出水量；

45、dw(i,j)＝v·δt·ds·h(i,j)    (8)

46、

47、其中，δt为计算时段长；dw(i,j)为单元格(i,j)流出水量，当流向为正上、正左、正右、正下时，dw(i,j)由公式(8)计算，当流向为左上、右上、左下、右下时，dw(i,j)由公式(9)计算；

48、s42、基于霍顿产流公式，计算单元格渗漏水量；

49、fp＝[fc+(f0-fc)e-βt]h(i,j)    (10)

50、其中，fp为单元格渗漏水量；f0为干河床时最大下渗率，随着入渗过程进行，减小到近似常数的入渗率fc；参数β控制着入渗能力的下降速度；f0、fc和β根据实测入渗资料率定；

51、s43、基于单元格流出水量和单元格渗漏水量，根据水量平衡原理，计算单元格水位；

52、

53、其中，z(t)和z(t+δt)分别为计算时段初和时段末的单元格水位；∑dw′(i,j)为相邻单元格流入单元格(i,j)的水量。

54、优选的，步骤s5具体为，计算完毕所有单元格水位后，利用如下公式，对所有单元格逐个进行水量平衡修正，

55、z‘(t+δt)＝z(t+δt)+δz (12)

56、其中，δz为每个单元格流入扣除流出和下渗水量得到的单位面积不平衡水量；z‘(t+δt)为平衡修正后的水位。

57、优选的，步骤s6具体为，完成当前时段所有单元格的水流流速、水位计算后，将相关计算结果作为下一时段的初始值，返回步骤s2。

58、本发明的有益效果是：1、本发明提供的河道枯水流量计算方法，解决了河道枯水流量小水深或干河床情况下的演算不稳定问题。一方面，方法保持了一维度水动力学方法的简单实用，在保证水流计算效率和计算精度的基础上，解决了以往水动力学方法在小水深、干河床情况下的不适用或不稳定问题，另一方面，方法能实现河道水流的空间二维演算，可更为准确地刻画出天然河道小水深情况下横向流速、水位的空间差异。2、本发明方法实现了枯水期河道水流的稳定演算，可为枯水河道补水调度和生态保护等工程实践提供技术支持。3、本发明方法提供了一种权衡计算效率和计算精度的水流演算新思路，可推广于相关类似工程实践中，尤其是枯水期补水调度、生态调度、水资源管理等问题。