化肥既是重要的农业生产要素, 又是导致农业面源污染的主要因素之一。中国化肥施用强度较大, 利用效率总体较低^[1], 对粮食增产的边际贡献率不断下降^[2], 单位耕地面积化肥施用量远超国际安全使用水平^[3], 严重威胁到中国农业经济的可持续发展^[4-5]。过量的化肥投入不仅增加了农业生产的经济成本, 还带来了严重的生态环境负担^[6]。探讨农业发展与化肥面源污染之间时空演化关系, 可以更好地推动化肥零增长行动, 促进中国农业高质量绿色发展。

环境库兹涅茨曲线(Environment Kuznets Curve, EKC)是经济增长与环境关系实证研究的重要理论。Grossman等^[7]提出环境与经济之间可能存在类似库兹涅茨曲线的“倒U型”关系。Antle等^[8]和McConnell^[9]较早探讨了农业面源污染与经济发展间的关系, 认为农业面源污染也有EKC的变化特征。相关研究表明, 农业面源污染EKC可能呈现为“倒U型”^[10-12]、“U型”^[13-15]、“N型”或“倒N型”^[16-18]以及线性增长^[19-20]等多种演变趋势。然而, 以上研究存在两方面的不足: 1)用化肥投入量来表征农业面源污染, 忽视了农业面源污染是化肥施用后, 流失的养分进入水体才构成污染的事实^[21], 从而无法反映化肥面源污染实际状况。2)基于同质性假设, 只关注国家或地区的整体状况, 忽略了区域间异质性的影响。各地区经济发展水平、地理环境、资源禀赋、劳动力水平等差异都可能对当地农业面源污染带来不同影响^[22], 使得基于同质性假设的EKC拟合结果不符合区域特质。部分研究按照东、中、西部分组^[23]、省份单独检验^[24]、数理方法选取分类变量^[22]等方式进行分组EKC拟合, 试图克服同质性问题。然而, 东、中、西部的划分方法不能合理反映化肥面源污染的地区特质^[25], 而各省份分别检验的方法则忽视了面源污染在区域间的流动性和分散性^[26], 以及省份间产业流动和劳动力转移等因素导致的空间上的相关性, 进而影响估计结果的准确性。数理方法分组虽有利于运用数据的内生特性优势, 但指标选取往往具有主观性且假定一个区域的环境污染与相邻区域在地理空间上不存在空间相关性^[18], 从而可能使拟合结果与实际情况产生偏差。

本文利用全国1997—2018年的省际面板数据, 首先对各省化肥面源污染排放强度进行核算, 更加准确地反映各省历年化肥面源污染状况。在此基础上, 利用ArcGIS空间统计模型对其进行空间自相关分析和热点分析, 综合时空格局演变分析结果对各省份进行分组。此种分组方式既考虑了不同区域化肥面源污染的实际情况, 又考虑了区域的空间相关性, 更加符合各区域的实际。基于分组, 拟合各组化肥面源污染排放强度与人均农业生产总值之间EKC关系的时间路径, 考察空间异质性条件下各组时间路径的差异性, 最后对时间路径的未来演变趋势进行预测, 试图为化肥减量化和农业面源污染防治工作提供一定参考。

1 研究方法与数据来源 1.1 研究方法 1.1.1 面源污染排放强度核算

当前相关研究多采用单元调查评估方法^[27-28]和清单分析法^[29]进行大尺度区域的化肥面源污染核算, 然而此类方法没有考虑土地类型、自然条件(如气候)等因素的差异对化肥流失产生的异质性影响, 利用化肥施用强度和施用比例确定的化肥流失率也可能导致估计结果不够准确^[30]。因此, 本文在单元调查评估法和清单分析法的基础上, 结合流失系数法进行化肥面源污染的核算。参考《第一次全国污染源普查—农业污染源—肥料流失系数手册》的6大分区, 综合考虑各省气候、地形、土壤等主要影响因素, 将全国化肥面源污染源分为6个区域类别, 各类别的流失系数主要运用“一污普”后发布的《农业污染源肥料流失系数手册》中各区域各模式流失系数的均值, 并通过广泛的文献调研对各类区域的流失系数进行确认和校正, 以使流失系数更符合各区域的实际情况, 保证流失系数的科学性和准确性, 流失系数见表 1。

