基于栅格与云模型的风暴潮洪水风险模拟评估方法——以珠海市香洲区为例

孙 海，嵇文捷，郑雅芝

（1.中国海洋大学工程学院，山东青岛 266100；2.中国海洋大学海洋发展研究院，山东青岛 266100）

引言

风暴潮灾害在我国是最常见和破坏性最大的自然灾害之一。风暴潮洪水涌入沿海城市，会造成多方面的负面影响，如堤坝桥梁等基础设施被冲毁，城市道路被淹没导致交通中断，港口码头仓库中的库存物资因水位上升被浸泡损失严重，海岸防护工程受损巨大［1］等。同时，沿海地区人口和经济高度聚集［2］、地面沉降［3］等因素亦会加剧风暴潮灾害的威胁，最终造成大量的人员伤亡和财产损失。风暴潮已经成为沿海地区经济发展的一个严重制约因素［4，5］。

准确高效地进行风暴潮洪水模拟及风险评估可以最大限度地减轻风暴潮灾害损失［6］，模拟风暴潮洪水的演进对沿海城市灾害应对有着重要的指导意义。目前有许多成熟的二维水动力学模型可以用于洪水演进模拟，其中常用的有DHI的MIKE 21，SOBEX，JFLOW及UIM等［7］。上述模型可以得到较为精确的洪水随时间空间变化的信息，如在不同时刻、不同地点的水深、流速等。但这类模型对数据的要求较高，需要地形、糙率、水位边界条件等原始数据，且输入模型前需进行数据格式转换，转换过程会产生一定的误差和丢失精度。此外，传统的二维水动力学模型求解过程复杂，操作和计算耗时长，无法适应实时应用。元胞自动机（cellular automata，CA）是一种时间、空间、状态均离散，空间相互作用和时间因果关系为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力，在模拟地表径流和洪水事件等方面具有较强的优势，能很好地符合洪水波的运动变化规律，可便捷高效地与GIS系统集成。胡蓓蓓等［8］提出基于栅格图像的种子蔓延算法，该方法能够较好模拟有源淹没的洪水。刘坤等［9］提出基于栅格的区域生长算法，实现了对堰塞湖溃坝后淹没区的预测。上述元胞自动机模型可以直接利用DEM地形数据，通过输入起始位置和水位近似求解，计算规则简单，计算效率高。弊端是模拟过程不涉及水动力学计算，而是按照较为简单的转换规则更新元胞状态，无法得出精确的洪水随时空变化的信息，不能很好地反应水流实际运动规律。考虑到元胞自动机的计算效率优势和水动力学模型的准确性，一些学者尝试将二者结合起来。Guidolin［10］团队利用水位权重和曼宁公式的网格动力学方法模拟了考虑降水的河道洪水淹没，Dottori和Todini［11，12］建立了基于扩散方程的二维元胞自动机模型。上述研究表明元胞自动机洪水模拟的准确性得到了有效提高，初步实现了对洪水随时空变化过程的模拟。

近年来越来越多的案例表明，如何应对和处置洪水灾害造成的风险，也是抢险救灾、确定抢救重点、制定应急方案的重要部分。进行洪水风险评估需要综合考虑各种因素，包括灾害强度与区域特征。部分学者通过传统的洪水风险评估方法，例如层次分析法［13］、模糊聚类模型［14］、信息扩散理论［15］等，开展了洪水灾害危险程度的区划研究，研制出基于可靠性分析的模糊风险计算模型。上述研究，将洪水风险区划的模糊性、随机性，灾情数据的不完备性和研究区之间的关联性分别进行了研究。但是，模糊性、随机性、不完备性以及关联性都是洪水灾害风险评估中的一部分，需予以综合研究。云模型是解决这一问题的有效手段，该模型能兼顾等级概念的模糊性与随机性，能从实际数据分布中抽取等级概念，实现不同层次上的分析与综合，能更好模拟洪水灾害风险等级的概念。目前，基于云模型进行风暴潮风险的研究相对较少，本文拟根据灾害风险系统的定义以及洪水灾害的特点，以云模型模拟区域洪水灾害风险的不确定性，以期为区域洪水灾害风险评估提供新思路。

