基于元胞自動機(jī)的城市地表徑流流向優(yōu)化算法

何夢男，陳 誠，李 港，胡智華，何湖濱，陳求穩(wěn)，6

(1.南京水利科學(xué)研究院生態(tài)環(huán)境研究所，江蘇 南京 210029； 2.四川大學(xué)水利水電學(xué)院，四川 成都 610044；3.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098； 4.重慶交通大學(xué)河海學(xué)院，重慶 400074；5.浙江省交通規(guī)劃設(shè)計研究院有限公司，浙江 杭州 310030；6.水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210029)

近年來，快速的城市化進(jìn)程在改變地表下滲過程的同時增加了各種復(fù)雜陡、緩地形，使得城市地表水力梯度組成紛繁復(fù)雜。城市水文模型可為城市內(nèi)澇預(yù)警和面源污染防控提供重要的決策工具，其發(fā)展需要適應(yīng)城市地表復(fù)雜變化的特征，從而實現(xiàn)城市降雨徑流過程的精確模擬[1-3]。其中，水流流向是城市水文模型中表征地表徑流的重要參數(shù)，描述了水流在不同地表地貌單元之間從高到低的流動路線。流向算法是眾多城市水文模型中計算水流流向的基本方法，影響著匯水面積、地形指數(shù)等參數(shù)的計算[4]，在城市水文模型研究中起著重要作用。

目前流向算法主要分為單流向(single flow direction，SFD)算法和多流向(multiple flow direction，MFD)算法，均假設(shè)局部水力梯度可用局部坡度估計[4-5]。單流向算法是將水力梯度最大的相鄰單元作為水流的下游單元，具有計算簡單、效率高等優(yōu)勢[6]。單流向算法往往趨向于形成收斂的直線，而實際情況中水流有可能是漫散流態(tài)，因而單流向算法一般適合模擬水流在陡坡面上的流動，不適合模擬緩坡面上的漫散流動。而多流向算法的方法學(xué)原理主要基于“越陡下游方向獲取的水量越多”和“動力系統(tǒng)傾向于往平衡狀態(tài)發(fā)展”兩種假設(shè)[7]，將水量按照分配權(quán)重p分配至不同的下游單元中。水量分配權(quán)重是指水量分配到各個網(wǎng)格的權(quán)重，用于描述水流從當(dāng)前時刻到下一時刻的運動過程。有關(guān)學(xué)者基于水量分配權(quán)重的不同取值，提出了FD8、FMFD、MFDfg等多流向算法[8-9]。雖然多流向算法可以彌補單流向算法模擬緩坡漫散流的缺陷，但目前多流向算法中水量分配權(quán)重多為固定值，忽略了地形特征對流場劃分的影響，在坡度較陡的情況下模擬結(jié)果往往較差[10]。因此，如何兼顧單流向算法與多流向算法的優(yōu)點，構(gòu)建適用于城市地表復(fù)雜水力梯度的徑流優(yōu)化算法對于城市水文過程模擬具有重要意義。

元胞自動機(jī)(Cellular Automata，CA)是大量個體通過局部聯(lián)系組成的時間、空間以及狀態(tài)均離散的數(shù)學(xué)方法。在水文學(xué)和土壤侵蝕研究領(lǐng)域，通常被用于快速洪澇[11-13]和泥沙沖淤演變模擬[13]，已被實踐證明是一種適合復(fù)雜系統(tǒng)模擬的離散化數(shù)學(xué)工具，不僅能避免微分方程求解，還可有效減少數(shù)值計算的復(fù)雜度。

本文在考慮城市地表復(fù)雜水力梯度特征的前提下，提出一種基于元胞自動機(jī)的城市徑流流向優(yōu)化算法，通過計算元胞鄰域各單元水位值標(biāo)準(zhǔn)差，得到直接反映局部水力梯度的指標(biāo)，并設(shè)計流態(tài)判斷規(guī)則以彌補其他多流向算法的不足。

圖1 Moore型鄰域Fig.1 Moore type neighborhood

1 材料與方法

1.1 元胞自動機(jī)模型

元胞自動機(jī)模型主要包括元胞空間、鄰域、狀態(tài)和規(guī)則四部分，其中元胞規(guī)則是元胞自動機(jī)的核心，主要是流態(tài)判斷規(guī)則，定義了元胞從當(dāng)前時刻狀態(tài)向下一時刻演變的法則，決定了元胞自動機(jī)模型的演化方向。元胞鄰域采用典型的Moore型鄰域，如圖1所示。元胞狀態(tài)除了元胞索引坐標(biāo)(x，y)、邊長l以及時間步長t以外，還包括表征產(chǎn)匯流過程的地表特征、徑流兩部分狀態(tài)，見表1。

表1 元胞狀態(tài)參數(shù)