本文核算的化肥面源污染单元为磷肥、氮肥与复合肥。污染物指标采用总氮(TN)、硝氮(${\rm{NO}}_{\rm{3}}^ - {\rm{ - N}}$)、氨氮(${\rm{NH}}_{\rm{4}}^ + {\rm{ - N}}$)、总磷(TP)和可溶性总磷(DTP)^[31]。运用化肥流失系数法, 核算氮的总排放量和磷的总排放量, 进而核算化肥面源污染排放总量及强度。计算公式为:

$ E = \sum {{E_{ij}} = } \sum {{C_i} \times {\eta _{ij}}} $

(1)

$ I=E/A $

(2)

式中: E为化肥面源污染总排放量; E_ij为污染单元i产生并流失进入水体的第j种污染物的量; C_i是污染单元i的指标统计数(氮肥的统计数为N的量, 磷肥统计数为P₂O₅的量, 复合肥按32%折算P₂O₅的量, 按40%折算TN的量^[32]); η_ij为污染单元i中第j种污染物的流失系数; A为各省农作物播种面积^{[22, 27]}; I(kg·hm-²)为单位面积的污染排放量, 即化肥面源污染排放强度。

1.1.2 空间统计模型

1) 全局空间自相关。化肥面源污染排放强度在全局空间上的相关程度用Global Moran’s I测度, 其含义是空间相邻省份的化肥面源污染排放强度在整体区域内空间相关的总体趋势^[33]。在给定的显著水平下, 若Moran’s I值为负, 则污染排放强度呈显著的空间离散趋势; 若Moran’s I值为正, 则污染排放强度呈显著的空间集聚趋势^[34]。Global Moran’s I的计算公式为:

$I = \frac{{\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{W_{ij}}{\rm{(}}{X_i} - \bar X} {\rm{)(}}{X_j} - \bar X{\rm{)}}} }}{{{S^2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{W_{ij}}} } }}$

(3)

式中: n为省份个数, W_ij为空间权重矩阵(两省相邻为1, 不相邻为0), X_i和X_j为i省和j省的观测值, $\overline X $为观测值的平均值, S²为观测值的方差。

2) 局部空间自相关。局部自相关特征用热点分析的G_i^*指数进行测算, 用于识别区域内的冷点区和热点区。若G_i^*指数显著为正, 则i省周围污染排放强度相对较高, 属于热点区; 若G_i^*指数显著为负, 则i省周围污染排放强度相对较低, 属于冷点区; 若G_i^*指数不显著, 则i省不是热(冷)点区。计算公式为:

${G_i}^*{\rm{(}}d{\rm{)}} = {{\sum\limits_{j = 1}^n {{W_{ij}}} {\rm{(}}d{\rm{)}}{X_j}} \mathord{\left/ {\vphantom {{\sum\limits_{j = 1}^n {{W_{ij}}} {\rm{(}}d{\rm{)}}{X_j}} {\sum\limits_{j = 1}^n {{X_j}} }}} \right. } {\sum\limits_{j = 1}^n {{X_j}} }}$

(4)

式中: W_ij(d)是根据距离规则定义的空间权重。

1.1.3 EKC模型设定

EKC模型通常先设为三次曲线, 若不显著, 则再设为二次曲线, 若二次曲线仍不显著, 则拟合成线性模型^[35]。因此, EKC曲线可能表现为“N型” “倒N型” “U型” “倒U型”或者直线型。为减少异方差, 本文对各指标进行对数化处理, 设定化肥面源污染排放的EKC模型为:

