基于遺傳算法的SWMM模型參數(shù)率定研究

馬彥斌，盛 旺，李江云, 代文江

(武漢大學(xué)，武漢 430000)

0 引 言

美國環(huán)境保護署(EPA)開發(fā)的SWMM(Storm Water Management Model)降雨徑流模型，廣泛應(yīng)用于對城市區(qū)域降雨徑流水量、水質(zhì)的動態(tài)模擬[1]。SWMM模型參數(shù)繁多，且部分參數(shù)或者具有不確定性，或者為概化結(jié)果，不具實測意義，無法直接測量，一般通過收集可直接觀測數(shù)據(jù)進行模型參數(shù)率定。即通過尋找一系列適合的模型參數(shù)，使模型預(yù)測結(jié)果盡可能地接近監(jiān)測數(shù)據(jù)[2]。

早期關(guān)于參數(shù)率定方法較少，以人工試錯法和單參數(shù)敏感性[3]分析法為主，由于人工試錯法效率低下，主觀性強。隨著計算機技術(shù)的發(fā)展，一系列自動尋優(yōu)算法[4,5]被運用在降雨徑流模型的參數(shù)識別中。采用自動尋優(yōu)算法時，如果直接對模型中所有參數(shù)進行率定分析，模型計算量會隨著參數(shù)維數(shù)的增加而呈指數(shù)增加，陷入“維數(shù)災(zāi)難”[6]。因此，本文將按率定參數(shù)選擇和模型參數(shù)率定兩個步驟進行研究。

1 研究方法

首先采用敏感性分析法確定敏感參數(shù)，模型參數(shù)的率定僅針對敏感性參數(shù)進行，借此達到降維效果，這在一定程度上降低了模型率定的復(fù)雜性，提高了模型率定的效率[7]；其次，利用收集到的研究區(qū)域?qū)崪y數(shù)據(jù)，基于MATLAB中的遺傳算法模塊，對初步構(gòu)建的城市水文模型進行參數(shù)率定。

1.1 研究區(qū)域概況

研究區(qū)域位于中山市火炬開發(fā)區(qū)張家邊某居民小區(qū)。該小區(qū)于2007年建成，整體綠化率較高，區(qū)域總面積約為6.14 hm2。如圖1所示，根據(jù)研究區(qū)域影像圖，采用遙感軟件ENVI5.1解析可得到下墊面分析結(jié)果。各下墊面類型中，路面所占比例最高，占研究區(qū)域總面積的43.97%，其次為建筑屋面以及綠地，分別占31.92%和24.11%。按照下墊面透水能力進行用地類型劃分，可以把除去綠地以外的下墊面類型劃分為不透水表面，則該區(qū)域的整體不透水表面所占總面積比例為75.89%。

圖1 研究區(qū)域下墊面分析

模型概化的過程主要包括模型中的節(jié)點、管線以及子匯水區(qū)的概化過程。模型概化后共包括21個節(jié)點，其中包括20個檢查井和1個排放口，包括20根管渠以及27個子匯水區(qū)，子匯水區(qū)包括建筑屋面12個，廣場2個，市政路面3個，小區(qū)路面10個。模型概化示意圖及監(jiān)測位置如圖2所示。

圖2 研究區(qū)域模型概化示意圖

1.2 參數(shù)選擇及取值

SWMM模型對降雨-徑流過程的模擬主要包括地表產(chǎn)匯流和管道輸水兩個過程。本研究采用霍頓(Horton)滲透模型來描述透水地表產(chǎn)流過程，采用非線性水庫法描述地表匯流過程，采用圣維南(Saint-Venant)方程進行管渠水流計算。

根據(jù)參數(shù)獲取的方式，SWMM模型中參數(shù)分為可測量參數(shù)和不可測量參數(shù)。本文中可測量參數(shù)有節(jié)點底標高(Invert Elevation)、節(jié)點埋深(Max Depth)、子匯水區(qū)面積(Area)、管道長度(Length)、管道直徑(Diameter)和子匯水區(qū)不透水比例(Pct-Imperious)，這6個參數(shù)通過資料收集或下墊面解析方法進行確定，在模型率定時固定可測量參數(shù)。不可測量參數(shù)則只能通過經(jīng)驗估計或者實測率定而進行獲取，參考SWMM模型手冊及相關(guān)文獻[8,9]，表1給出了模型中不可測量參數(shù)的取值范圍。

表1 SWMM不可測量參數(shù)及其取值范圍

參數(shù)寬度修正因子Kwidth和坡度修正因子Kslope分別用于計算子匯水區(qū)的特征寬度及特征坡度：

(1)

式中：Area為子匯水區(qū)面積。

Slope=KslopeSlope0

(2)

式中：Slope0可以通過原有的地形圖獲取，作為識別特征坡度的初始取值。

1.3 Morris全局敏感性分析法

