基于SWAT模型對(duì)洙趙新河流域徑流的模擬研究

蘆昌興， 王甲榮,2， 宮雪亮， 孫秀玲

(1.山東大學(xué) 土建與水利學(xué)院， 山東 濟(jì)南250061； 2.河海大學(xué) 水文水資源學(xué)院， 江蘇 南京 210098)

1 研究背景

水文模型可以用來(lái)探索復(fù)雜水文循環(huán)機(jī)理和過(guò)程，也可以解決許多實(shí)際水文問(wèn)題[1]。SWAT模型作為可以進(jìn)行大尺度、長(zhǎng)時(shí)期模擬分布式水文模型，是目前解決流域水文模擬、環(huán)境評(píng)價(jià)等問(wèn)題的主要手段之一[2-3]。目前SWAT模型在國(guó)內(nèi)外已大量應(yīng)用于徑流和非點(diǎn)源污染的模擬，Tamm等[4]用SWAT模型模擬土地利用和氣候變化對(duì)波羅的海東部地區(qū)水資源的影響；Yesuf等[5]利用SWAT模型對(duì)埃塞俄比亞的Maybar流域進(jìn)行模擬，結(jié)果表明SWAT模型在熱帶流域的適用性良好；趙杰等[6]模擬了烏魯木齊河流域徑流，模擬結(jié)果與觀測(cè)流量過(guò)程線擬合程度較好；喻曉等[7]以葫蘆河流域?yàn)檠芯繉?duì)象，研究了不同亞流域數(shù)量對(duì)徑流模擬的影響；賈靜[8]、劉鵬[9]、馬放等[10]利用SWAT模型分別對(duì)秦皇島地區(qū)、滇池流域、阿什河流域非點(diǎn)源污染進(jìn)行模擬，結(jié)果均表明模型適用性良好。

中國(guó)大氣同化驅(qū)動(dòng)數(shù)據(jù)集(China Meteorological Assimilation Driving Datasets for the SWAT model，CMADS)是孟現(xiàn)勇[11]建立的應(yīng)用于SWAT模型的氣象驅(qū)動(dòng)數(shù)據(jù)集，與國(guó)內(nèi)外NCAR/DOE、JRA-25及Princeton等數(shù)據(jù)集相比，CMADS數(shù)據(jù)可更好的對(duì)我國(guó)大尺度地表分量時(shí)空演變規(guī)律進(jìn)行細(xì)致化分析[12]。由于CMADS數(shù)據(jù)建立的時(shí)間較短，目前國(guó)內(nèi)將CMADS作為SWAT模型氣象驅(qū)動(dòng)數(shù)據(jù)集的研究較少，孟現(xiàn)勇等[12-14]在分別在精博河流域、黑河流域、瑪納斯流域利用CMADS數(shù)據(jù)集驅(qū)動(dòng)SWAT模型，其結(jié)果優(yōu)于使用傳統(tǒng)氣象站的模擬結(jié)果。

南四湖作為南水北調(diào)東線工程中重要的調(diào)蓄場(chǎng)所，其水量水質(zhì)直接影響東線工程運(yùn)行的成敗。洙趙新河是南四湖流域具有代表性的入湖河流，是魯西南跨菏澤、濟(jì)寧兩市的大型防洪、排澇、灌溉綜合性河道，也是山東省15條骨干河道之一，可作為研究湖西平原區(qū)的典型流域。本文選擇洙趙新河為研究區(qū)，引入CMADS作為SWAT模型氣象驅(qū)動(dòng)數(shù)據(jù)集，構(gòu)建洙趙新河流域日尺度分布式水文模型，結(jié)合SWAT-CUP模型對(duì)參數(shù)進(jìn)行率定、敏感性以及不確定性分析，討論SWAT分布式水文模型在洙趙新河流域徑流模擬的適用性，以期為洙趙新河乃至整個(gè)南四湖水資源綜合管理與決策提供科學(xué)依據(jù)。

2 研究區(qū)概況

