AnnAGNPS模型(Annualized Agricultural Non-point Source Pollution Model)是由美国农业部研发的广泛应用于农业非点源污染研究的分布式参数模型, 能以日为固定步长对流域一个时间段内每天累计的地表径流、土壤侵蚀、营养盐流失等进行连续模拟^[1-2]。国内外研究表明, 该模型机制复杂, 输入参数类型多样, 参数不容易直接获取, 同一参数应用于不同地区的模拟结果也不尽相同, 参数的不确定性对模拟结果有显著影响^[3]。因此, 如何使参数的取值更加符合实际情况成为模型可靠运行以及模拟结果真实可信的重要前提条件, 参数敏感性分析以及参数识别和率定成为分布式水文水质模型研究及应用的关键。采用定量或定性的方法, 评价模型参数的敏感性及重要性程度, 识别模型主要敏感性参数, 提高模型参数的量化能力, 降低模型参数维度及模拟结果的不确定性, 进而提高模型效率, 便于模型进一步推广应用^[4]。

纵观国内外对AnnAGNPS模型的应用和研究, 在模型机理与预测偏差^[5]、数据库构建^[6]、模型适应性检验^[7]和流域非点源管理措施效果模拟^[8]等方面的研究成果较多, 而专门针对参数敏感性分析的研究成果相对较少或不成系统。SARANGI等^[1]以St. Lucia流域为例, 研究表明径流曲线数CN值是影响径流量的重要敏感性参数, 土壤可蚀性因子K值是影响输沙量的重要敏感性参数。CHAHOR等^[9]以西班牙Latxaga流域为例, 发现RUSLE-P和冠层覆盖度是影响输沙量的重要敏感性参数。AnnAGNPS模型在我国不同地区应用时影响模拟结果的敏感性参数也存在差异性。以南方低山丘陵区为例, 桃溪河流域^[10]、九龙江流域^[7]、马槽河流域^[11]、四岭水库流域^[12]和章溪河流域^[13]等同属亚热带季风气候, 地形以低山丘陵为主, 土壤以黄红壤、水稻土为主, 植被覆盖度较高, 研究发现CN值是影响流域径流的最主要敏感参数, 此外, 土壤有效含水率是影响桃溪河流域径流量的重要敏感性参数; 田间持水量是影响四岭水库流域径流量的重要敏感性参数; 降雨侵蚀力和坡度因子是影响九龙江流域、马槽河流域产沙量的主要敏感性参数; 地表残留覆盖率、年根系生物量和土壤可蚀性因子是桃溪河流域产沙量的主要敏感性参数; 水土保持因子是章溪河流域产沙量的最主要敏感性参数; CN值、地形因子和耕作管理因子是影响四岭水库流域氮磷输出的重要敏感性参数; CN值和水土保持因子是影响章溪河流域氮磷输出的重要敏感性参数。这些研究主要集中在模型应用模拟方面, 而对模型参数敏感性分析缺乏系统性研究。与上述地区相比，黄土丘陵沟壑区气候干燥少雨, 地形支离破碎，沟壑纵横, 土质疏松, 植被覆盖度低, 这些因素导致模型在该地区应用时参数取值及对模型输出结果的敏感性与其他研究区相比会存在较大差异性。而在黄土丘陵沟壑区对模型参数进行系统的敏感性分析鲜有报道。因此基于前人研究, 结合岔口小流域实际情况, 筛选出适合流域敏感性分析的参数, 采用修正的摩尔斯分类筛选法对AnnAGNPS模型参数进行相对系统的敏感性分析, 以期识别出对模型模拟结果影响较大的敏感性参数, 以此作为模型参数优先调试和优化的对象, 尽可能使参数值准确、可靠, 有效缩小模型估计的不确定性范围。

1 研究数据与方法 1.1 研究区概况

岔口小流域位于山西省西南部永和县、石楼县和隰县交界地带, 涉及永和县4个行政村, 隰县和石楼各1个自然村, 总面积131.91 km², 地处北纬36°47'26"~36°57'14", 东经110°38'01"~110°50'02", 平均海拔1 184 m。研究区属于典型的暖温带半干旱大陆性季风气候区, 四季分明, 日照时间长, 多年平均气温8.6 ℃, 多年平均降水530.9 mm, 并多集中在7—9月, 约占全年降雨量的62.25%。流域内河流属季节性河流, 丰水期为历时很短的洪水径流且径流中携带大量泥沙和营养物质。研究区位于典型的黄土丘陵沟壑区, 流域内沟壑纵横, 导致水土流失严重, 生态环境脆弱。流域内共有17种土壤类别, 土壤质地以粉砂质壤土为主, 湿陷性较大。研究区主要土地利用类型为耕地、园地、林地、草地、水域和居民点, 其中林地面积占比较大。流域内产业结构单一, 以农业生产为主, 经济发展迟缓, 无较大城镇，无大型养殖场和工业企业, 因此可以不考虑点源污染。研究区示意图见图 1。

