1 引言(Introduction)

水环境是社会与经济系统健康稳定发展的重要基础, 随着经济方式的多样化发展, 城市化进程加快, 污染负荷超过水环境容量而导致的水质恶化已造成富营养化、酸化和缺氧等水环境问题, 从而制约了经济社会可持续发展(Xu et al., 2018).目前, 我国流域污染物排放总量控制是最为重要的流域管理措施之一, 对水环境污染负荷的估算和时空分析成为水环境健康诊断、预测、预警及控制和管理的重要基础.

污染负荷估算是对流域污染负荷时空变化进行定量分析的前提与基础.目前, 国内外学者常用输出系数模型(龙天渝等, 2008；Cai et al., 2018)估算每个计算单元(人、畜禽或单位土地面积)的污染负荷.该模型忽略了非点源污染复杂的迁移转化过程, 计算区域既可以是边界明确的流域, 也可以是不同等级的行政单元, 时间步长的设定比较灵活.虽然测算精度通常比实证模型(国外应用较多, 数据要求高, 可移植性差)和SWAT、HSPF、AnnAGNPS等机理模型低(耿润哲等, 2013；Fan et al., 2017), 但对尺度不敏感, 可移植性强, 可以在较大尺度和较长时间段对非点源污染负荷进行估算(刘庄等, 2015).同时, 因其结构简单和数据易获取的特点在国内得到广泛应用(薛利红等, 2009).目前, 国内在输出系数法方面的研究大多是对经典输出系数模型(Johnes, 1996)进行改进, 如排污系数法就是借鉴经典输出系数模型的思路, 将污染物输出系数替换为排放系数, 得到污染物排放量的一种污染负荷估算方法(熊昭昭等, 2018), 对点面源污染负荷都可进行估算, 具有参数较少、应用范围广、适用于缺乏实验条件的区域的优点.

在合理估算污染负荷的基础上, 对污染负荷进行时空变化分析可以直观反映流域内各污染源污染物在一定时空内的排放特征, 该变化特征既受土地利用方式和人类社会经济活动的影响, 也受各污染源之间相互作用的影响.因此, 如何探究流域污染负荷时空变化与土地利用(包括各农作物种植类型的改变)、社会经济间的关系成为解决流域水环境污染的关键.目前对污染负荷的时空分析多以源解析和直接分析污染强度的空间差异(胡芸芸等, 2015；谢经朝等, 2019)为主, 缺乏对影响污染负荷时空变化原因的定量分析.多元统计分析法具有强大的数据分析功能, 能够定量地探究数据之间的内在联系, 但其在污染负荷时空分布方面的应用研究较少, 更多地应用于水质时空变化方面的研究.例如, 如富天乙等(2014)和杜麦等(2017)运用因子分析和层次聚类分析及多元线性回归分析等方法分别对辽阳太子河和浐灞河的水质进行了评价, 充分利用水质指标数据的时间和空间特征, 识别出监测样本的内在差异性和相似性.张涛等(2017)运用方差分析法研究了水质的空间差异性和季节性变化.以上多元统计分析方法多用于对水质监测数据进行分析来反映流域受污染状况, 但水质监测数据往往较难获取, 以污染负荷为基础数据进行多元统计分析可以弥补水质监测数据较难获取这一不足, 也可以达到真实反映流域受污染的情况和探究污染负荷时空变化驱动因素, 进而为水环境质量时空变化的内在机理分析提供数据支撑的目的.

沱江流域既是四川省工业和人口集中之地, 又是四川省最大的棉、蔗产地.近年来, 在人口与经济发展的双重压力下, 沱江流域的水质现状不容乐, 水环境问题逐渐显现(汪嘉杨等, 2017).目前, 对沱江流域的水污染研究主要集中在污染源解析(胡芸芸等, 2015)、污染物排放现状和防控(秦延文等, 2020)及水质时空变化特征分析(张明锦等, 2020)等方面, 缺乏污染负荷时空变化特征的研究.因此, 本文采用排污系数法估算沱江流域污染负荷, 并结合社会经济和土地利用方式对流域污染负荷的时空变化特征进行多元统计分析, 以期客观地反映沱江流域2007—2017年污染负荷变化的时空分布及变化特征与社会经济变化因素和土地利用方式之间的相互作用关系, 从而为该流域的水污染治理及污染防控提供可靠的数据支撑和理论依据.

