基于SWAT模型的合理亞流域劃分和徑流模擬

喻 曉，馬孝義，蔡朵朵

(西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，陜西 楊凌 712100)

0 引 言

SWAT模型作為可以進行大尺度、長時期模擬分布式水文模型，以擴展模塊的形式集成于ArcGIS軟件中，是目前解決流域水文模擬、環(huán)境評價等問題的主要手段之一[1,2]。它基于水量平衡、MUSLE方程等原理，能夠在不同土壤類型和土地利用的大尺度流域內(nèi)模擬流域產(chǎn)流、產(chǎn)沙以及營養(yǎng)物負荷等的變化情況。亞流域劃分作為構(gòu)建SWAT模型初始環(huán)節(jié)，通過給定不同的集水面積閾值，得到不同的亞流域劃分方案及模擬結(jié)果。閾值設(shè)置得越小，生成的水系越密集，劃分的亞流域數(shù)目越多，HRU越精細；但當(dāng)閾值小到一定程度后，亞流域數(shù)目太多影響運算效率，且容易出現(xiàn)程序劃分出錯，導(dǎo)致無法運行。因此，如何平衡模擬精度與模擬效率的關(guān)系、確定研究流域合理的亞流域劃分水平成為眾多學(xué)者研究的問題。Mamillapalli等[3]研究認(rèn)為SWAT模型產(chǎn)流預(yù)測隨亞流域數(shù)量的變化存在一個上限，超過這一上限，模擬結(jié)果趨于穩(wěn)定；Bingner[4]、FitzHugh[5]的研究表明SWAT模型對不同亞流域數(shù)量劃分的響應(yīng)，總是泥沙的敏感性強于徑流；張雪松等[6]認(rèn)為亞流域劃分?jǐn)?shù)量改變了模型輸入?yún)?shù)的空間集總程度，從而對流域模擬結(jié)果產(chǎn)生一定影響；胡連伍等[7]研究認(rèn)為對于豐樂河流域，徑流、泥沙、營養(yǎng)物的模擬在11～41個亞流域劃分對模型結(jié)果的影響比較穩(wěn)定。因此，不同的亞流域劃分會對模型模擬結(jié)果產(chǎn)生不同程度的響應(yīng)，針對不同流域合理的亞流域劃分方案的確定是保證模擬精度的前提。

葫蘆河流域作為北洛河最大支流，目前沒有合理亞流域劃分的相關(guān)研究，且其他流域的大多數(shù)相關(guān)研究僅是確定了研究流域亞流域劃分層次，并沒有進行模擬參數(shù)的率定和驗證，基于對徑流、泥沙、營養(yǎng)物影響分析的合理亞流域劃分方案的適用性未得到證明?；谏鲜鰡栴}，本文通過不同亞流域劃分對徑流量、泥沙負荷、營養(yǎng)物負荷影響分析的方法，旨在確定葫蘆河流域合理的亞流域劃分方案；并在此基礎(chǔ)上進行月徑流模擬參數(shù)的率定及驗證，證明其結(jié)果的可靠性，以期為葫蘆河流域的進一步的水文模擬以及泥沙、營養(yǎng)物模擬研究提供基礎(chǔ)。

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

1.1 流域概況

葫蘆河發(fā)源于甘肅省華池縣川家岔，在陜西省交口河附近匯入洛河，是北洛河最大支流。流域面積5 449 km2，干流全長235.3 km。流域地勢整體西北高東南低，河道平均比降0.24%，高程分布在766～1 748 m之間。流域年平均氣溫9.2 ℃，年蒸發(fā)量1 600～2 000 mm，多年平均徑流量1.80 億m3,多年平均降雨量539 mm，降水量年際變化大，季節(jié)分配不均，旱澇災(zāi)害頻繁，且交替出現(xiàn)。研究流域控制站點為張村驛水文站，本文研究區(qū)域即為張村驛站以上流域部分。

1.2 數(shù)據(jù)及預(yù)處理

