松涛水库上游流域土地利用类型变化对氮负荷影响

雷湘龄 叶长青 朱丽蓉 薛杨

（1.海南省林业科学研究院 海南海口 571100；2.海南大学生态与环境学院 海南海口 570228；3.海南文昌森林生态系统国家定位观测研究站 海南文昌 571300；4.海南大学旅游学院 海南海口 570228）

随着各国对水污染问题的重视，点源污染已逐步得到有效的控制和治理，非点源污染成为水污染的主要来源[1]，北京密云水库、安徽巢湖、天津于桥水库、上海淀山湖、云南洱海等水域，非点源污染比率超过了点源污染[2]。土地利用类型的变化使土地下垫面发生改变，不适当的土地利用方式和农业管理措施在降雨-径流驱动因子的作用下加剧了土壤侵蚀，化肥中过量的氮、磷污染物附着于泥沙颗粒及可溶性氮磷污染物进入水体，从而形成非点源污染[3-5]。

水文模型模拟是定量描述流域非点源污染过程的重要手段[6-7]。SWAT（Soil and Water Assessment Tool）模型被称为是以农业和森林为主的流域中较有前途的具有连续模拟能力的非点源污染模型。Anna等[8]使用SWAT模型对欧洲南部的多瑙河流域水质进行评估，发现流域内主河道的总氮和总磷主要来源于农作物种植施用的化肥；陈健等[9]和耿润哲等[10]的研究结果也表明，农业活动所使用的化肥、农药等带来的营养物流失是流域内发生非点源污染的主要原因。大量的研究表明，土地利用变化会对非点源污染产生影响，Sushil等[11]在印度莫罗尔流域的研究说明，采用等高耕种和种植植被过滤带能够减少泥沙和土壤中营养物质的流失；刘晶晶等[12]模拟湖北省清江流域宜都段2020年的土地利用情景，结果显示，在耕地减少、林地增加、城镇面积扩大的情景下，总氮和总磷负荷分别减少23.77%和29.53%；秦耀民等[13]对陕西黑河流域不同土地利用情景进行模拟，发现流域内非点源污染负荷随着林地面积的增加、耕地和未利用地的减少而减少；Lin等[14]在伊利湖西部流域的研究结果显示，采用合理的农业措施能够降低流域内地表水体中氮、磷等营养物的负荷。

近年来非点源污染成为当下研究的热点之一，国内已对我国西北部、南部、东部典型流域的非点源污染开展研究，但鲜见关于热带地区河流的研究报道。海南省位于我国的最南端，属于热带季风气候，松涛水库是海南省最大的水库。本研究以松涛水库上游流域 1975—2014年的气象数据和 1988、2002、2016年土地利用图构建SWAT模型，分析土地利用变化（砍伐热带雨林种植橡胶等经济林）对非点源污染氮负荷的影响，并识别造成非点源污染的“源-汇”区。

1 材料与方法

1.1 研究区概况

松涛水库上游流域位于海南省白沙县内（图1），流域内山地多，森林资源丰富，平均气温高，日照长、雨量丰沛，以台风带来的降雨为主；多年平均降雨量为 1 859 mm，年最大降雨量可达2 637.8 mm，存在明显的干湿两季，11月至次年4月为干季，降雨量占全年的15.9%，5—10月为湿季，占全年降雨量的84.1%。

图1 研究区域图

以流域内福才水文站点的水质数据进行模拟研究，控制的集水面积为 492.7 km2。根据 2012年海南省统计年鉴，白沙县一年施用氮肥8 688 t，耕地和水旱田面积共有268.38 km2，平均氮肥施用320 kg/hm2（折纯氮）。流域内有2个排污口，分别用于南开乡卫生院排污和生活用水排污，氨氮负荷分别为0.03和0.04 t一年（表1）；而根据福才站的实测数据，氨氮负荷每年约有5 t；由此可见，在该流域内水质污染以非点源为主。

表1 研究区内点源排放统计

1.2 方法

1.2.1 模型简介SWAT模型中，水文循环陆地阶段基于水量平衡的方程为：

式中：SWt表示土壤最终含水量（mm）；SWo表示第i天的土壤初始含水量（mm）；t表示时间（d）；Rday表示第i天的降水量（mm）；Qsurf表示地表径流总量，单位为mm/hm2；Ea表示第i天的蒸散发量（mm）；Wseep表示第i天离开土壤剖面底部的渗透水流和旁通水流的水量（mm）；Qgw表示第i天回归流的水量（mm）。

修正的通用土壤流失方程（USLE）如下：

