阿什河流域氮磷输出负荷及其与下垫面特征的关系研究

高晓曦，左德鹏，马广文，徐宗学，李佩君

(1.北京师范大学水科学研究院 城市水循环与海绵城市技术北京市重点实验室，北京 100875;2.中国环境监测总站 国家环境保护环境监测质量控制重点实验室，北京100012)

水是人类赖以生存的自然资源，水质好坏对人类生产生活至关重要。随着工业和农业的快速发展，对水资源的需求日益增加，而工业废水、生活污水乱排乱放，农药化肥肆意使用，使不同地区水资源及水生态环境均遭受到不同程度的污染和破坏，环境问题已经成为全社会关注的焦点，其中水环境问题尤其突出。

为解决这一问题，世界上许多国家最初都在点源污染防治上投入了大量人力和物力。但把水环境安全仅仅寄托在点源污染防治上，并不能从根本上解决问题[1,2]。越来越多国家逐渐开始关心非点源污染在水环境污染方面所扮演的角色。美国环保局2003年调查结果显示，农业非点源污染是美国河流和湖泊污染的最大污染源，是造成40%河流和湖泊水体水质变差、地下水体污染和湿地退化的主要原因[3]；欧洲国家农业非点源污染也是造成水体，尤其是地下水污染和地表水磷富集的主要原因，其中由农业非点源排放所造成污染占地表水污染总负荷的24%～71%[4]；我国大量研究结果同样表明[5,6]，农业非点源污染在非点源污染中占决定性作用。土壤中的氮磷等营养物质在降水冲刷和泥沙携带作用下进入河道，造成水生态环境恶化[7,8]。农业非点源污染形成过程随机性大、影响因子复杂多样，其分布范围广、潜伏周期长等特点导致研究和控制非点源污染难度较大[5,9]。

阿什河是松花江干流南岸支流，位于黑龙江省南部。近些年来，随着该地区经济快速增长和粮食产量不断增加，流域水体污染愈发严重，已经被《松花江流域水污染防治“十二五”规划》列为“哈尔滨市污染最为严重的河流”[10]。随着阿什河流域点源污染治理的稳步推进，非点源污染所引发的问题逐渐显现。阿什河流域是我国重要的粮食生产基地之一，化肥农药用量高，施用方式不合理，造成河水污染负荷较高。流域面临的农业生产压力和非点源污染风险，使得开展非点源污染研究迫在眉睫。

基于以上背景，本研究通过在阿什河流域构建SWAT模型，在水文相应单元尺度上分析了该地区的总磷、总氮的空间分布特征，重点分析了不同土地利用类型、不同土壤类型、不同坡度等下垫面特征上氮磷输出负荷的变化规律，以期为阿什河流域非点源污染控制和土地管理提供科学依据。

1 研究区概况

阿什河位于黑龙江省哈尔滨市辖区内(126°40′E～127°43′E，45°05′N～45°50′N)，属于松花江上游右岸一级支流。阿什河发源于尚志蜜蜂山，流经尚志、五常、阿城和哈尔滨等区县，于哈尔滨市东郊注入松花江[11]。流域面积3 545 km2，河流全长257 km。流域内地形起伏不大，上游以山地为主，中下游主要是平原和丘陵。阿什河流域属大陆性季风气候，具体表现为春季干燥多风，夏季降水集中，秋季降温迅速，冬季干旱漫长[12]。流域年内降水量分布不均，多年平均降水量565 mm，但6-9月降水占全年降水量的77%左右[13]。阿什河属山溪性河流，主要以降水补给为主，融雪和地下水补给为辅。流域内径流流量与地势高程呈正相关，上游地势高，降水量量大，径流量相应较大；而下游地势较低，降雨较少，使得径流量也相对较少[14]。阿什河流域是我国重要的粮食生产基地，流域内主要农作物有水稻、玉米和大豆，而农业生产大量施用化肥、农药所造成的非点源污染，导致阿什河阿城以下河段水质污染严重[15]。研究区地理位置如图1所示。

图1 阿什河流域位置示意图Fig.1 The location of the Ashi River basin

2 数据和方法

2.1 数据来源

