张水锋，张金池，庄家尧，王新猛，张思玉

(1.南京林业大学，南方现代林业协同创新中心，江苏省水土保持与生态修复重点实验室，南京林业大学林学院，江苏 南京 210037；2.南京森林警察学院信息技术学院，江苏 南京 210023)

非点源(NPS)污染已经成为全球性的环境问题之一[1-7]，在我国长江流域、淮河流域、太湖流域等，非点源污染一直是困扰流域管理的关键问题[8-12]。我国人均耕地相对较少[13]，根据2018年耕地面积测算,平均施肥量约419 kg/hm2，已远超世界平均水平120 kg/hm2[14-16]，种植高产品种过程中大量施用化肥农药导致的环境污染尤其是农业非点源污染已成为我国重要的污染源之一[17]。此外，《长江流域水土保持公告(2018年)》显示，长江流域水土流失面积3.47×105km2，占流域总面积的19.36%，流域年输沙量占全国的17%(8.31×107t)[18-19]。因此，我国流域的非点源污染与水土流失问题也不容忽视。

农业非点源污染物模型“AnnAGNPS v1.00”支持以天为基础的连续模拟和改进的处理方法[20]，用于预测径流、峰值流量、沉积物负荷和养分负荷。AnnAGNPS在国外已被广泛用于不同类型、尺度的流域研究，成功模拟了水文、沉积物和养分输送[2-6]，根据粒径级别计算的沉积物产量，以及营养(氮、磷和有机碳)浓度，可用于评估流域对农业管理的响应实践。其中，Chahor等[6]应用AnnAGNPS模型对地中海农业小流域的8个参数进行了敏感性评价，其标定误差小于1%，验证误差小于7%。边金云[21]采用修正的摩尔斯分类筛选法模拟并讨论了AnnAGNPS模型参数的敏感性；席庆[22]将AnnAGNPS模型应用在田中河流域，并采用扰动分析法进行了水文水质参数敏感性分析。多位学者选择采用差分灵敏度分析(DSA)方法对AnnAGNPS模型参数进行敏感性分析[23-27]。参数敏感性分析对模型运行具有重要意义，可以使参与流域模拟的重要参数能够有针对性地得到校准和改进，从而提高模型预测的准确性,已广泛应用于SWAT[27]、HSPF[28]等水文模型中。然而，应用AnnAGNPS模型从敏感性参数分析与适用性评价角度对长三角农业小流域的模拟研究较少，尤其对长江中下游水阳江水系控制区域以传统农业生产为主的小流域研究则更少见。

我国一直将以小流域为单元治理水土流失、改善生态环境、发展经济作为水土保持工作中一项基本策略[29]。本次研究区“沛桥河小流域”属于长江中下游水阳江水系，是典型的以传统农业为主的小流域，位于长三角地区的重要农业基地江苏南京高淳区内，且存在非点源污染等管理难题。随着近年来我国测土配方施肥以及农用地整治措施的逐步推进，区内农业作物种植造成的非点源污染和水土流失等问题得以逐步缓解，为此，本研究拟通过AnnAGNPS模型在沛桥河小流域的适用性研究，以期为长三角地区的流域综合管理工作提供数据支撑和科学依据，辅助流域生态-经济复合系统可持续发展的多目标决策。

1 材料与方法

1.1 研究区概况

沛桥河小流域位于江苏高淳(118°41′～119°21′E，31°13′～31°26′N)，高程为4～179 m。沛桥河小流域总面积34.85 km2，地处茅山山脉的余脉，山势平缓，为长三角区域典型的低山丘陵区，属于水阳江水系，土壤类型主要包括粗骨土、黄棕壤、红壤、石灰(岩)土、渗育型水稻土和潴育水稻土。水田和旱地各占地约37%和25%，为典型的以农业用地为主的小流域，属于北亚热带和中亚热带过渡季风气候区。受季风环流影响，四季分明，光照充足，无霜期长，每月平均气温下降5.6 ℃，2005—2018年多年平均降雨量1 296.1 mm，降雨主要集中在6—9月，占全年降雨的52.82%，且雨强较大，常发生短历时暴雨，其中7月降雨量最大，1—3月与11—12月降雨量较小。高淳区是华东长三角地区的重要农业基地[30]，截至2017年末，农用地总面积580.93 km2，占全区总土地面积的73.51%，粮食作物种植主要以小麦和水稻为主，占当季粮食作物种植面积比例的87%以上，而经济作物主要以油菜为主，2016—2018年年均种植达60×103hm2，面积大小位于小麦和水稻之间。该区也是江苏省第一个国家级生态区，农业快速发展的同时，作物种植导致的非点源污染与水土流失问题也不容忽视。

