基于社会-经济因子修正的沱江流域农业面源总磷污染负荷时空演变研究

姚婧，杨丽君，肖宇婷，樊敏*，谌书，刘云峰，王华伟，陈雯，邓粤，王美玲

（1.西南科技大学环境与资源学院，四川 绵阳 621010；2.成都市沱江流域投资发展集团有限公司，成都 611741）

农业面源污染是指在农村生活和农业生产活动中溶解的或固体的污染物（如生活垃圾、农田中的化肥和畜禽粪便等），在降雨和径流的冲刷下，从非特定的地域进入受纳水体而引起的污染，其具有形成过程复杂、随机性大、分布范围广以及潜伏期长等特点。随着工业和城市点源污染逐步得到治理和控制，农业面源污染对水环境的影响日益凸显，我国《第二次全国污染源普查公报》显示，2017年农业源总磷（TP）排放量为2.120×10t，占全国总TP 排放量的67.22%，农业面源TP 污染已成为导致流域水环境恶化的主要原因之一，同时也是制约农业经济可持续发展和农村水环境质量改善的重要原因。因此，农业面源污染负荷未来变化趋势、时空分布特征分析及其空间相关性探索研究对流域水环境质量改善具有重要意义，可以为环保部门制定相关政策提供重要科学依据。

目前，国内外学者针对面源污染开展了一系列研究，主要围绕污染负荷估算方法及源解析、时空分布特征分析、评价方法、污染类型划分及防治等内容。污染负荷估算方法从总体上可分为集总式水文模型和分布式模型。集总式水文模型主要包括水质水量法、输出系数法和平均浓度法等，其中应用较为广泛的输出系数模型计算原理简单，所需参数少，且具有一定的精度，适合于估算大中尺度流域面源污染的年负荷量；分布式模型主要包括ANSWERS、BASINS、SWAT和AnnAGNPS等，此类模型所需参数多，建模要求高，同时模拟精度也较高。在估算出污染负荷的基础上进行时空分布研究，有助于识别出研究区污染高负荷区域，从而可制定针对性治理政策。目前，污染负荷时空分布研究方法主要有多元统计分析法、GIS 技术、Kernel 密度估计法等。综合评价污染源对水环境的影响是进行污染物总量控制的重要依据，污染源评价方法大致可分为单一评价法、综合评价法和经济技术评价法3 类，其中，应用较广泛的是综合评价法中的等标污染负荷法，该方法简单、可操作性强，且具有较好的综合性。

沱江是长江上游的一级支流，其水质对整个长江流域水环境有着至关重要的作用。近年来，流域磷矿开采活动的不断开展、畜禽养殖规模的扩大以及农田农药化肥用量的增加，导致流域TP 污染问题非常突出。2016 年，沱江TP 平均浓度为0.274 mg·L，在长江一级支流中最高。对此，四川省开展了一系列磷污染排放攻坚工作，但成效有限，TP 污染问题依旧严峻。综合国内外和沱江流域污染负荷的研究现状，发现排污系数法在流域面源污染负荷估算中广泛应用，但该方法忽略了社会-经济因子时空差异对污染负荷估算结果的影响，并不能反映流域实际污染状况。因此，本文以沱江流域为研究区域，基于历史社会-经济因子在该流域28 个区县和四川省的时空变化特征及其比例关系，采用修正的排污系数法预测各污染源TP污染负荷，分析其时空变化趋势，识别主要污染源及污染区县，为相关部门精准治污提供科学依据，进而促进流域水环境持续健康发展。

1 研究区域概况

沱江位于四川省中部，发源于绵竹市清平镇，流经简阳市、资阳市、资中县、内江市、自贡市和富顺县等，全长627.4 km，流域面积为2.76×10km，其中，四川境内流域面积为2.55×10km，地理坐标为28°49'17″～31°42'03″ N，103°40'57″～105°42'50″ E（图1）。沱江流域是四川省人口密度最大、城市分布最密集、经济社会发展最好、工农业生产最发达的地区，承载了全省30.8%的经济总量和26.2%的人口，流域内工业、城镇生活、畜禽养殖和农业种植等各类污染源众多，粮食和经济作物种植力度大、畜禽养殖规模大、土地利用强度大等因素导致流域TP 超标严重，严重影响流域内人民的生产生活，制约着整个流域的社会经济发展。

图1 沱江流域空间地理图Figure 1 Spatial geographical map of Tuojiang River watershed

