基于SWAT與新安江模型的閩江建陽流域徑流模擬研究

肖 豪，周春輝，尚艷麗，嵇海祥，鄭何聲園，瞿思敏，石 朋

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210024；2.國網(wǎng)新源水電有限公司富春江水力發(fā)電廠，浙江 杭州 310000；3. 山東省泰安市水文中心，山東 泰安 271000；4. 水利部南京水利水文自動化研究所，江蘇 南京 210008)

0 引 言

閩江流域地處福建省北部，河長577 km，面積約6萬km2，流域的水資源變化深刻的影響著福建省的經(jīng)濟和環(huán)境。該流域?qū)賮啛釒Ъ撅L(fēng)濕潤氣候，雨量充沛，洪水頻發(fā)。近年來，流域內(nèi)的水循環(huán)和水資源在全球氣候變化和人類活動的影響下發(fā)生了重大的變化。實時精確的徑流模擬和流域尺度上的環(huán)境變化深刻影響著徑流時空分布的影響研究。流域水文模型是對流域的水文過程進行模擬分析，在洪水預(yù)報、防洪減災(zāi)、水資源開發(fā)利用、徑流對變化環(huán)境響應(yīng)研究等方面得到了廣泛的應(yīng)用[1]。目前，在我國應(yīng)用最廣泛的流域水文模型是新安江模型。根據(jù)研究實踐，新安江模型適用于我國南方濕潤地區(qū)[2]，被用在我國各大流域評估氣候變化和人類活動對徑流的影響[3-5]。美國農(nóng)業(yè)部研發(fā)的SWAT模型是一個大尺度分布式水文模型，自開發(fā)以來在國內(nèi)外得到了廣泛的認可和應(yīng)用[6-12]，被廣泛用于水流、泥沙、水質(zhì)模擬及變化環(huán)境的徑流響應(yīng)研究[13-16]。比如，朱炬明等[6]在雙橋流域探討了新安江模型、SWAT模型和BTOPMC模型在不同流量級下的日徑流模擬效果。劉柏君等[7]劃分青海省主要流域并構(gòu)建SWAT模型進行徑流模擬，分析了模型的適用性。林若蘭等[8]基于SWAT模型評估北江流域生態(tài)徑流調(diào)節(jié)效果并分析徑流調(diào)節(jié)量的時空變化。但是，在閩江流域，同時利用新安江模型和SWAT模型進行徑流模擬并基于模擬結(jié)果評估徑流時空變化特征的研究并不多見。

為了提高閩江流域徑流模擬精度并分析變化環(huán)境對徑流時空分布影響，本文以福建閩江上的建陽流域作為研究區(qū)域，利用空間數(shù)據(jù)和屬性數(shù)據(jù)構(gòu)建SWAT模型，研究SWAT模型在研究區(qū)的適用性和模擬精度，并與三水源新安江模型的日徑流模擬結(jié)果進行對比分析；再利用SWAT模型的輸出結(jié)果，分析研究區(qū)徑流的時空分布變化特征，為研究流域的防洪減災(zāi)和水資源規(guī)劃提供技術(shù)支持。

1 流域概況及資料

1.1 建陽流域概況

本文選擇福建閩江上游的建陽流域作為研究區(qū)(見圖1)。研究區(qū)流域面積4 848 km2，范圍為東經(jīng)117°31′～118°19′，北緯27°10′～28°5′。流域位于福建省北部，屬亞熱帶季風(fēng)濕潤氣候，降雨充足，多年平均降水量1 702 mm，降雨年際變化較大，空間分布不均勻，總體趨勢從北向南遞減，流域內(nèi)多年平均徑流系數(shù)為0.56。該流域植被較好，地下水埋深較淺，是典型的山丘區(qū)地形，洪水大多數(shù)是由暴雨形成[17-18]。流域內(nèi)有武夷山、興田等17個雨量站和建陽水文站(見圖1)。

1.2 流域資料情況

流域內(nèi)的水文資料包括建陽水文站的日徑流數(shù)據(jù)以及17個雨量站1993年～2000年的日降水?dāng)?shù)據(jù)；氣象資料有研究區(qū)內(nèi)及周邊的氣象站數(shù)據(jù)。新安江模型計算時選用了研究區(qū)內(nèi)的雨量站數(shù)據(jù)和水文站數(shù)據(jù)來率定和檢驗參數(shù)。

構(gòu)建SWAT模型需要輸入的數(shù)據(jù)除了屬性數(shù)據(jù)以外，還需要空間數(shù)據(jù)。此外，因為流域內(nèi)氣象站數(shù)據(jù)序列時間長短不一致。本文綜合考慮選取流域內(nèi)及周邊邵武、武夷山、浦城3個氣象站點的數(shù)據(jù)用于構(gòu)建SWAT模型的天氣發(fā)生器。研究區(qū)內(nèi)及周邊站點的分布如圖1所示。