2 数据与方法(Data and methods) 2.1 研究区概况

沱江是长江的一级支流, 发源于四川茂县九顶山南麓, 全长627.4 km, 流域面积为2.76×10⁴ km², 其中, 四川境内流域面积为2.55×10⁴ km², 是四川省面积最大的流域, 地理坐标为北纬28°49′17″~31°42′3″, 东经103°40′57″~105°42′50″(图 1), 流域多年平均降水量为1200 mm, 年径流量为351×10⁸ m³, 其中, 岷江补给约占33.4%, 是成都平原经济区、川西北生态经济区和川南经济区等主要经济区的重要水源地(汪嘉杨等, 2017).2007—2017年, 沱江流域的土地利用类型以耕地和林地为主, 分别约占整个流域面积的78%和14%, 且变化不大(耕地减少0.98%, 林地减少1.55%), 建设用地面积由759 km²增加到987 km², 增加了30%(图 2).沱江流域包含的28个区县隶属于成都市、德阳市、眉山市、资阳市、内江市、自贡市和泸州市等社会经济发达、人口密集的地级市, 地处成渝经济圈腹心地带(仇开莉等, 2013), 区位优势明显.

2.2 数据来源

在2017年8月四川省生态环境厅出台《沱江流域水污染防治规划(2017—2020年)》以前, 沱江流域水质长期处于劣Ⅲ类且有恶化的趋势.为了展示这种恶化趋势在时空尺度上的特征, 并结合社会经济和土地利用类型数据的可获取性, 本文特选取2007和2017年作为典型年份, 研究流域污染负荷的时空分布特征, 采用的数据来源及其特征如表 1所示.

2.3 研究方法 2.3.1 排污系数法

本研究采用排污系数法分别对沱江流域28个区县的城镇生活、畜禽养殖、农村生活垃圾、农村生活、农田固废类、农田径流和城市径流污染源中COD、NH₃-N、TN和TP的污染负荷进行计算.沱江流域土地利用类型以农田为主, 约占整个流域的78%(图 2), 近年来, 在经济利益的驱使下, 区域农田内各农作物种植类型变化显著, 但各农作物种植类型的耕作模式、管理状况等因素几乎保持不变.因此, 本文采用统一排污系数估算2007年和2017年沱江流域污染负荷, 排污系数的确定通过参考《第一次全国污染源普查产排污系数手册》(国务院第一次全国污染源普查领导小组办公室, 2009)和文献检索法(Bramer et al., 2018)确定.各污染源污染负荷的详细计算信息如表 2所示.

2.3.2 多元统计分析法

① Pearson相关分析法：如果变量之间存在显著关系, 则可以从变量之间的高度相关性来理解(Sezen et al., 2020).Pearson相关系数利用变量间的线性相关性来表示用户相似度, 其计算公式见式(2).

(2)

式中, sim(X, Y)为X、Y两个指标间的相似度, 即相关系数, i为某种指标的第i个样本, n为样本总数.农用化肥施用量、农业产值、农业人口是影响农业面源污染负荷的重要因素, 工业产值、非农业人口则是影响流域点源污染负荷的重要因素.因此, 本文选择上述5个社会经济指标与4种污染负荷进行相关性分析.

② 主成分分析：研究如何通过少数几个主成分来揭示多个变量间的内部结构, 即从原始变量中导出少数几个主成分, 通常数学上的处理就是将原来P个指标进行线性组合作为新的综合指标(林海明等, 2013), 使它们尽可能多地保留原始变量的信息, 且彼此间互不相关的一种多元统计方法.本文以沱江流域内各污染源为变量确定主成分个数, 从而确定聚类个数.

③ K-means聚类法：该方法属于快速聚类, 计算速度快, 适用于样本较多的研究, 对于n个数值变量参与快速聚类, n个变量组成一个n维空间, 每个样品是空间中的一个点, 最终按照要求聚类成K个类别.聚类前计算机随机产生初始的聚类中心, 然后计算机迭代新的聚类中心(武松, 2019).本文对流域内28个区县各污染源污染负荷进行K-means聚类, 揭示污染负荷时空分布的相似性和差异性.