2 数据与研究方法

2.1 数据来源

本研究所采用的相关基础统计数据主要源于统计年鉴和网络数据平台，详细信息见表1。

表1 数据来源信息Table 1 Data sources

2.2 农业面源核算基础数据预测

考虑到数据的可获取性，本研究采用GM（1，1）模型分3 部分进行相关社会-经济数据的预测：第一部分是分别将2011—2015年和2016—2017年沱江流域各区县旱地和水田面积、作物产量（稻谷、小麦、玉米、豆类、薯类和油料）作为GM（1，1）模型的校核和验证，预测2021—2025 年变化趋势；第二部分是分别将2011—2017年和2018—2019年沱江流域各区县农村人口数、畜禽年末出栏头数（生猪、牛、羊、家禽）和第一产业GDP 作为GM（1，1）模型的校核和验证，预测2021—2025 年变化趋势；第三部分是分别将2011—2017 年和2018—2019 年四川省农、牧业及第一产业GDP 作为GM（1，1）模型的校核和验证，预测2021—2025年变化趋势。

2.3 各区县TP污染负荷估算

本研究采用排污系数法，综合相关文献和沱江流域实际地理环境状况，确定各污染源TP 污染负荷产污系数，具体计算公式及参数值见表2。

表2 各污染源TP污染负荷计算所需参数Table 2 The required parameters for calculating TP pollution loads from diverse pollution sources

2.4 TP污染负荷预测修正

人口、农业GDP 和牧业GDP 分别是驱动农村生活类（农村生活污水和农村生活垃圾）、农业种植类（农田径流和农田固废）和畜禽养殖类污染负荷变化的主要因素。传统排污系数法估算污染负荷时忽略了这些因素，其估算结果可能会低估（高估）部分农业人口较多（少）或农牧业发达（落后）区县的污染负荷。因此，本文首先引入人口因子对农村生活污水和农村生活垃圾污染源产污系数进行修正，其次引入农业GDP 对农田径流和农田固废污染源产污系数进行修正，最后引入牧业GDP 对畜禽养殖污染源产污系数进行修正，进而得到修正的各污染源TP 污染负荷。沱江流域各区县人口数据可以直接获取，因此直接基于各区县人口数据计算各区县的人口修正系数α，其详细计算过程见图2和表3中的公式（6）。鉴于该流域28 个区县2011—2017 年农、牧业GDP 统计数据的可获取性及沱江流域涉及行政区域对四川省社会经济发展的贡献性（沱江流域以仅占全省3.5%的水资源支撑了全省30.8%的经济总量和26.2%的人口），首先，基于四川省2011—2017 年第一产业和农、牧业GDP 数据，采用GM（1，1）模型预测其2021—2025 年变化趋势，计算其农、牧业GDP 比例系数ω和Ψ，其次将2021—2025 年沱江流域各区县第一产业GDP 值乘以相应的农、牧业比例系数ω和Ψ，得到各区县2021—2025年农、牧业GDP 值，最后计算农业修正系数β和牧业修正系数γ，详细计算过程见图2 和表3中的公式（7）～公式（12）。

表3 各污染源污染负荷修正系数Table 3 The modified coefficients of pollution loads from various pollution sources

图2 TP污染负荷预测修正流程图Figure 2 Flow chart of prediction and correction of TP pollution loads

2.5 局部空间相关性分析

局部空间自相关可反映一个区域单元上的某种地理现象或某一属性值与邻近区域单元上同一现象或属性值的相关程度，该方法能够描述空间对象分布中所存在的局部特征差异，常用Local Moran’s I 衡量，具体计算公式为：

3 结果与分析

3.1 GM（1，1）模型及修正的排污系数法验证

本文将2016年和2017年沱江流域各区县耕地面积、作物产量实际数据与其预测数据分别进行比较（表4），同时将2018 年和2019 年沱江流域各区县农村人口、畜禽年末出栏头数和第一产业GDP 及四川省农、牧业和第一产业GDP 实际数据与其预测数据分别进行比较（表5），其最大误差低于10%，表明此方法具有较高的计算精度，计算结果可信，可以将该模型应用到2021—2025 年社会-经济数据预测和TP污染负荷修正。

表4 耕地面积-作物产量的GM（1，1）模型验证Table 4 GM（1，1）model validation of cultivated land area and crop yield