$\begin{gathered} \;\;\;\;\ln Iit = {\beta _{\rm{0}}} + {\beta _{\rm{1}}}{{\rm{(ln}}{G_{it}}{\rm{)}}^{\rm{3}}} + {\beta _{\rm{2}}}{{\rm{(}}\ln{G_{it}}{\rm{)}}^{\rm{2}}} + {\beta _{\rm{3}}}\ln{G_{it}} + {\beta _{\rm{4}}}{\rm{ln}}{P_{it}} + \\ {\beta _{\rm{5}}}\ln{M_{it}} + {\beta _{\rm{6}}}\ln{S_{it}} + {\beta _{\rm{7}}}\ln{R_{it}} + {\beta _{\rm{8}}}\ln{N_{it}} + {\beta _{\rm{9}}}\ln{U_{it}} + {\beta _{{\rm{10}}}}{\rm{ln}}{A_{it}} + {\mu _{it}} \\ \end{gathered} $

(5)

式中: I表示化肥面源污染排放强度; G表示人均农业产值, 作为解释变量, 用来表征农业经济发展水平; G(元∙人^–1)=农业总产值(狭义)/乡村总人口数(狭义农业总产值指种植业和其他农业的主产品和副产品的总产值, 主要包括粮、棉、油料、麻类、烟叶、蔬菜、药材、瓜类和其他农作物的种植, 以及茶园、桑园、果园的生产经营。狭义上的农业是化肥面源污染的主要来源, 因而此处用狭义的农业总产值计算)。其他指标是参考现有相关研究选取的控制变量: P表示作物种植结构, 用粮食播种面积与农作物总播种面积的比值表征, 作物种植结构差异可能致使化肥施用量的差异; M表示农业机械投入强度(kg∙hm-²), 用单位农作物播种面积的农业机械投入强度表征, 机械化水平提高可能会对化肥施用量产生正负两面的影响; S表示产业结构, 用农业产值占GDP的比重表征

$ 农业产值占 \rm{GDP}的比重=农业产值(狭义)/地区国内生产总值×100\% $

(6)

农业产值所占比重大小影响到化肥施用量的大小, 地区产业结构会影响当地的环境质量; R表示农村居民收入结构, 农村居民收入结构用农村居民人均工资性收入与农村居民人均纯收入的比值表征, 农村居民人均工资性收入比重的高低可能影响化肥施用量的大小; N表示农村居民人均纯收入, 农村居民收入的高低会影响化肥的投入; U表示城乡居民收入差距, 用城镇居民人均可支配收入与农村居民人均纯收入之比表征, 收入差距可能是环境退化的重要因素, 均衡的收入分配可能减少污染^[36]; A表示农业产业结构, 用狭义农业总产值占农林牧渔业总产值的比重表征, 以反映种养结构调整对化肥面源污染的影响; μ为随机误差项; i表示省份, t表示年份; 回归模型的估计系数β_n (n=1~10)表示弹性系数。

1.2 数据来源

本文在我国化肥面源污染排放强度核算的基础上, 对其进行时空格局演变分析, 并以时空特征为分类条件, 对空间异质性条件下化肥面源污染EKC的时间路径进行拟合与预测。研究数据为1997—2018年31省(市、自治区, 港澳台数据缺失, 以下均称省)共22年的面板数据, 来源于《中国农村统计年鉴》(1998—2019年)和《中国统计年鉴》(1998—2019年), 并用插值法补齐少量年份因年鉴统计口径变化造成的缺失值。为实现数据的可比性, 农业总产值通过农业产值指数折算为1997年不变价, 农村居民人均纯收入和城镇居民人均可支配收入运用消费者价格指数折算为1997年不变价。

2 化肥面源污染排放强度核算及其时空格局演变 2.1 化肥面源污染排放强度核算

