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

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

(武漢大學，武漢 430000)

0 引 言

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

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

1 研究方法

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

1.1 研究區域概況

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

圖1 研究區域下墊面分析

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

圖2 研究區域模型概化示意圖

1.2 參數選擇及取值

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

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

表1 SWMM不可測量參數及其取值范圍

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

(1)

式中：Area為子匯水區面積。

Slope=KslopeSlope0

(2)

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

1.3 Morris全局敏感性分析法

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

圖3 Morris設計原理

(3)

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

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

(4)

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

1.4 遺傳算法

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

2 結果分析

2.1 Morris敏感性分析法計算結果

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

表2 各參數敏感性分析計算結果

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

2.2 遺傳算法步驟及計算結果

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

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

圖4 遺傳算法率定參數框架設計

算法具體計算步驟如下：

(1)隨機產生初始種群。初始種群即父代種群，4個參數在各自的取值范圍內構建一個解空間。控制matlab在該解空間中隨機產生初始種群，本研究中所采用的初始種群數量為10個，即產生一個10行4列的初始參數矩陣，矩陣中每一行代表一個個體，也就是代表一組參數取值。

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

(3)由初始種群生成SWMM input文件。控制matlab將初始種群中的10組參數分別讀寫輸入input文件中，使程序中產生10個可用于調用計算的SWMM input文件。

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

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

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

(5)

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

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

調用遺傳算法工具箱函數recombin從產生的新種群中完成交叉，采用單點交叉xovsp函數，使用交叉概率Px=0.7執行交叉，并將執行交叉后的個體寫入新種群之中。

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

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

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

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

圖5 500次迭代后種群最優解與種群均值變化跟蹤圖

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

圖6 參數率定后計算結果與實測數據對比圖

3 結 語

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

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

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