段明印，李傳奇，韓典乘，楊幸子

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

0 引 言

城市化的進(jìn)程使得我國城市雨洪災(zāi)害問題日益嚴(yán)峻，因此我國在防治雨洪災(zāi)害的同時(shí)正在大力推行海綿城市的建設(shè)，暴雨管理模型(SWMM)可以模擬城市降雨地表徑流的產(chǎn)生和在排水管網(wǎng)中的輸送，能夠用來指導(dǎo)海綿城市的建設(shè)，因此在國內(nèi)得到了較為廣泛的運(yùn)用[1]。

準(zhǔn)確率定SWMM的參數(shù)是精確構(gòu)建SWMM模型、指導(dǎo)海綿城市建設(shè)的前提。但是SWMM參數(shù)較多，在率定時(shí)為提高效率可僅率定對輸出變量影響較為顯著的參數(shù)，那些對輸出變量影響較小的參數(shù)可直接取經(jīng)驗(yàn)值。因此，需要對SWMM參數(shù)進(jìn)行敏感性分析，識別出影響較大的參數(shù)[2]。參數(shù)敏感性分析方法分為局部敏感性分析方法和全局敏感性分析方法，局部分析方法只是分析單個(gè)參數(shù)對輸出變量的影響，因其計(jì)算速度快，得到了較廣泛的運(yùn)用，但是忽略了參數(shù)之間相互作用對輸出變量的影響，在“異參等效”的情況下會使一些敏感參數(shù)無法識別；全局分析方法可以分析參數(shù)之間相互作用對輸出變量的影響，能夠在“異參等效”的情況下識別出敏感參數(shù)，適用于參數(shù)眾多的模型[3]。

常見的全局敏感性分析方法有多元逐步回歸法[4]、Sobol法[5]、FAST法[6]、偏相關(guān)法、偏秩相關(guān)法等[3]。這些方法都是通過分析參數(shù)和輸出變量的線性相關(guān)關(guān)系來識別敏感參數(shù)，忽略了參數(shù)和輸出變量的非線性相關(guān)關(guān)系，可能導(dǎo)致某些敏感參數(shù)無法識別。因此，本文以某小區(qū)為例，選取能夠分析變量之間非線性關(guān)系的互信息法分析參數(shù)全局敏感性，將以互信息法分析的結(jié)果與能夠分析變量之間線性關(guān)系的偏秩相關(guān)法分析的結(jié)果相結(jié)合，精確識別所有敏感參數(shù)，為下一步高效率定打下基礎(chǔ)。

1 材料與方法

1.1 SWMM模型參數(shù)

在利用SWMM模型對興隆山校區(qū)進(jìn)行降雨徑流模擬時(shí)，需要輸入的參數(shù)共計(jì)14個(gè)。這14個(gè)參數(shù)中的管渠長度和不透水面積可由設(shè)計(jì)資料讀取或者直接測量，所以僅對余下的12個(gè)參數(shù)的敏感性進(jìn)行分析。其中子匯水區(qū)的面積、寬度、坡度這3個(gè)輸入?yún)?shù)具有顯著空間特性，輸入時(shí)存在測量誤差，所以定義3個(gè)修正系數(shù)：K-Area(面積修正系數(shù))、K-With(寬度修正系數(shù))和K-Slope(坡度修正系數(shù))，確定這3個(gè)修正系數(shù)的分布范圍，則SWMM模型的輸入值為修正系數(shù)和對應(yīng)參數(shù)的乘積[7]。依據(jù)SWMM使用手冊和國內(nèi)相關(guān)文獻(xiàn)[7-11]確定參數(shù)的范圍，依據(jù)研究區(qū)域內(nèi)的實(shí)測值確定參數(shù)初始值，參數(shù)分布及取值如表1所列。

表1 SWMM參數(shù)取值表Tab.1 The parameters table of the storm water management model

1.2 SWMM模型構(gòu)建

