官厅水库上游农业面源污染时空分异研究

郭慧

(清鉴环境与文化工作室，北京 102288)

农业面源污染是相对工业点源污染而言的水体污染类型。随着我国针对工业点源污染治理、工业供给侧改革的不断深入，农业面源污染对我国水环境造成的负面影响不断凸显。我国农业生产中化肥、农药等施用方式粗放、亩均投入量较高等问题，导致其对环境造成较大的负外部性影响[1]。农田种植、农村畜禽养殖，以及农村及城乡接合部的生活排污等，被认为是造成水体氮、磷富营养化的主要原因[2-3]。农业面源污染又因其分布与扩散方式分散，不确定性强等特点[4]，导致其治理难度大。这已成为未来我国水环境治理面临的焦点和棘手问题。

我国已有的针对流域农业面源污染时空分异的相关研究，多是基于数据分析与地理信息系统(GIS)技术应用，依靠研究区土地利用类型、农药使用量、化肥施用量、畜禽养殖量等，结合研究区其他空间信息数据进行的面源污染负荷研究。徐丽萍等(2011)以化肥投入浓度、单位面积农药施用量等衡量指标，在GIS技术的支持下，研究了新疆地区规模化农业活动中主要污染物在空间分布上的差异[5]。韩书成等(2018)建立了化肥、农药及畜牧养殖污染程度评价指标模型，对广东省湛江市典型年份的农业面源污染进行了空间化与数字化的处理，进而研究其空间分布格局[6]。也有针对小流域面源污染负荷，结合GIS与SWAT、HSPF等水文模型进行的面源污染负荷模拟研究[7-8]。以上研究仅针对农业生产与水环境关系的简单二维分析研究，缺少与宏观产业结构发展的多元统计分析。本文在研究方法上的创新之处在于，在应用GIS技术对农业生产与水质关系的时空分异进行研究的基础上，根据芬兰环境经济学家塔皮奥(Tapio)的弹性脱钩指数理论，研究区域环境压力，尤其水体水质情况与农业经济增长间的关系，以期探索一条环境与产业经济协调可持续发展的道路。

1 数据与方法

1.1 研究区概况

官厅水库作为中华人民共和国成立后建设的首座大型水库，于1954年5月竣工，设计库容41.6亿立方米，控制流域面积达43402m2。作为北京曾经的两大地表水源地之一，官厅水库为周边地区的民生与经济发展做出了巨大贡献。然而由于改革开放初期的工业与农业污染日趋严重，尤其是有机氯农药的大量使用以及重金属的超量排放对水质安全造成的巨大威胁，官厅水库于1997年退出北京市饮用水源地体系。其后的二十余年间，众多学者对官厅水库及其上游流域的水质进行长期研究，发现其水体的富营养化问题仍旧严重。氨氮、总磷、化学需氧量等水质指标仍旧在Ⅲ-Ⅳ类标准间徘徊，其结论均指向农业面源污染未能得到有效控制[9-12]。

官厅水库库区位于北京市延庆区与河北省张家口市怀来县交界处的永定河上游地区。官厅水库上游包括库区西部的洋河及桑干河两大支流。洋河由东洋河、西洋河、清水河等河流汇合形成，分别流经怀安、万全、宣化、怀来等县区；桑干河及其支流分别流经阳原、蔚县、宣化、涿鹿、怀来等县区。两条河流在怀来县汇合为永定河，继续向东行汇入官厅水库。另外库区东部的妫水河自北京市延庆区发源，并在延庆汇入官厅水库。

本研究的主要范围为借助ArcMap10.2版软件对数字高程模型(DEM)数据的分析，进行水库上游流域边界提取，并选取在河北省张家口市内的区域，具体包括张家口市的尚义县、怀安县、涿鹿县、怀来县、阳原县、蔚县、崇礼县(现为崇礼区)、宣化区、万全区、下花园区、张家口市辖区，共计7县4区。由于研究时期(2004年—2017年)，部分县区的行政区划有所调整，本研究以当时的县区区划为准，研究区内流域分布及行政区划分布，参见图1。

图1 研究区内流域分布及行政区划分布示意图

1.2 数据来源与初步处理

1.2.1 空间数据

