武海霞，李清雪，孫玉壯，冀雅珍

(1.河北省水資源高效利用工程技術研究中心，河北工程大學，河北 邯鄲 056038；2.河北省煤炭資源綜合開發(fā)與利用協(xié)同創(chuàng)新中心, 河北省資源勘探研究重點實驗室, 河北工程大學，河北 邯鄲 056038；3.山西省水利水電科學研究院，太原 030002)

0 引 言

目前，SWAT分布式水文模型在流域水文、生態(tài)和環(huán)境等領域得到廣泛應用。然而，模型變量和參數具有很大的隨機性和時空變異性，且這些變量和參數的變化會造成模擬結果較大的差異[1]。流域離散是分布式水文模型的基礎，子流域的數量和面積對模擬結果影響差異極為顯著，子流域離散方案產生的模擬偏差，可能會超過模型參數率定所產生的偏差，而且在地形復雜、坡度較大的流域，不同子流域數量使模擬精度存在不確定性[2]。關于子流域劃分對流域模擬結果的影響，國內外學者進行了相關研究。一些學者認為子流域劃分對流域產流模擬影響較小,而對產沙模擬影響較大[3,4]。郝芳華等和張雪松等的研究表明，當子流域劃分數量達到一定水平時，增加子流域劃分對模型輸出結果的影響較小[5,6]。部分學者通過進一步分析，認為在不同集水面積閾值條件下其模擬結果存在閾值適用范圍的上下限。Tripathi M和Jha M認為對流量、泥沙、NO3-N 和無機磷模擬的最適子流域面積占流域總面積的比例分別為4%、3%、2%和5%[7,8]。同時，基于SWAT模型的子流域劃分方案對模擬結果的影響和研究尺度有一定的關系。Soni等的研究表明SWAT模型的產流預測對子流域和水文響應單元的數量具有一定響應，當研究的尺度太大時，流域參數的集總會影響徑流曲線CN(Curve Number)值的分布，從而導致預測的不準確[9]。雖然國內外學者針對分布式水文模型子流域劃分方案進行了相關研究，但由于子流域的數量和面積對模擬結果影響差異極為顯著，模擬結果的不確定因素比較復雜，對于如何準確而有效地確定子流域劃分數量并沒有統(tǒng)一的結論和標準，而且在不同流域，結論并不完全一致。本次以大尺度復雜山區(qū)漳河上游流域為研究對象，分析了在不同水平年流域劃分方案對SWAT模型模擬結果的影響，對準確模擬該地區(qū)的水沙情況奠定一定基礎，可為該地區(qū)水資源水環(huán)境的進一步研究提供技術支撐。

1 流域概況及模型模擬

1.1 流域概況

本次研究選擇漳河上游流域作為研究對象,漳河上游位于海河流域西南部,具體位置見圖1。流域面積約17 730 km2，占海河流域總面積的6%。流域屬溫帶大陸性季風型氣候，四季分明。全年夏短冬長，多年平均氣溫7.4～10.3 ℃，無霜期170 d左右。流域多年平均降水量為576.9 mm，降水量總的分布趨勢是南部大于北部，迎風坡多于背風坡和盆地中心。流域內支流眾多，水系呈扇形分布，分東西兩大支流。東支為清漳，該地區(qū)為石質山區(qū)，山高谷深，巖石裸露，坡陡流急，含沙量?。晃髦闈嵴模搮^(qū)域為山丘盆地區(qū)，盆地內黃土覆蓋較厚，植被較差，水土流失嚴重，洪水挾帶泥沙較多。

漳河上游流域內以山區(qū)地貌為主，地形落差高達2 000多米，東鄰太行山脈，西鄰太岳山脈，中間是長治盆地。河流橫穿太行山脈，大地構造處于太行山新華夏系第三隆起帶，由一系列北東向平緩復式群皺組成，巖性地層出露主要有震旦系石英砂巖、泥巖，寒武系紫紅色頁巖，奧陶系石灰?guī)r，河谷兩巖階地由第四系積物組成。

圖1 漳河上游流域位置

1.2 模型數據和參數