1.2 AnnAGNPS模型构建

模型包括水文、土壤侵蚀和污染物迁移等模块; 基础数据库主要有气象、地形、土壤和土地利用数据库。流域地形数据选用分辨率为10 m的数字高程模型栅格数据。流域土地利用现状图参考第2次全国土地调查土地利用现状数据, 以2013年QuickBird影像为底图进行全流域野外实地调查生成。径流曲线数CN值参考取值条件和流域不同的土地利用方式、土壤类别及植被覆盖等确定; 流域内作物、田间管理和肥料数据主要通过入户调查和实地采样测得; 蓄水坝数据通过岔口小流域坝系工程设计文本和实地调查获得; 土壤数据由永和县土壤1:5万地图经矢量化得到; 土壤理化性质主要通过实地采样获得; 土壤可蚀性K值根据土壤粒径含量并参考诺莫图获得; 气象数据为2001—2013年永和县气象站逐日气象资料数据(露点温度、云覆盖、每日地面太阳辐射等通过经验公式间接获得^[14-15])和流域内设立的3个雨量监测站实测数据。根据流域实际情况, 最终确定临界源面积(CSA)和最小初始沟道长度(MSCL值)分别为5 hm²和100 m, 将流域划分为8 121个子流域(图 2)。将模型所需基础信息数据导入模型输入编辑器, 初步建成岔口小流域AnnAGNPS模型。

1.3 模型敏感性分析参数选取

在参数初步率定的基础上, 结合模型机制结构, 同时兼顾模型操作运算情况, 从可能对模型输出结果有影响的参数中筛选出与SCS-CN径流曲线数模型、RUSLE方程和营养盐输出模拟过程等有关的10个参数(表 1)进行相对系统的敏感性分析。

1.4 参数敏感性分析方法

对模型参数进行相对系统的敏感性分析可识别出对模型模拟结果影响较大的敏感参数, 在此基础上对模型进行优化处理, 使模型在实际应用中更加科学合理。摩尔斯分类筛选法因简单易懂, 计算量较小, 特别适于参数较多的复杂模型的敏感参数识别, 但仅能给出定性结果, 无法给出定量结果^[5]。笔者采用修正的摩尔斯分类筛选法, 假设各参数间相互独立, 选择1个参数为随机变量, 在其阈值范围内按固定比例取值, 并将不同取值代入模型从而得到不同的模拟结果, 通过基准值判断参数值改变对模型输出结果的影响程度^[16]。计算公式为

$ S = \sum\limits_{i = 0}^{n - 1} {\frac{{\left( {Q{}_{i + 1} - {Q_i}} \right)/{Q_0}}}{{\left( {{P_{i + 1}} - {P_i}} \right)/100}}} /\left( {n - 1} \right)。$

(1)

式(1)中, S为敏感性指标因子; Q_i和Q_i+1为第i次和i+1次模型模拟输出值; Q₀为参数率定后模型模拟基准值; P_i和P_i+1为第i次和第i+1次模型运算参数值相对于率定后参数值的变化百分率; n为模型运行次数。

基于模型各参数的阈值和基准值, 在基准值基础上以5%的步长对参数进行扰动, 每个参数值分别±5%、±10%、±15%和±20%, 利用模型的不同输出结果与参数的相对变化值来分析模型参数的敏感性。

2 分析与讨论

分析输入参数对模拟结果的影响, 可更加清晰和有针对性地描述各参数对模型输出结果的不同影响趋势, 为参数调试提供参考。10个输入参数变化对模型径流量、泥沙负荷、氮磷负荷、总有机碳负荷模拟结果的影响如图 3所示。

