夏曉玲,曾莉萍,劉 濤,3,王加敏,方 荻,2,張明祥

(1.貴州新氣象科技有限責(zé)任公司,貴州 貴陽 550002; 2.貴州省山地環(huán)境氣候研究所,貴州 貴陽 550002; 3.貴州省人工影響天氣辦公室,貴州 貴陽 550002)

0 引 言

貴州省地處云貴高原東斜坡地帶,氣候濕潤(rùn),降雨豐沛,河流眾多,水資源充足。但貴州省下墊面情況復(fù)雜,全省73%的面積為碳酸巖鹽[1],92.5%的面積為山地和丘陵,山間平壩區(qū)面積甚少[2],導(dǎo)致貴州水資源時(shí)空分布非常不均勻、水資源供需矛盾突出,難以有效利用水資源。隨著社會(huì)經(jīng)濟(jì)的快速發(fā)展,水資源的可持續(xù)利用面臨著更加巨大的挑戰(zhàn)[3-4]。水力發(fā)電作為水資源開發(fā)利用的重要手段之一,需要精細(xì)化的水文(徑流量)預(yù)報(bào)為其提高開發(fā)利用率提供支撐。隨著“3S”等科學(xué)技術(shù)的發(fā)展,Soil and Water Assessment Model(SWAT模型)成為運(yùn)用最為廣泛的分布式水文模型之一。該模型可以在下墊面復(fù)雜、降水時(shí)空分布不均情況下對(duì)流域水文過程進(jìn)行精細(xì)化模擬[5-7]。

2015年,Uniyal等[8]使用ArcSWAT模型研究氣候變化對(duì)印度拉尼河水平衡等的影響,發(fā)現(xiàn)21世紀(jì)末氣候條件變化會(huì)對(duì)該區(qū)域的徑流產(chǎn)生重大影響。2017年,Reshmidevi等[9]使用與ArcGIS集成的SWAT進(jìn)行水文模擬,預(yù)測(cè)了水資源壓力可能來源于河流流量、地下水補(bǔ)給量的減少以及灌溉需求的增加。中國(guó)針對(duì)分布式水文模型的研究發(fā)展迅速[10]。1997年,黃平等[11]總結(jié)了具有物理機(jī)制的分布式水文模型,構(gòu)建了流域三維動(dòng)態(tài)水文數(shù)值模型。2007年,熊立華等[12]利用拓?fù)潢P(guān)系表達(dá)流域水流的空間聚合與分散,建立了水文模型。2013年,于嵐嵐等[13]在大凌河流域進(jìn)行小流域降雨模擬,針對(duì)洪水傳播時(shí)間為固定值這個(gè)缺陷,采用變動(dòng)態(tài)存儲(chǔ)系數(shù)法演算洪水,提高了洪水模擬的合格率。2015年,朱悅璐等[14]將CanESM2氣候模式下的兩種情景與半分布式水文模型VIC相對(duì)接,分析未來2020s、2030s、2040s、2050s四個(gè)時(shí)期渭河徑流變化。2020～2021年間,劉君龍、劉飛、徐志等均運(yùn)用SWAT模型對(duì)不同流域的水文過程進(jìn)行了模擬[15-17]。陳長(zhǎng)征等[18]對(duì)SWAT模型的參數(shù)化優(yōu)化方法進(jìn)行了研究分析。國(guó)內(nèi)外研究大多從水文過程著手,通過優(yōu)化水文模型以提高徑流量或洪水預(yù)報(bào)。本文基于ArcGIS軟件平臺(tái),針對(duì)貴州省平寨流域,首次在喀斯特地貌特征明顯的區(qū)域應(yīng)用SWAT分布式水文模型開展徑流量預(yù)測(cè)研究,以期為類似流域研究提供參考。

1 研究區(qū)域概況