洙趙新河流域流經(jīng)菏澤、濟(jì)寧兩市，是魯西南地區(qū)具有防洪、排澇和灌溉功能的大型綜合性河道，其地理位置東臨南陽(yáng)湖，西靠黃河，南與萬(wàn)福河和東魚(yú)河搭界，北接梁濟(jì)運(yùn)河流域。洙趙新河流域包括洙趙新河干流以及鄆巨河、鄄鄆河、洙水河3條大的支流，干流起源于東明縣宋砦村，向東流經(jīng)8個(gè)縣(區(qū))于劉官屯村東入南陽(yáng)湖，全長(zhǎng)145.05 km，流域面積4 206 km2，其中菏澤市境內(nèi)4 119 km2，濟(jì)寧市境內(nèi)87 km2。該流域?qū)冱S泛沖積平原，地勢(shì)西高東低，西部地面高程57.40 m左右，東部濱湖地面高程為34.00 m左右，地面坡度在1/5000至1/12000之間。流域?qū)倥瘻貛啙駶?rùn)氣候區(qū)，多年平均氣溫為13.7℃，流域內(nèi)多年平均(1964-2009年)降水量為641 mm。

3 研究方法

3.1 適應(yīng)性評(píng)價(jià)指標(biāo)

本次研究采用NS系數(shù)和確定性系數(shù)R2來(lái)評(píng)價(jià)模型模擬的結(jié)果。研究者普遍認(rèn)為SWAT模型模擬滿意標(biāo)準(zhǔn)為R2>0.5，NS>0.5[15]；Moriasi等[16]認(rèn)為NS≥0.65(月尺度)就說(shuō)明模型的適用性是可以接受的。本文綜合選取的3種適用性評(píng)價(jià)指標(biāo)，參考陳軍鋒等[17]對(duì)梭磨河流域的研究，基于NS系數(shù)及R2的評(píng)價(jià)值可將模擬的結(jié)果劃分為4個(gè)等級(jí)，具體的方案見(jiàn)表1。

表1 徑流模擬評(píng)定等級(jí)標(biāo)準(zhǔn)值劃分表

3.2 不確定性分析算法

SWAT-CUP是獨(dú)立于SWAT模型的參數(shù)校準(zhǔn)分析軟件，其有運(yùn)算快、算法多、自動(dòng)繪圖等多個(gè)優(yōu)點(diǎn)，同時(shí)可應(yīng)用SWAT-CUP進(jìn)行參數(shù)的敏感性分析和不確定性分析。SWAT-CUP內(nèi)部鑲嵌了與SWAT鏈接的SUFI-2、PSO、GLUE、Parasol、MCMC算法程序。本次校準(zhǔn)采用Sufi-2算法對(duì)參數(shù)進(jìn)行自動(dòng)校準(zhǔn)，自動(dòng)校準(zhǔn)分為4步：設(shè)置初始參數(shù)、建立目標(biāo)函數(shù)、確定算法終止條件、計(jì)算評(píng)價(jià)指標(biāo)。

Sufi-2算法是一種反演建模法，參數(shù)的不確定性表示為范圍(均勻分布)，考慮了所有不確定性來(lái)源，例如驅(qū)動(dòng)變量(例如降水量)的不確定性、概念模型、參數(shù)和測(cè)量數(shù)據(jù)。參數(shù)中不確定性的傳播導(dǎo)致模型輸出變量的不確定性，表示為95%概率分布。這些是使用拉丁超立方采樣通過(guò)參數(shù)不確定性的傳播產(chǎn)生的輸出變量的累積分布的2.5%和97.5%水平計(jì)算的，被稱為95%預(yù)測(cè)不確定度，即95PPU，這些95PPU是隨機(jī)校準(zhǔn)方法中的模型輸出。重要的是沒(méi)有單個(gè)信號(hào)代表模型輸出，而是由某個(gè)參數(shù)范圍生成的95PPU表示的良好解決方案[18]。

4 構(gòu)建SWAT模型

4.1 數(shù)據(jù)預(yù)處理