模型所需的空間數據主要包括數字高程模型(DEM)、土壤數據、土地利用類型數據和氣象數據。DEM使用國際科學數據服務網站提供的分辨率為30 m的數據，依照研究區(qū)經緯度范圍下載，將多景的數據進行融合與拼接，使用數字化的研究區(qū)邊界裁切出邊界精準的漳河上游流域的DEM圖。土壤數據來源于世界糧農組織網站的HWSD(Harmonized World Soil Database)數據。數據格式為柵格數據，分辨率為1 000 m，投影為WGS1984，土壤分類系統(tǒng)是FAO90。由于模型采用的是美國土壤分類系統(tǒng)，需要重新建立研究區(qū)的土壤屬性數據庫。土地利用類型數據由海河流域南部 1∶100 萬土地利用圖片數字化得到，并根據土地利用的一級分類，轉化為 SWAT 模型能識別的模型代碼，最終劃分為 5種土地利用類型：耕地、林地、草地、水域和其他用地。河網水系圖來源于國家基礎地理中心提供的分辨率為30 m的數據。降雨、溫度、相對濕度、風速及太陽輻射等氣象數據來源于SWAT 官方網站及中國氣象科學數據共享服務網。模型所需基礎數據的詳細描述見表1。利用Arcgis軟件對DEM、土壤類型圖和土地利用圖作Albers等積圓錐投影統(tǒng)一地理坐標系WGS1984，并統(tǒng)一生成相同的柵格大小。

表1 模型所需基礎數據

根據郝芳華[5]等的研究，由于現有的確定模型參數的最優(yōu)化方法對實測資料的依賴性很大，只能反映模擬值與實測值的擬合程度，而不能揭示參數的物理意義，因此本次直接采用模型默認的參數值進行模擬。

1.3 子流域劃分方案

本次研究基于最小集水面積( Critical Source Area，CSA)閾值劃分出不同數量的子流域。最小面積閾值(CSA)是形成河流的最小集水面積，也可以表示為其與流域總面積的百分比，CSA越大子流域劃分數量越少。

為深入分析流域離散化程度對SWAT模型模擬結果的影響，本次研究在前人成果的基礎上，結合研究區(qū)特點選擇6種子流域劃分方法，即CSA分別為0.4%、0.5%、0.7%、1%、2%和10%，并通過確定土地利用和土壤類型面積閾值分別為 5% 和10%來確定子流域中 HRU 的數量。表 2 列出了6種流域離散化方案的特征參數。

表2 子流域劃分及流域參數

根據CSA值對漳河上游流域進行子流域劃分,獲得6種不同的子流域劃分結果見圖2。

圖2 流域不同子流域劃分層次圖

定義相對偏差為：

Re=(Vi-Vmin)/Vmin×100%

式中：Vi為不同子流域劃分模擬得到的模擬值；Vmin為子流域數目最小時得到的模擬結果。

本次研究選擇的實測資料來自于課題組其他成員的分析成果。王喆[10]利用漳河上游近40年的氣象資料進行面降雨量的趨勢性分析，得出1979-1995年為枯水年，1996年為豐水年，1997-2013年為平水年。本次研究在時間尺度上選擇1995、1996年和1997年的氣象資料作為枯水年，豐水年和平水年的輸入數據進行模擬結果的分析。

子流域劃分方案下年均徑流、泥沙及營養(yǎng)物在平水年的相對偏差(Re)見表3。

表3 不同子流域劃分下年徑流、泥沙及營養(yǎng)物在平水年的相對偏差(Re)

注：Re1、Re2、Re3、Re4分別為年均徑流、泥沙、有機氮和有機磷相對偏差值。

2 結果與討論

2.1 子流域劃分方案對地表徑流的影響

由圖2可看出，不同子流域劃分結果最直接的表現是各子流域內面積大小的變化及河網詳實程度的描述。由表2可看出，子流域劃分方案對徑流模擬結果有一定的影響，最大偏差值為11.23%，但和泥沙負荷相比而言，子流域劃分方案對徑流的影響并不大，尤其在子流域個數為58～111范圍內，影響很小，這與前人的研究結論一致。這主要是因為本文所選用的產流方程為徑流曲線數法，其中決定產流量的敏感性參數為徑流曲線數CN，而不同的劃分方法對流域面積加權CN值的影響比較小。

圖3 不同水平年徑流量隨子流域數量變化曲線

圖3顯示了在不同水平年徑流量隨子流域數量變化的趨勢，總體表現為先增大，然后趨于穩(wěn)定，最后呈減小趨勢。具體表現為當子流域數量從5增加到25時，年徑流量隨之增大，豐水年增加的幅度最大，而枯水年呈略微減小趨勢，表明這兩者之間的子流域劃分并不是最佳的子流域劃分方案。當子流域劃分數量進一步增加，子流域個數大于25，小于111，CSA值為流域總面積的0.5%～2%時，年徑流量趨于穩(wěn)定，且不同水平年趨于穩(wěn)定的范圍一致。當子流域劃分數量超過111個時，徑流量明顯減少，尤其在豐水年減小更明顯，而枯水年呈略微的增加趨勢。這說明在不同水平年徑流量隨子流域劃分數量變化的趨勢并不完全一致。造成這種情況的原因可能是過多的子流域劃分，形成了較多的過于狹長的子流域，使得模型出現了對水文過程的不真實模擬[11]?？梢?，在漳河上游流域，合理的子流域劃分水平對徑流模擬的影響并不大，但存在一個使徑流模擬趨于穩(wěn)定的子流域劃分的上下限，當低于下限，高于上限時，模擬結果受劃分數量的影響較大。

