媯水河流域水文模擬及參數(shù)不確定性分析

劉 昭，趙樹旗，劉培斌，王利軍，高曉薇

(1.北京工業(yè)大學(xué)建筑工程學(xué)院，北京 100124；2.北京市水利規(guī)劃設(shè)計研究院，北京 100048)

媯水河是北京市延慶區(qū)內(nèi)的主干河流，其干流及各支流是城鄉(xiāng)居民生活用水和工農(nóng)業(yè)生產(chǎn)用水的重要來源。近年來，由于人口增加、經(jīng)濟(jì)發(fā)展與人類活動增強(qiáng)以及持續(xù)的干旱氣候條件，使得水資源日益減少，水資源短缺嚴(yán)重，供需矛盾突出[1-2]。

水是自然界的一種重要載體。由于水的輸送作用，泥沙、污染物等隨地表徑流進(jìn)入河流。對流域地表徑流的模擬也是非點源污染研究、水質(zhì)水量科學(xué)管理、產(chǎn)沙量計算的重要前提。SWAT模型(Soil and Water Assessment Tool)是由美國農(nóng)業(yè)部(USDA)農(nóng)業(yè)研究中心(ARS)開發(fā)的基于物理過程，并能夠模擬不同土地利用和多種農(nóng)業(yè)管理措施對流域水、泥沙、營養(yǎng)物質(zhì)、殺蟲劑等輸送遷移影響的分布式流域水文模型[3]。

本文選取SWAT模型進(jìn)行媯水河流域地表月徑流量模擬，根據(jù)現(xiàn)有數(shù)據(jù)選取2014年、2015年為模擬期，2017年為驗證期。運用SWAT-CUP中的SUFI-2算法分析模型參數(shù)的敏感性、對所選參數(shù)進(jìn)行率定和驗證以及分析各參數(shù)的不確定性。通過以上分析，以期為研究SWAT模型在媯水河流域的適用程度和模型參數(shù)的敏感性及不確定性提供參考、為延慶區(qū)水資源管理及世園會、冬奧會水質(zhì)水量安全保障提供支持。

1 研究區(qū)概況與基本數(shù)據(jù)

1.1 研究區(qū)概況

媯水河流域見圖1。媯水河發(fā)源于延慶區(qū)大吉祥村南，屬海河流域永定河水系。河流穿過延慶山間盆地匯入官廳水庫，其自東向西的流向在我國西高東低的地形下顯得尤為特別。媯水河橫亙在八達(dá)嶺西北方向，譽稱“首都西北大門”。本文研究區(qū)選自媯水河谷家營斷面以上流域，控制流域面積732 km2。上游沿途有9條主要支流匯入；支流中僅古城河流量較大常年有基流，其余支流流量較小在枯水期時常斷流，屬季節(jié)性河流。東大橋水文站(東經(jīng)116°00′00″，北緯40°25′30″)為媯水河唯一的水文站，等級為區(qū)域代表站。媯水河流域近20年的年降雨量大致在200～450 mm之間，流域降雨較為缺乏，屬河流的枯水年段。

圖1 媯水河流域

表1 SWAT模型數(shù)據(jù)信息

1.2 基本數(shù)據(jù)

本研究所用的數(shù)據(jù)包括：研究區(qū)4個氣象站點(見圖1)的日數(shù)據(jù)以及東大橋水文站的月徑流量數(shù)據(jù)，DEM，2015年土地利用數(shù)據(jù)，土壤數(shù)據(jù)集及土壤屬性數(shù)據(jù)，具體如表1所示。

2 研究方法

2.1 SWA模擬

SWAT模型是流域尺度的分布式動態(tài)水文模型[4]。該模型能夠利用GIS和RS提供的空間信息，模擬流域中多種不同的水文物理過程，是流域尺度上的動態(tài)模擬模型[5-7]。本次流域水文模擬可分為以下幾個步驟：①加載媯水河流域的DEM數(shù)據(jù)、添加點源、水庫入水口等，同時河流將被生成，確定流域的出水口后可通過SWAT劃分各個子流域。②添加流域土地利用圖、土壤類型圖和坡度的疊加分析，其后計算水文響應(yīng)單元。③加載氣象數(shù)據(jù)以及點源污染、水庫入水口的相關(guān)參數(shù)。④啟用SWAT模擬運算。完成后需要分析模型參數(shù)的敏感性、率定和檢驗參數(shù)最后將分析參數(shù)的不確定性。從而使流域內(nèi)水文狀況得到更加準(zhǔn)確的模擬。

