基于偏秩相關(guān)-逐步回歸法的SWMM模型全局敏感性分析

李傳奇，崔佳偉，孫 策，段明印，馬夢蝶，楊幸子

(山東大學(xué)土建與水利學(xué)院，山東 濟南 250061)

0 引 言

近幾年，隨著我國經(jīng)濟社會的發(fā)展、氣候的變化以及城市化和工業(yè)化的快速推進，導(dǎo)致城市'熱島效應(yīng)'不斷加重和下墊面條件發(fā)生變化，造成城市洪澇災(zāi)害頻繁發(fā)生，人民的生命財產(chǎn)安全受到了嚴(yán)重的威脅[1]。利用水文模型可以有效地模擬城市洪澇過程，對城市洪水進行預(yù)報預(yù)警和管理，降低甚至避免城市洪水帶來的危害[2]。SWMM模型是由美國環(huán)保署研究開發(fā)的一種暴雨洪水管理模型，能夠動態(tài)的模擬研究區(qū)域內(nèi)的各種水文過程，包括降雨-徑流、地下水與排水管道的水量交換、污染物運移和地面積水等，在國內(nèi)用來管理城市洪水已經(jīng)得到了廣泛的認可和使用[3]。

與大多數(shù)的水文模型類似，SWMM模型由于其眾多的相關(guān)參數(shù)，在進行模型率定時工作量大、費時費力，而且不能保證每一個參數(shù)的精度[4]。這時我們就需要對SWMM模型的相關(guān)參數(shù)進行敏感性分析，得到各個參數(shù)的敏感性即對于模型輸出的影響程度，對于那些敏感性小的參數(shù)，在進行參數(shù)率定時可以取其經(jīng)驗值進行固定，高效快速地完成參數(shù)率定，提高模型的可靠性和穩(wěn)定性[5]。一般敏感性分析方法有兩種，分別是局部敏感性分析和全局敏感性分析[6]。局部敏感性分析是固定其他參數(shù)不變，在取值范圍內(nèi)對研究參數(shù)選取一系列的值帶入模型得到一系列的輸出再進行分析，這種方法沒能考慮到參數(shù)之間的交互作用對模型輸出的貢獻；全局敏感性分析通過全部參數(shù)在取值范圍內(nèi)的隨機抽樣形成一系列的輸入數(shù)組帶入模型得到一系列的輸出再進行分析，這種方法考慮了參數(shù)之間的聯(lián)合效應(yīng)對輸出的影響，結(jié)果穩(wěn)定可靠，能夠更加全方位的分析參數(shù)敏感性[7]。目前，常用的敏感性分析方法有基于修正的Morris篩選法[8]、互信息法[9]、偏相關(guān)法[10]、基于Sobol的指數(shù)法[11]、逐步回歸法[12]和偏秩相關(guān)法[13]等。其中，逐步回歸法和偏秩相關(guān)法都是全局敏感性分析方法，對參數(shù)的概率分布無特殊要求，適合SWMM這種參數(shù)眾多而且非線性的模型。

本文以山東大學(xué)千佛山校區(qū)為例，運用SWMM模型分析降雨時該校區(qū)的水文過程，針對該模型運行時的輸入?yún)?shù)，采用拉丁超立方對其在取值范圍內(nèi)進行隨機抽樣，利用逐步回歸法和偏秩相關(guān)法對參數(shù)的敏感性進行分析和量化，為模型進一步的高效率定和不確定性分析提高參考和指導(dǎo)，進一步提高模型的精度和可靠性。

1 材料與方法

1.1 SWMM模型及參數(shù)

SWMM模型自從被開發(fā)以來就被全世界的科學(xué)研究人員用來研究模擬雨洪，主要包括水文模塊和水質(zhì)模塊，應(yīng)用范圍廣泛，特別是城市區(qū)域內(nèi)降雨條件下的各種水文過程，管渠以及下水道中的水流流動等。在對山東大學(xué)千佛山校區(qū)進行模擬時，主要水文過程有管道輸、排水和地表徑流兩個方面，求解方法為非線性水庫模型和圣維南方程。

