生态系统服务是指人类可以直接或者间接从生态系统中获得的各种效益^[1], 通常被划分为供给、调节、支持和文化四类服务^[2]。目前, 由于人类活动的干扰和高强度的经济开发活动, 全球生态系统服务一直处于不断退化中。例如, 2005年的千年生态系统服务评估报告指出, 全球范围内有一半以上的生态系统服务出现明显退化^[3]。生态系统服务保护与管理已成为当前人们关注的热点, 但是, 由于生态系统空间分布结构的异质性和生态过程的复杂性, 各类生态系统服务间都或多或少地存在着此消彼长的权衡关系和共同增减的协同关系^[4-5], 对某种生态系统服务的保护和管理有可能导致另外一种生态系统服务受到损害^[6]。因此, 正确认识生态系统服务间的权衡与协同关系及其在不同空间尺度上的差异, 是规划生态功能区、完善生态保护建议的关键举措, 也是推进区域协调发展的必要步骤。

随着自然资源开发利用与生态保护间的矛盾日益凸显, 探究生态系统服务间的权衡与协同关系已成为生态系统服务研究的核心议题之一。目前, 国内外学者主要应用统计学方法^[7-8]、空间分析方法^[9-10]和情景模拟方法^[11-12]等开展研究, 内容多聚焦于权衡与协同关系的概念内涵、统计数量、时空特征和区域差异等^[13-14], 研究结果适用于农业生产^[15]、渔业管理^[16]、水源区管理^[17]等多方面。例如, 饶胜等构建草地生态系统服务权衡模型, 分析供给功能和防风固沙功能之间的权衡作用^[18]。Turner等在10 km×10 km的网格尺度上对丹麦的11种生态系统服务进行空间分析, 结果显示供应服务与文化服务、调节服务之间存在着明显的权衡倾向, 文化服务与调节服务具有形成协同效应的潜力^[19]。Qiao等在1km²像素尺度和县域尺度上对太湖流域8种生态系统服务间的权衡与协同关系进行研究, 结果显示在两种尺度上生态系统服务间的关系有所不同^[20]。傅伯杰和于丹丹通过对生态系统服务权衡与集成方法的深入研究, 指出生态系统服务间的消长和权衡具有尺度依赖性^[21]。诸多学者在生态系统服务权衡与协同方面的研究已经取得一些成果, 但是这些研究大多基于统计关系的数量分析, 在权衡与协同关系的动态变化分析和区域内部空间差异分析等方面的研究相对较少^[22]。

南四湖是重要的鸟禽以及水生动植物栖息地, 是山东省自然保护区, 在山东省具有“生物宝库”的美称。同时, 南四湖在南水北调东线输水干线上具有重要的调蓄作用^[23]。因此, 无论是在生态环境保护还是经济发展方面, 南四湖流域都具有不可忽视的重要地位。伴随着城市化水平飞速提高, 南四湖流域内工农业活动不断发展, 人类对生态系统的干扰程度不断增强, 生态环境的保护与发展成为流域管理活动中的一大挑战^[24-25]。因此, 本文以南四湖流域为研究对象, 基于1980、2000、2018年3期土地利用数据, 结合气象、土壤等多种数据, 利用InVEST模型模拟测算流域关键生态系统服务, 包括产水量、水质净化、土壤保持、碳存储、生物多样性维持(生境质量)和授粉等6种生态系统服务, 在ArcGIS 10.2中分析生态系统服务的时空变化和空间冷热点分布特征, 使用相关分析展现生态系统服务间的相互作用关系, 并运用主成分分析从子流域尺度和县域尺度分析权衡与协同关系的动态变化, 以期揭示南四湖流域生态系统服务的时空差异和生态系统服务权衡与协同关系在时间上的变化特征, 以及在不同的尺度上权衡与协同关系的区别。

1 研究地区与数据准备 1.1 研究区概况