以山東大學(xué)興隆山校區(qū)為研究區(qū)域，該區(qū)總面積為63.82 萬m2，其中綠地面積較大，綠地面積高達(dá)30.82 萬m2，建筑占地面積為9.2 萬m2，道路占地面積8.9 萬m2，塊石鋪裝占地面積1.9 萬m2，硬質(zhì)鋪裝面積8.7 萬m2，水體面積2.8 萬m2，足球場面積1.5 萬m2，研究區(qū)域概況如圖1所示。

圖1 研究區(qū)域概況Fig.1 The overview of study area

結(jié)合興隆山校區(qū)的地勢條件、地下管網(wǎng)資料、下墊面類型以及校區(qū)內(nèi)的河道等資料，構(gòu)建研究區(qū)的SWMM模型。SWMM模型包括45個(gè)子匯水區(qū)域，55個(gè)管道鉸點(diǎn)，55條雨水管道，3個(gè)排水口，興隆山SWMM模型如圖2所示。

圖2 研究區(qū)域SWMM模型Fig.2 Storm water management model of the study area

1.3 抽樣方法

常見的抽樣方法有簡單隨機(jī)抽樣、分層抽樣和整群抽樣，分層抽樣方法的原理是將樣本空間分成互不交叉的若干層，然后從各層中抽取樣本，相對于簡單隨機(jī)抽樣和整群抽樣，分層抽樣能夠使得樣本點(diǎn)均勻地分布在樣本空間內(nèi)，避免樣本集中的問題。因此本文選用拉丁超立方抽樣法在SWMM模型參數(shù)的分布空間內(nèi)進(jìn)行抽樣，拉丁超立方抽樣法屬于多維分層抽樣方法，能夠高效的抽取分布較為均勻的樣本[8]。

首先通過拉丁超立方抽樣法對有關(guān)降雨徑流的12個(gè)SWMM模型參數(shù)隨機(jī)抽樣1 000次，得到1 000組12個(gè)SWMM模型參數(shù)，然后通過MATLAB編程將這1 000組參數(shù)導(dǎo)入到SWMM模型的輸入文件中，生成1 000組新的SWMM輸入文件，利用SWMM模型將這1 000組輸入文件依次模擬，得到這1 000組參數(shù)對應(yīng)的SWMM輸出文件，提取這1 000組輸出文件中某排放口的3個(gè)重要水文變量：流量峰值、流量峰值發(fā)生時(shí)間、總產(chǎn)流量，因?yàn)檠芯繀^(qū)域主要通過1號排放口排放總產(chǎn)流，因此提取1號排放口的3個(gè)水文變量作為3個(gè)因變量，以12個(gè)輸入?yún)?shù)作為自變量，分別以偏秩相關(guān)法和互信息法分析自變量和因變量之間的關(guān)系。

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

偏秩相關(guān)分析是一種對變量分布不做要求的非參數(shù)統(tǒng)計(jì)分析方法，能夠?qū)Φ燃墧?shù)據(jù)進(jìn)行線性相關(guān)分析，在分析兩變量的相關(guān)關(guān)系時(shí)能夠控制其他變量的影響，因此能夠精確識別變量間的線性相關(guān)關(guān)系，可用做分析參數(shù)的全局敏感性[9]。

令輸出結(jié)果為Y，輸入?yún)?shù)為X1，X2，…，Xn，設(shè)這n+1個(gè)變量之間的Spearman秩相關(guān)系數(shù)r組成的矩陣為T，求矩陣T的逆矩陣得矩陣C。

(2)

則得某輸入?yún)?shù)Xi與輸出結(jié)果Y之間的偏秩相關(guān)系數(shù)PXi為：

(3)

式中：r為兩變量間Spearman秩相關(guān)系數(shù)；R2Xi為X1，X2,…,Xi-1,Xi+1,…,Xn線性回歸Xi的可決系數(shù)；RY2為它可決系數(shù)；Bi為X1，X2，…，Xn線性回歸Y的標(biāo)準(zhǔn)回歸系數(shù)；PXi為偏秩相關(guān)系數(shù)，它的絕對值大小為參數(shù)敏感性大小，它的正負(fù)表示呈正相關(guān)還是負(fù)相關(guān)。

1.5 互信息分析

