NAM降雨徑流模型的參數(shù)全局敏感性分析

趙然杭，伍 謀，王興菊，齊 真，周 璐，顧士升

（1.山東大學(xué) 土建與水利學(xué)院，山東 濟(jì)南 250061；2.山東省防汛抗旱物資儲(chǔ)備中心，山東 濟(jì)南 250013）

NAM模型由丹麥學(xué)者Nielsen和Hansen于1973年首次提出［1］，后經(jīng)丹麥水力研究所逐步完善，是用于模擬流域范圍內(nèi)由降雨產(chǎn)生的徑流過(guò)程的集總式概念性水文模型［2］。由于NAM模型結(jié)構(gòu)簡(jiǎn)單、模型參數(shù)較少，因此在我國(guó)的流域降雨徑流模擬中得到較為廣泛的應(yīng)用［3］。NAM模型不僅可以單獨(dú)用于流域降雨徑流模擬，也可以與MIKE11水動(dòng)力模型耦合進(jìn)行河道洪水過(guò)程模擬［4］。

模型概化使得模型與實(shí)際的物理過(guò)程差異較大，且NAM模型中很多參數(shù)難以通過(guò)實(shí)測(cè)方式獲取，需要通過(guò)模型率定獲得，導(dǎo)致參數(shù)值存在一定的不確定性，從而影響模擬預(yù)測(cè)結(jié)果［5-6］。參數(shù)敏感性分析是水文模型構(gòu)建與應(yīng)用的關(guān)鍵環(huán)節(jié)，其目的在于定性或定量評(píng)估模型參數(shù)對(duì)模擬結(jié)果的影響，確定參數(shù)的重要程度，從而識(shí)別敏感參數(shù)，以提高模型參數(shù)率定的效率，降低參數(shù)值的不確定性影響［7］，避免出現(xiàn)異參同果現(xiàn)象［8］。參數(shù)敏感性分析方法可根據(jù)作用范圍分為局部敏感性分析方法和全局敏感性分析方法，局部敏感性分析方法只分析單個(gè)參數(shù)對(duì)模型輸出結(jié)果的影響程度，而全局敏感性分析方法檢驗(yàn)多個(gè)參數(shù)及參數(shù)間相互作用對(duì)模型輸出結(jié)果產(chǎn)生的影響［9］。常見的全局敏感性分析方法有偏相關(guān)法［10］、Sobol法［11］、FAST法［12］等。

近年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)水文模型的敏感性進(jìn)行了相關(guān)研究。劉松等［13］聯(lián)合運(yùn)用Morris和Sobol法對(duì)三水源新安江模型進(jìn)行了全局輸出結(jié)果分析，研究了參數(shù)間的相關(guān)性對(duì)模型參數(shù)敏感性的影響；溫婭惠等［14］基于單目標(biāo)與多目標(biāo)GLUE方法研究了新安江模型參數(shù)的不確定性和敏感性，表明多目標(biāo)函數(shù)條件下的模型參數(shù)不確定性范圍較大；于永強(qiáng)等［15］基于LH-OAT敏感性分析方法分析了MIKE11模型參數(shù)的全局敏感性，表明降雨強(qiáng)度、雨型和評(píng)價(jià)目標(biāo)不同會(huì)影響模型參數(shù)的敏感性；Werkhoven等［16］對(duì)SAC-SMA模型參數(shù)進(jìn)行全局敏感性分析，優(yōu)先考慮敏感參數(shù)，降低了模型維數(shù)。Liu等［17］通過(guò)耦合Morris分析方法與NSDE算法對(duì)MIKE／NAM降雨徑流模型進(jìn)行參數(shù)敏感性分析與率定，篩選了模型敏感參數(shù)，提高了模型率定的效率；Zaghloul等［18］利用人工神經(jīng)網(wǎng)絡(luò)研究了SWMM模型參數(shù)的全局敏感性，并與重復(fù)試驗(yàn)結(jié)果進(jìn)行對(duì)比分析。目前對(duì)NAM模型的參數(shù)敏感性分析較少且主要針對(duì)其中部分參數(shù)的局部敏感性，如李磊等［19］利用擾動(dòng)分析法，分析了NAM模型中5個(gè)不同參數(shù)的局部敏感性；王振亞等［20］基于灰色系統(tǒng)關(guān)聯(lián)度分析方法研究了NAM模型中地表蓄水層最大含水量、淺層蓄水層最大含水量、地表徑流系數(shù)、壤中流系數(shù)、匯流時(shí)間常數(shù)對(duì)徑流深、洪峰流量和峰現(xiàn)時(shí)間模擬結(jié)果的影響；解恒燕等［21］基于擾動(dòng)分析法分析了NAM模型中9個(gè)主要參數(shù)的敏感性。以上研究表明，水文模型參數(shù)在高維空間通常會(huì)表現(xiàn)出較強(qiáng)的相關(guān)性［7］，針對(duì)NAM模型部分參數(shù)的局部敏感性分析，忽略了模型參數(shù)之間的相互作用對(duì)模擬結(jié)果的影響，因此有必要利用全局敏感性分析方法對(duì)NAM模型參數(shù)進(jìn)行敏感性分析。

