基于元胞自動機的實時降雨下流域洪水淹沒過程分析

蔡木良，支歡樂，李文歡，黃 河，劉 賢，蔣水華，張秀平

（1. 國網江西省電力有限公司電力科學研究院，江西 南昌 310058； 2. 南昌大學工程建設學院，江西 南昌 330031；3. 中鐵水利水電規劃設計集團有限公司，江西 南昌 330029； 4. 江西省水利科學院，江西 南昌 330029）

0 引 言

近年來，結合數字高程模型（DEM）、水動力模型和地理信息系統（GIS）推求洪水淹沒范圍是水利應用領域的研究熱點之一。強暴雨引起的具有威脅性的洪水是我國最為嚴重的自然災害之一，往往會對流域出口及下游緩坡及山間盆地人口密集地帶造成重大的人員傷亡和財產損失。目前我國尚未建立完善的基于降雨資料的洪水災害預評估體系。如何快速預警預報暴雨洪水災害是當前迫切需要解決的難題。

建立一個適用性及可推廣性強，并且具有實際物理意義的洪水淹沒分析模型是實現暴雨洪水災害預警預報的重要前提。傳統的預警預報方法根據臨界雨量等預警指標和結合風險區域地形地貌等特征進行綜合評斷［1］，然而該方法確定預警指標臨界值缺乏科學依據，無法獲得準確的暴雨洪水災害淹沒信息，給決策者制定有效的洪災防災措施帶來了困難。結合GIS和分布式水文模型是目前提倡的洪水災害預警預報方法，可考慮降雨資料的可獲得性以及預警預報技術的可推廣性。其中采用相對簡單并具有一定物理意義的水文模型［2］，是目前使用較為廣泛的方法之一［3-7］。如劉仁義［3］提出了有源和無源淹沒情形下基于種子蔓延的淹沒區計算方法，但是該方法運算效率及計算精度會受到種子點位置及數據精度等因素影響。馮麗麗等［7］盡管基于GIS 技術運用無源淹沒和有源淹沒方法確定洪水淹沒區范圍，但是不能考慮實時降雨因素進行洪災預警分析。綜上，雖然基于GIS技術的無源淹沒算法較為簡便，但是大多難以準確模擬實時降雨造成的洪水淹沒過程。另外常用的洪水淹沒分析方法如水動力學模型，雖然能夠模擬降雨和河道匯流過程得到相對準確的結果，但是需要的輸入參數較多，運算過程復雜，對數據精度要求較高。這些方法容易受到數據觀測與采集手段的限制，在無資料流域或無法確定河道上游斷面流量的情況下具有一定局限性。此外，上述方法建模過程繁瑣，且模型所需文件數據較多，數值求解微分方程的近似解時存在很大困難，導致實際工程應用非常有限［8-10］。

相較之下，元胞自動機為模擬洪水淹沒過程提供了一種有效的分析手段［8-12］。如王復生［8］運用元胞自動機模型模擬流域降雨徑流、產匯流和洪水淹沒的時空變化過程。Dai 等［12］考慮水文因素影響，構建了三層空間連接城市屋頂、地面和人工排水網絡，建立了基于元胞自動機的水文和NPS 污染模型。孫海等［13］提出了基于柵格水動力學的元胞自動機模型，獲得了與MIKE 21 模型相近的計算結果，并有效提高了計算效率。王威等［14］建立了基于“體積法”與元胞自動機的潰壩洪水演進模型，能較好地模擬城市地區潰壩洪水演進的時空特性。劉伊萌等［15］利用CADDIES-2D 洪水模型模擬了北京五環內城區內澇情況，研究發現CADDIES-2D 模型的模擬精度優于Flood Area模型，但稍遜于BUW 模型。綜上，相比于無源淹沒分析方法，基于元胞自動機的洪水淹沒分析模型可以充分考慮流域下墊面條件。相比于水動力學模型，該模型在保證計算精度的前提下需要的模型輸入參數更少，同時元胞自動機受數據的影響更小，不僅可以計算獲得流域內降雨產匯流量，而且可以動態輸出徑流空間變化過程，因而具有較大的應用前景［16］。

