SWAT模型在清江流域上游径流模拟中的应用

谭 珉，孙 毅*，胡 洋，宋 健，杨松元

(1.湖北民族大学 林学园艺学院，湖北 恩施 445000;2.湖北省地质局 第二地质大队，湖北 恩施 445000；3.恩施州水文资源勘测局，湖北 恩施 445000)

径流作为水循环中的基本环节和重要组成部分，同时也是水量平衡的基本要素，它的变化可反映区域水文过程与水资源量的变化[1].径流模拟作为揭示径流对环境变化、人为活动影响响应规律的研究基础，现已成为水文过程模拟研究的重要环节之一[2].

目前，随着人类活动对区域环境影响日益增强，下垫面环境变化显著，传统的概念性水文模型对径流形成的空间异质性缺少考虑，致使在进行水文模拟过程中存在不足.而SWAT(soil and water assessment tool)模型可以在各种土地利用、土壤类别以及管理模式下进行模拟，并通过强大的内部物理运行机制加以运算，自20世纪90年代开发以来便在不同国家与地区开展相关应用研究，应用研究内容包括水文模拟、水沙模拟、养分输移模拟、环境变化下的水文效应等方面[3-4].SWAT模型在清江流域的应用主要包括非点源污染模拟[5-6]、环境变化对径流的响应[7-8]、旱灾监测[9]和径流模拟[10]，这些研究多以大尺度为研究对象，对中小尺度研究较少.

清江上游是清江的发源地，其水量水质发生变化一方面会影响当地的水资源环境，另一方面也会对中下游地区水资源环境造成连锁反应.开展SWAT模型在清江上游流域的应用研究，不仅对了解流域径流特点具有现实意义，而且还可以进一步丰富SWAT模型在清江流域的研究内容.本文以恩施出口断面控制的清江上游流域为研究对象构建SWAT模型，使用SUFI-2算法分析模型的敏感性参数，并进行参数校正与验证，评估SWAT模型在清江上游的适应性，为清江上游水资源调配提供技术支持，同时也为后续研究打下基础.

1 研究区概况

清江上游流经湖北省恩施州境内，干流自西向东分别流经利川市和恩施市，恩施市之后为清江中游[11].河段全长153 km，流域面积2 928 km2,流域内岩溶地貌分布广泛，地下暗河、溶洞、伏流发育.流域四周高、中间低，总体上从西部向东部延伸.清江上游地处中亚热带季风气候区，年平均降水量1 466 mm，年平均流量464 m3/s，降雨集中在4-9月，汛期多出现在6-7月，1月降水量最小，7月降水量最大.清江上游流域支流主要包括长偏河、云龙河、龙桥河等.

流域内有9个雨量站、2个水文站，其中以流域出口的恩施水文站作为径流模拟站点.

2 研究方法

2.1 SWAT模型

SWAT模型是美国农业部20世纪90年代针对复杂大中尺度流域开发的分布式水文模型.该模型采用模块化结构进行模型组成与模型开发，由亚流域模块和汇流运算模块组成.亚流域模块包括8个部分：水文过程、气候、沉积物、土壤、植物生长、养分、农药和管理方式.汇流运算模块由2个部分组成：河道径流运算与储水体汇流运算[3,12].径流模拟主要在亚流域模块中的水文过程部分中进行，水文过程部分估算地表径流主要包括两种方法：SCS曲线数法与Green & Ampt下渗法，其计算公式如下[13].

1) SCS曲线数法

(1)

式中Qsurf表示累计径流量，mm；Rday表示某天的雨深，mm；Ia表示最初损失量，包括产流前的地表填洼量、植被截留量和下渗量，mm；S表示滞留参数，mm.滞留参数的定义为：S=25.4(1 000/CN-10)，其中CN表示某天的曲线数.

2) Green & Ampt下渗法

(2)

式中finf,t表示t时刻的下渗率，mm/h；Ke表示有效渗透系数，mm/h；ψwf表示潮湿锋处的基质势，mm；Δθv表示潮湿锋处土壤体积含水量的变化量，mm/mm；Finf,t表示t时刻的累计下渗量，mm.

2.2 SUFI-2算法

SUFI-2算法(sequential uncertainty fitting version 2)是内置于SWAT-CUP(calibration and uncertainty programs)软件中的一种参数估计最优化方法.它具有运行效率高、运行次数少、模拟精度高[14]等优点，在模型参数敏感性分析、率定以及水文模型不确定性应用过程中使用广泛.