南四湖由微山湖、昭阳湖、独山湖、南阳湖等相连湖组成, 是山东省的第一大湖, 也是中国大型淡水湖泊之一。流域地处山东省西南部, 位于34°15′—34°53′N, 114°50′—117°48′E, 总面积约3.14万km², 行政上涉及山东省、江苏省、安徽省和河南省。南四湖流域地势东高西低, 湖区是流域内地势最低的区域, 流域内的所有河流呈辐射状汇入湖区(图 1)。南四湖流域具有丰富的矿产资源和水产资源, 盛产鱼、虾、蟹等水产品, 是山东省重要的淡水渔业基地。流域内以种植小麦、水稻、玉米等粮食作物为主, 是鲁西南的鱼米之乡, 同时也是全国重要的商品粮基地。流域内林地、草地和湖区等地区生物多样性较高, 在山东省具有较高的生态地位。近年来, 作为典型的农业流域, 南四湖水土流失加剧, 非点源污染严重, 水环境质量恶化, 且动植物栖息地丧失, 生境质量下降, 生物多样降低, 流域各类生态系统服务功能退化。根据南四湖流域生态环境特征, 本研究选择流域关键生态系统服务展开研究, 包括产水量、水质净化、土壤保持、碳存储、生物多样性维持(生境质量)和授粉等6种区域关键生态系统服务, 探讨南四湖流域各类生态系统服务及其权衡与协同关系, 以期为流域资源配置以及生态环境保护政策制定提供重要参考。

1.2 数据来源

