劉凌雪，敖天其,2，黎小東，胡 正，胡富昶，李孟芮

(1.四川大學(xué)水利水電學(xué)院, 成都 610065; 2.四川大學(xué)水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 成都 610065)

隨著水文氣象資料和下墊面信息的愈加豐富及計(jì)算機(jī)科學(xué)的發(fā)展，分布式水文物理模型成為當(dāng)今水文科學(xué)研究的前沿方向之一[1]。與傳統(tǒng)的集總式模型相比，分布式水文模型在擁有眾多優(yōu)點(diǎn)的同時(shí)也面臨著因參數(shù)過(guò)多而難以率定的問(wèn)題。對(duì)于足夠復(fù)雜且具有高維度和高計(jì)算需求的分布式水文模型而言，伴隨著參數(shù)相互作用的過(guò)度參數(shù)化導(dǎo)致模型參數(shù)不是唯一可識(shí)別的，Beven[2]將此現(xiàn)象稱為“異參同效”現(xiàn)象，即不同的參數(shù)集會(huì)導(dǎo)致相同或相似的模擬結(jié)果。避免過(guò)度參數(shù)化或不可識(shí)別性的一種可能的方式是減少參數(shù)個(gè)數(shù)[3]。

模型敏感性分析(SA)是過(guò)濾模型關(guān)鍵參數(shù)的有效工具[4]，它旨在從不敏感的參數(shù)中識(shí)別出最敏感的參數(shù)以減小參數(shù)維數(shù)[5]，降低參數(shù)率定的不確定性，提高模型的參數(shù)優(yōu)化效率。目前，SA方法已廣泛用于許多領(lǐng)域的實(shí)際研究中：Collins和Avissar[6]采用傅立葉振幅靈敏度測(cè)試(FAST)方法評(píng)估LAID模型中參數(shù)對(duì)顯熱和潛熱的影響情況；Werkhoven等[7]將全局SA方法作為篩選工具來(lái)降低參數(shù)維數(shù)，進(jìn)行參數(shù)的多目標(biāo)優(yōu)化率定，并將其成功應(yīng)用于SAC-SMA模型；Sun等[8]考慮了多目標(biāo)敏感性分析和多目標(biāo)條件的參數(shù)優(yōu)化，將Morris篩選法與非支配排序差分進(jìn)化(NSDE)算法耦合并用于MIKE/NAM這一降雨徑流模型的參數(shù)敏感性分析和率定中。

BTOPMC模型[9]屬分布式水文物理模型，有5個(gè)可以定量反映流域物理特征的參數(shù)，主要涉及每個(gè)網(wǎng)格單元的土壤類型，植被和土地利用(土地覆蓋)。盡管該模型的參數(shù)具有明確的物理意義，但由于測(cè)量困難，參數(shù)值往往通過(guò)參數(shù)率定來(lái)確定。基于BTOPMC模型參數(shù)的物理意義及以往的模擬經(jīng)驗(yàn)，在不同的研究區(qū)，由于流域物理特征及類型的不同，模型參數(shù)需要分為不同類型的土壤、植被、地形以及地表覆蓋情況來(lái)率定[10]，這就使得BTOPMC的參數(shù)率定過(guò)程通常需要較長(zhǎng)時(shí)間，且由于“異參同效”現(xiàn)象使得參數(shù)最優(yōu)值也存在較大的不確定性。本文分別采用定性分析方法MOAT[11]和定量分析方法DGSM[12]對(duì)BTOPMC模型的參數(shù)進(jìn)行敏感性分析，篩選較敏感參數(shù)以減少需要被率定的模型參數(shù)個(gè)數(shù)，從而達(dá)到減少參數(shù)率定值的不確定性、提高參數(shù)率定效率的目的。

1 BTOPMC模型