平寨水庫壩址在三岔河中游木底河平寨附近,處于貴州省六枝特區(qū)與織金縣交界處(圖1)。平寨水庫壩高162.7 m,淹沒區(qū)面積為15.10 km2,正常蓄水位1 331 m,死水位1 305 m,平均水深50 m,校核洪水位1 333.29 m,總庫容量為10.89億m3,調(diào)節(jié)庫容4.48億m3,電站總裝機(jī)容量146.6 MW。流域總面積3 492 km2,占三岔河流域面積的48.2%。

圖1 平寨流域

平寨流域范圍內(nèi)現(xiàn)有數(shù)據(jù)較為完整的水文站為陽長(zhǎng)水文站,1993年1月1日該水文站開始記錄數(shù)據(jù)。陽長(zhǎng)水文站由畢節(jié)地區(qū)水文水資源局設(shè)立,控制集水面積為2 696 km2,在納雍縣、陽長(zhǎng)鎮(zhèn)、新陽長(zhǎng)大橋下游約50 m的位置。氣象站點(diǎn)選用水城、汪家寨、松林坡、新房、牛場(chǎng)共5個(gè)站點(diǎn)。

2 數(shù)據(jù)與方法

2.1 建模數(shù)據(jù)來源

建立SWAT模型所需要的基礎(chǔ)數(shù)據(jù)庫主要分為兩類:空間數(shù)據(jù)和屬性數(shù)據(jù)。

(1) 空間數(shù)據(jù)包括高程數(shù)據(jù)(DEM)、土地利用類型數(shù)據(jù)、河流水系分布數(shù)據(jù)、土壤類型分布數(shù)據(jù)等。其中DEM高程圖和土地利用數(shù)據(jù)的精度為30 m×30 m,格式為Raster;土地利用數(shù)據(jù)來源是2015年全國(guó)范圍的遙感影像(基于Landsat),分辨率為1 km,分為14類;河流水系圖精度為1∶500 000,格式為Shapefile。以上所有空間數(shù)據(jù)需具有統(tǒng)一的投影坐標(biāo),本文投影采用WGS_1984_UTM_Zone_48。

(2) 屬性數(shù)據(jù)主要是氣象和水文的觀測(cè)資料,氣象資料:氣溫(最低、最高)、降水量、日照時(shí)數(shù)、平均相對(duì)濕度、平均風(fēng)速等;時(shí)間范圍為2009年1月1日至2018年12月31日,時(shí)間分辨率為逐日數(shù)據(jù)或逐月數(shù)據(jù);水文資料為陽長(zhǎng)水文站逐日/逐月的流量觀測(cè),逐日水文資料時(shí)間范圍為2017年1月1日至2019年12月31日;逐月水文資料時(shí)間范圍為2009年1月至2018年12月。

2.2 預(yù)報(bào)、分析數(shù)據(jù)來源

氣象數(shù)據(jù)包括平寨流域關(guān)鍵氣象站點(diǎn)(水城、汪家寨、松林坡、新房、牛場(chǎng))2021年6～9月每月9,19,29日的未來8 d降水預(yù)報(bào)數(shù)據(jù)。實(shí)況數(shù)據(jù)包括2021年6～9月平寨水庫日入庫流量數(shù)據(jù)。

2.3 研究方法

2.3.1 模型評(píng)價(jià)指標(biāo)

本文的評(píng)價(jià)指標(biāo)有確定系數(shù)R2、納什系數(shù)NSE、偏差百分率PBIAS。公式如下:

(1)

(2)

(3)

式中:Qobs,Qobsave,Qsim,Qsimave分別為觀測(cè)值、觀測(cè)值的平均值、SWAT模型的模擬值、模擬值的平均值,m3/s;n為數(shù)據(jù)長(zhǎng)度,此處為年份。確定系數(shù)和納什系數(shù)取值范圍在0～1,值越接近1,說明模型的模擬效果越好。確定系數(shù)大于等于0.5時(shí),認(rèn)為模擬效果可以被接受。納什系數(shù)大于等于0.75時(shí),表示模型的模擬效果很好;當(dāng)NSE為0.36～0.75時(shí),說明模型的模擬結(jié)果基本達(dá)要求;如果納什系數(shù)不大于0.36,說明模擬結(jié)果達(dá)不到要求。偏差百分率越接近0說明模擬效果越好,絕對(duì)值在10%以內(nèi)表示可以接受。