本研究使用的原始數(shù)據(jù)包括數(shù)字高程模型(DEM)數(shù)據(jù)、土地利用數(shù)據(jù)、土壤類型數(shù)據(jù)等空間數(shù)據(jù)和氣象數(shù)據(jù)及水文數(shù)據(jù)等屬性數(shù)據(jù)，具體各類型數(shù)據(jù)的來源及用途如表1所示。

表1 SWAT建模數(shù)據(jù)及數(shù)據(jù)用途Tab.1 Data and data purpose of SWAT model

為得到滿足建模需求的數(shù)據(jù)形式，在收集到各類數(shù)據(jù)之后，需要對數(shù)據(jù)進行必要的預(yù)處理。將原始DEM數(shù)據(jù)統(tǒng)一投影WGS_1984_UTM_Zone_49N，將土地利用數(shù)據(jù)經(jīng)過統(tǒng)一投影、合并、剪切、重分類后得到SWAT模型所需的圖形，通過《中華人民共和國土壤圖》剪切得到流域土壤圖，其土壤分類系統(tǒng)為FAO-90，投影為WGS1984。

1.3 亞流域劃分方案

根據(jù)GIS軟件建議的集水面積閾值范圍(2 712～542 470hm2)按照圖1所示步驟進行不同的亞流域和HRU劃分[8]，分析其對模型產(chǎn)流、產(chǎn)沙、營養(yǎng)物負荷模擬結(jié)果的影響，以期確定葫蘆河流域合理的亞流域劃分水平使子流域個數(shù)從1開始增加，按照控制變量的原則劃分HRU時均根據(jù)土地利用、土壤和坡度所占的百分比分別為10%、5%、10%進行，具體亞流域劃分?jǐn)?shù)量見表2所示，其所得到的不同亞流域劃分如圖2所示。

采用相對誤差Re來表示由不同亞流域數(shù)量所引起的變化，其計算公式為：

圖1 GIS劃分亞流域流程圖Fig.1 Sub-watershed delineation flow diagram by GIS

圖2 不同亞流域劃分圖Fig.2 Different sub-watershed delineation

(1)

式中：Vi代表由不同亞流域劃分所得到的模擬值；Vmin代表最小集水面積閾值的模擬值。

1.4 徑流模擬

泥沙和營養(yǎng)物的遷移與徑流過程具有緊密聯(lián)系，徑流模擬是泥沙，氮、磷等營養(yǎng)物模擬的基礎(chǔ)，確保徑流模擬的準(zhǔn)確至關(guān)重要。根據(jù)確定的合理亞流域劃分構(gòu)建葫蘆河流域SWAT模型，進行月尺度下的徑流模擬。利用SWAT-CUP程序的全局敏感性分析方法(global)進行參數(shù)敏感性分析，采用PSO算法進行率定及驗證。選取張村驛水文站2007-2009年的數(shù)據(jù)為率定期，2010-2012年數(shù)據(jù)為驗證期，將2006年設(shè)置為1年的緩沖期。

表2 不同亞流域劃分參數(shù)值Tab.2 Parameter of different sub-watershed delineation

選用決定系數(shù)(R2)和納什系數(shù)(Ens)評價模擬過程與實測過程之間的擬合程度，其計算公式分別為：

(3)

2 結(jié)果與討論

2.1 合理亞流域劃分

分別對不同亞流域劃分下模型徑流、泥沙、營養(yǎng)物模擬結(jié)果進行分析，確定葫蘆河流域合理的亞流域劃分。

(1)亞流域劃分對徑流的影響。不同亞流域劃分引起的月均徑流量的變化如圖3所示。對葫蘆河流域來說，產(chǎn)流量在亞流域數(shù)量為1～17個和17～37個兩個變化過程時都有一個緩慢增加又減小的階段，產(chǎn)流在37個亞流域劃分后處于穩(wěn)定，繼續(xù)細分亞流域數(shù)量相對誤差均在1%以內(nèi)，說明對其再進行更詳細的劃分并不能明顯提高徑流模擬精度。因此，存在使流域徑流模擬趨于穩(wěn)定的亞流域劃分水平，這與李曼曼[10]、張召喜[11]、馬放[12]的研究結(jié)論一致。