1.2 数据来源与预处理

AnnAGNPS模型对流域过程的模拟需要大量的流域数据为支撑，指标包括地形数据、土壤数据、土地利用数据、气象数据和土地管理数据等5大类。模型模拟精度与所获取指标数据的确定性和准确性紧密相关，为使参数指标能更好地反映流域特征且保证模拟的精度，本研究做了大量的实地和资料调研，确定各参数指标的数据可得性和可用性。

数字高程模型(DEM，分辨率30 m)、土地管理数据(包括研究区农用地种植概况、农用地利用方式，以及农作物耕作和施肥管理)由高淳区自然资源部门提供。土地利用数据则基于南京市自然资源部门提供的2016年沛桥河小流域1 m×1 m分辨率原始遥感影像，通过ArcGIS、ENVI与BigMap软件进行协同监督分类，获取沛桥河小流域土地利用类型图。气象数据包括7类逐日数据：24 h降雨量值、气温值(最高、最低、平均和露点)、蒸散发、风速值、地表土壤反射率，以及云量和太阳辐射量值等，主要来源于南京市气象部门提供的高淳气象观测站(118°54′E，31°20′N)资料。其中，太阳辐射数据根据高淳气象观测站提供的监测点经纬度和日值日照时数等数据，利用Angtrom-Prescott方程[31]计算获取；云量计算是根据金波尔提出的云量与总辐射经验关系式[32]；地表土壤反射率采取闫胜军[33]的经验公式。

模型所需的土壤理化数据包括土壤质地、土壤比质量、饱和导水率、田间持水量、凋萎点、土壤可侵蚀因子K值、土壤水文组类别，以及有机碳、氮、磷含量与pH等，主要通过实地采样与实验室分析获取。采样主要通过小流域内不同土壤类型区域坐标，在各区确定不同土地利用类型中的3个重复样地，按照“S”形取样方法分3层各取土壤样品6个，共取土壤样品270个、环刀156个，带回实验室分析。土壤质地根据农业行业土壤检测标准NY/T 1121.3—2006利用比重计法测定，土壤容重、饱和导水率、田间持水量的测定采用环刀法，土壤有机碳含量采用德国Elementar Macro元素分析仪测定，土壤pH采用IQ150土壤原位pH计测定，土壤N、P、K等养分含量通过TPY-6A土壤养分速测仪测定，用于计算流域输出养分负荷总量。凋萎点采用SPAW软件中的Soil Water Characteristics模块[33]进行计算，土壤可蚀性因子K值则根据Williams等[34]提出的EPIC(Erosion-Productivity Impact Calculator)模型中的土壤可蚀性K值的估算方法。参照美国农业部自然资源保护局(Natural Resources Conservation Service，NRCS)通过的土壤最小下渗率的范围将土壤水文类型组划分为A、B、C和D等4类[35]，其中的最小下渗率计算采用车振海[36]设计的经验公式。

1.3 最佳水文响应单元划分

AnnAGNPS模型通过临界源区面积(critical source area，CSA，式中记为ACS)与最小源区沟道长度(minimum source channel length，MSCL，式中记为LMSC)的取值来划分水文响应单元(cells)和自动切取流域河段(reaches)水系，计算子流域的面积、坡度与坡长等各项特征参数。经测试，CSA和MSCL的取值越小，流域河段越密集，水文响应单元数量也越多，对地表的概化程度越详细。然而并不是越小越详细效果越好，以过小的响应单元划分出的沟道并不是真正的沟道，即容易失真[37]，同时AnnAGNPS的计算工作量和时间会大量增加。CSA与MSCL的取值主要依据研究区的地形和土地利用的特点来确定，需尽可能以较少的单元反映实际的地形、土地覆被等，从而保证模拟的精度。如表1所示，本研究参照了国内[22-25, 33]和国外[38-39]对其取值的研究，结合研究区的面积大小与地形特征，将本研究测试的取值范围设置为5≤ACS≤50(增量取5)，20≤LMSC≤120(增量取10)，分别对取值组合进行测试。

表1 CSA与MSCL取值参考

1.4 参数敏感性分析方法

径流曲线数(CN)作为一个能综合反映降雨前流域特征的综合参数，其值取决于流域土壤类型、植被覆盖类型、水文状况、土壤湿润状况以及农业耕作措施等。依据以往的研究经验[21-26,38-41]，CN值被一致验证为影响流域径流、泥沙，以及氮与磷产量最敏感的参数，且可通过AnnAGNPS模型提供的技术手册(TR-55)对研究区的CN值进行人工取值率定，达到模型校准的目的。

