氮磷净化能力是水质净化功能的重要体现, 对制定水源保护规划, 落实陆域-水域综合保护, 防治水体污染, 保障区域用水安全具有重要意义^[1-2]。丹江口库区作为南水北调中线工程重要水源地, 其水质状况对工程具有举足轻重的作用, 不仅关系到库区及沿线地区人民生活、生产用水安全, 更影响着周边地区社会和经济的发展。

目前常用的水质评价方法有单因子法、综合指数法、灰色评价法、水质指数法、多变量分析法等^[3], 这些方法多用来评价水质状况, 无法针对水质净化功能进行评估。随着3S等技术的发展, 生态系统服务模型的建立与使用成为研究热点, SWAT模型^[4]、HSPF模型^[5]和InVEST模型^[6]等被开发利用。InVEST模型以其简单、快速和强大的空间表达能力成为应用最广泛的工具, 该模型营养盐输送率(NDR)模块专门用于模拟氮、磷出口的水质净化功能研究, 并已在世界许多地区得到应用。Goldstein等^[7]研究了粮食作物、牧草情景及生物燃料等七种情境下对夏威夷岛屿水质净化功能空间分异情况；王大尚等^[8]利用InVEST模型对密云水库上游流域氮磷水质净化服务在内的诸多生态系统服务功能及居民福祉进行了定量评估和空间特征刻画；吴瑞等^[9]、韩会庆等^[10]利用模型对官厅水库、贵州省珠江流域氮、磷水质净化功能开展研究。

水质净化能力不仅受地形、地貌、降水等自然因素的影响, 人类活动引起的景观格局变化以及城镇化集聚、农业生产等均会影响污染物的迁移和转化, 从而改变进入河流污染物的数量, 对流域水质产生重要影响。研究表明景观格局变化是影响水质净化的关键因子^[11-12], 不同景观要素如景观类型、景观面积及景观特征均对水质净化能力产生显著影响^[13-14]。赵鹏等^[15]运用多变量分析方法, 分析流域、河岸带尺度的景观组成和空间格局指数与水质的关系, 结果表明流域尺度下景观指数对水质影响较大；李昆等^[16]利用3S技术和景观格局指数以7种空间尺度的湖滨带缓冲区为对象, 研究景观格局对洪湖水质变化的空间尺度效应；Oliverira^[17]等选择河岸缓冲区和专属贡献区, 分析了景观格局与五个选定水质指数之间的联系；Rajaei^[18]研究了塔简河流域景观格局与硝酸盐等污染物的相关性, 结果表明流域硝态氮与斑块密度、景观形态指数等景观格局指数呈正相关。王瑛^[19]研究了景观格局与水环境变化的响应, 发现通过人为措施合理改变景观结构对控制面源污染具有显著作用。因此, 探讨景观格局与水质净化能力的关系, 对于揭示库区景观格局变化, 合理规划景观减少面源污染的传递具有重要意义。

丹江口库区是南水北调中线工程重要水源地, 水环境状况不仅威胁当地生态环境和用水安全, 更直接关系到工程受水区的水量和水质。《丹江口汇水区土壤养分状况及农作物施肥量调查研究报告(2015)》显示, 库区污染源中农业面源污染占比高达70%, 库区坡耕地主要经济作物柑橘为面源污染主要来源, 化肥和农药的不合理使用带来较高的氮磷输出, 高强度集中降水带来的水土流失成为面源污染迁移的重要载体, 大量养分进入水体造成水质退化, 由此产生的一系列生态问题为学者所关注。Wang等^[20]研究了丹江口库区景观格局和土地利用变化对水生态时空变化的关系, 但由于数据量较少, 而未能准确反映生态系统服务的变化；张煦等^[21]基于水质理化指标对丹江口库区水质开展监测, 得出农业污染是影响水质安全的重要因素, 丁霞等^[22]对库区马尾松林水源涵养价值进行了研究, 缺乏水质时空动态变化。目前, 该区景观格局变化和氮磷净化能力及二者关系的研究罕有报道。基于此, 在县域尺度上以丹江口市为研究区域, 利用景观格局指数和InVEST模型分析2003—2018年景观格局变化和水质中氮磷净化能力, 并利用Pearson相关分析和RDA冗余分析研究二者的关系, 探讨引起氮磷净化能力变化因素, 以期为丹江口库区水源水质保护提供科学依据。

