基于Morris與Sobol法的SWMM模型參數(shù)敏感性分析

劉鵬霄 馬立山 唐中楠 王伊林

(河北建筑工程學(xué)院，河北 張家口 075000)

0 引 言

近年來，我國經(jīng)濟發(fā)展迅速，城市化進程加快，城市的熱島效應(yīng)隨著城市化的發(fā)展而加劇，城市內(nèi)澇頻發(fā)[1].在此種情況下，水文模型的精確性為預(yù)測洪澇災(zāi)害，提早進行防范提供了可能[2].SWMM模型是由美國環(huán)保署研究的一個水文水質(zhì)模型，它能夠?qū)τ诮涤陱漠a(chǎn)生開始到形成徑流等一系列過程進行模擬，在國內(nèi)的眾多水文模型中得到了廣泛的應(yīng)用[3].模型主要由結(jié)構(gòu)和參數(shù)兩部分組成，參數(shù)又可以分為純物理參數(shù)、純經(jīng)驗參數(shù)和半經(jīng)驗參數(shù)三類，對于模型而言，模型的參數(shù)選擇對于模型的精確程度有很大影響[4].通過分析不同研究區(qū)的模型參數(shù)對于輸出結(jié)果的影響程度，降低模型參數(shù)的不確定性，從而提高模型率定的效率與精度[5].

參數(shù)敏感度分析方法可以分為局部敏感性分析和全局敏感性分析兩類[6].局部敏感性分析主要用于評價參數(shù)對于輸出結(jié)果的局部影響，而全局敏感度分析則可以考慮到參數(shù)整體對于輸出結(jié)果的影響，注重參數(shù)之間的作用對于輸出的綜合效應(yīng)[7].局部靈敏度分析方法主要包括Morris篩選法、EFAST法、擾動分析法等，全局分析法包括Sobol法(方差分解法)、Extend FAST法、RSA方法(區(qū)域靈敏度分析方法)、GLUE法(普適似然度法)、RSMSobol法、逐步回歸法、方差分析法等[8].史蓉等[9]利用Morris篩選法將SWMM的模型參數(shù)排序；張勝杰等[10]利用修正Morris篩選法對不同雨強下的參數(shù)進行敏感性分析，得出雨強與下滲參數(shù)關(guān)系；常曉棟等[11]利用Sobol法進行敏感度分析，認為此方法適應(yīng)于城市水文模型敏感性相關(guān)研究；張質(zhì)明等[12]將Sobol法應(yīng)用于WASP模型得出模型參數(shù)之間耦合性強結(jié)論.就目前情況看，在對于模型參數(shù)進行敏感度分析時，一般只運用某一種方法，考慮單個參數(shù)影響，而未能將參數(shù)組合影響考慮進去，造成分析結(jié)果的偏差.本文選取總產(chǎn)流和峰值流量兩相作為目標函數(shù)，利用Morris法和Sobol法兩種局部與全局相結(jié)合的敏感度分析方法，對于研究區(qū)的SWMM模型進行分析.其研究結(jié)果有利于今后模型敏感度分析方法的選取，進一步提升模型的精度.

1 研究材料及方法

1.1 研究區(qū)域概況

研究區(qū)亞熱帶季風氣候，日照時間長，歷史平均氣溫15.43℃，年平均降水1170mm，因為屬亞熱帶季風氣候，受東南沿海暖濕氣流影響，降雨本身存在明顯季節(jié)性變化差異，夏季降雨占全年總降雨接近50%，全年的5月至9月份是暴雨高發(fā)期，這一段時間易形成洪澇災(zāi)害，并且持續(xù)時間長，危害性大.

1.2 研究區(qū)模型構(gòu)建