由于本次模拟涉及的相关土壤理化基础数据均采用了野外取样与实验室测定的方式获取，并非通过查询全国土壤普查数据获取，模型中土壤理化数据相关的输入参数相对准确。因此，本研究参数敏感性分析以“径流曲线数(CN)”为主，主要用于率定AnnAGNPS模型。同时，由于模型输入参数较多，参考以往研究[21-26,38-41]，遴选了与土壤、作物管理相关且可通过人工生态修复方式改善的10个参数(土壤可蚀性因子、田间持水量、饱和导水率、凋萎点与土壤有机质含量，以及作物冠层覆盖度、作物留茬30%、肥料施用量、作物管理因子与水土保持因子)进行敏感性分析。

本研究采用DSA方法进行参数敏感性分析，DSA是一个有限差分数值逼近方法[40]。如式(1)所示，y0作为x0输入模型的输出值，在保持其他参数值不变的情况下，其中的某个参数从初始值以±Δx为步长上下增减，x1=x0-Δx，x2=x0+Δx对应的输出值y1和y2。指数符号表示模型对输入参数的变化是否有响应，即参数的增加是否导致输出变量的增加，参数的减少是否导致变量的减少，或者效果相反。DSA方法根据参数的敏感性指数取绝对值|S|，将其敏感性程度划分为4个等级：可忽略(Ⅰ为0.00～0.05)、中等(Ⅱ为≥0.05～0.20)、高(Ⅲ为≥0.20～1.00)、非常高(Ⅳ为≥1.00)。

(1)

式中：S为敏感性指数；x0为参数初始值；y0、y1、y2为参数x0、x1、x2情况下对应的模型输出值；Δx为变化步长，10%×x0。

1.5 模型率定及评价方法

根据以上的敏感参数分析方法，获取对模型模拟结果最敏感的参数，采用试错法人工调参校准模型[41]。针对模拟结果与观测数据进行对比分析，以2008—2014年作为率定期，2015—2018年作为验证期。采用Legates等[42]在应用AnnAGNPS模型评估流域泥沙和氮磷负荷中的建议，通过相关系数(R2)、效率系数(E)以及常用的相对误差(ER)来综合评价模型的适用性，对每个评价指标的局限性和适宜性做进一步的阐述：相关系数(R2)的取值范围为0.0(模拟精度低，不良模型)至1.0(模拟精度高，完善模型)，如式(2)所示；效率系数(E)如式(3)所示，它的范围从负无穷(模拟精度低，不良模型)到1.0(模拟精度高，完善模型)；相对误差(ER)如式(4)所示。

R2=

(2)

(3)

(4)

式中：Oi为观测数据；Pi为相应的模拟数据；O为观测数据集的均值；P为模拟数据集的均值；i为降雨与径流事件的编号；N为数据对的总数。

Moriasi等[43]与Saleh等[44]通过对常用水文模型的模拟应用与对比分析，认为R2和E大于0.5时模拟效果尚可，大于0.65则表示模拟的精度很高；ER取值可为正，亦可为负，在±40%之内为可接受的范围，若其绝对值小于15%则可满足研究需要[21-24]。

2 结果与分析

2.1 最佳水文响应单元划分分析

流域地形越平坦，土地利用差异越小，则CSA与MSCL的取值越大。鉴于研究区相对地势平坦，土地利用类型以农用地为主，也考虑到模型运算过程的时效性，CSA和MSCL不能取值太小。但其取值也不能过大，过大的取值会使划分结果概化程度大，无法反映小流域下垫面的真实利用场景属性，因此，根据设置的取值范围分别取值进行组合测试。结果表明，随着CSA值和MSCL值增大，地块单元越大，划分结果概化程度越大。根据小流域土地利用类型图中林地、草地、耕地的面积大小，经多次分组测试后，当CSA取值5、10或20，MSCL取值50或100时(共6个组合，如图1所示)，最小水文响应单元的划分结果更符合研究区的土地利用现状特征。

图1 不同CSA与MSCL取值组合的最小水文响应 单元离散分布图Fig.1 Discrete distribution of cells with different values of CSA and MSCL

根据土地利用类型图和DEM的分辨率水平，并结合野外实地调研对土地利用现状特征的深入理解，当取值组合为“ACS=20，LMSC=50”时，沛桥河小流域的117个水文响应单元和50个河段，更符合沛桥河小流域下垫面的地块和河段特征，其具体分布如图2所示，各断面之间的河段长度为42.43～3 137.94 m。

