生态系统服务是指生态系统为人类提供的福祉, 是人类生存与发展的基础^[1—2], 土地利用变化可以显著改变生态系统的模式和过程, 导致生态系统服务的变化^[3]。景观格局是由不同的土地利用类型在空间上呈现的分布与组合形式, 是各种生态过程在不同尺度上作用的结果^[4—5]。随着经济、人口的不断增长, 人类对于生态系统的利用和消耗持续增加, 带来了严峻的生态环境问题^[6], 导致生态系统服务不断退化^[1]。从景观角度考虑人与生态环境之间的相互作用, 能够为区域生态环境保护提供有力保障, 切实有效地促进区域可持续性发展^[7]。尤其是在人类活动与自然生态环境矛盾突出的快速城市化区域中, 明确景观格局与生态服务之间的关系, 可以为深入研究生态系统服务时空演变对景观格局的响应机制提供理论基础, 从而在景观格局优化角度对生态系统服务的提升提供新的方法途径。

景观格局与生态服务之间的耦合关系研究是生态系统服务管理以及景观可持续规划等跨学科领域的热点问题^[8]。目前, 大部分学者重点关注土地利用如何对生态系统服务产生影响, 分析自然生态系统服务价值动态变化与土地利用之间的关系, 例如：Terefe等分析了埃塞俄比亚森林地区的土地利用变化和生态系统服务价值变化, 发现林地面积占比的减小, 会使得该地区水量供应等重要生态系统服务功能丧失^[9]；张发等探究了博斯腾湖流域土地利用变化与生态系统服务价值时空演变, 发现水域、草甸和草原的面积变化是流域生态系统服务价值变化的主要影响因素^[10]。也有学者关注景观格局指数与生态系统服务价值之间的关系, 例如：Ayinuer等研究了艾比湖流域生态系统服务价值与景观格局变化之间的相关性, 发现景观多样性指数与生态系统服务价值有显著的正相关关系^[11], Hou等研究了西安市生态系统服务价值与景观格局变化之间的关系, 发现聚集度指数对生态系统服务价值有正向影响^[12]。其中, 生态系统服务价值是通过对不同土地利用类型进行赋值的价值当量法来计算, 赋值主要依据专家知识经验, 具有较大的不确定性^[13]。因此, 为了更加客观地分析景观格局与生态系统服务之间的定量关系, 本研究旨在：(1)从多时空尺度, 量化北京市景观组成、景观配置指数以及多种生态系统服务；(2)从景观和类型两个水平, 探究多种景观格局指数如何影响生态系统服务。

2016年公布的《新城市议程》中提到, 到本世纪中叶, 世界城市人口预计将增加一倍, 随着而来的城市化问题对于自然环境的可持续性带来巨大的挑战^[14]。尤其对于高度城市化地区, 北京市1980—2018年常住人口由904万人增长到2154万人, 地区经济总值由139亿元增长到33106亿元^[15]。在经济和人口不断发展的过程中, 北京市建设用地不断扩张, 景观格局发生剧烈变化^[16], 并且引发了一系列的生态环境问题^[17], 生态系统服务也随之受到威胁^[18], 研究北京市景观格局对生态系统服务的影响对于区域景观生态的可持续发展尤为重要。基于此, 本研究旨在：(1)定量分析北京市1980—2018年景观格局时空演变特征；(2)定量评估北京市1980—2018年生态系统服务时空演变规律；(3)探究北京市景观组成与景观配置指数对多种生态系统服务的影响。

1 研究区概况

北京市位于华北平原北部(39.4°—41.6°N、115.7°—117.4°E), 下辖16个区县, 面积为16410km², 是全国的政治、经济和文化中心。整体地势西北高、东南低, 其中, 西、北和东北三面环山, 山区面积约占总面积62%, 东南部相对平坦^[19]。平均海拔43.5m, 年均降水506mm^[20], 属于暖温带半湿润半干旱季风气候。2018年全市常住人口2154万人, 地区生产总值33106亿元, 城镇化率为86%^[15]。作为中国特大城市之一, 北京市1980—2018年经历了快速的城市化过程, 也伴随着巨大的景观格局变化和生态系统服务退化^[21]。