由图 3可知, CN值与模型模拟输出结果呈明显正相关, 即模拟结果随着CN值的增大而相应的增加。对于模拟结果, 径流量的变化幅度最大, 是由于模型采用SCS-CN径流曲线数模型计算径流量, CN值的大小直接影响模拟径流量的大小。氮磷输出中, 吸附态氮磷较溶解态氮磷变化幅度小是因为CN值通过影响模型水文模块的运算进而影响流域径流量的变化, 溶解态污染物输出随着径流量的变化而改变, 吸附态氮磷输出主要受土壤侵蚀模块影响。饱和导水率(SHC)对模型模拟结果均有不同程度的影响。总磷和总氮随着SHC值的增加而增加, 不同形态的氮磷变化趋势也不同, SHC值的增加会导致溶解态氮减少, 吸附态氮磷、溶解态磷增加，总体而言该参数主要影响吸附态氮。

化肥施用量(FERT)主要影响模型氮磷输出, 对其他模拟结果均无影响。氮磷流失量随着施肥量的增加而增加, 两者基本上呈正比关系。与吸附态氮磷输出相比, 溶解态氮磷对该参数更敏感, 尤其是溶解态氮。降雨侵蚀因子(R)对模型输出结果无影响, 因此可以判定其不是模型的敏感因子, 这一结论与李家科等^[17]研究结果相一致。

土壤可蚀性因子(K)对模型径流量、溶解态氮输出均无影响, 对其他输出结果均有较大影响, 并呈明显的正相关特征。K值是RUSLE方程中计算水土流失的重要因子之一, 是评价土壤对侵蚀影响作用的因子, 代表土壤对侵蚀应力分离和搬运作用的抵抗能力其主要影响流域泥沙输出量。地形因子(LS)对模型径流量、溶解态氮的输出均无影响, 对模型其他输出结果均有较大影响, 并呈明显的正相关特征, 对吸附态氮磷的影响大于溶解态氮磷。流域属于典型的黄土丘陵沟壑区, 流域内植被稀疏，沟壑纵横，土质疏松, 地形坡度越大, 雨水对土壤侵蚀越严重, 流域内泥沙流失也越严重, 从而导致吸附态氮磷输出量增加, 总氮、总磷污染物负荷必然随着坡度的增加而增加。

耕作管理因子(C)和水土保持因子(P)对径流量和溶解态氮输出无影响, 而对其他输出结果有较大影响, 并呈明显的正相关特征。2个参数对泥沙负荷和有机碳负荷影响最大, 其次是吸附态氮磷。田间持水量(FC)对模型模拟结果均有不同程度的影响, 且呈明显负相关特征, 该参数对径流量影响最明显, 其次是有机碳负荷, 再次是泥沙负荷。田间持水量主要通过SCS-CN径流曲线数模型影响地表产流, 随FC值的增加, 径流量和泥沙逐渐减少, 导致氮磷等污染物负荷输出减少。沟道曼宁系数(MA)对模型模拟结果的影响趋势与田间持水量相同, 呈明显负相关特征。流域内河道粗糙度随着沟道曼宁系数的增加而增加, 直接导致径流阻力增大，流速变缓, 使得泥沙输出量减少, 同时也会减少吸附态污染物输出。

岔口小流域属于典型的黄土丘陵沟壑区, 地形、气候、土壤和土地利用等特征明显, 使得在该地区应用时所选取的AnnAGNPS模型参数具有地域性、特殊性和复杂性等特点。研究区坡耕地坡度较大, 土质疏松, 以细沟侵蚀为主; 其他耕作地类坡度较小, 以细沟和细沟间侵蚀为主。因此, 在模型参数设定时, 除考虑流域潜沟侵蚀和层蚀外, 增加了因耕作不同而引起的短暂切沟侵蚀数据。土壤以地带性褐土为主, 结构松散, 湿陷性较大, 抗侵蚀能力差, K值远大于其他地区。降雨历时短且强度高的次降雨是造成流域水土和氮磷流失的主要途径, 次降雨中E_I30与径流和泥沙均呈较好的指数关系, 且拟合程度高^[18]; 流域施肥类型多样, 以复合肥为主, 有机肥较少, 可忽略不计。调查发现在7月农民以播撒方式对农作物进行追肥, 此时处于雨季, 导致肥料不能及时被农作物吸收而随径流和泥沙进入环境中造成污染, 因此模拟流域农业非点源污染时应重点考虑施肥量及方式的影响。流域内共有130座大中小型淤地坝, 模型淤地坝模块设计上具有局限性, 只考虑淤地坝对流域泥沙淤积的影响, 未考虑坝内死水位对径流的扣留和超过死水位时由卧管排出对泥沙和径流流速的影响。结合流域实际情况和开发人员Ron Bingner的建议, 将wetland和impound数据一起使用且共用1个沟道编号(reach id)模拟淤地坝对流域产流输沙及污染物的影响, 效果较好。

