基于SWAT 模型的沙塘川流域非点源氮磷污染特征及关键源区识别

张倩，巢世军，杨晓丽，李发娟，胡健

青海省环境科学研究设计院有限公司，西宁 810000

随着城市工业废水和生活污水等点源污染控制水平的提高，非点源污染的严重性日益凸显（吴蓓等，2007；马振邦等，2011；聂铁锋，2012）。非点源污染是指在工农业生产与人类生活中产生的污染物，如生活垃圾、农田中的化肥和农药等，在降雨和径流的冲刷下，从非特定的地点进入受纳水体而引起的污染（Olness，1995）。非点源污染具有形成复杂性、随机性、模糊性、污染负荷时空差异性显著等特点（李怀恩和李家科，2013；Xiang et al，2017），导致非点源污染无法进行集中处理和控制。有研究表明：在一些流域内，非点源污染已经超过点源污染，成为水环境质量下降的主要原因（李春林等，2013；罗倩等，2014）。因此，开展重点流域非点源氮磷污染特征分析对水环境质量管理具有重要意义（耿润哲等，2016；王萌等，2018）。

SWAT 模型作为一种发展成熟的分布式流域水文模型，因其强大的功能、先进的模型结构和高效的计算性能（Hosseini et al，2011；王晓朋等，2015；周帅等，2019），已经被国内外学者广泛应用于流域水文过程模拟，气候变化与土地利用对径流的响应以及非点源污染负荷估算及形成机制探讨，环境变化及农业管理措施对水文水质的影响等研究中（李明涛，2014；Yasarer et al，2016；徐畅等，2019），其在国内外各个流域均表现出适用性（乔卫芳等，2013；刘洪延，2019；袁远等，2020）。

湟水流域位于青藏高原与黄土高原的生态过渡地带，承载着青海省62%的人口和64%的GDP，是全省人口最密集、开发程度最高、经济最发达的地区，也是国家战略格局中兰西国家级重点开发区（兰州—西宁城市群）的核心组成部分，在区域经济社会发展和生态环境保护中具有重要地位（李云成等，2017）。近年来，通过对湟水流域涉水企业和生活污水的大力治理，流域整体水质有了显著的改善和提升，但随着时间推移改善效果越来越弱，说明非点源污染对水体污染的贡献率在逐渐增大。青海省的非点源污染缺乏科学性，治理实践具有一定的盲目性，主要原因是对非点源污染关键源区和污染负荷空间分布不清楚。沙塘川位于青海省东北部海东市互助土族自治县（以下简称互助县）西部，属黄河流域，为黄河干流的二级支流、湟水河一级支流，是青海省东部重要的农业产区，流域内工业点源较少，主要为互助县第一污水处理厂及区域内白酒制造等企业，面源以种植业、农村生活、养殖业和天然面源为主，是一个研究非点源污染的理想流域。因此本文以沙塘川流域为研究对象，分析了流域内非点源氮磷污染特征，并且对关键源区进行了识别，旨在探索一套适用于高原高海拔地区非点源污染关键源区识别的技术方法。

1 材料与方法

1.1 研究区概况

沙塘川流域（图1）主体位于青海省东北部海东市互助县境内西部，流域北侧为祁连山支脉达板山，南侧小部分流域位于西宁市城东区韵家口镇，地理坐标为北纬36°15′—36°57′和东经101°79′—102°21′，流域面积1113.51 km2。流域地势北高南低，海拔在2184—4333 m，北部为高山区，海拔4000 m 以上的达坂山横贯流域北部，中部为脑山地区，南部为浅山丘陵区，中下游河段河谷阶地为川水谷地。气候类型为高原半干旱大陆性气候，其特点为：降雨量少而集中，蒸发量大，日温差量大，冻土期长，无霜期短，紫外线强，冬长夏短，春秋相连。流域内总人口26.03 万人，其中西宁市人口1.34 万人，海东市互助县人口24.69 万人。农作物以粮食和油料作物为主，总播种面积37034 hm2。2016年末大牲畜存栏数34301 头，其中牛29900 头，占大牲畜存栏数的87.17%，其次为骡、马、驴；羊、猪、鸡存栏数分别为117256 只、116481 头、193676只*互助土族自治县统计局.2016.青海省互助土族自治县国民经济统计资料.内部资料.西宁市城东区统计局.2016.西宁市城东区统计年鉴.内部资料.。

图1 研究区地理位置示意图Fig.1 Schematic diagram of the geographical location of the study area

1.2 研究方法及数据来源

1.2.1 SWAT 模型数据库的建立