洙趙新河流域SWAT模型模擬需要5類資料，其中數(shù)字高程數(shù)據(jù)(DEM)、土地利用和土壤類型數(shù)據(jù)用于模型的構(gòu)建，氣象數(shù)據(jù)資料(包括降水、溫度、風(fēng)速、太陽(yáng)輻射、相對(duì)濕度)用于驅(qū)動(dòng)模型，實(shí)測(cè)徑流資料用于模型參數(shù)率定，各類型數(shù)據(jù)的來(lái)源及用途見(jiàn)表2。

研究區(qū)的數(shù)字高程數(shù)據(jù)(DEM)是空間分辨率為3 s×3 s(90 m×90 m)的SRTM_DEM，該數(shù)據(jù)獲取自CGIAR-CSISRTM90數(shù)據(jù)庫(kù)，在ArcGIS中經(jīng)格式轉(zhuǎn)換、掩膜提取出SWAT模型所需的南四湖流域的DEM數(shù)據(jù)，數(shù)據(jù)如圖1所示。

表2 SWAT模型所需數(shù)據(jù)來(lái)源及用途

圖1 洙趙新河流域數(shù)字高程數(shù)據(jù)圖

土壤數(shù)據(jù)來(lái)自于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心的1∶100萬(wàn)土壤數(shù)據(jù)，空間分辨率為1 000 m×1 000 m。該數(shù)據(jù)格式同土地利用類型。土壤分類系統(tǒng)采用FAO-90系統(tǒng)。根據(jù)流域內(nèi)現(xiàn)有的100中土壤類進(jìn)行合并重新分類成31中土壤類型、11種土組(除水體、沙丘流沙外)，如圖2所示。

圖2 洙趙新河流域土壤類型分布圖

土地利用數(shù)據(jù)選取寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心的中國(guó)地區(qū)土地覆蓋綜合數(shù)據(jù)集[19]，數(shù)據(jù)集坐標(biāo)系為D_Krasovsky_1940，空間分辨率為1 000 m×1 000 m，數(shù)據(jù)格式為GRID，經(jīng)過(guò)投影變換成統(tǒng)一的投影坐標(biāo)系(WGS_1984_UTM_Zone_50N)，掩膜提取研究區(qū)的土地利用數(shù)據(jù)，并把原有的24種類型重分類成10種，如圖3所示。

流域內(nèi)氣象站點(diǎn)相對(duì)較少，且氣象數(shù)據(jù)不易獲得，本次研究采用寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心發(fā)布的SWAT模型CMADSV1.1數(shù)據(jù)集[11]。數(shù)據(jù)時(shí)間序列為2008-2014年，空間分辨率為1/3°，采用數(shù)據(jù)集站點(diǎn)共9個(gè)，具體站點(diǎn)分布見(jiàn)圖4。

圖3 洙趙新河流域土地利用類型分布圖

圖4 洙趙新河流域CMADS數(shù)據(jù)集站點(diǎn)分布圖

4.2 模型構(gòu)建

利用Mask對(duì)洙趙新河流域的DEM、土壤類型、土地利用數(shù)據(jù)進(jìn)行提取，并進(jìn)行流域坡度分析。根據(jù)對(duì)子流域劃分的研究分析及模型的推薦值，取洙趙新河流域的最小匯水面積的閾值為1 000 hm2，并基于DEM數(shù)據(jù)運(yùn)用“burn-in”算法生成流域模型中的的流域河網(wǎng)，運(yùn)用流域勾勒的八流點(diǎn)法則計(jì)算水流貢獻(xiàn)單元數(shù)，產(chǎn)生匯流柵格計(jì)算河道出水口，共劃分子流域個(gè)數(shù)57個(gè)，并進(jìn)行流域坡長(zhǎng)、坡度等特征參數(shù)計(jì)算，流域水系、子流域劃分分別見(jiàn)圖5、6。