2.3.2 模型不確定性分析

運(yùn)用SWAT-CUP軟件對(duì)模型參數(shù)進(jìn)行率定。SWAT-CUP是獨(dú)立于SWAT與ArcGIS,專門開發(fā)的計(jì)算機(jī)程序。本文使用SWAT-CUP軟件中SUFI2算法對(duì)平寨流域SWAT模擬徑流量結(jié)果進(jìn)行參數(shù)自動(dòng)校準(zhǔn),獲取12個(gè)影響最大的參數(shù),并代回模型反算結(jié)果進(jìn)行驗(yàn)證。

2.3.3 實(shí)時(shí)入庫流量修正

目前常用的入庫流量計(jì)算方法有水量平衡法、區(qū)間推流法、入庫控制站代表法等,其中以水量平衡法應(yīng)用最廣。該算法基于水量平衡原理,用時(shí)段內(nèi)入庫水量減去出庫水量得到庫內(nèi)水量變化值,利用庫容變化和出庫流量反推計(jì)算入庫流量。但該算法對(duì)壩上水位數(shù)據(jù)的敏感度較高,特別是水庫面積較大,壩上水位的小幅波動(dòng)會(huì)導(dǎo)致庫容差大幅變動(dòng),從而影響入庫流量計(jì)算,使入庫流量呈現(xiàn)鋸齒狀震蕩,無法準(zhǔn)確反映實(shí)際情況,有些水庫在枯水期甚至出現(xiàn)計(jì)算入庫流量為“負(fù)”值的現(xiàn)象。

用m次多項(xiàng)式來擬合得到的實(shí)驗(yàn)數(shù)據(jù),設(shè)擬合多項(xiàng)式為

Y(t)=a0+a1t+a2t2+…+amtm

(4)

用最小乘法來確定方程(4)中的待定系數(shù),令

(5)

=φ(a0,a1,…,am)

(6)

為使φ(a0,a1,…,am)達(dá)到最小,將它分別對(duì)ak(k=0,1,…,m)求偏導(dǎo)數(shù),并令其為0,可得方程組:

(7)

該方程組稱為正規(guī)方程組。當(dāng)n=2(5個(gè)節(jié)點(diǎn)),m=3時(shí),得到具體的正規(guī)方程組,由此解出a0,a1,ai,ait代入式(4),并令t=0,1,-1,2,-2,得到五點(diǎn)三次平滑公式:

(8)

(9)

(10)

(11)

(12)

該算法要求節(jié)點(diǎn)個(gè)數(shù)為k≥5,當(dāng)節(jié)點(diǎn)個(gè)數(shù)多于5時(shí),為對(duì)稱起見,除在兩端分別用式(8)、(9)、(10)、(11)、(12)外,其余都用式(10)進(jìn)行平滑,這就相當(dāng)于在每個(gè)子區(qū)間上用不同的三次最小二乘多項(xiàng)式進(jìn)行平滑。對(duì)于平寨入庫流量運(yùn)用五點(diǎn)三次算法進(jìn)行平滑處理的結(jié)果,平滑后的流量數(shù)據(jù)走勢(shì)與原始數(shù)據(jù)基本相同,且洪水總量、洪峰流量與峰現(xiàn)時(shí)間洪水三要素與原始數(shù)據(jù)差異較小,能夠較好擬合原始流量數(shù)據(jù),具有較高的實(shí)際應(yīng)用價(jià)值。

2.3.4 水文氣象模式耦合

