基于Bootstrap法的水文模型參數(shù)不確定分析
——以伊通河流域為例

宋文博，盧文喜，董海彪，張 正，肖傳寧(.吉林大學環(huán)境與資源學院，長春 3006；.松遼流域水環(huán)境科學研究所，長春 3006)

流域水文循環(huán)過程及其演化規(guī)律通?？梢詰昧饔蛩哪Ｐ瓦M行模擬。為了將復雜的水文過程概念化和抽象化，流域水文模型通常包含大量數(shù)學方程和參數(shù)用以描述實際流域水文過程和流域特征[1,2]。在確定模型參數(shù)的過程中往往由于主觀和客觀的因素影響而存在著較大的不確定性，這會導致在模擬過程中將會出現(xiàn)“失真”現(xiàn)象，進而導致模型模擬結果的不確定性。水文模型中，參數(shù)的不確定性主要來源于2方面。首先，由于下墊面的空間異質(zhì)性將使參數(shù)值存在空間異質(zhì)性，這進而導致了參數(shù)的不確定性[3,4]；其次，模型中的參數(shù)值常使用觀測數(shù)據(jù)進行率定來確定，但是模型率定獲取參數(shù)值受到諸多因素的影響，例如參數(shù)之間的相互影響、參數(shù)的敏感性差異、時間和空間尺度以及模型殘差的統(tǒng)計特征[5]，這些因素使模型參數(shù)值與其真實值“所謂”的等價，導致不確定性問題始終貫穿于整個模型模擬的過程中[6]。因此，需要對參數(shù)的不確定性進行現(xiàn)實評估，并正確評價參數(shù)的不確定性對模型模擬的不確定性影響程度以及模型模擬結果的不確定性范圍。

國內(nèi)外學者在模型參數(shù)不確定分析方面開展了大量研究，很多方法被應用于分析參數(shù)的不確定性問題，常用的方法如GLUE方法[9]、馬爾科夫鏈蒙特卡洛方法[10]、貝葉斯方法[11]等。Bootstrap法是一種非參數(shù)的統(tǒng)計方法，有別于其他不確定分析方法，它從實際數(shù)據(jù)出發(fā)，不需要對觀察數(shù)據(jù)的分布作任何假設，簡單易行，在許多領域都具有廣泛的應用。本文運用SWAT模型對伊通河流域的水文過程進行模擬，通過參數(shù)敏感性分析，得到模型敏感性參數(shù)。與以往不同的是，本文嘗試應用Bootstrap法獲取模型敏感參數(shù)的邊緣分布，以95%置信區(qū)間對其不確定性范圍進行描述，并結合變異系數(shù)對參數(shù)的不確定性程度進行評估。在此基礎上，進一步分析了模型參數(shù)的不確定性對模型模擬結果不確定性的貢獻。研究成果可以為Bootstrap法運用于模型參數(shù)不確定性分析以及SWAT模型在伊通河流域水文模擬應用方面研究提供借鑒與參考。

1 材料與方法

1.1 Bootstrap法

Bootstrap法是由美國Stanford大學統(tǒng)計學教授Efron最早提出的，它是利用再抽樣技術來評估不確定性的方法[12]。Bootstrap法通過計算機模擬來替代對偏差、方差和其他統(tǒng)計量的復雜而不精確的近似方法[13,14]。這種方法不必對未知分布做出人為假設，通過對原始數(shù)據(jù)進行再抽樣來獲取未知分布。因此Bootstrap法是一種對獨立同分布的中小樣本進行統(tǒng)計推斷，并能在統(tǒng)計信息不充分的條件下提高推斷的方法[15]。然而統(tǒng)計數(shù)據(jù)在很多情況下并不能滿足獨立同分布，而是存在一定的相依結構，如時間序列數(shù)據(jù)等，在這種情況下，如果繼續(xù)應用獨立同分布的Bootstrap，則很可能會導致失敗。由此衍生出一系列帶有相依結構數(shù)據(jù)的Bootstrap方法，如區(qū)間Bootstrap法和基于模型的Bootstrap法[16]。

(1)在ε(t)中進行有放回抽樣，獲得新的殘差系列ε*(t)。

(4)重復步驟(1)-(3)，本研究進行再抽樣1 000次。

(5)應用Bootstrap法獲取 的1 000個Bootstrap估計量(θ*i1,…,θ*i 1 000)，確定θi在置信水平α的置信區(qū)間[θ*i 1 000(α/2),…,θ*i 1 000(α/2)]。

