基于元胞自動機的城市雨洪模型研究

李繼政，李傳奇，王復生，張焱煒

(山東大學土建與水利學院，濟南 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)模型產流規則中參數較多，本文選取經驗參數進行計算，未對其進行敏感性分析，一些參數優化方法如人工神經網絡等可以適用于多參數的優化計算，這對于提高元胞自動機模型精度來說也是一種不錯的方法。