BTOPMC是由日本山梨大學(xué)[9]在TOPMODEL的概念基礎(chǔ)上加入了洼地消除、數(shù)字河網(wǎng)生成、馬斯京根-康奇法等子模型而開(kāi)發(fā)的一個(gè)分布式水文物理模型。BTOPMC模型的結(jié)構(gòu)組成見(jiàn)圖1[13]。目前，該模型在河網(wǎng)技術(shù)、自然子流域劃分、產(chǎn)匯流模型參數(shù)以網(wǎng)格為單位的空間分布給定方法、匯流精度等方面取得了明顯進(jìn)步，涉及的主要內(nèi)容有[10]：①發(fā)現(xiàn)了按數(shù)字高程模型(DEM)中的洼地的空間位置計(jì)算高程微增量，從而消除洼地的自動(dòng)解析新技術(shù)；②研究了Pfafstertter流域自動(dòng)分割法用于分布式水文模型的問(wèn)題并提出改進(jìn)措施；③提出了用Pfafstertter法分割流域時(shí)的穩(wěn)定分割程度的判別方法；④開(kāi)發(fā)了解決馬斯京跟匯流方法中的負(fù)出流問(wèn)題的近似處理方法和時(shí)空步長(zhǎng)實(shí)時(shí)調(diào)節(jié)方法；⑤提出了將TOPMODEL用作大流域的產(chǎn)流模型時(shí)，以各網(wǎng)格(而不是由網(wǎng)格組成的子流域)中土壤、植被等流域物理特性給定個(gè)模型參數(shù)的新模式；⑥提出了用Pareto解析法減少模型參數(shù)優(yōu)化過(guò)程及建立參數(shù)轉(zhuǎn)移函數(shù)過(guò)程中的不確定性的必要性。

圖1 BTOPMC模型結(jié)構(gòu)組成Fig.1 The structure of BTOPMC model

BTOPMC模型具有5個(gè)需要被率定的參數(shù)?；贐TOPMC模型參數(shù)的物理意義及以往的模擬經(jīng)驗(yàn)，在不同的研究區(qū)，由于流域物理特征及類型的不同，模型參數(shù)需要分為不同類型的土壤、植被、地形以及地表覆蓋情況來(lái)率定[10]。根據(jù)所使用的數(shù)據(jù)源[14]，本研究中模型參數(shù)的具體分類見(jiàn)表1，其中“標(biāo)識(shí)”列中的Xi為對(duì)應(yīng)于不同參數(shù)的在本研究中需要被率定的變量。

表1 BTOPMC模型參數(shù)表[10]Tab.1 Model parameters of BTOPMC[10]

BTOPMC模型中用于模型精度評(píng)價(jià)的指標(biāo)為Nash-Sutcliffe效率(NSE)[14]。NSE的數(shù)值越大，模擬結(jié)果越好，其計(jì)算公式如下：

NSE=

(1)

式中：N為總的時(shí)間步長(zhǎng)；Qsim和Qobs分別為時(shí)間步長(zhǎng)t所對(duì)應(yīng)的模擬流量和實(shí)測(cè)流量；Qav是整個(gè)時(shí)期觀測(cè)流量的平均值。

2 研究區(qū)域及資料

富士川流域位于日本本州島中央，是日本三大險(xiǎn)峻河流之一，屬一級(jí)河川，研究區(qū)范圍示意圖見(jiàn)圖2。河流主干流長(zhǎng)度約128 km，流域控制面積約786.2 km2。全流域大概88%被森林覆蓋，流域內(nèi)最高點(diǎn)是著名的富士山。全流域平均降雨量約1 550 mm，年平均氣溫在4～18 ℃之間。年內(nèi)，4月的降雨較3月有明顯增多；6月，氣壓配置成為“梅雨型”，直到7月上、中旬梅雨結(jié)束。當(dāng)梅雨氣壓停滯日本上空時(shí)，若有南方熱帶低氣壓出現(xiàn)，舌狀的暖氣流就會(huì)進(jìn)入日本上空，“濕舌”與“梅雨”相遇便會(huì)導(dǎo)致流域內(nèi)出現(xiàn)暴雨天氣。8月末、9月上旬至9月末、10月上旬是易降雨的季節(jié)，而且這一時(shí)期也是臺(tái)風(fēng)多發(fā)期，降大雨的可能性很大。

圖2 研究區(qū)范圍示意圖Fig.2 Schematic diagram of the study area

本文選取富士川流域1991年9月的一場(chǎng)降雨過(guò)程，7個(gè)測(cè)站的降雨資料均來(lái)自氣象資料自動(dòng)集取系統(tǒng)(AMeDAS)的每小時(shí)降雨觀測(cè)記錄，基本無(wú)數(shù)據(jù)丟失。年潛在蒸散發(fā)量采用Penman-Monteith方法估計(jì)為1 800 mm。對(duì)于洪水的模擬，潛在蒸散發(fā)量被假定為5.5 mm/h。流域的植被覆蓋圖來(lái)自IGBP第二版(USGS)，土壤圖來(lái)自FAO，且植被和土壤的分布情況是對(duì)每一個(gè)網(wǎng)格給出的[14]。