SWMM模型的相關(guān)水力水文參數(shù)共計14個，這其中管道的長度和不透水區(qū)域面積比率一般為研究區(qū)域的固定值，目前可以通過GIS、RS、DEM或者直接測量得到固定數(shù)值且誤差不大，對其值進行大幅度的變化毫無意義，則只對其余12個參數(shù)進行敏感性分析。其中匯水區(qū)的坡度、面積和寬度3個參數(shù)由于其空間特性和測量時存在不可避免的誤差，定義3個修正因子Pct-Area、K-With、和K-Slope來概化這3個參數(shù)[14]。通過查找用戶手冊和相關(guān)文獻以及研究區(qū)域的實地情形，得出各個研究參數(shù)的取值范圍見表1。

表1 SWMM模型各研究參數(shù)取值范圍及意義

1.2 逐步回歸法

逐步回歸法是一種利用最小二乘法建立因變量與自變量的多元線性回歸方程，再求解回歸系數(shù)的方法，假設(shè)因變量為Y,則它與一系列自變量X1,X2,…,Xn之間的回歸方程可以表述如下：

y=a0+a1X1+a2x2+…+anXn+ε

(1)

式中：y為因變量Y的回歸值；a0,a1,a2,…,an為各個自變量所對應(yīng)的回歸系數(shù)；ε為誤差項，表示除了自變量以外的因素對因變量Y回歸值的影響。

在式(1)中，不是所有的自變量X1,X2,…,Xn都對因變量 有顯著的貢獻，如果能將影響程度不明顯的自變量從回歸方程中去掉，不但能減少計算量還能提高回歸方程的預(yù)測水平，這時就可以利用逐步回歸法，原理是：將自變量X1,X2,…,Xn一個一個逐步引入到回歸方程中，具體的做法是該變量通過了給定偏 統(tǒng)計量顯著性檢驗的情況下引入方程，此時便計算方程中所有存在的自變量的偏回歸平方和，并進行排序，對偏回歸平方和最下的自變量在進行顯著性檢驗，如果不顯著則從回歸方程中去除，這樣直到?jīng)]有新的自變量引入和老的自變量去除，得到一個新的回歸方程。此時每一次回歸分析的R2代表引入的該自變量對于回歸方程的解釋程度，標(biāo)準(zhǔn)回歸系數(shù)SRC的絕對值大小表示參數(shù)的敏感性程度，正、負符號表示該自變量與因變量的正、負相關(guān)關(guān)系[15]。

1.3 偏秩相關(guān)法

偏秩相關(guān)法是一種利用'等級位差'來進行分析的全局敏感性分析方法，適用條件為各參數(shù)有相同維度的隨機變化。首先要計算各參數(shù)的秩相關(guān)系數(shù)，通過抽樣方法生成 維自變量 ，帶入模型得到 維輸出變量 ,則組成的輸入、輸出變量矩陣如(2)所示。

(2)

變量Y,X1,X2,X3,…,Xn中，任意兩個變量之間的秩相關(guān)系數(shù)為:

(3)

式中：Ri為n維變量Y/Xk中第i個變量的秩(該變量在 維變量從大到小排序后所在的位次)；Qi為n維變量Xm中第i個變量的秩。則得到這些輸入、輸出變量間的秩相關(guān)系數(shù)矩陣P和P的逆矩陣C為：

(4)

(5)

式中：偏秩相關(guān)系數(shù)|PXi|的值越大，說明該參數(shù)的敏感性越大即對模型輸出的影響程度，PXi的符號代表該參數(shù)與輸出變量之間的正、負相關(guān)關(guān)系[16]。

1.4 抽樣方法

1.5 研究區(qū)域概況

山東大學(xué)千佛山校區(qū)南苑位于山東省濟南市歷下區(qū)，總面積為79 887.53 m2，是季風(fēng)氣候明顯區(qū)域之一。寒冷晴朗，雨雪稀少，多偏北風(fēng)。夏季受熱帶、副熱帶海洋氣團影響，盛行來自海洋的暖濕氣流，天氣炎熱，雨量充沛，光照充足，多偏南風(fēng)。春季和秋季是冬季轉(zhuǎn)夏季、夏季轉(zhuǎn)冬季的過渡季節(jié)，風(fēng)向多變；年平均降水量672.7 mm，其中7、8月份最多，占全年的55%以上，三面環(huán)山的地理情況讓熱空氣無法擴散，使得區(qū)域內(nèi)夏季高溫濕熱。該研究區(qū)域由于特殊的地形關(guān)系，無其他徑流流入，是一片相對獨立的區(qū)域。其中，不透水區(qū)域包括混凝土地面、屋頂和運動場地(混凝土和塑膠)，透水區(qū)域包括花壇、草坪等綠地，土地利用情況見圖1。