本文所使用的SWAT 模型基于ArcGIS 10.2 平台，充分利用GIS 强大的空间分析能力和地理数据的管理能力，对模型运行所需的地理参数进行统一的采集、提取和存储，通过建立查询表，将地理空间数据和非空间属性数据链接起来，本文所涉及地理空间数据及属性数据如表1、表2 所示。该模型自开发以来，在世界各国得到了广泛的应用，并取得了较好的模拟效果（吕乐婷等，2014；Abbaspourk et al，2015；Meaurio et al，2015）。模型的输入数据，包括气象、土壤属性数据、土地利用及作物管理数据和点源排放数据、非点源输入或产生数据。

表1 流域地理空间数据库Tab.1 Geospatial database for the watershed

表2 流域属性数据库Tab.2 Attribute database for the watershed

1.2.2 SWAT 模型的建立

（1）河网的生成以及子流域的划分

基于SWAT 模型的流域划分模块（Watershed Delineation）来进行河网的提取及子流域的划分，将沙塘川流域划分为19 个子流域（图2）。

图2 沙塘川流域SWAT 子流域划分结果Fig.2 Results of SWAT sub-watershed delineation in Shatangchuan watershed

（2）土地利用、土壤的加载和水文响应单元的划分

在划分水文响应单元之前，首先利用Arc SWAT 的土地利用、土壤和坡度定义模块（Land Use/Soils/Slope Definition）加载土地利用与土壤图，并选择对应的索引关系表来进行空间叠加分析。其次，采用多水文响应单元划分方法，设置土地利用阈值、土壤面积阈值和坡度等级，本次模拟的土地利用阈值均设为5%，土壤面积阈值设为10%，坡度等级阈值设置为10%。

（3）模型运行时期设定

根据收集到的流域相关资料，选取2008—2016年进行模拟，其中，2008—2010年为预热期，2011—2014年为率定期，2015—2016年为验证期。

1.2.3 率定与验证

模型适用性评价参数采用决定系数（R2）和纳什效率系数（NS）来衡量。其计算公式分别为：

式中：Qm, i为实测径流量；Qp, i为模拟流量；Qm,avg为多年实测平均流量；Qp,avg为多年模拟平均流量；n为实测时间年份。

从R2可以得出实测值和模拟值之间的相关程度。理论上，R2的取值范围在0—1，假如R2越接近于1，说明相关性越高，若越接近于0，说明相关性越小。率定期与验证期沙塘川实测与模拟月径流总氮、总磷负荷变化如图3、图4 所示，通过参数调整，使得模拟值与实测值趋于一致。由表3 可知，月径流量率定期和验证期R2达到0.84以上，NS 达到0.68 以上，日径流量率定期和验证期R2达到0.74 以上，NS 达到0.66 以上，说明模型对沙塘川流域水文过程的模拟效果较好。月径流总氮负荷率定期和验证期R2达到0.65 以上，NS达到0.50 以上，月径流总磷负荷率定期和验证期R2达到0.72 以上，NS 达到0.65 以上，说明模型对沙塘川流域氮磷负荷的模拟效果较好，研究所建立的SWAT 模型能够较为真实地反映沙塘川流域的水文和污染物迁移转化过程，模拟结果可以满足研究需要。

表3 流域径流量及径流总氮、总磷负荷模拟评价指标Tab.3 Simulated evaluation indexes of runoff and total nitrogen and total phosphorus load in the watershed

图3 率定期、验证期沙塘川实测与模拟月径流总氮负荷变化Fig.3 Measured and simulated monthly runoff total nitrogen load changes in Shatangchuan in calibration periods and validation periods

图4 率定期、验证期沙塘川实测与模拟月径流总磷负荷变化Fig.4 Measured and simulated monthly runoff total phosphorus load changes in Shatangchuan in calibration periods and validation periods

2 结果与分析

2.1 非点源氮磷污染负荷时空分布特征

利用构建好的沙塘川流域SWAT 模型，模拟全流域内的非点源氮磷负荷，并以2015 和2016年平均值为基准，得到不同时间、不同来源、不同区域非点源氮磷进入河流的负荷量。

2.1.1 非点源氮磷污染负荷时间分布特征

沙塘川流域进入河流的非点源总氮污染负荷为236.77 t·a−1，非点源总磷污染负荷为51.02 t·a−1。根据湟水河民和水文站逐月径流量变化图（图5），将流域水期分为丰水期（7—10月）、平水期（3—6月和11月）、枯水期（12月—次年2月）。从各月污染负荷变化情况（表4）来看，非点源总氮和总磷负荷均呈现丰水期（7—10月）＞平水期（3—6月和11月）＞枯水期（12月—次年2月）的规律，这与陈俊宏等（2021-06-06）研究结论一致。总氮在各月负荷占比相对较为平均，最小月负荷占比为2.06%（2月），最大月负荷占比为17.11%（10月），枯水期、平水期和丰水期负荷占比分别为10.33%、33.34%和56.33%。总磷在各月的负荷占比差异非常大，最小月负荷占比为0%（1月），最大月负荷占比为41.18%（8月），枯水期、平水期和丰水期负荷占比分别为0.12%、13.74%和85.91%，这也说明了丰水期雨季更多的持续强降水对非点源总磷入河负荷的影响更高。

