一种含行洪区河网结构的洪水模拟计算方法与流程

：
1.本发明涉及对洪水的形成和发展过程进行模拟和预测，特别涉及一种含行洪区河网结构的洪水模拟计算方法。


背景技术：

2.洪水模拟是以一系列控制方程为基础，通过赋予初始条件和边界条件，模拟河网中水流运动的一项数值技术。由于洪水灾害在全球范围内经常发生，且危害巨大，因此对洪水的形成和发展过程进行模拟和预测，在防洪减灾领域占有重要的地位。
3.行洪区是一类比较特殊的河网单元，兼具消峰、排洪两种功能。“消峰”指行洪区运用前期可以储存、迟滞一部分洪水；“排洪”指行洪区运用的中后期可以扩大河道断面，使所在河段的泄流能力大幅提升。淮河洪水的特点是峰高量大、持续时间长。为弥补河道蓄泄能力的不足，将干流两侧的洼地开辟为行洪区成为保障防洪安全的重要措施。据统计，沿淮行洪区在1950年至2021年之间，累计启用151处次，为降低干流洪峰水位，缩短高水位持续时间，保证淮北大堤等重要堤防的防洪安全发挥了重要作用。
4.随着流域经济的发展和河道整治工程的持续推进，流域内的行洪区数量一直在增减变化。截止2021年底，淮河王家坝至洪泽湖段共有行洪区8处，总面积达738km2，滞洪总库容为39.76亿m3。
5.由于行洪区蓄泄作用突出，运用频繁。研究含有行洪区河网结构的洪水模拟计算方法，准确快速地模拟洪水在主河道和行洪区之间的滞泄和吐纳过程，为决策者提供合理可行的调度方案，是近年来防洪减灾领域研究的热点问题。
6.现有技术对含行洪区河网结构水动力的模拟，主要有一维、平面二维(下文简称为二维)和三维。河网地区模拟范围较广，三维模型计算量大难以实用，故常用一维、二维以及一、二维耦合模型进行洪水推演。下面分别介绍它们的技术背景：
7.1)一维河网水动力数学模型
8.河流具有明显的一维特征，垂向尺度(水深)和横向尺度(河宽)比纵向尺度(河长)要小很多，因此国内外学者广泛采用一维模型计算河道水流的演进和运动规律。该模型的理论基础为saint-venant方程组，该方程组由连续方程和动量方程组成，见式1-2。
9.连续方程：
[0010][0011]
动量方程:
[0012][0013]
式中t为时间坐标，x为河道沿程坐标，q为流量，z为水位，a为过水断面的面积，b为水面宽度，k为流量模数，g为重力加速度，q为旁侧入流流量。
[0014]
上述方程的离散主要有显格式和隐格式两种。由于隐式差分格式稳定性好，能允
许较大的时间步长，所以在一维水动力数值模拟中应用最为广泛。离散方程组的系数矩阵具有明显的带状特点，求解时可采用直接解法、分级解法及松弛迭代法。
[0015]
直接法即直接求解由断面方程和边界方程构建的方程组。这是早期河网计算中常用的方法，需要很大的计算机存贮量和计算时长。分级解法基本思想是先将未知数集中在汊点，待汊点水位求出后，再分河段求解各断面的水力参数。松弛迭代法是将河网分解为一系列单一河道，各河道耦合汊点处的流量先给预估值，再用松弛迭代法校正，并逐步逼近精确值。目前，分级解法和松弛迭代法应用较多，而直接解法已逐步淡出主流应用市场。
[0016]
2)二维水动力模型
[0017]
对于河网中拥有宽广水域的湖泊、行洪区、河口而言，水流的速度和方向都是时变的，一维模型很难刻画细节，此时可采用二维模型。将navi－stocks方程沿水深进行积分，利用libnitz公式，即可得到沿水深平均的平面二维流动的基本方程式：
[0018]
连续性方程:
[0019][0020]
动量方程:
[0021][0022][0023]
式中:h为水深，ζ为自由面相对高度,h＝h ζ。u，v为笛卡尔坐标下x，y方向的深度平均水速分量,γ
φ
是有效粘性系数:γ
φ
＝η
t
η。
[0024]
fhu和fhv为考虑地球自转引起的科里奥利力(coriolis)力的作用。f为科里奥利力系数,f＝2ωsinψ,ω为地面自转角速度,ψ为当地纬度。
[0025]
τ
bx
和τ
by
为考虑底部阻力影响产生的附加项。
[0026][0027]cf
为底部摩擦系数:cf＝n2g/h
1/3
，n为底部糙率。
[0028]
与为考虑表面风力影响产生的附加项。
[0029]
ρa为大气密度,cw为风应力系数,w为水面上10m高处的风速值,β为风向与x坐标轴的夹角。
[0030]
上述方程的求解一般常用有限差分法、有限元法和有限体积法。有限差分法是比较经典的方法，其基本思想是将各微分方程中的微分项离散成微小矩形网格上各邻近节点的差商，得到一个以各节点上的函数值为未知变量的代数方程，然后进行求解；有限元法最先应用于固体力学，20世纪60年代开始在流体力学中有所应用。它吸收了有限差分法中离
散处理思想，同时采用了变分计算中选择逼近函数及对任意形状(三角形或四边形)的许多微小单元进行积分处理的合理方法，因而具有很广泛的适应性，特别适合于几何边界比较复杂的问题；有限体积法把计算区域离散为若干点，以这些点为中心，把整个计算区域分为若干互相连接但不重叠的控制体，这一点与有限元法是类似的。不同之处在于它使用守恒型的积分控制方程，在整个计算区域，不存在守恒量的误差。
[0031]
2)一、二维耦合水动力模型
[0032]
应用一维水动力模型可以较准确地模拟出河道内各断面的水位和流量变化情况，但在处理宽浅水域的对象时精度不高；应用二维水动力模型可以较准确地模拟出湖泊或行洪区流速的空间变化，但计算复杂耗时。将两种模型相结合，一维为二维提供来流量信息，二维为一维提供模拟区域的水位信息，通过互提边界就形成了一、二维耦合模型。
[0033]
一、二维耦合模型实际上是对大规模水域进行整体模拟的重要方法，它充分发挥了两种模型的优势并避免各自的不足。
[0034]
但是，针对一维模型、二维模型和一、二维耦合模型，这三类方法依旧存在缺点，具体缺点如下：
[0035]
(1)全一维模拟法主要有两种技术思路：一种是直接采用实测断面计算；另一种则是采用行洪区的概化断面进行计算。
[0036]
直接采用实测断面计算行洪区的蓄泄能力比较自然，大部分天然河道的洪水计算都是这么处理的，但行洪区有其特殊性——水面宽度远大于河道，不同区域的流速分布存在明显的差异，正是这种差异给实测断面计算带来了一些困难。根据行洪区二维模型的计算结果，水流进入行洪区后一般会在上下口门(或进退洪闸)及其连线形成高流速带，特别是口门附近水流收缩，流速增大明显。而在其他区域，行洪区的流速很低，甚至是零流速。所以貌似整体过流的行洪区，实际只有一部分起输水作用。如果仅仅是断面流速不均的问题，还可以通过引入动量校正系数进行修正，但行洪区内部还会出现回流，这就使得直接采用实测断面进行计算，难以反映行洪区大面积水域真实的水动力特性，模拟精度难以提高。
[0037]
概化计算较好地避免了断面流速分布不均的问题，因为无流速或回流的那部分区域已经通过滞洪节点反映，概化处理的断面只要反映行洪能力即可，经过验证，这种操作方式是可行的。最大的难点在于寻找与实际行洪能力一致的断面型式。首先该断面不能太宽，否则与滞洪节点的库容叠加会高估行洪区的蓄水能力，另一方面，该断面也不能太窄，否则泄水能力与原型不相似。除断面宽度外，需要确定的参数还包括断面深泓、边坡系数等。这些只能通过反复试算确定，耗时漫长，而且主观性很强，没有一定之规。不同的人处理同一行洪区，可能会得出不同的结果，也可能花费很长时间也得不到合适的概化断面。由于此方法的主观性和经验性太强，目前没有标准化的试算流程。
[0038]
这两种方法还有一个共同的缺点：基于断面地形和断面间距的计算模式仅能得到水流因子的断面平均值等有限信息，不能估算行洪区不同区域的淹没风险，功能上存在一定的局限性。
[0039]
(2)整体二维模型是所有方法中精度最高的，但此类方法存在两个问题。一是计算开销太大。随着研究范围扩大，网格密度增加，模型计算的速度显著变慢。目前为止，尚未有淮河所有行洪区同时计算的文献报道。二是对口门(或闸门)的处理不够灵活。行洪区在运用过程中，口门冲宽变深直至最终稳定是一个逐步发展的过程，行洪区的进退洪闸也是逐
步开启，以避免冲刷破坏。实际上，现有大多数二维模型对水工建筑物的处理都比较简单，不够灵活。
[0040]
3)为克服全一维模型以及整体二维模型遇到的技术困难，许多学者分别采用一、二维模型计算狭窄的河网和行洪区大面积水域。其中一、二维模型通过连接关系衔接在一起而实现同步计算。这类模型已广泛应用于规划设计领域，但在实时洪水计算中还存在一些先天不足：一是模型边界条件复杂，一维、二维模块以显式方式衔接，稳定性欠佳；二是计算量仍然较大，难以及时为防洪决策提供信息。