表5 其他社会-经济数据的GM（1，1）模型验证Table 5 GM（1，1）model validation of other socio-economic data

此外，本文比较了基于传统排污系数法和修正排污系数法的TP 污染负荷估算结果，进一步验证了本研究提出的基于社会-经济因子修正TP 污染负荷计算方法的可行性。本研究首先基于2011—2016 年沱江流域各区县农业人口、耕地面积、作物产量和畜禽年末出栏头数等数据，采用GM（1，1）预测其2017 年变化趋势，其次分别采用传统排污系数法和修正排污系数法计算不同污染源TP 污染负荷，并将其与前人研究结果进行对比分析，结果见表6。各污染源TP污染负荷预测值与前人研究结果相对误差保持在5%以内，农村生活污水、农村生活垃圾、农田径流和农田固废TP污染负荷修正值与前人研究结果相对误差保持在10%左右，修正结果在可接受范围内。畜禽养殖TP污染负荷修正值与前人研究结果相对误差保持在15%左右，修正结果偏高。统计年鉴中未反映出各区县畜禽死亡淘汰数，这部分产出粪尿未计入粪尿总量中，导致估算结果偏小，因此，本研究修正后的畜禽养殖TP污染负荷可大体反映该流域未来污染状况。同时，污染负荷估算参数应是动态数据，而传统排污系数法估算污染负荷时取静态值，其估算值未能如实反映各地区污染物排放量的时空差异。修正的TP 污染负荷与未修正TP 污染负荷的误差是由基于四川省农、牧业GDP 与第一产业GDP 比例系数估算得到的沱江流域各区县农、牧业GDP 数据与其实际值之间的差异导致的，鉴于沱江流域涉及行政区域对四川省社会经济发展的贡献性，四川省农、牧业GDP 数据仍可用于沱江流域各区县农、牧业GDP 的估算。综上分析，修正排污系数法更为合理，其估算值可客观地反映流域未来污染状况。

表6 排污系数法与修正的排污系数法验证Table 6 The verification of corrected and uncorrected pollution discharge coefficient method

3.2 TP污染负荷时间变化分析

如图3a 所示，2021—2025 年，沱江流域总TP 污染负荷由12 194.92 t 增加至12 461.26 t，将增加266.34 t，增幅为2.18%；如图3b～图3d所示，农村生活污水、农村生活垃圾和畜禽养殖的TP 污染负荷均呈逐年递减趋势，其将分别减少94.24、77.90 t 和86.52 t，降幅分别为16.92%、16.92%和1.10%；如图3e 和图3f 所示，农田径流和农田固废的TP 污染负荷均呈逐年递增趋势，将分别增加74.52 t和450.49 t，增幅分别为4.82%和19.64%。各污染源TP污染负荷的绝对变化量将表现为农田固废＞农村生活污水=农村生活垃圾＞农田径流＞畜禽养殖。2021—2025 年沱江流域各污染源TP 污染负荷贡献率如图4 所示，畜禽养殖贡献率将最高，在60%以上，将成为沱江流域主要污染源，农村生活垃圾贡献率将最低，低于4%。各污染源TP污染负荷贡献率将表现为畜禽养殖＞农田固废＞农田径流＞农村生活污水＞农村生活垃圾。

图3 2021—2025年各污染源TP污染负荷时序变化特征Figure 3 Temporal variation of TP pollution loads from diverse pollution sources from 2021 to 2025

图4 各污染源TP污染负荷贡献率Figure 4 Contribution rate of TP pollution loads from diverse pollution sources to summed TP pollution loads from all pollution sources across whole Tuojiang River watershed

3.3 TP污染负荷空间演变分析

本研究选取2017 年和2025 年各污染源TP 污染负荷数据，探究空间分布演变规律（图5）。农村生活污水和农村生活垃圾TP污染负荷空间分布特征较为相似（图5a和图5b）。2017年，农村生活污水、农村生活垃圾TP 污染负荷较高的区县均位于中游的仁寿县、资中县和安岳县以及下游的泸县，2025 年将位于上游的新都区和龙泉驿区。2017 年，农田径流TP 污染负荷较高的区县位于中游的安岳县，2025 年将位于中游的安岳县和下游的荣县和江阳区（图5c）。2017 年中游安岳县的农田固废TP 污染负荷较高，2025 年中游的安岳县和东兴区以及下游的泸县和荣县的农田固废TP 污染负荷将较高（图5d）。2017 年，畜禽养殖TP污染负荷较高的区县位于中游的安岳县和下游的泸县，2025 年将位于中游的仁寿县和安岳县以及下游的荣县、泸县和江阳区（图5e）。