图2 沛桥河小流域的水文响应单元分布图(A)和 河道网络分布图(B)Fig.2 Cells distribution of Peiqiao River watershed(A) and channel network distribution(B)

2.2 参数敏感性分析

参数敏感性分析通过将“径流曲线数(CN)”以及遴选的其他10个参数对径流、泥沙、总氮与总磷输出的变化率代入式(1)，计算结果如表2所示。由表2可以看出，径流、泥沙、总氮与总磷负荷的输出均对CN的敏感性相对最高，且呈明显正相关关系，泥沙输出的敏感等级为Ⅲ，而对径流、总氮与总磷负荷输出影响的敏感程度达到Ⅳ级。除CN外，径流深仅对凋萎点、田间持水量和作物冠层覆盖度的变化表现敏感，且等级为Ⅲ，而对其他参数的敏感性可忽略；泥沙对土壤可蚀性因子(K)、水土保持因子(P)与作物管理因子(C)的敏感性等级均为Ⅲ，对其他参数则不敏感；总氮、总磷负荷对土壤有机质含量的敏感等级均为Ⅲ，总磷对田间持水量、肥料施用量的敏感性等级也达到Ⅲ；泥沙负荷对凋萎点不敏感，而总氮、总磷负荷却表现为负相关的中等敏感(Ⅱ)。

表2 所选参数敏感性指数值与等级

2.3 模型的率定与验证

2.3.1 径流、泥沙率定与验证

1)基流计算。基流长期补充河川的径流，是水文计算中不可忽略的重要部分。本研究从方法的客观程度和基流特征出发，选择应用广泛且效果较好的数字滤波法[45]中单参数数字滤波方程[46](式5)：

(5)

式中：qt，t时刻地表径流量，m3；Qt，t时刻实测径流量,m3；bt，t时刻的基流量,m3；t，时间，d；β，滤波参数，一般取值为0.900～0.950[47]。

本研究中β取值为0.900、0.920、0.935、0.950，进行4次过滤，将2008—2018年的观测径流和AnnAGNPS径流模拟数据代入式5，获得2008—2018年逐年的基流占总径流比例(BFI)值如表3所示。通过不同β取值计算得到多年平均BFI分别是0.434、0.407、0.391和0.367。其中，当滤波参数为 0.935，计算后基流量比前两次取值分别减少4.31%、6.52%，且反映水文序列变异程度的离散系数Cv分别减少了0.08和0.14，基流相对平稳。因此，本研究的滤波参数β取值0.935。

表3 沛桥河小流域2008—2018年基流占总径流比例计算结果

2)率定与验证。由于长江下游流域基流受降雨影响较大，按照降雨量年平均值把径流时间分为汛期4—9月(降雨量>100 mm)和非汛期10月至次年3月(40 mm<降雨量≤100 mm)两个区间，进行分别率定。通过基流处理后，对月尺度径流、泥沙负荷按汛期和非汛期两个阶段进行了率定与验证。由于CN值是对径流深、泥沙、总氮与总磷的模拟效果敏感性程度均最高的参数，因此本研究主要通过采用试错法人工调参CN值校准模型。模型校准后，林地CN取值范围是30～77、草地CN取值范围是32～82、农用地的CN取值范围是63～91。

以Moriasi等[43]与Saleh等[44]的研究为参考标准，该模型在沛桥河小流域对月尺度的模拟结果可靠，且模拟效果较好。对泥沙按照汛期和非汛期分别进行率定与验证分析，结果如表4所示，汛期和非汛期的相关系数(R2)、效率系数(E)均分别大于0.87和0.83，相对误差(ER)均小于11%，且汛期的率定与验证效果均高于非汛期。因此，模型的径流模拟效果虽略低于泥沙，但整体模拟效果可靠。

表4 沛桥河小流域汛期与非汛期径流、泥沙、总氮、总磷模拟评价结果

验证期(2015—2018年)汛期与非汛期之径流、泥沙的模拟与观测值随时间的分布趋势基本一致(图3)。

基于图3中的径流和泥沙数据进行Pearson分析可知，泥沙与径流随时间变化的趋势在汛期较为相似，两者模拟值之间的Pearson相关系数达到0.916，而观测值之间的Pearson相关系数达到0.894，均在0.01水平(双侧)上显著相关；而非汛期泥沙与径流的模拟值和观测值之间Pearson相关系数分别为0.666、0.701，均小于汛期，但也均在0.01水平(双侧)上显著相关。验证期(2015—2018年)模型模拟的径流输出、泥沙负荷与观测值的相对误差(ER)均值分别为-11.50%、9.53%，因此，模型对泥沙的模拟效果略优于径流，与表4的模拟结果相一致。