综上所述, 由于流域内各种污染物主要以径流冲刷和土壤侵蚀为输出载体, 所以模型模拟结果的准确性与径流和土壤侵蚀密切相关, K值和LS值在一定时间内发生变化的可能性极小。因此, 控制和治理岔口小流域水土流失及非点源污染时应着重考虑径流曲线数、植被覆盖、耕作管理方式合理性和水土保持情况。

3 计算结果

利用修正的摩尔斯分类筛选法分别计算10个参数的敏感性指标因子(S值), 结果如表 2所示。

由表 2可知, 对模型模拟输出结果最为敏感的参数是CN值, 其S值远远高于其他参数对模型模拟结果的影响。

对径流量影响最大的是CN值, S值为10.30;其次是FC值, S值为-5.44;再次是MA值, S值为-0.53, 其他参数对径流影响较弱或者无影响。对泥沙量影响最大的是CN值, S值为7.12;其次是K值, S值为1.7, 再次是C值, S值为1.67。此外, P值和LS值对泥沙输出量也相当敏感。这与RUSLE方程中产沙量主要影响因子是吻合的。

对总氮和总磷影响最大的是CN值, S值分别为5.35和5.02;其次是K值, S值分别为1.41和1.60；再次是C值, S值分别为1.32和1.58。对泥沙负荷和氮磷负荷影响最敏感的参数主要是CN值、K值、P值、C值和LS值。泥沙负荷与氮磷负荷的影响因子较一致, 说明氮磷负荷主要取决于泥沙的流失。在CN值保持不变的条件下, 对吸附态氮、磷和溶解态磷最敏感的都是C值, S值分别为1.68、1.58和1.57。根据流域2007—2013年氮磷污染物负荷空间分布(图 4)可知, 氮磷流失空间分布并不具有一致性, 但区域性较强, 氮流失主要发生在沟谷两侧, 这与径流冲刷有关; 磷流失主要发生在植被覆盖度较低的地区, 这与泥沙侵蚀有关, 说明CN值对总氮的影响大于总磷, 产沙因子对总磷的影响均大于总氮。

对总有机碳影响最大的是CN值, S值为7.70;其次是FC值, S值为-1.80;再次是P值, S值为1.67。土壤有机碳流失主要以泥沙为承载体被带走, 土壤侵蚀造成了有机碳在泥沙中富集, 且富集比大于1^[19]。因此，对泥沙敏感程度高的参数对总有机碳的敏感度也高。对各所选参数敏感性进行排序, 结果如表 3所示。

由于流域地理位置的不同, 气候、地形、土壤、植被和土地利用等具有显著的地域性和空间异质性, 导致同一参数在不同流域敏感性也存在显著差异。R对模型输出结果无影响, 属于非敏感参数, 而在九龙江流域^[7]、马槽河流域^[11]、大沽河流域^[20]则是影响泥沙、氮磷的较高敏感因子, 原因可能在于黄土丘陵沟壑区水土流失受场次降雨强度影响较大, 受日降雨量影响小。贾志军等^[21]、王万忠^[22]认为在黄土高原E_I10或E_I30为R值的最佳表达式, 基于日降雨数据利用章文波等^[23]提出的计算公式得出的降雨侵蚀力为半月侵蚀力, 因此可能导致R值在研究区为非敏感参数。李海防等^[24]认为水土流失是一个受多种因素影响的复杂过程, 降雨侵蚀力不能独立反映土壤侵蚀过程, 要想准确描述土壤侵蚀过程, 还需综合考虑流域土地利用、植被覆盖、坡度、坡长、土壤可蚀性和水土保持措施。与其他研究区相比, 岔口流域化肥施用量只对溶解态氮高度敏感, 对其他输出结果敏感性较低或无影响, 而这个参数在中田河流域^[25]、九龙江流域^[7]则是影响氮磷的高敏感参数, 这主要与各研究区作物种类、施肥量及种类不同有关。对岔口小流域, 农民为了单纯追求作物产量, 施用含无机氮比例较高的复合肥, 导致对溶解态氮影响较大。沟道曼宁系数对流域径流影响较大, 而在其他研究区对径流并无影响, 这可能是因为沟道曼宁系数对日峰值流量影响显著, 岔口流域径流属于洪峰径流, 因此该参数对径流敏感。此外, 田间持水量对流域径流和泥沙高度敏感, 而在其他研究区敏感性适中或较弱, 这可能与流域土壤质地和理化性质有关, 流域土壤以粉砂质壤土为主, 田间持水量越大, 产流和产沙量越小。