研究所采用的基础数据包括地理空间数据云提供的90 m DEM高程数据、中国科学院资源环境科学数据中心提供的1∶10万土地利用栅格数据、1∶100万土壤矢量数据集以及1∶25万河网图，采用的这些基础地理信息数据已经过严格的质量检查，数据质量较为可靠，且这些基础数据已被广泛应用于水文模型建模中[16-18]，具有一定的代表性。采用SWAT模型中国大气同化驱动数据集(CMADS V1.0)以及中国气象数据共享网气象资料作为模型的气象驱动数据。水文资料及水质监测数据分别来源于水利部阿城水文站和哈尔滨市环境监测中心，时间序列长度为2010-2015年。研究区农业用地的化肥、农药施用量是根据哈尔滨统计年鉴以及根据实际情况折算得到。由于阿什河产沙量较低，该流域没有进行泥沙监测；选取的水质指标总氮总磷只有8个月的例行监测，监测时间为1月、2月及5-10月。

2.2 研究方法

2.2.1 SWAT模型介绍

SWAT模型是由美国农业研究部农业研究中心开发、具有很强物理机制的分布式水文模型。模型根据流域高程将流域划分为若干子流域，每一个子流域根据土地利用、土壤、坡度的不同进一步划分为不同的水文响应单元(HRU)。模型可用于模拟地表水和地下水水质水量、预测土地管理措施对不同土壤类型、土地利用方式和管理条件的大面积复杂流域的水文、泥沙和农业污染物的影响[19]。

SWAT模型的水文模拟主要是基于水量平衡，其水文循环基于的水量平衡方程为[20]：

(1)

式中:SWt为土壤最终含水量；SW0为土壤初始含水量；t为时间；Rd为第i天降水量；Qsurf为第i天地表径流；Ea为第i天蒸发量；Wseep为第i天土壤剖面底层的渗透量和测流量；Qgw为第i天的地下水出流量。

SWAT模型采用修正的MULSE方程来模拟土壤侵蚀过程。修正的MULSE方程用径流因子代替降水动能，可清楚地分离和输送泥沙的能量，并可用于单次暴雨事件。具体公式[20]为：

Y=11.8(Qpr)0.56KUSLECUSLEPUSLELSUSLECFRG

(2)

式中：Y为土壤侵蚀量，t；Q为地表径流，mm；pr为洪峰径流，m3/s；KUSLE为土壤侵蚀因子；CUSLE为植被覆盖和作物管理因子；PUSLE为保持措施因子；LSUSLE为地形因子；CFRG为粗碎快土壤成分计算因子。

土壤中的氮磷等营养物随地表径流和泥沙输移进入河道，而附着在泥沙上的营养物随泥沙进入地表径流的量为[21]：

(3)

式中：N0为随泥沙进入地表径流的营养物量；C0为地表10 mm土层中有机物浓度；Y为模拟步长内的泥沙产量；An为HRU面积；εN营养物富集系数。

而随地表径流进入主河道的营养物量表示为[21]：

(4)

2.2.2 适用性评价

SWAT涉及的目标函数类型一共有9种，本文采用目前国内外常用的两个评价指标—决定系数(R2)和纳什效率系数[22](ENS)作为参数优化的评价标准，其表达式如下所示：

(5)

(6)

其中R2表征模拟值和实测值变化趋势的一致性，R2越接近于1，表示两者的拟合程度越高；ENS表示实测值和模拟值的偏离程度，ENS越接近于1表示偏离程度越小。当ENS≤0.36时，认为模拟效果不好；当0.36≤ENS<0.75时表示模拟效果令人满意；当ENS≥0.75时，认为模拟效果好[23]。

3 结果与分析

3.1 模型构建

阿什河流域上游地形起伏较大，主要为山区，下游为适合农业耕种及居住的平原地带。由于平原区地形起伏不大，直接采用D8算法容易使得栅格水流流向一致而生成不合理的平行伪河道[24]。为了减少这种误差，采用Burn in算法导入流域实际数字水系图加以修正。基于阿什河流域90 m分辨率DEM高程数据，通过设置集水区阈值2 500 hm2，将流域划分为89个子流域。通过对研究区土地利用类型进行重分类和编码，将研究区划分为包括森林(FRSD)、灌丛(RNGB)、草地(PAST)、水田(RICE)等10种土地利用类型，详细的研究区土地利用见图5(a)。基于中国土壤数据库，查询和计算研究区内各类土壤的相关属性参数，研究区具体的土壤类型可见图5(b)。根据划分得到各子流域内土地利用、土壤类型及坡度范围，最终将研究区划分为1 255个水文响应单元(HRUs)。