筆者以小清河黃臺(tái)橋斷面以上321 km2為研究區(qū)域，采取拉丁超立方抽樣方法隨機(jī)抽取1 000組初始模型參數(shù)，利用互信息法和偏秩相關(guān)法對(duì)NAM模型參數(shù)進(jìn)行全局敏感性分析，研究NAM模型參數(shù)對(duì)模擬結(jié)果中的洪峰流量、峰現(xiàn)時(shí)間、徑流總量的影響，以進(jìn)一步提高模型模擬精度和減少建模工作量。

1 模型與方法

1.1 NAM模型

NAM模型是一個(gè)用于模擬自然流域降雨徑流過(guò)程的集總式概念性水文模型［22］。模型中影響流域水循環(huán)的土壤狀態(tài)被描述為一系列簡(jiǎn)化的數(shù)學(xué)語(yǔ)言形式，并將流域產(chǎn)匯流分為融雪蓄水層、地表蓄水層、淺層蓄水層、地下蓄水層4個(gè)相互關(guān)聯(lián)的蓄水層進(jìn)行模擬計(jì)算，每一層代表流域內(nèi)不同的物理元素。

NAM模型的參數(shù)及輸入量是定義在整個(gè)流域上的，即為流域的平均值［23］。模型要求的輸入數(shù)據(jù)有限，主要包括降雨量數(shù)據(jù)、蒸發(fā)量數(shù)據(jù)和模型參數(shù)。NAM模型的主要參數(shù)及其一般取值范圍見表1。

表1 NAM模型主要參數(shù)及其取值范圍

1.2 互信息法

互信息是指兩個(gè)及兩個(gè)以上變量之間共享的信息量［24］，用于度量隨機(jī)變量間的相互影響程度，即變量間的相關(guān)性大小，互信息越大，相關(guān)性越強(qiáng)。互信息法分析得出的R統(tǒng)計(jì)量可識(shí)別NAM模型參數(shù)與輸出的徑流總量、峰現(xiàn)時(shí)間、洪峰流量之間的非線性相關(guān)關(guān)系，用于NAM模型參數(shù)的全局敏感性分析。但R統(tǒng)計(jì)量只能識(shí)別參數(shù)敏感性的大小，無(wú)法識(shí)別NAM模型參數(shù)與輸出的徑流總量、峰現(xiàn)時(shí)間、洪峰流量之間的正負(fù)相關(guān)關(guān)系［25］。

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

式中：p（x）為離散型隨機(jī)變量X取值x的概率；p（y）為離散型隨機(jī)變量Y取值y的概率；p（x，y）為X取值x、Y取值y同時(shí)發(fā)生的概率。

1.3 偏秩相關(guān)法

偏秩相關(guān)分析是一種非參數(shù)統(tǒng)計(jì)分析方法，在等級(jí)位差的基礎(chǔ)上，通過(guò)控制其他因變量對(duì)輸出結(jié)果的影響來(lái)分析兩個(gè)變量之間的相關(guān)關(guān)系［26］，可準(zhǔn)確地反映自變量與因變量之間的線性關(guān)系，用于分析參數(shù)的全局敏感性。其P值（偏秩相關(guān)系數(shù)）的絕對(duì)值大小與計(jì)算值正負(fù)可分別識(shí)別參數(shù)敏感性的大小及輸入?yún)?shù)與輸出結(jié)果間線性關(guān)系的正負(fù)。

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