將降水預(yù)報(bào)產(chǎn)品作為輸入場(chǎng)驅(qū)動(dòng)水文模型,實(shí)現(xiàn)水文氣象模式耦合,進(jìn)而開展徑流量預(yù)報(bào),旨在延長(zhǎng)洪水預(yù)報(bào)的預(yù)見期。降水預(yù)報(bào)數(shù)據(jù)是貴州新氣象科技有限責(zé)任公司在結(jié)合貴州省氣象臺(tái)與ECWMF模式降水預(yù)報(bào)產(chǎn)品的基礎(chǔ)上制作的主客觀融合格點(diǎn)降水、氣溫、風(fēng)等產(chǎn)品,這些數(shù)據(jù)分辨率為5 km×5 km,屬于精細(xì)化產(chǎn)品。將這些數(shù)據(jù)插值至平寨流域范圍內(nèi)的5個(gè)氣象站點(diǎn)。雨量訂正方法為首先隨機(jī)選取一部分?jǐn)?shù)據(jù),運(yùn)用一元線性回歸方式建立訂正模型,然后運(yùn)用其余數(shù)據(jù),帶入訂正模型,檢驗(yàn)訂正后誤差百分率,再用消除極值方法,保證預(yù)報(bào)誤差在10%以內(nèi),作為訂正后的雨量預(yù)報(bào)數(shù)據(jù),和訂正前預(yù)報(bào)數(shù)據(jù)同時(shí)帶入水文模型對(duì)比分析。

3 結(jié)果分析

3.1 SWAT在平寨流域的模型建立與檢驗(yàn)

3.1.1 SWAT模型數(shù)據(jù)庫建立

獲取基礎(chǔ)數(shù)據(jù)后,根據(jù)不同地理信息數(shù)據(jù),平寨河流域被劃分出2 117個(gè)水文響應(yīng)單元、23個(gè)子流域。陽長(zhǎng)水文監(jiān)測(cè)站被劃分在第11個(gè)子流域,運(yùn)用陽長(zhǎng)水文站觀測(cè)的實(shí)測(cè)流量數(shù)據(jù)對(duì)模型進(jìn)行率定和驗(yàn)證,2017年日徑流值為預(yù)熱期,2018年為率定期,2019年為驗(yàn)證期;2009～2010年月徑流值為預(yù)熱期,2011～2014年為率定期,2015～2018年為驗(yàn)證期。圖2為平寨流域SWAT模型建模界面。

圖2 平寨流域SWAT模型建模界面展示

3.1.2 參數(shù)敏感性分析

由于陽長(zhǎng)站的水文數(shù)據(jù)時(shí)間長(zhǎng)度最長(zhǎng),數(shù)據(jù)質(zhì)量控制最好,因此參數(shù)率定選擇陽長(zhǎng)站,經(jīng)過SUFI2算法2 000次迭代計(jì)算常用的25個(gè)參數(shù)后,發(fā)現(xiàn)12個(gè)對(duì)陽長(zhǎng)站徑流量影響最大的參數(shù)(表1)。結(jié)果表明:徑流曲線數(shù)(CN2)對(duì)徑流影響最大,土壤有效含水量(SOL_AWC)、地下水滯后系數(shù)(GW_DELAY)、和基流α系數(shù)(ALPHA_BF)次之,淺層地下水徑流系數(shù)(GWQMN)影響最小。

表1 參數(shù)敏感性分析

3.1.3 率定及驗(yàn)證

借助SWAT-CUP軟件,運(yùn)用陽長(zhǎng)站實(shí)際觀測(cè)的徑流量對(duì)模擬日徑流量(2018年)和模擬月徑流量(2011～2014年)進(jìn)行率定,將最佳模擬結(jié)果的參數(shù)值代回SWAT模型,重新模擬陽長(zhǎng)站2019年日徑流量及2015～2018年月徑流量,并與實(shí)際觀測(cè)的徑流量做對(duì)比分析(表2)。可以看出,模擬日徑流量值校準(zhǔn)效果優(yōu)于模擬月徑流量值,日值與月值率定期R2和NSE均為0.8左右,驗(yàn)證期日值R2和NSE均大于0.7,而月值NSE小于0.7。率定期日值和月值的PBIAS絕對(duì)值均小于10%。整體而言,驗(yàn)證期評(píng)價(jià)指標(biāo)基本在誤差范圍內(nèi),因此SWAT模型對(duì)平寨流域的徑流過程模擬有一定的效果。

