基于改进SWAT模型的土地利用变化对流域设计洪水的影响研究

魏 冲，董晓华,2，成 洁，时 磊，刘 冀，喻 丹，杨芝辰，谭雪松

(1.三峡大学水利与环境学院，湖北 宜昌 443002；2.水资源安全保障湖北省协同创新中心，武汉 430072；3.深圳市水务规划设计院有限公司，深圳 518001；4.宜昌晟泰水电实业有限责任公司，湖北 宜昌 444299)

0 前 言

随着流域内城市的飞速发展，该流域的土地利用会有所改变。人类在流域内进行诸如流域河道整治、水库塘坝以及跨流域调水等水利工程建设的各种活动，必然会对流域下垫面产生一定的影响[1]，使土地利用/覆被状况发生变化，从而改变流域蒸发截留、填洼下渗等产汇流条件[2-4]，导致一系列的水文生态效应。研究流域内土地利用变化对水文过程的影响对流域内城市的发展有重要指导意义。

国内外学者应用各种形式的方法对土地利用变化的响应进行研究，而应用最广的还是能考虑流域自然条件的分布式水文模型[5]，它可以体现出流域的各种物理条件对水文过程的影响，并可以结合遥感与GIS，灵活设置不同的土地利用情景，从事更精细的研究[6]。SWAT模型开发于20世纪90年代[7]，是近年来在土地利用变化的水文响应研究领域应用比较成熟的模型[8]。王强[9]等人在太湖上游西苕溪流域构建SWAT模型，应用GWR模型在空间上定量评估了土地利用变化对径流过程的影响；董晓华[10]等人在淮河上游构建SWAT模型，通过设定不同的土地利用模式来模拟未来土地利用情景下的水文响应；吴淼[11]等对潢川上游流域构建SWAT模型，模拟流域不同土地利用下最大几日设计暴雨洪水过程，研究结果表明：该流域土地利用变化对设计洪水过程的影响比较显著。但是洪水过程与径流过程相比，具有时间短、强度大和产流快等特性，需要更为精细的模型进行模拟。雷超桂[12]等人在奉化江皎口水库流域构建HEC-HMS水文模型，分析不同频率的暴雨洪水对下垫面变化的响应情况，指出暴雨洪水受土地利用变化的影响程度不一，对设计频率高的洪水影响更明显；毛慧慧[13]等人基于新安江-海河模型，分析了永定河上游土地利用变化对设计洪水的影响；冯平[14]等人对滦河流域下游构建流域水文模型，通过分析土地利用变化对洪水过程的影响，发现水土保持工程对洪水过程的峰值及总量都有一定削减作用。

城市化会影响该地区的产汇流过程，从而加大其受到洪涝灾害的风险[15]，了解下垫面变化对流域设计洪水的影响对城市防洪以及城市发展有重大意义。深圳市石岩水库建于1960年3月，同年设立石岩水库雨量站，距今已经有近60年的历史。近年来深圳市发展迅猛，城市化过程非常剧烈，水库修建之初使用的设计暴雨及设计洪水是否适用于现在，亟待研究。本文对石岩流域1988-2016年间四期的土地利用遥感图进行解译，构建改进后的SWAT模型，模拟该流域不同设计频率的洪水过程。以此对不同时期的土地利用情景下不同频率设计暴雨对应的设计洪水进行模拟，研究其变化情况，为城市防洪管理及水库管理提供理论支持。

1 研究区概况

石岩水库位于深圳市石岩镇境内，地理位置为东经113°54′，北纬22°41′，流域面积约44 km2(图1)，库容为3 120 万m3，水面面积约2.98 km2。库区位于亚热带常绿林海洋性季风气候区，气候特征为均温高，日照时间长，无霜期时间长，多雨，且雨量在年内不均。根据资料统计，该地年均气温22.4 ℃，最高温38.7 ℃，最低温0.2 ℃；年均日照时数2 120 h，年均蒸发量1 322 mm，流域年均降雨量为1 628.3 mm，平均每年4-9月占其年降雨量的85%左右。石岩水库有6条支流，石岩河为其干流。

图1 石岩流域概况图Fig.1 Overview diagram of Shiyan River basin

2 数 据

2.1 土地利用数据

