基于SWAT模型的布尔哈通河流域面源污染的变化研究

石金昊，朱卫红，田 乐，金日

（延边大学 地理与海洋科学学院/长白山湿地生态系统功能与生态安全重点实验室，吉林 延吉133002）

0 引 言

【研究意义】据2010年《第一次全国污染源普査公报》数据显示，面源污染已经成为我国水体的主要污染源，是我国湖泊水质恶化的主要原因。面源污染具有随机性、复杂的迁移转化机理、时空变化幅度大等特点[1]。面源污染是流域不同土地利用类型对降雨径流的综合响应，其产生、输出、迁移和转化与土地利用类型及其不同组合与空间格局特征密切相关[2]，流域景观的组成格局和空间格局变化也会影响面源污染物的产生、迁移和转化。降雨对地表的冲刷通过陆地汇流将污染物带入水体，土地利用和土地覆盖方式影响了地表各类污染物的分布和负荷，用地类型的改变也影响了天然的水文过程。因此，运用景观生态学原理探究景观格局变化与面源污染之间的关系成为面源污染控制的关键。

【研究进展】随着水文模型发展和应用的不断成熟，其中SWAT 模型被越来越多地用来评估流域尺度下的面源污染，分析其时空分布进而识别关键污染区域和关键污染期，也用于分析和评价污染控制管理措施对水环境的影响[3]。Uriarte 等[4]研究了1977—2011年流域土地利用变化对水质的影响，发现河流水质恶化是因为城市化和牧场的发展。Kibena 等[5]发现城市化进程的加快和农业用地的增加，使该地区的面源污染加剧。而李怀恩等[6]、宋林旭等[7]运用SWAT 模型模拟了多年径流、泥沙、氮磷等污染物浓度，表明不同的土地利用条件下产出的氮磷污染物浓度有较大差异，且植被覆盖的增加可以很好地降低面源污染负荷。Wang 等[8]、孙丽娜等[9]利用SWAT 模型和CLUE-S 模型相结合的方法，分析了长时期土地利用变化下面源污染的输出特征，并结合2 种模型模拟未来不同土地利用情景下的面源污染特征。【切入点】上述研究通常集中在流域内的土地利用类型的构成如何解释面源污染的变化，很少有学者研究土地利用模式（景观格局）的空间配置，例如斑块的大小、形状，边缘配置或景观内土地利用类型的斑块的空间连通性等。但土地利用类型指标只能粗略预测水环境的变化，大多数研究没有区分不同的景观格局与面源污染的关系，所以难以全面反映土地利用变化对面源污染的影响。【拟解决的关键问题】因此，基于SWAT模型分析布尔哈通河流域1986、1996、2006、2016年的面源污染时空分布特征，并结合实测的水文水质数据进行率定及验证。基于所构建的布尔哈通河流域SWAT 模型，将重分类后的土地利用数据代入SWAT模型，分别模拟1986、1996、2006 和2016年4 期土地利用现状下的面源污染时空分布，计算1986年和2016年布尔哈通河各子流域的景观格局指数，分析子流域面源污染和景观格局之间的关系，以期为该流域有效减少非点源污染提供科学依据。

1 研究区概况

布尔哈通河位于42°27′N—43°23′N，129°46′E—129°38′E 之间，是国际性河流图们江的重要支流之一（图1），源于安图县哈尔巴岭山脉东南麓，横穿延边朝鲜族自治州中部平原盆地，最终汇入嘎呀河。河流全长约172 km，流域面积约7 064 km2，多年平均径流量约63 198 万m3[10]。全流域内有耕地11.3 万hm2，耕地率16%，其中水田2.85 万hm2，占耕地面积的25.2%，且多集中在布尔哈通河的河谷盆地，是吉林省著名的水稻产区。由于开发历史较早，垦殖率高达20%～30%，河谷内的土地大多被开发成水田和建筑用地，河谷二侧坡地开垦的旱田、果园非常普遍，原有的森林资源面积急剧减少，破碎化程度不断加剧，水土流失面积占流域总面积的39%～49%[11]，加重流域面源污染风险。

图1 布尔哈通河流域DEMFig.1 DEM of Burhatong River Basin

图2 布尔哈通河流域1986—2016年4 期土地利用Fig.2 The land use map for the fourth period of 1986—2016 in the Burhatong River Basin

2 数据与方法

2.1 数据来源