2.2 SUFI-2敏感性分析方法

SUFI-2是利用拉丁超立方體抽樣法獲得參數(shù)值帶入模型進(jìn)行模擬[8]。在SUFI-2中，有兩類共3種方法得到參數(shù)的敏感性[9]。3種方法：①敏感性分析，通過每一次只改變一個參數(shù)的敏感性進(jìn)行計算，其他參數(shù)保持不變。②全局敏感性分析，這種方法是在模型參數(shù)率定的過程中計算下次率定時所需參數(shù)的敏感性，可在一次迭代后執(zhí)行。③分析dotty plots圖，其展示了全部參數(shù)的散點圖，目的是通過抽樣點分布圖像的展示來反映參數(shù)的敏感性。通過顯示每次模擬參數(shù)值在其取值范圍和NSE值范圍所構(gòu)成的坐標(biāo)系內(nèi)的分布，把參數(shù)的敏感性區(qū)間縮小到最符合要求的區(qū)間。

本文所用SUFI-2敏感性分析步驟：①結(jié)合已有文獻(xiàn)的經(jīng)驗論述以及參照SWAT-CUP使用手冊推薦選擇相對較大的參數(shù)范圍；②按照SUFI-2算法要求，設(shè)定參數(shù)范圍等文件并選擇目標(biāo)函數(shù)；③對參數(shù)的推薦值或推薦范圍進(jìn)行重復(fù)迭代，當(dāng)?shù)Y(jié)果符合要求時停止。

2.3 結(jié)果評價標(biāo)準(zhǔn)

本文選用確定性系數(shù)R2及其修正參數(shù)bR2、Nash-Sutcliffe模型效率系數(shù)NSE作為模型參數(shù)的率

表2 參數(shù)敏感性及參數(shù)率定分析結(jié)果

定及驗證過程中的結(jié)果評價標(biāo)準(zhǔn)。以上系數(shù)在SWAT-CUP軟件中通過自動求參得到：①對于R2，R2越靠近1說明模擬值與觀測值的變化過程具有較高相似性；R2偏離1程度越大，說明觀測值與模擬值的變化過程具有較低相似性。②對于bR2，通過其取值可分析觀測值與模擬值在數(shù)值上的大小差異及其變化情況。③對于NSE，NSE值越靠近1說明模擬值在數(shù)值大小上與觀測值具有較高相似性；NSE偏離1程度越大，說明模擬值與觀測值在數(shù)值大小上具有較低相似性。若NSE出現(xiàn)負(fù)值，說明使用模擬值的可信度低于使用觀測值的算術(shù)平均值。

確定性系數(shù)R2的計算公式為

(1)

式中，Qm·i為徑流觀測值；Qm為徑流觀測值平均值；Qs·i為徑流模擬值；Qs為徑流模擬值平均值。

bR2的計算公式為

(2)

式中，b為回歸系數(shù);R2為確定性系數(shù)。

效率系數(shù)NSE的計算公式為

(3)

式中，Qo·i為觀測值；Qp·i為模擬值；Qavg為觀測值平均值；n為觀測值個數(shù)。

2.4 不確定性分析方法

SUFI-2算法作為參數(shù)估計的最優(yōu)化方法，考慮了輸入數(shù)據(jù)、模型結(jié)構(gòu)、參數(shù)及觀測數(shù)據(jù)等的不確定性，并將其反映在率定后的參數(shù)范圍內(nèi)。參數(shù)率定后95%置信水平上的不確定性區(qū)間(95% prediction uncer-tainty， 95PPU)包含大多數(shù)觀測數(shù)據(jù)，模擬結(jié)果總的不確定性，在2.5%(L95PPU)和97.5%(U95PPU)上的累積分布得到，通過拉丁超立方采樣計算輸出[9]。

本文使用該算法的兩個指標(biāo)來評價參數(shù)的不確定性大?。孩貾-factor(95%不確定性區(qū)間內(nèi)觀測數(shù)據(jù)所占百分比)，該指標(biāo)表示95%置信區(qū)間內(nèi)觀測數(shù)據(jù)量。P-factor的取值范圍為0～100%，值接近0時表示區(qū)間內(nèi)只包含極少觀測數(shù)據(jù)值接近100時表示區(qū)間內(nèi)包含幾乎所有的觀測數(shù)據(jù)。②R-factor(95PPU上下限的平均距離與標(biāo)準(zhǔn)偏差所占比例)，該指標(biāo)表示95%置信區(qū)間內(nèi)樣本的集中程度。R-factor的取值范圍為0～+∞。綜上可得：當(dāng)P-factor趨于1同時R-factor趨于0時，模擬結(jié)果與觀測值具有高度的一致性；故，通常以此作為標(biāo)準(zhǔn)來評價模擬結(jié)果。

