于永強(qiáng)， 沙曉軍， 劉 俊， 周 宏， 班 超

(1. 河海大學(xué)水文水資源學(xué)院，南京 210098；2. 江蘇省江陰市水資源管理辦公室，江蘇 江陰 214400)

MIKE11模型是一款比較成熟的模型，通過(guò)耦合降雨徑流模塊和水動(dòng)力模塊實(shí)現(xiàn)徑流模擬，廣泛應(yīng)用于防洪、水資源保護(hù)及水利工程設(shè)計(jì)管理等方面。MIKE11模型已經(jīng)得到廣泛應(yīng)用，在使用中也產(chǎn)生了很多問(wèn)題和困惑。模型的參數(shù)過(guò)多，不能與現(xiàn)有資料很好地銜接。使得確定參數(shù)時(shí)帶來(lái)的誤差會(huì)在模型運(yùn)行中積累，影響模型精度[1]。

目前國(guó)內(nèi)外對(duì)MIKE11模型敏感性分析較少且主要針對(duì)單一參數(shù)局部敏感性[2]，較少考慮多種參數(shù)共同影響下的全局參數(shù)敏感性。水文模型敏感性分析方法主要有局部分析法、隨機(jī)OAT法、全局抽樣法。全局抽樣法可以提供客觀有效的全局解但計(jì)算量太大，不適合復(fù)雜模型。局部分析法簡(jiǎn)單高效但輸入輸出不呈線性關(guān)系，不滿足模型高效性和穩(wěn)定性要求。本次研究采用LH-OAT參數(shù)敏感性分析法，用LH抽樣法代替MC抽樣法得到改進(jìn)的OAT法。綜合LH法的穩(wěn)定高效性和OAT法的精確性，保證參數(shù)敏感性分析的魯棒性和有效性[3,4]。該方法在Brussels、Texas等流域，SWAT和SWMM等模型敏感性分析和參數(shù)率定中都獲得良好效果[5,6]。參數(shù)敏感性分析得出不同參數(shù)對(duì)模型輸出的影響程度，將高敏感性參數(shù)作為主要率定參數(shù)。控制模型計(jì)算的不確定性，提高模型率定效率和模擬精度，提高M(jìn)IKE11模型使用效率。

1 模型與方法

1.1 MIKE11模型

MIKE11水動(dòng)力學(xué)模塊的數(shù)學(xué)表達(dá)基于明渠不穩(wěn)定流隱式格式有限差分解，控制方程為一維非恒定流Saint-Venant方程組，見(jiàn)式(1)。方程組利用Abbott-Ionescu六點(diǎn)隱式有限差分格式求解，每一個(gè)網(wǎng)格點(diǎn)按順序交替計(jì)算水位或流量，離散后的線性方程組用追趕法求解。

(1)

式中：x、t分別為計(jì)算點(diǎn)空間和時(shí)間的坐標(biāo)；A為過(guò)水?dāng)嗝婷娣e；Q為過(guò)流流量；h為水位；q為旁側(cè)入流流量；C為謝才系數(shù)；R為水力半徑；α為動(dòng)量校正系數(shù)；g為重力加速度。

本次降雨徑流模塊采用NAM模型，NAM模型作為一種概念性的集總模型，用模型參數(shù)反映流域單元整體狀態(tài)，對(duì)產(chǎn)流過(guò)程進(jìn)行簡(jiǎn)化和定量模擬。NAM模型包含了積雪儲(chǔ)水層，地表儲(chǔ)水層，根區(qū)儲(chǔ)水層和地下儲(chǔ)水層4個(gè)不同且相互影響的儲(chǔ)水層，這4個(gè)儲(chǔ)水層反映了流域內(nèi)不同物理狀態(tài)，考慮了融雪、截留填洼、植物蒸騰、土壤根系蓄水和地下儲(chǔ)水等水量交換過(guò)程，其中積雪儲(chǔ)水層可不考慮。通過(guò)由上及下依次連續(xù)計(jì)算4個(gè)儲(chǔ)水層來(lái)實(shí)現(xiàn)降雨徑流模擬，最終計(jì)算結(jié)果以蒸發(fā)量、坡面流、壤中流和基流4部分輸出。

1.2 LH-OAT參數(shù)敏感性分析法

LH-OAT參數(shù)敏感性分析法以歸一化的指標(biāo)來(lái)體現(xiàn)參數(shù)變化對(duì)于模型輸出的影響程度，即參數(shù)敏感性。LH-OAT法兼顧了OAT法的精確性和LH抽樣法的健壯穩(wěn)定性，引入Latin-Hypercube抽樣算法，相對(duì)于Monte-Carlo抽樣法簡(jiǎn)便，同時(shí)也很好地兼顧了不同參數(shù)區(qū)間。