2.2 子流域劃分方案對泥沙的影響

由表3可看出，泥沙負荷對子流域劃分方案的敏感性強于徑流量，最大偏差值高達30.31%。由圖4可看出，在不同水平年泥沙負荷隨子流域劃分數量的增加逐步呈增加趨勢，這和年均徑流量的變化趨勢并不完全一致。當子流域數量由5變到25時,泥沙負荷有增加趨勢，當子流域個數在25～111范圍時，泥沙負荷趨于穩(wěn)定，這和年徑流量趨于穩(wěn)定的范圍一致。但泥沙和年均徑流在不同水平年隨子流域數量變化的趨勢并不完全相同，當子流域個數超過模擬結果趨于穩(wěn)定的上限值111個時，年均徑流量明顯減少，而泥沙量明顯增加，尤其在豐水年表現的更明顯。這可能是因為子流域個數改變引起流域內高程的空間變化，從而影響到流域的坡面坡長和坡度，而坡度和坡長參數(LS 因子)對于泥沙量是敏感因子，能顯著影響泥沙產量的預測。泥沙演算中的泥沙沉積和降解過程會影響產沙預測的結果，隨著子流域劃分數目的減小，流域描述簡化，河網密度減小，從而影響泥沙演算，并使預測的準確性降低[12]。

圖4 不同水平年泥沙隨子流域數量變化曲線

2.3 子流域劃分方案對營養(yǎng)物的影響

營養(yǎng)物氮、磷的變化原因比較復雜，除了受徑流、沉積負荷的影響外，還受到自身復雜計算算法的影響，而有機態(tài)養(yǎng)分的遷移主要是通過吸附在泥沙等固體上進行的。

隨著子流域劃分數目的變化，營養(yǎng)物有機氮和有機磷負荷發(fā)生明顯變化。由表3可看出，最粗略的子流域劃分和最精細劃分之間的相對誤差，最大偏差值分別為26.53%和45.30%，而且有機磷的相對偏差值總體都比有機氮的大。 由此可見，子流域劃分對營養(yǎng)物負荷的影響較大，這和前人的研究結果一致。

由圖5、6可看出，在不同水平年，隨著子流域劃分數量的增加，流域出口營養(yǎng)物有機氮、磷總體上呈現增加的趨勢，和泥沙負荷隨子流域劃分數量的變化趨勢類似，營養(yǎng)物負荷隨子流域數目變化趨于穩(wěn)定的范圍和泥沙負荷的范圍一致，都在25～111。

圖5 不同水平年有機氮隨子流域數量變化曲線

圖6 不同水平年有機磷隨子流域數量變化曲線

有機氮、磷負荷與泥沙負荷是成比例的，但是由于有機氮磷的河道傳輸與泥沙的河道傳輸算法的不同，所以出口有機氮磷與泥沙并不一定表現出共同的趨勢。徑流量和泥沙負荷都隨子流域劃分方案的不同而產生一定變化，且各自具有不同的變化趨勢和變化幅度(圖 3、圖 4)，因此營養(yǎng)物負荷量的變化原因很難歸結為某一簡單的原因。

2.4 最佳子流域劃分方案的確定