流域內(nèi)共有10類土地利用類型、15類土壤類型。流域內(nèi)土壤類型較為單一，92%以上為黃河沖積土(主要為石灰性沖積土、鹽化沖積土)；土地利用類型80%以上為耕地。流域地處平坦湖西地區(qū)，流域坡度在0～58.08°之間，平均坡度為0.711°，流域坡度97%以上在2°及以下，流域坡度定義采用Multiple Slope分二級(jí)定義(0～6°、6～58.08°)；流域的HRU劃分采用Multiple HRUs方法進(jìn)行，并通過(guò)子流域內(nèi)土壤、土地利用所占面積分析，設(shè)置土地利用、土壤、坡度閾值為1%、9%、1%(相對(duì)面積)，共劃分384個(gè)HRUs。流域驗(yàn)證測(cè)站梁山閘站所在子流域編號(hào)為53。

圖5 洙趙新河實(shí)際水系與模擬水系圖

圖6 洙趙新河流域子流域劃分圖

4.3 徑流量還原計(jì)算

洙趙新河流域徑流量的校準(zhǔn)采用梁山閘水文站的實(shí)測(cè)月徑流量數(shù)據(jù)，由于洙趙新河是1973年完成全流域開(kāi)挖的坡水河，梁山閘站的流量過(guò)程受閘門影響嚴(yán)重，本次流域模擬在不考慮閘門影響的情況，對(duì)梁山閘站天然流量過(guò)程進(jìn)行基于降水?dāng)?shù)據(jù)的還原計(jì)算。

(1)

Ri=M1Ri-1+M2Ri

(2)

式中：Pi為研究時(shí)間范圍內(nèi)第i時(shí)段內(nèi)月降水量，mm；Ri、Ri-1分別為研究時(shí)間范圍內(nèi)第i時(shí)段及i-1時(shí)段對(duì)應(yīng)的月徑流量,m3/s；α、M1、M2分別為洙趙新河流域全區(qū)的徑流系數(shù)及由梁山閘、魏樓閘及鄄鄆河入河口測(cè)站的降水、徑流的歷史資料修正的相鄰時(shí)段間的排水模數(shù)。公式(2)中的Ri為最終用于模型率定的還原徑流量。梁山閘站降水量與修正徑流量過(guò)程線對(duì)比圖、修正徑流量過(guò)程及梁山閘站實(shí)測(cè)徑流與修正徑流過(guò)程線對(duì)比圖見(jiàn)圖7、8。

對(duì)于湖西平原區(qū)控制面積較大的排水骨干河流，設(shè)計(jì)排水模數(shù)可采用經(jīng)驗(yàn)公式法或單位線法，但以經(jīng)驗(yàn)公式法較常用，本文采用經(jīng)驗(yàn)公式形式如下：

M=KRmFn

(3)

式中：M為設(shè)計(jì)排澇模數(shù)，m3/(s·km2)；K為綜合系數(shù)；R為設(shè)計(jì)凈雨深，mm；F為流域面積，km2；m為峰量指數(shù)；n為排澇模數(shù)隨面積增大而減小的指數(shù)。流域內(nèi)梁山閘控制站的不同頻率的排水模數(shù)如表3所示。

表3 洙趙新河梁山閘站斷面排水模數(shù)成果表

5 結(jié)果分析

研究區(qū)模擬時(shí)間選定2008-2014年，將各類數(shù)據(jù)向前延伸1年得到2007年的5類數(shù)據(jù)，作為預(yù)熱期；2008-2010年為校準(zhǔn)期，進(jìn)行流域徑流、泥沙相關(guān)參數(shù)的調(diào)整與校準(zhǔn)；2011-2014年作為驗(yàn)證期，以檢驗(yàn)建立模型的模擬效果。

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

本次模擬對(duì)常用的25個(gè)參數(shù)進(jìn)行率定，主要為地下水參數(shù)、土壤參數(shù)等。參數(shù)取值范圍參考SWAT模型數(shù)據(jù)庫(kù)，通過(guò)在SWAT-CUP中的Sufi-2算法進(jìn)行1000次迭代運(yùn)算，比較參數(shù)的敏感性，最終確定13個(gè)對(duì)徑流量的敏感性最大參數(shù)，參數(shù)模擬的最優(yōu)結(jié)果及敏感性排序如表4所示。P-Value 是t檢驗(yàn)值表對(duì)應(yīng)的P概率值。P-Value值越接近于0，參數(shù)敏感性越高。

