流域下垫面是大气与地面或水面的分界面^[1]。农田化肥、畜禽养殖和农村生活污染排放物迁移、转化最终进入河道的过程在很大程度上受下垫面径流条件影响^[2-4]。大量研究表明，由于各种面源污染源空间分布的差异性、污染排放的随机性以及下垫面径流过程的时空变异性，不同土地利用条件下径流过程以及面源污染负荷表现出明显差异性^[5]和不确定性^[6-8]。分布式水文和面源污染模型（如SWAT ^[9-10]、AGNPS ^[11]、CASC2D ^[12]和DWSM^[13-14]等）通常根据下垫面土壤性质、土地利用、地形等方面的差异，将汇流区划分为一系列基本水文单元，采用不同的参数描述对应下垫面中的水文过程和污染物迁移、转化、入河过程。然而，在流域尺度上难以获取详细的面源污染负荷及通量监测数据，通常仅能够通过若干控制断面的水文和水质监测数据进行参数率定和模型分析，模型中参数的敏感性在很大程度上受监测点位置、监测时间和长度、汇流区大小等多种因素的影响。了解下垫面对面源污染负荷的影响机制，对于分布式模型计算精度的提升以及面源污染监控都具有重要的意义。

笔者以涪江流域为研究对象，基于SWAT模型模拟流域水文和污染负荷，采用SWAT-Cup进行参数率定和模型校验，基于Sobol全局性方法^[15]分析流域不同下垫面参数对面源污染负荷的影响，以期为流域面源污染防控决策提供依据。

1 面源污染负荷监测及分析方法 1.1 基于SWAT模型的涪江流域面源污染通量过程模拟

涪江全长670 km，流域（29. 30° ~33. 05° N，103. 73°~106. 27° E）面积36 400 km²，位于三峡水库上游，四川省绵阳市、遂宁市和重庆市潼南区、铜梁区境内。海拔介于200~5 500 m之间，流域多年平均气温17. 8~18. 1℃，多年平均年降水量800~1 400 mm。流域内水系、水文站、气象站和水质自动监测站位置如图 1（a）所示，流域内共有平武、冶城、江油、涪江桥、梓潼、罗江、三台、潼川、射洪、南北堰和小河坝11个水文监测站，逐日测定河道流量，时间序列为2007年1月—2016年12月。流域内及周边有九寨沟、平武、安县、北川、江油、中江、青川、剑阁、苍溪、梓潼、三台、盐亭、西充、射洪、蓬溪、安岳、潼南和武胜18个气象站，测定每日的日均温度、最高温度、最低温度、日照时数、日均风速、日均湿度和日降雨量7个参数，时间序列为2006年1月— 2016年12月。流域内香山、老池和大安3个断面为四川省控水质断面，设置有自动水质监测站，测定氨氮、TP和高锰酸盐指数（I_Mn）3种污染物浓度，每日进行6次测定，3个站的水质监测时间分别为2012年1月—2016年12月、2014年1月—2016年12月和2006年1月—2016年12月。

流域30 m精度数字高程模型（DEM）如图 1（b）所示，土地利用分类如图 1（c）所示。流域内包括农田、农村生活、畜禽养殖和水土流失4类面源污染源，基于农村人口数量、农田类型及面积、畜禽种类及养殖数量、土壤类型、土地利用方式和地形数据，采用HJ 884—2018《污染源源强核算技术指南准则》核算面源污染物排放量。点源污染包括工业点源和城镇污水处理厂2类，均采用四川省环境监测总站环境统计数据。流域内点源分布和单位面积面源污染物排放量如图 1（d）所示，涪江流域土壤TN和TP含量分布分别如图 1（e）和图 1（f）所示。

选择流域内涪江上游、通河口、梓潼江、老池、大安及香山6个断面对应的汇流区，研究不同下垫面条件下面源污染负荷的构成及其变化机理。涪江上游断面为涪江江油市与平武县的交界断面，大安断面为涪江的一级支流安居河四川—重庆交界断面，老池断面为涪江干流绵阳市和遂宁市交界断面，香山断面为四川省和重庆市交界断面，梓潼江与通口河断面为区内典型支流出口断面，6个典型汇流区基本情况如表 1所示。

基于ArcGIS软件将流域基本信息进行模型展布，确定SWAT模型的输入信息后，以2005—2006年作为模型预热期，消除初始条件对模拟结果的影响，采用2007—2014年的水文流量以及3种污染物的浓度监测结果进行模型参数调整，采用SWAT-Cup进行参数率定。将2015—2016年作为模型验证期，采用率定后的参数进行模型校验。