表2 SWAT模型在平寨流域陽長(zhǎng)站率定期、驗(yàn)證期的模擬結(jié)果

平寨流域陽長(zhǎng)站率定期、驗(yàn)證期的水文觀測(cè)資料與模型模擬值對(duì)比結(jié)果(圖3～4)表明:從變化趨勢(shì)來看,模擬值與實(shí)測(cè)值基本一致;模擬的洪峰流量、洪峰值出現(xiàn)時(shí)間也和實(shí)際觀測(cè)的較為吻合。日值模擬中,在率定期2018年6月22日、7月11日、9月26日分別出現(xiàn)了3次洪峰流量過程,第一次和第三次過程模擬值低于實(shí)測(cè)值,第二次過程中模擬值較實(shí)測(cè)值高;2019年驗(yàn)證期,洪峰過程分別出現(xiàn)在6月18日、7月23日、9月9日,第一次過程模擬值與實(shí)測(cè)值較為接近,其他過程模擬值均小于實(shí)測(cè)值。月值模擬中,率定期2011年6月、2012年7月、2014年7月出現(xiàn)了3次洪峰流量,徑流量模擬值在第一次洪峰過程中顯著高于實(shí)測(cè)值,而在第二次和第三次過程中略小于實(shí)測(cè)值。整體而言,SWAT模型模擬出的平寨流域的徑流量可基本滿足日常應(yīng)用的要求。

圖3 平寨流域陽長(zhǎng)站逐日徑流量實(shí)測(cè)值與模擬值對(duì)比

圖4 平寨流域陽長(zhǎng)站逐月徑流量實(shí)測(cè)值與模擬值對(duì)比

3.2 耦合氣象預(yù)報(bào)數(shù)據(jù)預(yù)測(cè)入庫流量

在平寨流域SWAT模型中加入水城、汪家寨、松林坡、新房、牛場(chǎng)氣象站點(diǎn)2021年6～9月每月9,19,29 d的未來8 d降水預(yù)報(bào)數(shù)據(jù)。由于陽長(zhǎng)站和平寨水庫距離較為接近,且平寨水庫水位流量站建站時(shí)間較短且數(shù)據(jù)質(zhì)量不高,因此無法運(yùn)用平寨水庫流量站開展參數(shù)率定建模工作,本文嘗試運(yùn)用陽長(zhǎng)站點(diǎn)的參數(shù)率定結(jié)果,代入模型對(duì)平寨水庫的入庫流量進(jìn)行預(yù)測(cè)。

3.2.1 整體預(yù)報(bào)效果分析

針對(duì)2021年6～9月平寨水庫的SWAT模型入庫流量預(yù)報(bào),分別計(jì)算逐月的確定系數(shù)R2、NSE系數(shù)、百分比偏差PBIAS(表3),可以看出,2021年8月和9月的NSE系數(shù)為負(fù),其余時(shí)段的NSE系數(shù)為正,最高為2021年6月,NSE系數(shù)為0.45,說明耦合了氣象預(yù)報(bào)數(shù)據(jù)的SWAT模型在平寨流域的模擬結(jié)果在6～7月基本達(dá)到要求,但是8～9月模擬效果不佳。將降水預(yù)報(bào)偏差較大的日期進(jìn)行一定的訂正,運(yùn)用訂正后的降水模擬平寨水庫的入庫流量,可以看出6～9月的NSE系數(shù)和R2相比降水訂正前有明顯的提高,偏差有一定縮小,但整體表現(xiàn)沒有建模時(shí)期的效果好,最高納什系數(shù)僅為0.61,降水修正后8～9月的NSE系數(shù)由原先的負(fù)值變?yōu)檎?說明由于8～9月的降水預(yù)報(bào)偏差較大,造成了入庫流量預(yù)報(bào)效果不佳。

表3 2021年主汛期SWAT模型在平寨水庫預(yù)報(bào)效果

3.2.2 逐日入庫流量情況預(yù)報(bào)效果分析