本研究使用的南四湖流域1980、2000、2018年3期土地利用数据, 来自中国科学院资源环境科学与数据中心(http://www.resdc.cn/), 通过遥感影像解译获得, 分辨率为30 m。降水和气温数据来自中国气象数据网(http://data.cma.cn/), 选取南四湖流域及其周边的气象站点数据进行空间插值。土壤数据来自联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)构建的世界土壤数据库(HWSD)中的中国土壤数据集(1 ∶ 100万)。DEM数据来自美国马里兰大学全球土地覆被数据库(http://glcfapp.glcf.umd.edu/), 分辨率为30 m。县级行政区划矢量数据来自全国基础地理数据(1∶400万)。

2 研究方法 2.1 生态系统服务价值评估方法

本文使用InVEST模型和ArcGIS 10.2软件, 定量评估南四湖流域产水量、水质净化、土壤保持、碳存储、生物多样性维持和授粉六种生态系统服务。

2.1.1 产水量

InVEST模型的产水量模块是一种基于水量平衡并且结合气候、地形、土壤和植被类型等计算得出流域每个栅格产水量的估算方法^[26]。产水量越多, 水资源供给量越大。模型参数设置和校验参考课题组前期研究^[23], 产水模块计算的基本原理如下：

(1)

式中, Y_xj为土地利用类型j上栅格x的年产水量(mm)；AET_xj为土地利用类型j上栅格x的年实际蒸散量(mm)；P_x为栅格x上的年降水量(mm)。

2.1.2 水质净化

InVEST模型的水质净化模块是利用氮和磷的总量来表征水质状况, 通过评估生态系统中植被和土壤对氮磷营养物质的截留能力来评估其水质净化的能力^[27]。模型评估分为两个部分, 首先通过产水量模型计算年平均径流量；然后计算污染物截留量。模型参数设置参考模型帮助文件及相关研究文献[27-29], 水质净化模块计算的基本原理如下：

(2)

(3)

式中, ALV_x是栅格x的营养物输出值, HSS_x是栅格x的水文敏感性得分, pol_x是栅格x的输出系数, λ_x 是径流指数, λ_w表示区域内平均径流指数。

氮磷营养物输出越高, 水质净化能力越低, 为合理表征流域水质净化功能, 对模型模拟结果负向标准化后进行下一步分析。

2.1.3 土壤保持

InVEST模型的土壤保持模块根据地形地貌和气候条件评估潜在土壤流失量, 通过通用土壤流失方程(USLE)评估土壤保持能力。土壤保持量包含侵蚀减少量和泥沙持留量^[30-32], 通用土壤流失方程整合了土地利用、土壤性质、DEM、降水和气温等信息。模型参数设置参考模型帮助文件和前期研究^{[27, 30-31]}, 土壤保持模块计算的基本原理如下：

首先, 计算潜在土壤侵蚀量, 潜在土壤侵蚀量计算公式如下：

(4)

其次, 根据通用土壤流失方程计算得出生态系统中的土壤实际侵蚀量, 计算公式如下：

(5)

最后, 土壤保持量计算公式如下：

(6)

式中, R表示降雨侵蚀力因子, K表示土壤可蚀性因子, LS表示坡度坡长因子, C表示植被覆盖和管理因子, P表示水土保持措施因子。

土壤侵蚀越强, 泥沙输出越高, 土壤保持服务越低, 为合理表征流域土壤保持服务, 对模型模拟结果负向标准化后进行下一步分析。

2.1.4 碳存储

InVEST模型的碳存储模块根据不同的土地利用类型数据及其所对应的五大碳库的碳密度等数据计算不同时期和不同土地类型的碳储量^[33]。由于中国市场的木材衰减率等信息较难获取, 只考虑前四种基本碳库。模型参数设置参考模型帮助文件和前期研究^{[27, 33]}, 碳存储模块计算的基本原理如下：

(7)

式中, C_z为总碳储量, C_above为地上碳储量, C_below为地下碳储量, C_dead为死亡有机质碳储量, C_soil为土壤碳储量。各碳储量由碳密度与面积相乘所得。

2.1.5 生物多样性维持(生境质量)

InVEST模型的生境模块以土地利用信息和生物多样性胁迫因子为基础开展生境质量评价, 进而评估生物多样性维持功能^[34]。模型依据各生态系统类型对动植物的生境适宜度和人类干扰因子的威胁强度来模拟评估生境质量的空间分布。模型参数设置参考模型帮助文件和前期研究^{[27, 35]}, 生境模块计算的基本原理如下：

(8)

式中, Q_xj为土地利用类型j中栅格x的生境质量指数; H_j为土地利用类型j中的生境适宜度; D_xy^z为土地利用类型j中栅格x的生境退化度; z为模型默认参数；k为半饱和常数。

2.1.6 授粉

InVEST模型的授粉模块利用筑巢地和花蜜资源的可用性以及野生蜜蜂飞行距离的范围来获得在景观内每个栅格单元上蜜蜂筑巢的丰度指标, 再结合飞行距离信息来估算农业栅格单元内蜜蜂访问花卉的物种丰度指标^[36]。模型参数设置参考模型帮助文件和前期研究^{[27, 36]}, 授粉模块计算的基本原理如下：

由于传粉者丰度受到筑巢和花蜜资源的限制, 对于栅格x和物种β的物种丰度的指数, P_xβ是传粉者丰度觅食和筑巢的产物：

(9)

式中, N_j是土地利用类型j筑巢的适用性, F_j是土地利用类型j产生的花蜜资源的相对数量, D_mx是栅格m和x之间的欧氏距离, α_β是传粉者β的预期觅食距离。

根据上述方程的框架确定每个农作物栅格上的蜜蜂的相对丰度即为：

(10)

式中, P_xβ是传粉物种β在栅格x传粉者供应, D_oxβ是物种β从源栅格x到农业栅格o的距离, α_β是物种β的平均觅食距离。

2.2 生态系统服务权衡与协同关系评估方法

为了分析南四湖流域各类生态系统服务之间的相关性, 揭示其权衡与协同关系, 本研究参考太湖流域^[20]的研究方法, 在ArcGIS 10.2中使用“Create Random Points”工具随机创建1000个点, 然后提取每个点的生态系统服务值进行相关分析。如果生态系统服务间的相关性数值为负, 即为权衡关系, 反之则为协同关系^[37]。

从行政区划因素和自然地理性质两方面出发, 选取县域尺度和子流域尺度来探讨流域各类生态系统服务间的权衡与协同关系以及尺度效应。采用主成分分析的方法, 在县域尺度和子流域尺度进行定量分析, 以确定各类生态系统服务以及三大类生态系统服务之间的权衡与协同关系。

3 结果与分析 3.1 生态系统服务时空变化

南四湖流域6种生态系统服务的时空分布和不同土地利用类型的生态系统服务时空差异分别由图 2、图 3所示。结果显示产水量总体上处于增加的趋势, 1980—2018年南四湖流域的平均产水深度从101.83 mm增加到185.53 mm；从空间分布来看, 2000年开始流域南部城市化较快的县区产水量明显增加, 到2018年城市建设用地所对应的产水量高值区相较于其他地区更加显著, 洪水风险增加。在水质净化服务中, 1980—2018年氮、磷输出平均值分别由0.46 kg和0.02 kg增加到0.49 kg和0.03 kg, 水质净化服务降低；在空间上看, 2000年后流域南部氮、磷输出明显增多, 该区域根底比重大, 氮磷负荷较高, 故水质净化服务降低；2018年建设用地分布与氮、磷输出高值区一致, 建设用地氮、磷输出较多。泥沙输出量总体呈增加趋势, 近40年来单位面积泥沙输出量由0.01 t增加到0.05 t, 土壤保持功能逐渐下降；从空间上看, 建设用地分布的区域泥沙输出量明显较低, 原因在于城市扩张过程中建设用地面积扩大, 不透水地面的增加减少了泥沙输出。研究期内碳储量先增加后下降, 总体呈下降的趋势；从空间分布上看, 东部林地和草地碳存储能力最强, 耕地次之, 建设用地和水域依次降低。生境质量总体上处于不断降低的趋势, 1980—2018年从0.21降低到0.20, 说明整体上生物多样性维持服务较低并且处于逐渐恶化的趋势。相对授粉潜力在1980—2018年总体上处于下降趋势, 这与研究期内土地利用变化密切相关, 建设用地和水域的增加, 耕地减少, 导致相对授粉潜力下降；空间分布上, 东部林地和草地相对授粉潜力较高, 耕地次之, 建设用地和水域最低。

本研究制作风向玫瑰图表征不同土地利用类型上的各类生态系统服务功能。首先, 在ArcGIS 10.2中使用空间统计工具统计每种土地利用类型上6种生态系统服务的平均值, 由于6种生态系统服务值的单位和数量级不同, 在分析过程中具有较大困难, 因此, 选择极差标准化对生态系统服务值进行标准化处理。本文将研究重点放在不同生态系统服务之间的关系上, 以土地利用类型为基础, 将3个年份同种土地利用类型上的六种生态系统服务价值进行升序排列, 在序列基础上进行等级划分, 从0到1划分为20级。最后, 为方便结果的可分析性及可视性, 运用origin对各土地利用类型的生态系统服务数据制作南四湖流域不同土地利用类型生态系统服务风向玫瑰图(图 3)。在草地中, 授粉服务与土壤保持服务较高, 水质净化功能次之。在耕地中, 水质净化功能一直较高, 但是在研究期内耕地的泥沙输出量越来越大, 土壤保持功能逐渐下降。建设用地土壤保持与产水服务相对较高。在林地中, 一直存在较高的碳存储功能。当土地利用为水域时, 氮、磷输出和泥沙输出数值都低到接近于0, 因此水质净化和土壤保持功能呈现较高的状态, 生境质量也一直较高。未利用地中, 土壤保持功能一直较高, 产水功能在不断增强。

3.2 生态系统服务空间冷热点分析

对各种生态系统服务进行空间冷热点分析(图 4), 结果表明不同种生态系统服务的冷热点区域出现类似分布, 比如碳存储、生境质量和授粉服务, 研究期内东部丘陵山区一直为这3种服务的热点区域, 与东部的土地利用类型大多为林地和草地有密切关系, 湖区南部和北部小部分区域则为这些服务的冷点区域, 该区域多为耕地, 人类活动干扰强烈, 生态系统服务价值低。

1980年, 产水量服务的热点区域分布在研究区东部, 而2000年和2018年热点区域分布在研究区南部, 这与南部各县区城市化发展较快, 地面硬化导致的不透水地面增加有关, 研究区北部一直是产水量的冷点区域, 因为当地海拔较高, 山地上的各类植物对降水有一定的拦蓄作用。水质净化服务中, 1980年很少有热点区域, 水域的周围出现了冷点区域, 说明氮和磷输出较少。2000年和2018年流域南部出现热点区域, 与同时期城市扩张, 氮磷输出增多有关, 水质净化功能较差。冷点区域分布在研究区东北部, 原因是这些地区草地和林地的面积较大, 拦截氮磷输出, 水质净化功能较好。土壤保持服务中, 流域东北部一直为泥沙输出的热点区域, 因为当地海拔较高, 土壤流失严重, 泥沙输出较多, 土壤保持功能较差。土壤保持服务的冷点区域变化较大, 1980年冷点区域分布在研究区南部, 2000年和2018年流域北部兖州市、曲阜市出现少量的冷点区域, 与当地经济发展, 建设用地增多, 不透水地面土壤流失较少有关。

3.3 生态系统服务权衡与协同关系分析

通过相关分析, 研究6种生态系统服务间的关系及其强度(表 1)。1980—2018年, 产水量与水质净化、碳存储和生境质量均存在显著的负相关, 属于权衡关系, 其中产水量与水质净化的权衡关系更为突出, 相关系数绝对值均在0.5左右, 且随着时间的推移相关性不断增强, 到2018年产水量与水质净化相关系数高达-0.597。

水质净化与碳存储和生境质量之间存在显著的正相关, 其中水质净化与生境质量之间的相关系数最高值达到0.636, 具有较强的协同关系。碳存储和生境质量、碳存储和授粉之间存在显著的正相关, 且碳存储和授粉相关系数在2018年高达0.849, 说明碳存储和授粉服务是相互增益的, 山地丘陵区的花蜜可得性高且人类活动干扰较少, 并且在山地丘陵区林地和草地面积较大, 因此碳储量和授粉服务都较高。生境质量与授粉服务在1980年和2000年存在显著正相关, 但相关系数均较低。

综上所述, 研究区的供给服务(产水量)与部分调节服务(水质净化和碳存储)和部分支持服务(生境质量)之间存在此消彼长的权衡关系, 调节服务和支持服务之间(水质净化、土壤保持、碳存储、生境质量和授粉)都是相互增益的协同关系。因此, 当前生态管理策略制定的核心就是要综合平衡好流域供给服务与调节服务以及供给服务与支持服务间的权衡关系, 优化流域生态系统服务结构。

3.4 生态系统服务主成分分析

使用主成分分析来确定各类生态系统服务间以及各大类生态系统服务间的权衡和协同关系(表 2和表 3)。1980年, 在子流域尺度上前两个主成分占总贡献率的85.948%(特征根大于1), 第一主成分占贡献率的62.767%, 代表水质净化与碳存储、授粉、产水量之间的权衡关系, 第二主成分占总贡献率的23.181%, 代表了水质净化(N输出)和生境质量之间的协同作用。2000年, 第一主成分占总贡献率的56.965%, 代表水质净化、土壤保持、生境质量与产水量、碳存储之间的权衡作用, 第二主成分代表了碳存储与授粉服务之间的协同关系, 占总贡献率额外的25.5%。2018年, 第一主成分占总贡献率的58.767%, 代表了水质净化、生境质量与碳存储、授粉、产水量之间的权衡关系, 第二主成分占总贡献率的24.483%, 代表了存在于土壤保持与授粉服务之间的空间权衡作用。

由上可知, 1980—2018年, 各生态系统服务间权衡与协同效应不断变化, 但总的来说, 在子流域尺度上, 空间权衡关系总是发生在供给服务与调节服务以及供给服务与支持服务之间。

在县域尺度上, 1980—2018年, 第一主成分的贡献率都在40%以上, 且都代表了水质净化、生境质量和碳存储之间的协同效应。1980年产水量与土壤保持具有协同作用, 2000年, 新的协同作用出现在产水量、碳存储和授粉服务之间。1980年与2018年水质净化(P输出)和授粉服务之间都具有权衡关系。1980年和2018年, 供给服务、调节服务和支持服务之间都只存在相互增益的协同关系, 然而2000年, 并非只有协同关系, 在供给服务与调节服务、供给服务与支持服务之间存在权衡关系。各生态系统服务间的权衡与协同关系在县域尺度上与子流域尺度上具有较大差异, 但是2000年县域尺度上的权衡关系与子流域尺度上3大类生态系统服务之间的权衡关系一致。

4 讨论与结论 4.1 讨论

不同尺度生态系统服务的时空权衡与协同关系研究旨在充分认识和掌握各类生态系统服务之间的相互作用关系、驱动机制和尺度效应等, 寻求区域经济开发、自然资源利用与生态环境保护之间的平衡点, 从而实现经济发展与生态保护的双赢。本研究直接揭示了研究区关键生态系统服务之间在数量与空间上的权衡与协同关系, 并探讨其尺度效应, 但是各生态系统服务之间的相互作用机理方面还需要深入研究。

在本研究区, 产水量和与水质净化服务呈显著负相关, 与丹麦Raudsepp-Hearne等^[38]的研究结果相同, 其中产水量与水质净化的权衡作用最明显且随时间不断增强, 随着经济发展建设用地面积不断增加, 产水功能越来越强, 洪水风险在不断增加。而由于各个工厂企业的氮和磷输出较高, 城市建设用地所占区域的水质净化功能降低。水质净化与生境质量、碳存储与授粉之间都存在较强的协同关系, 在研究区林地和草地广布的地区, 不但具有较高的生物多样性, 还具有较高的花蜜可用性, 同时植物可以起到截留污染物的能力。对于土壤保持和碳存储的权衡与协同关系, 在不同研究区有不同的结论, 如戴路炜^[9]在多伦县的研究表明两者之间存在协同关系, 而Egoh等^[39]在南非的研究表明这两者之间存在权衡关系, 本研究中土壤保持和碳存储两者之间关系并不显著, 可能因为不同地区自然环境的差异, 导致研究结果不同。

各类生态系统服务间的权衡与协同关系在不同尺度上经常会有不同的体现, 孙泽祥等^[40]在呼包鄂榆地区的研究结果显示, 在城市群尺度上, 产水量和固碳服务表现出显著的权衡关系, 而在区域尺度上这两种服务间的权衡关系只体现在农业区。本研究中, 子流域尺度和县域尺度上, 生态系统服务间的关系也存在较大差异。子流域尺度上表现出的空间权衡关系与加拿大^[38]、太湖^[20]等研究结果一致, 权衡关系主要存在于供给服务与调节服务、供给服务与支持服务之间, 而本研究在县域尺度上主要表现出协同关系, 在1980年和2018年, 空间协同关系存在于供给服务、调节服务和支持服务之间。两个尺度上表现出差异的原因是不同尺度的每一个小单元上土地利用类型的构成和配置会不同, 在县域尺度上行政划分导致土地利用构成较为综合复杂, 而子流域尺度各单元则具有相似的自然地理特征, 可以从单元内部的土地利用构成角度更加详细地揭示各类生态系统服务间的作用机制。

4.2 结论

本研究根据南四湖流域自然地理特征、生态环境状况和社会经济特征, 选择了产水量、水质净化、土壤保持、碳储存、生物多样性维持(生境质量)、授粉等6类区域关键生态系统服务开展研究, 运用InVEST模型模拟评估了南四湖流域关键生态系统服务的空间格局和时间变化趋势, 在此基础上, 运用相关分析和主成分分析揭示了自改革开放以来(1980—2018年)这些生态系统服务的时空权衡关系。

(1) 模拟结果表明, 1980—2018年南四湖流域产水量上升, 洪水风向增加, 而水质净化服务、土壤保持服务、碳储存、生境质量、授粉服务总体呈下降趋势。

(2) 各类生态系统服务的高值区分布与流域土地利用分布格局有密切关系, 研究区东部林地和草地广布, 生物多样性较高, 人类活动干扰较小, 是碳存储、生境质量和授粉服务的高值区。城市建设用地分布的地区水泥硬化路面较多, 降雨的下渗量较小, 是产水量和土壤保持服务的高值区。

(3) 研究发现, 产水量和与除土壤保持外的调节服务显著负相关, 水质净化与生境质量、碳存储与授粉之间都存在较强的协同关系。而且随着时间的变化, 各类生态系统服务之间的相关性类型和强度在不断变化。

(4) 在子流域尺度和县域尺度上, 各类生态系统服务间的权衡与协同关系存在较大差异。在子流域尺度上, 权衡关系主要存在于供给服务(产水量)与调节服务(水质净化和土壤保持)、供给服务(产水量)与支持服务(生境质量)之间；在县域尺度上, 供给服务、调节服务与支持服务(产水量、水质净化、碳存储、生境质量和授粉)之间主要是空间协同关系。子流域和县域尺度上权衡与协同关系的差异, 为流域不同尺度的生态保护与管理提供了数据和理论支撑。