第i個(gè)輸入?yún)?shù)Xi與輸出結(jié)果Y之間的偏秩相關(guān)系數(shù)PXi為

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

選取小清河黃臺(tái)橋斷面以上321 km2為研究區(qū)域，該研究區(qū)域位于山東省濟(jì)南市，包括整個(gè)濟(jì)南市主城區(qū)和南部山區(qū)及部分西北郊區(qū)，屬溫帶季風(fēng)氣候區(qū)，降水主要集中于6—9月，多年平均降水量為671.1 mm。研究區(qū)域內(nèi)設(shè)有黃臺(tái)橋、東紅廟、燕子山、劉家莊、吳家鋪、紹而和興隆7個(gè)雨量站。

通過(guò)對(duì)研究區(qū)域內(nèi)7個(gè)雨量站2009—2013年的實(shí)測(cè)降雨資料及流量資料進(jìn)行分析研究，選用20120708、20130709和20130723場(chǎng)次暴雨洪水過(guò)程進(jìn)行參數(shù)率定，并以20110702、20120818和20130719場(chǎng)次暴雨為基準(zhǔn)進(jìn)行驗(yàn)證。根據(jù)率定期內(nèi)3場(chǎng)暴雨洪水過(guò)程得出流域率定參數(shù)，見表2，模擬計(jì)算結(jié)果見表3和圖1。

圖1 率定期與驗(yàn)證期洪水過(guò)程模擬結(jié)果

表2 NAM模型模擬參數(shù)率定結(jié)果

表3 各場(chǎng)次降雨洪水模擬結(jié)果

在3場(chǎng)率定洪水場(chǎng)次中，確定性系數(shù)R2分別為0.938、0.901和0.890，洪峰流量誤差分別為-1.67%、9.77%和-0.18%，峰現(xiàn)時(shí)間誤差分別為42、-24、60 min。在3場(chǎng)驗(yàn)證洪水場(chǎng)次中，確定性系數(shù)R2分別為0.915、0.903和0.883，洪峰流量誤差分別為-3.01%、-5.39%、-4.65%，峰現(xiàn)時(shí)間誤差分別為60、12、0 min。綜上，本模型對(duì)率定期和驗(yàn)證期共6場(chǎng)次暴雨洪水模擬計(jì)算的確定性系數(shù)均大于0.85，洪峰流量相對(duì)誤差均小于20%，峰現(xiàn)時(shí)間誤差均小于180 min，根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》可知，本模型的模擬精度達(dá)到良好水平，可以反映該流域流量的變化規(guī)律。說(shuō)明NAM模型在研究區(qū)域內(nèi)有較好的適用性，可用于模型參數(shù)的全局敏感性分析。

3 結(jié)果與分析

利用拉丁超立方抽樣方法在9個(gè)參數(shù)值各自取值范圍內(nèi)隨機(jī)抽取1 000組初始模型參數(shù)，通過(guò)MATLAB編寫腳本寫入NAM模型的輸入文件中，將1 000個(gè)輸入文件依次代入模型進(jìn)行模擬計(jì)算得到相應(yīng)的洪峰流量、峰現(xiàn)時(shí)間和徑流總量，再將參數(shù)值、模擬結(jié)果組成3個(gè)10×1 000的矩陣，利用互信息法和偏秩相關(guān)法對(duì)NAM模型進(jìn)行全局敏感性分析，敏感次序排第一和第二的參數(shù)即為敏感參數(shù)。

3.1 互信息分析

互信息分析方法得到的各參數(shù)對(duì)洪峰流量、峰現(xiàn)時(shí)間和徑流總量的敏感性分析結(jié)果見表4。

表4 互信息分析結(jié)果

根據(jù)表4中的R統(tǒng)計(jì)值可以判斷參數(shù)和輸出結(jié)果之間非線性關(guān)系的強(qiáng)弱，即參數(shù)敏感性的大小。