综上所述，本文提出了基于栅格水动力学的元胞自动机模型，将二维浅水动力学方程离散化并作为转化规则引入到元胞自动机模型，利用动量和能量守恒机制模拟风暴潮洪水演进过程，通过简化次要信息实现及时高效地模拟风暴潮洪水致灾过程。同时，在采用基于栅格水动力学的元胞自动机模型对研究区域进行风暴潮洪水仿真的基础上，引入基于云模型的多属性风险综合分析方法，综合考虑人口、经济、社会环境因素与评价指标的随机性和模糊性，结合模拟得出的淹没水深、淹没范围等风暴潮洪水灾害信息进行综合分析，弥补了灾害实时综合风险评估的不足。最后，对基于栅格水动力学的元胞自动机模型进行实时风险评估的有效性进行实例验证。文章通过基于栅格水动力学的元胞自动机综合风险度耦合计算模型，深化了对风暴潮灾害动力学与风险评估过程的研究，以期为风暴潮灾害的快速精确模拟及风险评估提供新的思路，形成一套高效的风暴潮灾害风险预警方法。

1 研究方法

1.1 浅水动力学方程离散化过程

元胞自动机由元胞空间、元胞状态、邻域及转化规则组成［16］。本文选用二维元胞空间，邻域采用Von Neuman型元胞网格，其洪水演进过程见图1。元胞自动机模型中每个元胞下一时刻的状态（水位、流量）由该元胞及其邻域元胞的上一时刻状态（水位、流量）决定［17］。元胞自动机最关键的部分就是定义转换规则。因此，为准确描述洪水演化过程，本文拟将二维浅水动力学方程离散化，并将其作为元胞自动机模型的转化规则引入，利用动量和能量守恒机制模拟风暴潮洪水演进过程，水位和流量的演化关系通过圣维南方程［18］进行控制，见公式（1）～（3）。

图1 Von Neuman型邻域Fig.1 Von Neumann type neighborhood

式中：h为水深，即水面离河底的距离；Z为水位，即自由水面相对于黄海基面或测站基面的高程；M、N分别表示x、y方向上的平均单宽流量；u、v分别表示在x、y方向上的平均流速，n为曼宁糙率，t为时间，g为重力加速度。动量方程的第1项反映局地加速度，第2、3项反映对流加速度，1、2、3均为惯性项；由于假设的风暴潮漫滩区域地形较为平坦的符号随涨落洪而变化，参考扩散波的简化方法［19］，可知公式（3）和式（2）同理。且由于第2、3项对流项为非线性项，在进行微分方程求解时会造成计算结果的振荡，对计算影响较小，因此本文将其省略［20］。对式（1）～（3）简化，进行显式差分，联立求解可得到和，见公式（4）～（6）。

其中，k为元胞模拟变化次数，Δt为时间步长，Δx、Δy定义了元胞大小，i，j定义了元胞的位置，代表第i行j列。初始化网格的水位值，初始水流速度（w为风暴潮的风速导致的初始流速）。通过反复利用式（4）～（6），就可得到每一时刻的元胞状态。元胞空间差分示意图，见图2。

图2 元胞空间差分示意图Fig.2 Cellular space difference diagram

显式差分需满足收敛性和稳定性条件，需对时间步长与距离步长进行限定，如公式（7）所示，c是波速，一般取

根据上述分析，基于元胞自动机的淹没模型算法过程如下。首先，将所研究的区域划分成若干大小一致的矩形单元（元胞），单元尺寸应在10～200 m之间。设置过小，单元水平比例尺接近垂直比例尺，无法用沿水深方向的平均值表示水力要素。设置太大，则会损失精度。为提高模拟效果，本文选择50 m作为单元大小。这些元胞组成了元胞空间，每个元胞上一时刻的状态决定了元胞的下一时刻的状态，同时这个变换过程受转化规则控制。元胞自动机模型网格图及淹没过程流程见图3、图4。

图3 元胞自动机模型网格图Fig.3 Grid chart of the CA model calculation

图4 元胞自动机模型计算流程图Fig.4 Flowchart of the CA model calculation

1.2 基于云模型的风暴潮灾害多属性综合风险分析方法

本文首先构建风暴潮灾害综合风险评价指标，常用的指标有：台风持续时间、环境影响、经济损失、人口伤亡、海洋工程损伤长度、社会影响、人均GDP、医疗机构数量等［21］。本文选取人口伤亡、经济损失、环境影响、社会影响四个权重较高的指标，能较好地反映灾害导致的经济、社会环境损失及受灾人口数。其次，建立基于云模型的风暴潮灾害的多属性综合风险分析方法，综合考虑指标的不确定性、模糊性和随机性，评价洪水灾害风险等级。

1.2.1 风暴潮灾害综合风险评价指标估算

风暴潮灾害导致的损失总体分为经济损失和非经济损失。经济损失分为直接经济损失和间接经济损失［22］。本文将受灾人口数及社会和环境损失归为非经济损失。