1 研究区概况

丹江口市位于湖北省西北部、鄂豫两省交界, 是南水北调中线工程水源地, 国家级重点生态功能区县, 素有“中国水都”之称, 水源区位于秦岭巴山之间, 汉江中上游, 市跨32°14′10″—32°58′10″N, 110°47′53″—110°34′47″E之间；属于北亚热带季风区的温暖半湿润气候, 四季分明, 降水分布不均, 年平均气温15.9 ℃, 年降水量约800 mm, 年蒸发量为1979.1 mm；河流以汉江为主, 自西向东横跨全境；土壤以黄棕壤和黄壤为主, 成土母质由石灰岩、片麻岩等发育而成, 质地疏松, 地貌类型多以低山丘陵为主；库区森林类型由人工林和次生林组成, 马尾松(Pinus massoniana)、栓皮栎(Quercus variabilis)、杉木(Cunninghamia lanceolata)等广泛分布；经济以农林牧渔业发展为主, 并广泛种植柑橘(Citrus reticulata Blanco)。

2 研究方法 2.1 数据来源及处理

InVEST模型水质净化模块所需数据包括研究区土地利用/覆被图、降水量、年潜在蒸散量、土壤深度、植物可利用含水量、DEM数据与流域及子流域、生物物理系数表、氮磷输出负荷数据以及流域累积阈值等数据。土地利用/覆被图通过遥感影像图解译得到, 选择2003年Landsat 7 ETM+、2008年Landsat-5 TM、2013年和2018年Landsat 8 OLI_TIRS影像, 分辨率为30 m；在ENVI软件中通过几何校正、融合镶嵌、裁剪处理, 将研究区分为林地、耕地、园地、建设用地、水域和未利用地6类, 见图 1。

降水数据通过中国气象数据网获得, 对研究区域及周边气象站点进行克里金插值分析。年潜在蒸散发利用中国地面气象观测报告, 收集研究区周边4气象站数据, 运用Hargreaves方程^[23]计算年潜在蒸发量。土壤深度根据森林资源二类调查结果, 采用克里金插值得到模型所需栅格图。植被可利用水含量可评估土壤为植被生长所储蓄的总水量, 也可确定田间持水量和萎蔫点的差值^[24], 计算公式选取周文佐等^[25]所提出的AWC模型进行估算, 将计算的结果进行克里金插值得到模型所需栅格图。

(1)

式中, SAN为土壤的沙粒含量(%)；SIL为土壤的粉粒含量(%)；CLA为土壤的黏粒含量(%)；C为土壤的有机质含量(%)。

DEM数据从地理空间数据云网站获取, 经镶嵌、裁剪、填洼处理, 使用水文分析工具和矢量数据生成流域和子流域, 每一级流域赋予唯一数值。生物物理系数表和氮磷输出负荷数据, 参考丹江口库区的调查和文献^{[22, 26-27]}及InVEST模型用户指南获得, 相关参数见表 1。

2.2 氮磷净化能力研究方法

InVEST模型水质净化模块主要基于植被土壤可通过存储、转换的方式截留径流中营养盐污染物的机理, 计算每个栅格中养分元素的量, TN、TP输出量越高表明氮磷净化能力越低, 其主要算法：

(2)

式中, ALV_x是像素x的调整负荷值, HSS_x是像素x计算方法的水文敏感度评分, pol_x是像素x的输出系数。计算得到TN、TP输出量后, 结合土地利用/覆被类型对污染物的移除效率来得到营养物保持(截留)量。

