基于SWAT模型的龙溪河泸县境内面源污染特征分析

马亚丽，白祖晖，敖天其

(1.甘肃农业大学水利水电学院，兰州 730000；2.甘肃省水利水电勘测设计研究院有限责任公司，兰州 730000；3.四川大学 水力学与山区河流开发保护国家重点实验室，成都 610065)

相关研究表明，SWAT模型在河灌区、流域面积仅为6.26 km2的重庆市万州区陈家沟小流域、柴河水库流域、三峡库区等多个流域[7-10]的面源污染负荷模拟研究都可以取得较好效果，因此，本文首次利用SWAT模型在龙溪河泸县境内的面源污染进行研究，探讨SWAT模型在该流域的适用性，并建立适合该流域的SWAT模型数据库及确定敏感性参数值，定量计算该流域的氮磷污染负荷，分析氮磷污染负荷特征。其成果将对该流域的面源污染的控制和管理、水源地的保护、城市饮水安全以及流域水资源管理有一定的理论和实践意义，可以作为龙溪河流域水资源保护、规划和面源污染控制、规划、决策的依据。

1 流域概况

龙溪河发源于四川省登东山，属于长江左岸一级支流，总干流长110 km，流域面积502 km2，途径泸县立石、云锦、玄滩、奇峰、兆雅、云龙等镇以及龙马潭区长安、石洞、鱼塘、罗汉等镇，在龙马潭区汇入长江，其中，龙溪河泸县境内干流长约69 km，流域面积242 km2，涉及人口14.60万人。多年平均流量为6.58 m3/s，平均深度为0.4 m，河宽40～60 m，河道地形起伏不大，平均比降为0.17%。研究区四季分明，具有亚热带湿润气候的特点。全县多年平均降水量为1 013.6 mm，降雨量最大的年份可达1 450.2 mm；多年平均风速1.2 m/s，一般以偏北风为主；多阴天，日照少；多年平均温度为17.7 ℃；海拔在316～376 m左右，地形起伏不大。

2 数据来源及处理

SWAT模型数据库包括空间数据库和属性数据库，空间数据包括高程地形图(DEM)、土地利用图、土壤类型图、水系图，属性数据包括土壤、气象、水文、水质、污染源、农业管理方式等，其中空间数据必须具有统一的投影坐标[11]，本研究采用Albers等积圆锥投影，相应投影参数见表1。由于泸县环保局委托编制《濑溪河流域(泸县境内)水污染评估项目成果报告》、《泸州市龙溪河生态环境泸县段(2013-2015年)保护方案》，由泸县环保局协调收集获得泸县相关资料。各种数据的精度、格式和来源具体见表2。

表1 研究区投影参数表Tab.1 The projection parameters of the study area

表2 模型所需数据分辨率、格式和来源Tab.2 The resolution, format, and source of the data required for the model

2.1 流域内子流域划分

SWAT模型基于DEM图以及已有的水系图生成最终水系以及划分子流域。利用SWAT模型进行流域离散化，首先根据DEM生成的河网分布和出水口的位置，将整个流域划分为若干个子流域，其中，最小河道集水面积阈值设置得越小，生成的水系就越详细，子流域划分数目就越多，然后根据每个子流域内的土地利用、土壤和坡度属性叠加划分出一个或多个水文响应单元。HRU是SWAT模型模拟的基本单位，假定认为具有统一的水文效应。本研究中设定最小河道集水面积阈值400 hm2，以龙溪河在泸县的出口断面作为流域断面划分17个子流域，具体子流域划分见图1。

图1 子流域划分图Fig.1 The sub-watershed division diagram

2.2 各子流域中水文响应单元(HRUs)分配

根据土地利用类型、土壤类型和坡度的分类，将每个子流域划分为一个或多个水文响应单元(HRUs)，HRU作为SWAT模型模拟的基本单位，计算每个HRU的径流量、输沙量和营养物质污染负荷。为适应SWAT模型非点源模拟中土地利用类型的分类，根据流域内原来的分类，土地利用类型重新分类编码[12]，土壤类型图与土地利用类型图的处理方法相似。将重分类的土地利用数据、土壤数据以及坡度分类进行叠加，本研究中设置土地利用、土壤和坡度的阈值分别为20%、25%、15%，最终将流域划分为65个水文响应单元。

2.3 模型参数的写入及编辑

建立各类气象信息测站位置表和实测数据表，写入各类气象数据，创建模型所需的各种输入文件，包括：结构文件(.fig)，土壤文件(.soil)、气象文件(.wgn)，子流域文件(.sub)、水文响应单元文件(.hru)、主河道文件(.rte)、地下水文件(.gw)、水利用文件(.wus)、农业管理文件(.mgt)、土壤化学文件(.chm)、池塘数据文件(.pnd)和河流水质文件(.swq)，输入完毕后可以通过模型界面修改编辑。