1.2 SWAT模型

SWAT(Soil and Water Assessment Tool)模型是20世紀90年代由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心(USDA-ARS)在SWRRB[19]基礎上開發(fā)的基于過程的連續(xù)性分布式流域水文模型，它主要用于氣候變化和管理措施對流域水文過程和非點源污染物的影響研究[20]。SWAT模型結構上分為水文過程、土壤侵蝕和污染負荷3個子模型[21]。

其中水文過程子模型采用水量平衡原理，表達式為：

(1)

式中：SWt為土壤最終含水量，mm；SW0為土壤前期含水量，mm；Rday為第i天降水量，mm；Qsurf為第i天的地表徑流，mm；Ea為第i天的蒸發(fā)量，mm；Wseep為第i天存在于土壤剖面地層的滲透量和測流量，mm；Qgw為第i天地下水含量，mm。

在描述地表降水徑流時模型采用SCS曲線數(shù)法，其表達式為：

(2)

式中：P為一次性降水總量，mm；Q為徑流量，m3/s；Iα為產(chǎn)生地表徑流前的降雨損失，mm；F為產(chǎn)生地表徑流后的降雨損失，mm；S為流域當時的可能最大滯留量，mm。

1.3 研究區(qū)概況

伊通河流域面積8 840 km2，地處東經(jīng)124°30′～125°45′，北緯43°0′～44°45′。伊通河發(fā)源于吉林省伊通縣，全長342.5 km，多年平均徑流總量5.4 億m3，由南向北貫穿長春市城區(qū)，流經(jīng)伊通縣、長春市、德惠市、農(nóng)安縣，在農(nóng)安縣境內(nèi)的靠山鎮(zhèn)附近匯入飲馬河后流入第二松花江，流域地理位置見圖1。研究區(qū)屬于中溫帶大陸性半濕潤季風氣候類型，年平均氣溫4.8 ℃，年平均降水量583.5 mm，多集中在6-9月份，占全年降水量的80%以上。研究區(qū)包括旱田、水田、林地、草地、水域、城鎮(zhèn)6種土地利用類型，以旱田為主。土壤類型有暗棕壤、白漿土、黑土、黑鈣土、草甸土、水稻土、新積土、沼澤土、風沙土9種。

圖1 伊通河流域地理位置與相關測站分布Fig.1 Watershed location and meteorological stations distribution

1.4 模型數(shù)據(jù)準備

SWAT模型所需數(shù)據(jù)由地形、土壤、土地利用、水文和氣象數(shù)據(jù)組成。本研究選用的伊通河流域地形數(shù)據(jù)來源于分辨率為30 m數(shù)字高程模型(DEM)柵格數(shù)據(jù)；土地利用數(shù)據(jù)為2008年伊通河流域1∶50萬土地利用類型圖；土壤數(shù)據(jù)為吉林省1∶50萬土壤類型圖矢量化獲得；氣象數(shù)據(jù)為2008-2012年的長春、農(nóng)安、雙陽、煙筒山、四平5個氣象站的逐日氣象資料(包括降水量、最高和最低氣溫、相對濕度、太陽輻射、風速)；本文采用伊通水文站2009-2012年月徑流數(shù)據(jù)進行模型參數(shù)的率定與檢驗。該水文系列包含了流域豐(2010年)、平(2011、2012年)、枯(2009年)水期的徑流，具有代表性和可靠性。其中應用2009-2010年月徑流數(shù)據(jù)進行模型參數(shù)率定，2011-2012年數(shù)據(jù)進行驗證，數(shù)據(jù)來源及說明見表1。

表1 數(shù)據(jù)來源及說明Tab.1 Source of data and instructions

1.5 子流域劃分及水文響應單元(HRU)生成

SWAT模型基于流域DEM提取水系，通過設置子流域集水面積閾值和流域出口位置進行子流域的劃分和子流域的參數(shù)計算。本文設置子流域集水面積閾值為15 km2，并將伊通河流域劃分為31個子流域，見圖1。在子流域劃分的基礎上，根據(jù)不同的土地利用類型和土壤類型的組合，在每一個子流域內(nèi)進一步劃分水文響應單元(Hydrologic Response Unit，HRU)。HRU的劃分是根據(jù)不同土地利用類型和土壤類型組合將子流域劃分為多個HRUs，為提高計算效率設置土地利用和土壤面積閾值，忽略低于閾值面積的土地利用類型和土壤類型。根據(jù)實際情況，本文設置土地利用面積閾值為4%，土壤面積閾值為10%，共劃分為495個HURs。

