国家尺度种植业面源污染负荷估算方法研究

邱捷，佘冬立，夏永秋

（1.河海大学农业科学与工程学院，江苏 南京 210098；2. 江苏常熟农田生态系统国家野外观测研究站，中国科学院南京土壤研究所，江苏 南京 210008）

农业生产活动中产生的溶解态或固态污染物，在降雨径流的冲刷和淋溶作用下进入受纳水体，形成农业面源污染[1]。20世纪70年代以来，中国水环境质量急剧恶化，面源污染日趋成为水体富营养化的主要原因。有关研究表明，农业源氮磷污染对太湖的贡献率分别达到77%和33.4%[2]，对巢湖贡献率分别为53%和42%[3]。然而，面源污染因其形成过程受区域地形、气候、土地利用和植被覆盖等多种因素影响，具有随机性、分布范围广和产生机理复杂等特点[4]，这些特点使得面源污染在监测、模拟与控制方面面临巨大的困难。

空间尺度上，面源污染估算可用分为田块尺度、流域尺度和区域尺度。田间尺度种植业面源污染负荷估算研究大多是以典型田块为基础，开展典型田块的径流动态特征观测，基于管理措施如施肥、灌溉等田块理化性质与作物生长、降雨强度等进行估算[5]。流域尺度上种植业面源污染负荷估算是基于代表性流域为研究对象，开展流域尺度的河流断面与流域出口水质水量观测，基于流域土地利用、景观格局、土壤、管理和气象等因子空间变异进行估算[6]。区域尺度上种植业面源污染负荷估算是通过调查区域典型水体水质时空变异，将区域划分成网格、行政单元和流域，基于种植制度、气候条件和人类活动等因子空间变异进行估算[7]。因此，区域尺度是掌握整个地区或国家种植业源污染负荷概况、检验种植业面源污染防治措施成效、预测面源污染发展趋势的重要尺度[8]。

目前，面源污染的量化研究大多基于田块尺度和流域尺度，由于面源污染负荷有很强的时空变异性，将小尺度估算结果放大到区域污染负荷估算有很大的不确定性[9]。随着面源污染损失规律研究不断深入，越来越多全国尺度上的估算方法出现。从早期经典输出系数法到考虑土地利用变化、水文等因素的改进系数法，从经验统计模型到复杂机理模型。本文简要综述了国家尺度种植业面源污染负荷的估算方法，结合研究结果分析了我国种植业面源污染防治中存在的若干问题，提出了面源污染发展的优先研究领域和重要方向，可为种植业面源污染控制和管理的定量化研究提供参考。

1 估算方法

1.1 输出系数法

输出系数法是指利用一个单一的系数或系数集来表征区域种植业面源污染负荷强度。20世纪70年代初期，美国和加拿大在研究土地利用—营养负荷—湖泊富营养化关系的过程中，首次提出并应用输出系数法，这种方法为人们估算面源污染提供了一种新的途径[10-11]。但早期的输出系数法土地利用分类简单，精度不高，只能用于估算土地利用类型较为单一区域内的面源污染。针对这一不足，Johnes[12]对不同作物类型提出不同的输出系数，并利用它对英国Windrush流域的氮磷损失量进行估算。这种改进方法丰富了输出系数模型的内容，提高了模型对土地利用状况变化的灵敏度，因此被广泛应用于面源污染负荷研究。在我国，朱兆良[13]开展了长期定位试验，借鉴Johnes经典输出系数法，建立了我国1979年和1998年农业生态系统氮素损失的输出系数，其中氮素淋失系数约为化肥氮投入量的2%，径流系数约为化肥氮投入量的5%。

1.2 改进系数法

传统的输出系数法假定农田输出系数固定不变，这种假设导致估算结果和现实状况差异很大[14]。为提高估算精度，胡玉婷[15]在输出系数模型中将农田区分为水田和旱地，对两种典型土地利用方式的径流和淋洗损失采用不同的输出系数，其中水田氮素淋失率和径流率分别为投入化肥氮的3.39%和5.15%，旱地氮素淋失率和径流率为投入化肥氮的9.97%和5.15%。第一、二次全国污染源普查[16-17]进一步将37种土地利用方式的径流和淋失区分开，在全国不同地区选取代表性地形地貌，设置数百个地表径流和地下淋溶监测试验点，从而得出不同土地利用类型的径流和淋失系数。上述输出系数模型都不考虑面源污染负荷的滞后性，为此，Worrall和Burt[18]提出考虑土地利用变化对污染物输出系数影响的氮流失动态模型，减小了预测误差。面源污染除随水迁移之外，还随泥沙进行迁移，尤其是易被泥沙吸附的磷素[19]。由于泥沙输移率在迁移过程中逐渐减小，泥沙颗粒中吸附态磷在泥沙输移过程时无法全部进入受纳水体。为此，Soranno等[20]提出考虑磷迁移损失的改进输出系数法来预测和评价面源磷负荷。

