基于元胞自动机的城市雨洪模型研究

李继政，李传奇，王复生，张焱炜

(山东大学土建与水利学院，济南 250014)

0 引 言

洪涝是一种复杂的水文现象，并且具有频发、危害大的特点。随着我国城市化进程的不断推进，城市不透水面积增大，原有水文循环遭到破坏，加剧了城市洪涝灾害的问题，逢雨必涝已成常态。而洪涝模拟技术是研究城市洪水、防洪减灾的重要手段[1]，可以在洪灾发生前对其进行预测，并对其造成的影响进行评价。因此，如何构建精准的模型对城市雨洪进行快速响应具有重要意义。

目前，城市雨洪模型可分为3类：水文模型、水动力模型、简化模型[2]。不同模型有其各自优缺点，SWMM模型是代表性水文模型，但其对数据要求较高，应用存在局限性；水动力模型运用微分方程计算水流运动，因此精确度较高，但运算速率低，难以满足快速响应的要求[3]。而其他一些简化模型如元胞自动机(Cellular Automata, CA)模型由于其较低的数据要求以及快速运算速率已成为水文模型研究领域的热点。

元胞自动机起源于20世纪50年代初，在70年代迅速发展，目前已成功应用于各种自然现象和物理系统的模拟，如城市扩张[4]、土地利用变化[5,6]、交通流仿真[7]等。近些年诸多学者将其应用到洪水模拟当中，逯琳[8]详细阐述了利用元胞自动机来模拟暴雨洪涝过程，并进行了模型不确定性分析和参数优化。赖泽辉等[9]进行了城市地表积水的元胞自动机模拟，证明了CA模型能够很好的应用于城市内涝模拟。Prasetya[10]等利用元胞自动机进行了Solo River流域的洪水预测分析。Gibson[11]等基于元胞自动机对城市洪水模拟的精度与计算效率作了详细研究，结果表明CA模型与行业标准软件相比，处理速度有很大提升。

目前，元胞自动机在水文领域应用主要为自然流域，城市雨洪模拟相对较少[9]。在同时考虑坡面汇流及河道汇流情况时，元胞自动机网格化的汇流方法会受到一定限制，模拟结果可能存在较大误差[12]。对于城市区域，由于地面覆盖物的多样化以及排水设施的布置，CA模型将元胞集特征化的方法具有很好的适用性，其关键问题在于如何全面考虑各类土地利用类型并精确制定模型的产汇流规则。本文以济南主城区为研究区域，针对CA模型的演化规则，对各个产流因素分别确定较为精确的模型，采用曼宁方程结合最小差异法制定汇流规则，以求在CA模型的高效运行速率下也能够精准地模拟出城市各地的淹没水深，对城市洪水预测及防洪预警进行快速响应。

1 模型构建

元胞自动机是一种时间、空间都离散的局部网格动力学模型。它将系统分割成一个个小的元胞，每个元胞的状态受其相邻元胞状态的影响，通过制定合理的演化规则来模拟复杂系统的变化。

1.1 CA模型组成

元胞自动机由元胞空间、邻域关系、元胞状态、演化规则4部分组成。元胞空间即为整个研究区域网格点的集合，可为一维、二维或多维，元胞形状一般为三角形、正方形或六边形，不同的格网有其各自的优缺点[8]。本研究采用二维正方形元胞，方便与DEM栅格相对应。元胞的邻域类型通常有：冯诺依曼型、Moore型、扩展Moore型，不同的邻域类型决定了元胞之间的传递关系及影响范围，本研究采用Moore型邻域关系，即每个元胞的邻域由其周围8个元胞组成。演化规则是构建CA模型的核心，建立合理的演化规则是CA模型取得成效的关键，在CA模型中，这些规则决定了元胞状态的变化，城市洪涝受到社会因素和人类干预影响很大，外力或外部因素如何反映在演化规则中，怎么决定地表径流方向，是基于DEM的分布式降雨径流模型中的一个十分重要和关键的问题。CA模型基本结构如图1所示。

图1 CA模型结构图Fig.1 structure diagram of CA model

城市雨洪CA模型演化规则主要包括产流规则和汇流规则。对于产流规则，需要尽可能全面地考虑产流因素，如降雨量、植被截流量、土壤下渗量、雨水口排放量等，以及城市地表覆盖特征，如道路、绿地、房屋、高架桥等。汇流规则主要为确定水流方向和流量，需要全面考虑到特殊情况下的解决方案，如时间步长小于流动时间的情况，计算的分配水量高于中心元胞水量的情况等。

1.2 产流规则