對(duì)于洪峰流量，9個(gè)參數(shù)的R統(tǒng)計(jì)值處于0.162～0.379之間，說(shuō)明各參數(shù)與洪峰流量的非線性關(guān)系較弱。其中，參數(shù)CK12的R統(tǒng)計(jì)值最大，為0.379，說(shuō)明與其他參數(shù)相比，CK12的參數(shù)敏感性較強(qiáng)。對(duì)于峰現(xiàn)時(shí)間，C K12的R統(tǒng)計(jì)值最大，為0.683，說(shuō)明CK12為最敏感參數(shù)，與峰現(xiàn)時(shí)間的非線性關(guān)系強(qiáng)于其他參數(shù)；第二敏感參數(shù)是TOF，R統(tǒng)計(jì)值為0.310。對(duì)于徑流總量，敏感參數(shù)分別為TOF和CK12，R統(tǒng)計(jì)值分別為0.434和0.426。

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

偏秩相關(guān)法得到的洪峰流量、峰現(xiàn)時(shí)間和徑流總量對(duì)各參數(shù)的敏感性分析結(jié)果見表5。

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

由表5中P值的絕對(duì)值大小可以得出參數(shù)和輸出結(jié)果之間線性相關(guān)關(guān)系的強(qiáng)弱，即參數(shù)的敏感性大小。P值正負(fù)可以反映參數(shù)和輸出結(jié)果之間線性相關(guān)關(guān)系是呈正相關(guān)還是負(fù)相關(guān)。

對(duì)于洪峰流量，最敏感的參數(shù)是TOF，P值為-0.598；第二敏感參數(shù)是CK12，P值為-0.418，說(shuō)明參數(shù)TOF和C K12與洪峰流量的線性相關(guān)關(guān)系強(qiáng)于其他參數(shù)，均為負(fù)相關(guān)關(guān)系。

對(duì)于峰現(xiàn)時(shí)間，參數(shù)CK12的P值最大，為0.817，說(shuō)明CK12為最敏感參數(shù)，與峰現(xiàn)時(shí)間的線性相關(guān)關(guān)系強(qiáng)于其他參數(shù)，呈正相關(guān)，在率定過(guò)程中可優(yōu)先率定參數(shù)C K12。

對(duì)于徑流總量，最敏感的參數(shù)是TOF，P值為-0.594；第二敏感參數(shù)是CK12，P值為-0.434，說(shuō)明TOF和CK12與徑流總量的線性關(guān)系強(qiáng)于其他參數(shù)，均為負(fù)相關(guān)關(guān)系。

3.3 結(jié)果對(duì)比分析

將偏秩相關(guān)分析的結(jié)果和互信息分析的結(jié)果進(jìn)行對(duì)比發(fā)現(xiàn)，對(duì)于洪峰流量，兩種方法求得的敏感參數(shù)都為TOF和CK12。在互信息分析結(jié)果中，參數(shù)CK12的R統(tǒng)計(jì)值為0.379，TOF的R統(tǒng)計(jì)值為0.282，說(shuō)明參數(shù)CK12的敏感性強(qiáng)于TOF；在偏秩相關(guān)分析結(jié)果中，參數(shù)CK12的P值為-0.418，TOF的P值為-0.598，說(shuō)明參數(shù)CK12的敏感性弱于TOF。

對(duì)于峰現(xiàn)時(shí)間，兩種方法求得的最敏感參數(shù)都是CK12，互信息法求得的R統(tǒng)計(jì)值為0.683，偏秩相關(guān)法求得的P值為0.817，說(shuō)明C K12與峰現(xiàn)時(shí)間的線性關(guān)系更為顯著。

對(duì)于徑流總量，兩種方法求得的最敏感參數(shù)都是TOF，互信息法求得的R統(tǒng)計(jì)值為0.434，偏秩相關(guān)法求得的P值為-0.594，說(shuō)明TOF與徑流總量的線性關(guān)系更為顯著。

3.4 討 論