3 模型参数率定与验证

模型参数的率定，从空间上来说，流域从上到下，从支流到干流；从时间上来说，时段从粗到细；从校准的内容来说，先径流后污染负荷的顺序，依次校准。SWAT模型可以进行人工校准也可以进行自动校准，本研究区采用两种方法相结合进行参数校准。在确定模型的敏感性参数之后，利用模型自带的自动校准分析方法SCE-UA (Shuffled Complex Evolution Algorithm )来确定这些敏感性参数的初始取值，将模型自动校准优化后的参数代入模型模拟，比较分析模拟值与实测值，如果3个校准指标满足要求，采用自动校准之后敏感性参数值；如果校准指标不满足要求，则需要对敏感性参数采用人工校准的方法进行相应的调整，直到评价指标在允许的范围之内。

采用龙溪河泸县出口监测断面的2006-2010年实测流量、负荷数据，以2005年的数据作为预热，2006-2008年数据进行率定，2009-2010年的数据进行验证。

3.1 参数敏感性分析

采用模型自带的敏感性分析模块对径流和氮、磷负荷进行敏感性分析，选取几个排序在前的参数作为调整参数，确定敏感性参数的最优值，从而得到很好的模拟结果，并且节省工作量。敏感性参数排序见表3。

3.2 径流的调试与验证

评价模型结果是否能够被接受，不存在统一的标准，一般认为，当R2>0.6，Ens>0.5，Re≤0.3%时，模型是准确的[13]。月均径流率定及验证结果见表4，模拟值与实测值月均流量对比如图2所示。径流模拟结果可以被接受，引起模拟值与实际值存在的差异的原因是多方面的，包括流域的空间差异性分布所带来的参数难确定性，以及模型结构本身的模拟误差和数据的准确性等原因。

表3 敏感性参数表Tab.3 The sensitivity parameters table

表4 径流模拟率定及验证结果Tab.4 The runoff simulation results of calibration and verification

3.3 氮、磷污染负荷的调试与验证

与径流模拟相比，氮磷负荷的模拟影响因素更多，不确定性将更大，因此，氮磷负荷模拟的精度要求低于径流模拟的要求。氮磷负荷率定及验证结果见表5，模拟值与实测值TN月日均负荷对比见图3，TP月日均负荷对比见图4。氮磷负荷模拟值与实测值差异，与径流模拟的精确程度密不可分，同时污染负荷的迁移转化过程所引起的不确定性也是造成模拟值与实测值不一致的重要因素。

图2 月均流量对比图Fig.2 The comparison of average monthly flow

图3 TN月日均负荷对比图Fig.3 The comparison figure of TN monthly average load

图4 TP月日均负荷对比图Fig.4 The comparison figure of TP monthly average load

表5 氮磷负荷模拟率定及验证结果Tab.5 The calibration and verification results of nitrogen and phosphorus load simulation

4 面源污染的特征分析

4.1 面源污染时间变化规律

依据研究区3个气象站的降雨数据以及泸县1960-2010年多年平均降雨资料，将一年分成3个时期进行分析，5-9月作为丰水期，占年降水量的73%，12-3月作为枯水期，占年降水量的9%，其他月份作为平水期。从5年的产量均值，可以看出产水量、TN产量、TP产量都呈现丰水期显著多于平水期，平水期多于枯水期，丰水期产量比重大约为70%～90%，枯水期产量比重在5%以内，可见，丰水期是面源污染产生的关键时期，决定了面源污染的严重程度。这与面源污染的产生机理有关，面源污染是以径流冲刷为主要的动力来源，降雨量大的月份造成的土壤侵蚀就严重，随之携带的氮磷负荷量就大。具体见表6。

表6 研究区各时期产量情况表Tab.6 The production situation in each period of the study area

4.2 面源污染空间分布规律

从空间角度分析面源污染的分布情况，子流域5、9、11、13的TN输出系数最大，为5 599～5 903 kg/km2，且东部地区产量比西部地区大；子流域12、13、14的TP输出系数最大，为184～197 kg/km2，且下游地区比上游地区大。综合来看，5、11、13三个子流域氮磷负荷都较为严重的，4号子流域最轻，TN负荷东部地区较西部地区更严重，TP负荷下游比上游更严重。这不仅与子流域的下垫面差异性有关，还与各子流域污染排放情况密切相关。具体见图5、图6所示。

图5 年均TN产量分布Fig.5 The annual yield distribution of TN

图6 年均TP产量分布Fig.6 The annual yield distribution of TP

