韩江流域面源污染及与景观格局的关系

郑 宇，程香菊，王兆礼,2，赖成光,2

(1.华南理工大学土木与交通学院，广东 广州 510641； 2.华南理工大学亚热带建筑科学国家重点实验室，广东 广州 510641)

面源污染导致河流湖泊等水体出现富营养化和生物多样性损失，并威胁人们的用水安全[1-2]。我国面源污染尤其是农业面源污染已经成为一个重要的水环境问题[3-4]。面源污染的产生和发展与地形、气象、土壤、土地等诸多因素有关[5-6]，而流域景观的组成格局和空间格局变化会改变这些因素，进而影响面源污染物的产生、迁移和转化[7-9]。因此，以景观生态学原理为指导，研究流域面源污染与景观变化之间的关系，有助于对景观的组成格局和空间格局进行有效的规划和管理，从而改善流域的面源污染。近年来，流域水文水质模型结合数理统计分析的方法被越来越多地应用于研究面源污染与景观格局的关系[10]，但大部分研究重点关注景观的组成格局(例如土地利用面积比例)与面源污染的关系，而忽视了景观的空间格局(例如景观的斑块密度、形状指数等)与面源污染的关系，因而难以全面反映景观对面源污染的影响[11]。

韩江流域是广东省内的第二大流域，从20世纪80年代开始，随着经济的发展，流域内的景观发生了较大变化，对生态环境造成了深刻影响。韩江的上游梅江水质日趋恶化，在某些河段存在富营养化现象，部分省界断面水质已经严重恶化，威胁下游潮州、汕头等市人民的用水安全。本文基于水土评价工具(soil and water assessment tool, SWAT)模型构建面源污染模型，定量模拟和分析韩江流域面源污染的时空分布特征，识别关键污染地区，并应用多元线性回归、冗余分析等方法，探讨流域面源污染与景观格局之间的关系，以期为流域水环境管理提供依据。

1 数据与方法

1.1 研究区概况

以韩江流域潮安站以上部分为研究区，在SWAT模型中处理划分为48个子流域(图1)，总面积约29 077 km2。韩江流域拥有800多万人口，流域内的多年平均降雨量为1 600 mm，主要土地利用类型是林地。流域内地形变化较大，高程在20～1 500 m之间，坡度在11.7°～20.0°之间。近年来，由于城市化的快速发展和水利工程的建设，流域的水质逐渐恶化[12]。

1.2 数据来源

本文所用的主要数据见表1，其中，依照全国土地利用分类标准，土地利用类型数据被重新划分为水田、旱地、草地、林地和城镇5个主要类型；气象数据包括日平均降雨、最高和最低气温、相对湿度、风速及日照时长等；水文水质数据包括径流量、含沙量、氨氮质量浓度等的监测数据。除了表1所列数据，还根据研究区内12个县的统计年鉴等资料统计了各县的人口、禽畜与水产养殖数量、耕地面积和化肥用量等数据，按照第一次全国污染源普查得到的排污系数成果，估算出研究区的生活污水、水产养殖、禽畜养殖和无机化肥施用产生的污染负荷并添加到模型中，以提高模拟的精度。模型中的农业管理措施通过走访和文献调研获得，例如研究区的主要作物被设定为双季水稻和花生，早稻于每年的4月种植，7月收割；晚稻于每年的8月种植，11月收割；花生于每年的3月种植，7月收割。

图1 韩江流域的地理位置、水系、站点与子流域分布

1.3 流域面源污染模拟

SWAT是美国农业部开发的流域尺度的时间上连续、空间上呈半分布式的模型，可用于模拟土地利用、土地管理措施、气候变化等因素对流域水文水质的影响。本文利用地形、气象、土地利用、土壤、管理措施等数据构建SWAT模型，用SWAT-CUP(SWAT calibration/uncertainty program)校准和不确定性分析程序，以径流量、含沙量和氨氮质量浓度为目标对所建模型进行校准和验证，并选用决定系数R2、纳什系数ENS和相对误差ER等指标来评价模型的拟合效果[13]，选用fp因子和fr因子来评价模型的不确定性[14]。用校准后的模型分别加载3期的土地利用数据，输出各子流域出口处1990—2014年的总氮、总磷、有机氮、有机磷、氨氮、无机磷等污染物的多年平均质量浓度，以进行进一步的分析。

表1 研究采用的主要数据

1.4 景观格局分析