在抽樣時(shí)，第一步是把各個(gè)參數(shù)區(qū)間等分為n個(gè)子區(qū)間，保證在抽樣時(shí)每個(gè)參數(shù)子區(qū)間被抽到的概率相等，都為1/n；第二步，將m個(gè)參數(shù)進(jìn)行隨機(jī)抽樣，使每個(gè)子區(qū)間都抽出且只抽出一個(gè)隨機(jī)參數(shù)；最后，將各個(gè)隨機(jī)參數(shù)任意組合，得到n組隨機(jī)參數(shù)組合帶入模型計(jì)算。模型計(jì)算結(jié)果的差異取決于所有參數(shù)的變化，無(wú)法確定每個(gè)參數(shù)對(duì)于影響模型輸出結(jié)果的貢獻(xiàn)程度。因此每次對(duì)n組隨機(jī)參數(shù)序列中的一個(gè)參數(shù)進(jìn)行微小擾動(dòng)fk，擾動(dòng)幅度分別取±10%，±20%和±30%[7,8]，以減小不同隨機(jī)參數(shù)組合下擾動(dòng)幅度大小對(duì)于參數(shù)敏感性的影響。綜上可得到n(6m+1)組參數(shù)序列，依次帶入模型計(jì)算得到結(jié)果。

為了方便比較敏感性差異，對(duì)參數(shù)敏感性進(jìn)行歸一化處理得到相對(duì)敏感度。本次研究考慮了6種擾動(dòng)幅度，改進(jìn)的敏感度計(jì)算公式如下：

(2)

式中：RSij為第i個(gè)參數(shù)的第j組參數(shù)序列的相對(duì)敏感度；fk為第k次施加的微小擾動(dòng)；M為模型函數(shù)；M[e1j,…,eij(1+fk),…,emj]為對(duì)第j組參數(shù)序列中第i個(gè)參數(shù)施加擾動(dòng)fk后模型輸出結(jié)果。

將得到的RSij求均值，見(jiàn)下式：

(3)

得到歸一化的敏感度RSi參照表1確定參數(shù)敏感性。

表1 敏感性取值表Tab.1 Sensitivity values Table

2 研究區(qū)域和模型構(gòu)建

本次研究區(qū)域?yàn)榻K省張家港市中南部相對(duì)獨(dú)立的河網(wǎng)區(qū)域，位于長(zhǎng)江下游西南岸，江蘇省東南部，共230.2 km2。張家港市年平均降水量為1 025.7 mm，地區(qū)多年平均蒸發(fā)量1 395.7 mm。張家港市河流屬長(zhǎng)江流域太湖水系，河、港、套、塘、浦縱橫貫通，交織成網(wǎng)。長(zhǎng)江環(huán)繞西北、北和東北面，張家港河、二干河與太湖相通。

研究區(qū)域共有大小河道3 557條，總長(zhǎng)1 326.99 km。本次研究概化了市級(jí)及以上河道15條，總長(zhǎng)131.48 km，鎮(zhèn)級(jí)河道114條，總長(zhǎng)285.25 km，水系現(xiàn)狀圖見(jiàn)圖1。邊界條件包含內(nèi)河邊界和外河邊界，由太湖流域水文動(dòng)力模型HOHY2得出奚浦塘、曲塘、鹽鐵塘等8個(gè)外河邊界；由張家港市河網(wǎng)水利模型得出二干河、三干河、四干河等11個(gè)內(nèi)河邊界。根據(jù)張家港市防洪排澇形勢(shì)，確定張家港市中心城區(qū)防洪標(biāo)準(zhǔn)為100年一遇，河道排澇標(biāo)準(zhǔn)為20年一遇。研究區(qū)域匯水面積較小，匯水歷時(shí)較短，選擇暴雨最大控制時(shí)段為24 h，降雨輸入選擇了100年一遇雨型1、20年一遇雨型1和20年一遇雨型2的設(shè)計(jì)暴雨過(guò)程，分別記為降雨過(guò)程A、B和C，見(jiàn)圖2?？偨Y(jié)了MIKE11的10個(gè)重要參數(shù)，其中RR模塊相關(guān)參數(shù)取值范圍采用用戶手冊(cè)推薦范圍，HD模塊相關(guān)參數(shù)取值范圍由天然流域河道情況決定，見(jiàn)表2。對(duì)無(wú)斷面資料的小型河道概化為規(guī)范化梯形斷面，無(wú)詳細(xì)坡度資料地區(qū)采用平均坡度。初始水位和流量采用張家港市各河道常水位及相應(yīng)流量。