图5 民和水文站月均径流量变化图（2011—2016年）Fig.5 Variation diagram of monthly average runoff of Minhe hydrological station (2011—2016)

表4 沙塘川流域非点源总氮、总磷负荷及月变化情况Tab.4 Non-point source total nitrogen and total phosphorus load and monthly changes in Shatangchuan watershed

2.1.2 非点源氮磷污染负荷空间分布特征

沙塘川流域共划分了19 个子流域，每个子流域内进入河流的非点源总氮和总磷负荷量见图6。根据图6，从非点源总氮负荷来看，沙塘川流域从上游至下游非点源总氮负荷逐渐增大，在中上游东部流域非点源总氮负荷明显大于西部流域。流域内非点源总氮负荷贡献最高的地区位于塘川镇下游和韵家口镇的19 号子流域，非点源总氮负荷达到42.8 t·a−1；其次为位于沙塘川中下游塘川镇、东山乡、西山乡、东沟乡的12、15、16和17 号子流域，非点源总氮负荷贡献达到21.6—25.8 t·a−1；位于沙塘川中上游林川乡东部、东和乡、威远镇6、10、11 和13 号子流域总氮负荷稍低，为11.4—17.1 t·a−1。非点源总氮负荷最低的区域主要分布在沙塘川流域上游西部的南门峡镇、林川乡西部、台子乡和威远镇西部的1—5 号子流域和7—9 号子流域，这些地区的非点源总氮负荷贡献均低于10 t·a−1，且多在2.1—6.6 t·a−1。

图6 沙塘川流域进入河流的非点源总氮、总磷负荷分布示意图Fig.6 Schematic diagram of the distribution of non-point source total nitrogen and total phosphorus loads entering the river in the Shatangchuan watershed

从非点源总磷负荷来看，上游东部和下游地区是最主要的负荷区。流域内非点源总磷负荷贡献最高的地区位于上游东和乡、中游威远镇东部的11 号子流域，非点源总磷负荷达到8.15 t·a−1；其次为位于中游地区东沟乡和威远镇中部的12、13 号子流域，以及位于下游塘川镇的19 号子流域、上游南门峡镇和林川乡西部的1、2、4、5 号子流域，这些地区非点源总磷负荷在2.88—5.55 t·a−1；非点源总磷负荷最低的区域位于西山乡和东山乡的16、17 号子流域，以及台子乡北部的7、8 号子流域，非点源总磷负荷贡献在0.005—0.54 t·a−1。

根据谷歌地球在线地图及遥感土地利用分类结果，非点源总氮负荷贡献较高的区域主要分布在土地利用类型为旱地的农用地集中区，这与李明龙等（2021）的研究结果相对符合，非点源总磷负荷贡献较大区域的土地利用类型以有林地、灌木林和高覆盖度草地为主。

2.1.3 非点源氮磷流失强度

根据沙塘川流域内单位面积上非点源总氮和总磷流失强度，来识别区域内非点源总氮、总磷排放关键源区。根据图7，沙塘川流域非点源总氮流失强度较高的区域集中在中下游地区，其中威远镇中部的13 号子流域和台子乡北部的8 号子流域非点源总氮流失强度最高，分别达到2.42 kg·hm−2·a−1和1.90 kg·hm−2·a−1；其次为塘川镇和东山乡的15、16 和19 号子流域，非点源总氮流失强度在1.53—1.71 kg·hm−2·a−1；流失强度最低的区域集中在南门峡镇、林川乡西部的1、2、4、5 号子流域和东和乡、巴扎藏族乡的11 号子流域，这些区域非点源总氮流失强度在0.24—0.44 kg·hm−2·a−1。

图7 沙塘川流域非点源氮磷流失强度分布图Fig.7 Distribution map of non-point source nitrogen and phosphorus loss intensity in Shatangchuan watershed