互信息是指兩個(gè)變量之間或者多個(gè)變量之間共享的信息量，能夠度量變量之間的相關(guān)性，相關(guān)性越強(qiáng)，互信息越大。基于互信息的R統(tǒng)計(jì)量可以用來分析參數(shù)的全局敏感性，R統(tǒng)計(jì)量能夠精確識別參數(shù)和輸出結(jié)果之間的非線性相關(guān)關(guān)系，R統(tǒng)計(jì)量的大小代表參數(shù)的敏感性大小，但R統(tǒng)計(jì)量只能識別參數(shù)敏感性的大小，無法識別參數(shù)和輸出結(jié)果之間是正相關(guān)還是負(fù)相關(guān)[10]。

兩離散型隨機(jī)變量X與Y之間的互信息定義為：

(5)

R統(tǒng)計(jì)量表示為：

(6)

式中：p(x)為離散型隨機(jī)變量X取值x的概率；p(y)為離散型隨機(jī)變量Y取值y的概率；p(x,y)為X取值x、Y取值y同時(shí)發(fā)生的概率；I(x,y)為變量X、Y的互信息；b為互信息的量綱。本文互信息的量綱采用自然常數(shù)e，將SWMM模型的輸入?yún)?shù)視為離散型隨機(jī)變量X，將SWMM模型輸出結(jié)果視為離散型隨機(jī)變量Y，求解變量X與Y之間的互信息以及R統(tǒng)計(jì)量，分析參數(shù)對輸出結(jié)果影響的大小。在計(jì)算時(shí)為提高計(jì)算效率，將變量X分為i等份，將變量Y分為i等份，用pi表示變量X第i份中的數(shù)量占變量X總數(shù)量的概率，同時(shí)用pj表示變量Y第j份中的數(shù)量占變量Y總數(shù)量的概率，將分成i等份的變量X與分成j等份的Y組成聯(lián)立表，求得變量X和變量Y的聯(lián)合概率，用pij表示變量X在第i份取值、自變量Y在第j份取值同時(shí)發(fā)生的概率，則互信息的公式變?yōu)椋?/p>

(7)

式中：i、j為輸入?yún)?shù)和輸出結(jié)果分成的等份數(shù)，本文都分為10等份；

1.6 降雨資料

降雨資料采用興隆山校區(qū)內(nèi)雨量站(測站編碼：41822450)測得的降雨數(shù)據(jù)，選取20120708場次的降雨，該場次降雨歷時(shí)為13個(gè)小時(shí)，降雨量為47.5 mm，雨峰系數(shù)約為0.4，降雨強(qiáng)度過程如圖3所示。

2 結(jié)果與分析

2.1 初步分析

為初步分析SWMM模型的12參數(shù)和3個(gè)水文輸出結(jié)果的關(guān)系，分別畫出這1000組12個(gè)參數(shù)和3個(gè)輸出結(jié)果的殘差散點(diǎn)圖，其中散點(diǎn)圖的橫坐標(biāo)為某一參數(shù)與控制參數(shù)回歸分析得到的殘差，縱坐標(biāo)為輸出結(jié)果與控制參數(shù)回歸分析得到的殘差。各個(gè)參數(shù)和峰值流量的關(guān)系如圖4所示，觀察發(fā)現(xiàn)Manning-N與峰值流量的殘差散點(diǎn)圖幾乎在一條一直線上，具有高度的線性關(guān)系，呈負(fù)相關(guān)。各個(gè)參數(shù)和峰值時(shí)間的關(guān)系如圖5所示，觀察發(fā)現(xiàn)N-perv和峰值時(shí)間具有較強(qiáng)的線性關(guān)系，呈負(fù) 相關(guān)。各個(gè)參數(shù)和總產(chǎn)流量的關(guān)系如圖6所示，觀察發(fā)現(xiàn)Manning-N、Min-Rate和總徑流量具有較強(qiáng)的負(fù)相關(guān)關(guān)系，其中Manning-N和總產(chǎn)流量的負(fù)相關(guān)性強(qiáng)于Min-Rate。

圖3 研究區(qū)域20120708場次降雨強(qiáng)度過程Fig.3 Rainfall intensity processes of label 20120708 in the study area