依据氮磷净化能力结果, 利用水域营养类型的划分标准, 结合《生态功能区划技术暂行规程》, 对研究区的净化能力重要性进行分类, 共分为5级(表 2), 以此表现区域氮磷净化能力的强弱。

2.3 景观格局指数选取

景观格局指数能有效反映生态系统与空间格局的关系, 可在斑块水平、斑块类型水平和景观水平3个层次上评估景观格局变化、分析景观与水质变化的关系等^[28-29]。基于前人的基础^[30-31], 本文在土地利用/覆被图的基础上利用ArcGIS和Fragstats 4.2软件, 从斑块类型和景观水平层次开展研究, 选取表征面积、聚集度、形状和多样性测量等8个具有代表性和广泛应用的景观指标(表 3)。

2.4 数据处理与统计分析

采用统计分析软件SPSS 24.0和Canoco 5.0软件, 对不同时期TN、TP输出量与景观格局指数进行Pearson相关分析和冗余分析(Redundancy analysis, RDA), 探究氮磷净化能力与景观格局的相关关系, 利用RDA冗余分析识别各个因子贡献大小, 揭示在景观尺度下氮磷净化能力变化的主要驱动。

3 结果与分析 3.1 丹江口市景观格局特征及变化

选择斑块类型和景观水平对丹江口市的景观格局开展研究, 结果表明(表 4), 丹江口市水域面积增加, 斑块类型比和最大斑块指数总体呈上升趋势, 水域占研究区面积的比例持续上升；林地面积逐年增加, 从40.80%增长至50.84%, 最大斑块指数先增后减表明受人类活动干扰较大；耕地与园地面积逐渐减少, 耕地、园地的斑块类型比和最大斑块指数总体呈下降趋势, 但依然占比较大, 表明当地对耕地和园地的经营管理强度依然较大；建设用地先减后增, 总体新增面积约2553 hm², 城市及镇村的发展和移民搬迁政策执行, 都加剧了建设用地变化；未利用地主要包括裸地等, 基本保持持平状态, 因此合理开发利用裸地将成为下阶段土地利用的重点。

景观水平上, 2003—2018年景观形状指数不断减少, 蔓延度指数和聚集度指数整体呈增长趋势, 说明形状与相同面积的规则图形之间偏离程度逐年减小, 形状复杂程度降低, 斑块类型间形成了良好的连接性, 景观聚集程度提升, 空间分布趋向集中, 景观的异质性减弱；散布与并列指数增加了5.58, 香农多样性指数变幅较小, 说明各斑块类型趋于规则, 呈均衡趋势分布。究其原因, 水库蓄水、移民搬迁、退耕还林和生态文明建设等项目的实施, 社会经济发展以及人类不断活动, 都加速了景观格局的变化。

3.2 氮磷净化能力时空变化与等级划分

利用InVEST模型对氮磷净化能力开展研究, 其中TN、TP输出总量越大氮磷净化能力越差。结果表明(表 5、图 2)TN、TP单位面积输出量最大值均出现在2013年, 分别为1.625 kg hm^-2 a^-1和0.129 kg hm^-2 a^-1, 总输出量分别为776.979 t, 68.163 t；但整体上看2003年TN、TP输出量最高分别为899.224 t、77.308 t, 输出量的最大值分别为1.559 kg hm^-2 a^-1和0.124 kg hm^-2 a^-1；2008年TN、TP的输出量开始降低, 输出量的最大值分别为1.583 kg hm^-2 a^-1和0.125 kg hm^-2 a^-1；2018年TN、TP的输出量最低, 总量为672.139 t、60.802 t, 其中最大值分别为1.588 kg hm^-2 a^-1和0.126 kg hm^-2 a^-1。