1.2 流向優(yōu)化算法

基于元胞自動機(jī)的城市徑流流向優(yōu)化算法(the optimal urban flow direction algorithm based on Cellular Automata，UCA)流程如圖2所示，主要分為3個步驟：

圖2 UCA算法流程Fig.2 Flowchart of UCA algorithm

步驟1轉(zhuǎn)移時間計算。比較中間元胞與鄰域元胞的水位值，排除水位值高于中間元胞的鄰域元胞，將剩下的鄰域元胞作為可轉(zhuǎn)移水量的下游元胞，再根據(jù)式(1)計算中間元胞水量全部轉(zhuǎn)移到各個下游元胞所需時間：

(1)

步驟2流態(tài)判斷。設(shè)置流態(tài)判斷規(guī)則對漫散流態(tài)與匯聚流態(tài)進(jìn)行定量區(qū)分，主要根據(jù)9個元胞高程(水深)值標(biāo)準(zhǔn)差(SD)與標(biāo)準(zhǔn)差閾值(SD0)的大小進(jìn)行判斷：

(2)

(3)

步驟3轉(zhuǎn)移水量計算。先根據(jù)式(4)(5)確定各個下游元胞的水量分配比例di，再根據(jù)式(6)對比每個元胞轉(zhuǎn)移時間與模型時間步長的大小，計算時間步長內(nèi)轉(zhuǎn)移至各個下游元胞中的水量Ri。

(4)

(5)

(6)

(7)

1.3 流向算法試驗

圖3 兩種坡面示意圖Fig.3 Illustrations of two slopes

為篩選漫散流態(tài)和匯聚流態(tài)時各自水量p的最優(yōu)值，考慮如圖3所示的陡、緩坡面流問題：x方向和y方向長度分別為500 m和100 m，坡度β分別為0.20和0.05。穩(wěn)定、均勻的超滲降雨落于坡面，降雨強度I為100 mm/h，降雨歷時為25 min，模擬時長為60 min，坡面水力糙率采用曼寧系數(shù)表示，取值為0.01。假定初始時刻坡面處于干燥狀態(tài)，l=1.0 m，采用UCA算法求解坡腳出口處的流量過程線，并利用羅志等[15]提出的坡面流解析解(analytical solution)作為參考進(jìn)行定量評價。

為定量評價UCA算法在不同地形條件下的性能，采用不同算法模擬計算了城市中常見的凸坡、凹坡、山脊和直坡上的單位匯流面積(specific catchment area，SCA)。單位匯流面積指單位長度等高線上的匯流面積，是定量描述地表水流分布的重要水文參數(shù)，通?？梢詫Ρ萐CA的模擬值與理論值對流向算法進(jìn)行定量評價。采用l=1.0 m及4種人工數(shù)學(xué)曲面[5-6，9](橢球面、反橢球面、馬鞍面和斜面)分別代表城市中常見的凸坡、凹坡、山脊和直坡。人工數(shù)學(xué)曲面的優(yōu)點是可以通過數(shù)學(xué)推導(dǎo)獲得匯流面積的理論值，其定義及對應(yīng)的SCA理論值如圖4所示。參與對比的算法有單流向D8算法、固定水量分配權(quán)重的FMFD算法、UCA算法和基于差異最小化的4+4N算法。

圖4 人工數(shù)學(xué)曲面及其單位匯流面積理論值Fig.4 Artificial mathematical surface and its theoretical value of SCA

采用均方根誤差(root mean square error，RMSE)、納什效率系數(shù)(Nash-sutcliffe efficiency coefficient，NSE)和平均絕對誤差(mean absolute error，MAE)對模擬結(jié)果進(jìn)行評價。

2 結(jié)果與分析

2.1 水量分配權(quán)重篩選

在坡度較緩時，Quinn等[16-17]認(rèn)為水量p=1.0可以模擬水流完全漫散的情況，F(xiàn)reeman[17]通過對比p=1.0、1.1和1.25情況下的模擬效果，認(rèn)為p=1.1可作為緩坡漫散流態(tài)情況下最優(yōu)水量分配權(quán)重。因此，設(shè)置水量p=1.1、2.0、5.0、8.0、10.0等5種工況下的坡面流模擬試驗。不同水量分配權(quán)重下出口流量模擬統(tǒng)計結(jié)果見表2。由表2可以看到，不同水量分配權(quán)重出口流量陡坡條件下的納什效率系數(shù)都很好，接近于1，而均方根誤差有明顯的區(qū)別。按照均方根誤差越接近于0，納什效率系數(shù)越接近于1的選擇依據(jù)，可看出在城市地區(qū)陡坡匯聚流態(tài)和緩坡漫散流態(tài)的最優(yōu)水量分配權(quán)重分別是2.0和1.1。

