改进的暴雨洪水淹没模拟算法及其在淮河流域的应用研究

卢燕宇 邓汗青, 田红, 何冬燕, 戴娟,

（1安徽省气象科学研究所，合肥 230031；2 安徽省大气科学与卫星遥感重点实验室，合肥 230031；3 安徽省气候中心，合肥 230031）

0 引言

由局地强降水造成的中小河流突发性洪水频繁发生，易造成人员伤亡。由于大部分中小河流站网偏稀，缺少必要的应急监测手段，预报方案和能力有待健全与完善，加上中小河流源短流急，洪水具有暴雨强度大、历时短、难预报、难预防的特点，针对暴雨诱发中小河流洪水和山洪地质灾害的预报和防御成为目前防洪减灾工作中突出的难点。充分利用暴雨后的短暂时间及时快速地做出洪水或洪涝预警决策，做好人员转移，是防范洪涝灾害保证人民生命安全的有效途径。因此基于降水监测结果，开展暴雨洪水演进模拟可以提前识别可能遭受灾害影响的承灾体信息，有助于做出合理预判和决策，能够更加有效地发挥预警实效。基于此，本研究将在现有暴雨洪水淹没算法的基础上进行改进，使之与现实更为贴近，将有助于发挥预警信息的最大效益，推动暴雨灾害风险实时动态评估的业务化进程，促进气象灾害风险管理理念和方法的发展与创新。

IPCC《管理极端事件和灾害风险，推进气候变化适应》特别报告（SREX）中指出构成灾害风险有3个要素，即致灾危险性、暴露度和脆弱性，只有将三要素有机结合才能合理表征灾害风险。就暴雨灾害风险而言，致灾危险性可由暴雨致灾临界阈值来反映，脆弱性体现承灾体对暴雨灾害强度的响应，暴露度则是连接危险性和脆弱性的关键环节，只有将空中降水转化为地面径流、洪水，才能判断出可能受灾害影响的承灾体，因而需要借助暴雨洪水淹没模拟来识别暴雨灾害的风险暴露度[1,2]。暴雨灾害的危害与洪水的淹没范围和水深直接相关，20世纪90年代以来，基于GIS空间分析技术，进行给定淹没高程下的洪水淹没分析已有不少研究[3-6]。但是目前大多研究仍然以静态分析为主，尚不能快速提供不同的入流量、不同时相的淹没状况和水深分布。由于水动力模型能够实现洪水演进的动态模拟，可以比较准确反映淹没范围、淹没深度及其历时特征，因而成为当前研究的一个热点方向[7-11]。而水动力学模型研究的重要问题是如何在充分考虑洪水演进的物理机制的同时高效快速地实现洪水演进的动态模拟，以及如何将洪水模拟结果与社会经济等数据相匹配结合，以评估灾害影响。当前，将水动力学模型与GIS技术相结合[1,7]，为这一问题的解决提供了思路，但仍然处于探索阶段，并且目前关于洪水淹没的模拟多直接以河流水文特征量为起始值，忽略了与灾害发生的根源——降水之间联系的分析，在这些方面仍有待进一步的研究。

针对风险预警评估业务服务的迫切需求，以及现有暴雨洪水淹没算法存在的问题，本研究将融合降水产流模型和二维水动力模型，改进淹没模拟算法，使之反映前期降水和下垫面要素空间变化对暴雨洪水形成及其演进的影响；通过优化降水空间插值，开展滚动化叠加运算，改进算法的输入和运行模式，使之体现降水时空变化对暴雨灾害风险分布的影响。选择淮河流域典型地区完成算法的检验和应用分析，从而为实时气象防灾减灾提供技术支撑。

1 资料与方法

1.1 研究思路与方案

开展暴雨洪水淹没模拟的关键是要抓住降水→径流→淹没三个环节，在这三个环节中需解决以下关键问题：1）降水精细化时空分布，2）降水向地表径流的转化，3）洪水的演进动态。在以往的暴雨洪水研究中这三方面存在一定的脱节，在开展洪水淹没模拟时往往仅考虑地形因素对暴雨洪水演进的影响，对降水产流过程以及降水时空分布变化的影响考虑不足，因而影响了暴雨洪水淹没模拟的仿真效果。因此本研究将着眼于现状问题与迫切需求，通过产流系数将降雨产流模型（SCS）与二维水动力模型（Floodarea）相融合（图1），使暴雨洪水淹没算法能够体现前期降水和下垫面特性的影响，提升仿真模拟效果；改进降水格网化方案和模型运算方式，实现滚动叠加运算，更好地体现降水时空分布变化对暴雨灾害风险分布的影响；最后基于改进的模型算法，在淮河流域典型区内对模型算法和参数进行检验和优化，实现暴雨洪水淹没的分布式动态模拟。

