自20世纪60年代以来，随着对点源污染治理的日益加强，非点源尤其是农业非点源逐步成为水体恶化的主要污染源^[1]。由于非点源污染具有随机性、复杂性、分散性、滞后性和难监测的特点，对其治理存在很大困难^[2]。

随着技术的发展和研究的深入，应用机理模型对非点源污染进行模拟已成为非点源污染研究的常用方法^[3]。目前，常用的非点源污染模型有ANSWERS(Areal Nonpoint Source Watershed Environment Response Simulation)^[4]、SWAT(Soil and Water Assessment Tool)^[5]和AnnAGNPS(Annualized Agricultural Non-point Source Pollution Model)^[6]等。耦合地理信息系统(GIS)与非点源污染机理模型，可以模拟非点源污染形成过程、识别关键源区、分析污染负荷时空分布，相关结果可为制定非点源管理方案提供参考，成为非点源污染研究的重要手段^[7]。AnnAGNPS模型作为连续、分布式机理模型，能够连续模拟指定时间段内每天及累积的径流、泥沙、氮磷等营养物的形成及迁移转化过程，在流域农业非点源污染研究中得到广泛应用。Suttles等^[8]应用AnnAGNPS模型模拟估算了南乔治亚岛上游流域的泥沙量和非点源污染负荷；Polyakov等^[9]利用该模型对夏威夷Kauai岛流域径流量和泥沙侵蚀进行了模拟研究。该模型在美国的Neal Smith国家野生动物避难所^[10]、澳大利亚的Currency Creek流域^[11]、西班牙Olive Orchard小流域^[12]、哥伦比亚的Cauca流域^[13]也得到了应用。国内学者则在三峡库区^[14]、丹江口库区^[15]、太湖流域^[16-17]、辽河源头小流域^[18]、黑河流域^[19]、胶东半岛大沽河流域^[20]等地应用AnnAGNPS进行了研究。东北平原是我国粮食主产区之一，农业面源污染严重，污染的特征、机制和影响因素等问题亟待深入研究。然而，目前还没有AnnAGNPS在东北平原区应用情况的报道。

本研究采用AnnAGNPS模型对苇子沟流域的农业非点源污染负荷进行了模拟，通过实地监测数据对模拟结果进行了校准，在此基础上，分析了苇子沟流域农业非点源污染负荷的时空变化特征。研究结果将对AnnAGNPS模型在我国东北平原小流域的应用提供参考。

1 材料与方法 1.1 研究区概况

研究区域位于饮马河下游，属长春市农安县境内的支流苇子沟小流域(图 1)，上起源头高家店镇高家店村，下至万金塔乡苇子沟村，全长约15 km，地理范围为44°43′~44°37′N，125°26′~125°14′E，流域面积69 km²。苇子沟小流域年平均气温4.7 ℃，无霜期145 d，年均降水量507.7 mm，多集中在6—8月，四季分明，属中温带大陆性气候。该地区位于松辽平原，地势平坦，一般海拔250~270 m，土地利用方式以旱作耕地为主，主要种植玉米，兼少量西瓜和葡萄。

苇子沟流经区域包括高家店、德胜唐、九德号、万金塔、苇子沟等14个行政村，经调查统计，该区域2015年人口为14 490人。

1.2 试验设计与方法

采用实地地表径流监测方法，在苇子沟流域出口处利用ICSO-6712便携式自动水质采样器同步监测每场次降雨过程的地表径流量和氮磷输出浓度。样品的采集频率为降雨前采集一次水样，在降雨开始后的前2 h内每30 min采集水样，此后每1 h采集水样，直至出口径流量基本恢复正常。每次采样1000 mL，采集的水样24 h内送到实验室检测总氮、总磷。所有水质指标均按照标准方法^[21]进行测定。采样器每5 min自动记录流速、水位。场次降雨径流的污染负荷(L)计算公式^[22]为：

式中：C_j和C_j+1为某场降雨径流中第j次和第j+1次监测的污染物浓度，mg·L^-1；q_j和q_j+1为某场降雨径流中第j次和第j+1次监测的径流量，m³·s^-1；Δt为相邻两次采样间隔时间，s；n为某场降雨径流过程中的采样次数。

1.3 AnnAGNPS模型介绍

AnnAGNPS模型由美国农业部自然资源保护局和农业保护局于20世纪90年代在AGNPS模型的基础上研发，是用于模拟评估流域地表径流、泥沙侵蚀和氮磷等营养物流失的连续分布式参数模型^[23]。AnnAGNPS模型由数据输入和编辑模块、年污染物负荷计算模块、数据输出和显示模块3部分组成。