4 模型模拟验证

在参数敏感性分析的基础上, 再次对模型进行微调, 采用经过校准之后的参数对岔口小流域进行模拟验证, 用径流、泥沙和氮磷等观测数据与模型模拟结果相比较来验证模型在流域的模拟精度。观测数据为2013年7月5场暴雨实测值, 采用相对误差D_V、决定系数R²和Nash-Suttcliffe模拟效率系数Ens对模型模拟结果进行验证。其中D_V值越小，拟合度就越高, R²越大则表示实测值与模拟值的相关关系越好, Ens值越大表明模型模拟效果越理想。观测值与模拟值比较结果如表 4~6所示。

径流模拟结果显示, 模型D_V为-4.7%, R²和Ens均为0.99, 模拟结果很理想, 模拟精度符合研究要求, 表明AnnAGNPS模型对流域径流模拟的可行性。

泥沙模拟结果显示, 模型D_V为181.83%, R²和Ens分别为0.43和-0.28, 模拟结果不理想。原因在于:流域土质疏松, 雨季易发生滑坡, 滑坡具有不确定性, 难以监测且模型缺乏滑坡模块; 模型对淤地坝模拟不准确; 泥沙测算采用RUSLE方程, 该方程通常用来测算长时间序列泥沙侵蚀量, 对日泥沙侵蚀量测算精度不高。

氮磷模拟结果显示, 模型D_V分别为152.64%和-1 629.62%, R²和Ens分别为0.04、0.25和-0.28、-0.21, 模拟结果不理想。原因在于:氮磷等营养物的流失与土壤侵蚀有关, 泥沙预测不理想，直接导致氮磷预测精度不高; 流域坡改梯工程实施的影响。

模型预测精度为径流＞泥沙＞氮＞磷, 参考其他研究成果^[26-27]可知模型对磷输出预测误差均偏大且不稳定。模型可从以下2方面进行改进:(1)历时短、降雨强度高的场次降雨是造成黄土丘陵沟壑区水土流失的重要因素, AnnAGNPS模型中采用的RUSLE方程适合长期土壤侵蚀量计算, 对场次降雨泥沙量模拟精度较差。针对黄土丘陵沟壑土壤侵蚀的特点, 应选用适合该地区的土壤侵蚀模型, 研究发现改进版通用MUSLE模型^[28]和Si模型^[29]对黄土丘陵沟壑区次洪或日尺度产沙量模拟精度较高, 因此采用多模型相结合的方法是提高黄土丘陵沟壑土壤侵蚀量预测精度的重要途径。(2)模型impound模块设计具有局限性, 只考虑淤地坝对流域泥沙淤积的影响, 未考虑坝内死水位对径流的扣留和超过死水位时由卧管排出对泥沙和径流流速的影响。结合流域实际情况, 将wetland和impound数据一起使用且共用一个reach id模拟淤地坝对流域产流输沙及污染物的影响, 模拟效果比两者单独使用好, 可为类似地区应用提供参考。虽然模型对泥沙、氮磷等污染物模拟误差较大, 但从总体上还能反映出流域泥沙侵蚀和氮磷流失的趋势, 对流域治理水土流失和非点源污染具有一定的参考价值。

5 结论

采用定量方法评价模型主要参数的敏感性程度, 得到不同参数对模型模拟输出结果的影响。总体而言, CN值对模型模拟结果有显著影响, 且呈明显正相关特征; 对模型径流输出影响最大的参数主要是CN值和FC值; 对模型泥沙、营养盐输出影响最大的参数主要是CN值、K值、C值、LS值和P值。结合流域实际情况应着重考虑CN值、植被覆盖和水土保持情况对模型输出结果的影响。

氮磷污染物负荷空间分布不具一致性, 但区域性明显。土壤氮磷库即施肥量是氮磷污染物的主要来源; 总氮流失主要与径流冲刷有关, 总磷流失主要与泥沙侵蚀有关。

通过模型拟合分析可知:径流预测精度较高, 泥沙和氮磷等污染物预测精度较低, 精度较低的原因与滑坡、坡改梯工程、模型内部的局限性有关。模型总体上能反映流域水土流失及非点源污染趋势, 对流域水土流失及非点源污染的防治具有一定的参考价值。

在对模型参数进行敏感性分析的基础上, 还需进行不确定性研究, 因为敏感度低但不确定性高的参数对模型模拟结果的影响可能大于敏感度高但不确定性低的参数, 这也是敏感性分析的不足之处。