图1 技术路线和算法流程Fig. 1 Technical route and algorithm flow

1.2 降水产流模型

构建降水量与地表径流之间关系的关键问题是建立方便有效的降水产流模型。在对土壤入渗产流机理研究的基础上，许多学者提出了不同的估算坡面产流过程模型，其中美国农业部根据美国气候特征及多年水文径流资料所研发的SCS-CN模型由于其结构简单、所需参数较少、模拟结果准确度较高，而被广泛应用于场次降雨地表产流及其过程的预测中。国内外一系列分布式—半分布式生态水文模型，如SWAT，EPIC，CREAMS等，均采用 SCS-CN模型预测径流量。

SCS-CN模型包含一个水分平衡方程和两个基本假设，表达如下：

水分平衡方程：

比例相等假设：

降雨初损量与潜在滞蓄量正比关系假设：

式中，P为降雨量（mm）；Q为径流深（mm）；F为累计入渗量（mm）；S为流域最大蓄水能力（mm）；Ia为初损量（mm）；λ为初损率，无量纲。

联立方程（1）～（3）可以得到计算径流深的表达式：

S可由CN系数计算得到，二者有以下统计关系：

该模型的主要参数径流曲线数CN是反映流域前期土壤湿润程度（Antecedent moisture condition，AMC）、坡度、土壤类型和土地利用现状等综合特性的参数。SCS模型对CN的敏感性很高，有研究指出，CN取值±10%的变化可导致计算径流量-45%～55%的变化。可见，CN值的确定对降雨径流量的准确预测非常重要。除CN系数外，初损率λ以往通常取标准值0.2，但近年来研究发现初损率随不同区域取值变异性也较大，其值变化对模拟效果有较为明显的影响。

1.3 二维水动力模型

洪涝灾害的危害与洪水的淹没范围和水深直接相关，因此，确定灾害范围和程度可通过模拟洪水演进及其水文特征来实现。

本研究引入了德国Geomer公司研制的基于GIS的水动力模型——FloodArea模型。其原理是充分利用GIS栅格数据在水文-水动力学建模上的优势，实现GIS与水文-水动力学模型的数据融合。模型以栅格为基本单元，淹没模拟基于二维非恒定流水动力学模型，用Manning-Strickler公式计算每个栅格与周围栅格之间的水量交换。用Manning-Stricker公式计算每个栅格单元与周围8个单元之间的洪水流量，FloodArea模型计算洪水汇流示意图见图2，相邻单元的水流宽度被认为是相等的，位于对角线的单元，以不同的长度算法来计算；图中R为相邻单元的栅格距离，R为对角线单元的栅格距离，阴影部分指栅格面积，箭头指水流方向。水流的淹没深度为淹没水位高程和地面高程之间差值，由下式表示：

淹没过程中的水流方向由地形坡向所决定，地形坡向反映了斜坡所面对的方向，坡向指地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角，表征该点高程值改变量的最大变化方向，计算公式如下：

其中，α为地形坡度。FloodArea模型有3种基本的淹没情景：漫顶、溃口以及暴雨，能够以给定水位，给定流量或给定面雨量3种方式进入模型，并可根据水文过程线进行实时调整，可视化表达流向、流速和淹没水深等水文要素的时空物理场，为洪水淹没风险动态模拟提供了有效工具。

图2 FloodArea汇流计算原理示意[12]Fig. 2 Schematic illustration of the confluence calculation in FloodArea model

1.4 研究区概况及资料

以淮河流域中的史河上游为例开展算法改进实验和应用研究。研究样区如图3所示，位于大别山腹地金寨县境内，以黄泥庄水文站作为研究区流域的出水口断面，获取了2008—2010年多场洪水过程的水文数据，用于计算实际径流深。研究河段属山区型河道，水浅流急，支流较多；流域面积约为800 km2，海拔高度悬殊大，研究区内高差超过1200 m，土壤类型以黄棕壤为主，植被覆盖类型主要为林地。研究区及周边布有多个气象观测站，可获取逐时雨量观测值用来计算流域面雨量。