AnnAGNPS作为分布式物理模型包括水文子模型、土壤侵蚀子模型和污染物迁移子模型三个子模型^[24]。地表径流量的计算依据美国土壤保持局提出的SCS曲线法(Soil Conservation Service Curve Number)。峰值径流量用Smith和Willam得出的经验相关方程计算；在土壤侵蚀子模型中，通过修正的通用土壤流失方程(RUSLE)计算各集水单元泥沙侵蚀量，由水文几何通用土壤流失方程(HUSLE)计算各集水单元的侵蚀物向沟道的输移率。污染物迁移子模型则采用CREAM模型和饲育场模型模拟计算溶解态和吸附态的氮磷等污染物质的迁移转化过程。

1.4 模型数据库构建 1.4.1 地形数据

采用经资源三号卫星图像遥感解译生成的数字高程模型(DEM)为数据源，并利用地形参数模块(TOPAGNPS)自动划分集水单元(Cell)和排水沟道(Reach)，生成集水单元文件(AnnAGNPS_cell.csv)和沟道参数文件(AnnAGNPS_reach.csv)。集水单元数据包括排水面积、高程、地面坡度、方位、LS因子等，沟道数据包括沟道长度、高程、坡度等。集水单元的划分和沟道网络的提取通过临界源面积(CSA)和最小沟道长度(MSCL)来确定。钟科元等^[25]的研究表明合理的集水单元划分影响模型模拟的精度。本研究根据苇子沟流域的实际地表状况，将CSA和MSCL分别取8 hm²和130 m，并依据此数值确定了740个集水单元和300个排水沟道数据(图 1)。

1.4.2 气象数据

模型所需的气象参数包括日均降雨量、日最高和最低气温、露点温度、日辐射值和风速。气象数据采用中国气象数据网中国地面气候资料日值数据集和中国辐射日值数据集，采用农安气象站点2009—2015年的逐日气象数据(降水、日最高气温、日最低气温、日均气温、风速、日照时数、相对湿度)，其中日辐射值的计算方法参考童成立等^[26]的研究。采用Wischmeier经验法^[27]计算降雨侵蚀力因子(R)，十年一遇降雨侵蚀因子(EI₁₀)的计算则采用Renard等^[28]的方法。露点温度根据日均气温和相对湿度计算^[29]，公式如下：

式中：A=0.198 0，B=0.001 7，C=0.840 0，均为系数；t_d为露点温度，℃；t为日均气温，℃；U为相对湿度，%。

1.4.3 土壤数据

模型需要的土壤数据分为空间数据和属性数据两类。空间数据采用世界土壤数据库中的中国土壤数据集，数据分辨率为1:100万。流域内土壤类型为石灰性黑钙土和潜育褐色土(图 1)。用GPS定位采集18个土壤样品，分析各类型土壤的pH值、总氮、总磷和有机碳含量。土壤有机氮、无机氮以及有机磷、无机磷的含量分别按照总氮的96%、4%和总磷的15%、85%的比例计算^[30]。根据对数正态分布方法^[18]将土壤粒径的FAO90制转换为美国制标准。利用土壤水分特性软件SPAW(Soil Plant Atmosphere Water)确定各类型土壤的容重、田间持水量、萎蔫系数、饱和导水率和饱和度。

1.4.4 土地利用数据

模型所需土地利用数据采用中国生态系统评估数据库提供的全国土地利用图(分辨率为90 m×90 m)。研究区主要土地利用类型为旱地，占82.61%，居住地占10.96%，林地占5.51%，其他占0.92%(图 1)。

1.4.5 作物及管理数据

流域主要种植玉米，少数村民种植西瓜和葡萄，考虑玉米种植规模，将作物概括为玉米。通过实地调查获得了农作物的种植、施肥、杀虫、收割等相关措施的操作时间和操作方式，以及肥料的养分含量等属性信息。作物参数主要参考模型自带的作物参数数据库(Crop.xls)和美国农业手册703号。

2 结果与讨论 2.1 模型校准

