基于SWAT模型的张家口清水河流域土地利用情景变化对径流影响研究

2019-06-26刘亭亭谢建治

水土保持研究 2019年4期

王磊, 刘亭亭, 谢建治

(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年耕地面积持续增加，造成了径流的不断增加。