随着工业点源污染控制水平的提高，非点源污染已成为水环境污染的主要因素，引起水体富营养化、 水体生态功能降低和饮用水源污染等问题，尤其是土地利用/覆盖变化导致的为甚^[1,2,3]. 不同土地利用类型的土壤侵蚀和营养物迁移特征存在显著差异^[4,5]，因此其对流域非点源污染的贡献也各有不同. 在定量化制定流域管理措施时，为了寻求环境效益、 经济效益最大化，可以通过优化各种土地利用类型大小及其空间分布，从而实现同步控制非点源污染和增加经济产值^[6,7,8].

在定量估算流域非点源污染负荷方面，许多研究采用依靠水文驱动的机制模型(deterministic modeling)或者基于经验的输出系数法(export coefficient modeling)^[9,10]. 前者应用非点源污染发生的物理机制定量污染物负荷的产出及其对水体的影响，考虑水文过程、 土壤侵蚀过程和污染物迁移过程等^[11,12]； 后者不考虑非点源污染发生的具体过程，而是根据不同单位源(如不同土地利用类型)负荷排放量与输出系数的乘积之和来计算^[13,14,15]. 但前者要求较多依赖建模区域的输入数据，人力财力成本高，在建模区域外的推广和应用受到限制^[16,17]； 后者则没有考虑单位源和受纳水体间的水力联系程度，不能提供连续或者逐月逐年负荷估算值，无法反映随降水径流变化而变化的负荷输出^[16,18,19]，而且由于不合理“借用”输出系数、 缺乏有效的率定和验证方式，导致估算结果较大的不确定性和风险^[9,17,20]，因此两种模型在应用中均存在一定限制.

区域营养盐管理模型(regional nutrient management，ReNuMa)^[21,22]是一种水文驱动的半经验模型，相比机制模型，其所需参数少、 成本低，易于向建模区域外推广； 相对于输出系数法，其含有连续降水径流传输模块，能够反映降水径流影响下的负荷输出，提供逐月逐年的负荷估算值，此外还能够分别给出溶解态、 颗粒态营养盐负荷及其来源，能够更为全面真实地反映流域负荷情况. 目前，ReNuMa已被用于美国东北部16个流域^[21]，中国长春的石头口门水库流域^[23,24]和天津于桥水库的沙河流域^[25]等的非点源污染定量估算. 本研究采用ReNuMa模型对新安江上游率水流域2000~2010年溶解态氮(dissolved nitrogen，DN)负荷进行定量估算和来源解析，不仅提供了流域内非点源DN负荷的来源和占比，确定了不同土地利用类型的DN负荷强度，还按照污染物负荷适量削减和农业经济产值最大化的原则，进行流域土地利用结构优化，以期对解决流域非点源污染问题和规划土地资源提供参考.

率水河是新安江、 富春江和钱塘江的正水源头，与横江交汇于安徽省黄山市中心城区屯溪区，全长148.2 km，流域面积约为1520 km²，流域内林地和耕地为主要土地利用类型. 该流域属于亚热带季风气候，年平均降水量为1752 mm，大多集中在4~8月之间. 流域范围中没有被记录到工业点源污染，污染物主要通过非点源排放的形式进入水体. 率水流域的地理位置关系如图 1所示.

图 1

Fig. 1

图 1 率水流域地理位置 Fig. 1 Location of Shuaishui watershed

ReNuMa模型的输入输出示意图如图 2所示. 其中，输入数据来源及处理步骤列于表 1.

图 2

Fig. 2

图 2 ReNuMa模型输入输出示意 Fig. 2 Input and output flow chart of ReNuMa

表 1(Table 1)