本文所使用的數(shù)據(jù)格式、來源及其用途如表1所示。

表1 原始數(shù)據(jù)說明

2 流域水文模型及評價指標(biāo)

2.1 新安江模型

新安江模型是河海大學(xué)趙人俊教授等人于1973年提出并逐漸完善起來的分布式概念性降雨徑流水文模型[19]。該模型結(jié)構(gòu)簡單，可以分為蒸散發(fā)計算、產(chǎn)流計算、分水源計算和匯流計算4個層次。蒸散發(fā)計算采用3層蒸發(fā)模式，產(chǎn)流計算采用蓄滿產(chǎn)流，分水源計算采用自由水蓄水庫，利用線性水庫計算坡地匯流，分段馬斯京根法計算河網(wǎng)匯流。模型的參數(shù)具有明確的物理意義。計算采用的方法和參數(shù)見表2。具體計算公式可參考文獻[20]。

表2 新安江模型計算方法與參數(shù)

2.2 SWAT模型

SWAT(Soil and Water Assessment Tool)模型是美國農(nóng)業(yè)部研發(fā)的一個具有很強物理機制的分布式水文模型[21]。該模型基于地理信息系統(tǒng)(GIS)并利用遙感(RS)和GIS提供的研究流域的空間信息，疊加分析土壤類型和土地利用方式，模擬復(fù)雜流域的水文物理過程從而分析下墊面在流域水循環(huán)的作用。模型包含氣象、水文、泥沙輸移、植被生長、營養(yǎng)物質(zhì)、污染物和農(nóng)作物管理等模塊[22-23]，本研究主要應(yīng)用的是模型的水文模擬模塊。

模型計算包括地表徑流、下滲計算、蒸散發(fā)計算、河道匯流計算。地表徑流量采用徑流曲線數(shù)法( SCS runoff curve number)；下滲采用土壤蓄水演算技術(shù)；蒸散發(fā)分為土壤水蒸發(fā)和植物蒸騰；河道匯流計算使用變動存儲系數(shù)模型或Muskingum方法，計算模塊和部分參數(shù)見表3。

表3 SWAT模型計算模塊和參數(shù)

2.3 模擬評價指標(biāo)

本文選用納什效率系數(shù)NS、相對誤差RE和決定性系數(shù)R2來評價模型模擬精度。

(1)納什效率系數(shù)NS，表征模型模擬徑流過程與實測徑流過程之間擬合程度。NS取值范圍為負無窮到1；越接近于1，則認為模擬的流量過程越接近觀測值。一般NS≥0.6時，模擬結(jié)果可以接受。計算公式為

(1)

(2)相對誤差RE，反應(yīng)了徑流總量模擬值和實測值之間的吻合程度。一般認為RE在±20%以內(nèi)，模擬結(jié)果具有可信度。計算公式為

(2)

(3)決定性系數(shù)R2，用于體現(xiàn)模型模擬值與實測值的相關(guān)性，R2越接近1，表明模型模擬值與實測值的擬合程度越高，一般R2>0.6可認為擬合效果較好。計算公式為

(3)

3 參數(shù)率定與模擬結(jié)果分析

兩個模型均采用1993年～1997年的日資料進行模型參數(shù)率定，1998年～2000年的日資料進行模型檢驗。

3.1 新安江模型參數(shù)率定

由建陽流域的地理位置、降雨特性和植被條件等因素可以看出，研究區(qū)的降雨徑流關(guān)系滿足蓄滿產(chǎn)流條件，故使用新安江模型進行徑流模擬。采用系統(tǒng)微分響應(yīng)方法對模型參數(shù)進行率定[24-28]，結(jié)果見表4。

表4 建陽流域新安江模型參數(shù)率定結(jié)果

3.2 SWAT模型參數(shù)率定

(1)構(gòu)建SWAT模型。對空間數(shù)據(jù)進行重分類，重分類結(jié)果如圖2a、2b。土壤屬性數(shù)據(jù)參考世界土壤數(shù)據(jù)庫HWSD中的數(shù)據(jù)。SWAT模型根據(jù)DEM、河網(wǎng)數(shù)據(jù)，綜合選取10 000 hm2為閾值將研究區(qū)劃分為31個子流域，如圖2c所示。其中，流域出口位于第 28號子流域。然后，根據(jù)坡度、土壤類型數(shù)據(jù)、土地利用數(shù)據(jù)并分別取10%的閾值進行疊加分析，將流域劃分為171個水文響應(yīng)單元(Hydrologic response units， HRUs)，再根據(jù)氣象數(shù)據(jù)構(gòu)建天氣發(fā)生器，進行模擬計算。

(2)參數(shù)率定。SWAT模型參數(shù)共有200多個，為提高模型參數(shù)率定的工作效率，本文選取了11個參數(shù)進行率定。參數(shù)的率定基于SWAT-CUP(SWAT-Calibration and Uncertainty Programs)軟件，采用SUFI-2算法進行模型參數(shù)率定，率定結(jié)果見表5。