3 结果(Results) 3.1 4种污染负荷时间变化特征

由图 3可知, 2007—2017年, 沱江流域COD、NH₃-N、TN和TP的污染负荷均有所增加, 分别增加了87556.54、3743.45、13807.68和1637.14 t, 其中, 增长量最多的是COD, 最少的是TP.从增长率来看, COD的增长率最高, 达到了24.16%, TN的增长率为20.75%.

3.2 4种污染负荷空间变化特征

由图 4可知, 2007、2017年沱江中游大部分区县(简阳市、仁寿县和雁江区等)和上游的旌阳区、彭州市及下游的泸县、富顺县和龙马潭区的4种污染负荷均较大, 污染较为严重, 而上游的什邡市、青白江区和下游的大安区、沿滩区和自流井区的4种污染负荷较小, 污染较轻.2007—2017年, 仅乐至县的NH₃-N污染负荷明显减少, 减少到1400 t以下；其余污染负荷发生明显变化的区县均为增加, 如中游的雁江区、资中县和上游的旌阳区、新都区的COD污染负荷达到20000 t以上, 中游的东兴区、金堂县和上游的彭州市、龙泉驿区及下游的富顺县、泸县的COD污染负荷达到16000 t以上；上游的绵竹市、新都区和龙泉驿区的NH₃-N污染负荷达到1000 t以上, 上游的什邡市和中游的内江市市中区及下游的沿滩区、龙马潭区的NH₃-N污染负荷达到800 t以上；中游的雁江区、资中县和下游的泸县的TN污染负荷达到3800 t以上, 中游的东兴区、金堂县和下游的富顺县的TN污染负荷达到3000 t以上, 中游的威远县、隆昌县和下游的龙马潭区、江阳区及上游的广汉市、新都区和龙泉驿区的TN污染负荷达到2300 t以上, 中游的内江市市中区和下游的沿滩区及上游的什邡市的TN污染负荷达到1600 t以上；上游的旌阳区和下游的泸县的TP污染负荷达到550 t以上, 中游的荣县和上游的新都区的TP污染负荷达到350 t以上, 中游的内江市市中区及上游的什邡市、龙泉驿区的TP污染负荷达到250 t以上.此外, 由图 4和图 2对比可知, 流域上游的绵竹市、什邡市和青白江区及下游的大安区、沿滩区、贡井区和自流井区的污染负荷较少, 这些区县的土地利用类型主要为林地和草地, 而其他区域污染负荷增加量较多, 污染严重区域的土地利用类型主要以农田和建设用地(城市密集区域)为主.

由图 5a可知, 流域内仅中游乐至县的COD污染负荷减少308.6 t, COD污染负荷增加较多的区县有位于上游的新都区、龙泉驿区、旌阳区、广汉市和中游的雁江区、资中县及下游的泸县、龙马潭区、江阳区和沿滩区, 增加量最多的是龙泉驿区, 达到7988.15 t.由图 5b可知, NH₃-N污染负荷减少的有金堂县、仁寿县、乐至县和简阳市, 这4个区县集中在中游, 其中, 减少量最多的是乐至县, 达到133.45 t, 最少的是简阳市, 减少量为21.06 t；NH₃-N污染负荷增加较多的有上游的新都区、龙泉驿区、广汉市和中游的雁江区、资中县及下游的沿滩区、泸县、龙马潭区和江阳区, 增加量最多的是龙泉驿区, 达到469.83 t.由图 5c可知, 流域内所有区县的TN污染负荷均增加, 其中, 污染负荷增加较多的有上游的新都区、龙泉驿区和中游的雁江区、资中县及下游的泸县、龙马潭区, 增量最多的是龙马潭区, 为1556.08 t.由图 5d可知, TP污染负荷减少的仅有青白江区, 减少量为0.71 t；增加较多的有上游的什邡市和中游的仁寿县、雁江区和隆昌县及下游的泸县和龙马潭区, 其中, 增量最多的是泸县, 为181.93 t.由此可见, 上游的新都区、龙泉驿区和下游的泸县、龙马潭区的4种污染负荷增加量都较大.