1.3 多元回归法

尽管输出系数法经过不断的扩充和发展，但只是在土地利用方式上的细化，而农田氮磷流失过程十分复杂，只考虑施肥量单一影响因子，忽略气象条件、地形地貌、土壤性质和作物管理等众多因素，在估算上具有很大不确定性。基于此，杨旺鑫[21]通过收集不同区域关于农田氮磷损失的文献，提出氮磷流失量与土壤全氮/磷含量、粘粒含量、氮磷施用量和降雨量之间的多元回归模型。夏永秋等[22]收集了我国大量农田氮素径流和淋洗损失数据，应用逐步回归分析方法，提取影响种植业面源氮负荷量的气候、地理、作物管理和土壤等因子，构建基于氮肥用量、降雨量和土壤粘粒含量的多元回归模型。结果表明，与仅考虑化肥施用量的输出系数模型相比，增加降雨量因子或土壤粘粒含量因子后的多元回归模型可大大提高国家尺度种植业面源污染负荷估算的准确性。

1.4 贝叶斯递归回归树模型

由于区域尺度上面源污染负荷变异大，影响因素众多，通过一个或几个系数的模型很难解释面源污染负荷的空间变异。贝叶斯递归树的思想是，将大问题分解为小问题来求解，然后再将小问题分解为小小问题。这样一层一层地分解，直到大变异被分解得足够小，在每一个小的变异范围内，都可以通过简单的因子进行充分解释，不用继续递归分解为止[23]。基于这种思路，周丰和周怀成[24]提出“多输入—单输出”的国家尺度种植业面源污染负荷贝叶斯递归回归树模型（Bayesian Recursive Regression Tree, BRRT）（式1），用来弥补输出系数模型和多元回归模型的不足。

式中：p(T|X,Y)为面源负荷的响应决策树，X为影响因素，Y为面源负荷，θ为模型参数，p(T)和p(θ|T)分别为决策树和参数的先验分布，p(Y|X,θ,T)为似然函数。

BRRT原理是通过随机搜索和局部贪婪搜索，以似然函数值最大为目标，建立因变量和自变量之间的时间或空间线性分段函数。Hou等[25]和高硕硕等[26]基于大量站点农田氮素负荷与影响因素，利用BRRT分别拟合氮素径流和淋溶损失的分段函数模型（PKU-NLEACH和PKU-NRUNOFF）。

1.5 过程模型模拟法

上述全局或者分段回归模型，都是建立面源污染负荷与主要影响因素之间的经验关系，仅代表区域内面源污染负荷的平均水平，缺乏对污染物负荷的机理表达[27]。为克服这一方法的局限性，Fu等[28]在大量田间观测和调研数据的基础上，将土壤—水界面的氮磷动态过程纳入面源污染负荷过程模型（Rice Paddy Runoff Module, RPR）中，并得到全国主要水稻产区11个地点127场降水径流事件田间观测的验证。RPR模型主要包括水分动态平衡、氮磷传输过程和灌溉施肥管理三个部分，采用水文循环、生物地球化学过程和作物物候学的耦合方程进行建模，用于模拟网格单元或农田尺度上的氮磷日径流通量，最后根据模型参数如气候、土壤、土地覆盖、移栽和收获日期、肥料施用率以及施肥和灌溉方案等数据在全国尺度上的分布，模拟了我国中稻、早稻和晚稻不同种植制度的稻田氮磷径流流失时空格局和减排潜力。

2 估算结果与分析

基于上述方法，众多研究者在全国尺度上对我国种植业面源污染负荷估算进行了研究（表1）。结果表明，我国总氮径流损失估算结果在0.30～2.40 Tg之间，总氮淋洗损失为0.36～2.03 Tg，总磷径流损失为3.5～6.37 Mg，由于总磷易被土壤吸附，总磷淋洗损失一般很低，很少有关于全国尺度总磷淋洗损失的估算。面源污染排放过程较为复杂，受地形、地貌、气候和水文等因素影响，同时，估算方法、数据来源和影响因素选择等的不尽相同，导致核算结果之间存在较大差异。

2.1 全国尺度种植业氮素径流损失比较