1.6 模型參數(shù)的選擇

SWAT模型中包含眾多參數(shù)，但并非所有參數(shù)對模型的輸出結果都具有顯著影響。因此本研究采用LH-OAT(Latin Hypercube-One foctor At a Time)方法進行參數(shù)敏感性分析，同時在參考相關研究經(jīng)驗的基礎上，篩選出CN2、ESCO、ALPHA_BF、SOL_AWC、SMTMP 5個對模型輸出較為敏感的參數(shù)進行不確定分析。由于SWAT模型為分布式水文模型，某些參數(shù)在不同的子流域上具有不同的參數(shù)取值，并需要在各自子流域上單獨率定。因此為簡化參數(shù)的率定過程并提高工作效率，同時保持參數(shù)在空間上的分布性，參數(shù)以a_、v_集成參數(shù)[22](aggregate parameters)的形式進行表示，其中a_表示在給定值域范圍內(nèi)，原值基礎上加減這個范圍內(nèi)的值，v_表示用給定值域范圍內(nèi)的值替換原值，各參數(shù)表示形式與取值范圍見表2。采用SWAT模型自動率定模塊中最小二乘法進行模型參數(shù)的率定和驗證，并應用決定系數(shù)R2、Nash-Sutteliffe效率系數(shù)Ens和相對誤差Re對模擬結果進行評價，當徑流模擬值與實測值相對誤差在±20%內(nèi)，Ens>0.5，R2>0.6，則符合SWAT模型模擬要求，可用于流域模擬。

表2 模型參數(shù)及取值范圍Tab.2 Parameter selection and its scope

2 結果與分析

2.1 模型模擬結果

根據(jù)前述模型評價標準，模型月徑流量模擬結果良好，見表3和圖2。其中率定期R2為0.93，Ens為0.89；驗證期R2為0.82，Ens為0.80。相對誤差Re率定期與驗證期分別為7.21%和10.38%。模型模擬結果符合SWAT模型模擬要求。

表3 模型率定與驗證結果Tab.3 model calibration and verification

2.2 參數(shù)不確定分析

應用基于模型的Bootstrap法生成1 000組SWAT模型敏感參數(shù)的Bootstrap估計值，分別對每個參數(shù)1000個自助估計值進行有序排列，并分別以2.5%和97.5%的閾值作為上限和下限，獲取敏感參數(shù)95%的置信區(qū)間。5個敏感參數(shù)的95%置信區(qū)間以及邊緣分布見表4和圖3。

表4 參數(shù)95%置信區(qū)間及變異系數(shù)Tab.4

圖2 月徑流量擬合曲線Fig.2 Monthly runoff fitting curve

圖3 參數(shù)邊緣分布頻率直方圖Fig.3 Parameters of the marginal distribution frequency histograms

從圖3中可以看出，a_CN2.mgt、v_ESCO.hru、v_ALPHA_BF.gw、v_SMTMP.bsn近似服從正態(tài)分布，而a_SOL_AWC.sol呈負偏態(tài)分布。表4中給出了每一個參數(shù)的不確定性范圍，CN2在SWAT模型中表示徑流曲線數(shù)，是土壤滲透性、土地利用和土壤前期含水量的綜合反映，是影響徑流量的重要參數(shù)，其不確定區(qū)間為[1.176,2.634]。ESCO在SWAT模型中表示土壤蒸發(fā)補償系數(shù)，該參數(shù)值越小說明土壤蒸發(fā)量越大，其參數(shù)不確定區(qū)間為[0.343,0.392]。ALPHA_BF為基流系數(shù)，反映了地下徑流對降水補給的響應，不確定性區(qū)間為[0.217,0.353]。SMTMP是雪融最低氣溫，對寒冷季節(jié)融雪具有重要影響，進而影響地表徑流，它的不確定性范圍為[1.332,2.017]。SOL_AWC為土壤可利用水量，其不確定區(qū)間為[0.024,0.080]。從中可以發(fā)現(xiàn)在進行bootstrap不確定性分析之后，這些參數(shù)的取值范圍都較初始的取值范圍小。