3 研究方法

近年來(lái)，基于實(shí)驗(yàn)設(shè)計(jì)(DoE)的全局SA方法十分受歡迎，原因在于它們?cè)诒３钟?jì)算效率的同時(shí)還可定量表征全局敏感度。典型的基于DoE的SA過(guò)程主要有兩個(gè)步驟：首先，使用所選擇的取樣方法在可行參數(shù)空間內(nèi)生成一組樣本參數(shù)；然后，利用SA方法獲得不同參數(shù)變化對(duì)模型輸出的影響情況。目前，有許多常用的DoE方法，如蒙特卡洛(MC)、拉丁超立方體(LH)[15]、對(duì)稱拉丁超立方體(SLH)等。

3.1 取樣方法

本文初步選取了4種DoE方法：MC是一種最常用的DoE方法，它依賴于重復(fù)隨機(jī)采樣來(lái)獲得未知概率實(shí)體的指定分布函數(shù)的近似數(shù)值，其斂散性依賴于獨(dú)立的隨機(jī)參數(shù)的個(gè)數(shù)；準(zhǔn)蒙特卡洛(QMC)是確定性的MC方法，也被稱為低偏差序列，它是一組填充樣本區(qū)域的“有效”點(diǎn)[16]；LH是一種根據(jù)多維分布生成合理的參數(shù)集合分布的統(tǒng)計(jì)方法，它適合于任意維數(shù)的抽樣設(shè)計(jì)；SLH與LH均具有“充滿空間”性質(zhì)，理論上來(lái)講，LH的試驗(yàn)點(diǎn)帶有隨機(jī)性，抽樣表現(xiàn)不穩(wěn)定，SLH則能夠產(chǎn)生分布性更好的試驗(yàn)點(diǎn)。

對(duì)于BTOPMC，我們使用300個(gè)采樣點(diǎn)來(lái)評(píng)估這4種方法的空間填充屬性，分別重復(fù)運(yùn)行MC、LH和SLH 6次，因QMC是一個(gè)確定性序列所以只運(yùn)行一次，計(jì)算空間填充性能好壞的度量數(shù)值：修正的L2偏差(MD2)，中心化L2偏差(CD2)，對(duì)稱化L2偏差(SD2)和可卷的L2偏差(WD2)。對(duì)于這些度量，數(shù)值越小，DoE算法的空間填充性能越好[17]。4種DoE算法的比較結(jié)果見(jiàn)圖3。對(duì)比分析之后，我們選擇LH作為后續(xù)BTOPMC敏感性分析的DoE方法。

3.2 SA方法

Morris One At A Time(MOAT)是一種簡(jiǎn)單、高效的參數(shù)篩選方法，通過(guò)少量的模型計(jì)算就可以獲得模型參數(shù)的定性排序，尤其適合參數(shù)較多的復(fù)雜模型[11]。該方法雖然可能將不重要參數(shù)判斷為重要參數(shù)，但不會(huì)出現(xiàn)將重要參數(shù)判為不重要參數(shù)的錯(cuò)誤[18]。每個(gè)參數(shù)的總體效果可以利用所采取樣本的對(duì)應(yīng)參數(shù)梯度的均值μ和標(biāo)準(zhǔn)差δ來(lái)近似估計(jì)(梯度是指目標(biāo)函數(shù)值的差值與樣本參數(shù)值的差值之比)。均值μ的大小表征參數(shù)的敏感度，而標(biāo)準(zhǔn)差δ表征參數(shù)之間相互作用的程度。Campolongo等[19]提出了修正的均值μ*，以解決部分參數(shù)對(duì)均值μ存在負(fù)效應(yīng)的問(wèn)題，公式如下：

(2)

di=[f(X1,…,Xi-1,Xi+Δ,Xi+1,…,Xn)-f(X)]/Δ

(3)

其中X=(X1,X2,…,Xn)是一個(gè)試驗(yàn)范圍內(nèi)的任意選定值，是一個(gè)n維p級(jí)的網(wǎng)格，Xi可以從{0,1/(p-1),2/(p-1),…,1}中取值；Δ是1/(p-1)預(yù)定的倍數(shù)，當(dāng)p為偶數(shù)時(shí)，Δ通常取p/[2×(p-1)]。

