基于非均匀网格的高效高精度洪涝过程模拟

侯精明，王俊珲*，同 玉，张兆安，郭敏鹏，苏 锋，高徐军

(1.西安理工大学 西北旱区生态水利国家重点实验室，陕西 西安 710048；2.中国电建集团 西北勘测设计研究院有限公司，陕西 西安 710065)

在全球气候变暖背景下，极端暴雨发生的频率和强度逐渐增加[1]，每年因极端暴雨洪涝灾害造成了巨大的人员伤亡和财产损失，严重制约经济社会的持续健康发展[2-3]。2016年汛期中国暴雨洪涝灾害频发，长江流域发生特大洪水，全国受灾人口达3 282万，直接经济损失约1 470亿元，与2000年以来同期均值相比，直接经济损失偏多51%[4]，水安全形式十分严峻。这就需要利用高精度数值模型研究洪水产生、演进及消退过程，为政府决策部门提供行之有效的科学依据[5]。

数值模型在预测地表水流动及相关过程中起着重要作用，其本质是对非线性微分方程组进行求解，其核心在于解的离散。解的离散程度会产生奇异性的解，奇异性使数值模拟受到阻碍。而网格的生成主导解的离散过程，因此一个高质量的计算网格对模型的模拟精度和计算效率和有着很大的影响。目前，通常使用两种网格对计算域进行离散，即结构网格与非结构网格。非结构网格在处理复杂边界问题时具有明显的优势，但非结构化网格关联的数据结构更复杂，控制方程的离散化更麻烦[6]。结构化网格结构简单，数值实现方便，且计算精度更高，已被众多浅水水流模型采用[7-8]。

然而，在实际的洪水模拟过程中，计算域内往往包含地形变化显著的狭窄河道区域及面积广阔的泛洪区[9]，且模拟通常涉及较大的研究区域。在均匀网格上，大范围的高精度数值模拟不可避免会涉及大量计算节点，导致计算效率降低，而非均匀网格可以很好适用于具有复杂地形特征的区域。对地形变化明显以及相对重要的区域提供精细网格，而其余区域则由粗网格离散，从而避免非必要的计算，减少计算成本。

近年来，随着计算机技术的发展，网格划分技术趋于成熟，周建中等[9]将自适应网格技术应用于溃坝洪水演进模拟中，计算效率有一定提高，然而自适应网格可能无法同时满足质量守恒和C-property[10]。质量守恒影响水量，进而影响水位。C-property通常反映模型保持所谓的良好平衡状态的能力，即通量项和坡源项之间是平衡的，从而不会出现假动量[6]；Liang[11]提出了非均匀网格方法以模拟洪水演进过程，模拟精度高，但研究区域土地利用单一，且算例多为理想算例，不适应于土地利用复杂（常为5种不同土地利用类型）、实际大规模范围的城市发达地区洪涝模拟。基于此，本文提出一种新的网格划分方法，分别对不同土地利用设置不同划分准则，以适应实际大规模区域的洪涝模拟。此外，静态非均匀网格可同时满足质量守恒与C-property[6]，对保证模型的数值稳定性具有重要意义。

本文采用一种具有层次拓扑关系的、结构化的静态非均匀网格进行浅层水流模拟的网格优化，以高程梯度变化为判定变量调整计算域内网格密度，在地势平坦区域适当加大网格尺寸，在水流模式受复杂地形特征影响的区域内适当加密网格并进行高分辨率计算。与基于GPU加速技术的高分辨数值模型相结合，从而在保持解的精度基础上，降低计算成本，提高模型的计算效率，以提升其在实际应用中的性能。

1 非均匀网格划分

本文采用递归层次的结构网格（图1），网格大小以网格划分层级level反映，（记为lev），叶网格（如网格（i-1，j）和（i，j））划分层级lev为0，其大小为Δx0和Δy0。对叶网格进行细化可以得到子网格，子网格大小为Δx=Δx0/2lev，Δy=Δy0/2lev。为降低模型复杂性，网格与其邻居网格大小满足2∶1规则，不允许任何网格拥有大于或小于两倍的邻居，采用网格之间的相对关系来确定存储邻居网格信息[6]。

