氮既是农业面源的主要污染物(McIsaac et al., 1991;Alberts et al., 1981;de Wit et al., 1999;Bardgett et al., 2001; Sharpley et al., 1987)，也是影响农作物产量的重要营养元素(Vitousek et al., 2009).随着我国农业生产集约化程度的不断提高，肥料投入量逐渐加大，农田氮素施用量(合成化肥+复合肥)在2008年达到4037万t，约占全球的41%，施N强度为232 kg · hm^-2，且农田泥沙流失速率(10~100 Mg · hm^-2 · a^-1)远高于其他透水地面(牧草地、天然草地、森林，0~0.5 Mg · hm^-2 · a^-1)，降雨产生的农田地表氮流失量逐渐增加(Quinton et al., 2010)，导致我国的江、河、湖、库及近海水域都面临着富营养化的威胁(黄文钰等，1998；赵亮等，2002；Yan et al., 2003)，并引起了许多水环境问题(Smith et al., 2001；陈志凡等，2006).

TN排放负荷的精确估算是解决水环境问题的基础，近年来，国内外许多学者对氮素污染的排放负荷及其与各个影响因素之间的响应关系展开研究(Liu et al., 2013; Ng Kee Kwong et al., 2002; Goolsby et al., 2000; Chen et al., 2004; Lewis et al., 2003).从估算方法上看，无论是早先提出的简单统计模型方法(郝芳华等，2006)，还是后期的机理模拟组合模型，如SWAT、HSPF、USLE等都有特定的适用范围.在大尺度区域中，由于各个地区的自然条件及社会发展状况不同，简单统计模型参数缺乏相应的空间分异性，精确度低，难以对污染产生过程机理进行深入的分析与研究；机理模拟模型对数据规模和精度要求过高，不太适合大尺度区域模拟.从活动数据的空间分辨率来看，目前还没有利用县域尺度活动数据估计我国氮地表流失负荷的研究；此外，还需要全面评估活动数据、流失率和参数的不确定性，以此估计氮流失负荷的可靠性.

已有研究表明，苏南太湖流域氮对地表水的污染负荷量高达2.5×10⁴ t · a^-1，占氮素化肥施用量的16.8%.太湖和巢湖面源污染物对TN的贡献率达到59%和63%(朱兆良等，2005)，其中，农业面源污染更突出，已成为水体富营养化最主要的污染源(张永春等，2006).因此，有必要对该地区农田TN排放负荷进行估算.考虑到区域的完整性，本文选取长江三角洲周围“三省一市”(上海、江苏、浙江和安徽)的农田为研究区域.以往的研究主要是对农田氮肥的流失做点位尺度或有限数量的田间试验来估算流失率，本研究则基于这些试验数据对2008年农田肥料氮素的地表流失率及其各个影响因子的驱动机制进行探索，结合GIS技术估算其肥料氮素的地表流失负荷，以期为该区域农业面源污染的控制和综合治理提供数据支撑.

本文研究区域为长江三角洲地区，选择上海、江苏、浙江和安徽共225个县(区、市)的农田，面积35万km²，其中，耕地面积15.4万km².数据空间分辨率为县级单元，基准年为2008年，农用肥料为合成化肥(氮肥和复合肥).

农田TN流失负荷的计算公式如下：

式中，R为农田TN流失负荷(Gg · a^-1，以N计)；NA为农田肥料施用量(折纯N，Gg · a^-1)，NA=NS + a·NC，其中，NS和NC分别为合成氮肥(折纯N)和复合肥施用量(Gg · a^-1)；a为复合肥中含N量，设为30%；η_i为农田TN流失率；RA为常规施肥TN流失量(kg · hm^-2 · a^-1)，其中，CK为不施肥情况下TN流失量(kg · hm^-2 · a^-1)；R_NA为年单位面积NA值(kg · hm^-2 · a^-1).

由于η_i的影响要素众多，且具有较强的空间分异性，因此，本文采用基于全局随机搜索算法原理的贝叶斯递回归树模型(Bayesian Recursive Regression Tree，BRRT，周丰等，2010)，建立影响要素不同水平下我国农田氮素净流失比例η_i的响应关系.在多次重启全局随机搜索和局部贪婪搜索下，使得趋于常数，其中，p(T|X，η)为考虑流失率的空间分异性的响应决策树，X为影响因素，η为流失率，为模型参数，p(T)和p(|T)为决策树和参数的先验分布，p(η|X，，T)为似然函数，具体原理参考文献(周丰等，2010；郭怀成等，2012).为简化计算，仅利用BRRT的全局分类功能，操作中所有变量并不进入节点回归拟合，具体参数设置包括：20次重启(经验而定)，每次10万次迭代(之后迭代对似然函数值改进水平<1%)，决策树先验分布的参数α和β及参数c和λ分别为0.95、1、1和0.1173(为敏感性分析最佳组合).BRRT的训练和验证样本皆来源于公开发表文献，最终得到不同叶节点的TN流失率，即不同分区的流失水平及不确定性范围.