该算法的原理为：首先确定一组较大的参数值范围，然后不断进行迭代运算，逐渐缩小参数值范围，最终使模拟值不断接近实测值[15].该算法通过不确定性P因子和R因子对模型的不确定性结果进行评估，其中P因子的取值范围为0到1，R因子的取值范围为0到+∞，最终以P因子接近1和R因子接近0的程度来判断模型校正效果.

SUFI-2算法步骤为：① 选择目标函数；② 确定初始参数范围和计算次数；③ 选择目标参数；④ 进行拉丁超立方抽样；⑤ 对抽样组合进行评价，确定模型敏感性参数；⑥ 参数敏感性分析.

2.3 模型模拟精度评价指标

采用国内外学者认可度较高的Nash效率系数(Nash-Sutcliffe efficiency,Ens)、相关系数(Correlation Coefficient，R2)、百分比偏差(Percent bias，PBIAS)3个评价指标进行模拟结果精度评价.

1) Nash效率系数(Ens) 该评价指标反映模拟值与实际测值两者的匹配程度，计算公式如下：

(3)

2) 相关系数(R2) 该评价指标反映模拟值与实际测值两者的吻合状况，计算公式如下：

(4)

表1 SWAT模型模拟结果评价Tab.1 Evaluation of SWAT model simulation results

3) 百分比偏差(Percent bias,PBIAS) 该指标反映模拟值与实际测值两者的偏差程度，计算公式如下：

(5)

式中Qobs为实际流量测值，m3/s，Qsim为模拟流量值，m3/s，n为实际测值个数.当PBIAS=0时，表示模拟值与实际测值一致；当PBIAS值为正值，说明模拟值偏小；当PBIAS值为负值，说明模拟值偏大.

参考前人对SWAT模型模拟结果的评价体系[16-17]，将Nash效率系数Ens、相关系数R2和百分比偏差PBIAS分为4个等级，分级情况如表1所示.

3 SWAT模型构建

3.1 基础数据库构建

输入到SWAT模型的主要原始数据包括研究区地形参数数据、土地利用类型数据、土壤种类数据、气象水文数据以及与流域相关的其他地理辅助数据.清江上游流域SWAT模型数据来源、格式及尺度如表2所示.

表2 研究区SWAT模型输入数据Tab.2 Input data of SWAT model in the study area

3.1.1 数字高程模型(DEM) 本文使用的数字高程数据来源于NASA EARTHDATA官网，数据尺度为12.5 m，数据格式为GRID.运用GIS软件对原始DEM数据进行图幅拼接、DEM数据裁剪、投影转换、填洼处理等一系列处理.

3.1.2 土地利用数据 土地利用数据来源于全球30 m精细地表覆盖产品[18]，分辨率为30 m，它是由我国研制的高精度全球地表覆盖数据，目前已发布GlobeLand30 2000、GlobeLand30 2010和GlobeLand30 2020三个版本，GlobeLand30数据共包括10个一级土地利用类型，分别是耕地、林地、草地、灌木、湿地、水体、苔原、人造地表、裸地、冰川和永久积雪.本文采用GlobeLand30 2000数据作为研究时段(1998-2005年)的模型输入数据.利用GIS软件对原始数据进行裁剪和投影，然后根据SWAT模型代码重新对土地利用类型进行分类，最终划分为5种土地利用类型.研究区各土地利用类型面积、占比和代码如表3所示.

表3 研究区土地利用类型分类表 Tab.3 Classification of land use types in the study area

3.1.3 土壤数据 土壤数据来源于联合国粮农组织发布的1∶100万世界土壤数据库(harmonized world soil database,HWSD)土壤类型数据集，数据格式为GRID，数据尺度为1 000 m的网格数据.因为HWSD数据与SWAT模型采用的美国农业部简化标准基本相同，所以不需转换土壤粒度，保证了土壤数据的准确性，降低建立SWAT模型土壤数据库的难度[19].本文输入SWAT模型的绝大部分土壤物理参数数据可从HWSD数据库中直接读取；其余参数如土壤水分参数中的土层结构、土壤湿密度、土层有效含水量和饱和水力传导系数使用土壤水特性软件SPAW(Soil Plant Air Water)中的Soil Characteristics模块计算；土壤侵蚀力因子利用Williams提出的修正通用土壤流失方程计算，对于难以获取计算的参数采用模型默认值.

运用GIS软件对HWSD数据进行裁剪、投影转换得到15种土壤类型，然后查询HWSD数据库属性表，将土壤物理属性相同的土壤类型进行归类，以面积最大的一类土壤作为该类土壤的代表性土壤进行重分类，减少模型计算量，最终重分类为8种土壤类型，重分类结果如表4所示.

表4 研究区土壤重分类结果Tab.4 Results of soil reclassification in the study area