（1）经济损失估算

直接经济损失根据受灾区域的土地类型进行分类估算，土地利用类型主要包括耕地、林地、水域、城市建成区和开发区五类，其中，耕地、林地及水域三类土地利用类型主要考虑减产造成的经济损失，城市建成区主要考虑家庭财产损失、房屋和基础公共设施的破坏，开发区主要考虑房屋和基础公共设施的破坏，其计算公式为：

其中，L os s直为风暴潮灾害导致的直接经济损失，Cos t为各土地利用类型单位面积的重置成本，Area为各土地利用类型的受灾面积，Vulnerab ility为各土地利用类型的潮灾损失率。由于间接经济损失影响因素众多，难以精确计算，本文采用系数法对其进行估算，公式为：

其中，L os s间为风暴潮灾害间接经济损失，B为折算系数。对于不同行业折算系数B的推荐取值区间为：农业15⁃30%，工业16%⁃35%［23］，由于风暴潮灾害属于较为严重的自然灾害，本文研究区域为珠海市香洲区，对经济造成的损失主要体现在工业方面而非农业方面，故本文折算系数取工业部门推荐取值的平均水平，取B=0.25。同时，通过研究区域历史资料和其他区域的风暴潮灾害损失率，可确定研究区域的灾损率Vulnerab i lity，见表1。

表1 各土地利用类型单位面积价值及不同水深对应的灾损率Table 1 The unit area value of each land use type and the damage rate corresponding to different water depth

土地利用类型根据国家标准GB/T 21010-2017土地利用现状分类确定，其中城市建成区包括住宅用地、商服用地、公共管理与公共服务用地，开发区包括工矿仓储用地、交通运输用地、其他土地。

（2）受灾人口数量估算

受灾人口估算模型采用Jonkman提出的基于洪水强度的脆弱性曲线公式［24］：

其中，y为洪灾伤亡率，x为淹没水深。本文通过土地利用类型提取城市建成区和开发区用地，结合人口统计资料及不同区域人口密度分布规律得到区域人口空间分布，根据洪灾伤亡率估算伤亡人口。

（3）社会与环境受灾程度估算

社会与环境受灾程度是指洪水灾害对城市政治、文化、生态环境等带来的损失，由于影响因素复杂，难以计算或用货币衡量。本文根据文献［25］中定义的社会与环境风险及风险影响指标参考值，估算潮灾带来的社会环境影响，见式（11）：

式中，f为社会与环境影响指数。其余系数均分为5个级别，其中N为风险人口系数，取值区间1.0-1.2、1.2-1.6、1.6-2.4、2.4-4.0、4.0-5.0，对应的风险人口数为1-10、10-103、103-105、105-107、107人以上；C为重要城市系数，取值1.0、1.3、1.6、2.0-3.0、4.0-5.0，对应散户、乡村、乡镇、县级市至地级市、直辖市或省会至首都；I为公共设施系数，取值1.0、1.2、1.5、1.7、2.0，对应设施重要程度为一般、一般重要、市级重要、省级重要、国家级重要；h为文物古迹系数，取值1.0、1.2、1.5、2.0、2.5，对应文物古迹珍稀程度为一般、县级、省市级、国家级、世界级；R为河道形态系数，取值1.0、1.3、1.6、2.0-3.0、4.0-5.0，对应受破坏规模为河道受轻微破坏、一般河流受一定破坏、大江大河受一定破坏、一般河流至大江大河受严重破坏，一般河流至大江大河改道；l为生物生活环境系数，取值1.0、1.2、1.5、1.7、2.0，对应丧失栖息地动植物珍稀程度为一般、较有价值、较珍贵、稀有、世界级濒临灭绝；L为人文景观系数，取值1.0、1.2、1.5、1.7、2.0，对应受破坏景观级别为自然、市级、省级、国家级、世界级；P为污染工业系数，取值1.0、1.2、1.6、2.0-3.0、4.0，对应污染工业种类为基本无污染工业、一般化工厂或农药厂、较大规模化工厂或农药厂、大规模化工厂或农药厂至剧毒化工厂、核电站或核储库。

1.2.2 基于云模型的灾害风险多属性评价方法

参考《海洋灾害等级标准》（DB21T 1715-2009）及洪水淹没风险等级划分［26］，风暴潮灾害影响程度可划分为Ⅰ、Ⅱ、Ⅲ、Ⅳ四级，分别表示特别重大、重大、较重、一般，其划分标准，见表2。