Morris法是一種基于篩選分析的全局敏感性分析方法[8]。其計算量較少，易于操作，并且該方法不會出現(xiàn)將重要參數(shù)遺漏的錯誤[9]，可用來篩選出那些不敏感的參數(shù)，從而保留敏感參數(shù)，這種定性篩選的作用正適用于本研究的研究目的。因此，本文采用Morris法對SWMM模型進行敏感性分析。一個典型的Morris設(shè)計圖如圖3所示，由圖可看出Morris法是利用OAT(one factor at a time)的采樣方法，即每次都只改變樣本中的某一個參數(shù)。生成的樣本在實際模擬計算的過程中映射到參數(shù)范圍的分布之中，運行模型計算，設(shè)Y為模型的輸出結(jié)果，則各個參數(shù)的基效應(yīng)(Elemenary Effect，EE)為：

圖3 Morris設(shè)計原理

(3)

式中：di(j)表示第i個參數(shù)第j組參數(shù)的基效應(yīng)，j=1,2，…，r(r為重復(fù)次數(shù))；Δ為第i個參數(shù)的擾動幅度。

通常表征Morris法的計算指標為“基效應(yīng)的平均值”，公式如下：

(4)

參數(shù)的平均值表征參數(shù)的敏感性，平均值越大表明該參數(shù)對結(jié)果輸出的影響越大，即越敏感，從而可以確定參數(shù)敏感性排序。

1.4 遺傳算法

遺傳算法是一種隨機搜索和優(yōu)化的方法，通過模擬大自然的生物進化過程，在種群不斷進化過程中逐步淘汰不適應(yīng)個體，而迭代中的種群進行不斷交叉、重組、變異又能產(chǎn)生新的個體。遺傳算法最重要的3個操作參數(shù)是：選擇(Selection)、交叉(Crossover)、變異(Mutation)，可利用MATLAB中遺傳算法工具箱函數(shù)完成該工作，本文中將運用MATLAB中遺傳算法工具箱函數(shù)進行程序編寫。

2 結(jié)果分析

2.1 Morris敏感性分析法計算結(jié)果

采用2016年7月10日實測降雨(73.6 mm)作為模型輸入的邊界條件，同時，本研究選取在城市降雨徑流模擬中具有重要意義的3個輸出變量：總徑流量、峰值流量、峰現(xiàn)時間。敏感性計算結(jié)果如表2所示。

表2 各參數(shù)敏感性分析計算結(jié)果

對總徑流量而言，敏感性排名前4的參數(shù)分別是Max Rate、Min Rate、Decay、N-Imperv，其他參數(shù)對徑流總量的影響較小，這可能是因為在此場降雨條件下，管道曼寧系數(shù)Manning-N的取值基本對總徑流量沒有影響。對于流量峰值而言，敏感性排名前4的參數(shù)依次是Manning-N、N-Imperv、Kwidth、N-Perv，其他參數(shù)對峰值流量的影響均不大；對于峰值時間而言，其敏感性排序結(jié)果基本和峰值流量類似，敏感性排名前4的參數(shù)依次是Manning-N、N-Imperv、Kwidth、S-Imperv。

2.2 遺傳算法步驟及計算結(jié)果

由前述敏感性分析可知，選取4個在對總徑流量、峰值流量和峰值時間敏感性較高的參數(shù)作為率定參數(shù)，這4個參數(shù)是：Max Rate、Min rate、Manning-N、N-Imperv，參數(shù)的取值范圍見表1。

利用MATLAB遺傳算法進行參數(shù)率定的主要框架圖如圖4所示。

圖4 遺傳算法率定參數(shù)框架設(shè)計

算法具體計算步驟如下：

(1)隨機產(chǎn)生初始種群。初始種群即父代種群，4個參數(shù)在各自的取值范圍內(nèi)構(gòu)建一個解空間?？刂苖atlab在該解空間中隨機產(chǎn)生初始種群，本研究中所采用的初始種群數(shù)量為10個，即產(chǎn)生一個10行4列的初始參數(shù)矩陣，矩陣中每一行代表一個個體，也就是代表一組參數(shù)取值。

(2)輸入監(jiān)測數(shù)據(jù)。將該降雨場次下監(jiān)測水深數(shù)據(jù)整理成可計算的數(shù)列變量輸入matlab之中。

(3)由初始種群生成SWMM input文件?？刂苖atlab將初始種群中的10組參數(shù)分別讀寫輸入input文件中，使程序中產(chǎn)生10個可用于調(diào)用計算的SWMM input文件。

(4)在matlab中調(diào)用SWMM模型動態(tài)鏈接庫(DLL)進行模型計算。由初始種群產(chǎn)生的10個input文件經(jīng)過調(diào)用計算后，分別在相應(yīng)的存儲路徑上產(chǎn)生10個report文件和10個output文件。其中report報告文件是對整個模型模擬過程的總結(jié)，以文本方式保存，output結(jié)果文件是在整個模擬過程中，各個元素在每個計算步長時的計算結(jié)果信息，例如節(jié)點水位、管道流量、子匯水區(qū)流量等，以二進制方式保存。

(5)讀取SWMM report文件將所需計算結(jié)果存入matlab之中。讀取模型計算后report文件中對應(yīng)計算結(jié)果值，提取與監(jiān)測數(shù)據(jù)時刻對應(yīng)的結(jié)果，將其整理成可計算的數(shù)列變量輸入matlab之中。

