基于SWAT模型的張家口清水河流域土地利用情景變化對徑流影響研究

王 磊, 劉亭亭, 謝建治

(1.北京市農林科學院 植物營養與資源研究所, 北京 100097; 2.河北農業大學 資源與環境學院, 河北 保定 071000)

1990年前后人們開始認識到土地利用/覆被變化對于研究全球變化特征與可持續發展戰略上占據了舉足輕重的地位[1]，因此“國際地圈與生物圈計劃”(IGBP)[2]與“全球環境變化的人文因素計劃”(IHDP)在1995年首先共同提出了土地利用/覆被變化(land use and land cover change, LUCC)這一跨越性的科學研究項目[3]，而土地利用變化下的水文響應研究一直是全世界的焦點問題之一。

大量的研究表明土地利用變化可直接導致水文要素發生改變。Santillan等[4]對菲律賓南部流域研究發現土地利用/覆被變化與地表徑流具有直接關系。楊宏偉等[5]的研究結果表明林地面積與年徑流量呈負相關關系；一年內枯水期的徑流量大，豐水期徑流量小。Cuo等[6]研究結果表明流域內建設用地的大量增加導致年均徑流量、總徑流量和階段性徑流量都發生了變化。Guo等[7]研究果表明，鄱陽湖流域林地面積的大量增加導致雨季流量減小，旱季流量增大，從而減少了干旱、洪澇災害的發生。

目前，利用3S(Remote Senescing，Geographical Information System，Global Positioning System)技術支持的非點源污染模型可較好地評估土地利用變化對流域徑流的影響，其中SWAT(Soil and Water Assessment Tool)模型的有效性已通過大量的研究證明。Nosetto等[8]使用SWAT模型分析南半球水循環的影響因素，發現土地利用變化是水循環影響的重要因子。Prokop等[9]發現引起喜馬拉雅山1930—2010年水文變化的重要因素是土地利用變化。羅巧等[10]應用SWAT模型研究不同土地利用方式對徑流的影響程度，結果表明林地和草地面積與徑流成反比，而耕地和建設用地面積與徑流成正比。王學等[11]在白馬河流域應用SWAT模型分析了不同土地利用/覆被情景下的水文響應單元，結果表明不同的土地利用類型對流域產流的影響不同。陳引珍等[12]應用SWAT 模型對長江上游清港河流域進行了研究，設置了兩種土地利用情景，得出年徑流的變化與年降雨量的關系并不密切，而平均月徑流的變化，則隨著降雨量的不斷增大而增大。李勇[13]基于SWAT模型分析計算了汾河水庫流域不同時期土地利用變化對徑流的貢獻率，結果表明汾河水庫流域土地利用變化對汾河徑流有一定的影響，雖然影響程度不大，卻是徑流變化的一個不可忽略的影響因素。

清水河是洋河的重要支流，上游地區位于河北省張家口市山西段造山運動的沉降帶。本文選取張家口清水河上游流域作為研究區域，建立SWAT模型，就該模型在研究區域的適用性進行評價，分析其主要土地利用類型變化對徑流的影響系數，了解張家口清水河流域土地利用—徑流機制，為該地區水資源保護利用和水土保持工作提供科學依據。

1 研究區概況

考慮到模型所需資料及流域監測數據的可獲得性，本文將研究區域集中在張家口市崇禮區清水河流域上游。地理位置位于東徑114°27′—115°30′，北緯40°46′—41°17′。東西長64 km，南北寬57 km，總面積為2 334.1 km2。崇禮區是2022年冬季奧林匹克運動會主辦城市之一，研究區內將有數座為冬奧會服務的大型室外滑雪場落成并投入使用。氣候屬東亞大陸性季風氣候中溫帶亞干旱區。區域內土壤類型一般從山地到河谷依次為：山區上部為栗鈣土、下部為淡栗鈣土、山腳溝谷(包括洪積扇、殘積臺地)為淡栗鈣土性土。植被資源主要分為暖溫帶落葉闊葉林區域和溫帶草原區域兩種。該區域面臨水資源匱乏，西溝區域土壤侵蝕嚴重等問題。

2 研究方法

2.1 SWAT模型模擬法