技术实现要素：

[0041]
本发明提出一种适用于含行洪区河网结构的洪水模拟新方法。该方法在保证计算高效的同时，克服现有一维概化计算方案流程不清，试算繁琐、难以应用的缺点。
[0042]
为了达到上述目的，本发明所采用的技术方案为：
[0043]
一种含行洪区河网结构的洪水模拟计算方法，主要包含9个步骤，如下：
[0044]
步骤1，运用gis工具构建行洪区数字高程模型dem，计算行洪区水位-面积-库容曲线；
[0045]
步骤2，基于dem构建行洪区平面二维水动力数学模型，并结合土地类型给单元糙率赋初值；
[0046]
步骤3，构建干流河道的一维水动力数学模型，并给滩槽糙率赋初值；
[0047]
步骤4，耦合一、二维模型，实现河道与行洪区的整体计算；
[0048]
步骤5，通过一、二维水动力数学模型的计算，得出行洪区正常行洪阶段的泄流能力；
[0049]
步骤6：基于步骤2、步骤5成果，构建行洪区一维概化模型；
[0050]
步骤7：通过水工建筑物节点，将行洪区一维概化模型与河道一维模型进行连接，构建含行洪区河网结构的一维水动力模型；
[0051]
步骤8：通过全一维模型计算，获取各节点水位、流量及蓄量信息；
[0052]
步骤9：以步骤8得出行洪区进出流量为点源边界，计算行洪区的洪水淹没要素。
[0053]
所述步骤1主要过程如下：1)由dwg格式的地形图数据转换成arcgis识别的shp格式文件，利用3d分析模块，由矢量等高线、点文件构造tin文件；2)根据地形图的精度，将tin文件插值生成对应分辨率的dem栅格文件；3)基于dem构建数字高程模型，进行不同蓄水位淹没区域的计算；4)统计得出不同水位下的淹没面积，点绘水位-面积-库容曲线。
[0054]
利用实测水文资料开展一、二维耦合模型的验证，若监测点水位的试算值与实测值之间误差足够小，满足用户设定的精度要求，则进入下一步，否则需要重新调整参数，直至满足要求。
[0055]
本发明与现有技术相比的优点如下：
[0056]
(1)与原有的一维模型相比，本发明通过蓄洪节点和行洪节点的设置，首次实现了蓄洪能力与行洪能力的彻底分离，概化流程清晰、人为干扰少，大大降低了人工试算的主观性。同时，以一维计算的成果作为边界，可以单独驱动二维模型，获取更多的洪水淹没信息，成果更加丰富。
[0057]
(2)与原有的二维或一、二维耦合模型相比，本发明采用先粗后细的模型结构，既
能对多方案的洪水风险进行快速动态分析，也方便对重点调度方案进行细化，对防汛应急指挥更具实际意义。
附图说明:
[0058]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图；
[0059]
图1为含行洪区河网结构的洪水模拟计算方法步骤示意图。
[0060]
图2为含行洪区结构河网的一维概化图。
[0061]
图3为先粗后细的模型结构。
[0062]
图4为一、二维模型耦合结构。
[0063]
图5为一、二维耦合模型的边界衔接方法。
[0064]
图6为行洪区一维断面概化。
[0065]
图7为花园湖行洪区库容及泄流能力概化。
[0066]
图8淮河干流水位及花园湖行洪区进、退洪流量比较(水位过程)。
[0067]
图9淮河干流水位及花园湖行洪区进、退洪流量比较(进洪闸流量)。
[0068]
图10淮河干流水位及花园湖行洪区进、退洪流量比较(退洪闸流量)。
[0069]
图11花园湖行洪区进洪2天后淹没水深分布。
具体实施方式:
[0070]
下面结合附图和实施例对本发明进一步说明。
[0071]
本发明提出一种适用于含行洪区河网结构的洪水模拟新方法。该方法在保证计算高效的同时，克服现有一维概化计算方案流程不清，试算繁琐、难以应用的缺点。
[0072]
如图2所示，该一维概化方式针对不同阶段的洪水传播规律，设计了相应的计算方法：当处于充蓄阶段时，洪水由退洪闸倒流至滞蓄节点，这一反向进洪的过程大部分时间内属于动力波范畴，采用完整的圣维南方程(式1-2)描述；随着蓄水位抬升，行洪区进入正常行洪阶段，采用单值化的蓄量与出流关系曲线q＝q(z)，可以近似模拟正常行洪阶段的泄水能力；当水位回落，比降减少，蓄量与出流的关系已比较散乱，此时重新使用完整的动力波方法进行洪水演算。经过上述处理，行洪区的行洪、蓄洪能力不仅可以单独概化，互不干扰。而且正常行洪阶段的q＝q(z)曲线，不需要反复试算就可得出，计算比较简单。
[0073]
整体一维的模型计算后，对于个别的优选调度方案可能需要知道行洪区内部更详细的洪水淹没要素，包括淹没范围、水深、流速、历时等。这就可以将可将一维模型计算出的进、退洪闸流量概化为点源，单独搭建二维水动力数学模型进行模拟，从而实现不同维度模型的单向耦合。这种先粗后细的模型结构(图3)，既有利于多方案的快速寻优，也方便对重点调度方案进行细化。
[0074]
如图1示，本发明提出的含行洪区河网结构的洪水模拟算法主要包含9个步骤，各步原理及实施方法如下：
[0075]
步骤1：构建行洪区数字高程模型计算库容曲线
[0076]
在本步骤内，要先利用行洪区实测散点数据计算行洪区的库容曲线，实现蓄洪能力与行洪能力的分离。
[0077]
主要过程如下：1)由dwg格式的地形图数据转换成arcgis识别的shp格式文件，利用3d分析模块，由矢量等高线(点)文件构造tin文件；2)根据地形图的精度，将tin文件插值生成对应分辨率的dem栅格文件；3)基于dem构建数字高程模型，进行不同蓄水位淹没区域的计算；4)统计得出不同水位下的淹没面积，点绘水位-面积-库容曲线。
[0078]
步骤2：基于步骤1得到的dem，构建行洪区二维水动力模型
[0079]
二维水动力模型的实质就是求解前述浅水方程，即式3～5。步骤1中得到的数字高程模型dem可直接作为模型的计算网格。除此之外，模型搭建还需要给定相关参数，其中比较敏感参数是床面糙率，它直接决定了计算结果的合理性。
[0080]
二维模型的糙率与床面粗糙程度、植被生长、区内涉水建筑物分布等多种因素有关。设定时，可结合天地图等高清卫星图片，依据不同的地表类型划分为若干糙率分区，先通过表1估计行洪区糙率的初值，待耦合一维模型后，再根据实测水位线继续调整，以最终确定行洪区的糙率分布。
[0081]
表1天然河道糙率n值
[0082]
[0083]
[0084][0085]
步骤3：构建干流河道的一维水动力数学模型
[0086]
一维水动力模型的实质就是求解圣维南方程组，即式1～2。方程的求解以水位为隐函数的三级解法，详细步骤可参阅文献：张二俊,张东生,李挺.河网非恒定流三级联合算法[j].华东水利学院学报,1982,10(1):1-13.
[0087]
一维模型采用地形剖分技术，计算的是每个断面水力参数的平均值。除了断面间距外，比较重要的参数依然是河道糙率的布设。一般而言，含行洪区的河道都处于中下游平原区，拥有宽阔的滩地和窄深的主槽，这两部分的糙率是有区别的，需要结合洪中枯不同量级的实测洪水资料，分滩槽、分水位级进行调试率定。
[0088]
步骤4：耦合一、二维模型，实现河道与行洪区的整体计算
[0089]
基于步骤2与步骤3的成果，构建干流一维、行洪区二维水动力数学模型，结构如图4。其中，一维模型还包括与行洪区连接河段，它控制了向行洪区分洪的进、退洪闸门。根据闸孔开度和上下游水位，可将过闸流态分为自由孔流、自由堰流、淹没孔流、淹没堰流四种，每种类型采用对应的水力学公式进行计算。
[0090]
一维与二维的耦合采用显式连接的方式：流量从一维模型的最后一个断面取出，作为源项进入二维模型中。二维模型计算出的水位再返回一维模型作为边界，供下一次循环使用。具体的操作方法见图5。
[0091]
图中预测器半步长的计算，采用了忽略迁移项的简化动量方程，见式(7)。
[0092][0093]
利用实测水文资料开展一、二维耦合模型的验证，若监测点水位的试算值与实测值之间误差足够小，满足用户设定的精度要求，则进入下一步，否则需要重新调整参数，直至满足要求。
[0094]
步骤5：通过一、二维水动力数学模型的计算，得出行洪区正常行洪阶段的泄流能力
[0095]
以步骤4构建的一、二维模型为计算平台，开展洪水的模拟计算。计算过程中，所有行洪区均按现行的调度规则开启使用。为保证行洪区参与行洪，且有较大的分洪比例，上游来水和区间洪水组成可选用设计洪水条件。计算完成后，记录退洪闸计算出流过程，以及进洪闸行洪区侧与退洪闸行洪区侧的计算水位过程。点绘正常行洪阶段行洪区泄流量(退洪闸计算出流量)与代表水位(进洪闸和退洪闸行洪区侧平均值)散点图，并运用最小二乘法
拟合曲线，得到q＝q(z)。
[0096]
步骤6：基于步骤2、步骤5成果，构建行洪区一维概化模型
[0097]
模型的构建思路如图2所示。步骤2提供行洪区概化中蓄洪节点的库容曲线，步骤5提供了行洪区概化中正常行洪时的泄流能力曲线q＝q(z)。行洪区的一维断面可简单概化为由底部窄深三角形和顶部宽浅矩形组成的断面，如图6所示，h1、h2分别对应行洪区的堤顶高程和区域最低点高程，h3是为了防止干河床演进设置的窄缝，取值可略低于淮河干流的起涨水位；b为行洪区的行洪宽度，一般可取进洪闸宽度的2～3倍。
[0098]
步骤7：通过水工建筑物节点，将行洪区一维概化模型与河道一维模型进行连接，构建含行洪区河网结构的一维水动力模型
[0099]
利用行洪区进、退洪闸，将行洪区一维概化模型与步骤3的淮河干流一维模型联合求解，得到含行洪区河网结构全一维水动力数学模型。
[0100]
步骤8：通过全一维模型计算，获取各节点水位、流量及蓄量信息
[0101]
基于步骤7构建的一维概化模型开展洪水的模拟计算。行洪区的启用时机可以采用人机交互或规则调度等方式进行，计算完成后可到各种调度方案下河网各断面的水位、流量过程，以及行洪区的蓄洪、行洪过程。
[0102]
步骤9：以步骤8得出行洪区进出流量为点源边界，计算行洪区的洪水淹没要素
[0103]
为进一步细化个别调度方案，可将步骤8得出的行洪区进、退洪闸流量作为点源边界，单独运行步骤2构建的水动力数学模型，这样就可以获取行洪区内部更详细的水力信息，如范围、水深、流速、历时等淹没要素。
[0104]
具体实施例：
[0105]
选择淮河干流蚌埠至小柳巷段作为本发明的试验区。该河段处于蚌埠以下，全长104km，河道左侧分布有一处花园湖行洪区(编号为8)。该行洪区口门变闸门的提升工程已先期获得批复，并于2018年实施完成。此外，自1952年实施高低水分流工程后，几乎没有支流入汇，河网结构比较简单，加之首末节点均设有水文站控制，观测资料丰富，是测试行洪区降维概化算法比较理想的河段。
[0106]
以visual c#2019为开发环境，建立研究河段整体一维模型。淮河干流一维模型的断面间距控制在500m以内，主槽和滩地糙率取值分别为0.021～0.024、0.035～0.039；图7给出了行洪区库容和泄流能力的概化成果，其中滞蓄节点的库容曲线通过分析dem数据得出，水位与出流关系曲线q＝q(z)利用最小二乘法拟合得出，拟合采用的原始数据取自一、二维耦合模型的计算结果。行洪区概化断面的参数取值为：b＝2000m,h1＝20.5m,h2＝16.5m,h3＝14.0m。
[0107]
该模型的精度采用如下方式进行检验：
[0108]
一是通过对实测洪水的反演，分析干流河道的糙率取值是否恰当。表2为2003年、2020年沿程各站最高水位计算值与实测值的对比，图8为2020年主要站点计算与实测水位过程的比较。由图表可知，模拟的洪水过程与实测过程峰谷对应，吻合甚好，5个站点的水位误差都不超过0.10m，nash效率系数均在0.91以上。
[0109]
二是对行洪区运用进行数值实验，并与一、二维耦合模型的成果相对照，以评估降维概化方式是否合理。由于花园湖行洪区近60年来没有启用记录，根据现行的控制运用办法，假设当临淮关达到19.75m时，启用花园湖行洪区。图9、10给出了该情景下两种模型计算
出的行洪区进、退洪过程。对比可知，行洪区开始分洪的时间相同，分洪的流量过程也基本一致，进、退洪闸最大流量的相对误差仅为0.6％和0.3％。验证结果表明，在行洪区蓄泄能力方面，概化一维模型与一、二维耦合模型的模拟精度相当。
[0110]
表2主要控制站水位洪峰计算值与实测值对比
[0111][0112]
以一维模型计算出的进、退洪流量为驱动项，单独运行花园湖二维模型，就可以得出洪水在行洪区内的演进情况，图11给出了行洪区进洪2天后的淹没水深分布情况。
[0113]
本发明所述的实施示例仅仅是对本发明的优选实施方式进行的描述，并非对本发明构思和范围进行限定，在不脱离本发明设计思想的前提下，本领域中工程技术人员对本发明的技术方案做出的各种变型和改进，均应落入本发明的保护范围，本发明请求保护的技术内容，已经全部记载在权利要求书中。