根据公式(1)和(2)对我国1997—2018年各省化肥面源污染排放强度进行核算, 各分区核算结果如图 1所示。我国化肥面源污染排放强度省际差异较大。北京、江苏、上海、浙江、湖北等经济发达省份排放强度较大, 22年间大多数年份在5 kg∙hm-²以上, 超过其他省份的平均强度。青海、西藏、黑龙江等省份排放强度较小。各省22年间化肥面源污染排放强度的走势各异, 表现出不同升降趋势, 表明我国化肥面源污染排放强度时空差异较大。各区域内各省份排放强度差异较小, 因此高(低)强度区域可能呈现一定的集聚特征。所以有必要进一步分析化肥面源污染排放强度的总体和局部时空格局演变。

2.2 省际化肥面源污染排放强度时空格局演变特征 2.2.1 总体时空格局演化特征

运用ArcGIS 10.2软件计算1997—2018年各省域化肥面源污染排放强度的Global Moran’s I指数, 结果如表 2所示。1997—2018年间, 我国省域化肥面源污染排放强度的Moran’s I指数均通过了显著性检验, 且均大于0.307, Z得分均大于4.352, 说明我国化肥污染排放强度呈正向空间自相关性, 化肥面源污染在区域上一直处于较高的集聚状态, 区域间存在较强的相关性和依赖性。Moran’s I指数22年间呈现小幅波动变化, 1997—2006年总体较为平稳, 2006年开始呈现出波动下降的趋势, 说明我国化肥面源污染排放强度空间格局总体较稳定, 空间自相关性呈现出先平稳后下降的特征。空间自相关性降低表明区域与周边地区变化趋势的一致性相对减弱, 应当注意区域间的异质性特征和空间溢出效应的影响。

2.2.2 局部时空格局演化特征

利用ArcGIS 10.2软件的空间统计工具计算化肥面源污染排放强度省域空间的G_i^*指数。依据自然断点法将G_i^*值分为7类:排放强度冷点(99%置信)、排放强度冷点(95%置信)、排放强度冷点(90%置信)、排放强度热点(99%置信)、排放强度热点(95%置信)、排放强度热点(90%置信)和不显著。部分年份空间格局如表 3所示(限于篇幅, 本文只展现时空格局变化较大的1997年、2003年、2007年、2009年、2014年和2018年的空间格局分布情况)。分析可知:

1) 化肥面源污染排放强度空间格局总体较稳定。热点区主要集中分布在华中和华南地区, 其中长江中下游地区尤其显著, 常年处于99%的置信热点。冷点区主要集中在西部地区的青海、甘肃和西藏, 东北的黑龙江也在2003年后长期呈现为冷点状态。化肥面源污染热(冷)点不显著的区域主要集中在我国西南部以及东北和西北的部分省份。

2) 化肥面源污染排放强度空间格局在各阶段呈现出局部变动, 黄淮海地区呈现出热点程度减弱的变化。1997—2002年, 全国化肥面源污染空间格局较为稳定, 各省份热(冷)点的属性未发生大的变动,只有少量省份在置信度上表现出一定的波动。2003年和2004年, 黑龙江在90%置信冷点和不显著间波动, 2005年后则持续呈现为90%的置信冷点。2007年全国空间格局发生较大变化, 京津冀地区首次从热点区消失, 山东的热点显著度也出现下降, 黄淮海部分地区迎来化肥面源污染的拐点。2009年广东首次进入热点区的版图, 至今虽有个别年份呈现为不显著, 但总体呈现为热点状态。2014年后, 部分热点省份发生波动, 如2014年, 广东由90%置信热点变为不显著, 江西由99%置信热点变为95%置信热点, 山东由95%置信热点变为90%置信热点; 2017年, 山东由90%置信热点进而变为不显著, 河南由长期的99%置信热点变为95%置信热点。从全国冷热点时空格局的变化来看, 近年来对于化肥面源污染的治理取得一定成效, 部分区域的热点程度有所下降, 尤其是黄淮海地区的效果较为明显。近年来, 高强度的减量化措施使我国化肥施用量短期减少, 但目前农业生产方式对于化肥的依赖性以及政策执行的波动性还是会使局部地区化肥面源污染状况产生波动。