图3 研究样区概况、地形、水文站和气象站分布Fig. 3 Overview of the study area, topography and distribution of hydrological and meteorological stations

下垫面信息包括数字地面高程（DEM）、土地利用类型等，DEM数据主要来自于1:5万地理信息数据，根据研究区范围进行裁剪得到，数据分辨率为25 m，土地利用类型来自于GLC30土地类型数据库，分辨率为30 m。气象数据来自于全省自动气象站的逐小时降水数据集。水文数据主要来自于安徽省水文遥测信息网。

2 结果与讨论

2.1 降水空间插值方案的优选

为使气象数据与模型输入相匹配，获取精细化的降水时空分布特征，采用空间插值的方法将离散的站点数据插值成细网格化的面状数据[13]。插值方法将考虑反距离权重（IDW）、自然邻域、克里金（Kriging）和薄盘样条 等方法，同时根据降水空间分布特点，采用海拔高度、地形起伏度、坡度、坡向等地形因子作为协同变量，构建空间插值模型。采用交叉验证（Cross-validation）的方法来比对不同插值方法和模型的实际效果，从中选择最优化处理方案，作为台站降水数据的空间拓展和网格化推算的方法。

选择近年来多次降水过程进行插值方法对比分析，其中包括了热带气旋、梅雨锋等不同天气系统所导致的强降水过程，具有一定的代表性。采用交叉验证的平均相对误差、均方根误差来衡量模拟效果，由于研究重点关注暴雨，因而进一步引入纳什效率系数来考察对于极值的拟合效果。从图4可以看出就平均相对误差来看，IDW和RBF 法拟合值总体比观测值偏小，其他3种方法则略偏高，其中以Co-kriging方法误差最小；对于均方根误差来看，RBF法的误差最大，而Kriging法的估计误差最小；纳什效率系数越接近于1说明拟合结果与观测值越接近，可以看出同样RBF法的估计偏差最大，而Kriging法则与观测结果最为接近。综合3种指标可以看出以Kriging法的拟合效果最佳。

图4 不同降水插值方法交叉验证的平均相对误差（a）、均方根误差（b）和纳什效率系数（c）Fig. 4 Mean relative error (a), root mean square error(b) and Nash efficiency factor (c) for cross-validation of different precipitation interpolation methods

2.2 降水产流模型的参数分析与应用

2.2.1CN系数计算

SCS模型中CN系数反映了流域下垫面的产流能力，对后续降水产流及淹没模拟计算具有重要意义，CN值与土壤类型、土地利用类型以及前期土壤湿润程度密切相关。根据研究区土壤类型和土地利用类型数据，通过查表可以得到中等湿润程度下CN值分布（图5），再通过计算研究区前5 d累计面雨量即可获取实时的CN系数。

由于SCS模型构建和参数化过程中主要依托于美国大量缓坡（坡度5%左右）降雨产流资料发展起来，而研究区地形复杂，海拔起伏大，在实际应用时还应进一步考虑坡度对CN系数的影响。根据已有研究结果，采用Huang等[14]的坡度修正公式（见式（8））对CN系数进行修正，使之适用于起伏山地的降雨产流模拟。

图5 研究区土壤类型（a）、土地利用类型（b）和AMC II 下的CN值（c）分布Fig. 5 Distribution of soil type (a), land use type (b) and CN values (c) under AMCII in the study area

式中：CNII为中等湿润程度下的CN系数，slp为平均坡度（%），CNIIs为坡度修正后的CN系数。

基于DEM计算了研究区的坡度分布，并采用公式（8）对CN系数进行了修正，确定了AMCII条件下的CN值分布结果（图6），作为SCS模型应用的输入参数。

图6 研究区坡度（a）及其修正后的CN系数（b）分布Fig. 6 Distribution of slope (a) and its corrected CN coefficient (b) in the study area

2.2.2 初损率λ优化

当确定了AMC II条件下的CN系数后，只需要计算API判定实时AMC条件，根据SCS模型公式可知仅有一个未知参数即初损率λ，在以往的SCS模型应用中，一般取λ=0.2，但该值是根据美国土壤水文观测资料推算而来，是否适用于本地区还需进一步研究，这里选择2008年的多次典型降水过程分析λ的取值对模拟结果的影响，并确定优化后的参数值。