图1 结构化非均匀网格图Fig.1 Structured non-uniform grid

生成精确反映区域地形特征的优化的非均匀网格，一般遵循6个步骤：

1）定义一个初始均匀细网格来反映所考虑区域的地形数据的最高分辨率，网格大小为Δx，如图2所示；在初始精细网格的基础上设置虚拟背景粗网格，其大小为初始均匀细网格的整数倍，如图3所示。一般粗化水平lev分为2级，最大粗化水平levm由式（1）计算：

图2 根据地形生成的初始精细网格Fig.2 Initial fine grid generated from topographic data

图3 定义levm=2的背景粗网格Fig.3 Defined background coarse grid with levm=2

式中，Δx、Δx0分别为初始网格和虚拟背景粗网格的单元大小。

2）对于任意背景粗网格，计算其内部所有初始细网格的高程梯度（高程的1阶导数）。对于初始细网格（i，j），其x和y方向的高程梯度与河床高程平均梯度分别为：

式中，zb为河床高程，Δx、Δy为初始细网格的大小，Grad_x、Grad_y分别为x和y方向的高程梯度，Grad为河床高程平均梯度。

3）对于任意背景粗网格，计算其内部所有初始细网格的高程梯度变化（高程的2阶导数），对于初始细网格（i,j），其x、y方向的高程梯度变化与河床高程平均梯度变化分别为：

式中，Gvar_x、Gvar_y分别为x和y方向的高程梯度变化，Gvar为河床高程平均梯度变化。

4）对于任意背景粗网格，取背景粗网格中所有初始精细网格的高程梯度变化最大值Gvar_max作为该背景粗网格的坡度变化值GvarΦ；取背景粗网格中所有初始精细网格对应的土地利用类型数量最多者作为该背景粗网格的土地利用类型。

5）根据研究区域内不同土地利用类型（即：建设用地、道路、裸地、林地和草地）重要程度设置不同阈值优化网格（若仅有一种土地利用，则仅设一种阈值）。对于任一土地利用类型，设置两层阈值grcr1和grcr2，即网格最高粗化水平为2级，当背景粗网格的梯度变化值GvarΦ大于阈值grcr1，初始精细网格第一次粗化；当背景粗网格的梯度变化值GvarΦ大于阈值grcr2，初始精细网格第二次粗化。否则，网格不需粗化。此外，可根据需求特定某区域，对该区域网格单独细化。

6）引入2∶1规则保证网格过渡平稳，即网格大小必须为邻居网格的2倍、1倍或1/2倍。在此基础上，定义网格顺序编号、连接规则等[6]。最终，生成优化后的非均匀笛卡尔网格，该网格邻居信息完全由简单代数关系决定，避免使用显式的数据结构来存储邻居信息[12]，从而使存储需求最小化。网格优化情况如图4所示。

图4 优化后的结构化非均匀网格Fig.4 Optimized structured non-uniform grid

2 数值求解

2.1 控制方程

本文利用高精度数值模型模拟研究区域内的洪涝灾害过程，模型所使用的控制方程为平面2维浅水方程，忽略运动黏性项、紊流黏性项、风应力和科氏力的2维非线性浅水方程矢量形式如下：

式中：x、y、t分别为笛卡尔空间坐标与时间坐标；h为水深；qx、qy为单宽流量；u、v表示x、y方向上流速；F、G表示x、y方向上通量矢量；S为源项矢量，包括降雨或下渗源项i、底坡源项及摩阻力源项；zb为河床底面高程；Cf为河床糙率系数，Cf=gn2/h1/3，n为曼宁系数。

2.2 数值求解

2.2.1 数值离散

采用中心格式的有限体积法，将变量存于网格中心，将式（4）在控制体 Ω上积分得：

式中，Г为控制体边界，n为垂直于面元dГ的单位向量，F(q)n为相应界面的通量。在笛卡尔坐标系下，式（7）中的曲面积分项可以转化为：