对于全国尺度种植业氮素径流损失，朱兆良[13]通过输出系数模型，估算出全国种植业氮素径流负荷为1.2 Tg。胡玉婷[15]基于改进的输出系数法对全国224组农田氮素损失数据进行统计分析，得到氮素径流损失量为1.31 Tg。Wang等[29]根据地表水系统各个流域农田输入的空间分布，采用改进输出系数法计算农田径流，得到2017年全国总氮径流输出量为2.4 Tg。Cui等[30]对462篇已发表文献和3 374个观测值组成的数据集进行分析，建立氮素损失与施氮量之间的指数模型，估算全国稻田氮素径流损失量为0.17 Tg。Sun等[31]通过分析全国232个监测点的测量数据，获取2004年至2017年我国氮肥的流失系数和流失量，研究发现，地表径流氮肥排放量从0.33 Tg降低至0.30 Tg。夏永秋等[22]应用逐步回归方法对全国187组文献数据进行多元回归分析，模拟得到2011年全国种植业氮素径流损失量为0.96 Tg。杨旺鑫[21]采用多元回归模型对2011年我国各省水田和旱地化肥使用的氮素损失进行研究，模拟全国水田旱地总氮径流总损失合计2.37 Tg。Hou等[25]基于BRRT算法校验后的经验模型对环境因素作用于氮素径流损失进行模拟，估算2008年中国农田总氮径流量为1.09 Tg。Fu等[28]建立我国稻田施肥、灌排和耕作制度数据集，通过开发的RPR模型估算全国稻作系统中氮素径流流失量为0.27 Tg。

表1 国家尺度种植业面源污染负荷估算方法及结果对比Table 1 Estimation methods and results comparison of national non-point source pollution loading from planting system

2.2 全国尺度种植业氮素淋洗损失比较

对于全国尺度种植业氮素淋洗损失，朱兆良[13]基于输出系数法，估算出全国种植业面源氮素淋洗负荷为0.5 Tg。胡玉婷[15]区分稻田和旱地的损失系数，得到全国氮素淋溶损失总量为2.03 Tg。Sun等[31]通过改进输出系数模型对全国氮素淋洗损失进行估算，结果表明，淋溶损失的氮肥排放量从2004年的0.38 Tg降低至2017年的0.36 Tg。夏永秋等[22]构建氮素损失及其影响因素的多元回归模型，计算得到2011年全国种植业氮素淋洗损失量为1.01 Tg。杨旺鑫[21]利用农田实验数据对我国2011年种植业源氮素损失进行多元回归分析，估算我国种植业总氮淋洗损失为1.80 Tg。高硕硕和周丰[26]通过BRRT模型检验并量化了我国农田氮素淋溶损失与施肥量的关系，计算出2008年我国农田上因施肥量引起的总氮淋溶量为0.88 Tg。第一次全国污染源普查[16]对2007年我国内地31个省市自治区种植业氮素流失潜力进行估算，得到全国农田氮素损失量为0.71 Tg。第二次全国污染源普查[17]得到2017年我国种植业氮素总损失量为0.72 Tg。

2.3 全国尺度种植业磷素径流损失比较

对于全国尺度种植业磷素径流损失，Sun等[32]通过改进输出系数法对全国尺度种植业磷素径流损失进行估算，2004年、2008年和2013年我国种植业磷素损失量分别为3.5 Mg、3.8 Mg和4.6 Mg。杨旺鑫[21]基于多元回归法估算2011年我国水田和旱地磷素径流总损失合计6.37 Mg。Fu等[28]采用过程模型模拟法对中国稻田总磷径流损失的空间格局和缓解潜力进行模拟研究，估算全国稻作系统中总磷径流流失量为1.70 Mg。第一次全国污染源普查[16]估算2007年我国种植业总磷损失量为37 Mg，第二次全国污染源普查[17]得到2017年我国种植业总磷损失量为76 Mg。

3 估算方法中存在的问题

3.1 缺乏沿程消纳过程模拟

种植业面源污染的成污过程是农田中的氮磷营养元素以径流、淋洗等形式到达地表、溪流、沟渠、池塘、下水道和地下水等，经输移到达河流、水库、湖泊和海洋，成为污染[33]。一般来说，田块产生的面源污染物无法全部到达流域出口，污染物在向下游输移过程中由于土壤植被的截留、泥沙沉淀、地下渗漏和微生物分解等过程逐步消纳[34-36]。现有国家尺度面源污染负荷模型，由于数据源大多基于田块流失量的监测，对污染物伴随径流产生与汇集时向流域出口迁移过程中的损失考虑不足，未能综合流域地形地貌、水文气候、植被覆盖和土地利用等因素的影响，因此计算出的结果严格上仅为流域面源污染物发生量, 并非流域出口的实际污染负荷。Khadam和Kaluarachchi[37]提出考虑水文因素和流域损失的改进输出系数模型，并将此模型应用于美国华盛顿州，发现该模型模拟准确性明显提高。因此在核算我国种植业面源污染总氮和总磷负荷量的研究中，应对区域污染物在输移过程中的损失加以考虑，进一步提高核算结果准确度。

3.2 未区分背景排放和肥料排放