圖3 月均徑流量隨亞流域數(shù)量的變化Fig.3 The variation of the average monthly runoff with the number of sub-watershed

模型選用SCS徑流曲線數(shù)法進行產(chǎn)流模擬，決定產(chǎn)流量的敏感性參數(shù)為徑流曲線數(shù)CN，CN值的大小變化直接影響了產(chǎn)流量的變化。為此我們將研究不同亞流域劃分與流域面積加權(quán)CN值的關(guān)系進行分析，如圖4所示，發(fā)現(xiàn)在亞流域數(shù)量從1～95個的過程中，CN的取值變化程序很小，保持在58.84～59.35。由此可知亞流域的不同劃分并不能引起CN值的劇烈變化，因而引起的徑流變化也相對較小。

圖4 CN值隨亞流域數(shù)量的變化Fig.4 The variation of CN with the number of sub-watershed

(2)亞流域劃分對泥沙負荷的影響。不同亞流域劃分引起月均泥沙負荷變化如圖5所示。流域亞流域數(shù)量由1增加至29個的過程中，產(chǎn)沙量呈現(xiàn)先顯著增加后減小的趨勢，在11個時處于最大值，37個亞流域劃分后產(chǎn)沙量趨于穩(wěn)定。

圖5 月均泥沙負荷隨亞流域數(shù)量的變化Fig.5 The variation of the average monthly sediment with the number of sub-watershed

SWAT模型模擬泥沙負荷變化受亞流域數(shù)量的影響，分析其原因有兩點：一是亞流域的不同劃分引起河網(wǎng)密度的變化，隨著亞流域數(shù)量的增加，生成的河網(wǎng)變得簡單化、河網(wǎng)密度減小，定義的河段以及所控制的流域面積被新生成的河網(wǎng)代替，由此模擬的泥沙演算發(fā)生變化；二是SWAT模型模擬泥沙生成量采用修正的通用土壤流失方程(MUSLE方程)[13]：

msed=11.8 (QsurfqpeakAhru)0.56

KUSLECUSLEPUSLELSUSLE

(4)

式中：msed為土壤侵蝕量，t；Qsurf為地表徑流，mm/h；qpeak為洪峰流量，m3/s；Ahru為水文響應(yīng)單元(HRU)的面積，hm2；KUSLE為土壤侵蝕因子；CUSLE為植被覆蓋和管理因子；PUSLE為保持措施因子；LSUSLE為地形因子；CFRG為粗碎屑因子。

MUSLE方程中PUSLE對所有的HRU是常量[6]。因此主要探討KUSLE、CUSLE、LSUSLE對亞流域劃分的影響，如圖6所示。地形因子LS作為一個復(fù)合因子是將坡度因子S和坡長因子L結(jié)合的一個參數(shù)表達，圖6(c)表示了亞流域劃分對坡面坡長和坡面坡度的影響。由圖6(a)，圖6(b)可以看出，KUSLE、CUSLE的值幾乎不受亞流域數(shù)量變化的影響，坡度和坡長隨亞流域數(shù)量的增加分別呈現(xiàn)增大和減小的趨勢，因而引起泥沙負荷發(fā)生變化。

圖6 土壤侵蝕因子、植被覆蓋和管理因子、坡面坡長、坡度隨亞流域數(shù)量的變化Fig.6 The variation of KUSLE、CUSLE、overland slope and slope length with the number of sub-watershed