圖1 千佛山校區(qū)南院土地利用情況

對研究區(qū)域進行SWMM模型數(shù)字概化，其中子匯水區(qū)域24個，節(jié)點23個，9條為地下管線，14條為道路，3個雨水排放口。該研究區(qū)域的概化圖見圖2。可以看出，該研究區(qū)域四季分明，降雨量充沛且土地利用類型覆蓋完整、子匯水區(qū)域復(fù)雜，對于季風(fēng)氣候下區(qū)域的SWMM雨洪模型具有一定的代表性。

圖2 研究區(qū)域概化圖

1.6 設(shè)計降雨

本研究的設(shè)計降雨計算基于最新的濟南市降雨強度計算公式，我國降雨中雙峰雨型和多峰雨型較少，大多為單峰雨型，并且大部分地區(qū)的雨峰系數(shù)在0.3～0.5之間，這種情況與芝加哥雨型非常相似。由此本研究利用芝加哥降雨模型推求降雨過程線具有一定的代表性，得到的設(shè)計降雨重現(xiàn)期為1 a，降雨量52.76 mm，降雨強度時間間隔為2 min，雨峰系數(shù)為0.4。設(shè)計降雨過程線如圖2所示，濟南降雨強度計算公式如(6)。

(6)

式中：q表示降雨強度即每分鐘的降雨量，mm/min；P表示降雨重現(xiàn)期，a；t表示降雨歷時，min。

圖3 研究區(qū)域設(shè)計降雨過程線

2 結(jié)果與分析

利用matlab程序?qū)崿F(xiàn)拉丁超立方抽樣，需要將SWMM模型的12個研究參數(shù)取值范圍的上限和下限分別帶入程序計算得到1 000組輸入變量組合，再通過matlab的for循環(huán)語句寫入SWMM開源的輸入文件，最后將1 000個輸入文件依次代入模型進行計算得到3個1 000組重要的水文輸出變量：總產(chǎn)流量、峰值流量和峰值時間。之后，再將輸入、輸出變量組成3個 的矩陣，應(yīng)用逐步回歸法和偏秩相關(guān)法對其進行全局敏感性分析。

2.1 逐步回歸法分析

利用matlab實現(xiàn)逐步回歸法的基本程序，之后對輸入、輸出的變量矩陣的進行敏感性分析，結(jié)果如表2所示。對于峰值流量來說，N-Imperv和k-Width的|SRC|值相差無幾，稍大于其他參數(shù)，且引入這兩個參數(shù)時，R2的增量也相差無幾，但整體上來看，N-Imperv的敏感性要稍大于K-Width；在引入所有參數(shù)后，R2的值達到了0.929，說明峰值流量與各個參數(shù)之間的線性相關(guān)關(guān)系很好。對于峰值時間來說，N-Imperv的|SRC|值和引入方程時R2的增量都大于其他的參數(shù)，是對峰值時間最敏感、起決定性作用的參數(shù)，并且是正相關(guān)關(guān)系；在引入所有參數(shù)后，R2的值達到了0.884，說明峰值時間與各個參數(shù)之間的線性相關(guān)關(guān)系很好，其中，剔除了S-Perv、Max.Infil.Rate和Min.Infil.Rate這3個不敏感參數(shù)。對于總產(chǎn)流來說，Per-Area的|SRC|值和引入方程時R2的增量都大于其他的參數(shù)，是對總產(chǎn)流最敏感、起決定性作用的參數(shù)，并且是正相關(guān)關(guān)系；在引入所有參數(shù)后，R2的值達到了0.952，說明總產(chǎn)流與各個參數(shù)之間的線性相關(guān)關(guān)系很好，其中，剔除了N-Imperv和Manning-N這兩個不敏感參數(shù)。

表2 逐步回歸法分析結(jié)果

同時，不透水區(qū)的曼寧系數(shù)N-Imperv對于排放口的峰值流量和時間都是最敏感的參數(shù)、但是對于總產(chǎn)流來說卻是不敏感的參數(shù)，分析可知，N-Imperv只會影響水流在不透水地面上的流動狀態(tài)，并不會影響水流的總量。這個結(jié)果也間接說明了逐步回歸法分析的準(zhǔn)確性。

2.2 偏秩相關(guān)法分析

