基于SWAT模型的鄂西犟河流域非点源污染模拟

陈俊宏， 张利华*， 马永明

(1.中国地质大学(武汉) 地理与信息工程学院， 武汉 430074； 2.昭通学院地理科学与旅游学院， 云南 昭通 657000)

对于流域内水质污染的研究，采用定量化的水文模型计算和评估，已成为有效量化的手段[1]，各种分布式水文模型如WEPPO模型、AGNPS模型、EPIC模型和SWAT模型等大量出现，在国内外的研究中应用十分广泛[2]，但各类模型通常都需要输入大量的参数，但目前的水文站点分布尚有不足，如我国就存在水文站点多分布在经济发达的中东部地区的状况[3]，许多山区地区缺少资料或者资料不足，各类模型难以应用，治理受到制约，因此无资料地区的水文研究成为水文预报研究中的重点和难点.

目前，无资料地区的研究多应用区域化方法，如Narbondo等[4]以乌拉圭国内的13个流域为例研究了区域化方法，开发了基于相似性概念的流域分区方法，建立了降雨-径流模型参数与流域物理属性的函数关系，以获取无资料流域的参数，从而进行模拟.而在区域化方法中最易使用的是基于水文相似性的参数移植法[5]，Dumedah等[6]在加拿大安大略省南部的8个流域应用数据同化进行了径流模拟，根据得到的结果，他们认为模型参数值具有高度的共同性，可以将另一个流域的参数应用到给定的流域，丁洋等[7]对水文资料缺失的湖北省松滋市小南海流域应用了参数移植法进行研究，将具有相似性的洞庭湖等流域的水文参数移植到小南海流域，得到了流域的径流模拟结果，Cho等[8]将美国乔治亚州西南部沿海平原流域的经过了率定优化的参数移植到同一地理区域的其他五个流域，不再进行参数校准，得到了较好的模拟结果，樊琨[9]对甘肃省泾河上游区、三水河等流域应用了参数移植法，模拟结果良好.

SWAT模型可应用范围广，在国内外的各个流域均表现出适用性，可用于非点源污染分析[10].有研究应用模型参数移植，实现了缺资料地区的SWAT模型构建[11-13]，SWAT模型在犟河流域的上级流域已经表现出了适用性[14-16]，陈昊荣[17]在汉江流域成功应用了参数移植法构建了SWAT模型.犟河流域位于南水北调工程核心区，在调水后曾表现出水质污染严重状态[18]，亟待解决水环境治理问题，但是由于水文资料缺乏，流域内的相关研究较少，因此基于SWAT模型，应用参数移植法进行犟河流域的SWAT模型构建和非点源污染模拟，为该流域污染防治提供研究依据.

1 材料和方法

1.1 研究区概况

犟河是长江流域汉江右岸的二级支流，同时也是堵河的一级支流，其起于十堰市茅箭区大川镇，穿越南水北调核心十堰市，在黄龙镇东湾村与堵河交汇，再汇入汉江至丹江口水库，作为十堰市的纳污河，长期接受十堰市城市污染排放.犟河干流全长50.2 km，流域范围为32°29′21″N～32°42′37″N，110°31′38″E～110°43′26″E，流域面积316 km2，流域地势上南高北低，西南高、东北低，整体海拔139～1 441 m，为亚热带季风气候，雨热同期，季节变化明显，生态基流小，在秋冬季枯水期流量较小.流域地理位置见图1，其位于上级流域堵河流域的东北角.

图1 研究区概况图 Fig.1 Geographical location

1.2 模型简介

SWAT模型是由1 013 个中间变量、701 个方程组成的综合的、分布式的模型体系[19]，经过了多年的开发和应用，已经有了比较成熟的集成软件体系如ArcSWAT、SWAT-CUP方便研究人员进行分析.SWAT模型构建简单，输入高程数据、土壤数据、土地利用数据、气象数据等基础地理数据后，即可初步构建流域的SWAT模型，根据研究方向，再输入其他管理数据，并对模型进行率定和验证，便得到模拟结果.

1.3 数据处理