参数优化的依据为确定系数（R2）、相对误差（PBIAS）和效率系数（NSE）。使用确定系数R2评价模型追踪观测值变化的准确程度，越接近于1表明模型越能准确反映变化趋势；利用相对误差PBIAS评价模拟结果高于或低于观测值的平均趋势，该值大于0表示模型低估产流量，小于0表示模型高估产流量，越接近于0表明偏差越小；采用NSE表示模拟结果与观测值的契合程度，越接近于1表明模型越准确。

从表1可以看出模拟径流深基本上随着λ的增加而减少，这也可以用λ的物理意义来解释，即λ增大，降水损失量增加，转化为径流的比例相应减少。对比历次过程，模拟径流深基本上是当初损率λ为0.05～0.15范围内与实测结果较为接近。进一步通过R2、PBIAS和NSE来看，当λ=0.10时SCS模型模拟效果表现最佳（图7），因此选择0.10作为研究区的λ的优化取值。

表1 典型洪水过程实测与模拟径流深Table 1 Observed and simulated runoff depths for typical flood processes

2.2.3 SCS 模型应用及检验

在确定了SCS模型的CN系数和λ后，进一步建立不同条件下CN系数的概化转换公式（图9 ），以中等湿润条件B类土壤的CN系数作为模型输入参数，然后结合API的分布结果进行CN系数的实时计算，从而实现降水产流量的实时分布式模拟（图8）。

基于模型应用流程，采用率定好的模型参数集，以降雨量作为输入实时模拟了2010年多次降水过程的产流量（表2），模拟与实测的相关系数在0.9以上，相对误差约为-15%，模拟值略小于实测。总的来看，率定后的SCS模型能够较好地适应研究区，可以实时反映降水产流关系，从而为后续暴雨洪水淹没模拟提供可靠的输入前提。

2.3 二维水动力模型的参数分析与应用

在二维水动力模型FloodArea中影响模型效果的主要控制参数为Manning-Sticker公式中的地面糙率，根据已有研究表明，糙率系数与土地利用类型有密切关系，参考相关研究及标准规范[15]，按不同土地利用类型设定了糙率系数（表3）。该套系数已在大别山区和淮河王家坝等地区经过了验证[11,16]，能够较好地反映实测淹没水深的时空动态分布。

在以往暴雨洪水淹没模拟中，多以过程模拟和事后评估为主，降水输入模型时仅体现降水随时间变化，不考虑降水空间分布的动态变化，即在一次过程模拟中，暴雨中心位置常常假定不变。本研究选择史河上游2010年7月11日暴雨过程为例研究FloodArea模型在强降水致洪淹没分析中的应用，基于实况观测资料采用Kriging方法计算了该日逐小时面雨量，经统计24 h累计面雨量超过80 mm，流域出口断面黄泥庄水文站水位上涨接近1 m。基于1:5万地理信息数据和糙率信息，将降水空间分布和过程雨量带入到FloodArea模型进行洪水淹没动态分析，同时为了考察降水输入方式对模拟结果的影响，我们进一步改变了模型输入方式，采用滚动迭代输入将逐小时的降水空间分布代入模型，进行叠加运算。

图7 不同初损率模拟结果的确定系数（a）、相对误差（b）和效率系数（c）Fig. 7 Determination factor (a), relative error (b) and efficiency factor (c) for different initial loss rate simulations

从图9可以看出，在整个过程中降水中心的改变十分明显，即使是比较接近的时刻降水的空间分布仍有较大变化。而如果在模型输入时采用固定的降水空间分布将不可避免地带来较大的偏差，例如在本例中过程累计雨量主要集中在流域的南部，而实际上在过程中降水中心时常在南北摆动。

对比两种输入方式的模拟结果（图略），可以看出淹没水深的空间分布具有一定差异，主要表现为相比于默认输入方式，滚动模拟结果流域北部淹没范围更大，而南部相对淹没范围较小，水深较浅。进一步采用黄泥庄水文站的观测数据进行分析（图10），可以看出采用默认的固定降水输入模拟的水深变化略滞后于实测，而滚动模拟结果则与实测较为同步，这可能是由于过程累计降水中心距流域出口断面较远，采用固定的降水空间分布使得汇水时间增长，从而流域出口的水位变化与实际相比存在滞后。