以江西省宜春市袁州流域作為研究區域，首先基于GIS 技術生成數字高程模型格網數據，結合流域下墊面條件，求解洪水淹沒范圍及水深，進而建立基于元胞自動機的洪水淹沒分析模型。將實時預報的未來24 h 降水資料作為模型輸入，發生最大暴雨洪水災害的區域作為模型輸出。同時，基于無源淹沒分析方法和MIKE SHE 的坡面流模型計算結果驗證所建模型的有效性。

1 洪水淹沒過程分析方法

在引入基于元胞自動機的洪水淹沒分析模型之前，首先簡要介紹洪水淹沒過程分析的無源淹沒洪水分析方法。

1.1 無源淹沒洪水分析方法

無源淹沒分析方法基本原理是運用GIS 技術遍歷DEM 上每個點的高程值，一旦該點的高程值小于給定洪水水位，即將其加入至淹沒區域。該方法僅將區域總降水量作為輸入，適用于大面積區域降水情況下洪水淹沒分析，計算相對簡單。無源淹沒分析方法將復雜的洪泛區域地形劃分為若干個相等的面積單元，以點（x，y）代表單元所在的地理位置，各單元淹沒水位和地面高程各不相同，可得單元網格對應的淹沒水深D（x，y）為：

式中：H（x，y）為單元網格對應的淹沒水位；E（x，y）為單元網格對應的地面高程值。

考慮到區域下墊面的不均勻性，通過分析地表降雨經過蒸散發、植被截留、土壤下滲等一系列損耗后產生的徑流量過程，可以得到單元網格處由降雨產生的徑流深R（x，y）為：

式中：α（x，y）為單元網格處的徑流系數；P（x，y）為單元網格處的降雨深度。

同時，根據無源淹沒基本理論，由下式可確保洪水來水量與洪水淹沒范圍內總水量相等：

式中：i=1，2，…，n；Ri（x，y）和Di（x，y）分別為第i單元網格處徑流水深和淹沒水深；Δσ為單元網格面積。

1.2 基于元胞自動機的洪水淹沒分析模型

元胞自動機（Cellular Automata）是將系統劃分為規則的元胞網格，每個網格處于有限狀態并受相鄰元胞狀態的影響［8］。在每個模擬步中，所有元胞狀態都根據指定的轉換規則進行同步更新。本文建立的基于元胞自動機的洪水淹沒分析模型，將流域產匯流和洪水淹沒時空分布過程分為3 個階段：①雨水降落到地表時經植被攔截、土壤非飽和帶吸收作用后，產生凈雨；②凈雨經流域坡面、飽和帶、非飽和帶作用，匯聚形成徑流；③最終形成洪水［17，18］。該模型主要由元胞空間、元胞鄰域、元胞狀態和轉換規則這四部分構成，按照既定規則進行動態模擬，主要特點如下：

（1）采用最常見的二維正方形元胞，以確保元胞空間劃分與研究區域地形的DEM格網劃分一致，便于數據提取、轉換、計算和數據前處理以及結果可視化輸出。

（2）元胞鄰域是指中心元胞周邊的其他元胞，與中心元胞有著直接的演化關系。本文研究中鄰域類型選擇摩爾型（Moore），定義每個元胞的左上、上、右上、左、右、左下、下和右下這8 個相鄰元胞為其鄰域元胞（見圖1），其中藍色元胞為中心元胞，相鄰的紅色元胞為鄰域元胞。洪水在360°方向上的演進過程皆為連續，且只與相鄰的水量單位有關。

圖1 摩爾型鄰域關系Fig.1 Moorish neighborhood relationship

（3）元胞狀態是洪水演進過程中的狀態量和參數，需要儲存在元胞空間內進行相應的計算［9］。根據其各自屬性和特征，可分為靜態量集和動態量集。表1給出了需要采用的元胞狀態量參數。