图5 2017年和2025年各污染源TP污染负荷空间分布Figure 5 Spatial distribution of TP pollution loads from diverse pollution sources in 2017 and 2025

3.4 TP污染负荷空间集聚效应分析

为进一步揭示社会-经济因子对TP 污染负荷空间分异格局的影响，本研究将修正前后各污染源TP污染负荷局部自相关集聚图进行对比，结果如图6 所示。农村生活污水和农村生活垃圾TP污染负荷修正后较高、较低污染风险区县数量将增加，高污染风险区县将消失（图6a 和图6b），修正后较高污染风险区为新都区，低污染风险区为青白江区，较低污染风险区为泸县。农田径流TP 污染负荷修正后较高、较低污染风险区县数量均将减少（图6c），修正后较高污染风险区为资中县，较低污染风险区为新都区和广汉市。农田固废TP污染负荷修正后较高污染风险区县数量将减少，较低污染风险区县数量将增加（图6d），修正后较高污染风险区为雁江区、资中县和安岳县，较低污染风险区为什邡市、彭州市、广汉市、新都区、青白江区和旌阳区。畜禽养殖TP污染负荷修正后较高、较低污染风险区县均将减少（图6e），修正后低污染风险区为龙马潭区，较低污染风险区为自流井区。

图6 各污染源TP污染负荷修正前后局部自相关集聚对比图Figure 6 Contrast of local autocorrelation cluster diagram of uncorrected and corrected TP pollution loads from diverse pollution sources

4 讨论

4.1 沱江流域TP污染负荷时空变化特征

2021—2025 年，沱江流域TP 污染负荷呈逐年增加趋势，其中，农村生活类和畜禽养殖类TP污染负荷呈逐年减少趋势，农业种植类TP 污染负荷呈逐年增加趋势（图3）。农村生活类TP 污染负荷的减少与人口的减少、生活污水处理设施的完善及生活垃圾的集中处理有关。沱江流域畜禽养殖业逐渐由散养向专业化养殖转变，产污系数较大的畜禽数量将逐年减少，畜禽粪便的处理效率不断提升，畜禽排泄物对环境的污染日趋减缓。农业种植源污染负荷的增加主要与农药化肥的不合理施用有关。此外，畜禽养殖、农田径流和农田固废将是主要污染源。胡芸芸等也发现沱江流域农业面源污染主要来自于畜禽养殖和种植类污染，杜明等研究发现沱江流域TP超标严重，畜禽养殖、农田径流和农村生活是其主要污染源。长期粗放的农业生产方式、畜禽养殖废水废物的肆意排放和农村生活废水的低处理率促使沱江流域发展为生产型污染，即种植业生产与畜禽污染。未来水环境治理应侧重于科学有效地使用化肥农药，减少不合理的用量，提高农业生产的效率和效益，同时推动畜禽粪便、农作物秸秆、农膜资源化利用，加强种养结合和农牧循环。农村生活污水和农村生活垃圾TP 污染负荷的空间差异主要与各地区人口有关。城市化进程的不断推进使得隶属于成都市的新都区和龙泉驿区以及彭州市人口增加，生活用水量增大，TP 污染负荷增加（图5a 和图5b）。农田径流和农田固废的TP污染负荷主要与化肥和农药的不合理使用以及农业生产固体废弃物的低效率处理有关。减征或免征农业税政策的实施使得农户种粮积极性高涨，粮食及经济作物产量增加，同时农药化肥使用量和农作物固体废弃物增加，使得TP污染负荷将增加，这种变化在沱江中游地区尤为明显（图5c 和图5d）。畜禽养殖TP污染负荷的差异性分布与其各地区畜禽种类和数量密切相关（图5e）。

4.2 各污染源TP污染负荷空间集聚特征