研究區(qū)域面積為21.55km2.不透水區(qū)面積16.94km2，占總面積的78.60%.根據(jù)研究區(qū)域研究時間段的衛(wèi)星影像、數(shù)字高程數(shù)據(jù)、土地利用現(xiàn)狀和現(xiàn)狀排水管網(wǎng),利用GIS將研究區(qū)域面積、不透水比例、坡度、下墊面情況和管網(wǎng)信息輸入到SWMM模型中，構(gòu)建雨水徑流模型，將整個研究區(qū)域劃分為26個子匯水區(qū)、27個節(jié)點和17個排放口.子匯水區(qū)主要通過高程圖和坡度圖，與具體的主干道分布位置進行大的集水區(qū)的劃分，經(jīng)過人工調(diào)整后的子匯水區(qū)劃分，輸入降雨時間序列，添加雨量計，完成模型構(gòu)建，如圖1所示.

圖1 SWMM模型圖

1.3 模型參數(shù)與目標函數(shù)

1.3.1 模型參數(shù)

模型采用了Horton模型與非線性水庫模型來模擬降雨徑流過程，Horton模型需要參數(shù)主要包括：最大入滲速率、最小入滲速率、衰減速率常數(shù)、干燥時間、最大可能滲透深度[13].非線性水庫模型所需參數(shù)主要包括匯水面積、不透水百分比、寬度、坡度、不透水區(qū)曼寧系數(shù)、透水區(qū)曼寧系數(shù)、不透水區(qū)積水深度、不發(fā)生積水百分率[14].其中，面積、寬度、坡度、不透水區(qū)域所占百分數(shù)、透水區(qū)不發(fā)生積水所占百分數(shù)等根據(jù)地理信息系統(tǒng)獲取，排水管道的管網(wǎng)數(shù)據(jù)從排水資料和施工圖中獲取.根據(jù)現(xiàn)有資料的具體情況，現(xiàn)選擇敏感性分析的參數(shù)共8項，其具體的物理意義和取值范圍見下表：

表1 SWMM模型參數(shù)取值范圍

1.3.2 目標函數(shù)

根據(jù)模型所涉及到的主要產(chǎn)匯流參數(shù)，將能表征徑流特征的峰值流量和總產(chǎn)流兩項提出，并結(jié)合水文相關(guān)規(guī)范，最終確定目標函數(shù)選擇為峰值流量絕對百分比誤差和總產(chǎn)流的平衡誤差系數(shù)，其具體的計算方法如下表所示：

表2 目標函數(shù)

1.4 參數(shù)敏感性分析方法

1.4.1 Morris法

Morris篩選法是由Max D.Morris首次在1991年提出的，它作為局部敏感性分析方法關(guān)注單個參數(shù)對于模型模擬效果的影響，操作性強、簡單快捷、被廣泛應(yīng)用[15].Morris篩選法主要思想是改變單一變量參數(shù)xi而其他參數(shù)不變后，在參數(shù)允許的閾值內(nèi)隨機改變xi取值，通過其對于輸出值的影響來判斷參數(shù)敏感性的一種方法[16].本次研究以10%為固定步長對于各個水文水力參數(shù)進行擾動，取值分別為±10%、±20%、±30%.在保證某一參數(shù)變化時控制其他參數(shù)不變，探究單一參數(shù)對于目標值的影響，選取目標值為總產(chǎn)流和峰值流量兩項重要表征項相關(guān)形式，具體式如下：

(1)

其中：SN：敏感性判別因子；Yi：模型第i次運行輸出值；Yi+1：模型第i次運行輸出值；Y0：參數(shù)調(diào)整后的計算結(jié)果初始值；Pi：第i次模型運行參數(shù)值相對于校準后參數(shù)值的變化百分率；Pi+1：第i+1次模型運行參數(shù)值相對于校準后參數(shù)值的變化百分率.

1.4.2 Sobol法

基于蒙特卡羅法的Sobol法是一種定量全局敏感度分析方法，它能夠求解高度非線性模型中的參數(shù)間的相互作用所產(chǎn)生的靈敏度[17].它主要是利用了目標函數(shù)的方差可以分解為單個參數(shù)的方差和各個參數(shù)之間相互作用的方差的思想，來表征各個參數(shù)的獨立或相互關(guān)系靈敏度.Sobol法中，假設(shè)模型的函數(shù)形式為y=f(x)，其中參數(shù)x可以由xi(i=1,2，…，n)的n維離散表示.若f(x)本身可積，xi服從[0,1]上的均勻分布，則f(x)可以表示為：