表5 建陽流域SWAT模型參數(shù)率定結(jié)果

3.3 模擬結(jié)果比較與分析

兩個模型率定期和檢驗期模擬結(jié)果詳見表6和圖3～4。從日徑流模擬結(jié)果可以看出，新安江模型和SWAT模型在建陽流域的模擬效果都較好，僅有個別年份的相對誤差較大。其中率定期，新安江模型的計算結(jié)果比SWAT模型的計算結(jié)果好一些，檢驗期新安江型計算結(jié)果的相對誤差略大于SWAT模型計算的相對誤差。從圖3和圖4可以看出，率定期和檢驗期兩種模型模擬與實測的徑流過程擬合較好，納什效率系數(shù)較高。率定期，新安江模型平均納什效率系數(shù)0.78，SWAT模型平均納什效率系數(shù)0.73；檢驗期，新安江模型平均納什效率系數(shù)0.79，SWAT模型平均納什效率系數(shù)0.73，均大于0.6。

表6 新安江模型與SWAT模型計算結(jié)果

4 徑流時空分布變化特征分析

4.1 徑流時空變化分析方法

由日徑流模擬結(jié)果可以看出，SWAT模型在建陽流域模擬效果較好。因此，利用SWAT模型輸出不同年份各個子流域上的產(chǎn)水量，得出研究區(qū)徑流的空間分布，并結(jié)合Sen坡度和空間變異系數(shù)Cv來研究徑流的空間分布變化。

(1)坡度Sen[29]估計法，用于確定樣本系列變化趨勢程度的大小，可用β表示趨勢變化情況。計算公式為

(4)

式中，Rj、Ri為樣本系列中任選兩個年份的徑流深數(shù)據(jù)；i和j為相應(yīng)數(shù)據(jù)的序號。

(2)空間變差系數(shù)Cv[30]，表現(xiàn)變量在空間分布上的離散程度。Cv越大，表明變量空間分布越不均勻，差異較大。計算公式為

(5)

4.2 徑流時空分布變化

計算研究區(qū)內(nèi)17個雨量站1993年～2000年的多年平均降雨量，采用反距離權(quán)重插值法，繪制多年平均降雨量空間分布圖(見圖5a)，再利用降雨插值數(shù)據(jù)進行分區(qū)統(tǒng)計，得到每個子流域的多年平均降水量(見圖5b)。

從圖5可以看出，多年平均降雨量呈現(xiàn)從西北向東南逐級遞減的態(tài)勢，多年平均降水量最大的區(qū)域位于流域西部；多年平均降水量最小的區(qū)域位于東南部。

利用SWAT模型輸出子流域逐日產(chǎn)流量，繪制各個子流域多年平均徑流量分布圖，見圖5c。對比圖5b和5c可以看出，徑流的分布規(guī)律與降雨的空間分布規(guī)律相似，空間分布差異明顯，總體上呈現(xiàn)從西北向東南遞減的趨勢。

利用不同年份的降雨量和下墊面資料模擬得到不同時期建陽流域各子流域徑流空間分布圖(見圖6)。

從圖6中可以看出，不同年份下徑流的空間分布略有不同，但整體上流域東南部的徑流深相對較小，西北部年徑流深較大。

由坡度計算可得，流域內(nèi)82.7%的區(qū)域年徑流深表現(xiàn)為較大程度的增加趨勢，僅有17.3%的區(qū)域年徑流深表現(xiàn)為減少趨勢，表現(xiàn)為減少趨勢的子流域為8號和9號子流域，9號子流域為多年平均徑流深最大的子流域，說明研究區(qū)徑流深最大值有減少表現(xiàn)，但流域整體徑流深表現(xiàn)為增加。

變差系數(shù)的變化如圖7所示。從圖7中可以看出，徑流空間變差系數(shù)呈減少趨勢，表明徑流的空間差異性在減少。

5 結(jié) 論

本文將SWAT模型與新安江模型應(yīng)用于建陽流域日徑流模擬中，并利用SWAT模型輸出結(jié)果研究徑流的時空分布變化，得到如下結(jié)論：

(1)從日徑流模擬結(jié)果來看，兩種模型在建陽流域日徑流模擬表現(xiàn)都比較好，平均納什效率系數(shù)均大于0.7，平均相對誤差均小于5%；從而表明采用這兩個模型對建陽流域的日徑流進行模擬具有一定的可靠性，可用于建陽流域徑流模擬。

(2)流域徑流空間分布與降雨空間分布相似，差異較明顯，表現(xiàn)為由西北部向東南部遞減的趨勢。

(3)多年平均徑流量最大的子流域的徑流有減少趨勢，但研究區(qū)總體徑流呈增加趨勢，空間分布差異性呈減小趨勢。