SWAT模型是在30 a前由美國農業部(USDA)和農業研究中心(ARS)針對較大尺度流域開發的具有很強物理機制的長時間的分布式水文模型[14]。SWAT模型綜合考慮了自然因素和社會因素，能夠較好地模擬地表徑流的過程[15]，其吸取了CREAMS(Chemic-als Runoff,and Erosion from Agricultural Management Systems)，EPIC(Environmental Policy-Integrated Climate)等模型的主要特征，主要形成于SWRRB(Simulator for Water Resources in Rural Basins)[16]。主要運用在基礎數據缺乏的區域，通過模擬預測流域內水、泥沙及農業管理措施所產生的響應分析，用于幫助進行水資源管理[17]。

SWAT模型主要由水文過程子模型、土壤侵蝕子模型(Epic)、污染負荷子模型3部分組成，其中水文過程子模型包括水文循環陸地階段、水文循環演算和水文循環的匯流階段[18]。

水文主要包括降水、地下徑流、入滲、疏散損失等過程[19]，SWAT模型采用的水量平衡方程[20]如下：

(1)

式中：HStn為第n天的土壤含水量(mm)；HSn為第n天土壤的最初含水量(mm)；t為時間(d)；Rn為第n天的降雨量(mm)；Qn為第n天的地表徑流量(mm)；En為第n天的蒸發蒸騰量(mm)；Wn為第n天土壤剖面底層的滲透量和側流量(mm)；Tn為第n天的地下水回流量(mm)。

SWAT模型采用修正之后的MUSLE模型來模擬由降水和徑流產生的土壤侵蝕[21]。方程如下：

SED=11.8(Qsurf·Qpeak·AREA)0.56·K·C·P·LS·CFPG

(2)

式中：SED表示產沙量(t)；Qsurf表示地表徑流總量(mm/hm2)；Qpeak表示洪峰流量(m3/s)；AREA表示HRU面積(hm2)；K表示水土流失方程中土壤可侵蝕因子[0.013 (t·m2·h)/(m3·t·cm)]；C表示土地覆蓋與管理措施因子；P表示水土保持措施因子；LS表示地形因子；CFRG為粗糙因子。

本文選取決定性系數R2和納什效率系數NSE來評價模型的模擬精度。

決定性系數R2表示模型模擬值與實測值的變化趨勢。計算公式[22]如下：

(3)

納什效率系數NSE表示模型的總體效率，NSE越高，模型的可信度越高，計算公式[23]如下：

(4)

2.2 單一土地利用類型動態函數

本文采用單一土地利用動態度來表示土地利用類型變化速率。單一土地利用類型動態度是指在一定時間范圍內研究區內一種土地利用類型速度變化，計算公式為：

(5)

式中：K表示研究時段內某一土地利用類型動態度；Ua表示研究期初某一種土地利用類型的數量；Ub表示研究期末某一種土地利用類型的數量；T表示研究時段長，當T的時段設定為年時，K值為該研究區某種土地利用類型年變化率。

2.3 多目標決策函數法

通過MATLAB(matrix & laboratory)軟件建立多元一次線性方程組求解每個單元的某一土地利用類型對流域徑流的影響[26]。

年平均徑流量是指時間步長內不同土地利用類型對于徑流量影響之和，不同土地利用類型對于徑流量的影響將其量化，可正可負，但意義不同，當為正時，代表土地利用類型與徑流量變化呈正相關，為負時代表土地利用類型與徑流量變化方向呈負相關。

年平均徑流量Xn的計算方程如下：

Xn=Ma×An+Mb×Fn+Mc×Gn

(6)

式中：Xn為年平均徑流量；An為林地面積；Fn為草地面積；Gn為耕地面積；Ma，Mb，Mc分別為林地、草地、耕地土地利用類型單位面積(km2)影響徑流量，單位均為m3/(km2·h)；n為不同年份下的土地利用情況。

將不同土地利用分布情景下的徑流量帶入建立如下矩陣，計算可得到不同土地利用類型對徑流的影響系數。

(7)

3 模型應用

3.1 流域基礎數據庫的建立