對(duì)于MOAT方法而言，μ*越大，參數(shù)越敏感。本文中，我們使用修正的平均值μ*作為MOAT方法的敏感性度量。此外，MOAT[11]有特定的取樣方法―Morris，所取樣本大小通常設(shè)置為n+1的倍數(shù)，其中n是參數(shù)個(gè)數(shù)。由表1可得，本研究中n取19。

圖3 四種DoE方法對(duì)應(yīng)的度量值比較Fig.3 Comparison of metrics for the four DoE methods

DGSM[12]是一種新的基于導(dǎo)數(shù)的全局SA方法，該方法隸屬于微分法，它是建立在求解響應(yīng)函數(shù)對(duì)輸入變量在多個(gè)點(diǎn)處的偏導(dǎo)數(shù)的基礎(chǔ)上，可以看作是局部敏感性在全局范圍內(nèi)的擴(kuò)展。對(duì)于涉及參數(shù)較多的高維模型而言，DGSM進(jìn)行數(shù)值評(píng)估所需的計(jì)算時(shí)間比廣泛應(yīng)用的Sobol敏感度指數(shù)的計(jì)算時(shí)間低許多個(gè)數(shù)量級(jí)。在本研究中，我們將計(jì)算所得的各參數(shù)敏感性指數(shù)占所有參數(shù)敏感性指數(shù)的比例作為該方法的敏感性度量。各參數(shù)敏感性指數(shù)的確定公式如下：

STi=vi/π2D

(4)

4 結(jié)果與討論

4.1 定性敏感性分析

根據(jù)第3.2節(jié)中關(guān)于Morris取樣方法的介紹，取300個(gè)參數(shù)樣本對(duì)BTOPMC進(jìn)行定性參數(shù)篩選，結(jié)果見(jiàn)圖4(a)，可以發(fā)現(xiàn)X8，X9，XC，X7，X6，X4，XB，X3，X5，XF，XE，XA的μ*值(每個(gè)Xi的箱線圖中紅色線的對(duì)應(yīng)值)明顯大于0，可以確定為敏感變量，涉及的參數(shù)類型有m(m)，Srmax(m)，Sbar0(m)和n0。此外，LH抽樣也可用于MOAT進(jìn)行定性SA分析，樣本數(shù)量同樣取300，分析結(jié)果見(jiàn)圖4(b)，我們可以觀察到的敏感變量有X9，X8，X6，XE，X7，X3，XF，涉及到參數(shù)的類型有m(m)，Srmax(m)和n0。綜合分析，兩種取樣方法所對(duì)應(yīng)的MOAT結(jié)果是具有一定程度上的一致性的，由于采樣不確定性，它們之間所存在的差異是合乎情理的。

根據(jù)MOAT方法雖然可能將不重要參數(shù)判斷為重要參數(shù)，但不會(huì)出現(xiàn)將重要參數(shù)判為不重要參數(shù)的錯(cuò)誤[18]的這一分析特征，將m(m)，Srmax(m)和n0確定為BTOPMC的敏感參數(shù)類型，其中Srmax(m)中所包含的X6，X7，X8，X9和m(m)中所包含的X3，X4，X5均可被確定為需要被率定的敏感變量，n0中XE和XF的敏感性較為突出，同樣被認(rèn)定為敏感變量。

圖4 BTOPMC模型參數(shù)敏感性的定性分析結(jié)果Fig.4 Qualitative sensitively analysis results of the parameters in the BTOPMC model

4.2 定量敏感性分析

為了驗(yàn)證MOAT方法對(duì)敏感參數(shù)的篩選結(jié)果，并量化每個(gè)敏感性參數(shù)的貢獻(xiàn)率，我們利用DGSM這一種定量SA方法對(duì)BTOPMC的參數(shù)敏感性進(jìn)行量化分析。本研究將LH作為DGSM的取樣方法。由于定量SA方法通常比定性SA方法需要更多的樣本，所以本文將樣本數(shù)量取為500，分析結(jié)果分別見(jiàn)圖5。可明顯看出，X9，X8，X7，X6的敏感性較強(qiáng)，此外，X3，X4.X5,XE，XF的總效應(yīng)也很突出，對(duì)應(yīng)到參數(shù)類別，我們可以將Srmax(m)，m(m)和n0作為BTOPMC模型的敏感參數(shù)類型。此外，雖然XA，XB，XC，XD對(duì)于SA結(jié)果有一定程度上的影響，但相比較而言影響較小，在此不將其作為敏感性參數(shù)考慮。