表1 元胞的狀態量參數Tab.1 Parameters for state variable set of a cell

（4）制定元胞轉換規則是構建元胞自動機模型的核心環節。元胞自動機模型中的元胞轉換規則主要有產流規則和匯流規則。這些規則主要用于計算水流的流動方向和進行水量分配。水流方向與中心元胞與其鄰域元胞之間的水位差有關，如果某中心元胞的水位最低，則不進行水量分配。先計算中心元胞及8個領域元胞的水位平均值：

式中：H0為中心元胞水位；Hi為未被去除的鄰域元胞的水位；m為未被去除的鄰域元胞的數目。將水位大于平均值（Hi＞ave）的鄰域元胞去除，計算剩余鄰域元胞與中心元胞的水位平均值，并繼續剔除水位大于平均值的鄰域元胞，重復上述步驟直到沒有鄰域元胞被剔除，即剩余的鄰域元胞水位都低于平均值。接著，再選擇剩余的鄰域元胞作為該中心元胞水量分配的對象，使中心元胞與剩余鄰域元胞具有相同的水位。按照最小差異算法，中心元胞向剩余鄰域元胞進行水量分配，直至達到平均水位。元胞水流流速采用以下曼寧公式推求：

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

據此，可進一步得洪水由中心元胞流向鄰域元胞的流動時間為：

式中：L為中心元胞與鄰域元胞的距離。在匯流模擬過程中，時間步長t的大小決定了流量值。

鄰域元胞的實際流量q為：

在進行水量分配時，需要注意以下兩點：①當計算所得的分配水量大于該元胞水深時，分配水量需要按照該元胞的水深進行計算；②當中心元胞與鄰域元胞水位相同時，則不再進行水量分配。為保證模型運行效率和穩定性，通常將時間步長t設置為小于大部分洪水流動時間的值，以防止水流在某一個時間步內穿過元胞。每個元胞將水量分配給其鄰域元胞的同時，也會獲取來自其他鄰域元胞的洪水，因此每個時間步長的元胞洪水量Qi計算公式為：

式中：qi表示中心元胞在該時間步長內的洪水流入量，來自其周邊8 個鄰域元胞；qt表示中心元胞在該時間步長內的洪水流出量；Qi-1表示通過上個時間步計算的元胞洪水量。

2 案例分析

下面依托江西省宜春市袁州流域，采用所建的基于元胞自動機的洪水淹沒分析模型模擬洪水淹沒過程。江西省宜春市袁州地處江西省西部，東鄰鋼城新余，南接古城吉安，西毗煤都萍鄉，西北連湖南腹地。DEM 提取的流域數字特征包括單元格網的流向、匯流路徑、河網間的拓撲結構等。江西省宜春市袁州流域DEM高程圖如圖2所示。

圖2 江西省宜春市袁州流域圖Fig.2 Diagram of Yuanzhou watershed in Yichun City，Jiangxi Province

圖3 給出了以2021 年5 月28 日12 時至5 月29 日12 時的總降雨量生成的江西省預報降雨數據分布圖。由圖3 可知，該時段內袁州流域地區總降雨量高達78 mm，極有可能出現因短歷時強降雨誘發的洪水災害。研究流域采用30 m×30 m 格網，以2021年5月28日12時為計算初始時刻，采用30 m×30 m 分辨率的DEM高程格網數據和降雨格網數據，其中高程格網數據和降雨格網數據保持一致。

圖3 江西省2021年5月28日降雨數據分布圖Fig.3 Distribution of rainfall data in Jiangxi province on May 28， 2021

2.1 模型計算結果