式中， Δx与 Δy分别为笛卡尔坐标系下的网格单元（i，j）的左右和上下边长，FE、FW、FN和FS分别为网格单元（i，j）的东西南北4个界面的界面通量。

不同于均匀网格，非均匀网格的单元尺寸存在差异，如图5所示。单元ic的东界面的通量由FE=FE1+FE2计算得出，其中FE1、FE2是通过两个较小单元的西界面的中点的通量。考虑以i为中心的较小单元，通量计算可能涉及位置w、i和e处的流量变量。由于守恒的流量变量h，qx和qy存储在单元中心，因此，i处的流量变量可直接获得，w和e处流量变量可通过自然邻居内插法获得。

图5 非均匀网格中通过控制体积边界的流量Fig.5 Fluxes through the control volume boundary under the non-uniform grid

对界面通量进行计算采用近似Godunov方法的HLLC近似Riemann求解[13]，该格式在处理干单元时的功能要优于其他格式，能较好地处理复杂地形下干湿交替变化，具有较强的激波捕捉能力。其界面通量求解过程为：

式中：SL、SR、SM分别为计算单元左、右两侧及中间的波速，具体计算见式(10)。当SL≥0和SR≤0时，计算单元界面的通量值分别由左、右两侧的水力要素确定；当0≤SL≤SM和SM≤0≤SR时，计算单元界面的中间对流通量值由HLLC近似Riemann解给出，通量值详细计算见参考文献[13]。

式中，uL、uR分别为计算单元左、右两侧的流速，hL、hR分别为计算单元左、右两侧的水深，中间变量h*与u*由下式计算。

2.2.2 2阶数值重构

考虑到采用基于HLLC近似Riemann求解界面通量时在空间上仅具有1阶精度，为使数值解的空间精度提高到2阶，同时为了避免数值重构后，在水面梯度大的地方出现数值振荡的现象，模型采用TVDMUSCL方法进行数值重构[11]，表达式为：

其中：

式中，U为计算变量，φ为限制器函数，本模型采用Min mod限制器，其函数表达式为：

此外，为保证数值解整体上提高到2阶精度，同时维持数值解的稳定性，对时间步长采用两阶Runge-Kutta格式[14]来获得2阶精度。干湿动边界处理方面，水深容差为0.000 001 m以区别干湿网格单元[15]。底坡项使用模型作者提出的底坡通量法处理，摩阻项的计算使用分裂点隐式法以提高稳定性[16]，计算过程中，Courant取0.5。

2.2.3 GPU加速方法

考虑到洪水过程模拟尺度较大，模型利用GPU并行计算技术，通过CUDA语言，将求解的浅水方程各项在GPU上运行，在空间上实现大规模计算各网格水力要素信息；采用C++语言，将数据的读写、变量初始化在CPU上运行。从而，在不降低精度条件下实现高速运算[17]。图6为GPU计算的流程图。图7为基于非均匀网格的数值模型进行洪水模拟流程图。

图6 GPU计算流程图Fig.6 Flowchart of the computation using GPU

图7 基于非均匀网格的洪涝模拟流程图Fig.7 Flowchart of flood simulation with non-uniform grids

3 非均匀网格在洪涝模拟中应用

3.1 3个驼峰溃坝波模拟

本文选择3个驼峰溃坝波模拟的理想算例来验证模型稳定性。该溃坝问题被广泛应用于测试模型在2维复杂地形下处理非恒定流的能力[18]，具有一定的代表性。研究区域长75 m，宽30 m，糙率n=0.018 s/m1/3，底部高程定义如下：

初始均匀网格分辨率为0.5 m，共计150×60个方形单元，以此为基础进行非均匀网格划分，网格最高粗化水平为2级，设置两层阈值：grcr1=0.001，grcr2=0.011。经非均匀网格系统划分，优化后的非均匀网格数为3 216，研究区域地形及非均匀网格分布如图8和9所示。

图8 3个驼峰区域地形图Fig.8 Topographic map of three humps

图9 3个驼峰地形非均匀网格分布图Fig.9 Non-uniform grid based on three hump terrain