圖4 參數(shù)和峰值流量殘差散點(diǎn)圖Fig.4 Residual scatter diagram of parameters and peak flow

圖5 參數(shù)和峰值時(shí)間殘差散點(diǎn)圖Fig.5 Residual scatter diagram of parameters and peak time

圖6 參數(shù)和總產(chǎn)流量殘差散點(diǎn)圖Fig.6 Residual scatter diagram of parameters and total runoff volume

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

通過MATLAB編程利用偏秩相關(guān)分析方法分析參數(shù)和輸出結(jié)果的相關(guān)性，分析結(jié)果如表2所列。由表2可以得出參數(shù)和輸出結(jié)果之間線性相關(guān)關(guān)系的強(qiáng)弱，即參數(shù)的敏感性大小。對峰值流量來說，最敏感的參數(shù)是Manning-N，它與峰值流量的偏秩相關(guān)系數(shù)P高達(dá)-0.996，具有高度的線性關(guān)系；第二敏感參數(shù)是N-perv，P值的絕對值位0.670，排名前兩位的敏感參數(shù)和峰值流量都是負(fù)相關(guān)關(guān)系，這是因?yàn)楣芮屯杆乇淼牟诼氏禂?shù)增大會使滯蓄雨水量增大，從而使得峰值流量減?。黄渌麉?shù)的敏感性相對較小，說明Manning-N和N-perv對峰值流量起著決定性作用。對峰值時(shí)間來說，最敏感的參數(shù)是N-perv，P值為-0.417，呈負(fù)相關(guān)，其他參數(shù)的敏感性都較小。對總產(chǎn)流量來說，最敏感的參數(shù)是Manning-N，P值高達(dá)-0.936，具有高度的線性關(guān)系，第二敏感參數(shù)是Min-Rate，P值為-0.838，其與總產(chǎn)流量也具有較強(qiáng)的線性關(guān)系；Manning-N、Min-Rate與總產(chǎn)流量都是負(fù)相關(guān)關(guān)系，Manning-N和Min-Rate的增大將會使總產(chǎn)流量減小，這是因?yàn)樽钚B透率增加將會使雨水下滲量增大，管渠糙率增大會使管渠滯蓄雨水量增大，這些都會使總產(chǎn)流量減小。

表2 偏秩相關(guān)分析結(jié)果Tab.2 Result of partial rank correlation analysis

2.3 互信息分析

通過MATLAB編程利用互信息相關(guān)原理分析參數(shù)和輸出結(jié)果之間的非線性相關(guān)關(guān)系，結(jié)果如表3所列。根據(jù)表3中的R統(tǒng)計(jì)值可以判斷參數(shù)和輸出結(jié)果之間非線性關(guān)系的強(qiáng)弱，及參數(shù)敏感性的大小。對峰值流量來說，最敏感的參數(shù)是Manning-N，R值高達(dá)0.886，這說明Manning-N和峰值流量之間具有較強(qiáng)的非線性相關(guān)關(guān)系；其他參數(shù)的R值較小，說明其他參數(shù)和峰值流量的非線性關(guān)系較弱。對峰值時(shí)間來說，N-perv是最敏感的參數(shù)，R值為0.457，說明N-perv和峰值時(shí)間的非線性關(guān)系較弱；第二敏感參數(shù)是K-Width，R值僅為0.337，對峰值時(shí)間影響較小，其他參數(shù)對峰值時(shí)間的影響也都較小。對總產(chǎn)流量來說，Manning-N是最敏感的參數(shù)，R值為0.574；第二敏感參數(shù)是Min-Rate，R值僅為0.372，說明其對總產(chǎn)流量的影響較小。