图3 2015—2018年沛桥河小流域月径流、泥沙、总氮、总磷的观测值与模拟值对比Fig.3 Comparisons of actual and simulation runoff, sand,TN and TP in Peiqiao River watershed from 2015 to 2018

2.3.2 总氮、总磷率定与验证

以2008—2014年作为率定期，2015—2018年作为验证期，对总氮、总磷负荷进行率定和验证，结果如表4所示：2015—2018年的汛期与非汛期，总氮的相关系数(R2)与效率系数(E)均大于0.80，相对误差(ER)均控制在10%以内；而总磷的相关系数(R2)与效率系数(E)均大于0.78，相对误差(ER)的绝对值均控制在11%以内，其模拟结果的精度总体略小于总氮。此外，非汛期的模拟精度略高于汛期，这不同于径流、泥沙的模拟结果。

验证期2015—2018年AnnAGNPS模型模拟的年度总氮、总磷负荷与其观测值随时间的分布趋势基本一致(图3)。基于图3中总磷和总氮数据，经Pearson相关性分析，总氮、总磷负荷的模拟值与观测值均在0.01水平(双侧)上显著相关。其中，汛期模拟值之间的相关系数达0.933，观测值之间的相关系数为0.897；而非汛期的相关系数小于汛期，分别为0.773和0.789。验证期2015—2018年AnnAGNPS模型模拟的年度总氮、总磷负荷与其观测值相比，相对误差(ER)的年度均值分别为-6.88%和-7.99%，均可满足研究需要。因此，对总氮的模拟效果略优于总磷，与表4的模拟结果相一致。

3 讨 论

选择长三角区域以传统农业为主的沛桥河小流域为研究对象，全面构建了模型基础数据库，尤其针对AnnAGNPS模型模拟所用的关键土壤理化基础数据，采用了野外取样与实验室测定的方式，避免了因查询获取的全国土壤普查数据尺度较大对小流域尺度非点源污染模拟的影响，相对提高了模型模拟的精度。

本研究中，AnnAGNPS模型在沛桥河小流域的最佳水文响应参数CSA与MSCL分别取值20与50为宜，这与文献[23-25, 33]的结果有差异，分析认为是由于各研究区地形的特点不同导致：文献中的研究区主要位于黄土丘陵沟壑区、西南岩溶地区、福建晋江流域，研究区面积均较大(130～2 581 km2)，且平均高程多在800 m以上，地势起伏均较大，相互之间地貌类型差异也较大；而本研究区面积与高程均较小，相对地势平坦，土地利用类型以农用地为主。

此外，由于前人对日尺度径流模拟得出一致的分析结果[45-46]，即基于24 h降雨数据测定的降雨数据通常达不到24 h，相对短暂，日尺度径流的模拟某种程度上延长了降雨时间，因此会削弱对降雨带来的效果，模拟结果一般不够理想，若要提高日尺度的模拟精度，则需要如小时降雨数据等更精确的气象数据。因此，本研究根据所获取的气象数据选择对2008—2018年汛期与非汛期的月尺度泥沙、总氮与总磷负荷进行了率定与验证。其中，AnnAGNPS模型对沛桥河小流域总氮负荷的模拟效果优于对总磷的模拟效果，这与边金云[21]、席庆[22]、梁丽营等[41]的研究结果相一致。

本研究的参数敏感性分析中，径流曲线数(CN)是对径流深、泥沙、总氮与总磷输出模拟总体影响最高的参数，除泥沙的敏感性等级为Ⅲ外，均达到最高级Ⅳ，这与以往的研究[21-26,38-41]相一致。由此可见，模型率定应以径流曲线数(CN)为主，土壤可蚀性因子(K)、水土保持因子(P)、作物管理因子(C)、土壤有机质含量、作物冠层覆盖度、田间持水量、肥料施用量对模型输出的影响较大，应以实验室精确测定为佳，而不能仅依靠中国土壤数据库中查询的大尺度数据。在长三角区域内以传统农业为主的农业小流域综合管理工作中，可通过人工生态修复措施调控这些参数以达到优化小流域生态环境的目标。

综上所述，本研究通过对AnnAGNPS模型的径流、泥沙、总氮与总磷输出模拟进行了参数敏感性分析与适用性评价，结果相对可靠，且模拟效果较好，可为长三角农业小流域非点源污染与水土流失管理决策提供依据。