基于柵格與云模型的風暴潮洪水風險模擬評估方法
——以珠海市香洲區為例

孫 海，嵇文捷，鄭雅芝

（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）、計算時間多個方面對比后，證明了模型在保證精度的前提下可有效提高模擬效率，能夠實現利用較少的參數和地形數據近實時地模擬風暴潮洪水演進過程。其次，引入云模型進行了風暴潮災害的多屬性風險分析，可在仿真基礎上實時評估災害綜合動態風險及風險的空間分布。最后，以珠海市香洲區作為研究區，對臺風天鴿導致的風暴潮洪水進行了模擬，分析得到風暴潮洪水淹沒范圍、淹沒水深，并對災害損失風險進行了綜合評估，研究成果可為風暴潮災害實時預警、風險分析和制定防災減災措施提供一定的技術支持。