3 結(jié)果與分析

將每個(gè)模型參數(shù)取值范圍分成8層進(jìn)行LH抽樣，得到8組隨機(jī)參數(shù)。在對(duì)參數(shù)敏感性進(jìn)行OAT分析時(shí)，分別施加±10%、±20%和±30%的微小擾動(dòng)，因此模型針對(duì)一場(chǎng)降雨過(guò)程需運(yùn)行共8×(1+10×6)=536次。本次研究選取不同參數(shù)條件下的徑流量和洪峰流量作為敏感度評(píng)價(jià)指標(biāo)，對(duì)模型輸出結(jié)果進(jìn)行歸一化處理后得到不同雨型和不同降雨強(qiáng)度情況下的參數(shù)全局敏感度，見(jiàn)表3。

圖1 張家港市水系現(xiàn)狀圖Fig.1 Current situation of water system in Zhangjiagang City

圖2 設(shè)計(jì)降雨過(guò)程Fig.2 Design rainfall process

參數(shù)名稱物理意義取值范圍Umax地表儲(chǔ)水層最大含水量/mm10～25Lmax根區(qū)儲(chǔ)水層最大含水量/mm50～250CQOF坡面流產(chǎn)流系數(shù)0．01～0．90CKIF壤中流匯流時(shí)間常數(shù)/h500～100CK12坡面流匯流時(shí)間常數(shù)/h3～48TOF坡面流根區(qū)土壤含水量閾值0～0．7TIF壤中流根區(qū)土壤含水量閾值0～0．7TG地下水補(bǔ)充的根區(qū)閾值0～0．7CKBF基流計(jì)算時(shí)間常數(shù)/h500～5000nManning河道糙率0．02～0．04

以徑流量作為評(píng)價(jià)目標(biāo)，從模型參數(shù)上來(lái)看，MIKE11模型參數(shù)靈敏度針對(duì)不同降雨過(guò)程，會(huì)產(chǎn)生相應(yīng)波動(dòng)。然而在圖3(A1～C1)中，敏感度都為敏感的參數(shù)為L(zhǎng)max、CQOF和TOF，三者都與坡面流產(chǎn)流關(guān)系緊密。其次，TG和CK12在3場(chǎng)暴雨過(guò)程中都為敏感參數(shù)，但在100年一遇和20年一遇時(shí)，TG敏感度分別為敏感和一般敏感；CK12敏感度分別為一般敏感和敏感。其他參數(shù)如TIF、CKBF、Umax等對(duì)徑流量影響甚微。從降雨過(guò)程來(lái)看，雨強(qiáng)對(duì)參數(shù)敏感性影響較大，CQOF和TG隨著雨強(qiáng)增大敏感性提高，CK12隨著雨強(qiáng)增大敏感性降低。雷達(dá)圖B1和C1基本擬合，可見(jiàn)雨型對(duì)于參數(shù)敏感性影響較小。

以洪峰流量作為評(píng)價(jià)目標(biāo)，從模型參數(shù)上來(lái)看，河道糙率n-Manning是最為敏感的，在3場(chǎng)暴雨過(guò)程中敏感度都達(dá)到0.6以上。從降雨過(guò)程來(lái)看，當(dāng)降雨強(qiáng)度增大時(shí)，主要敏感參數(shù)n-Manning的敏感度下降25%，且次要敏感參數(shù)明顯增多，給模型率定增加了難度。對(duì)比雷達(dá)圖B2和C2可知，當(dāng)雨峰靠前時(shí)，Lmax與CK12由不敏感參數(shù)變?yōu)槊舾袇?shù)。圖3中，A～C分別表降雨過(guò)程A、B和C；1為徑流量；2為洪峰流量。

圖3 不同降雨的參數(shù)全局靈敏度對(duì)比Fig.3 Comparison of the parameter sensitivity in different rainfall

綜上所述，降雨強(qiáng)度、雨型和評(píng)價(jià)目標(biāo)的變化會(huì)改變MIKE11模型參數(shù)敏感性。降雨強(qiáng)度對(duì)洪峰流量參數(shù)敏感性的影響大于對(duì)徑流量參數(shù)敏感性的影響。雨型對(duì)徑流量參數(shù)敏感性基本無(wú)影響，而對(duì)洪峰流量參數(shù)敏感性有一定影響。以上結(jié)果與Drechsler和Lei的相關(guān)研究結(jié)果[9,10]相映證。

