周跃华，岳志春，，潘汀超，Rizwan Qadir

(1.宁夏回族自治区水旱灾害防御中心，宁夏 银川 750001；2.天津大学水利工程仿真与安全国家重点实验室，天津 300072)

0 引 言

黄河宁夏青铜峡到石嘴山河段洪水灾害易发。20世纪青铜峡水文站有记载的洪峰流量大于5 000 m3/s的洪水共5次，洪灾严重威胁了沿河两岸广大人民群众的生命财产安全。

洪水演进数值模拟是开展防洪保护区洪水风险分析的主要方法。王扬等[1]基于一维、二维水力学耦合模型，建立了韩江南北堤防洪保护区溃坝洪水演进数学模型。陈平等[2]构建二维非恒定流数学模型模拟大陆泽及宁晋泊蓄滞洪区50年一遇洪水演进过程。袁水龙等[3]耦合分布式水文模型MIKE SHE和一维水动力模型MIKE 11研究淤地坝建设对黄土高原小流域暴雨洪水过程的影响。苑希民等[4]针对山洪沟防洪标准偏低、遭遇暴雨洪水时溃漫堤风险较大且较难准确预测的问题，系统构建溃漫堤洪水多维耦联数值仿真模型。李大鸣等[5]以二维非恒定流控制方程为理论基础，采用有限体积法，将改进的水量平衡模式应用于山区洪水演进计算，建立了洪水演进数学模型。郭立兵等[6]构建山洪沟道防冲建筑物壅水风险评估模型，能够较好地模拟宁夏北武当沟山洪演进过程和防冲建筑物壅水风险。田福昌等[7]建立山洪沟道溃堤洪水演进一、二维水动力耦合数值模拟模型，分析评估溃堤山洪淹没风险。针对防洪保护区河道洪水的漫滩、溃堤及其与内涝叠加耦合预测的问题，张靖雨等[8]构建区域内涝模型叠加联用的一、二维水动力耦合模型，模拟中运河南片防洪保护区内洪水演进。苑希民等[9]建立了模拟河道和灌区洪水演进过程的漫溃堤洪水全二维水动力模型。戚蓝等[10]基于高精度DEM数据建立三亚市主城区暴雨、洪水与风暴潮多元耦合精细化洪涝分析模型。

针对狭长河段计算区域采用二维水动力模型模拟计算效率低的问题，将河道及防洪保护区的洪水演进概化为一维流动，既可以反映漫溢洪水的特征，又可以减少网格数量提高模型效率。本文在前人研究成果的基础上，采用水力学方法对黄河青石段建立一维非均匀流水动力模型，以陶乐防洪保护区为研究对象，对陶乐计算区段河道进行网格加密处理，计算各断面最高水位，连成水面线；同时结合保护区地形计算保护区淹没水深，为漫溢洪水演进仿真模型构建、防洪保护区洪水风险分析提供一定的参考。

1 模型原理

1.1 河道一维水动力学模型

河道水流流动可采用基于Saint-Venant方程[11-13]构建的一维水力学模型模拟。Saint-Venant方程组为

(1)

(2)

式中，Q为流量；A为断面面积；q为旁侧入流；t为时间步长；x为空间步长；g为重力加速度；β为流速分布系数；Z为水面高程；Sf为摩阻损失。

1.2 离散求解

采用 Abbott 六点隐式差分格式来求解Saint-Venant方程组，即

(3)

(4)

式中，h为水位点；Q为流量点；j为计算断面号；n为时间步。Abbott六点隐式差分格式具体求解过程见参考文献[14]，该格式具有计算稳定、精度高、可靠性强等优点。

1.3 漫溢洪水风险评估方法

对于滩地漫溢洪水淹没的淹没范围和水深分布，可利用数值模型计算出的洪水高程作为洪水水位开展淹没分析[15]。具体来说就是基于河道一维水动力模型计算洪水演进过程，然后依据一维水动力模型模拟获得的各断面最高水位开展漫溢洪水影响评估。

首先，利用地理信息技术，在河道实测断面的基础上加密，从DEM数据提取河岸高程，将河道断面线向保护区外延，形成加密的保护区断面线，保护区断面线与河道断面线一一对应。继而从保护区DEM数据提取河道和滩区离散点高程，整理为河道一维水动力模型计算需要的断面数据，同一个断面线上的离散点最高水位通过一维水动力模型求解。