(3)亞流域劃分對營養(yǎng)物負荷的影響。隨著亞流域劃分?jǐn)?shù)目的變化，各營養(yǎng)物的負荷發(fā)生明顯變化，如圖7所示。營養(yǎng)物的變化隨著亞流域數(shù)量的增加均先增大后減小且最后在37個亞流域劃分后趨于穩(wěn)定。硝氮與徑流類似，7個亞流域劃分后相對誤差均小于10%，隨亞流域數(shù)量變化幅度不大，這可能是由于硝氮主要隨地表徑流等遷移的緣故。有機氮、有機磷、總氮、總磷和泥沙的變化趨勢基本一致，在亞流域數(shù)量為1～11個時逐漸增加，而后出現(xiàn)降低趨勢，且37個亞流域劃分后趨于穩(wěn)定，分析其原因認(rèn)為有機態(tài)氮、磷主要是吸附在泥沙顆粒上進行運移，泥沙的流失量變化一定程度上可以反映有機態(tài)氮、磷負荷的變化，總氮、總磷的變化受有機態(tài)氮、磷的影響較大。與Wang G[14]等一些學(xué)者的研究類似，營養(yǎng)物負荷變化趨勢表現(xiàn)出與徑流或泥沙相似。

圖7 月均有機氮、有機磷、硝氮、氨氮、總氮、總磷負荷隨亞流域數(shù)量的變化Fig.7 The variation of organic nitrogen, organic phosphorus, nitrate nitrogen, ammonia nitrogen, total nitrogen and total phosphorus with the number of sub-watershed

根據(jù)本文研究的不同亞流域劃分情況下徑流、泥沙、營養(yǎng)物的變化情況可知，葫蘆河流域最小集水面積閾值為10 000 hm2時，既有較高的模擬精度又能滿足模擬效率，此時亞流域劃分?jǐn)?shù)量為37個左右。

2.2 模型徑流模擬

(1)參數(shù)敏感性評價?；?7個合理的亞流域數(shù)劃分，建立月尺度的SWAT模型，在SWAT-CUP程序中進行參數(shù)敏感性分析，最終選取15個敏感性較高的參數(shù)，結(jié)果如表3所示。

由表3可知，SCS徑流曲線數(shù)(CN2)是葫蘆河流SWAT模型最為敏感的參數(shù)，主要影響徑流大小。土壤容重(SOL_BD)、冠層最大出水量(CANMX)等參數(shù)也是較強的敏感性，在調(diào)參時應(yīng)予以重點考慮。

(2)模型率定。根據(jù)參數(shù)敏感性分析結(jié)果，運用PSO算法進行率定期徑流數(shù)據(jù)的迭代計算，本文進行4次迭代，每次迭代模擬500次，率定期的模擬值與實測值基本吻合，比較準(zhǔn)確地捕捉到峰值，如圖8所示。率定期決定系數(shù)R2為0.84，納什系數(shù)Ens為0.76。由圖8可看出，除個別月份的峰值模擬有偏差外，實測徑流過程和模擬徑流過程起伏基本一致，說明SWAT模型的模擬結(jié)果可代表葫蘆河流域2007-2009年的徑流變化規(guī)律。

(3)模型驗證。在模型率定的基礎(chǔ)上，對2010-2012年的徑流進行驗證，如圖9所示。由圖9可以看出根據(jù)率定期參數(shù)進行驗證期徑流過程模擬，模擬值與實測值的起伏變化基本吻合，決定系數(shù)R2和納什系數(shù)Ens分別為0.81和0.62，驗證期模擬結(jié)果符合模擬精度，表明基于本文方法確定的亞流域劃分可以用于葫蘆河流域的SWAT模型的徑流模擬。

表3 參數(shù)敏感性分析結(jié)果Tab.3 Result of parameter sensitivity analysis

圖8 率定期模擬結(jié)果Fig.8 The simulation result on calibration

圖9 驗證期模擬結(jié)果Fig.9 The simulation result on validation

3 結(jié) 語

本文應(yīng)用SWAT2012，分析了葫蘆河流域不同亞流域劃分對模型結(jié)果的影響，主要得到以下結(jié)論。

(1)葫蘆河流域最粗略亞流域劃分到較精細劃分的過程，產(chǎn)流、產(chǎn)沙以及營養(yǎng)物負荷的變化均呈現(xiàn)先顯著上升后下降，最后趨于穩(wěn)定的變化規(guī)律，說明較低的亞流域劃分水平對模型結(jié)果是不穩(wěn)定的，綜合模擬精度和模擬效率，確定葫蘆河流域合理的亞流域劃分為37個。