①土地利用数据：本研究基于1986、1996、2006年和2016年布尔哈通河流域的Landsat系列卫星遥感影像，先对遥感影像数据做处理工作，主要包括几何校正、波段融合、波段合成、影像拼接与裁剪等。在ArcGIS10.1 和Cognition902 软件的利用下，运用面向对象分类方法结合人工目视解译方法进行土地利用信息提取工作。从TM 图像中把土地利用类型分为7 类：农田、林地、草地、建设用地（住宅、商业和工业用地）、湿地、水域、裸地（图2）；②为了划分研究所需要的流域，利用NASA 提供的ASTER DEM（分辨率30 m）；③土壤数据来源于寒区旱区科学数据中心，基于世界土壤数据库中国土壤数据集（V1.1）；④气象数据来自中科院水利水电科学研究院的SWAT 模型中国大气同化驱动数据集[12]（CMADS，http://www.cmads.org/），时间尺度为2008—2017年，数据包含了SWAT 模型所需的降水、气温、风速、太阳辐射、相对湿度；⑤水文数据来自榆树川站点月径流数据，时间尺度2010—2013年。月水质监测数据（总磷），时间尺度2011—2013年，该研究仅模拟流域总磷的时空分布。

2.2 分析方法

2.2.1 流域面源污染模拟

利用上述数据构建SWAT 模型后，采用SWAT-CUP 进行校准和不确定性分析，以月径流量、月总磷量为目标对模型进行校准和验证，为判断模型的模拟效果，分别选取决定系数R2和纳什效率系数NSE评价模型精度[13]。校准模型之后，将1986、1996、2006年的土地利用数据分别带入模型，模拟4 期土地利用下的流域面源污染特征。

2.2.2 景观格局分析

将1986年和2016年的土地利用数据由矢量数据转化为栅格数据，然后按照SWAT 划分出的子流域分别裁剪成44 份栅格数据，最后输入Fragstats 4.2 软件运行，计算景观格局指数。本研究在斑块类型水平上选取PLAND（Percentage of Landscape）和PD（Patch Densit）来表征景观格局的组成特征，这2 个指数是确定子流域优势斑块和分析子流域不同类型斑块破碎度的重要依据。在景观水平上选取表征斑块形状复杂程度的LSI指数、表征斑块聚集程度和连通性的CONTAG指数及COHESION指数、表征子流域景观异质性的SHDI指数，这几个简单且通用的景观度量，能比较全面地反映子流域景观空间特征。

2.2.3 统计分析

冗余分析（RDA）现在被广泛用于确定环境因素和生态指标之间的关系，以44 个子流域为样本，使用CANOCO5.0软件分析面源污染负荷和景观格局的关系，其中箭头的长度代表该变量被排序图解释的程度，箭头越长影响程度越高，当各景观指标箭头与水质指标箭头之间夹角小于90°，二者关系为正相关，大于90°为负相关，等于90°则不存在相关性[14]。

3 结果与分析

3.1 流域景观格局变化

林地作为布尔哈通河流域的主要土地利用类型，2016年的林地面积占流域面积的73%以上，并且从1986—2016年林地的面积有着小幅度的增加。由于开垦时间较早，农业用地在流域内的面积占比仅次于林地，在1986年达到21.45%，但到2016年缩减到了20.57%。受城市化进程和经济发展的需要，建设用地从1986—2016年呈不断扩张的趋势，其面积占比增长了0.92%，其余土地利用类型，仅有小幅变动。从1986—2016年，面积占比最大的林地有5 503.45 hm2被开垦成了农田（表1），而大量的湿地被侵占成了农田。因为2002 退耕还林工程在吉林省全面启动，到了2016年，共有6 718.72 hm2的农田转换为林地，其次有4 577.50 hm2的农田变为了建设用地。1986—2016年，共有1 182.28 hm2的建设用地转换成了农田，农田成为建设用地的主要转变方向。

表1 1986—2016年土地利用转移矩阵Table 1 1986—2016 Land use transfer matrix

本研究主要选取了林地、农田、建设用地这3 种面积占比前3 的用地类型来分析其斑块破碎化程度的改变，PD是表征流域斑块破碎化程度的重要指标，由表2 可知，1986—2016年，建设用地的PD平均值基本保持平衡，林地和农田的PD平均值有着上升的趋势，其中农田的PD平均值上升幅度较大，破碎化程度较高。景观水平上表征子流域斑块聚集程度的CONTAG指数平均值呈下降趋势，而斑块连通性指数COHESION平均值基本维持不变。LSI和SHDI分别表征流域斑块形状复杂程度和景观异质性的指数，1986—2016年，二者的平均值都呈上升趋势。