將偏秩相關(guān)分析的結(jié)果和互信息分析的結(jié)果進(jìn)行對比發(fā)現(xiàn)：兩種方法所求的對峰值流量最敏感的參數(shù)都是Manning-N，但是用偏秩相關(guān)法算出的敏感性大小為0.997，大于用互信息法算出的敏感性大小，說明Manning-N和峰值流量的線性關(guān)系強(qiáng)于非線性關(guān)系；對峰值時(shí)間而言，兩種方法求得的最敏感參數(shù)都是N-perv，但是互信息求得的R值大于偏秩相關(guān)法求得的P值絕對值，說明N-perv與峰值時(shí)間的非線性關(guān)系強(qiáng)于線性關(guān)系；對總產(chǎn)流量來說，兩種方法求得的最敏感參數(shù)都是Manning-N，且偏秩相關(guān)法求得的P值絕對值大于互信息法求得的R值，說明Manning-N和總產(chǎn)流量的線性關(guān)系強(qiáng)于非線性關(guān)系；兩種方法識別出的較敏感參數(shù)存在差異，這是因?yàn)槠认嚓P(guān)法是通過分析線性相關(guān)關(guān)系來識別敏感參數(shù)，互信息法是通過分析非線性相關(guān)關(guān)系來識別敏感參數(shù)，因此二者識別的敏感參數(shù)存在差異，但是這兩種方法可以共同使用，從而全面識別敏感參數(shù)，避免只用一種方法而導(dǎo)致的疏漏。

3 結(jié) 語

(1)偏秩相關(guān)分析方法能夠通過分析參數(shù)和輸出結(jié)果之間的線性關(guān)系來識別敏感參數(shù)，互信息法能夠通過分析參數(shù)和輸出結(jié)果之間的非線性關(guān)系來識別敏感參數(shù)，這兩種方法都是全局敏感新分析方法，將兩種方法結(jié)合使用能夠精確識別所有的敏感參數(shù)。

(2)通過兩種分析方法得知，Manning-N是對峰值流量影響最大的參數(shù)，Manning-N的增大將會使峰值流量減??；N-perv是對峰值時(shí)間影響最大的參數(shù)，N-perv的增大將會使峰值時(shí)間提前；Manning-N是對總產(chǎn)流量影響最大的參數(shù)，Manning-N的增大將會使總產(chǎn)流量減小。

(3)通過敏感性分析，精確識別所有敏感參數(shù)，能夠在率定過程中提高效率。

□

[1] 官奕宏,呂 謀,王 燦,等. 低影響開發(fā)技術(shù)的雨洪控制效果及水質(zhì)影響分析——基于SWMM模型[J]. 中國農(nóng)村水利水電,2017,411(1):84-87.

[2] 朱嘉祺,徐向陽,何 爽. 基于LH-OAT的SWMM模型參數(shù)敏感性分析[J]. 中國農(nóng)村水利水電,2014,377(3):21-26.

[3] 宋曉猛,張建云,占車生,等. 水文模型參數(shù)敏感性分析方法評述[J]. 水利水電科技進(jìn)展,2015,35(6):105-112.

[4] ZHOU N, PIERRE J W, TRUDNOWSKI D. A stepwise regression method for estimating dominant electromechanical modes[J]. IEEE Transactions on power systems. 2012,27(2):1 051-1 059.

[5] 王海霞,李 昱,周惠成,等. 基于Sobol’敏感性分析法的水庫多目標(biāo)調(diào)度模型求解質(zhì)量與效率[J]. 水電能源科學(xué),2016,34(6):43-47.

[6] 任啟偉,陳洋波,舒曉娟. 基于Extend FAST方法的新安江模型參數(shù)全局敏感性分析[J]. 中山大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,49(3):127-134.

[7] 王浩昌,杜鵬飛,趙冬泉,等. 城市降雨徑流模型參數(shù)全局靈敏度分析[J]. 中國環(huán)境科學(xué),2008,28(8):725-729.

[8] Chen Y, Wen J, Cheng S. Probabilistic load flow method based on nataf transformation and latin hypercube sampling[J]. IEEE Transactions on sustainable energy, 2013,4(2):294-301.

[9] 曹 淵,韓 峰,王鐵良,等. 爆炸容器力學(xué)安全的敏感度分析方法[J]. 科技導(dǎo)報(bào),2013,31(3):28-32.

[10] 熊劍智. 城市雨洪模型參數(shù)敏感性分析與率定[D]. 濟(jì)南:山東大學(xué),2016.

[11] 侯改娟. 綠色建筑與小區(qū)低影響開發(fā)雨水系統(tǒng)模型研究[D]. 重慶:重慶大學(xué),2014.