(2)葫蘆河流域亞流域劃分對徑流影響的敏感性要低于泥沙和營養(yǎng)物，徑流與硝氮的變化幅度相差不大，7個亞流域劃分后相對誤差均小于10%，泥沙與有機氮、有機磷、氨氮、總氮、總磷的變化趨勢基本一致。產(chǎn)流變化主要受CN值影響，坡面長度和坡面坡度的變化對泥沙的影響較大，硝氮和有機氮、磷分別與徑流和泥沙的變化趨勢基本一致，總氮、總磷的變化主要受有機氮、有機磷的變化影響。

(3)基于37個亞流域劃分的徑流模擬，滿足模擬精度要求，通過分析不同亞流域劃分對徑流、泥沙、營養(yǎng)物的不同響應(yīng)確定亞流域劃分的方法可行。SWAT模型在葫蘆河流域有良好的適用性，可為今后葫蘆河流域進一步的水文模擬研究提供參考。

□

[1] 徐宗學(xué)，程 磊.分布式水文模型研究與應(yīng)用進展[J]. 水利學(xué)報，2010，41(9):1 009-1 017.

[2] 江凈超，朱阿興，秦承志，等.分布式水文模型軟件系統(tǒng)研究綜述[J]. 地理科學(xué)進展，2014，33(8):1 090-1 100.

[3] Mamillapalli S, Srinivasan R, Arnold J G, et al. Effect of spatial variability on basin scale modeling[EB/OL]. http:∥www.ncgia.ucsb.edu., 1996.

[4] Binger R L. Effect of watershed subdivision on simulation runoff and fine sediment yield [J]. Transaction of the ASAE, 1997,40(5):1 329-1 335.

[5] FitzHugh T W, Mackay D S. Impacts of subwatershed partitioning in modeled source-and transport-limitied sediment yields in an agricultural nonpoint source pollution model [J]. Journal of Soil and Water Conservation, 2001,56(2):137-147.

[6] 張雪松，郝芳華，程紅光，等. 亞流域劃分對分布式水文模型模擬結(jié)果的影響[J].水利學(xué)報，2004，(7):119-123.

[7] 胡連伍，王學(xué)軍，羅定貴，等. 不同子流域劃分對徑流、泥沙、營養(yǎng)物模擬的影響-豐樂河流域個例研究[J]. 水科學(xué)進展，2007，18(2):235-240.

[8] 黎云云，暢建霞，金文婷. 基于SWAT模型的渭河流域分區(qū)徑流模擬研究[J]. 西北農(nóng)林科技大學(xué)學(xué)報(自然科學(xué)版)，2017，45(4):1-9.

[9] D N Moriasi, J G Arnold, M W Van Liew, et, al. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations[J]. Transaction of the ASAE, 2007,50(3):885-900.

[10] 李曼曼，韓會玲，劉曉英，等. SWAT模型最佳子流域劃分方案研究-以云南洱海流域為例[J].中國農(nóng)業(yè)氣象，2012，33(2)：185-189.

[11] 張召喜，羅春燕，張敬鎖，等.子流域劃分對農(nóng)業(yè)面源污染模擬結(jié)果的影響[J].農(nóng)業(yè)環(huán)境科學(xué)學(xué)報，2012，31(10)：1 986-1 993.

[12] 馬 放，姜曉峰，王 立，等.基于SWAT模型的亞流域劃分方法研究[J].中國給水排水，2015，31(7)：53-57.

[13] Williams J R. Sediment routing for agricultural watersheds[J]. Water Resources Bulletin, 1975, 11(5): 965-974.

[14] Wang G, Chen L, Huang Q, et, al. The influence od watershed subdivision level on model assessment and identification of non-point source priority management areas[J]. Ecological Engineering, 2016,87(3):110-119.