根据国内外对AnnAGNPS模型参数敏感性的研究结果^[31-32]，在水文模块校准时参考田耀武等^[33]的方法反复调试径流曲线数CN值，使径流量模拟值与实测值最大程度地接近。氮、磷输出量主要受土壤可侵蚀因子K、水土保持因子P、作物管理因子C、化肥施用量、土壤中氮、磷本底值和河道曼宁系数的影响，调整上述敏感因子，使得总氮、总磷模拟值最大限度接近其实测值。表 1为2016年6—7月的两次降雨事件模拟实验结果，其模拟偏差计算公式如下：

式中：D_V为模拟偏差(D_V值越趋向于0，则拟合度越好)；V为模型模拟值；V′为实测值。

由表 1可知，AnnAGNPS对小型降雨的模拟结果偏低，与Lenzi等^[34]、邹桂红等^[35]的研究结论一致。可能原因为研究区为平原，平均坡度较低，且土地利用类型主要为旱田，土壤含砂量大，土壤渗透性较大，所以小降雨事件时产生的地表径流值偏低。两场降雨事件的径流模拟结果与实测结果拟合度一般，可能与流域类型为平原有关。两场降雨事件总氮的模拟偏差分别为16.16%、15.48%，平均模拟偏差为15.82%，表明模型对总氮模拟效果较好，能较真实反映流域内总氮的变化情况。总磷的模拟结果较差，其中7月25日场次降雨事件总磷模拟结果与实测结果偏差达到53.74%。已有研究表明，日尺度的污染负荷模拟精度低于月、年尺度^[35-36]。此外，部分研究中，场次降雨事件总磷偏差大，如在黄金良等^[37]应用AnnAGNPS模型对九龙江小流域的氮磷负荷模拟中，有三次降雨事件总磷负荷的模拟偏差分别达到了111.04%、-80.19%和-76.36%；朱乾德等^[38]的研究中有两次降雨事件总磷负荷的模拟偏差分别达到40.09%和39.09%。因此，对于本研究而言，需要更多降雨事件的监测数据来验证总磷模拟的可靠性。总体而言，AnnAGNPS模型对总氮负荷的模拟精度高于对总磷负荷的模拟精度，与Novotny等^[39]对众多非点源污染的模型评估结果一致。可能原因在于，模型对磷元素的模拟只考虑了吸附态磷和溶解态磷，实际上磷在自然界中的循环属于沉积型循环，部分磷进入水体后会沉积在底泥中，而模型忽略了沉积部分的磷，故对磷的模拟值通常大于实测值。

2.2 模拟结果分析 2.2.1 氮磷污染负荷年内分布特点

本研究应用校准后的模型模拟了苇子沟流域2013—2014年的总氮、总磷负荷(图 2和图 3)。从年际变化来看，2013年总氮和总磷的年内变化与降雨量的年内变化一致。由图 2可见：总氮和总磷负荷的月流失量随月降雨量的变化而变化，呈现较强的正相关性。2013年总氮负荷量为26 377.63 kg，总磷负荷量为8 907.70 kg，总氮、总磷月产生强度与月降雨量Spearman相关系数分别是0.85(P=0.000 40)、0.84(P=0.000 40)。5—8月的降雨量分别为43.5、91.3、173.0、105.3 mm，占年总量(562.1 mm)的73.49%，对应期间的总氮负荷量分别为1 204.81、2 412.00、14 691.72、4 231.46 kg，占年总量的85.45%；总磷负荷量分别为382.95、967.00、4 880.93、1 401.80 kg，占年总量的85.69%。由图 3可知，2014年Spearman相关系数分别是0.71(P=0.009 8)、0.69(P=0.009 8)，5—8月的降雨量分别为92.7、73.6、52.2、43.9 mm，占年总量(331.9 mm)的79.06%；对应期间的总氮负荷量分别为730.15、255.06、965.70、126.52 kg，占年总量的84.24%；总磷负荷量分别为241.89、83.40、312.00、39.26 kg，占年总量的84.46%。这说明降雨是苇子沟流域非点源污染的主要驱动力。2014年7月降雨量较5、6月低，总氮、总磷负荷却高于5、6月，可能是因为7月场次降雨事件降雨量更大导致。

2.2.2 氮、磷污染负荷的年际分布特点