3.3 污染负荷与社会经济指标的相关性分析

本文以Pearson相关分析为研究手段, 选取沱江流域28个区县的农业人口、非农业人口、农用化肥施用量、工业产值和农业产值这5个社会经济指标, 分别与2007、2017年沱流域28个区县的4种污染负荷总量进行相关性分析, 结果如表 3所示.

由表 3可知, 2007—2017年, 农业人口与COD、NH₃-N、TN和TP污染负荷的相关性由高度相关变为强相关或中等程度相关, 相关性降低.而非农业人口与COD、NH₃-N、TN和TP污染负荷的相关系数分别由0.460、0.218、0.226和0.184变为0.953、0.938、0.881和0.871, 由中等强度相关或不相关变为高度相关.农用化肥施用量和农业总产值与污染负荷的相关性变化不大, 相关系数基本大于0.6, 为强相关, 工业产值与污染负荷的相关性则呈弱相关或不相关.

3.4 各污染负荷的污染源聚类分析

将流域内28个区县的COD、NH₃-N、TN和TP污染负荷的8种污染源分别采用主成分分析降维提取3个主成分, 分别可以提取总污染源信息的91.86%、84.81%、84.87%和84.09%, 均高于80%.因此, 本文运用K-means聚类将沱江流域内28个区县4种污染负荷的8种污染源聚为3类.

根据表 4和表 5中各污染源在各聚类中心的贡献大小来命名各污染源聚类类型, 结果表明, COD、NH₃-N和TN具有相同的污染源聚类类型, 即城镇污染源主导型、农村污染源主导型、生活-畜禽养殖污染源复合型.而TP的污染源聚类类型为农村-城镇生活污染源复合型、畜禽养殖-城市径流-工业污染源复合型、城镇生活-畜禽养殖污染源复合型.

本文进一步利用ArcGIS将聚类结果可视化, 不同污染负荷下各污染源的聚类空间分布变化特征如图 6所示, 4种污染负荷的污染源聚类结果均有明显的空间分布变化特征.由2017年的COD污染源聚类结果可知, 上游的绵竹市、什邡市、广汉市、青白江区和中游的金堂县、乐至县、威远县、内江市市中区及下游的部分区县(荣县、贡井区、自流井区等)为生活-畜禽养殖污染源复合型, 其余区县为城镇污染源主导型.上游的彭州市、新都区、龙泉驿区和中游的部分区县(简阳市、雁江区、仁寿县等)及下游的部分区县(隆昌县、泸县、富顺县等)的COD污染结构发生变化.

由2017年的NH₃-N污染源聚类结果可知, 上游的彭州市、新都区、旌阳区、龙泉驿区和中游的部分区县(金堂县、简阳市、仁寿县等)及下游的富顺县、泸县、江阳区为生活-畜禽养殖污染源复合型, 其余区县为城镇污染源主导型.上游的彭州市、新都区、龙泉驿区和中游的部分区县(金堂县、简阳市、仁寿县等)及下游的富顺县、泸县、江阳区的NH₃-N污染结构发生变化.

由2007年的TN污染源聚类结果可知, 上游的旌阳区和中游的部分区县(金堂县、简阳市、雁江区等)及下游的富顺县、泸县为生活-畜禽养殖污染源复合型, 其余区县为农村污染源主导型.2017年新增了城镇污染源主导型, 上游部分区县(彭州市、旌阳区、广汉市等)和中游的东兴区及下游的隆昌县、沿滩区、江阳区的TN污染结构发生变化.

由2017年的TP污染源聚类结果可知, 上游部分区县(绵竹市、什邡市、旌阳区等)和中游的金堂县、乐至县、威远县、东兴区及下游的荣县、贡井区、隆昌县、富顺县为畜禽养殖-城市径流-工业污染源复合型, 上游的新都区、青白江区、龙泉驿区和中游的内江市市中区及下游的部分区县(大安区、贡井区、沿滩区等)为城镇生活-畜禽养殖污染源复合型, 其余区县聚为农村-城镇生活污染源复合主导型.什邡市、雁江区、资中县、隆昌县、泸县的TP污染结构发生变化.