从上述方法原理可知，该模型具有2个特点：①具有更高的空间分异性，一则模型输入的地形、土壤、土地利用、气象、养分管理和农田氮素污染源具有空间分异性，其次，影响要素不同水平下的响应方程也不同；②响应方程是全局和局部寻优的结果，使得响应方程具有更高的可靠性和适应性.

各县(区、市)的氮肥、复合肥数据来自于48个地级行政区2009年社会经济统计年鉴或国民经济统计和社会发展统计公报，获取途径为国家图书馆和中国社会经济发展统计数据库.2008年长三角地区县域农田肥料(氮肥和复合肥)折纯N施用量见图 1.

图 1(Fig. 1)

图 1 研究区域、数据观测站及2008年县域农田肥料折纯N施用量分布 Fig. 1 Study area，data observation stations and agricultural fertilizers of counties in 2008

数据来源于维普期刊网和中国期刊网的公开文献，获得研究区文献35篇(2002—2011年)，共188组 TN地表流失率及其影响因素数据，其试验地点分布于整个研究区域(图 1).相应地，BRRT中因变量η_i为TN流失率，自变量x_i为其影响因素，包括降雨量(x₁，mm)、化肥施用量(x₂，kg · hm^-2 · a^-1)、土地类型(x₃，其中，0表示旱地，1表示水田)、TN(x₄，g · kg^-1)、pH(x₅)、有机质(x₆，g · kg^-1)、粘粒比重(x₇)，本文所有试验数据中坡度均为0.上述土壤理化性质数据来源于Harmonized World Soil Database(1 km×1 km)和ISRIC World Soil Information(5'×5')，年降雨量来源于中国气象科学数据共享服务网，各县的灌溉面积是基于2009年中国农业年鉴进行折算.所用二级分类的土地利用数据源于中国科学院地理与资源研究所(1 km×1 km，刘纪远等，2009).

数据处理：在ArcGIS平台上，对降雨量站点通过创建泰森多边形插值得到研究区域的降雨空间分布图，结合土地利用类型数据将每个县施氮量平均分配到该县的耕地上，最终得到1 km×1 km的各个影响因子栅格图；再根据得到的流失率驱动机理模型，利用栅格计算模块进行叠加分析，得到该区域农田氮素地表流失率及其流失负荷的空间分布图.

由于统计口径、区域差异性等原因，上述活动数据、径流率和模型相关参数是存在不确定性的.活动数据的偏差是根据县级数据总和与国家及FAO统计数多年(1978—2008年)数据差值的平均值来确定，相应的氮肥和复合肥的偏差分别为5.0%和0.6%；径流率和相关参数的不确定性范围则是根据本土实测值的样本而统计确定，选择R₅₀(25%和75%分位数的差值)为其不确定性水平.

采用文献中188个样本(148个样本用于校准，40个样本用于验证)校准和验证BRRT，结果表明，在校准和验证中的可决系数分别为0.820和0.744(图 2)，达到较为满意的模拟效果，可用于计算流失率.相应地，决策树形成17个响应方程(表 1).

图 2(Fig. 2)

图 2 长三角地区农田TN流失率的BRRT模型校准(a)与验证(b) Fig. 2 Calibration(a) and verification(b)of BRRT model for total nitrogen runoff rate of cropl and s in theYangtze River Delta area

表 1(Table 1)

表 1 BRRT模型相应的响应方程及适应条件 Table 1 Response equation and adaptation conditions of the BRRT model

表 1 BRRT模型相应的响应方程及适应条件 Table 1 Response equation and adaptation conditions of the BRRT model 编号 流失率 适用条件 1 y=0.00207±0.00113 x1≤72.5且x4≤1.8 2 y=0.01782±0.00455 x1≤72.5且x3=0且x4>1.8且x5>5.6 3 y=0.01782±0.00586 72.5< x1≤80且x3=1 4 y=0.00867±0.00445 80< x1≤140且x3=1 5 y=0.02700±0.00735 140< x1≤300.1且x3=1 6 y=0.00773±0.00418 300.1< x1≤470.3且x3=1 7 y=0.07015±0.00596 470.3< x1≤491且x3=1 8 y=0.02585±0.00818 491< x1≤957且x3=1 9 y=0.017859±0 x1>957且x3=1且x6≤29.54 10 y=0.06409±0.00860 x1>957且x3=1且x6>29.54 11 y=0.04238±0.00445 x1>72.5且x3=0且x4≤1.52且x7≤15.6% 12 y=0.00475±0.00290 x1>72.5且x3=0且x4≤1.52且x5≤6.99且 x7>15.6% 13 y=0.02374±0.00940 x1>72.5且x3=0且x4≤1.52且x5>6.99且 x7>15.6% 14 y=0.05201±0.01771 72.5< x1≤250.8且x3=0且1.52< x4≤2.88且x6≤38且x7>26.6% 15 y=0.07544±0.00873 x1>250.8且x3=0且1.52< x4≤2.88且x6≤38 16 y=0.01500±0.00262 x1>72.5且x3=0且1.52< x4≤2.88且x6>38 17 y=0.05201±0.00745 x1>72.5且x3=0且x1>2.88

