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

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

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

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

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

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

1 模型與方法

1.1 NAM模型

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

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

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

1.2 互信息法

互信息是指兩個及兩個以上變量之間共享的信息量［24］，用于度量隨機(jī)變量間的相互影響程度，即變量間的相關(guān)性大小，互信息越大，相關(guān)性越強(qiáng)?；バ畔⒎ǚ治龅贸龅腞統(tǒng)計量可識別NAM模型參數(shù)與輸出的徑流總量、峰現(xiàn)時間、洪峰流量之間的非線性相關(guān)關(guān)系，用于NAM模型參數(shù)的全局敏感性分析。但R統(tǒng)計量只能識別參數(shù)敏感性的大小，無法識別NAM模型參數(shù)與輸出的徑流總量、峰現(xiàn)時間、洪峰流量之間的正負(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同時發(fā)生的概率。

1.3 偏秩相關(guān)法

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

3.1 互信息分析

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

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

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

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

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

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

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

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

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

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

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

3.3 結(jié)果對比分析

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

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

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

3.4 討 論

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

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

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

4 結(jié) 論

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

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

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

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