对31个省22年热(冷)点发生的频次和显著水平进行统计, 据此将31个省份分为热点区、冷点区、非热(冷)点区3类区域, 结果如表 4所示。3类区域化肥面源污染状况的差异体现出了我国化肥面源污染在空间上的异质性。

3 化肥面源污染EKC分组模拟与预测 3.1 面板单位根检验和协整检验

面板数据兼具截面数据和时间序列数据的特点, 为避免伪回归, 在构建面板回归模型前需用单位根方法检验各个变量的平稳性和协整关系。常见的面板数据单位根检验有LLC、IPS和费雪式检验等方法, 本文运用以上方法对各组分别进行检验, 结果表明各种检验结论一致, 限于篇幅仅列出LLC检验结果。LLC检验结果(表 5)显示, 多数变量的原始序列存在单位根。因此对原始序列进行一阶差分, 结果显示各变量在一阶差分时为平稳, 序列为一阶单整序列I(1), 因此可以进行协整检验。

选取Kao检验方法对各组回归模型的二次曲线和三次曲线进行协整检验, 结果如表 6所示。各组ADF值均通过了显著性检验, 表明组内各个变量间的关系长期稳定, 可构建面板模型。

3.2 组间异方差、组内自相关和组间同期相关的检验

由于3类区域的分组数据均为长面板数据, 因此需要对各组数据进行组间异方差、组内自相关和组间同期相关的检验, 并根据结果确定合适的计量方法。检验结果(表 7)显示3类区域的二次和三次模型均存在组间异方差和组间同期相关; 而组内自相关检验表明, 热点区和非热(冷)点区的二次模型和三次模型均存在组内一阶自相关, 冷点区的二次模型和三次模型则不存在组内自相关。

3.3 化肥面源污染EKC空间分组拟合结果

为解决组间异方差、组内自相关和组间同期相关问题, 热点区和非热(冷)点区均采用全面FGLS方法进行估计, 该方法同时考虑了这3个因素, 能够保证估计的效率。对于不存在组内自相关的冷点区, 采用LSDV方法进行估计, 并使用组间异方差和组间同期相关的稳健标准误差来保证估计结果的稳健性。表 8报告了各组的回归估计结果。

回归模型结果显示, 热点区的三次曲线模型在10%显著水平未通过检验, 二次曲线模型在1%显著水平上通过检验, 故采用二次曲线形式; 冷点区的三次曲线形式模型在10%显著水平通过检验, 故采用三次曲线形式; 非热(冷)点区的三次曲线在5%显著水平通过检验, 故采用三次曲线形式。通过lnG的各次项的正负, 可以判断热点区的EKC形状为“倒U型”, 冷点区和非热(冷)点区为均“倒N型”。