式中：sed表示某天的产沙量（t）；Qsurf表示地表径流总量（mm/hm2）；qpeak表示洪峰流量（m3/s）；areahru表示水文响应单元HRU的面积（hm2）；KUSLE表示方程中的土壤可侵蚀因子0.013 t·m2·h/(m2·t·cm)；CUSLE表示土壤覆盖与农业管理措施因子；PUSLE表示水土措施保持因子；LSUSLE表示地形因子；CFRG表示粗糙度因子。

1.2.2 数据来源构建SWAT模型需要建立流域的空间数据库和属性数据库以保证模型模拟的精确性，构建SWAT模型所需数据来源见表2。

表2 构建模型所需的数据来源

1.2.3 模型的构建

1.2.3 .1 数字高程模型数据（DEM）采用SWAT模型对 DEM 进行填洼，计算河流流向，对流累积栅格进行预处理，通过设置上游汇水面积的阈值划分子流域，共划分16个子流域。

1.2.3 .2 土地利用图使用 ENVI软件的监督分类方法，结合谷歌地球影像图，采用人机交互的方式对30 m×30 m遥感影像进行解译，将研究区的土地利用类型分为七类：浆纸林、水域、橡胶林、耕地、居民区、果园、天然林。

1.2.3 .3 子流域和水文响应单元划分SWAT模型属于分布式的水文模型，通过离散化来体现流域的空间分布[15]。在子流域上叠加土地利用图、土壤类型图和坡度图，进一步划分水文响应单元。

1.2.4 参数的率定与验证使用 SWAT-CUP的SUFI-2算法对SWAT模拟输出的数据进行率定和验证。SWAT 模型中涉及径流与营养物负荷的参数众多，本研究通过查阅相关文献和 SWAT-CUP的多次迭代确定参数，利用 SUFI-2 的全局敏感性分析方法评价参数的敏感性，敏感性越高表示该参数对径流或营养物负荷的贡献率越大，选用的径流参数及其率定与验证见表3、4。

表3 SWAT模型径流敏感性参数及其率定值

以Nash-Sutcliffe效率系数Ens和确定性系数R2为评价指标，R2越趋近于 1，表明实测值与模拟值的拟合度越高[16]。当Ens>0.75时，模拟结果良好；当0.36≤Ens≤0.75时，结果令人满意，说明模型适用于该流域；当Ens<0.36，模拟效果较差，表明模型并不适用[17]。

表4 径流参数的率定与验证

1.2.4 .1 径流的率定和验证根据 Mann-Kendall法对年径流系数突变检验的结果，选取 1977—1996年为模型率定期（1975—1976年设置为模型预热期），1997—2014年为验证期。径流参数率定与验证的R2和Ens均达到0.90以上，说明SWAT模型在该流域内对径流具有较好的模拟效果。

1.2.4 .2 氨氮的率定和验证福才水质监测数据为2012—2014年，设置2012年为率定期，2013—2014年为验证期。氨氮率定与验证的Ens均在0.54以上，说明SWAT模型适用于该流域营养物负荷的模拟，见表5、6。

表5 SWAT模型氨氮敏感性参数及其率定值

1.2.4 .3 总氮的率定和验证根据福才水文站点的相关监测数据得到 2013年的总氮数据，于2014年分别对下游的南丰水文站点氨氮和总氮浓度数据进行监测，计算南丰站点氨氮占总氮的百分比，约为25%；以此比例估算福才站点2014年的总氮浓度，并对估算出的总氮数据进行率定与验证（表7、图2~4），以2013年为率定期，2014年为验证期。

图2 1988、2002、2016年土地利用下径流的率定与验证

表6 氨氮的率定与验证

表7 总氮的率定与验证

图3 1988、2002、2016年土地利用下氨氮的率定与验证

图4 1988、2002、2016年土地利用下总氮（TN）的率定与验证

2 结果与分析

2.1 土地利用类型的变化

研究区内土地利用类型以天然林和经济林为主，1988—2016年浆纸林、橡胶林和耕地的面积有所增加，橡胶林面积占比由1988年的9.74%增加至 2016年的 11.30%，浆纸林由 5.19%增加至9.49%；而天然林面积逐渐减少，由 80.15%降至64.63%（表8），全年总氨空间分布情况见图5。

图5 1988、2002、2016年土地利用下的全年总氮空间分布图

表8 土地利用类型及其面积比

2.2 总氮负荷的时空分布