3 結(jié)果分析

3.1 參數(shù)率定及敏感性分析

參數(shù)校準(zhǔn)的過程中SWAT內(nèi)置的調(diào)參功能存在耗時長、率定參數(shù)少并受調(diào)試者個人經(jīng)驗的影響等問題。所以本文通過SWAT-CUP進(jìn)行參數(shù)的自動率定。利用SUFI-2算法自動率定模型參數(shù)、分析及排序參數(shù)的敏感性、并對模擬的結(jié)果進(jìn)行分析和評價，這彌補(bǔ)了SWAT模型內(nèi)置校準(zhǔn)功能的缺陷。經(jīng)過多次迭代逐步滿足評價標(biāo)準(zhǔn)，得到模型參數(shù)率定及敏感性分析結(jié)果見圖2、表2。

圖2 率定期觀測值與最佳模擬值徑流過程線

3.2 模型結(jié)果分析和評價

本文采用2010年～2013年為模型預(yù)熱期，2014年～2015年為率定期，驗證期為2017年。經(jīng)過多輪迭代獲得模擬結(jié)果。在率定期中R2=0.65，bR2=0.56，NSE=0.61；在驗證期中R2=0.89，bR2=0.83，NSE=0.88，用前述評價標(biāo)準(zhǔn)評價模擬結(jié)果可知：R2值均大于0.6，bR2值均大于0.5，NSE值均大于0.6。迭代后結(jié)果顯示東大橋水文站月尺度模擬值與觀測值徑流過程線總體擬合結(jié)果良好，如圖2、3所示。SWAT-CUP軟件中通過自動求參所得系數(shù)值滿足評價標(biāo)準(zhǔn)，見表3。

圖3 驗證期觀測值與模擬值徑流過程線

表3 東大橋水文站月徑流量模擬結(jié)果

3.3 不確定性分析

通過利用前述方法分析模型參數(shù)的不確定性得到以下結(jié)果：①P-factor在率定期及驗證期取值分別為0.55和0.67；②R-factor在率定期及驗證期分別為0.24和0.17。③率定期部分參數(shù)500次模擬散點圖(見圖4)。不確定性分析結(jié)果見表4。

圖4 部分參數(shù)率定期500次模擬

表4 不確定性分析結(jié)果

圖4中橫坐標(biāo)代表各參數(shù)最終的取值范圍，縱坐標(biāo)代表NSE值[10]。由以上結(jié)果可得：①本次模擬率定期跟驗證期各自包含了95%置信區(qū)間55%和67%以上的觀測數(shù)據(jù)；②95%置信區(qū)間上下限的平均距離與標(biāo)準(zhǔn)偏差所占比例較小、寬度較窄、樣本的集中程度較大；③散點圖顯示抽樣點分布較為密集，取值點在0.5以上部分占絕大多數(shù)，說明模擬的置信水平較高，不確定性較小。

4 結(jié) 語

本文對媯水河流域的水文過程通過SWAT模型進(jìn)行了模擬研究，討論分析了模型參數(shù)的敏感性、對模型的參數(shù)進(jìn)行了率定及驗證并分析了模型參數(shù)的不確定性。通過以上研究總結(jié)結(jié)論如下：①本次模擬中模擬值與觀測值總體擬合程度良好，同時也存在率定及驗證年份偏少，參數(shù)擬合相對容易的缺點。②分析參數(shù)的敏感性可得：所選參數(shù)中對流域水文過程影響最大的為CN2(SCS徑流曲線數(shù))，此外ALPHA_BF(基流消退系數(shù))、GW_DELAY(地下水滯后時間)兩參數(shù)對流域水文過程也有較大影響，同時參數(shù)敏感性分析過程中也存在著SUFI-2算法在初次輸入?yún)?shù)時存在一定盲目性的缺點。③對于不確定性，不確定性只是綜合考慮，需要在后續(xù)的模型完善中深入考慮更多氣候變化與人類活動等因素影響并完善資料的收集。通過以上研究過程及結(jié)論以期為延慶區(qū)水資源管理及世園會、冬奧會水質(zhì)水量安全保障提供參考。