2007年COD、NH₃-N、TN的污染源聚类结果具有相似的聚类空间分布特征, 污染类型主要为农村污染源主导型和生活-畜禽养殖污染源复合型, 且农村污染源主导型主要分布在中游, 生活-畜禽养殖污染源复合型主要分布在上游和下游.而2017年只有COD和NH₃-N的污染源聚类结果表现为相似的聚类空间分布特征, 均表现为农村污染源主导型、生活-畜禽养殖污染源复合型大范围向城镇污染源主导型转变, TN的污染源聚类结果仅上游部分区县(旌阳区、广汉市、彭州市等)和中游的东兴区、隆昌县及下游的沿滩区、江阳区表现为农村污染源主导型、生活-畜禽养殖污染源复合型向城镇污染源主导型转变, 其余区县聚类类型未发生变化.2007年TP的污染源聚类结果则与以上3种污染负荷的聚类结果时空分布特征具有明显差异, 畜禽养殖-城市径流-工业污染源复合型在上、中、下游均有分布, 城镇生活-畜禽养殖污染源复合型则主要分布在上游和下游, 农村-城镇生活污染源复合主导型主要分布在中游.2017年TP的污染类型主要由城镇生活-畜禽养殖污染源复合型、农村-城镇生活污染源复合主导型向畜禽养殖-城市径流-工业污染源复合型转变.

4 讨论(Discussion) 4.1 流域污染负荷时空分布及变化特征的驱动因素

本研究通过排污系数法计算出沱江流域4种污染物的污染负荷, 2007—2017年, 沱江流域4种污染物的污染负荷时间上呈整体增加趋势, 空间上呈现异质性的变化特征.城市人口数量增多导致城镇生活污染负荷增加, 城镇生活污染源对沱江流域污染排放的贡献率提高, 使得非农业人口数量与污染负荷的相关性变强.并且城市化进程也会影响土地利用方式(Elizabeth et al., 2004), 2017年沱江流域建设用地面积相较于2007年增加了228 km²(图 2), 城市径流污染加剧.随着流域各城市进入城镇化、工业和农业现代化加速阶段(汪嘉杨等, 2017), 城镇生活污染负荷不断增加使得COD排放量变大, 其增长率达到24.16%.由于工业产值是影响流域点源污染负荷的主要因素之一, 随着沱江流域管网系统的完善和治理点源政策法规的颁布, 点源污染已得到较好的控制.另外, 沱江流域的农业发达, 农业面源污染是其主要污染源, 点源污染对总污染负荷的贡献较小, 故工业产值与该流域总污染负荷的相关性弱.

2007年和2017年, 沱江中游大部分区县(简阳市、仁寿县、雁江区等)和上游的旌阳区、广汉市、彭州市及下游的泸县、富顺县、龙马潭区的4种污染负荷均较大, 污染较为严重.原因是位于中游的区县均有大面积农业用地分布, 而氮、磷的负荷主要来自农业土地(Fan et al., 2015), 致使农业面源污染严重；位于上游和下游的区县则是工业发达和人口密集区, 旌阳区和广汉市所属的德阳市承担了三线建设时期的东电、二重、东汽等国家重要的工业项目建设, 工业污染严重, 且工业集中的地区有较多的建设用地, 城市径流污染严重(李彦旻等, 2020), 而富顺县和龙马潭区人口密集, 生活污染严重.上游的绵竹市、什邡市、青白江区及下游的大安区、沿滩区、贡井区、自流井区的污染负荷较少, 这些区县对应土地利用数据(图 2)中草地和林地分布区域.由于林地和草地对水体和土壤中的污染物有一定的吸收作用, 污染负荷与林地和草地呈负相关(孟晓云等, 2012；张彩玲等, 2017), 故该区域污染较轻.上游的新都区、龙泉驿区和中游的雁江区、资中县及下游的泸县、龙马潭区的4种污染负荷增加量都较大, 新都区和龙泉驿区都属于成都市, 社会经济发展迅速, 污染负荷增加较多, 而下游龙马潭区和泸县的城镇化率分别为增加了39.81%和31.08%, 城市化进程较快, 加剧了污染负荷的增加.