沙塘川流域非点源总磷流失强度空间分布相比总氮有着较为显著差异，非点源总磷流失强度较高的区域集中在流域中游以及上游部分区域。除塘川镇中部和北部村镇集中的两个面积小于4 km2的14与18号子流域流失强度高于1kg·hm−2·a−1外，威远镇中部和西部所在的9 号和13 号子流域非点源总磷流失强度相对最高，达到0.85—0.95 kg·hm−2·a−1；流失强度最低的区域集中在东山乡和西山乡的16、17 号子流域，以及南门峡镇南部、台子乡北部的3、7、8 号子流域和林川乡东部的6 号子流域，这些区域非点源总磷流失强度在0.001—0.17 kg·hm−2·a−1。

2.2 关键源区识别

关键源区识别是需要分析和识别出对流域内水环境整体状况有决定性影响的污染敏感区域（Howarth et al，2002；Sharpley et al，2003；Tripathi et al，2005）。本研究构建了非点源氮磷污染关键源区识别体系，将各子流域总氮、总磷负荷以及单位面积总氮、总磷流失强度按从高到低排序，并按表5 分别划分非点源总氮、总磷高负荷区、中负荷区和低负荷区，以及高流失强度区和低流失强度区，将高负荷区的全部区域，以及中负荷区和高流失强度区共有的区域作为非点源总氮、总磷的关键源区。

表5 沙塘川流域非点源总氮总磷负荷与流失强度分级一览表Tab.5 Classification of total nitrogen and phosphorus load and loss intensity of non-point sources in Shatangchuan watershed

根据沙塘川各子流域总氮和总磷负荷以及流失强度分级排序情况（图8），沙塘川流域关键源区识别结果如图9 所示。本研究确定沙塘川流域非点源总氮的关键源区为位于威远镇中部、塘川镇、东山乡、西山乡的13、15、16、17 和19 号子流域，这些子流域面积占全流域的30.8%，非点源总氮负荷贡献量占全流域的55.2%；非点源总磷的关键源区为11、12、13 和19 号子流域，这些子流域面积占全流域的38.9%，非点源总磷负荷贡献量占全流域的47.1%。

图8 沙塘川流域非点源总氮、总磷负荷及流失强度排序分级图Fig.8 Ranking map of total nitrogen,total phosphorus load and loss intensity of non-point sources in the Shatangchuan watershed

图9 沙塘川流域非点源总氮总磷关键源区识别结果Fig.9 Identification results of key source areas of non-point source total nitrogen and phosphorus in the Shatangchuan watershed

3 结论

本研究通过构建沙塘川流域SWAT 非点源污染模型，以验证期2015 和2016年的平均值为基础，定量估算了流域非点源氮磷污染负荷，并对其时空分布特征和来源进行了分析。在此基础上，通过构建关键源区识别体系，对沙塘川流域非点源总氮、总磷进行了关键源区识别，为后续制定流域非点源高效管控措施提供了基础。

（1）沙塘川进入河流的非点源总氮污染负荷为236.77 t·a−1，非点源总磷污染负荷为51.02 t·a−1，非点源总氮和总磷负荷均呈现丰水期＞平水期＞枯水期的规律。

（2）非点源总氮负荷从沙塘川流域上游至下游逐渐增大，在中上游流域东部流域非点源总氮负荷明显大于西部流域。流域内非点源总氮负荷贡献最高的地区位于塘川镇下游和韵家口镇的19 号子流域，负荷达到42.8 t·a−1，非点源总氮负荷最低的区域为1—5 号子流域和7—9 号子流域，非点源总氮负荷贡献均低于10 t·a−1，且多在2.1—6.6 t·a−1。总体来看，沙塘川流域总氮贡献较高的区域主要分布在土地利用类型为旱地的农用地集中区，总氮贡献较低的区域主要分布在南门峡水库及沙塘川支流城镇区域。

（3）从非点源总磷负荷的空间分布来看，总磷贡献较大区域主要集中在沙塘川上游东部地区和下游地区，该区域土地利用类型以有林地、灌木林和高覆盖度草地为主，畜禽养殖是总磷的主要来源。

（4）沙塘川流域非点源总氮的关键源区为13、15、16、17 和19 号子流域，这些子流域面积占全流域的30.8%，非点源总氮负荷贡献量占全流域的55.2%；非点源总磷的关键源区为11、12、13 和19 号子流域，这些子流域面积占全流域的38.9%，非点源总磷负荷贡献量占全流域的47.1%。

4 展望

本研究通过对沙塘川流域非点源氮磷污染负荷定量分析及关键源区识别，初步形成了一套适用于青海东部地区的非点源污染关键源区识别技术方法。下一步将在继续开展流域高分辨率非点源污染清单研究和非点源污染高效管控措施优化筛选研究的基础上，将流域非点源污染高效管控措施进一步推广到全湟水流域，为湟水流域乃至青海省非点源污染管控工作提供科学支撑。