進一步探討模型參數(shù)的不確定性，本研究采用變異系數(shù)Cv值對各個參數(shù)的不確定性進行評估。變異系數(shù)為標準差與均值的比，是用來描述數(shù)據(jù)離散程度的指標，變異系數(shù)越大，表明數(shù)據(jù)離散程度越大，即不確定性程度越高[23,24]。由表4可知，a_CN2.mgt、v_ESCO.hru相比于其他參數(shù)具有較高不確定性，變異系數(shù)分別為21.03%和16.43%，而v_ALPHA_BF.gw、v_SMTMP.bsn、a_SOL_AWC.sol參數(shù)不確定性程度較低，變異系數(shù)在10%左右。參數(shù)按不確定性程度高低進行排序為：a_CN2.mgt>v_ESCO.hru>v_ALPHA_BF.gw>v_SMTMP.bsn>a_SOL_AWC.sol。

2.3 模型模擬不確定分析

參數(shù)的不確定性會帶來模型模擬結果的不確定性，因此為探究參數(shù)的不確定性對模型模擬的不確定性的影響程度，本研究計算了月徑流量模擬的95%置信區(qū)間。若落入該置信區(qū)間的觀測值越多，則表明參數(shù)的不確定性對模型模擬的不確定性貢獻越大，若置信區(qū)間內(nèi)包含的觀測值很少，則說明參數(shù)的不確定性對模型模擬的不確定性貢獻較小。圖4是月徑流模擬的不確定性范圍，對結果進行統(tǒng)計發(fā)現(xiàn)58%的觀測流量落入到了這個置信區(qū)間內(nèi)。這表明參數(shù)的不確定性對模型模擬結果的不確定性具有較大貢獻，在所有不確定性因素中占據(jù)較大比例。其次，從圖4中可以看出，徑流量大的月份具有較大不確定性范圍，相反徑流量小的月份不確定性范圍較小，說明模擬結果的不確定性與徑流量的大小有關，這與程曉光[25]等人得到的結論一致。同時發(fā)現(xiàn)在每年11月份到第2年3月份期間僅8.3%觀測流量落入該不確定區(qū)間內(nèi)，表明這段時間內(nèi)模型參數(shù)的不確定性對模型模擬結果的不確定性貢獻較小。這段時期伊通河流域正值降雪期，流域內(nèi)徑流主要通過融雪產(chǎn)生。SWAT模型融雪模塊是利用溫度因子方法來實現(xiàn)[26]，然而僅由溫度不能作為模擬融雪過程的主導因素[27]，這會給模擬結果帶來較高的不確定性。因此研究區(qū)降雪期模型模擬的不確定性可能主要來自于模型結構的不確定性影響。

圖4 研究區(qū)2009-2012年月徑流模擬的不確定性范圍Fig.4 Uncertainty range of monthly runoff simulation in study area from 2009 to 2012

3 結論與建議

(1)Bootstrap法作為一種非參數(shù)統(tǒng)計方法，可以在不對模型參數(shù)作任何分布假設情況下，通過實測數(shù)據(jù)來獲取模型參數(shù)的未知分布，為水文模型參數(shù)的不確定分析提供了新的思路。

(2)本文應用Bootstrap法獲得了模型5個敏感性參數(shù)的邊緣分布以及這些參數(shù)的不確定性范圍，結果表明在進行Bootstrap不確定性分析之后，參數(shù)的取值范圍都較初始的取值范圍小，其中SOL_AWC呈負偏態(tài)分布，其余4個參數(shù)近似服從正態(tài)分布，并通過對參數(shù)的變異系數(shù)計算，表明CN2和ESCO 2個參數(shù)具有較高的不確定性。研究成果說明Bootstrap方法可以運用于水文模型參數(shù)的不確定分析研究。

(3)進一步分析了參數(shù)的不確定性對模型模擬結果的不確定性貢獻，結果表明參數(shù)的不確定性對模型模擬的不確定性具有較大貢獻。此外計算結果還顯示研究區(qū)降雪期模型模擬的不確定性受參數(shù)的不確定性貢獻較小，可能主要受模型結構的不確定性影響，建議進一步對模型結構的不確定性進行分析。

□

[1] 衛(wèi)曉婧．流域水文模型不確定性研究進展[J]．中國水運，2008，8(7)：165-167．

[2] 衛(wèi)曉婧，熊立華．改進的GLUE方法在水文模型不確定性研究中的應用[J]．水利水電快報，2008，29(6)：23-25．