其次，对计算区域进行二维网格剖分，由各离散点的最高水位数据内插得到各网格水深数据。河道一维水动力模型的计算结果是各个断面上离散点的最高水位，难以直观体现洪水的淹没范围和水深分布。需要通过空间插值方法，将各个断面上离散点的最高水位插值生成连续的洪水水位高程面。自然邻点插值方法是一种基于空间自相关性的面积权重线性内插法[16]。其原理是对于一组泰森多边形，当在数据集中加入一个新的数据点时，就会修改这些泰森多边形。邻点的权重将决定待插点的权重，待插点的权重和目标泰森多边形成比例，即

(5)

式中，h(x)为待插点x处的值；αi为自然邻点xi的权重；h(xi)为自然邻点xi处的值。采用自然邻点插值法对计算网格进行插值，具有计算稳定、精度高、速度快等优点，能够快速准确地反映洪水在保护区内最高水位的分布情况，适用于洪水风险应急分析。

最后，将插值得到的各网格洪水水面高程与研究区的地面高程相减得到洪水淹没水深分布，即

h=hw-hg(hw>hg)

(6)

式中，h为淹没水深；hw为水面高程插值结果；hg表示地面高程。利用GIS平台生成保护区淹没水深图。

2 应用实例

2.1 研究区域概况

黄河是世界著名的天然多沙河流，全长约5 464 km，流域面积约789.6 km2，黄河宁夏段自中卫市南长滩入境至石嘴山市麻黄沟出境，河段长度为397 km，区间流域面积为520 km2，穿越11个市、县(区)，地理位置如图1所示。

图1 黄河宁夏段地理位置示意

随着黄河宁夏段上游水库的建成使用和气候条件变化，该河段水沙关系发生变异并表现为非协调性变化。20世纪80年代以前，黄河宁夏段多年冲淤基本平衡，主河槽能够保持一定的泄洪输沙能力，随后在人类活动与自然气候变化的双重影响下，黄河宁夏段上游来水量、汛期输沙量和造床洪峰流量均大幅度减少，中常洪水过程加长使泥沙主要淤积在主河槽内，使得主河槽淤积萎缩日趋严重，平滩流量逐年减少，主河槽输水输沙能力严重下降，滩地漫溢淹没现象频发，威胁两岸人民生命财产安全。陶乐防洪保护区位于黄河宁夏段下游，面积约182km2，保护区西靠黄河，东临鄂尔多斯台地，地形狭长，类似于峡谷地区河流防洪保护区。

本研究需要收集的数据主要包括研究区域内的基础地理资料、历史洪水灾害资料、河道断面资料、1∶10 000地形图(DLG)、1∶10 000数字高程模型(DEM)、土地利用图或近期遥感影像资料、水利普查数据、黄河大断面数据、沟渠资料等。对地形图进行配准后，叠加保护区的行政边界、水系、防洪工程等图层得到洪水风险图的底图。

2.2 模型构建

采用水力学方法对陶乐防洪保护区进行洪水分析计算，建模流程如图2所示。

图2 建模流程

2.2.1河道一维水动力学模型构建

为了保证计算精度，一维水动力学模型计算断面在2012年实测断面的基础上按间距200 m进行加密。将河道断面线向保护区外延，形成间距为200 m左右的保护区断面线，保护区断面线与河道断面线一一对应，通过建立一维水动力学模型，计算河道各断面最高水位，从而求得保护区各断面最高水位。陶乐保护区断面线的位置如图3所示。

图3 陶乐保护区断面位置示意

根据资料统计，近40年来，宁夏出现具有代表性的洪水年份为1981年和2012年。1981年洪水，下河沿站测得最大洪峰流量为5 880 m3/s，青铜峡站测得最大洪峰流量为6 040 m3/s，石嘴山站测得最大洪峰流量为5 660 m3/s；2012年第3号洪水，下河沿站最大洪峰流量3 520 m3/s，石嘴山站最大洪峰流量3 400 m3/s，青铜峡最大洪峰流量为3 070 m3/s，洪水持续71天。2012年第3号洪水主要是受黄河上游持续降雨及水库下泄等影响，加之近年来黄河宁夏段大洪水极少，河道萎缩淤积严重，导致此次洪水漫滩几率增加，滩地淤积，主槽冲刷。综合考虑河道变化及上游水库调度的影响，本研究选择2012年第3号洪水用于率定模型。