1) 热点区。本组包含上海、江苏、浙江、河南、安徽、湖北等10个省份, 这些省份主要是我国经济发达地区或农业大省, 化肥面源污染强度普遍较大。本组EKC的时间路径呈现为“倒U型”曲线, 说明伴随农业发展, 化肥面源污染表现为“升—降”趋势。该曲线拐点对应的人均农业生产总值G分别为5 889.933元。从热点区10省的情况看, 江苏(2011年)、山东(2016年)、福建(2017年)、河南(2017年)、河北(2017年)5个省份已经跨过了“倒U型”曲线的拐点, 进入了“随着农业经济增长, 化肥面源污染减弱”的下降阶段。从时间上看, 多数省份在2017年左右跨过拐点, 其原因除经济不断增长促使化肥面源污染到达拐点外, 国家粮食安全战略下粮食主产区近年粮食播种面积比重的回升、农村居民人均纯收入的提高和农村居民对农业收入依赖的降低都可能促进了拐点的到来。此外, 2015年以来, 国家颁布的一系列化肥减量化政策, 如《到2020年化肥使用量零增长行动方案》以及近年中央一号文件对农业面源污染重视程度的不断加强也可能促进了拐点的到来。安徽、湖南、江西、上海和浙江则均尚未到达拐点, 仍处于“随着农业经济增长, 化肥面源污染增强”的阶段。安徽、湖南和江西作为我国重要的产粮区, 对于化肥的依赖程度较高, 难以在短期摆脱对化肥的依赖。上海和浙江作为粮食主销区, 较低的种植业比重反而致使化肥施用的集约化程度更低, 发展种养循环的动力更小, 从而未充分发挥有机肥代替化肥等方面的作用。热点区农业机械投入强度与化肥面源污染程度呈正相关, 表明科技和机械水平虽然在提高农业产值方面发挥了一定作用, 但是对于提高化肥利用效率和控制化肥流失率的作用尚未明显体现, 应该充分发挥农业新技术的作用, 如生态循环农业、配方施肥和精准施肥等。

2) 冷点区。本组包含甘肃、青海、西藏、黑龙江4个省份, 主要分布在我国西部和东北。本组EKC的时间路径为“倒N型”曲线, 说明伴随农业发展, 化肥面源污染表现为“降—升—降”趋势。该曲线的拐点对应的人均农业生产总值G分别为851.500元和4 433.746元。冷点区4省均在1997年以前就跨过了“倒N型”曲线的第1个拐点, 已经跨过第2个拐点的省份为甘肃(2017年)和黑龙江(2006年), 表明这两个省份目前处于曲线下降阶段, 即“随着农业经济增长, 化肥面源污染减弱”。青海和西藏则还未跨过第2个拐点, 仍处于曲线的上升阶段, 且表现出短期内不会达到拐点的趋势。对于青海和西藏而言, 其本身较低的农业经济增长水平促使其距离拐点较远。黑龙江农业产值高但化肥面源污染排放强度却一直较低, 其原因在于东北黑土资源肥沃、供肥能力强, 且多年来耕地质量普遍不断提高^[37], 化肥流失率低, 从而使黑龙江较早跨过拐点进入曲线下降阶段。

3) 非热(冷)点区。本组主要包括广西、云南、四川、吉林、辽宁、陕西等17个省份。本组EKC的时间路径为“倒N型”曲线, 说明伴随农业发展, 化肥面源污染表现为“降—升—降”趋势。该曲线两个拐点对应的人均农业生产总值G分别为1 490.698元和9 405.031元。从非热(冷)点区17省的情况看, 多数省份在1997年以前就跨过第1个拐点, 西南和华北部分省份则在2000—2010年间跨过第1个拐点, 进入了“随着农业经济增长, 化肥面源污染增强”的上升阶段。然而, 当前非热(冷)点区跨过第2个拐点的省份只有海南(2016年)和新疆(2016年)。表明当前非热(冷)点区的多数省份均处于曲线的上升阶段, 即“随着农业经济增长, 化肥面源污染增强”。按其曲线趋势, 非热(冷)点区各省的拐点相对热点区较晚, 且由于空间格局上近年空间自相关有所减弱, 因此应当注意非热(冷)点区和热点区之间空间相关性的变化, 讨论其相互间的影响。

3.4 化肥面源污染EKC时间路径分组预测