分析結(jié)果表明植被蒸騰補(bǔ)償系數(shù)、河岸調(diào)蓄的基流α因子、淺層地下水徑流系數(shù)的P-Value小于0.1，說(shuō)明對(duì)徑流量的敏感性較高。

5.2 模擬結(jié)果評(píng)價(jià)

洙趙新河流域的徑流量校準(zhǔn)期、驗(yàn)證期的模擬年系列值與實(shí)測(cè)年系列值對(duì)比圖及相關(guān)關(guān)系圖見(jiàn)圖9、10。

由洙趙新河流域校準(zhǔn)期、驗(yàn)證期的實(shí)測(cè)與模擬月徑流量對(duì)比圖可知，確定性系數(shù)在0.80以上，模擬效果基本滿意；洙趙新河流域的模擬與實(shí)測(cè)月徑流量過(guò)程的各評(píng)價(jià)指標(biāo)值均滿足模擬精度要求，模擬等級(jí)屬于乙級(jí)及以上，具體見(jiàn)表5所示。由圖表可知，模擬流域的實(shí)測(cè)月徑流系列值與模擬系列值趨于一致，模擬等級(jí)較高，豐水期的峰量控制較好，平枯水期及偏枯年份的模擬一般，模型整體模擬效果較好，滿足研究需要。

圖7梁山閘站降水量與修正徑流量過(guò)程圖8梁山閘站實(shí)測(cè)徑流與修正徑流過(guò)程線對(duì)比圖

表4 參數(shù)率定結(jié)果及敏感性排序

圖9 洙趙新河流域校準(zhǔn)期實(shí)測(cè)與模擬月徑流量對(duì)比圖及相關(guān)關(guān)系圖

圖10 洙趙新河流域驗(yàn)證期實(shí)測(cè)與模擬月徑流量對(duì)比圖及相關(guān)關(guān)系圖

時(shí)期月均值/(m3·s-1)實(shí)測(cè)值模擬值R2NS等級(jí)校準(zhǔn)期25.34325.1760.87630.8886乙級(jí)驗(yàn)證期15.15415.2570.92710.9267甲級(jí)

6 結(jié) 論

本文用SWAT模型模擬了洙趙新河流域徑流情況，并結(jié)合Sufi-2方法率定模型參數(shù)并進(jìn)行敏感性分析，得到以下主要結(jié)論：

(1)在如洙趙新河流域這樣氣象站點(diǎn)較少的地區(qū)，引入CMADS氣象數(shù)據(jù)集作為流域水文模型的驅(qū)動(dòng)數(shù)據(jù)集，有效解決氣象數(shù)據(jù)難獲取、流域氣象站點(diǎn)少的問(wèn)題，在模型精確度、空間驅(qū)動(dòng)分辨率上得到提升。

(2)根據(jù)敏感性分析結(jié)果發(fā)現(xiàn)水文相關(guān)參數(shù)與徑流有不同程度的相關(guān)性，其中植被蒸騰補(bǔ)償系數(shù)、河岸調(diào)蓄的基流α因子、淺層地下水徑流系數(shù)、基流α因子、主河道有效滲透系數(shù)、土壤蒸發(fā)補(bǔ)償系數(shù)、土壤表面到層底的深度、地下水滯后系數(shù)、土壤有效含水量、地下水再蒸發(fā)系數(shù)等參數(shù)對(duì)本研究區(qū)徑流模擬影響相對(duì)明顯。

(3)SWAT模型對(duì)豐水年份模擬效率高，對(duì)于某些枯水年模擬不夠理想，并且對(duì)徑流量特征呈“雙峰型”的時(shí)段模擬與實(shí)測(cè)徑流量相差較大?？傮w來(lái)看，SWAT模型模擬較為理想，模擬月徑流時(shí)，率定期R2和NS值均大于0.85，驗(yàn)證期超過(guò)了0.9，表明SWAT模型在洙趙新河流域適用性較好。