总的来看，FloodArea模型对水深的动态变化模拟效果较好，能够较为准确地反映水位的涨落情况，但两种输入方式下水深均显著高于实际，这主要是由于FloodArea模型在模拟过程中假定降水全部转换为径流参与汇流和淹没分析，而实际情况则是在降水产流过程中，有部分降水以下渗等方式损耗，产流量要低于降水量，因此在模拟中还需要进一步考虑产流系数，将降水转化为产流降水才能与实际情况更为接近。

图8 SCS模型实时应用流程示意Fig. 8 Schematic representation of the SCS model real-time application flow

2.4 模型耦合及算法验证

2.4.1 耦合思路与算法改进

基于精细的格网化降水量，通过引入SCS降水产流模型，综合前期降水量、土壤和土地利用类型、地形因子等参数来计算历次降水的产流系数，利用产流系数来对实时降水进行修正，并计算每个格点上实际产生的地表径流量，之后将实际产流量代入到FloodArea中运算模拟暴雨洪水的动态演进情况。即以“产流系数”这个关键环节来将SCS降水产流模型与FloodArea洪水淹没模型相融合，形成分布式的暴雨洪水淹没算法，能够充分反映流域内前期降水和下垫面要素空间变化对暴雨洪水形成及其演进的影响。

此外，在以往暴雨洪水淹没模拟中，多以均匀雨量输入模型，不考虑降水时空变化的影响，本文进一步采用滚动迭代的方式进行逐小时模拟以体现降水空间分布随时间的变化，以实现更好的仿真效果。

图9 累计降水分布（a）和降水空间格局动态变化（b）Fig. 9 Cumulative precipitation distribution (a) and spatial pattern dynamics of precipitation (b)

图10 不同降水输入方式下的淹没模拟结果与观测值对比Fig. 10 Comparison of inundation simulation results with observations for different precipitation inputs

表2 典型暴雨过程的观测与模拟径流深Table 2 Observed and simulated runoff depths for typical storm events

表3 不同土地利用类型对应的糙率系数Table 3 Roughness coefficients corresponding to different land use types

2.4.2 改进后算法模拟效果分析

为进一步验证改进算法的模拟效果，仍然采用2010年7月11日史河上游的模拟案例，从图11可以看出，在考虑了降水产流和损耗后，模拟的黄泥庄水深变化与实际更为接近，动态趋势吻合度更好，并且系统偏差也明显缩小，改进后的算法对暴雨洪水淹没模拟误差比之前缩小约50%。

3 结论

针对现有暴雨洪水淹没算法存在的问题，本文在降水空间插值、降水产流模型参数优化、二维水动力模型本地化、模型耦合及应用等方面开展了系统分析，融合SCS模型和FloodArea模型改进淹没模拟算法，使之实现更好的仿真效果。主要结论包括：

1）交叉验证结果表明以Kriging方法插值效果最佳，因此选用Kriging法作为降水空间插值和面雨量计算的优选方法。

2）基于下垫面信息实现了SCS模型的参数化，采用地形坡度对CN系数进行了修正；并通过典型个例模拟及率定，对SCS模型中初损率值进行了优化，率定后的SCS模型能够较好地反映降水产流关系，从而为后续暴雨洪水淹没模拟提供可靠的输入前提。

3）根据土地利用类型确定了FloodArea模型的地表糙率系数，降水输入方式对洪水动态模拟结果有一定影响，改进后的滚动模拟方式可以实现更好的仿真效果。

4）以降水产流过程为纽带实现了SCS模型与FloodArea模型的耦合，形成分布式的暴雨洪水淹没算法，来反映流域内前期降水和下垫面要素空间变化对暴雨洪水形成及其演进的影响，改进后算法对模拟效果有显著提升作用，模拟误差缩小约50%。

总的来说，通过参数优化、模型耦合等方法对现有暴雨洪水淹没算法进行了改进，使之与现实更为贴近。由于本研究重点关注降水落到地面后的产汇流过程及其可能带来的洪涝灾害，因此主要以降水实况监测作为输入来研究算法模拟效果。未来为进一步延长预见期，还有必要将预报降水与模拟算法相结合，探讨误差来源，优化匹配策略，以更好地与风险预警服务业务相衔接，在发布灾害性天气预报后，如果达到或超过临界条件，进一步开展暴雨洪水淹没模拟分析，指出灾害高风险区及可能灾损，为政府和群众提供直观的暴雨灾害风险信息，有助于发挥预警信息的最大效益。

图11 改进后算法的模拟效果对比Fig. 11 Comparison of the simulation effects of the improved algorithm