降雨是导致地表土壤养分流失的外部动力条件[18]。雨滴落下冲刷地表，使得土壤中的氮、磷等营养元素流失[19]，在降雨-径流驱动因子的作用下进入水体，从而产生了非点源污染。松涛水库上游流域 5—10月为湿季，降雨量占全年的84.1%，其总氮输出约占全年 85%；干季降雨量占全年的 15.9%，总氮输出占比约为 15%。有研究表明，山东省小清河流域6—9月的降雨量占全年的70%，在该时间段内非点源污染负荷的输出占比为全年的一半以上[20]；清水河非点源污染物总氮、总磷与月均降雨量、产水量的变化趋势基本一致，当降雨量达到峰值时，总氮、总磷输出量也到达最大值[21]。上述研究均与本文的研究结论相符，因此在该流域内非点源污染负荷以湿季输出为主。

降雨产生径流，而径流在迁移的过程同时还受到土地利用方式的影响，因此不同的土地利用方式下非点源污染物的输出存在较大的差异[22]。采用SWAT模型分别模拟1988、2002、2016年三期土地利用下总氮负荷的变化，结果表明，1988—2016年流域全年氮负荷量表现为显著的增加趋势；与 1988年相比，2002年流域内全年总氮负荷增加了114.6%，2016年增加了171.65%，在此期间流域内天然林面积减少，耕地、浆纸林、橡胶林的面积增加。经济林面积的增加意味着施肥量的增加，但作物对氮肥当季的利用率仅为30%~35%[23]，剩余的营养元素附着于土壤颗粒，当发生降雨时，随径流迁移入水体，从而造成了流域内非点源污染氮负荷量的增加。

由总氮负荷的空间分布（图5）可知，1988年土地利用下的总氮负荷量最大的是7和16号子流域，2002年是 3、5、7、8、16号子流域的总氮输出量最大；相比于 1988年，2002年属于总氮负荷最大梯度等级的子流域增加了3个，其中3和8号子流域的耕地面积增加，5号子流域的橡胶面积增加；16个子流域中总氮输出增长幅度最大的是1号子流域，涨幅为 223.58%，其耕地面积增加。2016年3、5、6、7、8、16号子流域的总氮输出量最大，相比于1988年，总氮负荷最大梯度等级的子流域增加了4个，增长幅度最大的是8号子流域，涨幅为407.66%，其浆纸林和耕地面积大幅度增加。总氮输出量高的子流域内均以耕地和经济林为主，1988—2016年间天然林面积减少，由人工开垦种植的耕地、经济林面积增加；土地下垫面的改变，施肥量的增加，造成氮负荷量的增加。

2.3 土地利用与总氮负荷的关系

土地利用类型的单位面积氮负荷量从高到低依次为：耕地>果园>橡胶林>浆纸林>天然林。不同土地利用类型的氮负荷量与植被的覆盖度相关，土壤中氮元素的流失是地表径流与土壤中氮素相互作用的结果，当植被覆盖度增加，土壤颗粒与地表径流充分作用，使水流速度减慢，同时使得随着径流迁移的粗颗粒沉淀；反之，泥沙中细颗粒含量增加，所以植被覆盖度的增加有利于泥沙中细颗粒的富集，从而导致泥沙中全氮含量增加[24]。在研究区内，耕地的单位面积氮负荷量最大，达到7.70 kg/hm2，在同等施肥条件下，相比于果园、橡胶林和浆纸林，耕地的植被覆盖率最低；其次是果园，而天然林在当地生长多年，植被覆盖率高，无人工干预不施用化肥，因此单位面积氮负荷最低，为1.76 kg/hm2，其土壤中的氮营养物质主要来源于生物固氮和微生物[25]。

土地利用类型对全年总氮输出的贡献率由大到小依次为：天然林>橡胶林>浆纸林>耕地>果园。天然林的单位面积氮负荷量最小，但其面积占比最大，因此输出氮负荷量的贡献率最大。1988—2016年，不同土地利用类型的氮负荷贡献率与其面积变化相关，耕地、橡胶林、浆纸林的面积呈增加的趋势，氮负荷贡献率随之增加，如1988年的耕地面积为25.56 km2，2016年为46.75 km2，其贡献率从3.37%涨到了15.25%（表9）；1988—2016年浆纸林面积增加了21.19 km2，氮负荷贡献率增加了一倍；相反，天然林面积减少了76.43 km2，贡献率从1988年的58.4%下降至2016年的34.33%。

对1988、2002、2016年研究区内各子流域不同土地利用类型的面积与其输出的总氮量进行Pearson相关性分析，见表8。1988年天然林与总氮负荷呈极显著相关；2002年果园与总氮负荷显著相关，浆纸林、天然林与总氮负荷呈极显著相关；2016年耕地、橡胶林与总氮负荷显著相关，浆纸林、天然林与总氮负荷呈极显著相关。天然林在流域内面积占比最大，因此在三期土地利用下均与总氮负荷呈极显著相关，但随着面积的减小，其相关系数逐年降低；相反，耕地、橡胶林、浆纸林的面积随年限有所增加，氮负荷量随之增加，其相关系数逐渐增大，甚至在2016年达到显著或极显著相关。