表2 风暴潮灾害综合评价指标赋值参考Table 2 Comprehensive evaluation index reference for storm surge disaster

公式（12）中，Ex ij为云滴在论域空间分布的期望，在本文是指标在论域中的中心点；En ij为熵，表示指标的不确定性度量，由指标的随机性和模糊性共同决定；He ij为超熵，是熵的不确定性的度量。

2 实例验证

2.1 研究区域

珠海市地处广东省东南海岸，人口密集，经济发达，由于台风经常袭击，风暴潮给该区域带来重大的经济损失和社会影响［28］。本文以2017年第13号台风天鸽为例，选择受灾严重的珠海市香洲区（22.06°N-22.30°N、113.41°E-113.60°E）作为研究区域，见图5。

图5 研究区域Fig.5 Study area

2017年8月23日12时50分前后，台风天鸽在广东省珠海市金湾区登陆，珠海站观测最大风暴潮增水达2.79 m，为1965年以来登陆珠江口的最强台风。据2017年中国海洋灾害公报统计，广东省码头损毁1.22 km，防波堤损毁240.18 km，海堤换损毁532.08 km，道路损毁0.02 km。数字高程数据取自30 m×30 m分辨率的GDEMDEM数据集，遥感图下载于Google Earth（2018）。

2.2 风暴潮洪水模拟与验证

为验证本文建立的栅格水动力学模型的有效性，拟将模型模拟的结果与传统二维水动力学模型的仿真结果进行水深与计算时间对比，对提出方法的性能进行验证。除此之外，我们还依据从社交数据集提取的台风天鸽登陆后香洲区多处实际淹没点在不同时刻的水深，与模拟水深进行对比分析，辅助验证模型的实际表现。目前流行的传统二维水动力学模型有MIKE 21、SOBEX、JFLOW、UIM等，它们的计算方法、时空离散化的形式和运算特点等也有所不同［7］，见表3。

表3 洪水淹没模型对比Table 3 Comparison of flood inundation models

据表3可知，MIKE 21在模拟洪水方面有可并行运算以及加入惯性项运算的优势，具有较好的时空离散化的表达形式和运算快速、精确等优点。因此，本文将对比基于栅格水动力学的元胞自动机模型与MIKE 21模型的模拟结果，从水深差值、统计指标（均方根误差及相关系数）和计算时间3个方面比较模型。用于对比实验的MIKE 21模型采用非结构化网格，糙率设定为固定值0.03，取0.05 m作为干湿边界的判断条件，水位边界条件取自三灶站与赤湾站的实测数据［29］。同时，据中国天气台风网播报，该地区受台风影响时长约为8 h，故本文模拟时间定为8 h（10：30-18：30）。

2.2.1 水深对比

选取t=240 min时刻的水深数据分别绘制洪水淹没图和水深绝对差值图，见图6，根据选取时刻点的淹没水深图和差值图，基于栅格水动力学的元胞自动机模型与MIKE 21模型模拟结果有良好的一致性，极少部分区域存在一些小的差异，且水深差绝对值不超过0.1 m。

图6 t=240 min淹没水深对比Fig.6 Comparison of simulated flooding at 240 min

选取5个时刻的水深数据，通过计算均方根误差（RMSE）和相关系数（R）验证结果，见表4。其计算公式见式（16）和式（17）。与表示CA模拟在元胞C i,j处的水深和平均水深，与表示MIKE 21模拟在元胞C i,j处的水深和平均水深。

表4 MIKE 21和基于栅格水动力学的元胞自动机模型结果统计指标对比Table 4 RMSE and R value of the MIKE 21 and CA model based on grid hydrodynamics at five moments

由表4知，在5个时刻两模型水深数据的均方根误差都在0.04 m以下，说明CA模型模拟的淹没范围及水深的精确度较高，与MIKE 21的结果一致性较好。此外，它们的相关系数均不低于0.93，这证明两模型的模拟结果相关性较高，异常性较低。

2.2.2 计算时间对比

计算时间是本文模型需要考虑的重要因素和评价模型效率的重要标准。本次模型对比分析均在台式服务机上运行，处理器为因特尔I5-9 400F CPU@2.90GHz，内存为16GB，操作系统为Windows 10。MIKE 21的模拟时间为3.27 min，CA模型的模拟时间为1.23 min。并且MIKE 21前期需要进行模拟区域的网格处理和边界条件、参数的设置，增加了模型模拟的复杂性和工作量，而CA模型仅需设置几个初始参数，模拟时间减少62.4%，有利于实现风暴潮洪水淹没过程的简单快速评估。

