王 磊, 劉亭亭, 謝建治

(1.北京市農(nóng)林科學(xué)院 植物營養(yǎng)與資源研究所, 北京 100097; 2.河北農(nóng)業(yè)大學(xué) 資源與環(huán)境學(xué)院, 河北 保定 071000)

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

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

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

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

1 研究區(qū)概況

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

2 研究方法

2.1 SWAT模型模擬法

SWAT模型是在30 a前由美國農(nóng)業(yè)部(USDA)和農(nóng)業(yè)研究中心(ARS)針對(duì)較大尺度流域開發(fā)的具有很強(qiáng)物理機(jī)制的長時(shí)間的分布式水文模型[14]。SWAT模型綜合考慮了自然因素和社會(huì)因素，能夠較好地模擬地表徑流的過程[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]。主要運(yùn)用在基礎(chǔ)數(shù)據(jù)缺乏的區(qū)域，通過模擬預(yù)測流域內(nèi)水、泥沙及農(nóng)業(yè)管理措施所產(chǎn)生的響應(yīng)分析，用于幫助進(jìn)行水資源管理[17]。

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

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

(1)

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

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

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

(2)

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

本文選取決定性系數(shù)R2和納什效率系數(shù)NSE來評(píng)價(jià)模型的模擬精度。

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

(3)

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

(4)

2.2 單一土地利用類型動(dòng)態(tài)函數(shù)

本文采用單一土地利用動(dòng)態(tài)度來表示土地利用類型變化速率。單一土地利用類型動(dòng)態(tài)度是指在一定時(shí)間范圍內(nèi)研究區(qū)內(nèi)一種土地利用類型速度變化，計(jì)算公式為：

(5)

式中：K表示研究時(shí)段內(nèi)某一土地利用類型動(dòng)態(tài)度；Ua表示研究期初某一種土地利用類型的數(shù)量；Ub表示研究期末某一種土地利用類型的數(shù)量；T表示研究時(shí)段長，當(dāng)T的時(shí)段設(shè)定為年時(shí)，K值為該研究區(qū)某種土地利用類型年變化率。

2.3 多目標(biāo)決策函數(shù)法

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

年平均徑流量是指時(shí)間步長內(nèi)不同土地利用類型對(duì)于徑流量影響之和，不同土地利用類型對(duì)于徑流量的影響將其量化，可正可負(fù)，但意義不同，當(dāng)為正時(shí)，代表土地利用類型與徑流量變化呈正相關(guān)，為負(fù)時(shí)代表土地利用類型與徑流量變化方向呈負(fù)相關(guān)。

年平均徑流量Xn的計(jì)算方程如下：

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

(6)

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

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

(7)

3 模型應(yīng)用

3.1 流域基礎(chǔ)數(shù)據(jù)庫的建立

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

圖1 研究區(qū)域土壤類型

圖2 不同時(shí)期研究區(qū)域土地利用類型

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

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

3.3 土地利用變化情景設(shè)置

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

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

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

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

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

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

表1 不同土地利用情景的面積統(tǒng)計(jì)

4 結(jié)果與分析

4.1 流域土地利用變化分析

根據(jù)2017年發(fā)布的《土地利用現(xiàn)狀分類》，并結(jié)合流域土地利用實(shí)際情況，制定張家口清水河流域的土地利用分類表，將流域土地利用類型劃分為林地、草地、耕地、建設(shè)用地與居民區(qū)、未利用地、水域6大類[34]。對(duì)1990年、2000年、2010年和2015年崇禮地區(qū)張家口清水河流域土地利用分類數(shù)據(jù)進(jìn)行統(tǒng)計(jì)計(jì)算(表2)。

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

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

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

4.2 模型校正、驗(yàn)證

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

分析清水河流域水文站點(diǎn)2008—2016年徑流實(shí)測數(shù)據(jù)，將模型預(yù)熱期設(shè)為2008—2009年，校準(zhǔn)期設(shè)定為2010—2013年，驗(yàn)證期設(shè)定為2014—2016年，本研究中模型校準(zhǔn)使用的手動(dòng)調(diào)錯(cuò)的方式，通過不斷調(diào)整的參數(shù)率定值，以期達(dá)到較好的模擬效果。得到的校正期和驗(yàn)證期的模擬效果見圖3。

圖3 校正期與驗(yàn)證期模擬與實(shí)測月徑流量數(shù)據(jù)對(duì)比

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

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

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

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

4.4 不同土地利用類型對(duì)流域產(chǎn)流的影響系數(shù)分析

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

表4 土地利用變化對(duì)徑流量的影響系數(shù)

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

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

5 結(jié) 論

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

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

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

張家口清水河流域的水資源及水土保持工作，要以保護(hù)生態(tài)作為基礎(chǔ)，調(diào)整優(yōu)化土地利用類型，科學(xué)合理布局，從而減少徑流，加大流域出水量。