降雨量扣除降雨过程的损失量即为形成径流的水量，降雨损失包括土壤下渗、植物截留、洼蓄、蒸发、雨水口排放等。流域产流计算方法通常有：蓄满产流法、径流系数法、下渗曲线法等[13]，这些方法通常侧重于某些产流因素(洼蓄、下渗等)的影响，适用于地面特征变化不大的流域区域。城市地表下垫面复杂多样，使其产流更加复杂，如何全面考虑不同下垫面的产流因素是构建城市雨洪模型的关键问题。对于城市暴雨来说历时较短，因此蒸发量一般可忽略不计，本研究采用式(1)来计算城市雨洪的产流量：

R=P-I-N-D

(1)

式中：R为产流量，mm；P为降雨量，mm；I为土壤下渗量，mm；N为植物截留量，mm；D为雨水口排放量，mm。

本文将城市土地利用类型划分为5种：城市绿地、林地、建筑物、水域、道路。不同的土地利用类型考虑不同的产流因素，具体规则为：①透水地表降雨损失为土壤下渗、植物截留，包括绿地、林地；②水域不产流；③含有雨水口的道路元胞考虑雨水口排放；④建筑物旁通常会有一定面积的绿地建设，而对于较低分辨率的土地利用类型无法考虑该部分的产流，因此本文将建筑物设为产流下垫面，与绿地有相同的产流因素，通过调整参数来体现建筑物的产流情况。采用上述产流规则基本能够兼顾城市区域不同下垫面下的重要产流因素，结合元胞自动机将研究区域网格化的特点，能够很好地构建城市雨洪模型。

1.2.1 土壤下渗模型

目前计算下渗量的方法主要有Horton方程、Philip方程及Green-Ampt方程[14]。大部分城市雨洪模型采用Horton方程计算土壤下渗量[15]，Horton公式反映了土壤下渗率随时间的变化关系，其表达式如式(2)所示。

It=fc+(f0-fc) e-ut

(2)

式中：It为t时刻的下渗率，mm/h；f0为初始下渗率，mm/h；fc为稳定下渗率，mm/h；u为衰减系数，1/h。

在应用到元胞自动机模型时，对每个时间步长进行积分，即可求得该时间步长内的土壤下渗量。

1.2.2 植物截留模型

影响林冠截留量的因素很多，如降雨量、降雨强度、林冠结构、温度等。为了分析植物截留特性，国内外学者建立了诸多不同的植物截留模型，这些模型主要基于Horton提出的截留机制[16]，即将截留量分为“枝叶吸附量”和因蒸发而增加的“附加截留量”。我国学者崔启武[17]和王彦辉[18]提出的植物截留模型已运用到许多研究中，可作为通用模型使用[16]。本研究采用崔启武模型来计算植物截留量，表达式如式(3)所示。

(3)

式中：N为截留量，mm；l0为林冠最大截留量，mm；C为林冠的特性函数，与本次降水无关，可用郁闭度表示；P为降雨量，mm。

该截留过程可以理解为：当降雨量较小时，截留率为常数，随着降雨量增大，截留量最终趋近于林冠的最大截留量。在应用到元胞自动机时，可用式(3)计算出每一时间步长的累计截留量，当前步长累计截留量减上一步长累计截留量即为本时间步长的截留量。

1.2.3 雨水口排放

除了下渗和植物截留外，城市暴雨洪水还需考虑雨水口的排放。通常雨水口沿着城市道路的边缘设置，雨水经过雨水口排放到城市的排水管网中。本研究采用式(4)[19]来计算雨水口的排放量：

(4)

式中：Q为孔口流量，m3/s；w为雨水篦孔口面积，m2；g为重力加速度，m/s2；c为孔口流量系数；h为雨水篦上水深，m；k为孔口阻塞系数。

本模型中雨水口与元胞相对应，含有雨水口的元胞需计算雨水口排放量，并且假设排水管网系统完全发挥作用，即不考虑排水管网的汇流情况。

1.3 汇流规则

目前，利用元胞自动机模拟洪水流动的汇流规则通常有两种：第一种为选取与中心元胞水位差最大的邻域元胞作为下游元胞[3,9]，即中心元胞只向一个元胞分配水量；第二种为基于动态系统向平衡条件发展原理的最小差异法[8]，中心元胞可流向低于其水位值的多个元胞来达到元胞之间的水位平衡。最小差异法已在熔岩流模拟、土壤侵蚀、洪水流动等模型中得到应用[19]，并显示出良好的效果。本文采用最小差异法结合曼宁公式来计算水流方向及元胞之间交换水量。最小差异法步骤如下：

(1)计算中心元胞及邻域8个元胞水位的平均值：