表2 1986—2016 布尔哈通河景观格局变化总体趋势Table 2 General trends of land use and landscape pattern changes in Burhatong River from 1986 to 2016

注 n 表示景观中斑块类型的数量。

3.2 流域面源污染时空分布特征

3.2.1 模型校准及验证

本研究通过SWAT-CUP 2012 进行参数的全局敏感性分析，最后筛选出敏感程度较大的23 个参数，利用SWAT-CUP 程序里的SUFI-2 参数估计最优化方法对榆树川站点2010—2013年的月径流量和月总磷量进行校准及验证。径流在校准期的R2和NSE 分别达0.86、0.81，在验证期的R2和NES 分别达0.72、0.72。总磷在校准期的R2和NSE 分别达0.62、0.57，在验证期的R2和NSE 分别达0.71、0.69。一般情况下，R2＞0.6，NSE＞0.5 时，模型能有较好的模拟效果。综合来看，本研究构建的布尔哈通河流域SWAT模型达到了模拟精度的要求，可以很好地模拟流域水文水质状况。

图3 1986—2016年4 期土地利用下子流域面源污染分布特征Fig.3 Distribution characteristics of non-point source pollution in each sub-basin under land use from 1986 to 2016

3.2.2 不同土地利用的面源污染时空变化特征

时间尺度上看，流域总磷流失负荷峰值主要发生在6—8月，降雨对总磷负荷的作用方式主要通过对径流和输沙量的影响产生，一般情况下降雨量和总磷流失量呈正相关。总磷负荷与径流量大致一致，表现为总磷负荷随径流量增加而增加。已有研究[15]表明，在夏季降雨量多，总磷负荷量也随之增多，汛期面源污染负荷量通常较高。空间尺度上看（图3），总磷的流失主要发生在流域北部和东部地区，集中在9、10、32、33、34 号子流域，在2016年，这5 个子流域的总磷流失负荷占总负荷的30%。1986、1996、2006、2016年的总磷负荷分别为74.04、73.78、82.50、128.31 t。1986—1996年，总磷流失量保持在相对稳定的状态，但从2006年以后，总磷流失负荷开始呈急剧上升趋势，这一差异的原因可能是进入快速城市化发展阶段后，城市的快速扩张和人类生活生产污水的排放直接加速了总磷流失。分析结果表明，建设用地和农田都与总磷流失负荷正相关（图4（b））。

3.2.3 景观格局与面源污染的关系

通过RDA 排序分析探讨1986年和2016年景观格局与非点源污染过程之间的关系。由图4 所知，林地与TP负相关，而建设用地和农田与TP正相关。1986年，景观指标对总磷污染的解释率高达86.43%（表3），总磷流失负荷贡献度较高的景观指标主要是农田、PD（林地）及PD（农田）（图4（a））。而在2016年，景观指标对总磷污染的解释率为77.60%（表3），总磷流失负荷贡献度排前二的景观指标变成了建设用地和PD（林地）（图4（b））。1986—2016年，农田的PD值呈上升趋势，并且农田的PD值与总磷负相关，表明农田斑块的破碎化在一定程度上减缓了面源污染的流失。林地在布尔哈通流域的占比达到了73%，尽管林地对改善面源污染有着积极作用，但林地的PD值与总磷正相关，因为林地斑块的破碎化会造成水土流失，从而加大面源污染流失风险。CONTAG与TP负相关，这表明破碎化严重且景观斑块较为分散的区域面源污染负荷更高。COHESIO与TP负相关，表明水质恶化更可能与流域内分散分布的土地利用有关，所以提高流域内土地利用斑块之间的连通性，减少斑块破碎化对养分流失改善有积极作用。SHDI 和LSI都与总磷正相关，表明流域内景观异质性的增加和形状的复杂化，能够导致污染负荷输出风险水平的增加。

表3 冗余分析各排序轴方差解释率Table 3 Total variance explained by the ordination axis

图4 1986 和2016年污染物流失量与景观指标的冗余分析Fig.4 Redundant analysis of pollutant loss and landscape pattern index in 1986 and 2016

4 讨论