进行流域边界提取的空间数据，来自“地理空间数据云”平台提供的免费数字高程模型(DEM)30米分辨率数据(1)数据来源：“地理空间数据云”平台，中国科学院计算机网络信息中心，http：//www.gscloud.cn。。通过AcrGIS10.2软件对DEM数据分别进行洼地填充、流向分析、流量分析、河网及河流分级分析，提取得到官厅水库流域的河网，对河网进行流域边界提取，与县区级行政区划数据结合，裁剪得到研究区域。因本研究侧重以县区行政区划进行的农业面源污染情况与农业经济比较，因此不再对研究区进行子流域划分。

1.2.2 研究区农业面源污染排放量的估算

根据研究区域的划定，本研究选取上述8县2区的与农业面源污染相关的化肥(包括氮肥、磷肥、复合肥)折纯施用量、农药施用量、畜禽养殖量等指标数据，以上原始数据均来自2005—2018年《张家口经济年鉴》(统计数据范围为2004—2017年)(2)数据来源：《张家口经济年鉴》，中国经济社会大数据研究平台，http：//data.cnki.net/yearbook/Single/N2019070141。。

针对农业面源污染排放量的估算，已有学者采取不同推算方法的研究，如韩书成等(2018)、徐丽萍等(2011)以单位面积的化肥、农药投入量即投入密度作为指标进行衡量[5，6]，亦有如吴义根等(2017)以统计数据结合产污强度系数与排污系数进行估算[13]。本研究参照后者的估算方法，结合研究区农业生产实际特征进行统计数据初步处理。首先，将各类农业面源污染源根据其污染方式、产污单元及排放清单进行分类统计(表1)；其次，分别根据《第一次全国污染源普查—农业污染源系数手册》(以下简称《手册》)中的流失系数指标选取原则，以2010年第一次全国污染源普查结果中，农业源排放占地表水体污染总负荷最高比例的三大指标——化学需氧量(COD)、总氮(TN)、总磷(TP)为主要考察对象；再次，结合《手册》中的产排污系数与流失系数，建立各污染类别的面源污染排放量与产污单元的关系式(式1-3)。

表1 面源污染的产污单元与排放清单分类

化肥类别污染物排放量：

Ei=∑iCFjLi

(1)

畜禽养殖类别污染物排放量：

Ei=∑iNjDjFPi

(2)

式中：Ei为农业面源污染物i的排放量，CFj为化肥类别中产污单元j的统计数，Li为污染物i的流失系数；Nj为产污单元j的饲养数量，Dj为产污单元j的饲养天数(其中猪、牛、禽类的饲养天数分别为180天，365天，55天)，FPi为产生污染物i的产污系数。通过上述公式，分别核算出7县4区各污染物的排放量。

农药的施用虽然不会对水质相应的CODCr、TN、TP等指标产生直接影响，但农药危害性对环境及生物体均造成不利影响。因此，本研究也将对农药类别污染物排放量进行时空分析。

农药类别污染物排放量：

式中：Pi为i地区农业面源污染中农药污染排放量，EPi为农药施用量的统计数，Li为《第一次全国污染源普查—农药流失系数手册》中，9大农药对应地区及土地类型的平均流失系数，计算结果为0.0033%。

1.2.3 研究区水质对照数据

本研究涉及的环境时序数据为水质数据，以此与研究区内第一产业生产总值进行弹性脱钩指数的分析。水质数据来自“青悦开放环境数据中心”抓取的生态环境部“全国主要流域重点断面水质自动检测周报”中公布的洋河—官厅水库入口断面，位于张家口八号桥点位2011年第25周至2017年第53周的数据(3)数据来源：全国主要流域重点断面水质自动检测周报，青悦开放环境数据中心，http：//data.epmap.org/eia/water。。考虑季节及丰水期与枯水期等因素，分别选取每年第25周至第52(或53)周的数据做进一步处理，选取水质指标为化学需氧量(CODMn，mg/L)，以及氨氮(NH3-N，mg/L)。选取原则是化学需氧量为农业源排放占地表水污染负荷最高比例的污染指标，生态环境部重点断面水质自动检测周报公布的水质指标数据中未有总氮与总磷指标数据(4)生态环境部公布的全国主要流域重点断面水质自动检测周报中的水质指标分别为：pH、DO、CODMn及NH3-N四个指标。，因此选择与其最接近的指标——氨氮作为研究对象。

1.2.4 经济统计数据

经济统计数据为研究各县区按当年价计算的第一产业生产总值，原始数据与水质时序数据的时间区间一致，来自2012—2018年的《张家口经济年鉴》(统计数据范围为2011—2017年)(5)数据来源：《张家口经济年鉴》，中国经济社会大数据研究平台，http：//data.cnki.net/yearbook/Single/N2019070141。。因考虑价格变动等因素对统计数值的影响，将原始数据换算为以2000年为基期的不变价用于后续计算。