根据上述数据处理方法估算得到1 km×1 km的农田肥料TN地表流失率与流失负荷(图 3).通过观察可知，流失率在空间上均存在着明显的变化(图 3a)，淮河北部及江苏东部沿海区域虽然降雨量不高，但土壤中TN含量高，超过1.52 g · kg^-1，致使部分区域流失率高达7.54%；长江两侧流失率相对来说普遍较低，最低处可达0.47%，原因是这部分农田土壤中TN含量相对较低；白色区域为非耕地，多集中在安徽的南部和浙江的大部分区域，长三角地区农田肥料TN平均流失率为3.36%.

图 3(Fig. 3)

图 3 长三角农田TN流失率(a)与流失负荷(b)空间格局 Fig. 3 Spatial pattern of total nitrogen runoff rate(a) and load(b)of cropl and s in the Yangtze River delta area

如图 3b所示，对于负荷总量而言，流失负荷较多的县(区、市)主要集中在淮河北部及江苏东部沿海区域，与流失率的空间分布相一致，这些区域正是我国粮食主产区，化肥施用强度相对较高；而流失负荷较少的地区主要分布浙江及安徽南部，且流失量大多在100 t · a^-1(以N计)以下，这片区域耕地较少，多为山地丘陵地区，化肥施用量也相对较低.研究区域农田肥料TN流失负荷总量为88.1 Gg · a^-1(71.9~104.4 Gg · a^-1，以N计)，平均达到7.91 kg · hm^-2(以N计)，江苏、安徽、浙江和上海分别贡献了54.3%、40.7%、3.7%、1.3%；从市域来看，流失负荷最大的市依次为盐城(江苏省，12.9 Gg · a^-1)、徐州(江苏省，10 Gg · a^-1)、阜阳(安徽省，7.7 Gg · a^-1)、亳州(安徽省，5.6 Gg · a^-1)，占研究区域农田肥料TN流失负荷的41%；从县域来看，流失负荷最大的依次是大丰市(江苏盐城市，3.4 Gg · a^-1)、沭阳市(江苏宿迁市，2.9 Gg · a^-1)、射阳县(江苏盐城市，2.8 Gg · a^-1)、东台市(江苏盐城市，2.3 Gg · a^-1)，占研究区域农田TN流失负荷的13%.

为了考察不同因素对研究区TN流失率的影响，计算了7个因素与流失率的Pearson相关系数(表 2).结果表明，农田中TN流失率与降雨量、土壤有机质含量及粘粒比重显著正相关(p<0.01)，与施氮量显著负相关(p<0.01).降雨量和土壤团聚体是决定农田TN流失率大小的重要因素，降雨强度越大，地表流失模数与径流模数增大，氮素流失越多(段小丽等，2012).土壤团聚体的富集又是流失泥沙富集颗粒态氮的主要原因(黄满湘等，2003)，而土壤有机质含量及粘粒比重又影响了土壤团聚体的形成(王清奎等，2005).流失率是TN流失量与施氮量的比值，因而施氮量的增加直接导致流失率的减小.土壤中养分流失的重要途径之一是随水分流失，由于水田的径流和积累径流量大，因而其TN流失大于旱地(于兴修等，2002).土壤TN含量的提高也有利于流失泥沙对N素养分的富集，因而与流失率呈正相关关系(段然等，2013；段永蕙等，2004；付伟章等，2005).由于土壤胶体是带负电荷的，随着土壤pH的增加，土壤胶体对NO^-₃-N的吸附能力下降，从而增大了表层土壤中NO^-₃-N的渗漏流失量，减小径流流失量(张国梁等，1998).因此，土壤pH与TN地表流失量呈负相关(张燕等，2012).农田地表TN的流失是一个复杂的综合过程，仅仅用这些简单的影响因素进行相关分析，很多现象难以解释，如水田与旱地的区分、土壤TN含量及pH值与TN流失率的关系，故本文将土地类型、TN及pH作为方程变量.