SWAT模型基礎數據庫的建立需要數據包括流域同一投影和空間坐標系統之下的DEM((Digital Elevation Model，90 m×90 m，來源于地理空間數據云)、土壤數據、坡向、坡度、土地利用數據、氣象水文數據。土壤數據與氣象數據由寒區旱區科學數據中心提供。根據研究區域土壤中砂土、黏土等已知量，利用SPAW(Soil-Plant-Air-Water)軟件計算得到土壤數據庫中所需參數和土壤類型分布圖(圖1)。氣象數據采用CMADS V1.1(The China Meteorological Assimilation Driving Datasets for the SWAT model V1.1)數據集[27-30]，本次研究區域空間范圍為40°46′—41°17′N，114°27′—115°30′E。土地利用數據來自地理空間數據云Landsat影像數據，共有4期，分別為1990年、2000年、2010年、2015年，分辨率均為30 m。在ArcGIS 10.2和ERDAS IMAGINE 9.3軟件平臺下，進行人工解譯，得到4期張家口清水河流域土地利用圖(圖2)。水文數據為張家口清水河流域水文站點2008—2016年的逐月徑流實測數據。

圖1 研究區域土壤類型

圖2 不同時期研究區域土地利用類型

3.2 子流域劃分及水文響應單元(HRU)的生成

HRU(Hydrologic Research Unit)是SWAT模型所特有的地方。在子流域的基礎上，根據土地利用數據、土壤類型和坡度，將流域劃分為具有獨特的土地利用和土壤組合的區域，使該模型能夠反映不同土地利用/覆蓋、土壤和其他水文條件的差異。根據DEM數據生成河網，并對研究流域總出口進行界定，以此為依據劃分子流域，并計算其參數。最后在集水面積閾值為42 km2的尺度上將研究區共劃分為23個子流域。為了保證模型的準確性，本文采用的面積閾值為土壤類型15%、土地利用10%、坡度類型10%[31]，生成99個水文響應單元。

3.3 土地利用變化情景設置

當模型精度達到要求后，通過兩種方法建立土地利用情景研究不同土地利用類型對徑流的影響，并將2015年土地利用設置為基礎情景(S0)[32-33]。第一種方法是極端土地利用法，旨在排除其他干擾因素，僅探討某一種土地利用類型對徑流的影響；第2種方法是基于實際情況，根據退耕還林還草政策建立土地利用情景。建立如下5種情景，面積見表1。

林地情景(S1)：保留建設用地及居民區和水域，將草地、耕地、未利用地全部設為林地；

草地情景(S2)：保留建設用地及居民區和水域，將林地、耕地、未利用地全部設為草地；

耕地情景(S3)：保留建設用地及居民區和水域，將草地、林地、未利用地全部設為耕地；

退耕還林情景(S4)：保留建設用地及居民區和水域、草地、林地、未利用地及25°以內的耕地，將流域內坡度大于25°的耕地設為林地。

退耕還草情景(S5)：保留建設用地及居民區和水域、草地、林地、未利用地及25°以內的耕地，將流域內坡度大于25°的耕地設為草地。

表1 不同土地利用情景的面積統計

4 結果與分析

4.1 流域土地利用變化分析

根據2017年發布的《土地利用現狀分類》，并結合流域土地利用實際情況，制定張家口清水河流域的土地利用分類表，將流域土地利用類型劃分為林地、草地、耕地、建設用地與居民區、未利用地、水域6大類[34]。對1990年、2000年、2010年和2015年崇禮地區張家口清水河流域土地利用分類數據進行統計計算(表2)。

由表2可知,林地、草地、耕地是崇禮地區張家口清水河流域的主要土地利用類型，這3種土地利用類型面積總和所占比例均較大，1990年為98.58%，2000年為98.58%，2010年為97.49%，2015年為96.14%；而建設用地及居民區、水域、未利用地所占比例均較小。

從不同年份不同土地利用類型的面積變化趨勢來看，林地、水域整體呈現下降趨勢，草地、耕地、建設用地及居民區、未利用地整體均呈現上升趨勢。就土地利用變化速率而言，在1990—2015年期間，各土地利用類型中未利用地年變化速率最快，這是因為25 a內土地受自然力(干旱、暴雨等)或人類不合理開發利用(不適當開墾、亂伐，不合理的種植制度和灌溉，不當使用化肥、農藥等)導致土地質量下降、生產力衰退，從而使未利用地面積不斷增大[35]。

表2 1990-2015年張家口清水河流域土地利用變化

4.2 模型校正、驗證

本文使用SWAT-CUP(SWAT-Calibration and Uncertainty Programs)軟件進行自動率定，然后結合率定結果進行手動調錯。最終篩選出對流域徑流模擬結果影響較大的參數包括CN2(徑流曲線)、SOL_AWC(土層有效含水量，mm/mm)、SOL_K(表層土壤飽和水力傳導率，mm/h)、CANMX(最大冠層截流量，mm)、ESCO(土壤蒸發補償因子)、Alpha_BF(基流α系數，d)、CH_K2(河道有效水力傳導度，mm/h)、GW-DELAY(地下水延遲時間，d)，其余參數影響較小[36-37]。

分析清水河流域水文站點2008—2016年徑流實測數據，將模型預熱期設為2008—2009年，校準期設定為2010—2013年，驗證期設定為2014—2016年，本研究中模型校準使用的手動調錯的方式，通過不斷調整的參數率定值，以期達到較好的模擬效果。得到的校正期和驗證期的模擬效果見圖3。

圖3 校正期與驗證期模擬與實測月徑流量數據對比

根據圖3，校準期和驗證期月徑流量模擬值與實測值變化趨勢相對一致，模擬效果較好。且NSE與決定性系數R2的值均高于0.6，說明SWAT模型在張家口清水河流域具有很好的適用性，可信度高。

4.3 土地利用情景下的徑流模擬

根據設定的5種土地利用情景，模擬張家口清水河流域徑流量結果(表3)。從表3可以看出，地表徑流的大小順序為S3>S0>S4>S5>S2>S1。模擬結果表明，與基礎情景(S0)相比，隨著地表草地和林地面積增加，流域年均徑流量減少，隨著耕地面積的增加，流域年均徑流量也隨著增加。這說明了耕地具有增加徑流的作用，林地和草地具有截留徑流的作用，且林地的截留作用更加明顯。

表3 不同土地利用情景下流域徑流變化

4.4 不同土地利用類型對流域產流的影響系數分析

根據前文可知，林地、草地、耕地是崇禮地區張家口清水河流域的主要土地利用類型，因此利用MATLAB R2017b進行計算，得到Ma，Mb，Mc的值即分別為林地、草地、耕地3種土地利用類型對流域徑流量的影響系數(表4)。

表4 土地利用變化對徑流量的影響系數

注：Δ1990—2000，Δ2000—2010，Δ2010—2015，Δ1990—2015分別表示1990—2000年、2000—2010年、2010—2015年、1990—2015年各土地利用類型對徑流總量的影響。

林地對徑流的影響系數為1.908 m3/(km2·h)，說明林地面積變化與徑流量變化呈負相關；研究期內，1990—2015年林地面積持續減少，被其他土地利用類型所占據，造成了徑流量的不斷增加。草地對徑流的影響系數為1.836 m3/(km2·h)，說明草地面積變化與徑流量變化呈負相關；研究期內，1990—2015年草地面積持續增加，造成了徑流量的不斷減少。耕地對徑流的影響系數為3.168 m3/(km2·h)，說明耕地具有促進徑流的作用，耕地面積變化與徑流量變化呈正相關；研究期內，1990—2015年耕地面積持續增加，造成了徑流的不斷增加。

5 結 論

(1) 建立了張家口清水河流域的月均徑流模擬模型，校準期和驗證期的決定性系數R2和NSE的值均高于0.6，說明SWAT模型在張家口清水河流域具有很好的適用性。

(2) 林地、草地、耕地是崇禮地區張家口清水河流域的主要土地利用類型，隨時間遷移，林地面積呈下降趨勢，草地、耕地均呈現上升趨勢，未利用地的年變化速率最快。

(3) 耕地對徑流的影響系數為正值，表明耕地具有增加徑流的作用；林地對徑流的影響系數為負值，表明林地和草地具有截留徑流的作用，且林地的截留作用更加明顯。

張家口清水河流域的水資源及水土保持工作，要以保護生態作為基礎，調整優化土地利用類型，科學合理布局，從而減少徑流，加大流域出水量。