基于大气-水动力集成模型的城市内涝模拟预报研究

陈光照，侯精明，同 玉，周 聂，高徐军，苏 锋，李继成，吕 鹏，杨 霄，张宏芳

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

近年来，随着城市化率的提高和气候变化的加剧，城市内涝灾害愈加频繁[1]。在下垫面硬化、排水能力不足、短历时强降水事件频发等因素的共同作用下，许多城市将遭受内涝灾害并造成严重损失[2]。政府间气候变化专门委员会(IPCC)的研究表明，1981年至2000年间，500 mm极端降雨的发生概率约为1%，至2100年这一概率将高达18%[3]。因此，开展城市内涝预报方法研究，对防涝减灾具有重要意义。

城市内涝预报是降低内涝灾害风险、减轻灾害损失的重要措施，但目前仍未开发出完全满足防灾减灾要求的预报系统[4]。影响预报效果的问题主要为预见期短及无法精确描述内涝变化过程。为了延长预见期，内涝预报可采用数值天气预报(NWP)模型生成的预报降雨产品作为地表内涝模型的输入数据[5]。例如，Li等[5]集成了柳溪河水文模型与WRF气象模型，将其用于南方流域大尺度洪涝预报。该方法延长了内涝预报的预见期，但只能提供网格分辨率为20km×20km的预报产品。Tian等[6]开发了耦合WRF模型、河北水文模型及3DVar数据同化模块，构建了大气-水文预报系统。结果表明，该预报系统对空间分布均匀的降雨有较好的预报效果，但对短历时强降雨预报误差较大。Habibi等[7]采用了由X波段气象雷达数据驱动的NWS水文分布模型。结果表明，该系统能较好地预测水位和河流的变化，但其时间误差较为严重。

常规水文模型无法精确模拟内涝变化过程，而求解浅水方程(SWEs)的全水动力模型数值则可解决这一问题[8]。基于动力波方法的水动力模型能够模拟水力要素在复杂地形上的变化过程，生成准确的内涝变化计算结果。本文提出集成数值气象模型和水动力模型构建内涝预报模型，并以此模拟预报沣西新城内涝淹没过程，研究结果将为提升防灾减灾能力提供参考借鉴。

1 研究区域

研究区位于西咸新区沣西新城，总面积22.5 km2。全年降雨量主要集中在7～9月。排水管网局部淤堵及结构设计不合理，限制了该地区的排水能力，当遇到暴雨时，内涝频发。该地区布设有微型气象站，可以气象站监测数据作为地面实况，对GRAPES_MESO预报数据进行评价。

本文采用无人机航测技术对研究区地形数据进行了测量。数字高程模型(DEM)和数字正射影像图(DOM)的网格分辨率为2m，约560万个网格单元，如图1、图2所示。根据正射影像图，采用最大似然分类法将网格单元划分为五种土地利用类型，如表1所示。曼宁系数参考相关标准和城市排水文献确定[9-10]。

图1 DEM数据Fig.1 DEM data

图2 土地利用数据Fig.2 Land use data

表1 土地利用参数Tab.1 Properties of the different land uses

2 大气-水动力集成模拟预报模型

本文集成GRAPES_MESO大气模型、二维水动力模型和预报降雨重构方法建立大气-水动力内涝预报模型，图3为模型工作流程图。集成模型每隔12 h将滚动更新的预报降雨数据输入基于回归分析方法的重构模块进行改进，而后将重构的预报降雨数据作为水动力模型的输入数据。为了评估重构法对预报降雨和内涝数据的改善效果，本文还利用气象站的实测降雨数据和原始预报数据驱动内涝模型，进行内涝淹没过程模拟计算，并以淹没面积和水深预报结果分析模型性能。

图3 城市内涝预报工作流程Fig.3 Flowchart for the urban inundation forecasting

2.1 GRAPES_MESO气象模型

全球/区域同化预报系统(GRAPES)是中国气象局自主开发的采用多级数据同化的数值天气预报模型。根据不同分辨率要求，选择区域或全球中期物理过程软件包，可形成全球中期天气预报系统(GRAPES_GFS)或区域中尺度天气预报系统(GRAPES_MESO)。城市内涝预报模型采用GRAPES_MESO 4.0版本预报降水数据。预报系统每天00:00和12:00对中国及周边地区的未来72h降雨数据进行滚动更新，时空分辨率为3h、10km。

2.2 预报降雨数据重构方法