對于子流域劃分方案對SWAT模擬結果的影響，多數學者認為對徑流量的影響較小，對泥沙和營養(yǎng)物的影響較大，但是如何準確有效地確定子流域劃分數目至今并沒有統(tǒng)一的標準。一般來講，地形復雜的山區(qū)，由于地形破碎，空間差異大，需要更詳細的子流域劃分，以此充分體現該地區(qū)局部空間的差異。但CSA值過小，劃分過細時會產生較多過于狹長或過小的虛假子流域，導致不真實的水文模擬，而當亞流域劃分數量較少時，對流域描述不夠充分導致模擬輸出結果不穩(wěn)定，難以達到理想的預測精度。本次研究在不同子流域劃分方案下的模擬結果分析表明，在漳河上游流域，無論是徑流、泥沙、還是營養(yǎng)物有機氮磷，對子流域劃分方案趨于穩(wěn)定的范圍一致，都在25～111之間。當子流域的數量從5個增加到25個時，年均徑流、泥沙和有機氮磷相對來說變化幅度較大，且在不同水平年表現出的變化趨勢并不一致，這表明此區(qū)間并不是最佳的子流域劃分方案。隨著子流域劃分數量的進一步增加，年均徑流、泥沙和有機氮磷逐漸趨于穩(wěn)定，當超過某一子流域劃分方案時，劃分的HRU數量急劇增加，計算單元劃分過于精細，可能導致計算結果失真，同時花費更長的計算時間，降低計算效率。因此，在漳河上游流域，存在上下限閾值111、25，當子流域個數低于下限閾值，對空間的描述不夠，難以滿足流域規(guī)劃與管理的要求，高于上限閾值，模擬結果出現失真。通過綜合權衡下墊面特征、與其他數據的協(xié)調及模型運行能力等因素的基礎上，根據模型模擬結果，最后確定漳河上游流域最佳子流域劃分個數為79個，CSA取值為12 411 hm2，僅占流域總面積的0.7%，這說明漳河上游流域內復雜的地理空間環(huán)境需要更加詳細的描述，這與實際需求相符。本次研究結果中子流域平均劃分面積占整個流域面積的比值與以往的研究結論并不一致，說明利用子流域面積閾值設定不同的子流域劃分方案存在地域差異，需要針對不同流域的實際地形、地貌、海拔等空間屬性對模型進行初步的徑流、泥沙和營養(yǎng)物質的敏感性分析，合理確定CSA的取值及子流域的劃分方案，從而減少分布式水文模型模擬的不確定性，為進一步更準確地模擬水文過程奠定基礎。

3 結 語

本文以中國海河流域西南部的漳河上游流域為研究對象，建立流域SWAT年尺度產流產沙模型，進一步分析不同CSA取值下，子流域劃分數量對徑流、泥沙及營養(yǎng)物模擬的影響。主要結論如下：

(1)從總體來看, 子流域劃分層次對模擬結果的影響存在上限下限兩個閾值25、111，介于兩個閾值之間模擬結果相對穩(wěn)定，超出上限閾值111，模擬結果出現失真；低于下限閾值25，對空間的描述不夠，難以滿足流域規(guī)劃與管理的要求。子流域劃分數量對年均徑流量影響較小，對泥沙和營養(yǎng)物影響較大。

(2)在不同水平年，模擬結果對子流域劃分數量的敏感性不一致，在豐水年最敏感，在枯水年最不敏感，建議在徑流量較大的地區(qū)，要合理選擇子流域劃分個數。

(3)綜合考慮模擬精度和模擬效率，漳河上游流域 SWAT 模擬的最合理子流域劃分水平為79個，CSA值為12 411 hm2，占整個流域面積的0.7%。

□

[1] Manillapalli S. Effect of spatial variability on basin scale modeling[C]∥Third International Conference Workshop on Integrating GIS and Environmental Modeling. Proceedings, Santa Fe, New Mexico, 1996.

[2] Kim J G, Park Y, Yoo D, et al. Development of a SWAT patch for better estimation of sediment yield in steep sloping watersheds [J]. Journal of the American Water Resources Association，2009，45(4):963-972.

[3] FitzHugh TW,Mackay D S. Impacts of input parameter spatial aggregation on an agricultural nonpoint source pollution model[J]. Journal of Hydrology, 2000,(236):35- 53.

[4] 盧文喜，伊燕平，張 蕾，等． 不同亞流域劃分數量對SWAT 模型模擬結果的影響[J]．水電能源科學，2010，28(10):23-25.

[5] 郝芳華，張雪松，程紅光，等．分布式水文模型亞流域合理劃分水平芻議[J]．水土保持學報，2003，12(4)：75-77．

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

[7] Tripathi M, Raghuwanshi N, Rao G. Effect of watershed subdivision on simulation of water balance components[J].Hydrological Processes,2005,20(5):1 137-1 151.

[8] Jha M, Gassman P W, Secchi S, et al. Effect of watershed subdivision on swat flow, sediment, and nutrient predictions[J].Journal of the American Water Resources Association,2007,40( 3):811-825.

[9] Soni M. Pradhanang, Russell D. B. Effects of critical source area on sediment yield and streamflow [J].Water and Environment Journal, 2014,28:222-232.

[10] 王 喆．氣候變化和人類活動對漳河上游地表水資源的影響研究[D]．河北邯鄲：河北工程大學城市建設學院，2015．

[11] 耿潤哲，李明濤，王曉燕，等. 基于SWAT模型的流域土地利用格局變化對面源污染的影響[J]．農業(yè)工程學報，2015,31(16):241-250.

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