(6)計算目標函數(shù)值。目標函數(shù)是用于評價計算模擬值與監(jiān)測數(shù)據(jù)值吻合程度的函數(shù)。本文中選用Nash-Sutcliffe效率系數(shù)作為目標函數(shù)，可以用作考察模型計算的準確程度[9]。Nash-Sutcliffe效率系數(shù)的計算公式為：

(5)

(7)分配初始種群適應(yīng)度函數(shù)。根據(jù)種群中各個個體的目標函數(shù)值，對初始種群適應(yīng)度進行排序，目標函數(shù)越大則個體適應(yīng)度值也越大，并將排序結(jié)果寫入對應(yīng)的個體適應(yīng)度值的列向量。

(8)選擇交叉變異產(chǎn)生子代。由上一步所得的個體適應(yīng)度列向量作為選擇依據(jù)，適應(yīng)度值越大的個體被選中作為下一代的概率越大，相反，適應(yīng)度值小的個體被選中的概率會很小。調(diào)用遺傳算法工具箱中選擇函數(shù)，從初始種群中選擇優(yōu)良個體，并將選擇的個體寫入下一代種群之中。

調(diào)用遺傳算法工具箱函數(shù)recombin從產(chǎn)生的新種群中完成交叉，采用單點交叉xovsp函數(shù)，使用交叉概率Px=0.7執(zhí)行交叉，并將執(zhí)行交叉后的個體寫入新種群之中。

調(diào)用遺傳算法工具箱函數(shù)mutbga對該種群進行實值變異，使用變異概率Pm=0.05執(zhí)行變異，并將執(zhí)行變異后的個體寫入新種群之中。

(9)控制matlab將所產(chǎn)生的子代種群(參數(shù)矩陣)輸入SWMM input文件中并保存。由初始種群經(jīng)過選擇交叉變異后產(chǎn)生子代參數(shù)矩陣，控制matlab進行文本文件的讀寫，將對應(yīng)的參數(shù)值修改并寫入input文件中，將input文件保存并用于下一步的計算。

(10)循環(huán)運行4～9步直至滿足迭代滿足要求，本文中采用最大迭代次數(shù)500為迭代結(jié)束條件。

毎代種群最優(yōu)解和種群均值變化跟蹤圖如圖5所示。實線代表每個迭代種群中目標函數(shù)的最大值，虛線表示迭代次數(shù)中種群均值的變化過程。從圖5中可以看出，初始種群的目標函數(shù)平均值在0.5以下，在迭代計算前120次左右時，種群最優(yōu)解和種群均值變化曲線有著顯著上升，為了防止目標函數(shù)解落入局部最優(yōu)，迭代繼續(xù)進行，在隨后的迭代計算中，即使有重組變異的過程使得種群目標函數(shù)均值處于波動狀態(tài)，但種群最優(yōu)解處于穩(wěn)定狀態(tài)，可以認為遺傳算法在此次尋優(yōu)的過程之中找到了全局最優(yōu)解。第500次迭代的結(jié)果為： Max Rate為50.25， Min Rate為12.61， Manning-N為0.018 2，N-Imperv為0.008 7。

圖5 500次迭代后種群最優(yōu)解與種群均值變化跟蹤圖

此時目標函數(shù)最大值為0.726 1，也就是此時經(jīng)計算后的Nash-Sutcliffe效率系數(shù)為0.726 1。當(dāng)計算結(jié)果NS值越接近于1，表示模擬值與觀測值越接近，說明模型的計算準確度越高。參考相關(guān)的文獻[9]可知，當(dāng)計算的NS效率系數(shù)大于0.7時，近似的認為模擬結(jié)果與實測結(jié)果的吻合程度較高，模型誤差在計算誤差允許的范圍內(nèi)。因此可以認為本次模型參數(shù)率定過程與實測結(jié)果吻合程度較高，模型計算結(jié)果與實測數(shù)據(jù)對比如圖6所示。

圖6 參數(shù)率定后計算結(jié)果與實測數(shù)據(jù)對比圖

3 結(jié) 語

(1)本文該場次降雨下參數(shù)敏感性分析結(jié)果表明，霍頓滲透參數(shù)對輸出目標總徑流量敏感性最大，管道曼寧系數(shù)對輸出目標峰值流量和峰值時間敏感性最大；

(2)選取4個敏感參數(shù)進行率定，遺傳算法結(jié)果表明，Max Rate為50.25，Min Rate為12.61，Manning-N為0.018 2，N-Imperv為0.008 7，該場降雨下目標函數(shù)NS值為0.726 1。

(3)通過MATLAB調(diào)用SWMM進行水力計算來實現(xiàn)模型自動率定功能，充分利用SWMM模型的水力計算和MATLAB的數(shù)據(jù)分析功能，實現(xiàn)模型自動率定。該方案及編程可以應(yīng)用于實際工程的參數(shù)率定過程，并得出與實測結(jié)果最接近的參數(shù)組合，提高了模型模擬的精度及穩(wěn)定性。