算例共分两部分：首先，验证和谐性，在整个计算域内设初始水位为1.875 m，1 000 s后，计算域内水位保持不变，表明模型具有良好的静水和谐性。其次，修改初始条件，在x=16 m处建坝，大坝上游初始水位为1.875 m，下游干河床。t=0 s时，大坝瞬时全溃，模拟运行300 s。此外，为验证水量平衡情况，地形四周设置为闭边界，且无下渗情况，初始水量为1 125 m3，模拟300 s后，水量基本保持不变，因此，水量平衡状况良好。图10为P1、P2和P3这3个点位基于均匀网格与非均匀网格模拟的水位情况。

由图10可看出，P1、P2和P3这3个测点基于非均匀网格模拟水位变化情况与基于均匀网格模拟水位情况吻合度较高，变化趋势一致，均方根误差RMSE分别为2.24%、2.23%、1.15%。因此，基于非均匀网格的数值模型具有较高的模拟精度。图11为基于非均匀网格模拟的t=6 s和12 s的2维与3维水深图。

图10 P1、P2、P3测点的水位模拟情况Fig.10 Simulated water level at gauge P1, P2, P3

图11 不同时刻的2维与3维水深图Fig.11 3D and 2D water depth maps at different times

考虑到地形关于直线y=15 m对称，理论上结果应保持对称性，图11呈现的水深模拟结果具有较好的对称性，符合地表水流运动规律。在溃坝水波演进过程中，水波与主峰及小峰之间不断接触碰撞，t=12 s水流绕过主峰在下游形成干湿界面，最终在经过多次碰撞及摩擦阻力的影响下，水流运动趋于平稳。

本算例证明，模型有效捕获由3个驼峰和区域边界相互作用的激波引起的复杂流动以及干湿界面及水位敏感区，且在干湿边界附近未发生明显的数值扩散，具有良好的和谐性。此外，基于非均匀网格的模拟计算可以提高模型模拟效率，模拟采用PC机，搭载显卡NVDIA RTX2070执行，基于均匀网格的GPU加速的数值模型运行26.95 s可完成300 s溃坝水流演进过程，采用基于非均匀网格进行计算，模型运行时间为10.43 s，计算效率约为均匀网格的2.7倍。

3.2 Toce溃坝波模拟

本模型为Toce河物理模型，该区域溃坝波洪水演进模拟[12]常用于验证水动力学数值模型的模拟精度和数值稳定性，具有一定典型性。该物理模型长约50 m，宽约10 m，中部有一空水库。初始均匀网格分辨率为0.05 m，共计约21万个方形单元，以初始均匀网格为基础进行非均匀网格优化，网格最高粗化水平为2级，设置两层阈值grcr1 = 0.09，grcr2 = 0.12。经非均匀网格系统划分，优化后的非均匀网格数为87 030，研究区域地形及非均匀网格分布如图12和13所示。

图12 Toce河区域地形图Fig.12 Topographic map of Toce River

图13 Toce地形非均匀网格分布图Fig.13 Non-uniform grid based on Toce

利用数值模型分别基于非均匀网格及均匀网格对溃坝波洪水演进过程进行数值模拟，干湿水深的网格单元判定条件为0.000 001 m，Courant设置为0.5，输入条件包括地形文件、入流文件等。左侧设为入流边界，右侧为出流开边界，其余设置为闭边界[12]。图14为该研究区域入流流量过程线。图15为P4、P9和P21这3个点位基于数值模拟及实测水位情况。此外，为验证水量平衡情况，除左侧入流边界外，其余边界设置为闭边界，研究区域无下渗，入流过程总水量约为18 m3，模拟180 s后，研究区域总水量基本保持不变，因此，水量平衡状况良好。

图14 入流流量过程Fig.14 Inflow flow process data

图15 P4、P9、P21测点基于模拟及实测水位情况Fig.15 Simulated and measured water level at gauge P4,P9, P21