2 研究数据与方法 2.1 数据来源

本文的数据来源包括：(1)遥感影像数据：分辨率为30m的Landsat-MSS遥感影像(1980年)、Landsat-TM遥感影像(1990年、2000年、2010年)、Landsat-OLI遥感影像(2018年)5期影像数据, 来源于美国地质勘查局网站(https://www.usgu.gov/)。(2)土地利用数据：分辨率为30m北京市1980年、1990年、2000年、2010年、2018年土地利用遥感监测数据, 来源于中国科学资源环境科学与数据中心(https://www.resdc.cn/)。(3)统计数据：1980、1990、2000、2010、2018年北京市各地区主粮作物单产、面积数据, 来源于北京市统计年鉴^[15], 中国农业年鉴^[22]。(4)DEM数据：分辨率为30m的ASTER GDEM数据, 数据来源于美国国家航天航空局网站(https://www.nasa.gov/)。

2.2 景观格局分析 2.2.1 景观类型变化

基于土地利用数据, 将北京市分为耕地、林地、草地、水体、建设用地、未利用地6类, 利用经过镶嵌, 裁剪, 几何校正等处理后的遥感影像数据, 采用监督分类方法对土地利用数据进行校正, 得到Kappa系数大于80%且空间分辨率为30m的土地利用分布图, 分析北京市1980—2018年不同土地利用动态变化趋势。空间叠加1980、2018年北京市土地利用分布图, 并利用转移矩阵表示景观类型转移情况。

2.2.2 景观指数变化

景观指数是指能够高度浓缩景观格局信息, 有效反映其结构组成和空间配置等方面特征的简单定量指标, 分为3个层次：斑块、类型和景观水平^[23]。各景观指数本身有一定关联, 许多指数之间存在冗余, 所以选取典型指数来表征景观格局^[24]。本研究通过阅读相关文献^{[16, 23]}, 从类型和景观水平2个尺度共选取8个指数：斑块密度、平均斑块面积、边缘密度、聚集度、景观形状指数、Shannon多样性、连通性、蔓延度(表 1), 利用Fragstats 4.2.1软件计算选取的景观指数。

2.3 生态系统服务评估 2.3.1 粮食供给

以北京市1980—2018年农田产量数据来量化北京市粮食供给服务, 统计1980、1990、2000、2010、2018年北京市各地区主粮作物单产、面积数据, 通过SPAM(Spatial Production Allocation Model)模型对多源数据格式进行一致化处理, 基于模拟的作物空间概率分布以及作物空间竞争优化过程, 将行政单元统计的作物产量值分配至像元尺度^[25]。

2.3.2 水质净化

水质净化服务功能主要考虑植被对水体中氮元素的净化能力, 使用InVEST模型中的水质净化模块进行评估。本文选择单位面积生态系统总氮输出的多少来表示水质净化服务的大小, 即：单位面积生态系统总氮输出越少, 水质净化服务越高。该模块计算公式如下：

式中, WP_x是指栅格单元x的水质净化服务功能, ALV_x是栅格单元x调节的负荷值, pol_x是栅格单元x的输出系数, HSS_x是水文敏感分值, λ_x是栅格单元x的径流系数, λ_w是集水区的平均径流指数, 是集水区所有栅格产水量总和。模型的输入参数参考Gao、Zheng、李屹峰等的研究^[26—28]。

2.3.3 空气净化

通过各类生态系统(例如森林、草地等)对总悬浮颗粒物TSP、可吸入颗粒物PM₁₀、二氧化硫SO₂以及氮氧化物NO_x等颗粒截留吸收量来表征空气净化生态系统服务。

式中, S_i (TSP, PM₁₀, SO₂, NO_x)表示植被种类i对颗粒物的截留吸收量；A_i表示植被种类i的面积；P_i表示植物种类i单位面积对颗粒物的截留吸收量, 各景观类型的颗粒截留吸收量参数参考Wu、Landuyt等的研究^[29—30]。

2.3.4 气体调节

生态系统通过增减大气中的二氧化碳等温室气体调节气候, 本研究使用碳储量来表征气体调节生态系统服务。对每一种景观类型, 采用InVEST模型对四种基本碳库进行估算, 把各类碳库的估测值用于土地利用空间分布图上, 从而生成整个景观的碳储量分布图, 计算公式如下：

式中, S_c为总的碳储量；C_above为地上部分碳储量；C_below为地下部分碳储量；C_soil为土壤部分碳储量；C_dead为死亡有机质的碳储量。各景观类型的生物质碳储量参考Guo、He等的研究^[31—32]。

2.3.5 生境质量

栖息地的质量取决于来自威胁的外部压力强度, 以及栖息地类型对每种威胁的相对敏感性, 采用InVEST模型中的生境质量模块进行评估, 结果用生境质量指数来表示, 计算公式如下：

式中, Q_xj为土地利用类型为j时栅格x的生境质量, H_j表示土地利用类型j的生境类型评估分数, D_xj则表示土地覆被类型j上栅格单元x的生境胁迫水平, z值表示所用到的尺度常量值, 通常是2.5, k值是使用的半饱和常数。根据Wu等人的研究和当地调查确定本模块中的参数^[33]。

2.3.6 休闲娱乐

通过Recreation Opportunity Spectrum(ROS)模型进行评估生态系统中的水体吸引力来表征休闲娱乐服务^[34], 结果由娱乐机会指数呈现, 计算公式如下：

式中, d表示距河岸的距离, α和K分别表示函数的大小和形状参数。

2.4 相关性分析

采用偏相关系数表示1980—2018年北京市生态系统服务与景观格局之间的相关性, 系数的大小表示两者关系的强弱程度, 系数为正表示两者呈正相关, 系数为负表示两者呈负相关。本研究利用R 3.6.2实现生态系统服务与景观类型面积以及景观格局指数之间的偏相关分析和显著性检验, 并进行可视化。

3 结果与分析 3.1 景观格局变化分析 3.1.1 土地利用变化

根据北京市1980—2018年不同土地利用变化情况(图 1), 北京市景观类型以林地为主, 主要分布在北京市西部和北部的山区, 而东南部的平原地区以耕地与建设用地为主, 建设用地扩张地区主要位于北京市东南部；北京市1980—2018年耕地面积大幅减少, 建设用地面积大幅增加, 其中, 耕地面积比例由37%下降到23%, 共减少2155km², 建设用地面积由8%增加到18%, 面积增加1598km²。此外, 草地、水体以及未利用地的面积变化并不明显。

忽略面积比例较小的未利用地, 由北京市1980—2018年土地利用转移情况表(表 2), 可以看出：1980—2018年, 北京市土地利用转移面积达3803km², 占总面积的23%, 表现出耕地的大量转出和建设用地、林地大量转入的特征。其中, 耕地主要转为建设用地, 转出为建设用地的面积占耕地总转出面积的64%；建设用地主要由耕地转入, 占建设用地总转入面积的87%；林地主要由耕地转入, 占林地总转入面积的62%。其他土地利用类型之间有不同程度的转换, 但转换率相对较低。

3.1.2 景观指数变化

由北京市1980—2018年景观水平指数(表 3)发现, 北京市1980—2018年整体景观更加破碎化。斑块密度不断增加, 斑块数目增多；景观形状指数不断增加, 景观形状趋于复杂化；蔓延度指数先增后减, 景观类型的聚集程度降低；连通度指数不断减小, 各景观类型之间连通性降低, 斑块空间分布趋于零散；Shannon多样性指数先增后减, 总体减少, 各景观类型面积分布越来越不均匀。

由北京市1980—2018年类型水平景观指数(表 4)发现, 耕地在1980—2000年斑块数目增多, 平均斑块面积减小, 斑块边缘密度增大, 边界形状复杂化, 斑块聚集度降低, 整体呈破碎化发展。2000—2018年, 耕地呈现不同的变化规律, 耕地斑块聚集性增强；林地、草地以及建设用地在1980—2018年斑块密度增加, 平均斑块面积减小, 边缘密度增加, 斑块聚集度降低, 呈现空间破碎化, 其中建设用地斑块密度和边缘密度均大幅增加, 平均斑块面积大幅减少, 聚集度降低, 斑块破碎化程度高。

3.2 生态系统服务变化 3.2.1 北京市1980—2018年生态系统服务时间变化

北京市1980—2018年生态系统服务变化(图 2)显示, 除空气净化外的其他生态系统服务均呈现下降趋势。1980—2018年北京市粮食生产、水质净化、气体调节和休闲娱乐服务先增后减, 总体呈现下降趋势, 生境质量服务持续减少。其中, 水质净化服务在1990—2000年大幅减少, 空气净化、气体调节、休闲娱乐服务在2010—2018年显著降低。

3.2.2 北京市1980—2018年生态系统服务空间变化

根据北京市1980—2018年生态系统服务空间变化(图 3)显示, 生态系统服务呈现显著的空间差异性。除粮食生产服务在东部、南部较高, 西部、北部较低, 其他服务均呈现西部、北部较高, 东部、南部较低, 中部最低的空间分布特征；1980—2018年北京市中部地区粮食生产服务明显减少, 空气净化、气体调节、生境质量、休闲娱乐服务在空间上呈现东部和南部地区减少的趋势, 尤其是在建设用地扩张地区呈现显著减少的变化特征。

3.3 北京市生态系统服务对景观格局的响应 3.3.1 生态系统服务对景观水平的响应分析

由生态系统服务与土地利用类型面积的偏相关分析(图 4)可知, 景观组成对生态系统服务有显著影响。耕地面积与粮食生产以及水质净化服务有正相关关系；林地、草地面积与空气净化、气体调节、生境质量、休闲娱乐服务均有正相关关系, 其中空气净化、气体调节、休闲娱乐等服务对林地的响应最为强烈；建设用地面积与所有的生态系统服务均有负相关关系, 其中生境质量服务对建设用地的响应最为强烈。

由生态系统服务与景观水平指数的偏相关分析(图 4)可知, 区域内景观破碎化程度以及多样性对生态系统服务有显著影响。斑块密度与区域内所有生态系统服务均有负相关关系, 斑块形状指数与区域内所有生态系统服务均有正相关关系；景观多样性与除粮食生产服务外的其他生态系统服务有正相关关系。其中, 斑块密度与生境质量服务有显著的负相关关系, 景观多样性与生境质量服务有显著的正相关关系。

3.3.2 生态系统服务对类型水平的响应分析

由生态系统服务与景观类型指数的偏相关分析(图 5)可知, 各景观类型的景观指数对生态系统服务的影响具有差异性。耕地的平均斑块面积与边缘密度与粮食生产服务有显著的正相关关系；林地的平均斑块面积和空间聚集度与空气净化、气体调节、生境质量、休闲娱乐服务有显著的正相关关系；建设用地边界密度和空间聚集度与空气净化、气体调节、生境质量、休闲娱乐服务有显著的负相关关系。

4 讨论

1980—2018年, 北京市中心城区建设面积不断向外扩张, 斑块呈现破碎化趋势, 其中受影响最大的是耕地, 耕地斑块变得破碎且分散^[35—36]。同时, 城市化发展过程对北京市的景观结构配置造成显著影响, 尤其是2000年之后, 受到人类活动的剧烈影响, 北京中心城区、东南部地区建设用地组分优势度明显, 导致北京市整体景观破碎化增加, 多样性降低^[16]。

北京市在过去38年间, 生态系统服务严重下降, 这主要是由于人类频繁的社会经济活动造成了自然生态系统的严重破坏^[37]。尤其是东南部地区, 建设用地占比不断增大, 因此生态系统服务明显下降。西部和北部地区, 山地面积占比较大, 受人类活动干扰较少, 植被生物量较高^[38], 提供的调节和文化生态系统服务相对更高。

土地利用变化是北京市生态系统服务变化的主要原因, 当林地等自然生态系统向建设用地等人工生态系统发生转变时, 会带来生态系统服务与生物多样性的减少^{[18, 39]}。林地面积的增加有助于提高生物物质以及木材生产的碳储存量, 并适度促进土壤中碳的储存^[40], 植被叶片的微观结构与去除空气中大小颗粒能力密切相关^[41], 因此, 增加林地面积对空气净化、气体调节服务有显著的正向效应；随着经济和人口的增长, 建设用地不断扩张, 自然生境遭到破坏, 导致生物多样性降低, 区域生境质量变差^[42]。

不同的景观配置会对北京市生态系统服务有显著影响。增加景观多样化与景观类型均衡分布, 有助于提高生态系统的稳定性^[11], 而破碎的景观会阻碍区域间的能量流动, 改变物质和营养循环过程, 对生态系统服务有显著的负面影响^[43]。耕地斑块破碎化不利于农田的规模化管理和农业资产的合理配置, 会严重影响粮食生产效率^[44—45]；林地斑块面积越大, 空间分布越聚集, 生态系统提供的空气净化、气体调节、生境质量服务则越多^[46]；不规则建设用地斑块数量越多, 空间分布越聚集, 造成自然生境斑块更加破碎, 连通性不断降低, 生境质量不断下降^[47]。

本文从景观组成与景观配置两个方面, 分析了北京市景观格局对多种生态系统服务的影响；未来的研究应基于具体的生态过程及其机制, 通过改变景观格局, 引起生态过程变化, 从而提升区域生态系统服务；针对特定区域进行景观优化, 在保证生态系统稳定性的同时为人类提供更多的福祉, 促进生态系统的可持续发展^[48—49]。

5 结论

(1) 1980—2018年, 北京市景观格局发生剧烈变化, 建设用地面积显著增加, 耕地面积显著减少；景观整体呈现破碎化趋势, 斑块异质性增加。多数生态系统服务呈下降趋势, 其中东南部地区的生态系统服务减少最为明显。

(2) 景观水平上, 增加林地、草地面积对空气净化、气体调节服务有显著的正向效应, 增加建设用地面积对生境质量服务有显著的负向效应；降低景观破碎化程度, 增加景观多样性对生态系统服务有正向效应；类型水平上, 增加耕地的平均斑块面积对粮食生产服务有显著的正向效应, 增加林地的聚集度对空气净化、气体调节服务有显著的正向效应, 而增加建设用地的聚集度对生境质量服务有显著的负向效应。

(3) 未来城市化发展过程中, 针对生态系统服务的空间异质性特征, 北京市不同区域应采取因地制宜的景观可持续规划策略。在顺义区、通州区等东部地区, 应重点提升粮食生产服务, 避免建设用地无序扩张与破碎化发展, 优先保护耕地资源, 增加平均斑块面积与空间聚集度；在市中心地区, 应重点改善调节服务, 减小建设用地的聚集度, 提高林地、草地的等绿色基础设施的景观组成比例与空间聚集度；在怀柔区、密云区、门头沟区等西北部地区, 应充分发挥其较高的调节和文化服务供给能力, 避免林地、水体等自然景观面积比例的减少, 从而促进区域的景观可持续发展。