根据模型输出文件得到苇子沟流域2009—2015年年际污染负荷情况(图 4)。结果表明，该地区多年平均降雨量为482.2 mm，年降雨量最大和最小年份分别是2010和2014年，降雨量分别为675.2 mm和331.9 mm。总氮、总磷污染负荷年际变化趋势一致，2010年总氮、总磷负荷量分别为79 015.29、25 018.02 kg；2014年则分别为2 466.16、801.01 kg。总氮、总磷多年平均负荷量分别为22 295.28、7 085.00 kg。总氮、总磷负荷随该地区降雨的年际变化而变化，主要原因一是降雨是氮磷流失的主要驱动因素，二是模型对不同年份进行模拟时只考虑气象数据随时间的变化，不考虑其他信息，如地形、土地利用方式、土壤等因素随时间的变化情况^[40]。经SPSS软件进行相关性检验，年降雨量与总氮负荷的Spearman相关系数为0.93(P=0.002 5)，与总磷负荷的相关系数为0.92(P=0.002 5)，均已达到极显著水平，表明流域氮磷污染负荷变化与年降雨量变化的紧密关系。

2.2.3 氮、磷污染负荷的空间分布特征

利用ArcGIS将模型输出的氮磷负荷数据库文件和集水单元图AnnAGNPS_SubWta.asc耦合，输出苇子沟流域2009—2015年总氮、总磷单位面积年均负荷量空间分布图(图 5)。结果表明，苇子沟流域总氮、总磷污染负荷在空间上呈现较强的区域性分布，总氮、总磷年均负荷量低(总氮0~3.00 kg·hm^-2·a^-1，总磷0~1.20 kg·hm^-2·a^-1)的集水单元主要分布在流域西北部，而总氮、总磷年均负荷量高(总氮5.00~13.06 kg·hm^-2·a^-1，总磷2.00~5.52 kg·hm^-2·a^-1)的集水单元主要分布在流域中下游。对比研究区河网图和土地利用图发现，总氮、总磷的流失主要集中在坡耕地和河道两侧等坡度较大区域，表明地形是影响总氮、总磷流失的重要因素。

总氮、总磷负荷在空间分布上有一定相似性，磷流失多的区域氮流失量往往也大，磷流失少的区域氮流失量通常也比较小。这与边金云等^[32]和李开明等^[41]的研究结果一致，可能原因是苇子沟流域营养物以颗粒态为主要赋存形式，而颗粒态营养物质往往附着在泥沙上，随着泥沙的迁移而迁移。

2.3 模拟的不确定性

本研究中AnnAGNPS模型模拟还存在一定的不确定性，主要表现为：

(1) AnnAGNPS假设每年肥料施用和农事管理措施不变，而实际中苇子沟流域农业为农户分散经营模式，耕作施肥方式不一，模型的假设不完全符合研究区域的实际情况。已有研究表明，肥料施用量以及肥料中氮磷比例等参数会影响模型的输出结果^[42]。

(2) 因野外监测和采样的条件限制，供模型校验的实测数据有限，使得模型的校验，特别是对总磷负荷模拟的校验不太理想。

(3) 降雨与径流、氮磷流失负荷的非线性关系，以及土地利用模式和管理措施的快速变化导致了定量化研究非点源污染的复杂性。

(4) 对各年污染负荷进行模拟时均使用了2010年的土地利用图，而现实中研究区域土地利用会随时间发生变化。通过对比2000年和2010年全国土地利用图发现，苇子沟流域土地利用在这10年间变化极小，而且研究区域为东北农村区域，受城市化影响很小。因此，本研究假设该区域2009—2015年土地利用相比2010年保持不变。

3 结论

(1) AnnAGNPS模型对总氮污染负荷量的模拟精度较高，能较好地反映该流域的总氮污染负荷实际情况，而对总磷污染负荷的模拟不确定性较大。

(2) 对2013、2014年污染负荷模拟结果表明，总氮、总磷流失趋势一致，且主要集中在5—8月，与降雨量呈较强的正相关关系。2009—2015年苇子沟流域总氮、总磷负荷年际变化分析表明，在年际尺度上，降雨总量是影响全年流域污染负荷量的重要因素。

(3) 总氮、总磷在空间分布上具有一定的相似性。流域西北部总氮、总磷单位面积负荷量总体较低，而流域中下游的总氮、总磷单位面积负荷量总体较高。旱地是氮磷输出的最主要来源。

(4) 模型对苇子沟流域非点源污染输出模拟具有不确定性，包括模型自身设计的局限性、输入参数的不确定性、模型的适用性以及校验数据有限。

(5) 总体而言，AnnAGNPS模型在苇子沟流域具有一定的适用性。本研究能够为AnnAGNPS模型在东北平原区小流域的应用提供参考。为了增加模型模拟的可靠性，需要更多的实测数据对模型进行进一步校准与验证。