1.3 研究方法

1.3.1 基于GIS的面源污染空间聚类分析

以各县区为单位，应用Excel软件对上述收集到的相关研究数据进行初步计算处理。以2005年、2010年及2015年作为研究代表年，应用SPSS数据分析软件分别对研究区内各县区的化肥、畜禽养殖及农药污染类别进行污染物排放量的K-means快速样本聚类分析。根据聚类分析结果进行上述污染类别的分级，并将处理后的数据导入ArcMap10.2软件中，在研究区的各县区行政区划内，根据分级结果进行展布，得到2005年、2010年、2015年研究区内面源污染空间格局图。

1.3.2 弹性脱钩指数的应用

在早期社会发展研究中，人们认为环境资源消耗及废物排放量是与经济总量“挂钩”的。20世纪末德国Wuppertal研究所的学者提出了“脱钩目标”，即提高国家的资源利用效率以实现资源消耗与经济增长的脱钩。为衡量脱钩目标的实现程度，以及识别经济发展与环境资源间关系的动态情况，于是由经济合作与发展组织(OECD)提出了脱钩因子模型[15]，随后Tapio提出了弹性脱钩指数模型。Tapio的弹性脱钩指数模型，采取时间维度的弹性分析方法相对客观的反映变量间的脱钩关系，克服了OCED脱钩因子模型的基期选择问题[16]。本文即采用Tapio的弹性脱钩指数模型，针对研究区内的环境与农业经济进行时空分异研究。

结合已有针对农业化学投入、废物排放、水体环境压力与经济脱钩关系的研究[14-18]，基于Tapio的弹性脱钩指数模型，建立研究区内国家重点断面——张家口八号桥化学需氧量及氨氮两个水质指标与农业产值间的时间序列模型(式4)：

式中：Et为第t年的脱钩状态；ΔPt与ΔGt分别为第t年各水质指标与农业经济产值的增长率；Pt与Gt分别为第t年各水质指标与农业经济产值；Pt-1与Gt-1为第t-1年各水质指标与农业经济产值。

Tapio(2005)在其道路交通与欧盟GDP增长脱钩程度的研究中，提出模型结果可以划分为8类状态分类[19]。本研究参考Tapio的模型应用与国内相关研究[15]，采用相似的判别与分析方法，如表2所示。

表2 农业经济与面源污染的脱钩状态分类

2 结果与讨论

2.1 农业面源污染时序特征分析

运用上述研究方法进行数据分析：对研究区内面源污染中的三个指标——总氮(TN)、总磷(TP)、化学需氧量(COD)，以及农药的排放总量做时序变化趋势分析，并将TN与TP两个污染指标对应的污染源(化肥投放及畜禽养殖)所占比例组合至序列中。对于统计数据计算得到图2至图5。

总氮排放与总磷排放在本研究中涉及农业化肥投入与畜禽养殖两部分，因这部分的年度变化趋势影响，导致TN与TP在时间分布上呈现相似的变化特征。自2004—2017年的14年中，TN与TP在研究区的排放总量均呈现缓慢上升趋势。其中2007—2010年，TN、TP呈现增势减缓且略有下降的情况，尤其化肥投入比例下降。2008年，张家口市开始推广测土配方施肥项目。按照原农业部《测土配方施肥技术规范》进行项目推进。以怀来县为代表的葡萄种植县区，开展“葡萄配餐工程”(针对葡萄种植的测土配方施肥技术)。化肥投入成本降低，生态环境压力得到缓解。据报道化肥施用量每亩降低约75公斤[20]。然而，2011年，TN、TP均呈现断崖式上升情况，且化肥投入比例急剧升高，这与吴义根(2017)等学者的研究结果一致，与当年我国粮食产量大幅增长有密切联系[13]。2011年，原环境保护部印发《关于进一步加强农村环境保护工作的意见》。同年底，原农业部出台《关于加强“十二五”中央农村环境保护专项资金管理的指导意见》，2012年即投入55亿元资金用于农村环保专项治理。从分析结果看出，2012年及其后大部分年份，研究区内TN与TP的排放总量大幅下降，且维持相对平缓的走势。

图2 农业面源污染总氮及污染类别比例特征

图3 农业面源污染总磷及污染类别比例特征

图4 农业面源污染化学需氧量排放特征

图5 农业面源污染农药排放特征