4.3 重要污染源识别

SWAT模型是基于流域水文地理特征的分布式水文模型，能够模拟出各子流域污染负荷总产量，但是污染确定污染来源，本文为了确定点源与面源污染哪个是主要污染来源，哪种面源污染类型比重较大，采用将排放量与模拟结果进行对比分析的方法。

4.3.1 主要污染源识别

从实际污染排放量与年均模拟值可以看出，面源污染是各子流域的主要污染排放类型，点源污染比重很小，且面源污染排放量大小与年均模拟值大小一一对应，排放量大的子流域相应的模拟值也大，因各子流域面源污染排放量明显高于点源污染排放量，是各个子流域最主要的污染类型，应该重点防治面源污染。具体见图7、图8。

图7 TN污染源类型识别Fig.7 the source type identification of TN

图8 TP污染源类型识别Fig.8 the source type identification of TP

4.3.2 主要面源污染类型识别

本文研究区考虑的面源污染类型包括农村生活污染、畜禽养殖污染、化肥污染，从TN的面源污染分布看，10、13、14子流域面源污染排放量最大，各面源污染排放比重分别为：农村生活污染50%、化肥污染40%、畜禽养殖污染10%；从TP面源污染分布看，10、13、14子流域面源污染排放量最大，相应的子流域年均模拟产值也最多，各面源污染排放比重分别为：化肥污染55%、农村生活污染23%、畜禽养殖污染22%，可见，化肥污染对面源污染产生影响最大，其次为农村生活污染和畜禽养殖污染。农村生活污染对TN模拟值影响最大，化肥污染对TP模拟值影响最大。具体见图9、图10。

图9 TN面源类型识别图Fig.9 the non-point source identification of TN

图10 TP面源类型识别图Fig.10 the non-point source identification of TP

4.4 污染关键区域识别

依据年均模拟TN、TP值分别分为5级，TN：Ⅰ级41～45 kg/hm2，Ⅱ级45～49 kg/hm2，Ⅲ级49～53 kg/hm2，Ⅳ级53～57 kg/hm2，Ⅴ级57～61 kg/hm2；TP：Ⅰ级1.0～1.2 kg/hm2，Ⅱ级1.2～1.4 kg/hm2，Ⅲ级1.4～1.6 kg/hm2，Ⅳ级1.6～1.8 kg/hm2，Ⅴ级1.8～2.0 kg/hm2。具体见图11、图12。

图11 TN年均产量分级Fig.11 The classification of TN annual output

图12 TP年均产量分级Fig.12 The classification of TP annual output

子流域3、9、11、13、15的TN负荷属于Ⅴ级，子流域产量最大；子流域11、12、13、14、17的TP负荷属于Ⅴ级，因此，11、13是非点源污染的关键区，是值得重点治理的区域。

5 结 语

本文首次以龙溪河泸县境内的面源污染为研究对象，以地理信息系统(GIS)和SWAT模型为工具，探讨了SWAT模型在研究区的适用性；定量计算研究区内面源污染负荷，并分析了其特征，主要结论如下。

(1)采用5年资料径流和水质监测数据，选取2000-2005年的数据作为预热期，2006-2008年数据进行率定，2009-2010年数据进行验证，依据相对误差(Re)、相关系数(R2)和Nash-Suttcliffe系数(Ens)这3个指标作为模型校准的评价指标，来表征模型模拟结果的优劣，结果表明，模拟值的精度在允许范围之内，敏感性参数的取值具有一定的可信性，所建立的SWAT模型在研究区具有适用性。

(2)年产水量、TN产量、TP产量三者之间存在着一定的相关关系，均呈现丰水期显著多于平水期，平水期多于枯水期，丰水期产量比重大约为70%～90%，枯水期产量比重在5%以内，可见丰水期是面源污染产生的关键时期，决定了非点源污染的严重程度。

(3)各子流域的污染排放量与模拟年均总值是基本对应的，并且面源污染是各个子流域最主要的污染类型。化肥污染与农村生活污染在面源污染中所占比重更大，影响也更大，对于TN模拟值影响最大的是农村生活污染，对TP模拟值影响最大的是化肥污染。

(4)从年均模拟值看，由于11、13是面源污染的关键区，位于云锦镇，是值得重点关注与治理的区域。

作为分布式水文模型代表的SWAT模型，尽管考虑了时空的异质性，将每个子流域划分为多个水文响应单元进行模拟，但模拟单元划分粗略决定输入数据向模拟单元尺度转换时参数集总程度的高低[14]，同时，输入数据的分辨率及精度、SWAT模型结构本身的不完善、模型参数难以准确获取都造成模拟结果的不确定性，需要后续进一步研究。