采用所建的基于元胞自動機的洪水淹沒分析模型計算降雨產流量時，考慮植被截留量和土壤下滲量這兩種降雨損耗，采取徑流系數表示降雨產生徑流的損耗過程。根據江西省暴雨洪水查算手冊［19］以及研究區域流域特征，取該研究區域的徑流系數為0.57。進行匯流計算時，需要輸入每個元胞DEM 的格網數據和產生的徑流量數據。本文基于DEM 格網數據中3×3的元胞矩陣進行計算，首先篩選出徑流流動方向及流量，求解單位計算步的中心元胞至領域元胞的流動水量，得到規定時間內洪水淹沒范圍及淹沒深度。將降雨發生前即模型開始計算時刻的研究區域淹沒面積視為0。為了更直觀地呈現降雨徑流過程中流域洪水淹沒及變化過程，計算過程中記錄各個時刻的元胞水位值，并以柵格數據的形式進行可視化輸出，進而可得到4 個不同時刻（10 min、30 min、1 h 和2 h）的宜春市袁州流域洪水淹沒數據，如圖4 所示。表2 給出了不同時刻對應的淹沒面積。另外，根據黃中發［20］和表3 基于淹沒水深的洪澇災害等級劃分，將研究區域可分為輕澇、中澇、重澇和特澇這4 個等級。

表2 宜春市袁州流域淹沒面積km2Tab.2 Inundation area of Yuan Zhou watershed in Yichun city

表3 洪澇災害等級劃分Tab.3 Classification of flood disaster levels

圖4 不同時刻的流域洪水淹沒水深結果Fig.4 Flood water inundation results for different durations using the established model

由表2 可知，這4 個時刻宜春市袁州流域的洪水總淹沒面積分別為136.28、40.91、33.68、33.68 km2。10 min時最大淹沒水深達4.93 m，淹沒區域分布較為密集，主要是由于此時洪水主要匯集在山谷溝壑等區域。隨著模擬時間的延長，輕澇淹沒面積和總淹沒面積逐漸降低，30 min 時最大淹沒水深逐漸增加至8.90 m，這是因為山谷溝壑等區域的洪水逐漸匯集到河道區域。由圖4（c）可知，1 h 時河道支流洪水逐漸減少，洪水貫通區逐漸增大，最大淹沒水深逐漸增加至10.18 m。此外，由圖4（c）和（d）可知，模擬時間超1 h 后，流域內淹沒范圍基本不變，最大淹沒水深逐漸增加至10.62 m 后逐漸減小。此外，洪水并非全部聚集在地勢低的區域，而是有條件地選擇行進路徑，即流入地形高程比其周圍淹沒水位更低的區域，符合水流運動的基本規律。綜上可知，所建模型能夠有效模擬洪水淹沒過程及時空分布特性。

2.2 模型驗證

本節分別采用無源淹沒方法和基于MIKE SHE 的坡面流模型計算結果驗證所建模型的有效性。首先，采用無源淹沒分析方法基于宜春市袁州流域的30 m×30 m 高程點的格網數據E（x，y）以及降雨格網數據P（x，y），通過式（2）計算單個格網的徑流深R（x，y），其中徑流系數取值同樣取0.57。進而根據洪水來水量與洪水淹沒范圍內總水量相等這一原則，聯立式（1）和式（3）推求研究區域淹沒水位。獲得宜春市袁州流域淹沒水位值之后，構建無源淹沒洪水區域輸出模型。接著依次對格網單元進行掃描，將格網單元的高程和給定洪水位進行比較，并把其中的淹沒單元進行連接，計算淹沒水深。最后獲得最終的洪水總淹沒面積為6.73 km2，最大淹沒水深為13.35 m（見圖5）。然而，利用無源淹沒洪水分析方法得出宜春市袁州流域淹沒災害等級最高區域主要集中在流域西南部地勢較低處，顯然與實際洪水災害特征不符。這是由于無源淹沒分析方法的基本原理是通過遍歷所有柵格單元，從中找出淹沒柵格，雖然實現過程較為簡單，但是不能考慮地域連通性，導致高程值低于淹沒水位H的柵格都被計入淹沒區域。

圖5 無源淹沒分析方法的淹沒水深結果Fig.5 Inundation water depth results using the non-source inundation analysis method