表2 不同水量分配權(quán)重下出口流量統(tǒng)計結(jié)果

2.2 匯流面積差異分析

圖5～圖8分別為凸坡、凹坡、山脊坡以及直坡面上不同流向算法求解的單位匯流面積與其理論值差異。對比4種算法在各坡面上匯流面積差異結(jié)果，可看到D8和4+4N算法單位匯流面積與理論值有差異的區(qū)域面積明顯大于FMFD算法，說明D8和4+4N算法的單位匯流面積模擬結(jié)果與理論值差異相對較大，而FMFD算法的單位匯流面積模擬結(jié)果與理論值差異小于D8和4+4N算法。同時，對比UCA算法和FMFD算法的各坡面匯流面積差異結(jié)果，可以看到UCA算法和FMFD算法在4種地形條件下模擬結(jié)果最接近，且在凹坡面上UCA算法的結(jié)果與理論值有差異的區(qū)域與FMFD算法相當(dāng)，而在凸坡面、山脊和直坡面上UCA算法的結(jié)果與理論值有差異的區(qū)域面積比FMFD算法更小，說明UCA算法在FMFD算法的基礎(chǔ)上進(jìn)一步實現(xiàn)了優(yōu)化。

圖5 凸坡面上不同流向算法求解的單位匯流面積與其理論值差異Fig.5 Difference between theoretical SCA and calculated SCA by different flow direction algorithms on convex surface

圖6 凹坡面上不同流向算法求解的單位匯流面積與其理論值差異Fig.6 Difference between theoretical SCA and calculated SCA by different flow direction algorithms on concave surface

圖7 山脊坡面上不同流向算法求解的單位匯流面積與其理論值差異Fig.7 Difference between theoretical SCA and calculated SCA by different flow direction algorithms on ridge

圖8 直坡面上不同流向算法求解的單位匯流面積與其理論值差異Fig.8 Difference between theoretical SCA and calculated SCA by different flow direction algorithms on plane surface

2.3 算法性能定量評價

圖9為不同流向算法求解的單位匯流面積與其理論值的誤差頻率分布。相比之下可以看到每個坡面上UCA算法的相對誤差分布最集中于0附近，相對誤差值在[-20，20]之間的元胞比例分別為99.84%、89.01%、99.89%和99.80%，在[-5，5]之間的元胞比例分別為89.61%、68.22%、84.72%和80.92%。FMFD算法次之，而D8算法和4+4N算法的誤差頻率相對分散?？傮w來看，UCA算法模擬結(jié)果與理論值的相對誤差最小。

圖9 不同算法求解的單位匯流面積與其理論值的誤差頻率分布Fig.9 Error frequency distribution of theoretical SCA and calculated SCA by different algorithms

表3 不同算法定量評價結(jié)果

為進(jìn)一步定量評價算法的性能，統(tǒng)計不同算法模擬單位匯流面積的均方根誤差和平均絕對誤差，結(jié)果見表3。從不同坡面的兩種誤差結(jié)果來看，4種算法在不同坡面上的誤差從大到小均呈現(xiàn)凹坡大于直坡大于山脊大于凸坡的趨勢。4種徑流算法在凸坡面的模擬效果均最好，而在凹坡面(反橢球面)上的均方根誤差最大，這是由于反橢球面底部的徑流過程屬于匯流填洼，流向算法對該復(fù)雜流向區(qū)域的概化本身就存在誤差[5]，再加上反橢球面底部等高線長度無窮小，誤差又被無窮小的等高線進(jìn)一步放大。

城市地形中常見的路面、屋面等多屬于直坡面[14，18]。從表3每種算法的誤差結(jié)果來看，D8算法在直坡面上的誤差均較大，而在凸坡面和山脊面上的誤差相對較小，說明D8算法適用于模擬凸坡、山脊這類陡坡匯聚的徑流過程，而不適用于模擬直坡這類緩坡漫散的徑流過程。4+4N算法相對于D8算法增加了在對角線方向分配水流的可能，因而在4種坡面上的誤差相比于D8算法有所減小。FMFD算法較D8和4+4N算法的誤差進(jìn)一步顯著減小，說明FMFD算法固定水量分配權(quán)重的方式可以有效彌補D8和4+4N算法不適用于模擬緩坡漫散流的缺陷。