将景观生态学的“源-汇”理论应用于非点源污染中，“源”是土壤养分流失的源头，“汇”则是能够接纳迁移物质，并对非点源污染具有滞缓作用的场所[26-27]。结合不同土地利用类型氮负荷的贡献率（表9）和Pearson相关分析（表10），研究区内耕地在1988和2002年与总氮负荷呈负相关，但随着面积的增加，2016年达到显著相关，其贡献率也逐年增加；橡胶林、浆纸林也随着面积的增加，相关系数和贡献率增加，因此耕地、橡胶林和浆纸林对非点源污染氮负荷起到“源”的作用。天然林由于在流域内面积占比最大，其贡献率最高且与总氮负荷呈极显著相关，而天然林具有涵养水源、保持水土等生态功能，因此起到“汇”的作用，滞缓非点源污染的发生。

表9 不同土地利用类型的总氮负荷

表10 土地利用类型的面积与总氮负荷的Pearson相关系数

3 讨论与结论

3.1 讨论

降雨是非点源污染产生的直接驱动力[28]，研究区内存在明显的干湿两季，湿季降雨量占全年的 84.1%，其总氮输出量占全年的 85%，总氮负荷与月降雨量的变化趋势基本一致，随着降雨量的增加，流域内总氮输出升高。

在降雨条件下，非点源污染的负荷强度还受到土地利用方式的影响。不同土地利用类型的单位面积氮负荷量从高到低依次为：耕地>果园>橡胶林>浆纸林>天然林，其与植被覆盖度和化肥施用量相关，因此耕地的单位面积氮负荷量最大，相反，天然林最小。同理，流域内的氮负荷变化和空间分布与土地利用的分布相关，1988—2016年天然林面积减少，耕地、橡胶林、浆纸林和果园种植面积增加，在此期间总氮负荷呈增加的趋势；相比于 1988年，在相同的划分条件下 2002年总氮负荷最大梯度等级的子流域增加了3个，2016年增加了4个，总氮增加的子流域均是由于耕地或其他经济林的种植面积增加了。不同土地利用总氮输出的贡献率从高到低依次为：天然林>橡胶林>浆纸林>耕地>果园，贡献率的大小与面积有关；Pearson相关性分析与贡献率分析的结果基本一致，天然林在流域内所占面积最大，氮负荷贡献率最高，与总氮负荷的关系均表现为极显著相关，但随着面积的减小，其相关系数和贡献率逐渐降低；耕地、橡胶林、浆纸林随着面积的增加，贡献率升高，与总氮负荷的相关系数逐渐增大，到了2016年达到显著或极显著相关。

当前非点源污染已成为水体污染的一大主要污染来源，而由农业造成的非点源贡献率已超过50%。松涛水库是海南省重点饮用水源保护区，其上游流域山多地广，森林资源丰富，以农业为主导产业，因此如何有效地防控非点源的发生已成为不容小觑的问题。本研究成果可为当地防治农业非点源污染提供参考。

3.2 结论

研究使用SWAT模型对松涛水库上游流域进行了非点源污染负荷的模拟，分析了土地利用变化对氮负荷的影响，结果如下：

（1）对以1988、2002和2016年土地利用数据为基础构建的SWAT模型进行率定与验证，对月径流率定和验证的纳什效率系数Ens均为0.90以上，氨氮和总氮的Ens分别为0.54和0.46以上，说明SWAT模型适用于该流域的非点源污染模拟。

（2）降雨是造成非点源污染的外部动力条件，湿季总氮负荷全年输出约占85%，干季总氮负荷全年约占15%。1988—2016年流域内氮负荷量呈增加的趋势，氮负荷的空间分布与土地利用类型相关，总氮输出量高的子流域内均以耕地和经济林为主。

（3）单位面积氮负荷量从高到低依次为：耕地>果园>橡胶林>浆纸林>天然林；土地利用类型对全年总氮输出的贡献率由大到小依次为：天然林>橡胶林>浆纸林>耕地>果园；天然林与氮负荷的关系为极显著相关，随着面积的减小，其Pearson相关系数逐渐减小；随着面积的增加，耕地、橡胶林、浆纸林与总氮负荷的Pearson相关系数逐渐增大，在2016年达到显著或极显著相关。