從逐日入庫流量預(yù)測(cè)和實(shí)況對(duì)比(圖5)可以看出,在2021年7月2日前后的洪峰過程,SWAT模型的模擬效果和實(shí)況流量基本一致,預(yù)測(cè)流量和出現(xiàn)的時(shí)間和實(shí)況吻合程度較高,預(yù)測(cè)的洪峰為224.7 m3/s,6月30日的實(shí)況流量為250.78 m3/s,隨后的幾天預(yù)測(cè)效果也比較理想。但是2021年6月9日前后的洪峰過程預(yù)測(cè),在時(shí)間上有一定的滯后,流量數(shù)值的預(yù)報(bào)和實(shí)況較為接近,7月27日前后以及8月30日前后小洪峰過程預(yù)測(cè)的入庫流量明顯小于實(shí)況。將降水訂正和實(shí)況對(duì)比,6～9月的預(yù)報(bào)流量走勢(shì)和洪峰過程均較訂正前更接近實(shí)測(cè)值,且8～9月的洪峰過程預(yù)報(bào)效果有明顯的提升,因此運(yùn)用SWAT模型耦合降水預(yù)報(bào)對(duì)平寨水庫的入庫流量進(jìn)行預(yù)測(cè)時(shí),在模型參數(shù)固定的條件下降水預(yù)報(bào)的準(zhǔn)確性起到關(guān)鍵作用。

圖5 平寨水庫降水訂正前后及入庫流量預(yù)報(bào)

4 結(jié) 論

本文以貴州省平寨流域作為研究區(qū),基于DEM高程資料、土地利用數(shù)據(jù)、土壤數(shù)據(jù)、氣象數(shù)據(jù)、陽長(zhǎng)水文站實(shí)測(cè)徑流量數(shù)據(jù)等,應(yīng)用SWAT模型,對(duì)其水文循環(huán)過程進(jìn)行了模擬預(yù)測(cè),具體結(jié)論如下。

(1) 利用實(shí)測(cè)徑流量資料對(duì)2017～2019年逐日徑流量和2009～2018年逐月徑流量進(jìn)行了率定及驗(yàn)證,其中日值模擬中設(shè)定2017年為預(yù)熱期,2018年為率定期,2019年為驗(yàn)證期,月值模擬中設(shè)定2009～2010年為預(yù)熱期,2011～2014年為率定期,2015～2018年為驗(yàn)證期。參數(shù)敏感性分析結(jié)果顯示CN2參數(shù)最敏感,其次為SOL_AWC;GWQMN對(duì)徑流量影響最小。

(2) 逐日數(shù)據(jù)率定期R2和NSE均為0.88,PBIAS小于10%,驗(yàn)證期R2和NSE均大于0.7,逐月數(shù)據(jù)率定期R2大于0.8,NSE為0.79,PBIAS為3.7%,驗(yàn)證期R2大于0.75,NSE為0.69,模擬效果達(dá)到要求,在平寨流域有很好的適用性。

(3) SWAT模型模擬的徑流量變化趨勢(shì)與實(shí)測(cè)徑流量變化趨勢(shì)基本一致。日值模擬的洪峰過程與實(shí)際過程相符,均能準(zhǔn)確模擬。月值模擬中,2016年峰值模擬效果最優(yōu)。

(4) 在平寨流域中SWAT模型耦合氣象預(yù)報(bào)數(shù)據(jù),模擬得到該時(shí)段徑流量變化,發(fā)現(xiàn)幾次洪峰過程預(yù)報(bào)有明顯的偏大的情況,洪峰出現(xiàn)的時(shí)間也有一定的錯(cuò)位。如果將降水訂正至和實(shí)況接近,這幾次洪峰過程預(yù)報(bào)降水會(huì)更接近實(shí)況,因此運(yùn)用SWAT數(shù)據(jù)耦合降水預(yù)報(bào)對(duì)平寨水庫的入庫流量進(jìn)行預(yù)測(cè)時(shí),降水預(yù)報(bào)的準(zhǔn)確性較為重要。