该预报降雨数据重构方法以GRAPES_MESO预报数据(数据来自国家气象科学数据中心http://data.cma.cn/)为自变量，气象站实测数据为因变量。利用SPSS统计分析软件的回归分析函数拟合修正公式。由于GRAPES_MESO的数值预报产品于2015年12月29日00:00首次发布，因此该方法选取的分析样本时间为2016年1月1日00:00至2019年1月1日24:00。在SPSS回归分析中，在筛去降雨量小于10 mm的场次降雨后，选取剩余降雨数据中的37个典型致涝降雨数据作为样本数据。为增强公式的适用性，各样本数据的降雨历时、雨强、累计降雨量都有所不同。此外，选择4个不同特征的降雨作为测试数据。

研究区第一和第三季度预报雨量小于实测数据，而第二季度预报雨量较大。由于第四季度的暴雨次数较少，因此不进行分析。本文选取线性方程、二次方程、三次方程、对数方程、指数方程、幂方程6种常用模型拟合公式。去除过度拟合和低相关性的函数，得到预测数据与实测数据之间的函数关系，如表2所示。第一季度符合对数函数，相关系数R2为0.81；第二和第三季度符合三次函数，R2分别为0.91和0.93。其中，xR为原始预测数据；xH为预测相对湿度；y表示重构值。

表2 修正公式Tab.2 Correction formulas

由于预报降雨数据的时间分辨率为3 h，当降雨持续时间小于3 h时，雨峰将变得不明显，致使降雨过程变得平缓，进而造成较大的误差。为了使重构后的降雨数据更符合实际情况，对采用修正公式(见表2)重构的降雨再进行重新分布。基于多年降雨数据的芝加哥雨型公式可以描述该地区的短期降雨过程，如式(1)所示。降雨将按芝加哥雨型每小时的雨量分配比例进行分配。根据公式计算结果，在第1小时分配降雨总量的10%，第2小时分配81%，第3小时分配9%。采用这种分配方法，可以在预报降雨时间小于3 h的情况下构建雨峰。

(1)

式中：q为暴雨强度；p为重现期；m为暴雨历时。

纳什效率系数(NSE)用于评估预测数据与实际情况的相似程度，ηNSE值越接近于1，则模拟结果越可靠。

(2)

2.3 二维水动力模型

2.3.1控制方程

模型控制方程为耦合水文过程的二维浅水方程(SWEs)，忽略运动黏性项、科氏力、风应力及紊流黏性项，其对应的二维非线性浅水方程守恒格式的矢量形式为：

(3)

其中，

式中：q为变量矢量，包括x和y方向的单宽流量qx和qy；u、v表示x和y方向的流速；t为时间；F和G分别为x、y方向的通量矢量；i为净雨率；S为源项矢量；zb为河床底面高程；Cf=gn2/h1/3为河床糙率系数，n为曼宁系数，g为重力加速度，h为水深。

2.3.2数值方法与GPU高性能计算技术

模型采用基于结构化网格的Godunov格式有限体积法求解SWEs。通过二阶显式Runge-Kutta法构造具有二阶时空精度的MUSCL型格式，从而确保物质守恒并有效解决不连续问题。针对模拟过程中可能产生急变流及非连续等复杂问题，模型选用HLLC近似黎曼求解器对单元界面上的质量以及动量的通量进行求解。底坡源项通过底坡通量法计算[11]。摩阻源项通过改进的显式方法计算[12]。模型代码使用C++和CUDA进行编程，从而能够在图形处理单元(GPU)实现并行计算，大大加快了计算速度。使用矩阵的方式将等待计算的单元分配给每个线程。每个GPU线程按顺序计算网格边界通量和源项。

3 结果及讨论

3.1 模型验证

以2016年8月25日00:36至8月25日12:27的实测降雨资料来验证模型的准确性。该场次降雨为双峰型，总历时为12h，降雨总量为97.2 mm，如图4所示。根据沣西新城管委会提供的监测资料，对其中4个内涝点(A、B、C、D)的模拟淹没面积进行分析，如图5所示。

图4 降雨过程线Fig.4 Rainfall process line

图5 内涝点分布图Fig.5 Inundation location map

表3为实测结果与模拟结果的对比。从表3中可以看出，模拟内涝面积与实测数据相吻合，相对误差最大为-4.2%，表明该模型精度较高，可用于大尺度复杂地形下的城市内涝过程模拟。

表3 实测结果与模拟结果比较(t=4 h)Tab.3 Comparison of measured results and simulated results (t=4 h)

3.2 预报降雨性能评估

以4场典型降雨(Ⅰ、Ⅱ、Ⅲ和Ⅳ)为例，对模型预报降雨性能进行评估。详细的降雨数据信息已在表4中列出。通过图6所示的原始预报降雨、重构预报降雨与气象站实测数据对比结果可以得出，重构数据的ηNSE值较原始预报数据的ηNSE值有了显著提高，说明该重构方法是合理的。但在原始预报误差较大的情况下，重构方法的修正能力也受到限制，存在着一定的不确定性，如Ⅰ、Ⅲ号降雨。

表4 典型降雨修正结果Tab.4 Corrected results of the typical rainfall events

图6 预报降雨过程线对比Fig.6 Comparison of forecast rainfall processes

3.3 淹没深度及面积预报性能评估

以2017年8月7日10:18发生的降雨致涝事件为例，该次模拟所用计算平台的CPU为IntelRCoreTMi7-8700，GPU为NVIDIA RTX2080。集成内涝预报模型耗时2.45小时完成了7小时内涝过程的模拟。由于GRAPES_MESO的预报降雨产品于2017年8月7日00:00发布，故此次内涝预报信息可以提前约9小时生成。内涝点模拟预报与现场实况的对比如图7所示，结果表明，预报模型可以准确预报出内涝位置。

以秦皇岛与开元路十字路口(图7中B点)为例，分析该区域淹没面积与最大水深的预报结果。图8、图9为采用三类降雨数据(实测、原始预报和重构预报数据)模拟预报的内涝过程。从图中可以看出，三类降雨数据模拟的面积和水深其变化趋势是一致的，但具体数值存在差异。由于GRAPES_MESO生成的原始预报数值高估了降雨量和持续时间，故用原始预报降雨量驱动预报系统得到的结果会稍大于用实测降雨模拟得到的预报结果，且峰现时间延迟了1小时；而利用重构数据模拟的预测结果更接近于实测结果，峰现时间也与实测降雨模拟结果一致，但1小时至2小时的模拟数值超过了使用原始预报降雨得出的数值，这说明重构方法仍存在一定的不确定性。

图7 2017年8月7日现场图片与模拟结果对比(t=6 h)Fig.7 Comparison of field pictures and simulation results on August 7, 2017 (t=6 h)

图8 不同输入降雨条件下的淹没面积变化过程Fig.8 Simulated inundation area by using different input data

图9 不同输入降雨条件下的最大水深变化过程Fig.9 Simulated max water depth by using different input data

图8、图9所示结果表明，使用重构降雨得出的内涝预报结果的相对误差在大多数时刻都较低。其预报面积和水深的平均相对误差分别为11.24%和3.75%；预报面积和水深的ηNSE分别为0.89和0.94。重构方法能够提升模型的预报性能，但重构方法改进的降雨数据存在一定的不确定性，可能会使预报结果出现偏差。同时，原始预测数据的时空分辨率不足也是影响预报准确性的另一因素。由于地形平坦，研究区域内的气候变化不大，因此，具有中尺度时空分辨率的预测数据可以在该区域使用。但是致涝强降雨具有局部分布的特点，应尽可能使用具有精细时空分辨率的预报降雨数据。

4 结 论

本文通过集成GRAPES_MESO数值天气预报模型与水动力模型，构建了城市内涝预报模型，并将该模型应用于沣西新城区域。

1) 针对GRAPES_MESO数值天气预报模型的预报产品存在不确定性的问题，基于回归分析方法，提出了一种重构预报降雨数据的方法。

2) 集成模型能够准确预报内涝点位置，预报淹没面积和水深的平均相对误差分别为11.24%和3.75%，ηNSE分别为0.89和0.94。

3) 在Ⅱ号典型降雨算例中，集成内涝预报模型耗时2.45小时完成了7小时内涝过程的模拟，内涝预报信息可于内涝发生前约9小时发布。

本文提出的城市内涝预报模型预见期较长、准确性较高，在城市防涝工作中具有巨大的应用潜力。内涝预报的准确性很大程度上取决于预报降雨的质量，因此，改进和优化数值天气预报模型仍然是提高内涝预报准确性最主要和最直接的手段。在下一步研究中，将开发天气预报模型降尺度方法，以便获取更高分辨率的降雨数据。