综合来看, 3类区域的曲线形态和拐点差异较大, 证明了空间异质性分析的必要性。截至2018年, 热点区仍有5个省份尚未跨过“倒U型”曲线的拐点, 处于曲线的上升阶段; 冷点区4个省份中还有两个未跨过“倒N型”曲线的第2个拐点; 非热(冷)点区仍有15个省份未跨过“倒N型”的第2个拐点。因此当前我国大部分省份正处于“随着农业经济增长, 化肥面源污染增强”的曲线上升阶段。在此情况下, 有必要预测这些省份EKC曲线的时间路径, 以根据其阶段特征和未来趋势采取相应的措施。本文运用趋势延伸法, 依据各省22年的人均农业产值年均增长率来预估下一拐点的出现年份, 结果如表 9所示。对于热点区省份, 除上海外, 其余省份可能在近年迎来拐点。对于上海而言, 其农业产业比重较低, 应当进一步加强产业结构的调整, 通过区域间产业转移等方式来寻求经济效益和生态效益的平衡, 获取更高的效益。对于冷点区, 青海和西藏距离拐点的时间较长, 且两省生态脆弱, 应进一步调整农业结构, 推广农业新技术的应用, 利用其畜牧业优势, 发展种养循环、生态农业, 使两省的拐点提前到来, 进入“倒N型”曲线的下降阶段。非热(冷)点区中除辽宁、吉林、内蒙古等省份外, 其他省份距离“倒N型”曲线的第2个拐点距离较远, 将较长期处于曲线上升阶段, 应当尤其引起重视。

4 讨论与结论 4.1 讨论

近年来, 化肥减量化和化肥面源污染防治成为我国农业工作的重点和研究的热点。随着农业面源污染对环境质量的影响越来越显著^[38], 基于区域异质性的化肥面源污染研究能使研究结果更加符合区域的实际情况, 从而更加准确地总结和预测我国化肥面源污染状况。当前有学者结合区域的实际情况对中国化肥面源污染环境风险时空变化进行分析, 发现各省差异较大, 地域联系较为明显, 呈现出聚集分布的特点^[39], 但并未探讨各区域差异的规律及原因。有学者对化肥施用及环境风险的时空变化特征进行了模拟, 得出了区域间差异化的结论, 然而并未考虑经济等因素对环境风险的影响^[40]。为研究区域异质性, 有学者通过面板门槛模型选取变量进行区域分组, 结合经济发展状况, 对化肥面源污染状况进行异质性条件下的EKC曲线检验^[22], 然而数理方法选取的分组变量可能存在主观性, 使得验证结果不能准确反映各区域的特质, 组间差异性也不明显。本文在我国化肥面源污染时空格局特征和相关性分析的基础上, 根据面源污染的空间异质性分组, 对化肥面源污染与农业经济发展之间的EKC时间路径进行模拟与趋势预测, 使研究结果更加符合区域的客观情况, 并得出了区域间差异化的结论, 其中一些现象值得进一步讨论。

1) 从化肥面源污染空间格局角度看, 我国化肥面源污染状况呈现出显著的空间正自相关性, 表明存在较强的空间集聚效应。近年来, 空间自相关性有所减弱, 长期的集聚格局可能发生一定变动, 不同区域间的流动性可能增强, 环境污染的空间溢出效应开始愈发显现^[41]。其原因可能是各区域间贸易^[42-43]和产业结构的调整和流动^[44-47]带来的污染物间的转移, 因此部分省份从热点区消失的部分原因可能是通过产业转移将污染物转移到了非热点区域。对于北京、上海等省市而言, 产业上的转移有可能提高区域的生产效率, 降低化肥面源污染的总体程度, 因此可以进一步适当推进农业产业的转移, 同时通过生态补偿、排污权交易等方式促进区域间的利益协调。从根本而言, 则应当进一步推广农业新技术的运用, 从根本上减少全国范围内的化肥面源污染。

2) 从EKC模型的分区模拟结果来看, 各区域在曲线形状和拐点上存在较大的差异。当前大多省份均处于曲线的上升阶段, 但距离拐点的时间则差异较大。热点区的多数省份距离拐点较近, 这表明热点区近年来化肥面源污染防控工作取得一定成效。但热点区迎来拐点后, 非热(冷)点区和冷点区却依然会在较长期内处于曲线上升阶段。其原因一方面可能是以上讨论的污染物空间上的转移, 另一方面则是非热(冷)点区的化肥面源污染防控工作亟待进一步加强, 应提高农业新技术的普及率, 推广有机肥和配方施肥。