3.1.4 气象、径流数据 输入SWAT模型的气象数据包含研究区的日降水量、日气温、日辐射、日相对湿度以及日风速，主要用于计算天气发生器各气象参数，其中降水、气温数据是模型的必要数据，可从气象站点获取.本文的气象资料主要分为两部分：一是研究区内利川、恩施两个气象站1998-2005年8年日气象资料，利用SWAT Weather程序计算各类气象资料，然后将其转换成模型需要的输入格式并构建用户气象数据库；二是研究区内气象站点少，为保证模拟精度，添加研究区内9个雨量站1998-2005年日降水数据.降水、气温数据格式根据模型需要的格式进行预处理.

按照SWAT-CUP软件径流观测数据输入格式对恩施水文站1998-2005年月径流量数据进行输入格式处理，为后期模型率定与验证提供参考.

3.2 子流域及水文响应单元(HRU)划分

加载前期处理好的DEM数据，参考前人研究成果[20-21],利用Burn in法添加流域实际河网进行河网提取与流域边界的确定，使用模型推荐的最优集水面积阈值并确定流域出水口、入水口和总出水口后，模型自动划分子流域，计算子流域地形参数.最终将清江上游流域划分为21个子流域.

在SWAT模型中添加之前处理好的土地利用及土壤数据并通过索引表链接到SWAT模型数据库,重新对土地利用及土壤数据进行分类.在重分类的土地利用数据与土壤数据的基础上，结合清江上游流域实际情况、模型运行效率以及参考相关研究成果[5,19-22]，分别定义土地利用类、土壤类和地形坡度类的阈值为5%、10%、10%，最终将清江上游流域划分为168个HRU.

3.3 参数敏感性分析

水文模型通过数学物理方程对现实水文过程进行描述与模拟，模型内部涉及到大量与流域空间属性或时间属性相关的参数，这些参数繁多且复杂[23].SWAT模型涵盖许多参数，对模型参数进行敏感性分析不仅可以了解研判各参数变化对模拟结果的影响深度，而且可以筛选出灵敏度高的模型参数，为后期进行针对性的参数校准奠定基础，减少参数校准过程中的不确定性和盲目性.本文使用SWAT-CUP软件中SUFI-2算法的全局敏感性分析进行模型参数的敏感性分析.

4 结果与分析

4.1 参数敏感性结果分析

根据模型运算原理及各参数含义，结合其他研究者的研究成果[16,24-25]，在SWAT-CUP软件中筛选与径流相关的参数进行全局敏感性分析，依照SUFI-2算法步骤筛选出研究区的模型敏感性参数，在拉丁超立方计算结果的基础上不断调整敏感性参数的范围，依次进行多次迭代运算，完成参数校准工作.最终敏感性参数结果如表5所示.

表5 参数敏感性排名及参数最佳取值结果Tab.5 Ranking of sensitivity of parameters and results of the best values of parameters

敏感性参数采用T检测法进行结果评估，t值(t-stat)代表参数灵敏性，该值绝对值越大，表明参数灵敏性越高；P值(P-value)代表参数的灵敏性显著水平，该值越接近于0，表明参数的灵敏性显著水平越高[26]，同时满足两者说明该参数越敏感.

由表5可知，选取的19个敏感性参数大致可以分为6种类型：地表径流参数(CN2、SURLAG)、地下径流参数(REVAPMN、GWQMN、GW_REVAP、ALPHA_BF、GW_DELAY和RCHRG_DP)、地形特征参数(HRU_SLP、SLSUBBSN)、下垫面特性参数(OV_N、EPCO、ESCO)、河道汇流参数(ALPHA_BNK、CH_K2、CH_N2)与土壤性质参数(SOL_BD(..)、SOL_AWC(..)、SOL_K(..)).

由表5可知，各敏感性参数对研究区的径流均产生不同差异的影响.一般认为，当t值(t-stat)的绝对值大于2时，表明该敏感性参数对流域径流灵敏性以及对径流影响的显著性水平较高[27]，由此可知对研究区径流模拟影响较为显著的5个参数分别是SCS径流曲线数(CN2)、河岸调蓄的基流α因子(ALPHA_BNK)、植物吸收补偿因子(EPCO)、主河道水力传导系数(CH_K2)和浅层含水层快速蓄水的阈值水深(REVAPMN).