从中国气象同化驱动数据集CMADS[20]中裁剪获得了研究区域的十堰站气象数据；从地理空间数据云(http：//www.gscloud.cn)，下载到了包含犟河流域的DEM数据，空间分辨率为30 m，裁剪得到的DEM数据见图2(a)；土地利用数据主要通过目视遥感影像解译获取，下载获取了包含研究区域的Landsat 8遥感影像，并结合GoogleEarth在线地图，在ENVI 5.3软件中进行遥感解译和裁剪等，得到了研究区域的土地利用类型数据(图2(b))；土壤数据参考姜晓峰等[21]论述的SWAT模型土壤数据库构建方法，基于HWSD数据库，建立研究区的土壤数据库(图2(c))，犟河流域土壤主要有四种类型：石灰性冲积土(FLc)、两种发育程度不同的弱发育淋溶土(LVh1、LVh2)、不饱和始成土(CMd).实测径流数据及实测水质数据收集自十堰市水文资源局，径流数据包括了竹山站和黄龙滩站两个水文站2014—2017年的数据，水质数据主要包括氨氮和总磷.其他管理数据如农业数据收集自十堰市农业农村局，包括张湾区各乡镇的农业化肥施用量数据；禽畜、水产养殖数据、人口生活污水、工业排放数据收集自十堰市统计年鉴及十堰市生态环境局.

此外，输入SWAT模型土壤数据库的部分参数使用了软件SPAW(Soil Plant Atmosphere Water)进行辅助计算，如土壤湿密度、土壤层有效持水量、饱和水利传导系数.数据库中USLE土壤可蚀性因子K，利用土壤粒径分布进行计算得到：

Kusle=fcsand×fcl-si×forgc×fhisand，

(1)

其中

fcsand=

fhisand=

式中，ms代表在美国制土壤粒径分类下，土壤中粒径为0.05～2.00 mm的砂砾的百分含量，msilt代表粒径为0.002～0.05 mm的粉砂的百分含量，mc代表粒径小于0.002 mm的粘粒的百分含量，Corg代表土壤中有机碳的含量.

图2 基础地理数据Fig.2 Basic geographic data

1.4 模型的率定和验证

分布式水文模型往往有多种需要结合实测资料进行率定的参数[22-23]，在缺乏资料的地区，为了提高模型模拟的精度，目前国内外研究者一般采用区域化方法[24-25]解决参数问题，参数移植法[26]是其中的一种.朱阿兴等[27]提出地理环境越相似，地理特征越相近.参数移植应用了相同的思想，Parajka等[28]对奥地利的320个流域进行了参数移植的试验，他们的结果显示，当流域的自然属性如地貌、土壤、土质等相似时，除了高寒地区和平坦流域，区域化处理之后的模拟的结果良好.Kokkonen等[29]研究参数移植时得到结论，如果流域水文特征相似，那么可以将有实测数据的流域的全部率定参数都移植到相似的流域.

犟河流域是堵河流域的子流域，两者位于相似气候区，有相似的降水、地貌等特征，因此，应用参数移植法，将堵河流域参数率定结果用于犟河流域进行犟河流域的模型构建，利用SWAT-CUP进行参数校准以降低模型不确定性[30].

基于实测径流数据，先进行参数敏感性分析，再对敏感性参数进行率定.敏感性分析采用GLUE算法，设定模拟次数5 000 次；率定采用SUFI-2算法，单次率定设定模拟次数500 次，迭代，以线性回归决定系数R2和纳什效率系数NSE最大化作为目标函数，参数选择参考牛利强[16]及戴枫勇[31].堵河流域的径流模拟以2014—2016年为率定期，以2017年为验证期，犟河流域的非点源污染模拟以2014—2015年为率定期，以2016—2017年对氨氮模拟结果进行验证，以2018年对总磷模拟结果进行验证.

2 结果和分析

2.1 径流模拟

SWAT模型中对模拟结果的评价指标，采用纳什效率系数NSE及线性回归决定系数R2[12， 32-33]，一般认为在进行径流模拟时，R2>0.6和NSE>0.5时，模拟的结果可以接受[12].堵河流域有黄龙滩、竹山两个水文站点，基于SUFI-2算法进行迭代模拟，每次设定模拟次数为500 次，迭代率定计算了12 次，R2和NSE不再增大，且满足精度要求，率定结束，最终得到的径流模拟率定参数值见表1.

表1 堵河流域径流模拟参数取值范围及最佳值

将经过率定的参数返回输入模型，进行2017年的径流模拟.以2017年的竹山站数据作为验证，计算得到的R2=0.575，径流模拟与实测值相关性见图3，径流模拟结果见图4.

图3 2017年径流模拟与实测相关性Fig.3 Correlation between runoff simulation and measurement in 2017

图4 月径流量模拟结果Fig.4 Monthly runoff simulation results

2.2 非点源污染模拟

堵河流域的径流模拟完成，基于参数移植法，将堵河流域的径流模拟参数输入到犟河流域的模型之中，并进行污染物的参数率定，率定方法与径流模拟相似，只在参数选取上有所不同.以犟河流域东湾桥监测断面的氨氮和总磷水质数据作为实测值，由于输入模型的的污染负荷为总量，单位为kg，需要将浓度值转换为总量，计算公式为：

L=F×c×10-3×t，

(2)

式中，L为污染负荷总量，单位为kg，F为断面径流量，单位为m3·s-1，t为时间，单位为秒，计算得到的数据，在率定时，输入SWAT-CUP，以2016—2017年数据进行氨氮的模拟验证，以2018年数据进行总磷的模拟验证，得到的率定参数值见表2，模拟结果精度见表3，氨氮模拟结果见图5，总磷模拟结果见图6.

表2 犟河流域污染物模拟参数取值范围及最佳值

表3 污染物模拟结果精度

图5 氨氮模拟Fig.5 Ammonia nitrogen simulate

图6 总磷模拟 Fig.6 Total phosphorus simulate

氨氮模拟，验证期的R2达到0.76，满足精度要求，效果较好；而总磷模拟，验证期的R2为0.43，效果不够好，但可以看到其年际变化与氨氮模拟结果相似，到2018年，污染量明显减少.氨氮在2017年汛期，9、10月出现了高值，可能是由于降水量增大，土壤侵蚀、对地表的冲刷等，更多的污染物汇入到了犟河之中.从总体上来看，犟河流域在2014年、2015年的污染物排放量较大，水质较差，在十堰市政府的逐渐重视和投资治理下，2018年后，水体中污染物含量已经有了明显降低.

图7 犟河流域污染负荷与径流量的月周期变化 Fig.7 Monthly variation of pollution load and runoff in the Jiang River Basin

2.3 氮磷污染时间分布特征分析

基于模拟结果进行氨氮和总磷的时间分布特征分析，从图7可以看出，犟河流域的氨氮量和总磷量与径流量随时间的变化趋势相似，在汛期，径流量大，氨氮和总磷污染量也相应较大，至非汛期，径流量和非点源污染负荷变小，其中氨氮的污染负荷最大值出现在2014年10月，总磷最大值出现在2017年9月，均对应于汛期.2017年至2018年全年氨氮污染总量下降了70.63%，下降幅度最大(图8)，2014年至2015年下降幅度为40.15%，2015年至2016年下降幅度为43.58%，从下降幅度来看，2018年治理效果良好，从氨氮减少的量来看，2014年至2015年，氨氮排放减少了45 738.27 kg，减少的总量最大.2016年至2017年，氨氮污染总量增加了24.73%，根据查到的《汉江流域气候公报(降水)》数据，2017年，十堰市汛期的降水量较常年降水量偏多60%，其中8月28日—10月18日，全市降水量达到了527.8 mm，达到历史同期的3.19倍，可能是降水造成了径流量增大，更多的污染物通过径流，汇入到了犟河之中.总体来看，经过多年时间，污染量每年以约40%的速率下降，犟河的污染得到了有效治理.

图8 犟河流域氨氮负荷年际变化Fig.8 Interannual variation of ammonia nitrogen load in Jiang River Basin

2.4 氮磷污染空间分布特征分析

利用SWAT模型模拟输出的氨氮和总磷污染量，以2014年和2018年的数据作为对比，计算到每个子流域的单位面积上的污染负荷量.计算结果见图9.对比结果显示，无论是氨氮还是总磷的单位面积负荷，高值都主要集中在1、4、5、7、9、10、13、14、15、17、18、25、29等子流域，根据谷歌地球在线地图及遥感土地利用分类结果，这些子流域正好对应于城市建筑群，是居民居住的主要区域.此外，结合犟河流域的河网分布(图1)，在财神沟和大西沟，自上游至下游，污染逐渐加重.犟河干流南岸污染高于北岸.同时，将2014年与2018年的氨氮及总磷单位面积负荷量在各子流域进行对比(图10)，氨氮的单位面积负荷在高值区有明显的降低，而总磷值略有降低.两者的模拟结果在空间分布上明显相近.

图9 犟河流域污染单位面积负荷的空间分布 Fig.9 Spatial distribution of pollution load per unit area in the stubborn River Basin

图10 犟河流域各子流域污染物单位面积负荷变化Fig.10 Variation of pollutant unit area load in each sub basin of the Jiang River Basin

3 结果与讨论

以缺少径流资料的犟河流域为例，基于SWAT模型，运用参数移植法探索了径流参数的移植.

1) 对犟河流域的上级流域——堵河流域进行了径流模拟，模拟期精度达到要求，验证期竹山站结果R2=0.62，并将堵河流域的径流模拟参数移植到犟河流域.

2) 对犟河流域的非点源污染进行了模拟，验证期非点源污染模拟结果(氨氮R2=0.76，总磷R2=0.43)，总磷的模拟结果精度偏低，但其和氨氮的模拟结果时空分布表现出相似的特征，表明SWAT模型在犟河流域是适用的.

3) 犟河流域经过多年治理，污染负荷量在时空两层次上均表现出下降的趋势.2017年至2018年治理效果表现最佳，全年氨氮污染总量下降幅度超过70%，下一步可以进一步关注模拟出现高值的4、5、7、9等子流域，做更精确的污染防治.

目前的参数移植，相似流域的选取主要是基于水文相似性，如汪银龙等[1]研究选择了在空间上相邻板涧河与毫清河流域，陈磊等[23]研究大宁河流域时，选择了流域内部的子流域进行对比，这是两种根据空间邻近性寻找水文相似流域的方法.犟河流域的模拟表明，无资料流域的水文模拟参数从空间上考虑，可以参考上级流域的参数，对其他缺乏资料区域的研究具有参考价值.