表 2(Table 2)

表 2 TN流失率及影响因素的相关性分析 Table 2 Correlation analysis of total nitrogen runoff rates and the corresponding factors

表 2 TN流失率及影响因素的相关性分析 Table 2 Correlation analysis of total nitrogen runoff rates and the corresponding factors 化肥施用量 土地类型 TN pH 有机质 粘粒比重 流失率 注：所有检验的样本量均为188;*p<0.05，**p <0.01水平(双侧). 降雨量 0.329** 0.070 0.055 -0.154* 0.008 -0.032 0.297** 化肥施用量 1 -0.022 -0.040 -0.242** 0.001 -0.112 -0.235** 土地类型 1 -0.210** 0.071 0.444** -0.047 0.121* TN 1 0.161* -0.044 -0.066 0.189* pH 1 0.173* 0.066 -0.131* 有机质 1 0.264** 0.265** 粘粒比重 1 0.248**

目前，有关地表TN流失的研究多集中在小流域，与本次的研究区不相匹配；而且本次研究中流失量的计算是将施肥农田TN流失量减去对照试验中不施肥农田TN流失量，所以与以往的地表TN流失研究略有差别.为达到分析验证的目的，本文采用2种方法：①TN流失量与流失率的空间形态对比验证；②纵向剖面数据验证分析.

TN流失负荷空间分布的合理性分析验证主要基于土壤流失过程是造成这类污染物发生迁移、输送直至产污的最直接动力这一思想而展开.因此，流失率高的地区其TN的流失水平也相应较高；反之，亦然.对照长三角地区农田TN流失率图(3a)与流失负荷图(3b)，总体而言，这两者在空间上存在很大的相似性.重点流失区，如淮河北部及江苏东部沿海区域，同时也是TN流失负荷高值区；轻微流失区，如长江两侧，其TN流失水平也较低.

分别对农田肥料TN流失率、流失量和施氮量空间图沿经线作纵向剖面划分及数据提取.剖面地理位置为东经120°，纬度27°45′~34°19′之间，从北向南依次跨过淮河、长江.在对剖面数据作去空白处理后，即剔除非农田区域，选择其中100个不连续单元数据网格(1 km×1 km)，其中，1~40号位于长江以南，41~80号位于长江与淮河之间，81~100号位于淮河以北(图 4).1~40号的N肥流失率最低，其施氮量也最少，故而其N肥流失量最少；41~80号施氮量最高，N肥流失率位于三者之间，流失量位于三者之间；81~100号N肥流失率最高，这是由于该地区土壤TN、有机质含量高及施氮量高所导致的；而41~80号N肥流失量却低于81~100号，可见淮河以北的N肥流失率决定了其流失量的高低.对上述100个网格单元的数据对比分析，发现不同气候带上网格单元的N肥流失状况符合其自然流失规律，即N肥流失率和施氮量是两大决定性因素.因此，剖面数据的合理性进一步证实了本研究估算结果的可靠性.

图 4(Fig. 4)

图 4 纵向剖面单位网格内施氮量、TN流失量及流失率 Fig. 4 Vertical profile of nitrogen application，runoff rate and load per unit grid

本估算是初步尝试，必然存在较大的不确定性水平：①仅考虑了农田合成氮肥和复合肥，缺少有机肥的贡献；②缺乏考虑种植模式、施肥方式、施肥季节降雨特点(雨量大小和持续时间)等的影响；③本估算方法仅用到BRRT的全局分类功能，没有真正建立基于全局分类的非线性响应方程(Nolan et al., 2006)，而这正是今后需要改进的研究要点.

基于高分辨率的活动数据、具有空间分异性的本土影响因素值，采用GIS方法，初步估算了2008年长三角地区农田肥料引起的TN流失率和流失负荷，分别为3.36%(3.09%~3.63%)和88.1Gg · a^-1(71.9~104.4 Gg · a^-1，以N计).流失率和流失负荷在空间上均存在着显著的分异性，在淮河北部及江苏东部沿海区域流失较为严重，长江两侧流失率相对来说普遍较低，浙江及安徽南部因耕地较少，因此，流失负荷也相对较小.其中，影响流失率的最主要因素是降雨量，其次是土壤有机质含量、土壤粘粒比重和施氮量，同时，土地类型、土壤TN含量及pH值也是关键的影响因素.肥料部分氮素的地表流失与降雨量、土壤有机质、TN含量、土地类型及粘粒比重正相关，与施氮量、pH负相关.因此，合理地选择施肥时间和科学地控制施肥数量是减少氮素径流流失的关键.