SCS径流曲线数是影响SWAT模型径流模拟精度的重要参数，表征下垫面的综合状况，决定了地表径流量大小.河岸调蓄的基流α因子又称河岸调蓄退水常数，表征河岸调蓄的退水曲线，反映河岸调蓄水量补给亚流域主要干流或支流的状况.植物吸收补偿因子是植物蒸散发需水量、蒸散发量与土壤有效水量的函数，该值会影响植物从土壤中吸收的水分含量.主河道水力传导系数又称浅层地下水再蒸发系数，它反映支流汇入主河道时地表径流的有效渗透率，对研究区水资源平衡计算以及基流量拟合状况有重要影响.浅层含水层快速蓄水的阈值水深与地下水水位有关，当浅层地下水的水位超过阈值水深，浅层蓄水层会再次发生蒸发.

4.2 模型参数率定与验证

大多模型敏感性参数运算早期的开始值为零，这会对SWAT模型输出的模拟结果造成影响.所以有必要设定模型的预热期，以合理预估模型敏感性参数的开始变量[26].在敏感性参数筛选好的基础上，利用研究区恩施水文站1998-2005年实测月径流数据进行参数校准与验证，其中设定1997年为预热期，1998-2002年为率定期，2003-2005年为验证期.

利用SWAT-CUP软件推荐的参数范围进行迭代运算，依据迭代运算的结果调整下一轮参数范围，依次类推，直至确定满足模拟精度的最佳参数值.本文进行3次迭代运算，每次运算500次，最终的最佳参数取值如表5所示.将迭代运算后的最佳参数取值带入模型，进行参数校准、验证，SWAT模型在率定期(1998-2002年)与验证期(2003-2005年)的模拟值与观测值对比结果及散点图如图1和图2所示.根据式(3)～(5)计算模拟精度，具体结果如表6所示.

图1 率定期(1998-2002年)实测径流与模拟径流结果对比及散点图Fig.1 Comparison of measured runoff and simulated runoff and scatter plot in rate period (1998-2002)

图2 验证期(2003-2005年)实测径流与模拟径流结果对比及散点图Fig.2 Comparison of measured runoff and simulated runoff and scatter plot during the verification period (2003-2005)

表6 月径流模拟评价指标结果统计表Tab.6 Statistical table of runoff simulation evaluation index results

由表6可知，在月径流模拟过程中，率定期(1998-2002年)与验证期(2003-2005年)的相关系数R2和Nash效率系数Ens均大于等于0.90，表明率定期与验证期的观测值与模拟值两者之间的线性相关性强，拟合效果好；率定期(1998-2002年)与验证期(2003-2005年)的百分比偏差PBIAS分别为0.48%和10%，两者均小于25%百分比偏差临界值，一方面表明观测值与模拟值两者的偏差程度小，另一方面说明率定期与验证期的月均模拟值均比月均实测值小.由表6可知，率定期和验证期的月均模拟值分别为80.85、63.14 m3/s，比观测值分别小了0.40、7.09 m3/s.

对于验证期(2003-2005年)Nash效率系数Ens和相关系数R2较高，百分比偏差PBIAS较低，由图2散点图可看出模拟径流值与观测径流值线性相关性好，散点分布在拟合直线的两侧，但部分点离拟合直线距离较远，从而导致验证期(2003-2005年)径流模拟值精确度较高，但准确度较低，造成这种现象的原因可能是构建模型的过程中土地利用输入数据为率定期(1998-2002年)期间的数据.

总体而言，SWAT模型在清江上游流域月径流模拟过程中，模拟流量过程曲线与实测流量过程曲线较为吻合，对流量峰值均有所响应，模拟等级高.说明SWAT模型能够很好再现清江上游流域月径流的实际变化过程，表明SWAT模型在清江上游流域径流模拟过程中具有良好的适应性，满足SWAT模型的模拟应用要求.

5 结论

径流模拟研究是水文过程模拟研究的重要环节之一，而SWAT模型作为径流模拟研究的重要工具，自开发以来便在世界不同地区、不同流域进行广泛应用.本文以恩施断面控制的清江上游流域为研究对象，整合研究区流域下垫面数据及气象水文数据构建研究区SWAT模型，结合SUFI-2算法进行敏感性参数分析及参数校正与验证,结论如下.

由全局敏感性分析可知，对研究区径流模拟影响较大的5个参数分别是SCS径流曲线数(CN2)、河岸调蓄的基流α因子(ALPHA_BNK)、植物吸收补偿因子(EPCO)、主河道水力传导系数(CH_K2)和浅层含水层快速蓄水的阈值水深(REVAPMN).

结合SUFI-2算法对模型参数进行校正与验证，结果显示：率定期(1998-2002年)相关系数R2、Nash效率系数Ens和百分比偏差PBIAS分别为0.96、0.96和0.48%；验证期(2003-2005年)相关系数R2、Nash效率系数Ens和百分比偏差PBIAS分别为0.90、0.91和10%，模型模拟结果良好，说明SWAT模型适用于清江上游流域月径流模拟.