模型计算时间步长为5 s，各断面水深初始值为起算时刻各断面水深，设定有利于计算稳定和结果精度。结果输出时间步长为1 h，方便结果展示。

对宁夏黄河2012年第3号洪水过程(2012年7月22日～2012年9月20日)进行模拟，提取相关断面水位流量数据，根据2012年大断面实测最高水位进行模型验证。上游入流边界采用青铜峡水文站2012年第3号洪水过程进行模拟水位过程曲线(见图4)；下游出流控制边界采用石嘴山水文站水位-流量关系曲线(见图5)。利用2012年遥感影像实测最高水位与模拟45个大断面最高水位对比，模型率定结果见图6，其中44个断面其中误差在0.2 m以内，满足洪水风险分析精度要求，计算精度保证率高达97.8%。

图4 青铜峡水文站2012年第3号洪水过程曲线

图5 石嘴山水文站水位—流量关系曲线

图6 断面率定结果

用模型对青铜峡水文站100年一遇洪水进行分析。上游入流边界为青铜峡水文站100年一遇流量过程曲线(见图7)；下游边界条件为石嘴山水位—流量关系曲线(见图4)。陶乐防洪保护区滩区初始视为无积水，按照河道断面水面线向外延伸确定淹没范围和淹没水深。

图7 青铜峡水文站100年一遇流量过程曲线

2.2.2防洪保护区水深内插

陶乐防洪保护区洪水风险分析依据青石段一维水动力模型计算获得的各断面最高水位开展。计算区二维淹没统计时采用的边界数据为一维模型河道水面线计算结果。根据水动力模型计算结果，利用自然邻点法插值计算得到各网格最高水位，将其与保护区地面高程叠加分析，得到洪水的淹没范围与淹没水深分布。

采用离散求解的方法求解陶乐计算区二维水深。将200 m左右间距的计算区断面线离散成50 m间距的离散点，同一个断面线上的离散点最高水位通过一维水动力模型求解。对陶乐计算区进行二维网格剖分，最大网格面积按1 000 m2进行控制，各网格水深数据由各离散点的最大水深数据内插得到。保护区网格如图8所示。各网格地面高程通过保护区地形数据提取得到，根据各网格最高水位和地面高程做差求得防洪保护区最大淹没水深分布，结果大于0的网格即被洪水淹没的区域，得到的每个网格的值即为淹没水深，得到淹没水深分布图层。将其与工作底图叠加和渲染后，绘制得到陶乐防洪保护区100年一遇洪水淹没水深分布图。

图8 陶乐计算区网格划分示意

2.3 计算结果分析

陶乐防洪保护区靠黄河侧无堤保护，为狭长的带状区域，地势较低的区段主要有红崖子乡河段上游部分、五堆子乡附近河段靠下游部分、陶乐镇附近河段靠下游部分、高仁镇河段上游部分，洪水分析计算结果显示，地势较低的部分淹没范围和淹没水深较大，符合实际情况。保护区内线状建筑物主要为省道S203，沿程高程均高于河道水位，挡水作用明显。洪水分析计算结果显示道路起到了明显挡水作用，如图9所示，计算结果符合实际情况，满足合理性要求。当陶乐计算区遭遇100年一遇洪水时，洪水淹没总人口数量约0.19万，主要影响月牙湖乡、高仁乡、陶乐镇、红崖子乡，总淹没面积约37.5 km2。

图9 陶乐防洪保护区100年一遇洪水淹没水深分布

3 结 论

基于水动力学理论构建了宁夏黄河青铜峡站至石嘴山站河道一维非恒定流水动力模型，采用所建模型计算得到陶乐防洪保护区黄河河道最高水面线，利用GIS平台按洪水淹没模拟分析结果确定陶乐计算区100年一遇滩地漫溢洪水的淹没水深及淹没范围。基于洪水分析计算结果，按照国家相关技术文件规定，绘制了标准、统一、规范的洪水淹没水深图。

研究成果表明，在宁夏黄河河床淤积抬高的情况下，遭遇洪水水位抬升造成漫滩现象风险较大，陶乐防洪保护区受灾较为严重的区域为芦苇、王家沟、小红柳滩、中滩村、下八顷村和黄泥岗等村庄。本文所建模型精度较高，分析结果可为黄河宁夏段整治工程规划设计及防汛指挥决策提供支持，文中采用的GIS与水动力模型融合方法可为其他区域淹没风险分析提供技术参考与应用借鉴。