[3] 杜新忠．流域水文模型的不確定性分析[D]．長沙：長沙理工大學，2011．

[4] 林凱榮，陳曉宏，江 濤，等．基于Beta-PERT分布的水文模型參數(shù)不確定性分析[J]．中山大學學報(自然科學版)，2010，49(4)：139-143．

[5] 李向陽．水文模型參數(shù)優(yōu)選及不確定性分析方法研究[D]．大連：大連理工大學，2005．

[6] 熊立華，萬 民．水文模型兩種不確定性研究方法的比較[J]．武漢大學學報(工學版)，2009，42(2)：137-142．

[9] 李 勝，梁忠民．GLUE方法分析新安江模型參數(shù)不確定性的應用研究[J]．東北水利水電，2006，24(2)：31-33，47．

[10] Kanso A，Chebbo G，Tassin B．Application of MCMC-GSA model calibration method to urban runoff quality modeling[J]．Reliability Engineering & System Safety，2006，91：1 398-1 405．

[11] K Engeland，Xu C Y ，Gottschalk L．Assessing uncertainties in a conceptual water balance model using Bayesian methodology[J]．Hydrological Sciences Journal，2005，50：45-63．

[12] 劉 偉，龍 瓊，陳 芳，等．Bootstrap方法的幾點思考[J]．飛行器測控學報，2007，26(5)：78-81．

[13] 鄧海軍，查亞兵．自助法中若干問題研究及其在命中精度評估中的應用[J]．2005，24(1)：59-63．

[14] Jhun M，Jeong H C．Application of Bootstrap Methods for Categorical Data Analysis[J]．Computational Statistics & Data Analysis，2000,(35)：83-91．

[15] 黃 瑋，馮蘊雯，呂震宙．基于Bootstrap方法的小子樣試驗評估方法研究[J]．機械科學與技術，2006，25(1)：31-35．

[16] 曾 翀．Bootstrap方法的歷史發(fā)展及其在金融風險管理和藥物動力學中的應用[D]．武漢：華中科技大學，2009．

[17] Domingo M，Leandro P，Laureano S．Bootstrap confidence regions in multinomial sampling[J]．Applied Mathematics and Computation，2004,(155)：295-315．

[18] 張金槐．精度評估中的自助方法[J]．飛行器測控學報，2000，19(3)：61-66．

[19] Arnold，J G，Williams，J R，Nicks，A D，et al. SWRRB: a basin scale simulation model for soil and water resources management[M]. College Station：Tex as A&M Press，1990.

[20] 尹雄銳，章光新，楊 帆，等. 東北半干旱地區(qū)流域分布式水文模擬——以洮兒河流域為例[J]. 吉林大學學報(地球科學版)，2011，41(1) ：137-144.

[21] 湯 潔，楊 巍，李昭陽，等. 遼河大伙房水庫匯水區(qū)農(nóng)業(yè)非點源污染入庫模擬[J]. 吉林大學學報(地球科學版)，2012，42(5)：1 462-1 468，1 476.

[22] Yang J，Reichert P，Abbaspour K C，et al．Hydrological modelling of the Chaohe basin in China：statistical model formulation and Bayesian inference[J]．Journal of Hydrological，2007,(340):167-182．

[23] Francos A, Elorza F J, Bouraoui F, et al. Sensitivity analysis of distributed environmental simulation models: understanding the model behaviour in hydrological studies at the catchment scale [J]. Reliability Engineering and System Safety, 2003,79:205-218.

[24] Zhang zheng，Lu wenxi，Chu haibo，et al. Uncertainty analysis of hydrological model parameters based on the bootstrap method：a case study of the SWAT model applied to the Dongliao River Watershed，Jinlin Province，Northeastern China[J]. Science China，2014，57(1)：219-229.

[25] 程曉光，張 靜，宮輝力. 半干旱半濕潤地區(qū)HSPF模型水文模擬及參數(shù)不確定性研究[J]. 環(huán)境科學學報，2014，34(12)：3 179-3 187.

[26] 孟現(xiàn)勇，吉曉楠，劉志輝，等. SWAT模型融雪模塊的改進與應用研究[J]. 自然資源學報，2014，29(3)：528-539.

[27] Marks D，Kimball J，Tingey D，et al. The sensitivity of snowmelt processes to climate conditions and forest cover during rain-on-snow：a case study of the 1996 Pacific Northwest flood [J]. Hydrological Processes，1998，12(10/11)：1 569-1 587.