由图15可看出，P4、P9、P21这3个测点基于非均匀网格模拟水位变化与均匀网格模拟情况及实测水位变化吻合度较高。与实测值相比，非均匀网格模拟水位情况的纳什效率系数分别为2.3%、1.0%和0.9%。因此，基于非均匀网格的数值模型可很好模拟溃坝波洪水演进情况。此外，由于入流边界在西侧，溃坝水流的洪水演进过程由实际地形的西侧逐渐经过P4、P9、P21测点，导致P4、P9、P21这3个测点实测水位在同一时刻数值不同，即峰值水位逐渐降低，峰现时间逐渐延后；而P4、P9这2测点地形较高且靠近入流边界，P21测点地势较平坦，因此，P4、P9测点水位与入流过程趋势基本相同，呈现先增后减趋势，而P21测点水位达到峰值后逐渐趋于平缓。图16为不同时刻洪水演进过程。

图16 t=40 s和120 s时溃坝水流演进过程分布Fig.16 Distribution of dam break flow evolution process at t=40 s and 120 s

此外，基于非均匀网格的模拟计算可以提高模型计算效率。模拟采用PC机，搭载显卡NVDIA RTX2070执行，基于均匀网格的GPU加速的数值模型需要耗时145 s才能完成洪水演进过程，相同条件下，基于非均匀网格则仅需约为64 s，计算效率约为基于均匀网格模拟的2.3倍。

3.3 城市内涝模拟

本文以陕西省西咸新区沣西新城为研究区域，降雨多集中在7至9月，且该区域土地利用类型多样、建筑物林立、下垫面情况复杂[19]，选此研究区域模拟城市内涝过程具有一定的代表性，图17为研究区域地理位置概况图。

图17 研究区域地理位置概况图Fig.17 Overview of the geographical location of the study area

利用机载激光雷达技术对研究区域地形进行测算，得到分辨率为1 m共计312万均匀网格单元的DEM数据，图18为研究区域数字高程图、正射影像图。

图18 研究区域数字高程及正射影像图Fig.18 DEM date and orthophoto image map in the area

考虑到计算区域范围较大，初始均匀网格较多，导致计算量大、耗时长、成本较高。故以初始化均匀网格为基础，采用优化的非均匀网格系统对计算域内不同的土地利用类型分别进行划分。网格最高粗化水平设置为2级，对于5种不同的土地利用类型（建设用地、道路、裸地、林地和草地）分别设置两层不同的阈值grcr1与grcr2。

对于高程相对较低的草地、裸地、道路设置同一种阈值，本文首先将阈值设置为grcr1=grcr1草地=grcr1裸地=grcr1道路=0.50，grcr2=grcr2草地=grcr2裸地=grcr2道路=0.55，当网格上的高程变化大于grcr1时，网格第1次粗化，当网格上的高程变化大于grcr2时，网格第2次粗化，达到最高粗化水平（叶网格水平）；对于林地、建设用地设置另一种阈值，首先将阈值设置为grcr1=grcr1林地=grcr1建设用地=0.65，grcr2=grcr2林地=grcr2建设用地=0.70，当网格上的高程变化大于grcr1时，网格粗化1次，当网格上的高程变化大于grcr2时，网格第2次粗化，达到最高粗化水平。经划分，非均匀网格数约为130万。

采用高分辨率数值模型分别基于初始均匀网格及该阈值设置条件下的非均匀计算网格对研究区域的内涝过程进行精细化模拟，输入条件包括地形文件、降雨数据、土地利用文件等，四周为开边界且无入流，计算过程中库朗数取0.5。其中，地形文件分别为初始均匀方形网格的DEM数据及非均匀计算网格数据文件；降雨数据为2016年8月25日西咸新区云谷10号楼气象站实测降雨，累计降雨量66 mm，降雨过程见图19。

图19 2016-08-25场次实测降雨Fig.19 2016-08-25 measured rainfall

根据研究区正射影像图，采用最大似然分类法分别对构建的初始均匀方形网格单元及非均匀网格单元划分成5种不同的土地利用（建设用地、道路、裸地、林地和草地），并制作成土地利用文件。不同土地利用的稳定下渗率及曼宁系数见参考文献[20]。