农村生活污水和农村生活垃圾TP污染负荷修正后，较高、较低污染风险区县数量将增加，高污染风险区县将消失（图6a 和图6b），这与修正的排污系数法基于每个区县人口与整个流域人口的比例重新分配污染负荷有关。修正后农村人口较多区县的TP污染负荷将增加，人口较少的区县TP 污染负荷将减少。农田径流TP污染负荷修正后，较高、较低污染风险区县数量均将减少，农田固废TP 污染负荷较高污染风险区县数量将减少，较低污染风险区县数量将增加（图6c和图6d），这与修正的排污系数法基于各区县农业GDP 与整个流域总农业GDP 的比例重新分配各区县的这两种污染源TP污染负荷有关，修正后农业GDP值较高区县的TP 污染负荷将增加，较少区县的TP 污染负荷将减少。畜禽养殖TP污染负荷修正后较高、较低污染风险区县均将减少（图6e），这与修正的排污系数法基于各区县牧业GDP 与整个流域牧业GDP 的比例重新分配各区县TP污染负荷有关，修正后牧业GDP较高区县的TP 污染负荷将增加，较少的TP 污染负荷将减少。

基于以上结论，就不同污染风险等级区县而言，较高污染风险区县及其周围区县的TP污染负荷均处于较高水平，因此，未来地方政府应采取一系列强有力的污染防治措施以减缓当地TP 污染，遏制高污染负荷区县的扩散，例如完善农村地区地下水管网设施，建立农村生活垃圾收集、运输和处理系统；优化农业种植模式，减少农药化肥的施用量；提高畜禽养殖废弃物的处理效率等。高污染风险区县TP污染负荷高，而其周围区县低，因此未来高污染风险区县可以经济补偿的方式将部分污染物转移到邻近污染负荷较低的区县。低污染风险区县TP 污染负荷低，易于受到邻近TP 污染负荷较高区县溢出效应的影响，因此今后决策者应加强低污染风险区县污染治理，密切关注邻近区县污染负荷发展动向，防止其污染负荷反向溢流。较低污染风险区县及其周围区县TP污染负荷较低，建议该地区实施“预测为主，治理为辅”的治理政策，使其污染负荷稳于较低水平。

综上所述，本研究提出的基于社会-经济因子修正的流域污染负荷预测框架，结构简单、计算方便，有效降低了各污染源产污系数时空变化对预测结果的影响，凸显出不同行政单元之间TP 污染负荷的空间异质性，可以推广到其他流域类似的污染负荷时空变化趋势研究，拓展了流域污染负荷预测的研究范畴。但本研究中排污系数法的参数取值部分是通过查阅文献获得，在今后的研究中，应尽量结合现场调研与试验，获取更为准确的参数值，以提高预测精度，从而为水环境污染的预防和管理提供正确方向。

5 结论

（1）沱江流域TP 污染负荷总体呈逐年递增的趋势，由2021年的12 194.92 t增加至2025年的12 461.26 t，将增加266.34 t，增幅为2.18%。农村生活污水、农村生活垃圾和畜禽养殖的TP污染负荷呈逐年递减的趋势，将分别减少94.24、77.90 t 和86.52 t，降幅分别为16.92%、16.92%和1.10%；农田径流和农田固废的TP 污染负荷呈逐年递增的趋势，将分别增加74.52 t和450.49 t，增幅为4.82%和19.64%。各污染源TP 污染负荷贡献率将表现为：畜禽养殖＞农田固废＞农田径流＞农村生活污水＞农村生活垃圾。

（2）2017 年，农村生活污水、农村生活垃圾TP 污染负荷较高的区县均位于中游的仁寿县、资中县和安岳县以及下游的泸县，2025 年，较高的区县将位于上游的新都区和龙泉驿区；2017年，农田径流TP污染负荷较高的区县位于中游的安岳县，2025 年，较高的区县将位于中游的安岳县和下游的荣县和江阳区；2017年中游安岳县的农田固废TP 污染负荷较高，2025 年中游的安岳县和东兴区以及下游的泸县和荣县的农田固废TP 污染负荷将较高；2017 年，畜禽养殖TP 污染负荷较高的区县位于中游的安岳县和下游的泸县，2025 年较高的区县将位于中游的仁寿县和安岳县以及下游的荣县、泸县和江阳区。

（3）2025 年修正前后各污染源TP 污染负荷集聚区域分布差异显著。农村生活污水和农村生活垃圾TP 污染负荷修正后较高、较低污染风险区县数量将增加，高污染风险区县将消失；农田径流TP污染负荷修正后较高、较低污染风险区县数量均将减少；农田固废TP污染负荷修正后较高污染风险区县数量将减少，较低污染风险区县数量将增加；畜禽养殖TP污染负荷修正后较高、较低污染风险区县均将减少。