2.2.3 水位辅助验证

国内外相关研究证明，社交媒体数据由于其自身的时间和地理空间属性，可以应用于灾害事件的实时监测和趋势预测。数据主要集中在暴雨暴发后的24 h内，本文从社交数据集获取150条相关信息，构建基于社交媒体的灾情信息网络验证集。其中，以保西路（22.16°N，113.48°E），洪湾港（22.16°N，113.46°E），东堤（22.14°N，113.54°E），琴海西路（22.13°N，113.46°E）4个特征点位置的照片及相关信息举例说明，见图7。

图7 基于灾情信息的网络验证集示例Fig.7 Example of network verification set based on disaster information

通过挖掘特征点的照片与相关灾情信息，当模拟时间t=140 min时，保西路实际水深约为0.45 m，CA模型模拟水位为0.51 m；t=330 min时，洪湾港实际水深约为0.55 m，CA模型模拟水位为0.50 m，其他特征点的水深也与实际情况基本相符，见图8，误差小于10%，一定程度上反映了模型的可靠性，结果可为风暴潮灾害管理提供洪水和水深数据参考。

图8 基于社交媒体特征数据点的淹没水深对比图Fig.8 Comparison of inundation depth based on social media feature points

2.3 基于云模型的风暴潮灾害多属性综合风险计算

根据土地利用类型分布和洪水最大淹没水深分布，见图9、10。利用公式（10）和公式（12）进行栅格计算，可以得到该地区风暴潮洪水灾害损失分布图和受灾人口分布图，见图11、12。

图9 土地利用类型分布图Fig.9 Distribution map of land use types

图11 经济损失风险图Fig.11 Risk map of economic losses

图10 最大淹没风险图Fig.10 Risk map of the maximum simulated flooding area

基于云模型的指标等级划分如表5所示。首先，确定指标权重系数并通过一致性检验，权重系数为：w=(0.560,0.250,0.096,0.096)。其次，根据云模型理论，借助正向云发生器算法计算生成每个评价指标的隶属度云图，构建灾害等级评价指标的隶属度矩阵R，见式（18），结合指标权重和隶属度矩阵，计算风暴潮灾害风险评价矩阵：B=w×R=(0,0.085,0.007,0.252)，将B归一化得B=(0,0.247,0.020,0.733)。最后，根据置信度准则，取置信度为0.7，可判断此次风暴潮灾害属于特别重大的自然灾害。该方法实现了评价指标与评语的不确定性映射，以此预判的灾害综合风险等级可为城市的风暴潮灾害风险预测和制定防灾策略提供科学依据。

表5 基于云模型的指标等级划分Table 5 Index classification based on cloud model

通过统计区域各土地类型对应不同水深的淹没面积和经济损失，可以得到总的直接经济损失为69.96亿元，间接经济损失为17.49亿元，合计87.45亿元，总伤亡人数约889人。

从经济损失风险图（图11）及受灾人口分布风险图（图12）可以看出此次台风事件在单位面积价值较高的地区、人口密度较大的地区，以及靠近河流、海洋地势较低的沿海地区会引发较高的损失风险，地形、位置、人口、单位面积价值都是导致风暴潮洪水损失的重要因素。同时，本文从权威媒体获取台风天鸽相关信息，与本文模型预测的风险区进行比较分析，辅助验证了风暴潮洪水对珠海市人员、经济、环境造成了巨大风险，结果与模型较为一致，见表6。

图12 人口伤亡分布图Fig.12 Distribution map of population casualties

表6 权威媒体相关信息提取及模型验证Table 6 Extraction of photo information in social media and verification

3 结论与展望

风暴潮灾害的发生具有显著的突发性和短历时性，为了实现实时高效的应急响应，本文首先建立了基于栅格水动力学的元胞自动机模型，通过与传统二维水动力学模型从淹没水深差值、均方根误差（≤0.04 m）、相关系数（≥0.93）、计算时间多个方面对比后，证明了模型在保证精度的前提下可有效提高模拟效率，能够实现利用较少的参数和地形数据近实时地模拟风暴潮洪水演进过程。其次，引入云模型进行了风暴潮灾害的多属性风险分析，可在仿真基础上实时评估灾害综合动态风险及风险的空间分布。最后，以珠海市香洲区作为研究区，对台风天鸽导致的风暴潮洪水进行了模拟，分析得到风暴潮洪水淹没范围、淹没水深，并对灾害损失风险进行了综合评估，研究成果可为风暴潮灾害实时预警、风险分析和制定防灾减灾措施提供一定的技术支持。