景观格局可以分为组成格局和空间格局，并用相应的景观指标来衡量[15-16]。本文子流域景观格局的组成特征用某一土地利用类型面积百分比PLS和最大斑块指数PLP来衡量，这两个指标是确定景观中优势元素的重要依据，也是决定景观中生物多样性和优势物种的重要因素。景观格局的空间特征用某一土地利用类型斑块密度DP和形状指数ILS来衡量，其中斑块密度反映了景观的破碎程度，形状指数反映了景观的形状复杂度。这些指标的计算公式如下：

(1)

(2)

(3)

(4)

式中：n为子流域中某种土地利用类型的斑块数量；ai为斑块i面积，m2；A为子流域面积，m2；ei为斑块i的边缘长度，m。

1.5 统计分析

在RStudio中，利用多元线性回归和冗余分析研究景观格局与面源污染的关系。多元线性回归能够反映单个因变量与多个自变量之间的关系，而冗余分析是多元线性回归的扩展，能够同时反映所有因变量与自变量间的关系。由于自变量数目较多，本文在多元线性回归中使用全子集回归对自变量集进行缩减，并根据回归诊断的结果选择最优的回归方程；在冗余分析中使用向前选择对自变量集进行缩减。

2 结果与讨论

2.1 流域的景观格局特征

在1990年和2010年土地利用背景下，整个流域内各土地利用类型的面积占流域总面积的百分比如表2所示。和其他土地利用类型相比，尽管城镇用地的面积比例是最小的，但是增长率却是正值且是最大的，说明研究区正在经历城市化的初始阶段，各种土地有转换为城镇用地的趋势。在2010年土地利用背景下，部分典型子流域中各土地利用类型的面积百分比PLS、最大斑块指数PLP、斑块密度DP、形状指数ILS如表3所示。可以发现，林地在每个子流域中都拥有最高的面积百分比和最大斑块指数，说明研究区内的土地利用类型以林地为主。人口密集的子流域，例如兴宁盆地附近的33、36、37、46号子流域，农业和城镇用地的面积百分比更高。水田的面积百分比一般大于旱地，说明流域内的主要农业活动是水稻种植。上游子流域的草地面积百分比普遍大于下游子流域。尽管草地的面积百分比一般大于水田，水田的最大斑块指数却常常大于草地，说明流域内草地的分布相比水田更为分散。城镇的斑块密度和形状指数相比其他土地利用类型更低，表明在人类活动的影响下，城镇斑块的破碎程度更小，形状更规则。

表2 各土地利用类型面积占全流域面积的百分比 %

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

在SWAT-CUP软件中通过全局敏感性分析筛选出较为敏感的重要参数，例如地表径流滞时、基流退水常数、曼宁系数、径流曲线数、不同土地利用类型的流失系数等，再用内置的SUFI-2(sequential uncertainty fitting version 2)算法进行校准和验证，结果如表4所示。在校准期和验证期，径流量、含沙量、氨氮质量浓度的决定系数R2和纳什效率系数ENS基本都超过0.7，相对误差ER基本都小于8.5%，说明模型的拟合效果良好。所有站点的fp因子都大于0.5，同时fr因子都比较小，表明模型具有较小的不确定性。因此，所建立的模型基本上能反映流域的面源污染规律。

应用SWAT模型，分别加载1990年和2010年的土地利用数据进行模拟，并输出各子流域出口处的总氮、总磷质量浓度以及它们的增量，结果如图2、图3所示。从空间上看，各子流域的污染负荷相差较大，总氮、总磷质量浓度呈现明显的空间差异性。位于西部宁江流域地区(28号及其上游的子流域)的总氮质量浓度较高，并在兴宁市及附近地区(36号子流域)达到最大值3.91 mg/L，这主要是因为该地区位于农业较为发达的宁江平原，农田径流、禽畜粪便、生活污水的产生量较大，产生了较大的总氮流失。流域东部的汀江流域地区的总磷质量浓度普遍较高，并在龙岩市永定区及附近地区(21号子流域附近)达到最大值2.15 mg/L，这可能是因为这些地区的草地和林地面积比重普遍较大，而草地相比其他土地利用类型具有更高的磷流失风险[17]。从1990—2010年，大部分子流域的总氮负荷呈现增加趋势。虽然在大部分子流域，农田、城镇等人类用地的面积比自然土地面积小，但是在整个子流域的出口处(45号子流域)，各污染物的质量浓度仍然较高且呈现上升趋势，这反映出流域内人类活动对水环境的影响大于自然过程。

表3 部分子流域2010年各土地利用类型的景观格局指标

表4 SWAT模型的拟合效果评价

(a) 1990年

(b) 2010年

(c) 增加量