1.2 下垫面条件对面源污染负荷的影响分析

下垫面对污染负荷的影响可以表现为3个方面：（1）面源污染物排放量，包括农田化肥TN和TP施用量、畜禽养殖和农村生活污染排放进入田间的污染物4个污染指标；（2）下垫面径流条件，包括土壤性质、地表坡度、河网密度以及土地利用类型4个参数；（3）土壤本底，包括土壤氮磷本底含量、有机质含量以及C/N比4个指标。

Sobol方法基于方差分解，通过将目标函数的总方差分解为单个参数的方差和参数间相互作用的方法，实现参数的全局性敏感性分析。对于所选择的3类12个下垫面指标参数中的任一参数，根据其分布特性（流域内下垫面径流参数和土壤本底参数空间分布基本符合对数正态分布，污染物排放量变化基本符合正态分布），采用Monte Carlo方法^[16-17]，基于参数调查值的均值和方差，随机产生2个不同的低差异随机序列：

$ \boldsymbol{M}_{\mathrm{A}}=\left[\begin{array}{cccc}{x_{11}} & {x_{12}} & {\cdots} & {x_{1 n}} \\ {x_{21}} & {x_{22}} & {\cdots} & {x_{2 n}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {x_{m 1}} & {x_{m 2}} & {\cdots} & {x_{m n}}\end{array}\right]，$

(1)

$ \boldsymbol{M}_{\mathrm{B}}=\left[\begin{array}{cccc}{x_{11}} & {x_{12}^{\prime}} & {\dots} & {x_{1 n}^{\prime}} \\ {x_{21}^{\prime}} & {x_{21}^{\prime}} & {\cdots} & {x_{2 n}^{\prime}} \\ {\vdots} & {\vdots} & {} & {\vdots} \\ {x_{m 1}} & {x_{m 2}} & {\cdots} & {x_{m n}}\end{array}\right]。$

(2)

式（1）~（2）中，m为矩阵行数，即抽样重复数；n为下垫面参数个数。M_A和M_B这2个随机序列具有不同的统计特征，反映下垫面状态的变化，利用矩阵M_B中的第i列数据（下垫面参数x_i的抽样值）替换矩阵M_A中的第i列数据（下垫面参数x_i的抽样值），形成一个新的联合矩阵M_Ci（i=1，2，3，…，n）：

$ \boldsymbol{M}_{c_{i}}=\left[\begin{array}{cccccc}{x_{11}} & {x_{12}} & {\cdots} & {x_{1 i}} & {\cdots} & {x_{1 n}} \\ {x_{21}} & {x_{22}} & {\cdots} & {x_{2 i}^{\prime}} & {\cdots} & {x_{2 n}} \\ {\vdots} & {\vdots} & {} & {\vdots} & {} & {\vdots} \\ {x_{m 1}} & {x_{m 2}} & {\cdots} & {x_{m i}} & {\cdots} & {x_{m n}}\end{array}\right]。$

(3)

对每一个下垫面参数都生成一个联合矩阵，最终形成n+2个矩阵，基于SWAT模型对所有参数矩阵下的水文及面源污染通量过程进行计算。下垫面参数x_i的一阶敏感性指数S_i和总敏感性指数S_Ti分别为

$ {S_i} = \frac{{\frac{1}{{m - 1}}\sum\limits_{k = 1}^m {y_B^{(k)}} y_{{c_i}}^{(k)} - \frac{1}{m}\sum\limits_{k = 1}^m {y_A^{(k)}} y_B^{(k)}}}{{V\left( {{y_A}} \right)}}，$

(4)

$ {S_{Ti}} = 1 - \frac{{\frac{1}{{m - 1}}\sum\limits_{k = 1}^m {y_A^{(k)}} {y_{c_i^{(k)}}} - {{\left( {\frac{1}{m}\sum\limits_{k = 1}^m {y_A^{(k)}} } \right)}^2}}}{{V\left( {{y_A}} \right)}}。$

(5)

式（4）~（5）中，y_A^(k)、y_B^(k)和y_Ci^(k)分别为矩阵M_A、M_B、M_Ci第k行下垫面参数对应的模型输出；V(y_A)为矩阵M_A输出的方差。

由式（4）可以看出，参数的一阶敏感性指数反映参数x_i变化对流域径流或污染负荷模拟结果变化方差所产生的贡献，因而参数的一阶敏感性指数越大，其敏感性越强。同样，总敏感性指数反映该参数与其他参数之间交互作用对径流或污染负荷模拟结果的敏感性，敏感性越强，则表明该参数与其他参数之间的交互作用对径流或负荷的影响程度越大。

2 结果与分析 2.1 基于SWAT模型的涪江流域污染负荷模拟

SWAT模型参数率定分为2个步骤：首先根据下垫面性质的实验结果以及相关性质确定参数的初始值；其次采用SWT-Cup确定参数的最终取值，使实测所有监测数据与模拟数据之间的误差达最低。SWAT-Cup分析表明，CN₂、ALPHA_BF、GW_DE- LAY、OV_N和CNCOEF等参数对径流过程具有显著影响（P < 0. 05），CH_K对于径流过程和污染物浓度均有明显影响，USLE_P、CMN、PPERCO和ERORGP等参数对NH₃、I_Mn和TP浓度的影响最为显著（P < 0. 05），涪江流域SWAT模型关键参数取值如表 2所示。

涪江流域全部监测断面的流量和污染物浓度误差标准均值如表 3所示。模型参数率定期和模型校验期涪江上游江油站、中游射洪站、下游小河坝站流量模拟和监测结果的比较，以及大安站、香山站、老池站NH₃、TP和I_Mn实测浓度和模拟浓度的比较见图 2。

采用Nash-Sutcliffe系数^[18]、相对均方根误差RRMSE、相对偏差F_B和相对总误差F_GE^[19]4个指标对流量以及NH₃、TP和I_Mn浓度模拟误差进行评价。SWAT应用于流域水文和污染负荷模拟的研究结果表明，流域尺度在10² km²级别时，Nash-Sut-cliffe系数普遍能够达到0. 58~0. 72；在10³ km²级别时，能够达到0. 51~0. 55；在10⁴ km²级别时，较难超过0. 55^{[5, 20-22]}。

由表 3可知，流量以及3种污染物浓度模拟结果的Nash-Sutcliffe系数普遍超过0. 5，具有较好的模拟效果。F_B用于描述计算值和模拟值之间的系统误差，若I_Mn的模拟值F_B/F_GE在±0. 3/0. 5范围内，则说明模拟值精度能够被接受^[22]。流量、NH₃和TP的模拟值F_B/F_GE在±0. 15/0. 3范围内，表明模拟精度良好。

2.2 不同下垫面的面源污染负荷特性分析

基于验证后的SWAT模型，将流域内的点源排放量全部置零，并将18个气象站的均值作为模型输入，以消除点源对河道污染通量的贡献以及流域降雨时空分布不均匀对面源污染通量的影响。图 3为零点源通量和流域无差别降雨情况下的面源NH₃、TP和I_Mn负荷（单位面积通量）比较，丰水期（6—9月）TP负荷在年总负荷中占44. 4%~73. 8%，平水期（4—5月、10月）占22. 9%~34. 2%，枯水期（11—2月）占5. 7%~12. 4%；丰水期I_Mn负荷占年负荷的74. 4%~92. 1%，平水期占18. 6%~29. 4%，枯水期占6. 1%~12. 4%，丰水期、平水期和枯水期NH₃负荷分别占年总负荷的53. 3%~71. 4%、24. 4%~32. 4%和15. 2%~22. 4%。

6个汇流区NH₃、TP和I_Mn负荷的比较如表 4所示。对比6个汇流区平均面源污染风险因子（表 1）和面源污染负荷（表 4）可知，涪江上游和通口河污染物排放量和下垫面径流条件均明显低于其他4个汇流区，大安站NH₃、I_Mn和TP这3种污染物负荷最小，梓潼江断面NH₃负荷明显超过其他断面，而香山和通口河断面的TP和I_Mn负荷值最大，且反映污染负荷变化幅度的标准差表现出随污染负荷增加而增大的趋势。相比其他汇流区，涪江上游不同水文期NH₃变化幅度最大，而I_Mn变化幅度最小，相同径流条件下，NH₃和I_Mn负荷变化表现出明显差异。这与I_Mn主要来源于农村生活和畜禽养殖，而NH₃则来源于多种污染源有关。

在下垫面以山丘区、自然林草为主的汇流区（如大安和涪江上游），下垫面径流条件中的坡度和坡长因子对面源污染负荷的影响最大。I_Mn负荷主要来源于土壤中有机质以及分散分布在流域各处的畜禽养殖和农村生活污染，下垫面径流条件中坡度和坡长对其的影响最为大。比较表 1和表 4可以看出，尽管老池汇流区的污染物排放量因子和下垫面径流条件均高于香山汇流区，然而老池断面汇流区内NH₃、TP和I_Mn这3种污染负荷及其变化幅度均在一定程度上小于香山汇流区，老池汇流区内农田（特别是水田区）面积比例明显增加，由于平原区农田对径流过程具有较强的滞留能力，很大程度上减少了径流入河量，因此降低了汇流区内的污染负荷，其他汇流区的情况基本相同。

在6个汇流区中，以农田为主的大安汇流区3种污染物负荷均为最小。此外，在以农田为主的汇流区，TP以悬移质携带入河为主要形式，由于农田对径流具有一定的调节能力，因此不同水文期TP负荷变化最小（图 3）。

涪江上游和通口河下垫面类型和面积占比较为一致，然而由于面源污染排放量因子差别明显，3种污染物负荷存在较大差异，面源污染排放量因子高的污染负荷变幅均显著大于面源污染排放量因子低的变幅。相比污染排放量因子，下垫面径流条件对于通量的影响更为明显。

2.3 下垫面对面源污染负荷影响敏感性分析

根据子流域内的农业区种植作物和化肥使用量的调查结果确定化肥中TN和TP施用量，根据农村畜禽养殖数量以及农村位置确定农村生活、畜禽养殖污染入田量。按照下垫面实际情况，基于GIS信息提取土壤渗透系数、平均坡度、河网密度、土地利用类型以及土壤本底信息。面源污染排放量、径流条件和土壤本底共12个参数的Sobol一阶和总敏感性指数如表 5所示。

NH₃、TP和I_Mn的一阶敏感性指数均值分别为0. 076、0. 054和0. 119，其中径流条件所占比例分别为71. 1%、88. 8%和55. 6%。单因子的敏感程度比较而言，径流条件的影响最大。面源污染物排放量、径流条件和土壤本底的总敏感性指数均值分别为0. 048、0. 161和0. 077，说明多因子协同条件下，径流条件仍为影响面源污染负荷的最敏感因子，而污染排放量因子的影响最小。已有研究也表明，区域尺度下主控流动途径对污染物迁移和入河过程的影响程度高于其他因素^[23]。

NH₃、TP和I_Mn这3种污染物的总敏感性指数平均值分别为0. 113、0. 072和0. 140，可以看出，TP对于下垫面的敏感性最弱，而I_Mn对下垫面的敏感性最强。这一结果表明，涪江流域TP来源于农田以及其他下垫面表层土壤，以悬移质携带的形式进入河道中，由于其来源和入河途径较为单一，因此在3种污染物中敏感性最弱，且径流条件的单因素敏感性和总敏感性在全部因子敏感性中占比均为最大。由于I_Mn来源高度分散，加上不同来源中径流入河条件存在差异，因此其敏感性最强，但其径流条件单因素敏感性和总敏感性在全部因子敏感性中占比均为最小。一些研究成果表明，流域的不同污染物通量过程之间及其与径流过程之间均存在明显差异，造成这种差异的主要原因在于下垫面径流过程与污染负荷之间的非线性响应过程^[24-25]。笔者对各种因子的敏感程度进行进一步分析，NH₃、TP和I_Mn污染物分别有5、3和6个因子的总敏感性指数超过0. 1，达到显著水平^[26]，多因子作用是造成污染负荷与径流过程非线性响应的重要原因。

3 结论

（1）基于SWAT模型模拟了涪江流域径流过程和面源污染负荷，采用SWAT-Cup进行了参数率定。采用Nash-Sutcliffe系数、相对均方根误差、相对偏差和相对总误差4种指标进行检验，结果表明模拟径流过程和污染负荷与实测值符合良好。

（2）各种下垫面条件下，NH₃、TP和I_Mn负荷均随着负荷均值的增加而增大，在以山丘区、自然林草为主的汇流区下垫面条件下，径流条件因子中的坡度和坡长对于面源污染负荷的影响最为明显，由于农田对于径流过程具有一定的调节能力，面源污染负荷随汇流区内农田面积比例增加而降低。

（3）采用Sobol方法分析下垫面中污染物排放量、径流条件和土壤本底对面源污染负荷的影响，结果表明3类因子中，径流条件对面源污染负荷的影响最为敏感，而排放量的影响最弱。对于3种污染物，下垫面对TP负荷的敏感性最弱，对I_Mn负荷的敏感性最强。

（4）涪江流域面源污染分析表明，径流条件对于河道污染负荷的影响最大，控制下垫面径流条件对于降低面源污染负荷具有积极作用。