3.2 模型适用性分析

基于SWAT模型，根据阿什河流域地理信息数据和属性数据进行研究区模型构建，使用SUFI-2算法对模型水量和水质相关参数进行率定。阿什河含沙量较少，该流域并没有进行泥沙方面的监测，考虑到一部分磷的流失主要是与泥沙的输移有关，因此在实际总磷参数的率定过程中，将泥沙参数和总磷参数统一率定。率定的总体原则是先进行水量参数的率定，然后再分别对总氮和总磷进行率定。本研究选取2008-2009年作为模型预热期，以降低初始条件的影响，尤其是初始土壤水的影响，但其模拟结果不参与模型评价，将2010-2013年作为率定期，2014、2015作为验证期。阿城站月径流、总氮、总磷在率定期和验证期的模拟效果如图2所示。模型适用性评价指标结果如表1所示。

从图2可知，在率定期和验证期，月径流模拟效果总体较总氮、总磷要好，月径流纳什效率系数(ENS)在率定期和验证期分别为0.70和0.67，确定系数(R2)则分别为0.62和0.60，拟合效果较好。总磷的模拟值和实测值变化趋势较为吻合，其中在2012年份出现一定的偏差，但月总磷纳什效率(ENS)在率定期和验证期分别为0.49和0.46，确定系数(R2)分别为0.47和0.45，拟合效果符合模拟的精度要求，因此可以根据校准好的模型进行进一步的分析。

图2 阿什河流域径流、总氮、总磷模拟值与实测值拟合图Fig.2 Fitting results of simulated and measured values of runoff, TN and TP in the Ashi River basin

表1 阿什河流域月径流模拟效果评价指标Tab.1 Evaluation index of monthly runoff, TN and TP simulation in the Ashi River basin

3.3 总氮总磷流失负荷空间分布

基于模型输出结果，在水文响应单元尺度上统计每一个水文响应单元在2010-2015年总氮总磷输出负荷大小，最终得到的阿什河流域多年总氮总磷年平均负荷空间分布如图3和图4所示。根据阿什河流域水体重污染现状，从空间分布上针对性的分析总氮总磷的输出规律。

图3 总氮负荷空间分布图Fig.3 The spatial distribution of TN load

由图3和图4可知，阿什河流域总氮总磷负荷空间分布都主要集中在下游的平原区以及上游的河岸两侧，总氮负荷最大可达到180.91 kg/hm2，总磷负荷最大的为4.59 kg/hm2。氮磷输出负荷较小的区域则主要分布在上游的林区，非点源污染空间变化差异较大。氮磷的空间分布特征与流域地形及植被覆盖类型有很大关系，在上游地形起伏较大、地表覆盖主要为林地的山区，总氮总磷负荷均比较小，只有在河流附近的农田区域氮磷负荷相对较高；在流域下游、适宜居住和耕种的平原区，总氮总磷负荷相对上游普遍偏高。而在下游，氮磷输出负荷高的区域集中在地表覆盖类型为水田的区域，尤其对于总磷；不同的是，土壤类型为草甸黑土的区域，总氮的污染负荷比总磷高。

阿什河流域作为我国重要的粮食生产基地，种植面积大、耕种强度高，使得该流域很容易产生农业非点源污染[25]。化肥施用强度高，而农作物实际吸收较少，多余的营养物则在降水的溶解和冲刷下进入河道，造成河流水质下降和富营养化。考虑到研究区不同区域所受的人类活动影响不同以及研究区本身所固有的自然特征变化，因此分析下垫面特征与氮磷输出负荷变化规律显得十分必要。

3.4 氮磷流失负荷与下垫面特征的关系