4.2 各污染负荷的污染源空间聚类特征

2007年COD、NH₃-N和TN的污染源聚类空间分布特征相似, 沱江中游地区主要以农村污染源主导型为主, 上、下游主要以生活-畜禽养殖污染源复合型为主.其原因是沱江中游有大面积农田分布, 农业面源污染严重, 而上游和下游是建筑用地分布广、人口密集区域, 生活污染严重, 畜禽养殖污染源是沱江流域农业面源污染的主要贡献源(胡芸芸等, 2015).聚类结果发生变化的区县多数集中在有农业用地和建设用地大面积分布的区域, 其社会经济发展迅速, 并且结合表 6分析可知, 聚类结果发生变化的大部分区县城镇化率或二、三产业占生产总值比的变化量较大, 如中游的乐至县、雁江区、东兴区等, 人口结构和经济结构的变化使污染结构发生变化, 聚类结果随之发生变化(杨耿等, 2018).个别区县的城镇化率或二、三产业占生产总值比的变化量虽然较大, 但其聚类结果未发生改变, 如下游的龙马潭区, 城镇化率增加了39.81%, 原因可能是这些地区人口基数小, 人口转化并未对污染负荷产生较大影响.未发生变化的区县对应土地利用数据(图 2)中林地和草地的分布区域, 有林地和草地大面积分布的区县人类活动小, 社会经济发展相对缓慢, 各污染源对污染负荷的贡献率变化小, 污染结构不易改变.

4.3 本研究对流域水环境管理的指导意义

本文提出了基于排污系数法-多元统计法-空间分析技术的研究框架, 探究了流域污染负荷的时空变化特征, 并揭示了流域时空变化特征与社会经济-土地利用之间的关系, 与之前的研究相比, 本研究在对点源、面源污染污染负荷都进行估算的基础上对流域污染负荷进行多元统计分析, 规避了水质监测数据难获取的问题, 并将多元统计分析结果赋予地理空间属性, 对研究区域的污染状况分析更加全面、直观.此外, 本研究估算污染负荷时的排污系数是基于研究区特殊的自然地理、社会经济状况, 其他类似流域估算污染负荷时可根据自身特点在本研究的基础上进行修正, 从而使估算结果更加接近实际.本研究采用的流域污染负荷时空分析方法能够更好地帮助决策者识别污染严重区域及污染负荷变化趋势, 为水环境污染管理提供参考, 同时也可为其他类似流域提供一个污染负荷时空分析思路.

5 结论(Conclusions)

1) 沱江流域4种污染物的污染负荷呈整体增加、局部减少的空间变化特征, 增长量和增长率最多的是COD, 分别为87556.54 t、24.16%.上游的龙泉驿区COD、NH₃-N的增加量最多, 分别为7988.15、469.83 t, TN和TP增加量最多的是下游的龙马潭区和泸县, 分别为1556.08、181.93 t.

2) 2007和2017年, 非农业人口与COD、NH₃-N、TN和TP污染负荷的相关系数分别由0.460、0.218、0.226和0.184变为0.953、0.938、0.881和0.871, 非农业人口对污染负荷贡献变大, 有大面积农业用地和建设用地分布的地区污染严重.

3) 4种污染负荷的污染源聚类结果均发生变化, 这与其他污染源类型向城镇污染源类型转变相关, 发生变化的区县多数集中在有农业用地和建设用地大面积分布的区域, 且大部分区县城镇化率或二、三产业占生产总值比的变化量较大, 如中游的乐至县、雁江区、东兴区等, 人口结构和经济结构的改变使污染结构发生变化, 污染源聚类结果随之发生变化, 聚类结果未发生变化的区县有大面积林地和草地分布.本研究对流域污染负荷的时空变化特征进行分析, 可为水环境污染治理和流域水功能规划管理提供有价值的参考信息, 也为建立和完善严格监管所有污染物排放的环境保护管理制度提供理论依据和数据支撑.