表 1 ReNuMa模型输入数据来源及处理 Table 1 Source and processing methods of ReNuMa model input data 数据名称 来源 处理 降水量 黄山市水文局提供的流域内13个降水站点2000~2010年逐日降水量 利用ESRIArcGis 9.3进行泰森多边形插值获得流域平均降水量 气温 中国气象科学数据共享服务网(http://cdc.cma.gov.cn)提供的景德镇、 衢州和杭州3个地面国际交换站2000~2010年日均气温 3个站点日均气温平均值为流域日均气温 土地利用 黄山市环保局提供的2010年《黄山市统计年鉴》 将统计年鉴提供的流域内各区县不同土地利用面积加和，得到流域内9种土地利用类型，包括水田、 谷物、 经济作物、 茶园、 果园、 林地和草地，根据《黄山市统计年鉴》2000~2010年期间流域内各土地利用类型面积改变不大，本研究不考虑此期间土地利用面积的改变导致的DN负荷改变 人口 地球系统科学数据平台(http://www.geodata.cn)2003年1km格网全国人口空间化数据 利用ESRIArcGis 9.3以河道周边500m为界划分缓冲区，分别计算缓冲区内外的人口数[25]，用于计算生活源DN负荷. 在ReNuMa中，假定缓冲区内外居民生活排放的N负荷分别由短循环系统(short circuit system)和普通系统(normal septic system)处理，具体算法参见文献[21,22] 施肥量 黄山市环保局提供的污染源普查数据和2010年《黄山市统计年鉴》 —

表 1 ReNuMa模型输入数据来源及处理 Table 1 Source and processing methods of ReNuMa model input data

模型的输入参数包括①定义产汇产流产沙过程的传输参数(transport parameter)和②定义营养盐(本研究中为DN)循环过程的营养盐参数(nutrient parameter)，参数初值根据文献资料或者实地监测获得，经过率定后用于正式模拟； 模型的输出结果包括：不同来源的径流深(streamflow,cm)和DN负荷量(DN,t). 采用Microsoft Excel中的规划求解加载项(Excel Solver Add-in)使模型参数获得率定值，此时径流和DN负荷的逐月观测值和模型模拟值差值的平方和最小. 径流观测值来自率水流域唯一的水文站月潭站点(黄山水利局提供)，根据流域出口的黄山林校断面(黄山市环境监测站提供)和负荷-流量相关关系^[26,27]得到的DN负荷用于模型率定和验证.

采用扰动分析法，即分别考察单一输入参数在相对其基线值高(Δx/x=+10%)和低(Δx/x=-10%)两种取值情形下的输出结果变化(Δy/y)与输入参数变化(Δx/x)的比值S，对ReNuMa模型输入参数进行敏感性分析^[21]. 根据敏感值的大小可以说明参数对模型结果影响的大小，敏感值越大则影响越大. 在率定过程中，优先率定敏感值高的参数，能够有效规避“同效异参”效应，提高率定效率和结果可靠性.

采用Nash-Suttcliffe效率系数Ens^[28]和决定系数R2评价模型模拟结果. 一般认为，Ens>0.5，R2>0.6时模型的模拟结果是可接受的^[24,29]，特别当Ens>0.75，R2>0.7时，表明模型模拟结果非常好^[30].

作为安徽省重要生态功能区，率水流域内的区 县工业产值比重很小，主要以农林经济为支柱产业. 根据《GB/T 21010-2007土地利用现状分类》和流域特点，本研究考虑农林用地[水田、非水稻谷物(谷物)、经济作物、 茶园、 果园、 林地和草地]的土地利用结构规划，建筑用地和其他类型土地不予考虑.

以农林经济产值最大化为目标，以耕地面积、 粮食产量、 农业用水量、 水土流失量、 DN负荷和农业可利用土地资源总量为约束条件，以不同土地利用类型的面积作为决策变量，建立土地利用结构优化模型，规划目标的求解过程在LINDO 6.1(Lindo System Inc，USA)软件中进行. 各约束条件中的系数列于表 2.

表 2(Table 2)

表 2 土地利用结构优化模型约束条件中的系数 Table 2 Coefficients in constraints of land use structure optimization 土地利用类型 产值系数/万元 ·(hm2 ·a)-1 粮食产量系数/万 t ·(hm2 ·a)-1 灌溉用水定额/m3 ·hm-2 土壤侵蚀模数/t ·(hm2 ·a)-1 水田 1.42 0.62 2100 13.72 谷物 0.94 0.39 450 39.87 经济作物 2.60 — 750 45.73 果园 13.31 — 900 10.88 茶园 2.00 — — 15.76 林地 0.36 — — 19.96 草地 — — — 26.94

表 2 土地利用结构优化模型约束条件中的系数 Table 2 Coefficients in constraints of land use structure optimization

目标函数： 式中，a_i为各土地利用类型的产值系数[万元 ·(hm² ·a)^-1],由《2011年安徽省统计年鉴》和《安徽省林业发展第十二个五年规划》得出，x_i为各土地利用类型的面积(hm²)，x₁~x₇分别指水田、谷物、经济作物、果园、茶园、林地和草地.

约束条件如下.

(1)耕地面积是水田、 谷物和经济作物面积之和，根据《黄山市土地利用总体规划(2006-2020年)》得出耕地面积(hm²)约束：

(2)根据《黄山市“十二五”农业和农村经济发展规划》，得出粮食产量(万t)约束：

b₁x₁+b₂x₂≥5090

(3)根据《黄山市“十二五”农业和农村经济发展规划》、《黄山市土地利用总体规划(2006-2020年)》、《黄山市农业产业化“635”转型增效行动规划纲要》以及当前实际面积，分别得出水田、 林地、 草地、 以蔬菜为代表的经济作物、 果园和茶园的面积约束空间，并且x_i≥0.

(4)根据(3)中的面积约束上限与相应的灌溉用水定额，得出农业用水量(亿m³)约束：

c₁x₁+c₃x₃+c₄x₄≤0.2045

(5)根据当前土地利用结构下的水土流失量，得出水土流失量(kt)约束：

(6)根据(3)中的面积约束上限与相应的DN负荷强度，得出DN负荷(t)约束：

(7)根据《黄山市土地利用总体规划(2006-2020年)》与当前实际面积，得出农林可利用土地资源总量(hm²)约束：

在率水流域，敏感性参数包括退水系数、 渗滤系数、 不同土地利用类型的径流CN值(表 3)，说明土壤的含水状况与不同土地利用类型的特征对径流深和DN负荷影响较大.

表 3(Table 3)

表 3 ReNuMa模型敏感性参数取值 Table 3 Calibrated parameters with higher sensitivity for model simulation 参数名称 参数敏感值S 参数取值 高(+10%) 低(-10%) 退水系数 0.1027 0.1226 0.0618 渗滤系数 -0.1092 -0.1116 0.0077 不同土地利用类型CN值 水田 5.4744 1.7605 90 谷物 0.5233 0.3439 75 经济作物 1.9451 1.1434 80 果园 0.0363 0.0238 75 茶园 0.2275 0.1860 50 林地 0.0216 0.0175 52 草地 0.0018 0.0014 60 水体 0.0000 0.0245 100 建筑用地 0.1197 0.0467 95

表 3 ReNuMa模型敏感性参数取值 Table 3 Calibrated parameters with higher sensitivity for model simulation

对于传输参数，采用2000~2006年的逐月径流监测数据进行率定，2007~2010年的逐月径流监测数据进行验证. 对于营养盐参数，以DN负荷监测数据相对完备的2007~2010年为率定期，2000~2006年为验证期. 从图 3中可以看出，无论是率定期还是验证期，径流深和DN负荷的Ens和R2都大于0.9，表明模拟的效果非常好，模型运算的结果能够反映流域真实情况.

图 3

Fig. 3

图 3 径流深和DN负荷的实测值和模拟值 Fig. 3 Observed and simulated streamflow and DN load (a)率定期和验证期的径流深； (b)率定期和验证期的DN负荷

在ReNuMa模型中，陆域产生的DN负荷被设定为通过4种途径进入受纳水体：①地表水来源(runoff flux)，指含氮污染物通过降水冲刷地表土壤，随地表径流进入河流； ②地下水来源(groundwater flux)，指含氮污染物下渗进入土壤后经由土壤浅层饱和区进入河流； ③生活源(septic system flux)，指居民的排泄物经过化粪系统处理后向河流排放； ④点源(point source flux)，指由于工业活动形成的含氮污染物向河流排放(图 2).

模型计算结果表明，率水流域2000~2010年均DN负荷为1.1×10³ t ·a^-1，年际变化较大，最大年份是最小年份的2.57倍(图 4). 根据2000~2010年均DN负荷计算结果，径流来源是流域中最主要的DN负荷来源，其中地下水来源的贡献(60.6%)明显高于地表水(20.1%)； 生活源贡献仅占19.3%； 由于研究流域中没有工业活动，因而点源贡献为0[图 5(a)]. 这说明，从陆域向水体迁移过程中与土地利用密切相关的这类DN负荷是率水流域最主要的非点源污染来源. 因此在定量化制定流域管理措施时，有必要细分各种土地利用类型对流域DN负荷的具体贡献，识别对DN负荷贡献最大的土地利用类型.

图 4

Fig. 4

图 4 2000~2010年ReNuMa估算的不同来源DN负荷量 Fig. 4 Estimated DN load of all sources by ReNuMa model during 2000-2010

图 5

Fig. 5

图 5 率水流域年均DN负荷不同传输形式和不同土地利用贡献率 Fig. 5 Contributions of annual average DN load sources in transport way and land uses (a)年均DN负荷不同传输形式贡献率; (b)年均DN负荷不同土地利用贡献率

ReNuMa模型结果中，仅直接给出了不同土地利用类型对地表水来源的DN负荷的贡献值，但由于模型认为地下水DN浓度是由不同土地利用类型的地表水DN浓度根据下渗(infiltration)权重加权平均得到的，因此本研究根据各土地利用类型的下渗权重计算出相应地下水来源的DN负荷，并将地表水与地下水来源的DN负荷之和定义为各土地利用类型对受纳水体贡献的DN负荷. 由各土地利用类型2000~2010的年均DN负荷可见，林地的DN负荷贡献最大(35.2%)，其后为经济作物(21.6%)、 水田(19.4%)[图 5(b)].

负荷强度是指单位面积上的污染负荷量，可以消除面积大小对负荷量分析的影响，是污染分析的重要指标之一^[15,34]. 率水流域DN负荷强度为(0.75±0.22)t ·km^-2，与Wang等^[35]在本流域所处地区的研究结果(0.16~3.37 t ·km^-2)较为一致. 从不同土地利用类型的DN负荷强度来看(表 4)，水田的负荷强度最大[28.60 kg ·(hm² ·a)^-1]，其后依次为经济作物[23.50 kg ·(hm² ·a)^-1]、 果园[22.99 kg ·(hm² ·a)^-1]和谷物[21.80 kg ·(hm² ·a)^-1]，而草地[2.80 kg ·(hm² ·a)^-1]和林地[2.71 kg ·(hm² ·a)^-1]的负荷强度最小，分别仅为水田的9.8%和9.5%.

表 4(Table 4)

表 4 各土地利用类型的年均DN负荷和DN负荷强度 Table 4 Annual average DN load and load intensity of all land uses in ReNuMa model 项目 土地面积 DN负荷 DN负荷强度 面积/hm2 次序 DN负荷/t ·a-1 次序 DN负荷强度/kg ·(hm2 ·a)-1 次序 水田 6097.27 5 174.36 3 28.60 1 谷物 3420.89 6 74.60 5 21.80 4 经济作物 8252.06 2 193.93 2 23.50 2 果园 197.68 8 4.55 7 22.99 3 茶园 6186.73 3 110.92 4 17.93 5 林地 116688.22 1 316.10 1 2.71 9 草地 6139.95 4 17.22 6 2.80 8 水体 157.63 9 1.60 9 10.14 6 建筑用地 973.14 7 4.00 8 4.11 7

表 4 各土地利用类型的年均DN负荷和DN负荷强度 Table 4 Annual average DN load and load intensity of all land uses in ReNuMa model

根据建立的土地利用结构优化模型，以2015年为规划年对流域内的农林用地土地利用结构进行优化(表 5). 从土地结构变化来看，规划年的耕地面积比现状减少1.6%，其中，水田面积增加6.1%而谷物面积减少19.1%； 果园面积增加13.6%，茶园面积增加21.3%； 林地和草地变化不明显. 总体而言，农林土地总面积增加0.7%.

表 5(Table 5)

表 5 土地利用结构优化前后状况对比 Table 5 Land use situation before and after structure optimization 项目 当前实际状况(2010年) 优化规划状况(2015年) 面积/hm2 面积比例/% DN负荷比例/% 经济收益比例/% 面积/hm2 面积比例/% DN负荷比例/% 经济收益比例/% 土地利用类型 水田 6097.27 4.1 19.6 9.6 6468.20 4.4 20.3 9.9 谷物 3420.89 2.3 8.4 3.6 2768.50 1.9 6.6 2.8 经济作物 8252.06 5.6 21.8 23.7 8252.06 5.6 21.3 23.0 果园 197.68 0.1 0.5 2.9 224.63 0.2 0.6 3.2 茶园 6186.73 4.2 12.4 13.7 7501.42 5.1 14.7 16.1 林地 116688.22 79.4 35.4 46.5 116688.22 78.8 34.6 45.1 草地 6139.95 4.2 1.9 ― 6136.57 4.1 1.9 — 约束条件限制量 耕地面积/hm2 17770.22 17488.76 粮食产量/万t 5114.45 5090.00 农业用水量/亿m3 0.2071 0.2122 水土流失量/kt 3191.57 3148.96 DN负荷量/t 890.42 910.95 农林土地总面积/hm2 146982.80 148039.60 农业总产值/万元 90341.46 93243.02

表 5 土地利用结构优化前后状况对比 Table 5 Land use situation before and after structure optimization

从污染和经济变化来看，规划年的DN负荷比现状增加2.3%，经济产值增加3.2%； 经济产值与DN负荷之比由2010年的101.46万元 ·t^-1上升到102.38万元 ·t^-1，增加0.9%.

ReNuMa模拟过程中仅需要较少的参数，模型率定和验证后模拟效果非常好(图 3)，能够反映流域真实情况. 由于ReNuMa模型含有连续降水径流传输模块，其计算结果不仅给出了逐年的地下水和地表水径流深的模拟值，还给出了与径流紧密相关的不同来源的DN负荷，反映了降水径流影响下的负荷输出. 因此应用ReNuMa模型计算流域非点源污染负荷既合理又简便.

本研究结果表明地下水来源的DN负荷大于地表水，这可能是因为该流域的地下水径流深(52.40~141.03 cm)大于地表水径流深(4.24~33.04 cm). 在孙晓庆等^[36]和黄秀梅等^[37]的研究中，也报道了进入地下水中的氮量大于进入地表水中的现象. 这可能是由于降水产流过程中地下水出水量大于地表水，而且氮肥经过土壤微生物转化成不被土壤胶体吸附的硝酸盐，很容易淋溶进入地下水，造成地下水DN负荷较大^[38]. Lerner等^[39]以及Güler等^[40]认为地下水硝酸盐污染严重与过多施用化肥、采用传统的灌溉、 耕作方式以及缺乏土地利用的系统规划有关.

本研究中，农业生产用地(水田、 谷物、 经济作物、 果园和茶园)以16.3%的面积占比贡献了62.2%的DN负荷，而林地以78.8%的面积占比仅贡献了35.2%的DN负荷，说明农业生产用地在很大程度上影响流域DN负荷的总量，是决定DN负荷的关键土地利用类型.

在本研究中，农业生产用地的负荷强度大于林地和草地，这一现象在国内各流域中普遍存在. 对汉江中下游流域^[41]、 太湖地区的西苕溪流域^[42]、 广州流溪河水库流域^[43]和福建晋江流域^[10]的研究都报道了类似结果. 这是因为林地和草地不施用化肥而且较少地受人为活动干扰，而农业生产用地施肥水平高而且受到人为活动频繁干扰. 结合土地利用类型的面积(表 4)来看，由于林地面积是所有土地利用类型中最大的，所以在林地的负荷强度是最小的情况下，林地贡献的DN负荷仍是最大的； 而果园的负荷强度为22.99 kg ·(hm² ·a)^-1，仅次于水田和经济作物，但由于面积占比仅为0.1%，所以果园贡献的DN负荷较小，是所有农业生产用地类型中最小的. 这说明非点源污染负荷量不仅与不同土地利用类型的负荷强度有关，还与各自的面积紧密相关. 所以，可以通过合理规划不同土地利用类型，实现非点源污染的有效控制^[44].

在土地利用结构优化的结果中，水田的面积增加而谷物面积减小，这是因为谷物相比于水田而言，产值与DN负荷强度之比和粮食产量系数都较小(表 2和表 4)，所以在粮食、 环境和经济约束条件综合考量下减少谷物种植面积. 果园和茶园的面积都有较大增加，说明作为流域内主要发展的两大特色农业产业的果园和茶园由于其产值系数高、 DN负荷强度小的优势而能够在土地面积有较大增加(表 2和表 4). 农林土地总面积略有增加而各土地利用类型有不同程度的变化，说明土地利用结构优化一方面较好地对未利用土地进行农业开发，另一方面对现有土地利用类型结构进行适当调整，从而支持特色农业产业的发展.

根据土地利用结构优化结果，发现规划年的DN负荷和经济产值都有所增加，经济产值与DN负荷之比也有小幅上升. 这说明在土地利用结构最优情况下，经济产值的增长依然伴随着负荷的增加，但经济产值的增幅大于DN负荷的增幅. 因此削减DN负荷可以通过降低DN负荷强度的方式实现. 特别是DN负荷强度高、 产值系数低的水田应当给予重点关注. 例如采用少耕免耕^[45,46]和秸秆留茬覆盖^[47]等保护性耕作方式^[48]来减少地表径流和土壤流失，以及测土施肥、 平衡施肥等合理施用化肥技术来提高化肥利用率^[49].

(1)农业生产用地对率水流域中DN负荷的贡献最大，表明人类影响下的农业生产活动是流域非点源污染的最主要来源. 不同土地利用类型的DN负荷贡献率排序为林地>经济作物>水田>茶园>谷物>草地>果园>建筑用地>水体； DN负荷强度的排序为水田>经济作物>果园>谷物>茶园>水体>建筑用地>草地>林地. 非点源污染负荷量不仅与不同土地利用类型的负荷强度有关，还与各自的面积紧密相关.

(2)对率水流域2015年的土地利用结构优化结果表明，在土地利用结构最优情况下，经济产值的增长依然伴随着负荷的增加，但经济产值的增幅大于DN负荷的增幅.