(5)

式中：H0为中心元胞水位，m；Hi为其邻域未去除元胞的水位，m;m为未去除的邻域元胞数量。

(2)将水位大于平均值的邻域元胞去除，即Hi>ave。

(3)计算剩余邻域元胞与中心元胞的水位平均值，并继续剔除水位大于平均值的邻域元胞。

(4)重复上述步骤直到没有邻域元胞被剔除。

(5)分配来自中心元胞的水量，使中心元胞与剩余邻域元胞具有相同的水位。

考虑洪水流动时间，采用曼宁公式计算元胞水流流速：

(6)

式中：V为流速，m/s；h为中心元胞水深，m；s为坡度；n为糙率。

洪水由中心元胞流向邻域元胞的时间为：

(7)

在设置时间步长t时，需大致计算出每个元胞的洪水流动时间，并使时间步长小于大部分的洪水流动时间，以保持模拟的稳定性。若t

2 实例应用

2.1 数据来源与处理

2007年7月18日，受高空槽、低涡等天气系统的影响，济南市出现了特大暴雨天气，最大点降雨量达到了180 mm。由于暴雨强度大，水流迅速汇集，导致济南市城区出现大面积积水，部分河道漫溢，低洼处房屋进水，给全市造成了巨大的经济损失。本文以济南主城区为研究区域(图2所示)，以“7.18”暴雨洪水为例，采用上述元胞自动机模型进行洪水模拟。

本研究所用降雨数据及实测水位数据来自2010年济南市水文年鉴，降雨数据采用3 h时段最大雨量，实测水位数据为各道路的洪痕水深，能够体现该次暴雨各地的最高水深。地形数据采用研究区域10 m分辨率的数字高程数据。土地利用类型数据来自地理空间数据云平台下载的Landsat8卫星遥感数据影响，经过ENVI和ArcGis工具的处理得到了研究区域30 m分辨率的土地利用类型数据，如图3所示，该数据包括了5种土地利用类型，即城市绿地、林地、建筑物、水域、道路。

图2 济南主城区影像Fig.2 Image map of Jinan main urban area

2.2 模拟结果与分析

模型参数选取经验值，依据来自SWMM用户手册及相关文献，模拟时间为3 h，时间步长设定为10 s。利用MATLAB工具对模型进行编程运算，最终得到整个研究区域不同时刻的淹没水深图，结果如图4所示。

选取5个站点的实测洪痕水位与模型模拟结果进行比较，结果如表1所示。

表1 实测水位与模拟水位对比Tab.1 Comparison of measured water level and simulated water level

图4 1、2、3 h时刻洪水淹没图Fig.4 Flood inundation map at 1、2、3 h

通过模型模拟结果可以发现，多数路段水深均超过0.5 m，部分路口和低洼地积水深度超过1 m，低洼处积水深度甚至可达3 m以上，与实际情况基本吻合。5个站点的模拟结果与实测数据相差不大，均方根误差为8.1 cm，相对误差基本保持在±20%以内，对于大面积高分辨率区域洪水模拟来说效果较为良好，后期可通过各种参数优化方法优化模型参数，进一步提高模拟精度。整个城区的洪水演化过程大致为：在降雨开始阶段，区域的淹没水深分布较为均匀，随着降雨及汇流的不断进行，水流总体向低洼处汇集，若降雨量较大，则会在低洼处形成积水，然后随着下渗和排水设施的作用淹没深度慢慢减小。

3 结 论

(1)与传统水文模型相比，元胞自动机没有复杂的数学物理方程，而是通过将研究区域划分成许多元胞，通过元胞之间的水量转化来模拟洪水演进，因此它对数据的要求较低，主要依赖高精度的DEM数据以及土地利用类型数据，在研究区域数据匮乏时该模型是一种合理且可靠的方法。并且该模型建模思路简单，结果表现形式直观，离散化时间和空间的特点使其可以模拟洪水的时空演进过程。

(2)5个站点的实测数据与模拟结果相对误差在±20%以内，对于城市这种具有复杂下垫面特征的区域来说，元胞自动机模型的模拟精度是可接受的。并且由于元胞自动机将研究区域网格化的特点，各个元胞的产汇流规则设置也较为灵活，针对城市各种下垫面可以进一步细化土地利用类型，对不同的土地利用类型设置多样化的演化规则，以求更加贴近实际情况。

(3)模型产流规则中参数较多，本文选取经验参数进行计算，未对其进行敏感性分析，一些参数优化方法如人工神经网络等可以适用于多参数的优化计算，这对于提高元胞自动机模型精度来说也是一种不错的方法。