利用模型模拟计算域从开始降雨到7.5 h的积水过程。可根据模拟结果获取到该区域不同时刻的内涝点位置分布、积水面积、积水深度等内涝情况。由模拟结果可知，t= 4 h时，该区域内涝积水量为最大值，图20为基于均匀网格与非均匀网格划分的模拟结果积水量最大时分布情况。图20中共标记出内涝积水量最大时的4处积水区（如图20中椭圆圈中区域所示），其内涝影响较为严重。

图20 基于非均匀与均匀网格的积水点位模拟情况Fig.20 Simulated inundation based on the uniform grid and the non-uniform grid

模拟过程中基于非均匀网格的积水位置与均匀网格及实际情况的发生位置大致吻合（4处）。将基于计算域划分的均匀网格的模拟积水程度与非均匀网格的模拟积水程度与实测数据[19]进行对比，对比结果如表1所示。其中，实测数据由陕西省沣西新城管委会实地测量获得。

由表1可知，基于结构化非均匀网格的模拟与在分辨率为1 m的均匀网格上生成的模拟及实际水位情况相匹配。与均匀网格及实际情况相比，整个计算域内4处积水位置的积水面积加权平均误差分别为2.94%、8.95%，表明基于非均匀模拟精度与均匀网格模拟精度及实际情况相似。此外，本次模拟在模拟采用PC机，搭载专业计算显卡NVDIA Tesla P100执行，均匀计算网格模拟共耗时约13 300 s，非均匀网格模拟共耗时约5 500 s。速度约为原来的2.4倍。可见，基于结构化非均匀网格的GPU加速的浅水流动模型在模拟内涝过程中可在保证准确度的基础上可有效提高模型的计算效率，缩短模拟时间。

表1 基于非均匀与均匀网格模拟水位及实际对比Tab. 1 Comparison of the field data and the simulated water level based on the non-uniform and uniform grid

4 结 论

计算域内地形的概化精度对模拟结果有着很大影响，本文建立了基于非均匀网格的GPU加速的数值模型以求解浅水方程。并将新模型应用到理想化地形的溃坝算例和实际大尺度地形的城市内涝算例中，得到以下结论：

1）采用具有层次拓扑关系的非均匀网格划分技术，以地形高程梯度变化为判定变量，根据研究区域不同土地利用类型制定不同划分准则以自动调整网格密度，对地形变化明显及水力要素敏感区域提供精细网格，而其余区域由粗网格离散，从而适应实际洪水模拟中研究区域复杂地形特征的需求。

2）对高分辨率数值模型采用GPU并行计算技术，在不降低模拟精度的条件下，大幅提升模型计算效率。并与非均匀网格划分模型相结合以模拟地表水流动过程。由于相对重要的区域仍被较准确描述，但网格数量大大减少，模拟速度在GPU加速的基础上进一步提高。因此，新模型可在保证模拟精度的同时，进一步提高模型的计算效率。

3）将新模型应用到实例中，在3个驼峰及Toce溃坝模拟中，模型可有效捕获地形和边界相互作用的激波引起的复杂流动以及干湿界面及水位敏感区，验证了模型稳定性，且模拟速度分别为基于均匀网格模拟的2.7与2.3倍；在沣西城市内涝模拟中，基于非均匀网格模拟的积水位置与均匀网格模拟位置及实测情况大致吻合，积水面积加权平均误差分别为2.94%、8.95%，速度约为均匀网格模拟的2.4倍，进一步证明基于非均匀网格模型的高效性。

本研究提出的非均匀网格方法可描述区域的复杂地形特征，并与高分辨率数值模型结合以可靠地求解浅水方程，运算速度在GPU加速基础上进一步提高。与基于均匀网格的模型相比，模型具有相似的求解精度，但计算成本较低，计算效率明显提高。因此，在实际浅水流动模拟中具有更好的潜力，可有效解决大规模洪水过程模拟遇到的精度低、效率慢等问题，适用于溃坝波洪水演进过程及土地利用复杂多样的发达地区的城市内涝过程模拟。