农田土壤氮磷流失是一个复杂的过程，施用化肥会显著增加氮磷流失量[38]。在不施肥情况下土壤本身也会因有机态氮磷矿化作用释放出未被作物吸收的无机态氮磷污染物，以径流或淋洗形式流失，称为背景排放[39-40]。许多学者在计算全国面源污染负荷时通常不区分背景排放和肥料排放，从而未能客观评价施肥对面源污染的贡献。倪玉雪[41]对旱地农田氮素径流和淋洗损失的研究表明，把握背景排放的影响因素规律对准确评估施入氮在农田径流和淋洗损失中的贡献有重要作用。土壤背景排放影响因素众多，其排放变异大。研究发现，土壤理化性状对土壤背景排放有较显著影响，作物生育期土壤水分、温度和养分矿化等情况的不同，导致土壤背景排放产生差异，其变化范围为11.4～12.9 Mg，进而使得背景排放对全国农田氮素排放总量的贡献率在20%～53%内变化[42-44]。此外，在气候、地理和土壤等多种因素影响下，部分肥料在固定—矿化等氮循环过程下转变为土壤氮库，以背景排放的形式流失。因而区分出土壤背景排放，把握背景排放和肥料排放的转化规律，为获得更为准确的土壤养分背景排放因子和化肥源排放因子提供依据，对开展全国尺度种植业面源污染负荷估算研究具有重要的意义。

3.3 未考虑南方一年多熟制

由于农业气候带的多样性，我国种植制度种类繁多，从一年一熟、两年三熟、一年两熟、两年五熟到一年三熟不等[45-47]。当前，全国尺度种植业面源污染估算大多基于田块尺度下一年一熟进行估算，对于南方一年多熟考虑较少，导致估算结果偏低。实际上，相比作物一季氮磷流失量，一年多熟种植制度下化肥投入量成倍数增加，污染物流失量也应成倍增加，因此用作物单季排放量不能准确描述污染全年排放规律。Xing[48]通过对已发表文献中数据进行整理分析，探讨旱地和水田氮素排放规律。结果表明旱地一年一熟、一年两熟和一年三熟排放量表现出明显差异，水田单季稻、双季稻和稻麦轮作排放量同样存在较大差别。Liu等[49]通过对我国亚热带丘陵区双季稻作系统径流过程进行1.5年连续高频监测，研究发现早稻季水田氮磷浓度较高，氮磷径流损失大于晚稻生长季。朱兆良[50]研究也同样发现，氮素淋洗过程是一种累进过程，在当季未被淋失的养分可在后续作物季释放出来。因此，在估算全国尺度种植业面源污染负荷时，应充分考虑种植制度的影响。

3.4 缺乏氮磷流失过程建模

由于强烈的时空变异特征，现有国家尺度面源污染负荷模型大都基于容易获取的气候、地理、作物管理和土壤等因子建立模型，缺乏对污染物流失过程的机理描述，使得模型预测性能较弱，只能对纳入模型的因子做简单的预测。此外，基于历史数据建立的统计模型仅限于参数分布在其建模时的参数范围内应用，一旦参数超过这个范围，模型将无法适用[51]。过程机理模型从土壤活性氮磷损失的过程和机理出发，阐述气候、地理、作物管理和土壤等因子的时空变异对面源氮磷损失的影响，得到用数学公式表达的物理、化学和生物过程[52]，能够科学地描述农田污染排放的时空差异性，并通过设置不同环境条件预测不同情景下区域污染损失情况[53]。然而过程机理模型所需参数较多，且要求详细明确，大量输入资料获取困难，参数选择和模拟精度上存在一定的问题，同时测试和相关费用较大，且必须通过长期和多点实测数据校验[54]，因而在我国区域尺度上应用存在一定的局限性。因此，如何根据研究区域特征模拟非点源产排污机理，建立方法简便但估算准确性较高且适合我国面源污染特点的预测模型是今后研究的重点。

4 结语

面源污染是导致水体和农田生态系统遭受破坏的重要因素。经过多年探索，种植业面源污染负荷研究已从定性化分析转向定量化核算，研究尺度由田块尺度和流域尺度扩展到国家尺度。纵观目前国家尺度面源污染负荷估算方法，已经形成输出系数法、改进输出系数法、多元回归法、贝叶斯递归回归树模型、过程模型模拟法等模型方法。各方法之间研究结果差异较大，具有很大的不确定性。

我国国家尺度种植业面源污染流失量估算方法总体上存在以下主要问题，源头排放未区分背景排放和肥料排放，也未考虑熟制，过程削减未充分考虑沿程消纳。针对上述问题，国家尺度种植业面源污染负荷模型应进一步明细各土地利用输出系数、区分背景排放和肥料排放、考虑我国熟制的区域差异、考虑田块到流域出口的沿程削减。为了同时兼顾模型的实用性和机理性，还应加强国家尺度面源污染半经验半机理模型研发，推进我国种植业面源污染负荷本土化模型开发和应用研究，为区域种植业面源污染控制提供科学依据。