通过构建布尔哈通河流域SWAT 模型发现，流域降雨和径流过程对面源污染流失有着重要影响，总磷流失量和径流量正相关，流域丰水期多出现在6—9月，总磷产生量峰值也主要在6—9月，所以该时期是非点源污染的关键期。农村和城市是非点源污染的主要来源，但在不同的历史发展阶段，二者对面源污染的影响也不同。基于4 期的流域土地利用数据，通过构建SWAT 模型分析因土地利用变化而导致的面源污染时空分布特征，发现受地形、农业生产及城市发展影响的9、10、32、33 和34 号子流域总磷污染负荷较高，属于流域污染的关键控制区。农田和建设用地是面源污染的主要来源，而森林对面源污染改善有着积极影响，这与前人研究的结果基本一致[16]。具有污染输出效应的建设用地和农田一直被视为面源污染的“源”景观，其中建设用地对水质的影响主要取决于流域内城市化水平的高低。大量的生活垃圾、工业污水的排放及农业生产使用的化肥、农药等会经过城市不透水表面进入河流。布尔哈通河中游到下游流域二岸分布城市、村屯、耕地等，来自城区的生活污水、农业污水及工业废水对流域面源污染有着负面影响。该流域传统的农业活动（即耕作或过度施肥）增加了土壤侵蚀的风险，并且过量的养分应用可能导致非点源污染经地表径流后流入到河流。此外，通过对流域的景观格局的进一步分析，发现了农田对面源污染贡献度降低的可能原因。有研究表明，农田斑块的破碎化程度越低，斑块廊道联通性越好，农田污染物容易形成面源污染[17]。尽管33 号子流域的建设用地面积最大，但33 号子流域的污染输出负荷比32 号和34 号小。已有研究表明，流域“源”和“汇”景观的空间位置不同，其对水环境的影响也不同[18]。33号子流域的建设用地和农田基本位于上游区域，下游区域则是大面积的森林，这对该流域的面源污染起到了更好的净化作用。

森林一直被认为与农业面源污染有着直接的联系，其植被覆盖可以减少土壤侵蚀，具有削减暴雨径流、减少水土流失、吸附污染物的作用，可以有效减少地表径流携带营养盐进入河流，对水质恶化在一定程度上有削减作用[19]，但往往局限于景观组成格局的角度分析得出该结论。以具有污染物消减效应的森林为优势斑块的流域景观格局中，其对景观破碎程度最为敏感[20]。林地景观越破碎，斑块越分散，越不利于抑制土壤侵蚀、径流等过程，从而造成污染元素的流失。本研究也通过冗余分析发现林地景观斑块破碎化的增加，会增大流域面源污染产生的风险，而林地作为布尔哈通河流域的主要用地类型，其景观空间格局变化引起的生态变化，更值得注意和警惕。

区域水质是多尺度环境因子的综合反映，不单单受土地利用的影响，景观格局的变化也被认为是影响区域水质的主要原因[21]，景观格局与流域非点源污染研究是探究景观格局与地表环境变化过程的典型模式。景观形状指数（LSI）、香浓多样性指数（SHDI）与总磷污染呈正相关，主要因为这些指标反映了子流域中斑块的多样性，这些多样性的增大又间接地体现了景观的异质性与人类活动干扰的增强，因此也增加了面源污染的风险，其他相关研究结果也表明这些指标与水质的恶化有密切联系[22]。斑块聚集度指数（CONTAG）反映的是景观的分离与散布程度，其高值反映了高度聚集。例如，当不同的土地利用最大限度地聚集时，CONTAG接近100，而当土地利用最大限度地分散时接近与0。斑块连通性指数（COHESIO）表示斑块内细胞的空间连通性和流域内土地利用的物理连通性，高蔓延度值和高连通性说明景观中的某种优势斑块类型高达集聚并且形成了良好的连接性，景观的破碎化程度较低，所以面源污染风险也会随之减小。有相关研究表明坡度主要通过影响产流量和产沙量进而影响磷素的流失总量[23]，而在布尔哈通河流域坡度较高区域，在雨水充沛且多暴雨出现的条件下，会加重水土流失，致使大量未被植物吸收利用的氮磷营养元素随水土流失进入河流，是造成流域面源污染的重要原因。

5 结论

1）SWAT 模型能较好地模拟布尔哈通河流域总磷污染时空分布，反映研究区的总磷污染流失负荷。

2）布尔哈通河流域面源污染防治的关键控制区域在坡度较高的北部区域和农田和城镇集中的东部区域。

3）通过冗余分析发现林地与总磷负荷负相关，农田和建设用地与总磷负荷正相关；林地斑块破碎化越高，流域面源污染流失越严重，农田斑块的破碎化在一定程度上减缓了面源污染的扩散。流域斑块之间的高度流通性和聚集性，对养分流失改善有积极作用。流域内景观异质性的增加和形状的复杂化，能够导致污染负荷输出风险的增加。对该区域进行土地利用规划与景观格局优化时，应从景观组成格局和空间格局二方面综合考虑，才能更好地改善流域面源污染。