(2)

式中共2n項相加和，將其寫作方差的形式，利用D表示，并且歸一化則可以得到如下公式：

(3)

式中各項即為一階、二階、總階靈敏度表達式：

(4)

(5)

(6)

式中：Si表征xi單獨作用時對于輸出項的靈敏度值，Sij表征xi和xj兩個參數(shù)共同作用時對于輸出項的靈敏度值，STi表征有xi參與的所有項對于輸出項的靈敏度值.

2 結(jié)果與分析

2.1 局部敏感度分析

利用修正Morris篩選法對研究區(qū)子匯水區(qū)以固定步長進行擾動，根據(jù)目標函數(shù)對各個參數(shù)的敏感度進行排序，結(jié)果表3所示：

表3 局部靈敏度結(jié)果排序

通過分析得知，對于峰值流量而言，敏感參數(shù)中敏感度大小依次為N-imperv、N-perv和Des-imperv，其中N-imperv為高敏感參數(shù).對于總產(chǎn)流而言，敏感參數(shù)則依照敏感度大小為N-imperv、Min.Infil.Rate和Des-imperv.對于此兩項目標函數(shù)，剩余參數(shù)均為不敏感參數(shù).

以峰值流量的絕對百分比誤差為目標函數(shù)時，N-imperv的敏感程度最高，這是因為研究區(qū)不透水區(qū)面積相對更廣，與不透水區(qū)徑流相關(guān)參數(shù)(N-imperv和Des-imperv)對于模型的產(chǎn)流過程取重要作用，相關(guān)的參數(shù)敏感程度高.對于總產(chǎn)流平衡誤差系數(shù)為目標函數(shù)時，由于下墊面的透水情況所造成的產(chǎn)流量的影響尤為明顯.本次采用的輸入降雨時間序列為20年一遇的實測點降雨數(shù)據(jù)，降雨強度較大的情況下，土壤中的含水量迅速增加，與下滲有關(guān)的參數(shù)敏感性相對較強.

兩目標函數(shù)下的敏感參數(shù)對比可知，N-imperv的敏感度排序未發(fā)生變化，但是其敏感度卻由高敏感參數(shù)變?yōu)槊舾袇?shù)，敏感度水平有所下降.Min.Infil.Rate和Des-imperv為敏感參數(shù)，其余參數(shù)項為不敏感參數(shù).Min.Infil.Rate由不敏感參數(shù)轉(zhuǎn)變?yōu)槊舾袇?shù)，敏感度有所上升.這表明同一參數(shù)在輸出結(jié)果不同時，它的敏感性也會表現(xiàn)出差異.

2.2 全局敏感度分析

本次研究的主要研究參數(shù)共8項，按照Sobol法步驟進行計算，利用拉丁超立方抽樣，抽取3000個基礎(chǔ)樣本，共2組，分別用矩陣A、B表示，將B矩陣的每一列替換為A矩陣中的每一列，組成新的矩陣，與原來A、B矩陣共同構(gòu)成(8+2)組樣本數(shù)矩陣.將得到的10組共30000個數(shù)據(jù)樣本帶入到SWMM模型中進行模擬，輸出共30000組狀態(tài)報告文件，按照要求提取相應(yīng)的輸出值，將其進行處理后作為目標函數(shù)，對其進行敏感度分析.

本文中主要求解模型的一階敏感度和總階敏感度，探究參數(shù)單獨作用時和與其他參數(shù)共同作用時對于目標函數(shù)的敏感度的影響.交互敏感度的計算方法為一階敏感度與總階敏感度差值，用以表示模型某一參數(shù)與其他參數(shù)的相互作用.

根據(jù)計算結(jié)果，將參數(shù)對應(yīng)于兩個目標函數(shù)的一階靈敏度(Si)、總階靈敏度(STi)與交互敏感度(STi-Si)匯總到表4和表5中.選用0.01作為敏感度的閾值，及敏感度值在0.01及以上，表示敏感參數(shù)，否則為敏感參數(shù).

表4 關(guān)于總產(chǎn)流的全局靈敏度結(jié)果分析