利用matlab程序?qū)崿F(xiàn)偏秩相關(guān)法的基本程序，之后將1 000 組SWMM模型的輸入、輸出數(shù)據(jù)帶入程序進行敏感性分析，得到如表3所示的分析結(jié)果。對于峰值流量來說，N-Imperv是最敏感的參數(shù)，偏秩相關(guān)系數(shù)P的值達到了-0.864，與峰值流量有較高的負線性相關(guān)關(guān)系；K-Width是第二敏感的參數(shù)，P值為0.853，與峰值流量有較高的正線性相關(guān)關(guān)系；是因為這兩個參數(shù)會影響水流在研究區(qū)域內(nèi)的蓄、流狀態(tài)，從而影響了峰值流量。對于峰值時間來說，N-Imperv和K-Width是敏感性排名前兩位的參數(shù)，P值分別達到了0.926和-0.822，與峰值時間有較高的正線性相關(guān)和負線性相關(guān)關(guān)系，同時也都通過影響水流在研究區(qū)域內(nèi)的蓄、流狀態(tài)，從而影響了峰值時間。對于總產(chǎn)流來說，最敏感的參數(shù)為Pct-Area，P值達到了0.919，與總產(chǎn)流有較高的正線性相關(guān)關(guān)系，這是因為，子匯水區(qū)面積比例因子的增大會增加計算降雨量，從而增大了總產(chǎn)流量。

表3 偏秩相關(guān)法分析結(jié)果

2.3 結(jié)果對比分析

將兩種敏感性分析方法結(jié)果畫在同一個柱狀圖里進行對比，見圖4、5、6。可以發(fā)現(xiàn)，對于峰值流量來說，兩種方法得到的最敏感的兩個參數(shù)都是N-Imperv和K-Width，并且所有參數(shù)的敏感性排序和正、負相關(guān)關(guān)系都相同；對于峰值時間來說，兩種方法得到的最敏感參數(shù)都是N-Imperv，逐步回歸法剔除了S-Perv、Max.Infil.Rate和Min.Infil.Rate 3個參數(shù)，在偏秩相關(guān)法里，這3個參數(shù)P值很小也是最不敏感的參數(shù)，并且其他參數(shù)的敏感性排序和正、負相關(guān)關(guān)系都相同；對于總產(chǎn)流來說，兩種方法得到的最敏感參數(shù)都是Pct-Area，逐步回歸法剔除了N-Imperv和Manning-N這兩個參數(shù)，在偏秩相關(guān)法里，這2個參數(shù)P值很小也是最不敏感的參數(shù)，并且其他參數(shù)的敏感性排序和正、負相關(guān)關(guān)系都相同。

圖4 峰值流量分析結(jié)果對比

圖5 峰值時間分析結(jié)果對比

圖6 總產(chǎn)流分析結(jié)果對比

3 結(jié) 語

(1)研究參數(shù)的輸入組合達到1000組，拉丁超立方的抽樣方法也保證了隨機抽樣的數(shù)據(jù)能均勻分布在取值范圍內(nèi)，數(shù)據(jù)具有一定的代表性，且逐步回歸法和偏秩相關(guān)法都是線性分析方法，兩種方法得出的結(jié)論完全相同，更說明了結(jié)論的合理性和可靠性以及SWMM模型建立的準(zhǔn)確性。

(2)通過兩種分析方法得知：N-Imperv是對峰值流量敏感性最大的參數(shù)，且它們之間是負相關(guān)關(guān)系，兩種方法得到的各個參數(shù)的敏感性排序相同；N-Imperv是對峰值時間敏感性最大的參數(shù)，隨著N-Imperv的增大峰值時間也會越來越長，兩種方法得到的各個參數(shù)的敏感性排序相同，并且敏感性最小的3個參數(shù)都是S-Perv、Max.Infil.Rate和Min.Infil.Rate，在進行模型率定時可以取固定值；Pct-Area是對總產(chǎn)流最為敏感的參數(shù)，與總產(chǎn)流呈正相關(guān)關(guān)系，兩種方法得到的各個參數(shù)的敏感性排序相同，并且敏感性最小的兩個參數(shù)都是N-Imperv和Manning-N。

(3)通過敏感性分析，精確識別了SWMM模型的敏感因素，為接下來不確定性分析和模型率定提供了 參考和指導(dǎo)，減少建模的工作量，提高模型率定的效率。