UCA算法在FMFD算法的基礎(chǔ)上考慮了局部地形特征變化，采用水量分配權(quán)重隨地表特征進(jìn)行變化的方式處理水流多流向分配問題。從表3還可以看到UCA算法在4種坡 面上的均方根誤差相比于FMFD算法進(jìn)一步減小，特別是對于直坡面上的漫流情況，UCA算法的均方根誤差為4.2 m2/m，相比于FMFD算法12.0 m2/m減小了近3倍，說明UCA算法不僅適合模擬城市平坦地形上的緩坡漫散流，還適合模擬城市局部地形變化較大坡面上的收斂匯聚流。從平均絕對誤差結(jié)果來看，UCA算法和FMFD算法在凹坡地形單位匯流面積模擬結(jié)果中，與其理論值的平均絕對誤差很接近，分別為15.0 m2/m和14.4 m2/m；而在凸坡、山脊和直坡上，UCA算法的單位匯流面積模擬結(jié)果與理論值的平均絕對誤差均最小，分別為2.8 m2/m、4.3 m2/m和10.5 m2/m。總體而言，UCA算法比FMFD算法更接近單位匯水面積理論值。

3 討 論

城市道路最大縱坡度受設(shè)計車速等因素限制，一般不低于8.0%；草地設(shè)計坡度通常較緩，一般在1.5%～12.0%；城市屋面常設(shè)為平坡面，坡度一般小于5.0%[18]。這些常見的下墊面表面徑流過程趨于漫散，而在坡度較大的坡面或不同下墊面之間過渡時，由于局部坡度較大，徑流過程趨向于匯聚收斂。準(zhǔn)確區(qū)分水流流態(tài)是進(jìn)行水量合理分配的前提，然而，目前研究領(lǐng)域?qū)τ诔鞘械乇硭μ荻鹊亩妇徧卣鳑]有明確的判斷標(biāo)準(zhǔn)。

本文在UCA算法中設(shè)計了流態(tài)判斷規(guī)則，提出流態(tài)判斷的標(biāo)準(zhǔn)差閾值。該流態(tài)判斷規(guī)則綜合考慮了元胞尺寸與鄰域內(nèi)9個元胞高程(水位)值對流態(tài)的影響，可以作為未來研究中城市地表水力梯度陡緩特征的判斷依據(jù)。采用UCA算法進(jìn)行不同水量分配權(quán)重下的坡面流模擬試驗，結(jié)果顯示UCA算法適用于城市地表緩坡漫散流態(tài)模擬的最優(yōu)水量分配權(quán)重pd=1.1，這與Freeman對緩坡漫散流態(tài)模擬的研究結(jié)論一致，說明UCA算法設(shè)計的流態(tài)判斷規(guī)則具有一定的合理性。同時，在陡坡匯聚流態(tài)模擬的研究中，Qin等[8，16，19-20]認(rèn)為在溝谷地形匯聚流動的水量分配權(quán)重pc=10較為合理，而UCA算法適用于城市地表陡坡匯聚流態(tài)模擬的最優(yōu)水量分配權(quán)重pc=2.0，說明相比于天然流域，城市地形不會出現(xiàn)類似于溝谷地形中的匯聚流態(tài)，即pc取值不會過大。

元胞自動機(jī)可以通過定義簡單的規(guī)則演繹復(fù)雜的動力系統(tǒng)，對于城市水文過程模擬具有獨特的算法優(yōu)勢[12-13]。UCA算法模擬效果相對于傳統(tǒng)流向算法更好，適用于分布式網(wǎng)格計算。未來研究可基于UCA算法開發(fā)適用于城市復(fù)雜地形變化的水文模型，為城市水利水務(wù)部門的洪澇預(yù)警、面源污染削減等提供更準(zhǔn)確的信息。本文的不足之處在于缺乏實測資料，算法評價選擇在4種具有代表性的人工數(shù)學(xué)曲面上開展，在實際應(yīng)用過程仍有值得改進(jìn)之處。未來研究可對真實城市地形展開模擬，并用實測匯流數(shù)據(jù)作為驗證條件來改進(jìn)算法。同時，地下管網(wǎng)匯流過程是城市徑流的重要過程，后續(xù)研究可以考慮在元胞自動機(jī)模型中增加地下管網(wǎng)匯流模塊，實現(xiàn)地表和地下管網(wǎng)匯流的耦合計算。

4 結(jié) 論

a.UCA算法中適用于城市陡坡匯聚流態(tài)和緩坡漫散流態(tài)模擬的最優(yōu)水量分配權(quán)重取值分別為2.0和1.1，可為后續(xù)多流向算法的水量分配權(quán)重取值提供參考。

b.UCA算法在FMFD算法的基礎(chǔ)上考慮了局部地形特征變化，采用水量分配權(quán)重隨地表特征進(jìn)行變化的方式處理水流多流向分配問題，其單位匯流面積模擬結(jié)果與理論值差異更小。

c.UCA算法兼顧了單流向與多流向算法的優(yōu)勢，既能模擬城市平坦地形上的緩坡漫散流，也能模擬城市水力梯度較大地形上的陡坡匯聚流，可為城市地表徑流過程模擬提供更優(yōu)的求解思路。