表5 關(guān)于峰值流量的全局靈敏度結(jié)果分析

對于總產(chǎn)流的一階敏感度而言，敏感參數(shù)排序由大至小依次為N-imperv、Min.Infil.Rate、Des-imperv、Des-perv，與局部敏感度分析相同，強降雨Rain1條件下，與下滲有關(guān)的參數(shù)Decay constant、Min.Infil.Rate為敏感參數(shù)，Max.Infil.Rate值同樣受影響，敏感程度上升.對于峰值流量的一階敏感度而言，敏感參數(shù)排序由大至小依次為N-imperv、N-perv、Des-perv、Min.Infil.Rate、Des-imperv、Decay constant.

對于總產(chǎn)流而言，一階敏感度中N-imperv的敏感性與其他參數(shù)相比較明顯要大，這說明N-imperv參數(shù)對于研究區(qū)域的總產(chǎn)流而言是十分重要的參數(shù)項，地表徑流受人類活動與下墊面影響較大.對于峰值流量的一階敏感度而言，除了N-imperv之外，N-perv的敏感度排序較高，這說明研究區(qū)的不透水區(qū)面積較大，但是透水區(qū)在對于峰值流量的影響仍然起了很大作用.

交互敏感度是總階敏感度與一階敏感度的差值，可以表征參數(shù)之間的耦合性.N-imperv、Min.Infil.Rate、N-perv、Des-perv的交互敏感度在10%左右，其值表示此獨立參數(shù)與其他參數(shù)的耦合性強.Des-perv和Decay constant的值相對較小，其表現(xiàn)出的與其他參數(shù)的耦合性更差一些，受到其他參數(shù)的影響較小，屬于相對獨立的參數(shù).

2.3 對比分析

表6 全局與局部靈敏度排序匯總

由此表可以看出，修正Morris篩選法作為局部敏感度分析法，對兩項目標函數(shù)的敏感度分析中，分別有三項參數(shù)被識別為敏感參數(shù)，其余為不敏感參數(shù).而Sobol法進行敏感度分析時，對于兩項目標函數(shù)的敏感參數(shù)項包括但是不僅限于原本的三項，而是有更多的敏感參數(shù)項被識別出來.由此可見，修正Morris篩選法對于參數(shù)的敏感性分析具有一定的借鑒意義，但是其對于模型總體敏感參數(shù)的識別并不突出，許多敏感參數(shù)被識別為不敏感參數(shù).在分析敏感度時，需要綜合考慮，不能作為直接依據(jù).

3 結(jié) 論

(1)根據(jù)研究區(qū)情況，利用局部Morris篩選法時，共篩選出N-imperv、N-perv、Des-imperv和Min.Infil.Rate 4項敏感參數(shù)，這些與地表徑流和下滲相關(guān)參數(shù)表現(xiàn)出的敏感性與所選目標函數(shù)相關(guān).利用全局Sobol法時，除去以上4項，同時將Des-perv和Decay constant判別為敏感參數(shù).兩種方法產(chǎn)生差異的原因主要是由于前者關(guān)注單個參數(shù)的影響，而后者則關(guān)注參數(shù)對整體的影響.Des-perv和Decay constant單獨作用時，對于研究區(qū)的兩項目標函數(shù)的影響較小，但在與其他參數(shù)共同作用時，對其他參數(shù)產(chǎn)生影響較大，間接影響最終結(jié)果，被標識為敏感參數(shù).

(2)局部修正Morris篩選法對于對于敏感參數(shù)項有一定的識別作用，全局Sobol法所需基礎(chǔ)參數(shù)大，通過對于參數(shù)進行敏感性分析，認為局部敏感性分析方法具有一定的參考性，但是卻會對于一些在全局敏感參數(shù)項識別有誤，因此在進行敏感性分析的時候，應(yīng)當全局與局部敏感分析綜合考慮.同時，敏感性分析時并未考慮降雨強度這一因素對于參數(shù)敏感性的影響，在今后的研究中，應(yīng)當加強此方面分析，為模型優(yōu)化提供更可行精確的研究方法.