采用假定最優(yōu)參數(shù)法分析參數(shù)敏感性對(duì)率定效率的影響。預(yù)先選定一套參數(shù)作為最優(yōu)參數(shù)系列，分別用兩套方案進(jìn)行率定：方案1，隨機(jī)選擇參數(shù)順序進(jìn)行率定；方案2，按照敏感度由大到小順序進(jìn)行率定。以率定輸出結(jié)果與最優(yōu)參數(shù)輸出結(jié)果偏離程度評(píng)價(jià)率定效率，結(jié)果見(jiàn)圖4。方案1率定17次，精度達(dá)到98.3%，方案2率定10次，精度達(dá)到99.6%。由此可見(jiàn)優(yōu)選率定參數(shù)可提高59%以上的率定效率。當(dāng)率定目標(biāo)為洪峰流量時(shí)，MIKE11模型參數(shù)敏感性受雨型、雨量影響大，給優(yōu)選率定參數(shù)帶來(lái)困難，降低率定效率。以洪峰流量為率定目標(biāo)時(shí)，不同重現(xiàn)期的，不同洪峰滯時(shí)的，以及多雨峰的降雨會(huì)產(chǎn)生不同的參數(shù)敏感性。通過(guò)分析和總結(jié)其變化規(guī)律，可以對(duì)MIKE11模型參數(shù)敏感性有更好的認(rèn)識(shí)，提高M(jìn)IKE11模型運(yùn)用水平。

圖4 率定效果對(duì)比Fig.4 Comparison of the calibration results

4 結(jié) 語(yǔ)

(1)利用LH-OAT法分析MIKE11模型參數(shù)全局敏感性。以徑流量為評(píng)價(jià)目標(biāo)時(shí)主要敏感參數(shù)由高到低依次為L(zhǎng)max、CQOF、TOF、TG和CK12；以洪峰流量為評(píng)價(jià)目標(biāo)時(shí)主要敏感參數(shù)為n-Manning。

(2)雨型、雨量和評(píng)價(jià)目標(biāo)變化會(huì)使參數(shù)敏感性產(chǎn)生波動(dòng)。以徑流量為評(píng)價(jià)目標(biāo)時(shí)，雨型對(duì)參數(shù)敏感性影響較小，當(dāng)雨量增大時(shí)，CK12敏感性降低，TG敏感性提高；以洪峰流量為評(píng)價(jià)目標(biāo)時(shí)，雨峰靠前會(huì)增加Lmax與CK12的敏感性，雨量增大會(huì)坦化各參數(shù)敏感性之間的差異。

(3)優(yōu)選率定參數(shù)可大幅度提高M(jìn)IKE11模型率定效率，應(yīng)根據(jù)模型實(shí)現(xiàn)功能(洪峰預(yù)警、洪量預(yù)報(bào))和輸入條件(雨量、雨型)選擇相應(yīng)參數(shù)率定次序。

□

[1] Griensven A V, Meixner T, Grunwald S, et al. A global sensitivity analysis tool for the parameters of multi-variable catchment models[J]. Journal of Hydrology, 2006,324(1-4):10-23.

[2] 許春東, 殷 丹. 基于MIKE11的河道糙率靈敏度分析[J]. 水電能源科學(xué), 2014,(11):101-103.

[3] 徐會(huì)軍, 陳洋波, 李晝陽(yáng), 等. 基于LH-OAT分布式水文模型參數(shù)敏感性分析[J]. 人民長(zhǎng)江, 2012,43(7):19-23.

[4] 田 雨, 雷曉輝, 蔣云鐘,等. 水文模型參數(shù)敏感性分析方法研究評(píng)述[J]. 水文, 2010,30(4):9-12.

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

[6] Leta O T, Nossent J, Griensven A V, et al. Incorporating rainfall uncertainty in a SWAT model: the river Zenne basin (Belgium) case study[C]∥ Egu General Assembly, 2013.

[7] Helton J C, Davis F J. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[J]. Reliability Engineering & System Safety, 2003,81(1):23-69.

[8] 黃金良, 杜鵬飛, 何萬(wàn)謙, 等. 城市降雨徑流模型的參數(shù)局部靈敏度分析[J]. 中國(guó)環(huán)境科學(xué), 2007,27(4):549-553.

[9] Martin Drechsler. Sensitivity analysis of complex models[J]. Biological Conservation, 1998,86(3):401-412.

[10] Lei J, Schilling W. Parameter uncertainty propagation analysis for urban rainfall runoff modeling [J]. Water Science and Technology, 1994,29(1-2):145-154.