从时间上看, 2003—2018年丹江口市TN、TP输出量呈减少趋势, 表明氮磷净化能力逐渐增强, 15年间TN的净化服务增强了25.3%, TP的净化服务增强了21.4%；空间上看, TN、TP输出量高的地区均分布在地势平坦和人口分布密集区。丹江口市农业以种植农作物和柑橘园为主, 农药、化肥的使用是氮、磷的主要来源, 有部分农业区距离水库较近, 因此濒临水库处应扩大林地面积, 进一步增强该区水质净化能力。

以2018年TN、TP模型输出结果为例, 进行氮磷净化能力重要性等级分类, 共分为5级(图 3)。整体上看丹江口市2018年氮磷水质净化重要性等级分布基本相似, 特别是氮磷净化能力强, 极重要区、高度重要区主要分布在库区、河流两岸和林地区域；净化能力稍弱, 一般重要区主要分布在城镇建设、农田、柑橘园等。

综上结果可知：氮磷净化能力与景观格局和土地利用类型有关, 随着时空变化丹江口市近15年氮磷能力不断增强, 较弱区域主要集中在平坦地区, 这些区域植被覆盖度较低, 对污染物的截留效用较弱, 且耕地和柑橘园化肥农药的使用增加了氮磷的输出量；同时这些区域城镇化建设明显, 人口和人类活动的增加也导致污染物排放递增, 从而导致净化能力较弱；净化能力强, 极重要、高度重要区主要分布在库区、河流两岸和林地等区域。

3.3 景观格局变化与氮磷净化能力的关系

由于InVEST模型水质净化模块主要与林地、耕地、园地和建设用地有关, 利用Pearson相关分析结果见表 6。在斑块类型水平上, 林地类型面积与氮磷输出量均呈显著负相关, 表明林地可以有效拦截氮磷进入水库；耕地、园地斑块类型面积和斑块类型比均与氮磷输出量呈正相关关系, 说明农业生产是影响水质的重要因素, 随着农业用地的减少, 污染源的输入逐渐下降, 其中耕地类型的最大斑块指数与氮磷输出呈负相关, 表明耕地最大斑块指数的增加对氮磷净化能力的提升起一定的促进作用；氮磷输出量与建设用地整体上呈中等程度相关, 因此城镇化发展中建设用地的增加也是影响水质变化的重要因素。

在景观水平上, 氮磷输出量与景观形状指数呈显著正相关, 说明TN、TP净化能力与景观形状的复杂度有关；与香农多样性指数呈中等程度正相关；氮磷输出量和景观聚集与分散度测量指标中蔓延度指数、聚集度指数呈负相关, 和散布与并列指数表现出显著负相关, 表明随着各类景观的连通性与聚集性的升高氮磷净化能力增强。

在景观格局与氮磷净化能力相关性分析的基础上, 通过Canoco 5.0软件进行解释变量的前向选择, 排序轴的特征值分别为0.843和0.109, 排序图可以很好地反映景观与氮磷净化之间的关系图 4表明氮磷输出量与林地类型, 蔓延度指数、散布与并列指数、聚集度指数呈负相关；与耕地、园地、建设用地类型, 景观形状指数、香农多样性指数呈正相关, 这与Pearson相关分析结果一致。由箭头连线的长短可看出土地利用类型依然以林地、耕地和园地对氮磷净化能力影响最大, 各景观格局指数带来的影响较为均衡。

4 讨论与结论 4.1 讨论

南水北调中线工程通水后, 为确保水质安全, 协调库区的区域社会经济可持续发展和水环境资源保护, 国务院先后批复了《丹江口库区及上游水污染防治和水土保持规范(2006年)》、《丹江口库区及上游水污染防治和水土保持“十二五”规划(2013年)》以及《丹江口库区及上游水污染防治和水土保持“十三五”规划(2017年)》, 随着规划的推行和实施, 面源污染氮磷等元素的输出、水质变化规律和变化驱动因子成为丹江口库区水环境研究的重点。孙玉君等^[32]采用单因子评价法, 对丹江口库区入库水系、水功能区以及水质指标进行了监测和分析；周颖^[33]借助GIS和RS等技术, 对丹江口库区胡家山小流域等流域水质进行监测, 探讨导致面源污染输出季节性差异的驱动因子, 得出林地、耕地和建设用地是影响小流域水质的主要因子；张博等^[34]研究发现丹江口库区林地对面源污染的输出有拦截和削减作用, 耕地则是水体氮磷污染的主要来源；丁飞霞^[35]以张家山小流域不同时期为研究对象, 分别从景观数量特征和景观格局特征探讨汇水区景观特征对小流域径流水质的影响。前人研究结果与本文研究结果相似, 且本文利用InVEST模型和景观格局结合的方法完善了丹江口库区氮磷净化时空变化研究, 探讨影响水质净化功能的因素。

