劉 昭，趙樹(shù)旗，劉培斌，王利軍，高曉薇

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

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

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

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

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

1.1 研究區(qū)概況

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

圖1 媯水河流域

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

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

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

2 研究方法

2.1 SWA模擬

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

2.2 SUFI-2敏感性分析方法

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

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

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

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

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

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

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

(1)

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

bR2的計(jì)算公式為

(2)

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

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

(3)

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

2.4 不確定性分析方法

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

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

3 結(jié)果分析

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

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

圖2 率定期觀測(cè)值與最佳模擬值徑流過(guò)程線

3.2 模型結(jié)果分析和評(píng)價(jià)

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

圖3 驗(yàn)證期觀測(cè)值與模擬值徑流過(guò)程線

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

3.3 不確定性分析

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

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

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

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

4 結(jié) 語(yǔ)

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