化学需氧量的污染总量仅通过畜禽养殖计算得到，因此在这里单独分析其变化特征。COD的变化趋势与TN与TP在2007—2010年的变化趋势相近，都出现增长较缓的态势。而在2011—2016年出现新一轮相对快速增长，直至2017年则大幅度下降。从TN与TP的排放总量趋势图中，也可看出2017年畜禽养殖的贡献占比有明显降低。2017年，张家口市环保局印发《张家口市环境保护局关于加强畜禽养殖污染防治工作的通知》，在畜禽养殖禁养区专项整治及规模化畜禽养殖小区(场)污染防治设施方面，于2017年至其后的几年取得了明显成效。此政策的建立与实施，对畜禽养殖及对流域COD、TN、TP等污染减排方面有显著贡献。

农药污染方面的时序变化特征同样呈现上升趋势。在2014年出现显著增长态势，2015年大幅回落。2014年原农业部开始“农药产品质量专项监督抽查工作”，以及“高毒农药定点经营示范项目”，河北省均作为示范省参与执行。2015年初，原农业部在全国开展“化肥、农药使用量零增长行动”，2015年的农药污染趋势回落应与此有紧密关联。

2.2 农业面源污染空间分布特征分析

农业面源污染的空间分析，可以借助地理空间数据直观地研究区域内面源污染的主要分布特征。叶延琼等(2013)通过对广东省不同年份的化肥、农药投入强度以及畜禽养殖的投放量进行数据分析，并运用GIS对不同污染源的投入强度进行分类及展布[22]。亦有不同学者利用GIS技术通过对农业面源污染的空间展布进行相关空间特征的分析[5-6]。本研究采用这种相对成熟且直观的分析研究方法，针对官厅水库上游研究区内各县区的污染物指标TN、TP与农药排放量，选取2005年、2010年、2015年作为研究代表期，利用SPSS26.0软件进行K-means快速样本聚类分析，将聚类结果导入ArcMap10.2软件，进行空间分布的可视化展布，得出图6、图7和图8各污染指标污染总量的空间分布图。

图6 2005年(a)、2010年(b)、2015年(c)总氮污染总量空间分布

图7 2005年(a)、2010年(b)、2015年(c)总磷污染总量空间分布

图8 2005年(a)、2010年(b)、2015年(c)农药污染总量空间分布

从图5中可以看出，聚类分析得出的第Ⅰ、Ⅱ、Ⅲ类，分别代表污染总量递增的三个分类，其颜色标识逐渐加深。研究区域内涿鹿县的总氮污染位于各区县之首，怀来县于2015年(图5-c)亦进入第Ⅲ类污染排放类别。研究区内Ⅱ、Ⅲ类县区个数也逐渐增加。涿鹿县与怀来县为张家口市历史悠久的葡萄种植大县，其化肥投入量亦远超其他县区，如涿鹿县2015年的氮肥、磷肥及复合肥施用量，分别为6560吨、3918吨及5327吨，与其相对的崇礼县三类化肥的施用量，分别为405吨、80吨及1053吨，差异大者相差40倍。

从图7中可以看出，总磷污染排放的空间分布中，研究区内2/3比例的县区于2010年均已达到第Ⅲ类污染类别(其中聚类Ⅰ中心的总磷流失为0.370吨/年，聚类Ⅱ中心的总磷流失为8.142吨/年，聚类Ⅲ中心的总磷流失为22.420吨/年)，仅张家口市辖区及下花园区位于第Ⅰ类别中。2015年怀安县、阳原县、蔚县与阳原县的总磷流失又有所好转(图7-c)，进入第Ⅱ类别。但各代表年的总磷污染空间分布，仍旧是涿鹿县与怀来县居最高类别。

从图8中可以看出，农药污染总量在研究各县区分布，仍旧为涿鹿县与怀来县居首，且随着时间发展，怀安县、阳原县及蔚县亦逐渐进入较高类别。

从以上三个污染指标的污染总量空间分布中可以了解到，研究区内张家口市辖区及下花园区的各污染排放总量均处于最低水平，主要与其城镇发展、农地耕地面积逐渐缩小有关；尚义县、崇礼县及万全区的污染总量在11个县区中处于较低水平，主要受气候条件影响，其农作物耕种与畜禽养殖量均较低；怀安县、阳原县、蔚县及宣化区的各污染总量有逐渐增加之势；而涿鹿县与怀来县因其气候适宜多果蔬种植与畜禽养殖，成为官厅水库及其上游流域的水体面源污染的重要来源。