對(duì)于洪峰流量，兩種方法求得的敏感參數(shù)均為CK12和TOF，但兩個(gè)參數(shù)的敏感次序不同，考慮到互信息法和偏秩相關(guān)法分別基于非線性和線性來(lái)識(shí)別敏感性參數(shù)，因此參數(shù)的敏感次序排列有所不同是不可避免的。對(duì)比其他學(xué)者對(duì)NAM模型參數(shù)敏感性分析結(jié)果，李磊等［19］的研究結(jié)果表明參數(shù)CK12對(duì)洪峰流量影響較大，與本文的結(jié)論一致。

對(duì)于峰現(xiàn)時(shí)間和徑流總量，互信息法和偏秩相關(guān)法求得的最敏感參數(shù)結(jié)果一致，峰現(xiàn)時(shí)間的最敏感參數(shù)均為CK12，R值為0.683，P值為0.817，說(shuō)明與峰現(xiàn)時(shí)間的線性關(guān)系較為顯著，呈正相關(guān)；徑流總量的最敏感參數(shù)均為TOF，R值為0.434，P值為-0.594，說(shuō)明與徑流總量的線性關(guān)系較為顯著，呈負(fù)相關(guān)。兩種方法的共同使用，避免了只使用單一方法的局限性，可全面地識(shí)別敏感參數(shù)。

對(duì)比其他學(xué)者對(duì)NAM模型參數(shù)敏感性分析結(jié)果，王振亞等［20］的研究表明參數(shù)CK12與峰現(xiàn)時(shí)間的關(guān)聯(lián)度最大，C QOF次之，本文中CQOF的敏感次序?yàn)?，其原因?yàn)榍罢呶纯紤]參數(shù)TOF對(duì)峰現(xiàn)時(shí)間的影響；解恒燕等［21］的研究表明參數(shù)Lmax和Umax對(duì)徑流總量的影響較為明顯，但其研究區(qū)域Fisher小流域地勢(shì)較為平坦且多為自然地貌，而本文研究區(qū)域小清河黃臺(tái)橋斷面以上321 km2流域，地勢(shì)南高北低且多為城市地形，說(shuō)明在不同下墊面情況的流域，對(duì)徑流總量影響最大的參數(shù)不同。

4 結(jié) 論

針對(duì)NAM模型參數(shù)，利用拉丁超立方抽樣方法抽取模型輸入樣本，使用互信息法和偏秩相關(guān)法進(jìn)行NAM模型參數(shù)全局敏感性分析，得到以下結(jié)論。

（1）通過(guò)拉丁超立方抽樣方法抽取了1 000組模型輸入樣本，保證了NAM模型參數(shù)值數(shù)據(jù)的代表性和均勻性?；バ畔⒎ㄍㄟ^(guò)分析參數(shù)與模擬結(jié)果間的非線性關(guān)系識(shí)別敏感參數(shù)，偏秩相關(guān)法通過(guò)分析參數(shù)與模擬結(jié)果間的線性關(guān)系識(shí)別敏感參數(shù)，且兩者均為全局敏感性分析方法，兩種方法共同使用可避免只使用單一方法存在的局限性，全面準(zhǔn)確地識(shí)別敏感參數(shù)，保證結(jié)論的合理性及可靠性。

（2）在兩種不同敏感性分析方法結(jié)果中，參數(shù)TG、TIF、CKIF、Umax和Lmax的R統(tǒng)計(jì)值和P值均較小，說(shuō)明這些參數(shù)為不敏感參數(shù)，在模型率定時(shí)可以根據(jù)經(jīng)驗(yàn)取固定值以提高率定效率。TOF和CK12是對(duì)洪峰流量敏感性最大的參數(shù)，且與洪峰流量呈負(fù)相關(guān)；CK12是對(duì)峰現(xiàn)時(shí)間敏感性最大的參數(shù)，隨著CK12增大，峰現(xiàn)時(shí)間會(huì)相應(yīng)延長(zhǎng)；TOF是對(duì)徑流總量最為敏感的參數(shù)，呈負(fù)相關(guān)，隨著TOF增大，徑流總量會(huì)相應(yīng)減小。

（3）通過(guò)全局敏感性分析，識(shí)別了NAM模型的敏感參數(shù)，為后續(xù)模型率定及耦合提供參考和指導(dǎo)，減少建模的工作量，提高模型率定的效率。