本研究所使用的土地利用图均来自对研究区域遥感影像数据进行解译。在土地利用/覆被分类的相关研究中，通常使用监督分类法和非监督分类法。本研究选用贝叶斯监督分类[16]，一种应用极广的监督分类方法。本文使用的遥感影像资料源于地理空间数据云提供的Landsat4-5和Landsat8的卫星数字产品，为减小云量和植被季节性变化给解译带来的误差，综合考虑这两方面的因素，选取成像时间为1988年11月24日、2002年11月7日，2010年12月23日和2016年2月7日的影像数据进行解译，此四期遥感影像的云量少，且植被处在相近生长周期。数据经过卫星系统自动校正后，选择2、3、4三个波段组合进行影像合成，再进行监督分类，得到1988、2002、2010和2016年的土地利用图。结果见图2。

图2 石岩流域土地利用图Fig.2 Landuse map of Shiyan River basin

2.2 水文气象数据

本研究选用气候预报系统CFSR(Climate Forecast System Reanalysis)提供的研究区1979-2010年共32 a的降雨、温度、太阳辐射、相对湿度等气象数据，计算天气发生器所需的参数。DEM数据从地理空间数据云获取，土壤数据来源于世界和谐土壤数据库(Harmonized World Soil Database-HWSD)1∶100万土壤数据(联合国粮农组织(Food and Agriculture Organization-FAO)和维也纳国际应用系统研究所(International Institute for Applied System Analysis-IIASA)所构建，土地利用数据选用通过遥感解译所得到的1988、2002、2010和2016年土地利用数据。

基于石岩水库雨量站1963-2010年的历史实测暴雨资料，采用P-III型频率曲线法进行适线分析，得到设计频率P=0.05%、0.1%、1%、2%、5%下的1、6、24、72 h年最大设计暴雨，结果见表1；再用推理公式求相应频率的设计洪水过程，用于率定模型参数。

表1 设计暴雨计算结果 mm

3 SWAT模型的构建

SWAT(Soil and Water Assessment Tool)模型是由美国农业部农业研究中心的研究团队于1994年研发的分布式水文模型[17]。该模型以日尺度为模拟步长，主要适用于大尺度流域的长序列径流模拟。

流域水文过程主要分为两个阶段：一是水文循环陆相阶段；二是水文循环汇流阶段。SWAT模型通过不同的土地利用、土壤类型及坡度组合将流域离散成若干个水文响应单元(HRU)，以此反映不同作物与土壤的蒸散发差异。模型分别计算每个HRU内的产流量，然后将每个子流域内所有HRU的产流量叠加，得到各个子流域的总产流量。再进行子流域坡面汇流，计算其主河道的水量，最后通过特征河长法或马斯京根法来进行汇流阶段的运算。

由于洪水事件的历时相对较短，需要更加精细的模拟步长来捕捉流量的瞬时变化，故以日为模拟步长的SWAT模型无法很好的应用于洪水模拟。本文对SWAT模型进行修改，建立研究区小时尺度的SWAT模型。

3.1 SWAT模型的改进

模型有两种方法来估算地表径流：SCS曲线法以及Green&Ampt下渗法，都是以日为默认的时间步长进行演算，故只能进行日尺度长序列的水文模拟。SCS曲线法是20世纪50年代在美国乡村广泛使用的经验模型，而Green&Ampt下渗法是在假设土壤剖面均质、前期水分在剖面中均匀分布的前提下连续模拟下渗过程的一个物理模型，具有一定的物理意义。

本文将Green&Ampt下渗法[18]作为计算日以下时间步长地表径流的方法，对模型进行改进。其计算公式见式(1)。

(1)

式中：Finf,t是给定时间步长下的累积下渗量，mm；Ke是有效渗透系数，mm/h；ψwf是湿润峰处的基质势，mm；Δθv是湿润峰处土壤体积含水量的变化量，mm/mm。

当降雨强度小于下渗率时，时段内的降雨全部渗入土壤，地表径流量为0；当降雨强度大于下渗率时，按土壤的下渗能力下渗，其余部分形成地表径流。对各个时段下渗雨量积分，得到累积下渗量，并计算地表径流量。

3.1.1 地表径流滞后模块的修改

由Green&Ampt下渗法计算的时段地面净雨通过一个滞后方程[式(2)]，将有部分净雨被滞留在HRU中，在后一时段再汇入河道中。该方程通过滞后系数surlag和汇流时间tconc计算滞留的地面净雨。而SWAT2005中的滞后方程基于日时间步长，并不适用于日以下时间步长。

(2)

(3)

式中：Δt是模拟步长，h；其余符号意义同前。

如图3可知，在汇流时间一定的情况下，当滞后系数取值从1～10，随着模拟步长Δt的增大，tconc/Δt逐渐减小，而汇流因子逐渐增大，汇入主河道的地表径流量也逐渐增大，符合“随着时间间隔的增加，将会有更多的地表径流汇入主河道”这一规律。

图3 汇流因子曲线图Fig.3 Confluence factor graph

3.1.2 河网汇流模块的修改

SWAT模型主河道汇流一般采用两种方法：特征河长法和马斯京根法。为方便计算，本文选用特征河长法，该方法将河段复杂的调蓄作用简化成线性关系(式(4))。

Vout,2=SC(Vin+Vstored,1)

(4)

式中：Vout,2表示时段末的出流量，m3；Vin表示时段内的入流量，m3；Vstored,1表示时段初的蓄水量，m3；SC表示蓄水系数。

蓄水系数与模拟时间步长之间存在式(5)的关系，且在河道演算过程中存在的诸如渗漏损失、蒸发损失及河岸调蓄等过程的模拟，模型均以日时间步长进行模拟。修改模型代码中的模拟步长，使其适用于任意Δt步长的模拟。

(5)

式中：TT表示水流的运动时间，h；Δt表示模拟时间步长，h。

修改前的模型，在时间空间上各存在三个循环层，即时间上的Δt→日→年和空间上的HRU→子流域→流域。修改后，取消原代码中时间上的年循环部分，将其替换成连续的场次循环。每场洪水的初始条件，包括初始土壤含水量、初始深层以及浅层含水层储量、初始河道储量，均是从模型的日模结果中提取。空间上，仍以HRU为最小计算单位，分别计算各个HRU的产水量。其中，冠层截留、地表径流、坡面汇流和河道流量演算以Δt为时间步长模拟，但是融雪、基流、侧流、蒸发以及含水层储水是以日为时间步长计算，然后平均到每个时间步长。随后汇总至子流域循环层进行坡面汇流，最后进入流域循环层进行河道流量演算。计算完每个Δt时间层上的循环后，再以同样的方法进行日循环和年循环，直至模型模拟期结束。

3.2 SWAT模型参数敏感性分析及自动率定方法

SWAT模型中常见的对径流模拟有影响的参数有26个，只有部分参数对模型的不确定性具有较大程度的影响，需对参数进行不确定性分析[19]，选出敏感性较强的参数来进行率定，以降低模型的不确定性。SWAT模型采用结合了LH抽样法的稳定性和OAT算法的精确性的LH-OAT算法[20]进行敏感性分析，可准确获取对模型结果影响较大的重要输入参数，提高模型适用范围。模型选用SCE-UA算法[21]率定模型参数，SCE-UA算法通过搜寻所有选定参数的可行范围，尽可能获取全局最优参数组。

本研究使用2010年土地利用和五种不同频率的设计暴雨及设计洪水作为模型输入，分析模型所有径流参数的敏感性，并对敏感性较强的参数进行率定。将2010年9月期间发生的实测小时降雨过程输入模型，模拟得到的小时径流过程聚合成日径流，与实测日径流数据进行对比，验证模型精度。(实测日径流数据由水库水位及取用水记录，根据水量平衡反推得到。)

敏感性排名前十的参数及其自动率定结果见表2，参数率定结果见表3，参数验证结果见表4和图4。

表2 敏感性参数及最优值Tab.2 Sensitive parameters and the optimal value

表3 模型率定结果Tab.3 Model calibration result

表4 模型验证结果Tab.4 Model validation result

图4 模型验证洪水过程Fig.4 The flood process of model validation

用于验证的三场实测洪水过程，第一场的模拟结果和实测结果相差较大，模拟的洪水过程有两天存在较大洪量，而实测洪水过程则只有一天存在较大洪量。由于实测洪水过程由水库水位数据反推得到，可能在实际过程中存在相应的泄洪措施，使得实测的径流过程在有降雨的情况下没有洪峰，模拟前由模型日模结果计算的初始条件可能存在些许偏差，影响到第一场暴雨的模拟结果，使结果产生较大误差；而第二场和第三场实测洪水过程的模拟精度很高，均合乎《水文情报预报规范》(GB/T 22482-2008)规定的要求。综合考虑，该模型模拟结果可以接受，模型可用于实际应用中。

4 结果及分析

4.1 土地利用变化情况

利用ArcGIS，统计流域各个时期各类土地利用的面积和比例，结果见表5。

由表5可知，林草地始终是石岩流域的主要土地利用类型，1988-2010年均超过流域一半的面积，但总体呈递减趋势，而城镇用地则呈持续递增趋势。1988-2016年间，林草地减少39%，耕地减少80%，城镇用地增加3 380%，未利用地增加93%。其中，1988-2002年间的土地利用变化最为剧烈，林草地减少27%，耕地减少61%，而城镇用地增加1 200%，未利用地增加440%；2002-2010年，林草地和耕地分别持续减少约10%和88%，城镇用地增加约114%，而未利用地减少30%；2010-2016年则相对平稳，只有城镇用地和未利用地发生了较大程度的变化。

表5 流域各期土地利用状况Tab.5 Basin periods of landuse status

耕地和林草地的持续减少，以及城镇用地的不断增加，标志着该流域城市化进程正迅猛发展，下垫面条件也随之发生剧烈变化。下垫面的改变直接导致洪水特性发生改变，适用于以前的设计洪水在下垫面发生改变后也会产生相应的变化。

4.2 不同土地利用的水文响应模拟结果

将五场不同频率的设计暴雨分别输入到不同土地利用条件的SWAT模型中,分别模拟四种不同时期土地利用下设计暴雨产生的洪水过程,并将之与对应的设计洪水过程进行对比，分析比较不同时期土地利用变化对不同频率的暴雨洪水的影响规律。模拟结果如图5。

图5 不同土地利用下设计洪水模拟结果Fig.5 The design flood simulation result under different landuse

由图5可知，同一频率的暴雨在四种不同时期土地利用情景下产生的洪水过程有明显差异。随时间的推进，同频率暴雨下产生的洪水总量与洪峰流量随之增加，且洪水起涨时间及峰现时间均呈逐渐提前的趋势；2010年与2016年情景下洪水过程线趋势相似，洪峰及洪量基本都接近设计洪水过程，洪水起涨阶段历时短且峰陡，退水阶段历时长且相对平缓；随着设计暴雨重现期的增大，2010年与2016年情景下洪水模拟过程逐渐趋近于实际设计洪水过程。

通过对不同年份的土地利用下暴雨洪水进行模拟，分析其结果，见表6和表7。

表6 不同土地利用情景下暴雨洪水模拟洪峰及变化量 万m3

表7 不同土地利用情景下暴雨洪水模拟洪量及变化量 万m3

根据表6和表7可得，同等频率的暴雨在1988-2016年四期土地利用情境下产生的洪水峰值呈逐次增大趋势。2002年模拟结果洪峰相对1988年模拟结果变化量与2016年模拟结果平均值相对于2010年模拟结果相近，分别为45.6和36 m3/s；而2010年模拟结果的洪峰相对于2002年则有巨大改变，其峰值平均相差205.6 m3/s。洪水重现期越短，土地利用变化下其洪峰和洪量的变化量越大。洪量与洪峰具有相同趋势，同等频率下暴雨在1988-2016年四期土地利用情境下产生的洪水洪量也呈逐步上升趋势。其中2002年模拟结果平均值相对1988年来说洪量仅增加了23.8 万m3，而2002-2010年则增加了634.2 万m3，2016年相对2010年也增加了84 万m3。

结合土地利用变化分析，1988-2002年与2002-2010年，林草地和耕地均大量减少，城镇用地均增加6 km2以上，但前者新增了6.6 km2的未利用地，后者则减少了2.4 km2；相比之下，后者的径流变化量远大于前者，其洪峰平均变化量为205.6 m3/s，洪量平均变化量为634.2 万m3；以此分析除却林草地和耕地有增大流域暴雨下渗量外，未利用地也对增大流域下渗量有积极作用。而2010-2016年间，林草地减少的速度变缓，耕地面积新增0.9 km2，城镇面积增加3.5 km2，同时未利用地减少2.8 km2，其间洪峰平均变化量为36 m3/s，洪量平均变化量为84 万m3，与分析结果相同。

5 结 语

本文分析了深圳市石岩流域1988-2016年间的土地利用变化情况，构建改进后的SWAT模型，模拟不同时期土地利用下设计暴雨的洪水过程，并分析其变化情况。主要结论如下：

(1)改进后的SWAT模型在石岩流域的洪水模拟中具有良好的适用性。其模拟精度均达到《水文情报预报规范》(GB/T 22482-2008)规定的要求，可用于该流域的洪水模拟研究。

(2)1988-2016年间，石岩流域下垫面发生了比较剧烈的变化。流域内林草地和耕地大量减少，城镇用地大量增加，且变化主要集中发生在1988-2010年间。

(3)随着石岩流域不同时期下垫面的变化，流域设计洪水过程也发生了剧烈的变化。流域内大量的林草地和耕地向城镇用地发展，改变了流域下垫面的渗透系数，进而影响其产汇流过程。在下垫面变化比较剧烈的1988-2010年间，各个频率的设计暴雨所产生的设计洪水洪峰均增加了一倍左右，洪量平均增加了77%，且设计频率越低的设计洪水所受到的影响大于设计频率较高的设计洪水，增加了其洪水风险。