水质净化功能和景观格局变化的相关性研究主要分为动态和静态的研究, 即通过时间和空间表征二者关系。朱兰保等^[36]基于GIS和RS利用多年的水质数据与景观格局进行分析, 研究发现景观类型是影响水质的重要因素, 但未明晰水质变化与景观格局变化的关系；刘怡娜等^[37]对长江流域景观格局变化与水质净化服务的关系开展研究, 表明河岸带景观格局的改变更能影响污染物的截留效率；伍恒赟等^[38]通过汇水单元和子流域划分, 研究不同区域景观格局变化与水质的联系, 从空间上分析了两者的相关性, 但时间期限较短。本文选择以时间变化序列研究方法, 从时空动态角度探讨水质净化中氮磷净化能力与景观格局的关系, 整体层面研究了丹江口市近15年氮磷净化能力、景观格局变化以及两者间的相关性。在斑块水平上选择林地、耕地、园地和建设用地, 着重探讨4种土地利用类型带来的影响, 在景观水平上选择聚集度、形状和多样性测量等指标, 进一步探究了各指标与水质净化功能的关系, 也表明了营养物质累积浓度的空间分布在不同景观类型上存在差异, 补充了丹江口市时空变化下景观格局与氮磷净化能力关系的空缺。

由于本文利用InVEST模型方法开展研究, 水质净化模块阈值累积量会随着季节性产生差异, 可能会引起氮、磷输出量存在误差；模型下一步应增加对其他污染物输出量的研究。

4.2 结论

丹江口市2003—2018年景观格局变化明显, 水域面积和林地面积不断增加, 耕地和园地面积逐渐减少, 建设用地面积整体上呈增长趋势；景观形状指数不断减少, 蔓延度指数和聚集度指数整体呈增长趋势, 形状复杂程度降低, 景观聚集程度提升, 散布与并列指数增加了5.58, 香农多样性指数变幅较小, 各斑块类型趋于规则, 呈均衡趋势分布。2003年、2008年、2013年和2018年氮磷净化能力呈持续增强趋势, 其中TN输出总量分别为899.224、801.481、776.979、672.149 t, TP输出总量分别为77.308、69.921、68.163、60.802 t, 15年间TN、TP的水质净化能力分别增强了25.3%和21.4%；对氮磷净化能力重要性等级划分发现极重要、高度重要区主要集中在库区、河流两岸和林地区域；氮磷净化能力稍弱, 一般重要区主要在城镇建设、农田、柑橘园。利用Pearson相关分析和RDA分析发现, 在斑块类型上林地面积与氮磷的输出呈显著负相关, 耕地、园地和建设用地与氮磷输出量呈正相关；在景观水平上氮磷输出量与景观形状指数呈显著正相关, 和散布与并列指数表现出显著负相关, 与蔓延度指数、聚集度指数呈负相关。

因此, 在丹江口库区氮磷输出量高的地区, 应选择种植喜氮、喜磷的作物, 合理施肥灌溉, 可以有效降低氮磷进入水体的风险；可调节景观格局结构, 合理规划土地利用, 增加林地面积, 建设河岸缓冲带充分滞缓径流；制定排污标准, 平衡水源保护与经济发展关系, 确保南水北调中线工程水源水质安全。