选取土地利用类型、土壤类型和坡度3个反映下垫面特征的要素来分别分析其与氮磷输出负荷间的关系。从图5可以看出，阿什河流域下垫面特征变化较为复杂。从土地利用类型上来看[图5(a)]，阿什河流域土地利用类型主要以耕地和林地为主。耕地主要分布在中下游，耕地面积占整个流域面积的46.42%，其中水田和旱地分别占全流域的6.58%和38.05%。林地次之，主要分布在流域上游，其面积约占全流域的44.59%，落叶阔叶林和灌丛则分别占流域面积的38.15%和1.02%。而草地、水体、聚落、荒漠所占面积较小，分别占流域面积的1.99%、0.98%、3.28%、0.17%。从空间分布上来看，水田主要以带状分布在中下游河岸两侧，在流域中游和河口附近分布有少量面积的建设用地。从土壤类型上来看[图5(b)]，流域的土壤类型主要有3种，分别是暗棕壤、黑土、草甸土，各占流域面积的37.51%、29.79%、21.87%，其余类型土壤则分布较少，从空间分布上看，暗棕壤主要分布在流域上游，黑土分布在流域下游，而草甸土则沿河道两侧分布，从上游延伸至下游。从土壤坡度上来看[图5(c)]，研究区主要坡度在9.42°以下，9.42°～ 17.98°(缓坡)以及大于17.98°(陡坡)主要分布在流域上游。流域整体上呈现上游陡峭，中下游平缓。

图5 研究区土地利用、土壤类型、坡度分布Fig.5 The distribution of land use, soil type and slope in Ashi River Basin

3.4.1 氮磷流失负荷与土地利用类型的关系

对研究区主要土地利用上的氮磷单位输出负荷及输出占比进行统计，统计结果如图6所示。从不同土地利用氮磷输出负荷方面来看，农业用地的总氮输出负荷最大，为50.383 kg/hm2，其次为林地、草地、灌丛和水体，分别为11.69、4.87、2.47和0.005 kg/hm2；对于磷来说，输出负荷最高的土地利用是农业用地和草地，分别为0.951和0.920 kg/hm2，其次为灌丛、湿地和森林，分别占0.483、0.370、0.265 kg/hm2。由此可见，耕地是造成非点源污染的主要土地类型，农业生产所施用的化肥、农药、杀虫剂等污染物经迁移进入河道，造成河流水质的下降。不同的是，草地总磷的输出负荷接近于农业用地，可能是畜牧业的发展有关，动物牲畜的粪便以及牧草施肥都会对草地的磷负荷产生影响。从不同土地利用类型氮磷输出总量占比来看，不同土地利用总氮总磷输出占比次序相似，输出占比最高的是农业用地、其次为林地、草地和湿地。林地虽然输出负荷比较低，但其所占流域面积很大，枯枝败叶以及动物尸体等经微生物分解的有机物都会成为氮磷营养物的来源。

图6 不同土地利用总氮总磷输出负荷及总量占比Fig.6 The yield of TN and TP and proportion of different land use

3.4.2 氮磷流失负荷与土壤类型的关系

对研究区的不同土壤类型上的氮磷输出负荷及占比进行统计，统计结果如图7所示。从不同土壤类型氮磷输出负荷来看，总氮输出负荷最大的为城区，达131.25 kg/hm2，其后为水稻土、黑土和草甸土，分别为65.22、39.81、39.18 kg/hm2，比较小的土类为沼泽土、暗棕壤、白浆土，依次是20.81、18.06和14.45 kg/hm2。对于总磷，不同土壤类型的输出负荷排序与总氮相似，最大的为城区，其后分别是水稻土、黑土、草甸土、白浆土、沼泽土以及暗棕壤，各输出负荷分别为1.34、0.87、0.83、0.76、0.50、0.46、0.37 kg/hm2。城区污染负荷大的原因可能与降水对城市地表的冲刷有关。城区地面含有许多污染物质，诸如城市垃圾、动物粪便、施工场地堆积物、空气沉降物等，同时城区地表相对于自然地表，不透水率增加而糙率减少，极易在降水的作用下引发污染物的冲刷和迁移而进入地表水体，最终对水体水质产生很大的影响[26,27]。不同土壤的氮磷输出负荷不同，除人工施加化肥等因素影响外，土壤的物理化学性质也不可忽略。研究区水稻土质地为黏土，表层土壤成碎块状结构，紧实，透水不良，有机质分解慢，养分不易释放，再加上其本身养分就比较充足，使得土壤固持能力有限，所以损失较为严重。研究区黑土土壤表层为黏土，虽呈团块状结构，但土质较松，相比水稻土，淋溶作用更明显。而草甸土与二者相比，腐蚀质层较厚，表层虽为黏土，但土质疏松，多根系，淋溶量以及土壤的吸附和吸收作用更强。不同土地利用磷负荷差异不同于氮负荷的原因可能与磷本身的物理性质有关。磷肥进入土壤后，较难溶解，且溶解的磷容易被土壤颗粒吸附而形成比较稳定的物质，不易被释放和移动[28]。从不同土壤类型氮磷输出占比来看，最高的依次均为黑土、草甸土和暗棕壤，三者累计占氮磷总输出分别为89.8%、90.3%。