3) 从控制变量的驱动作用来看, 3类区域农业机械强度均正向影响化肥面源污染, 这表明当前农业生产过程中, 机械水平的提高在提高农业生产效率的同时也增加了面源污染的风险, 因此在机械生产过程中应当提高精准施肥和配方施肥技术的运用。农村居民工资性收入的提高也使农业生产对于农民生计影响变小, 可能致使化肥施用更加粗放, 种养循环难度变大。此外, 相同因素在不同区域可能产生不同的影响, 该结论与崔鑫生等^[48]“同一因素对不同国家的作用存在着显著差异”的研究结论具有相似之处。粮食播种比例、产业结构的变化、种养结构的变化在不同区域会对化肥面源污染产生不同的影响。如畜牧业比重较大的冷点区, 种植业比重的下降明显地对化肥面源污染产生积极影响。

由于各区域在空间上的相关性, 因此在化肥面源污染治理中既要根据各区域的特性施策, 又要从整体上注重区域间的相互作用。对于热点区大部分省份, 应当进一步大力推广测土配方施肥, 精准施肥, 研发推广适用施肥设备, 机械深施, 提高化肥利用率。冷点区部分省份生态较为脆弱, 应实施保护性耕作, 增施有机肥, 推广水肥一体化等高效节水灌溉技术。非热(冷)点区应当注重种养循环, 推进有机肥的利用。从整体而言, 在化肥面源污染治理中, 应当注意区域间的相互影响, 尤其是注重产业结构调整可能带来的污染转移, 因此应当根据各区域省份的实际情况, 适当调整产业结构, 协调好经济效益和生态效益的平衡, 通过生态补偿、排污权交易等手段保障区域间发展均衡。

4.2 结论

本文在1997—2018年中国31个省(市、自治区)化肥面源污染排放强度核算的基础上, 运用ArcGIS空间自相关分析和热点分析空间统计模型, 分析了化肥面源污染排放强度时空格局演变特征, 并根据该特征, 将全国分为化肥面源污染排放强度的热点区、冷点区和非热(冷)点区3组, 在空间异质性条件下模拟和预测各组化肥面源污染和人均农业产值间的时间路径关系, 得出以下结论:

1) 我国化肥面源污染排放强度总体较高, 且省际差异较大。从时间上看, 各省化肥面源污染排放强度变化趋势各异, 表现出不同的升降趋势; 从空间上看, 空间自相关分析表明, 我国化肥面源污染排放强度表现出显著的空间正向自相关, 呈现为集聚模式, 集聚程度近年有所下降; 热点分析显示, 我国1997—2018年化肥面源污染排放强度时空格局演变较为稳定, 热点区主要分布在华中地区和华南地区, 黄淮海地区的热点程度下降较为明显, 冷点区主要集中在西部地区和东北的黑龙江。各区域间表现出一定的空间溢出效应。

2) 化肥面源污染排放强度时空格局演变特征分组的EKC模拟结果表明, 我国化肥面源污染与农业经济增长之间存在着显著的非线性EKC关系, 并且各组EKC时间路径曲线形状和拐点存在明显差异。热点区呈现为“倒U型”曲线、冷点区和非热(冷)点区呈现为“倒N型”曲线。我国大多数省份正处在化肥面源污染状况随农业经济发展而恶化的阶段, 且恶化的原因由于区域特质的不同存在差异, 亟需进一步因时因地推进各区域化肥面源污染防治, 并注意区域间的相互影响, 实现区域间的协同治理。

3) 化肥面源污染EKC的时间路径预测表明, 热点区上海、冷点区的青海和西藏以及非热(冷)点区的多数省份距离曲线拐点仍较远, 需因时因地施策, 促使拐点尽快到来。