图2 各子流域的总氮质量浓度模拟值

(a) 1990年

(b) 2010年

(c) 增加量

图3 各子流域的总磷质量浓度模拟值

2.3 流域的景观格局变化与面源污染变化的关系

以1990—2010年土地利用背景下各污染物质量浓度的增量为因变量，各景观格局指标的增量为自变量，建立多元线性回归方程如下：

yTN=0.536xDP5+0.093xPLP2

(5)

yON=-0.132xDP1+0.056xPLP2-0.003xILS5

(6)

yAN=-0.016xPLP5-0.015xILS3+
0.014xILS1-0.002xPLS3

(7)

yTP=-0.233xDP1+0.079xPLP2-0.005xILS5

(8)

yOP=-0.184xDP1+0.071xPLP2-0.004xILS5

(9)

yMP=-0.073xDP1+ 0.009xPLP2-0.002xPLS3

(10)

式中：yTN、yON、yAN、yTP、yOP、yMP分别为总氮、有机氮、氨氮、总磷、有机磷、无机磷质量浓度的增量，mg/L；xPLS、xPLP、xDP、xILS分别为面积百分比、最大斑块指数、斑块密度、形状指数的增量；数字下标1、2、3、4、5分别表示水田、旱地、林地、草地和城镇。

经过全子集回归的筛选，进入不同回归方程的自变量不同，说明景观格局对不同污染物的贡献程度各异。旱地最大斑块指数的增量与大部分污染物质量浓度的增量呈现正相关关系，水田斑块密度的增量与大部分污染物质量浓度的增量呈现负相关关系，林地形状指数和面积百分比的增量与氨氮质量浓度的增量呈现轻微的负相关关系。城镇斑块密度的增量与总氮质量浓度的增量呈现强正相关关系，表明城镇斑块可能是氮污染的主要贡献者；相反城镇的形状指数增量与有机氮、总磷、有机磷质量浓度增量呈现轻微的负相关关系；城镇的面积百分比和最大斑块指数没有进入多数回归方程，表明城镇的组成格局对面源污染的影响没有空间格局的影响大。

在多元线性回归的基础上，通过冗余分析来同时反映所有因变量与自变量间的关系。在冗余分析中，因变量的总方差被分解成为约束性和非约束性方差。由表5可见，前两个轴变量已经能解释81.7%的总方差，以这两个变量为横轴和纵轴，将各自变量、因变量和子流域绘制成图4所示的三序图。从图4可以看出，大部分子流域沿R2轴方向分布，而R2轴与城镇的xPLP、xDP和xILS形成的夹角较小，说明流域内各子流域之间的城镇化程度差异较大。进入冗余分析模型中的景观指标主要是城镇和旱地的景观指标，这说明污染物质量浓度的变化对城镇和农田景观指标的变化更敏感。城镇的最大斑块指数增量、斑块密度增量和形状指数增量以及旱地的形状指数增量、林地的面积百分比增量与各污染物质量浓度的增量呈现负相关关系，而旱地和草地的面积百分比增量与各污染物质量浓度的增量呈现正相关关系。

表5 冗余分析中各轴变量解释的方差比例 %

2.4 讨 论

相关研究表明，土地利用类型可以被划分成“源”景观和“汇”景观，前者对流域内的泥沙流失、污染物流失具有较大贡献[18]，而后者对泥沙流失和营养物流失的贡献较小甚至有削弱作用[19]。林地通常被认为是一种“汇”景观格局，因为它具有保持水土、净化水质的功能[20-21]。图4亦表明，林地的面积百分比与大部分污染物的质量浓度负相关。

由于削减径流、反硝化等作用，草地也常被认为是一种“汇”景观格局[22-23]，但也有研究指出草地增加了水环境污染的风险[24-25]。图4表明，草地面积百分比与污染物质量浓度呈正相关关系，说明研究区内的草地更像一种“源”景观格局。这主要是由以下原因导致的：①研究区的坡度较大，平均坡度达到16.6°，增加了营养物随土壤流失的风险[26-28]；②研究区受亚热带气候影响，降雨量较为充沛，多年平均降雨量超过1 400 mm，且80%集中在雨季，增加了土壤颗粒随着地表径流流失的风险[29]；③从1990—2010年，研究区有197.3 km2的林地被转化成草地，占了草地来源的78.5%，而草地对污染物的削减作用比林地弱，从而导致了径流中营养物通量的增加[30]。因此，尽管草地景观具有调节径流和削减污染的作用，依然不足以抗衡大坡度和强降雨的影响。本文得到的结果表明流域管理中考虑坡度的重要性，例如在坡度较大的地区退耕退草还林可能是降低流域面源污染的有效措施。为论证这一结论，设置了将10°以上地区的草地转换为林地的管理情景并进行模拟，结果如表6所示，该情景下流域出口处的多年平均总氮通量削减了1.82%，总磷通量削减了1.16%。