接著，采用MIKE SHE 的坡面流模型進行洪水淹沒過程分析。坡面流模塊主要采用擴散波對圣維南方程組進行近似分析，徑流系數同樣取值為0.57，糙率初始取0.025。MIKE SHE的坡面流模型運行時間與元胞自動機洪水淹沒模型相同。圖6給出了模擬得到的坡面流淹沒水深分布。由圖6可知，MIKE SHE的坡面流模型獲得的總淹沒面積達34.56 km2，與圖4（d）中歷時2 h 的總淹沒面積僅相差0.88 km2，較為吻合。圖6 和圖4（d）的淹沒水深分布基本一致表明，所建模型能夠準確模擬洪水淹沒過程。表4 進一步比較了3 種不同方法洪水淹沒結果。由表4可知，相較于無源淹沒分析方法，所建模型計算精度更高；相比于MIKE SHE 的坡面流模型，所建模型在滿足計算精度的同時，無需進行復雜的模型計算，更便于實際工程推廣應用。

表4 不同方法洪水淹沒結果的對比Tab.4 Comparison of flood inundation results obtained from different methods

圖6 MIKE SHE坡面流模型的淹沒水深計算結果Fig.6 Inundation water depth results using MIKE SHE overland flow model

2.3 參數敏感性分析

為探討輸入參數的影響，本節采用所建的元胞自動機洪水淹沒分析模型對糙率n進行敏感性分析，調查不同糙率n對研究區域不同時刻（30 min、1 h 和2 h）最大淹沒水深的影響。糙率n分別取0.01、0.025 和0.04，即是在基準值0.025 的基礎上分別增加和減少0.015，其他參數保持不變［20］。計算參數敏感性指標和H0.04分別表示糙率n取0.01、0.025和0.04的最大淹沒水深。參數敏感性指標α值越大，表示該參數對研究區域最大淹沒水深的影響越大。表5給出了糙率n的敏感性分析結果。由表5 可知，隨著模擬時間的延長，糙率n對研究區域最大淹沒水深的影響逐漸減小。因此，對于短時間洪水淹沒過程，糙率n的敏感性較大，而隨著模擬時間的延長，糙率n的敏感性逐漸減小。

表5 糙率n的敏感性分析結果Tab.5 Sensitivity analysis results of roughness n

3 結 論

依托江西省宜春市袁州流域，建立了基于元胞自動機洪水淹沒分析模型，進而基于實時降雨數據獲得了研究區域不同時刻的淹沒水深和淹沒范圍，并通過無源淹沒分析方法和MIKE SHE 的坡面流模型計算結果驗證了所建模型的有效性。主要結論如下：

（1）無源淹沒分析方法雖然計算簡單，但是不僅只能考慮總降雨量對研究區域的淹沒結果的影響，而且忽略了研究區域下墊面及地形連通性的影響，無法考慮降雨產匯流過程，獲得的洪水淹沒計算結果與實際情況不符，不利于洪災預警措施方案的制定。

（2）建立的基于元胞自動機的洪水淹沒分析模型可以充分考慮研究區域下墊面條件以及降雨的產匯流過程的影響，更加合理地模擬洪水淹沒過程，動態顯示不同時刻的洪水淹沒過程。相較于無源淹沒分析方法，所建模型具有更高的計算精度；相比于MIKE SHE 的坡面流模型，所建模型在滿足計算精度的同時，無需進行復雜的模型計算，更便于實際工程推廣應用，從而可為強降雨下洪水災害風險評估以及防汛搶險決策方案的制定提供了參考。

（3）所建模型的元胞大小和計算效率成反比，為保證計算效率，本文模型采用的元胞邊長較大（為18 m），在一定程度上會影響對真實地形條件的模擬精度。此外，本文模型只考慮了中心元胞向鄰域一個元胞進行分配水量。為更真實反映工程實際情況，關于中心元胞向鄰域多個元胞的水量分配問題還有待進一步研究。