2.3 脱钩状态的时空演变及其分析

将研究区由生态环境部发布的“全国主要流域重点断面水质自动检测周报”中位于张家口八号桥点位的2011年第25周至2017年第53周的数据，分别以化学需氧量(CODMn)及氨氮(NH3-N)两个指标，按照前述的数据处理与脱钩指数计算方法，计算得出农业产业生产总值与水质指标的增长率及弹性脱钩指数，如表3及图9所示。

由图9可以看出，研究区内农业面源污染与农业经济的脱钩状态演变：2011—2015年波动较大，2014—2017年呈现逐渐平稳的态势。整体呈现显著的“M”型走势。其中两个波峰分别出现在2013年及2015年期间，均表现为“扩张负脱钩”状态，说明农业经济缓慢增长，而面源污染大幅增长。而在这两个波峰年度前，即2012年与2014年则表现为“强脱钩”状态，即农业经济增长，面源污染下降。这本是环境状况与经济利好发展的趋势，在第二年却出现严重反弹，说明面源污染并未得到根本性改善，这与上述面源污染时序特征分析结果基本一致。2012年，我国投入55亿元资金用于农村环境保护的专项治理，在2013年可以看到其明显减排作用，但专项治理无法持续地改善农村面源污染问题，专项治理后反倒出现大规模反弹情况。

表3 第一产业生产总值与水质指标增长率的弹性脱钩状态

图9 农业面源污染与农业经济脱钩状态时序演变

2015年后，水质指标COD与NH3-N两个指标的增长率与对应的脱钩指数略有不同。COD与农业经济的关系，在2016年仍呈现扩张负脱钩状态。至2017年，由于农业经济增长率降低等原因，表现为“强负脱钩”状态，即农业经济衰退，COD指标对应的面源污染增长。NH3-N与农业经济的关系，在2016年呈现“增长连接”状态，即农业经济增长，面源污染同步增长。2017年则呈现“衰退脱钩”状态，说明农业经济衰退，面源污染同步衰退。

COD指标对应的面源污染来源，主要是畜禽养殖及农业固体废物等；NH3-N指标对应的面源污染来源，主要是畜禽养殖及农业化肥投入。从农业面源污染时序特征分析中可以看出，畜禽养殖所占比例明显减小，且NH3-N指标与农业经济脱钩状态也表明，畜禽养殖对水体造成的污染明显减小。因此，造成COD指标与农业经济呈现强脱钩状态的原因，可推断是由秸秆等农田固体废弃物引起，即秸秆等资源循环利用程度低导致。

3 结 论

本文通过对官厅水库上游流域面源污染的时序与空间分布特征的分析，并结合农业经济数据所开展的环境污染与经济发展的弹性脱钩状态分析，探讨了官厅水库面源污染贡献的时空分异特征，以及其与农业经济发展的关系，得出如下几点结论：

(1)在时序趋势方面，农业化肥投入呈现持续攀升的态势，研究区域内各县区的氮肥、磷肥及复合肥投入总量，从2004年的57422吨增长至2017年的274942吨，为研究初期的近5倍。而农业农药投入量与畜禽养殖量于2016年逐渐下降，这与农药及畜禽养殖污染防治相关政策的出台与有效实施有密切关系。

(2)在面源污染空间分布格局方面，研究区较为下游的涿鹿县与怀来县，为农业面源污染贡献量最大的地区。如涿鹿县的农业化肥投入量占研究区域内各县区的投入总量比例，从2004年的22.70%增加至2017年的24.37%。研究区内南部几个县区的总氮污染贡献量与农药投入量亦随时间有增长的趋势。张家口市随着城镇化发展，其市辖区、下花园区及宣化区则处于面源污染贡献比例较低的地区。

(3)在农业面源污染与经济脱钩状态方面，主要为相对波动的态势，且在部分强脱钩年度后面源污染情况进一步反弹，出现扩张负脱钩状态，甚至呈现出农业经济下行时期，面源污染仍旧加剧的局面。这与农业面源污染防治政策实施的有效性与持续性相关，对面源污染来源应做充分的研究，并针对主要污染源制定可行的防治政策。

本研究过程所采用的数据，大部分来自历年《张家口经济年鉴》，但经济年鉴中的统计数据，在统计过程中存在一定的局限性，无法切实反映污染源的排污情况。另外在研究数据处理方面，无法将所有可能的面源污染指标均纳入统计分析，在未来的研究中，还需将更多可能的污染类别与因子纳入分析体系中。