由于化肥的施用以及土壤的松动，农业用地是造成水环境恶化的重要因素之一[31]。式(5)～(10)和图4表明农业用地最大斑块指数的增量与污染物质量浓度的增量呈现正相关关系，但是，与水田相比，进入回归模型和冗余分析模型中的旱地指标更多，说明旱地与污染物质量浓度的相关性比水田与污染物质量浓度的相关性更强。这可能是由以下原因造成的：①农民倾向于在旱地施用更多肥料，导致农田径流中的氮磷等营养物质量浓度更高[32]；②旱地上的作物可能会产生更大的渗漏淋失[33]，导致旱地土壤中的有机质和总氮质量浓度比水田中的更高[34]；③水田具有更长的水力保留时间和缺氧环境，促进了反硝化作用，从而减少了总氮的产生。因此，调整农业结构，例如将旱地改为水田可能有助于减少肥料用量和土壤流失。为了证明这一结论，设置了将兴宁市的旱地改为水田的管理情景并进行模拟，结果如表6所示，该情景下流域出口处的多年平均总氮通量削减了4.63%，总磷通量削减了7.31%。需要注意的是，本文提出的管理措施主要是为了佐证通过数理统计方法得到的面源污染与景观指标之间的关系，因此在设定模拟情景时保持了模型数据库中各种土地利用类型的参数(如污染物去除率、城镇不透水面积等)不变，没有考虑设置过滤带、采取低影响开发措施、改善污水处理工艺、设置植草水道等具体工况的效果，导致模拟结果中对污染物的削减量有限；另外，虽然考虑到措施的可行性，只对局部地区施行了改造措施，但是没有进行详细的经济效益评估。如何结合本文关于面源污染与景观格局关系的结论，提出更加立体、有效、经济的流域管理方案，是一个有待进一步研究的课题。

除了组成格局，景观的空间格局也与水质相关。例如，水田的斑块密度与有机氮和各种形态磷的质量浓度负相关，旱地的形状指数与大部分污染物的质量浓度负相关。本文和前人的研究结果都表明，城镇景观的斑块密度和形状指数与各种形态的氮磷污染物质量浓度负相关[35]，这主要是因为城镇景观受人类活动影响较大，因而拥有更密集的空间分布和更规则的形状[36]。

式(5)～(10)和图4表明林地的空间格局指标与污染物质量浓度的相关性没有组成格局指标与污染物质量浓度的相关性强，而城镇景观则相反。考虑到林地和城镇分别是研究区每个子流域中面积最大和最小的土地利用类型，推测景观的组成格局和空间格局对水质的贡献程度与景观的面积有关。面积大的景观如林地和草地，其空间格局与水质的相关性更强；面积小的景观如农田和城镇，其组成格局与水质的相关性更强。

表6 不同管理情景下流域出口处污染物的多年平均通量和质量浓度

3 结 论

a. 研究区的土地利用以林地、草地为主，水田、旱地次之。流域处于城镇化的初始阶段，尽管城镇面积比例较小，但增长率最高。下游子流域的城镇化程度比上游高，呈现出向上游蔓延的趋势。在以后的研究中模拟未来土地利用情况或者进行流域管理规划时，应当重视城镇的影响。

b. 在当前土地利用策略和人类活动的影响下，研究区内的水质呈现恶化趋势。研究区西部的宁江流域地区是总氮污染的关键区，总氮质量浓度最大3.91 mg/L；东部龙岩市永定区的总磷质量浓度较高，最大值为2.15 mg/L。

c. 流域的面源污染与景观格局存在密切的关系。景观指标能解释81.7%的面源污染变化。景观组成格局指标与空间格局指标的作用与景观面积有关，景观面积越大，其组成格局可能比空间格局对水质的影响更大。林地的面积百分比与面源污染呈现负相关关系，其他景观则相反。城镇或农田等人类景观的斑块密度、形状指数等空间格局指标与污染物质量浓度具有较强的负相关性。旱地相比水田对研究区面源污染的贡献更大。

d. 在坡度较大的地区退草还林，可以削减1.82%的总氮和1.16%的总磷；将农业结构由旱地为主调整到以水田为主，可以削减4.63%的总氮和7.31%的总磷。