图7 不同土壤类型总氮总磷输出负荷及总量占比Fig.7 The yield of TN and TP and proportion of different soil type

3.4.3 氮磷失负荷与不同坡度的关系

采用自然断点法，将研究区坡度分为4个等级，并对研究区不同坡度等级的氮磷输出负荷进行统计，统计结果如图8所示。从不同坡度等级氮磷输出负荷大小来看，氮磷输出负荷随坡度的增加大致呈递减趋势，对于总氮，4个坡度等级总氮输出负荷分别为42.86、22.86、17.95、15.40 kg/hm2，对于总磷，分别为0.91、0.47、0.27、0.28 kg/hm2。原因是研究区耕地主要分布于坡度较低的平原区，耕地因农业耕作化肥施用量高，故单位负荷输出高。当坡度增加到一定程度，不宜种植业的发展，土壤中氮磷将由地表覆盖和自身土壤条件决定。坡度的增加使得地表所受的雨水条件将发生变化，坡度越高，流速越快，侵蚀作用也相应增加。因为磷易被土壤固定，而氮素容易被淋洗而被径流携带损失，所以氮素主要以径流流失为主，而磷容易吸附到泥沙上，随泥沙迁移而流失。张梦等[29]研究表明坡度较小时，磷素流失途径以径流流失为主，随着坡度的增加，磷素的流失途径以泥沙流失为主。这可能是坡度在18.74～33.07和大于33.07磷元素负荷量表现不同于氮负荷的主要原因。从不同坡度等级氮磷输出占比来看，占比最高的主要分布在坡度范围是7.16～18.76和0～7.16区域内，0～18.76坡度范围氮磷输出占总量的百分比依次是72.83%、78.15%，要占到全流域总流失量的绝大部分。

图8 不同坡度范围总氮总磷输出负荷及总量占比Fig.8 The yield of TN and TP and proportion of different slope

4 结 语

(1)本研究使用SWAT模型分析了阿什河流域氮磷输出负荷与下垫面特征的关系，结果表明，率定期径流的决定系数和纳什效率系数分别为0.62、0.70,验证期分别为0.60、0.67；总氮的率定期决定系数和纳什效率系数分别为0.50、0.67，验证期为0.49、0.75；总磷决定系数和纳什效率系数分别为0.46、0.49，验证期分别为0.45、0.47。除总磷部分年份模拟有一定偏差外，总体模拟结果较好，模型适用于流域非点源污染模拟。

(2)在空间尺度上，流域内总氮总磷输出负荷空间分布较为相似。总氮总磷单位输出负荷较高的区域分布在流域中下游主河道附近。在进行流域非点源污染的治理时，可首先对这部分敏感区域进行集中治理，可在一定程度上降低治理难度和提高治理成效。

(3)流域内不同土地利用氮输出负荷最高的是农业用地、其次为森林、草地等，而对于磷输出负荷，最大也为农业用地，之后则依次为草地、灌丛、湿地等；不同土壤类型单位输出负荷最高的是城区，然后为水稻土、黑土、草甸土等，不同土壤类型总磷输出负荷排序与总氮类似；研究区不同坡度范围的氮磷输出负荷随坡度的增加而呈下降趋势。

□