比較由MOAT方法所得的定性分析結(jié)果和由DGSM方法所得的定量分析結(jié)果，可以發(fā)現(xiàn)所篩選的BTOPMC的敏感參數(shù)基本一致。由于采樣方法及分析方法的不同，參數(shù)敏感程度的大小排序存在差異也是合乎情理的[20]。

圖5 BTOPMC模型參數(shù)敏感性的定量分析結(jié)果Fig.5 Quantitative sensitively analysis results of the parameters in the BTOPMC model

4.3 SA結(jié)果的物理解釋

BTOPMC中，T0和m是反映流域土壤類型影響的參數(shù)類型。作為T0的衰減因子，m的變化將導(dǎo)致T0明顯變化[10]，所以，當(dāng)m和T0發(fā)生相同的變化時(shí)，由m所引起的模型性能的變化是更加明顯的。也就是說(shuō)，在BTOPMC中m比T0更敏感。從表1可知，m中的X3，X4和X5分別是T0中X0，X1和X2的衰減因子，相應(yīng)地，X3，X4和X5更為敏感。

Srmax主要反映了植被和土地利用對(duì)模擬結(jié)果的影響情況。在根區(qū)的儲(chǔ)水值達(dá)到Srmax之前，水將在重力的作用下滲入不飽和區(qū)域。富士川流域面積的88%被植被覆蓋，Srmax變化對(duì)BTOPMC在該流域的模擬結(jié)果影響較大是合理的。此外，與深根植被區(qū)相比，淺根植被區(qū)的Srmax對(duì)地表徑流的影響更大，而灌溉區(qū)和不透水區(qū)的Srmax對(duì)模型模擬結(jié)果的影響最大，對(duì)應(yīng)到相應(yīng)參數(shù)，X9和X8比X7更敏感，X6的敏感性則相對(duì)較小。

n0即河道對(duì)水流的摩擦阻力，是BTOPMC匯流模型中唯一需要被率定的參數(shù)，受多種因素影響，如河岸或河床材料、河道結(jié)構(gòu)、河岸不規(guī)則程度、植被覆蓋情況等。因此，匯流模型通常對(duì)n0的變化比較敏感。此外，黏土、砂土和粉土的n0基本相同，且都小于河網(wǎng)和建筑區(qū)，因此XE和XF的變化對(duì)模型性能的影響要大于XG，XH和XI。

Sbar0是反映流域地形對(duì)模型模擬結(jié)果影響情況的參數(shù)類型，需要針對(duì)每個(gè)子流域進(jìn)行率定。與針對(duì)每個(gè)網(wǎng)格進(jìn)行率定的前4種參數(shù)類型相比，大面積均質(zhì)化可以解釋其在BTOPMC中的相對(duì)不敏感性。

當(dāng)然，由于SA方法和樣本選取方法的局限性、觀測(cè)誤差、模型本身的結(jié)構(gòu)不完善等原因，使得并非所有篩選結(jié)果都可以有精確的物理解釋[21]。

5 結(jié) 論

在本次研究中，采用MOAT這一定性SA方法和DGSM這一定量SA方法分別對(duì)BTOPMC模型的參數(shù)敏感性進(jìn)行定性和定量分析，將Srmax(m)，m(m)和n0確定為BTOPMC模型用于富士川流域時(shí)的敏感參數(shù)。

本文的敏感參數(shù)篩選結(jié)果與模型參數(shù)的物理解釋一致，這也證明了采用全局SA方法尋找復(fù)雜模型中重要參數(shù)的可行性。

雖然在富士川流域通過(guò)三種方法篩選出了敏感參數(shù)，但還需要在此基礎(chǔ)上對(duì)該模型參數(shù)的“異參同效”現(xiàn)象進(jìn)行分析，并進(jìn)行參數(shù)率定以進(jìn)一步明確該研究結(jié)果對(duì)BTOPMC參數(shù)優(yōu)化效率的影響情況。此外，通過(guò)減少參數(shù)個(gè)數(shù)來(lái)降低BTOPMC模型參數(shù)率定結(jié)果的不確定性和提高參數(shù)率定效率的普遍性需要在更多具有不同特征的流域上得到進(jìn)一步研究。