生态系统服务是人类从生态系统中所获得的惠益^[1], 包括生态系统向人类社会提供的人类生存所需的物质资源与服务, 与人类福祉密切相关^[2]。为有效地对生态系统进行管理并制定合理的政策, 生态系统服务价值核算成为一种重要手段^[3][4], 其可以识别区域重要的生态保护地^[5], 帮助理解生态系统服务之间的权衡和协同作用^[6][7]。随着GIS、RS等技术的发展及生态系统服务价值评估相关理论的不断完善, 生态系统服务价值的动态变化研究已成为当前研究趋势之一^[8-9]。但现有研究大多仅关注某个单一类型生态系统服务的价值评估^[10], 或对多种生态系统服务的总价值进行不同时间的对比^[11][12], 缺少对生态系统服务价值组成结构的探讨, 以及长时间序列的动态演变过程分析。

生态系统服务簇为此提供了一个可行的定量视角。生态系统服务簇是指一系列在空间和时间上重复出现的多种类型生态系统服务的组合^[13], 在研究各种生态系统服务之间的权衡与协同关系^[14]、区域主导生态系统服务的识别^[15][16]、生态功能区划^[17]、景观规划与管理^[18]等方面得到应用。目前生态系统服务簇的识别方法已较为成熟, 大多采用聚类的方法对多重生态系统服务进行分析, 得到的“类”即为研究区生态系统服务簇^[17][19]。借助生态系统服务簇的概念可以对区域不同生态系统服务的组合规律进行研究, 明确生态系统服务的价值构成, 识别出区域主导生态系统服务类型^[16], 有助于了解不同类型生态系统服务分布的关联特征^[17]。然而, 目前关于生态系统服务簇的研究多为静态分析, 对于生态系统服务簇的演变过程关注较少, 缺少长时间序列的研究。

本文提出一种基于生态系统服务簇的生态系统服务时空演变轨迹识别方法。这一方法从生态系统服务簇的视角对区域不同生态系统服务的组合规律进行研究, 使用K-means算法得到区域生态系统服务簇类型, 在此基础上对区域的生态系统服务簇进行长时间序列分析, 得到生态系统服务的时空演变轨迹, 并总结出生态系统服务时空演变的主导轨迹类型, 以期明确城市生态系统服务的演变特征, 为提出针对性的区域开发策略提供决策依据。

1 数据与方法 1.1 研究区概况

深圳市位于广东省中南沿海地区、珠江入海口东偏北, 所处纬度较低(113°43′—114°38′E, 22°24′—22°52′N), 属于亚热带季风气候区。截至2017年末, 全市下辖福田、罗湖、盐田、南山、宝安、龙岗、龙华、坪山、光明、大鹏10个区(图 1), 57个街道(2018年行政区划调整后, 深圳市辖74个街道), 土地总面积1997.47 km²；常住人口1252.83万人, 其中户籍人口434.72万人, 实际管理人口超过2000万, 城市化率100%。改革开放以来, 深圳在经济建设方面取得了辉煌的成就, 深圳的GDP从1979年的1.97亿元上升到2017年的2.25万亿元。但随着城市化的发展, 大面积的生态用地被用于城市建设, 原有生态过程随着土地利用类型的改变而消失, 城市生态系统的结构和功能均发生了较大变化。

1.2 数据来源

本文主要利用深圳市各街道的行政区划数据和1980、1990、2000、2010、2015年5个时期的土地利用遥感监测数据。行政区划数据来自深圳市规划和国土资源委员会2014年的街道级行政区划矢量数据。土地利用遥感监测数据使用中科院资源环境科学数据中心(http://www.resdc.cn/data.aspx?dataid=251)解译的中国土地利用现状遥感监测数据, 土地利用类型包括耕地、林地、草地、水域、居民地和未利用土地6个一级类型以及25个二级类型。之后根据研究区域自然条件及实际土地利用情况, 结合相关研究^[20]进行土地利用重分类, 最终将深圳市的土地利用生态系统分为水田、旱地、阔叶林、灌木林、灌草丛、水系、建设用地、裸地、湿地9类。

1.3 生态系统服务聚类

目前对生态系统服务的分类普遍采用二级分类体系^[20], 其中一级生态系统服务类型为供给服务、调节服务、支持服务、文化服务4个大类^[21]。具体而言, 供给服务是指生态系统提供食物和水等物质的服务, 包括食物生产、原料生产、水资源供给3个二级服务类型；调节服务是指调节生态环境的服务, 包括气体调节、气候调节、净化环境、水文调节4个二级服务类型；支持服务为其他生态系统服务提供前提条件, 包括土壤保持、维持养分循环、生物多样性3个二级服务类型；文化服务是人们从生态系统中所获得的非物质利益, 本研究重点考察美学景观。由此, 确定了生态系统服务的4个一级类型和11个二级类型。

本研究采用谢高地等^[20]于2015年提出的“中国二级生态系统服务价值当量因子表”对生态系统服务价值进行定量评估, 该表在以实物量计算方法为主的生态系统服务价值量评价结果的基础上, 结合了遥感影像数据对净初级生产力(Net Primary Productivity, NPP)和生物量进行模拟分析。本文基于研究区的实际情况, 假定建设用地的生态系统服务价值为0 ^[22-23], 最终确定符合深圳城市发展与用地特征的单位面积生态系统服务价值表(表 1)。以街道为生态系统服务价值的基本计算单元, 每个街道的生态系统服务价值计算公式如下：

(1)

式中, ESV_f为每个街道第f类生态系统服务价值；A_k为每个街道第k种类型生态系统的面积；VC_fk为每个街道第k种类型生态系统对应的第f类生态系统服务的单位面积价值。

完成研究区各年生态系统服务价值核算后采用K-means算法进行生态系统服务聚类。K-means算法是一种非监督的聚类方法^[24], 这种算法在输入给定的聚类个数k, 以及包含n个数据对象的数据集后, 可以输出满足同一聚类中的对象相似度最高、而不同类别的对象相似度最小的k个类别^[25]。K-means算法虽然存在对初始点敏感度低的问题^[26], 但因其聚类结构明确、聚类过程简单成为地理学空间格局研究中被广泛使用的聚类方法^[27-29]。

因此, 根据上文计算得出各街道每类生态系统服务价值, 求出每类生态系统服务价值占街道生态系统服务总价值的比重, 一个街道一年为一个样本, 五年即得到285(57×5=285)个样本, 每个样本包括11个数据。最后将样本输入K-means算法进行聚类, 进行1—20类逐步迭代聚类, 通过每种聚类情景下的组内离差平方和与组间离差平方和判断出最优的聚类类别数。可以看出当聚类数为6时, 组内离差与组间离差趋于稳定(图 2), 因此最佳聚类数为6。

1.4 生态系统服务簇时空演变轨迹挖掘

在最优聚类类别数目下, 得到57个街道5个时期的生态系统服务簇类型。由此, 每个研究单元将会得到5个类别标签, 这些标签反映出每个研究单元生态系统服务簇类型的时间变化轨迹。本研究采用层次聚类分析(Hierarchical Clustering Analysis, HCA)进行生态系统服务簇变化轨迹的聚类, 具体通过计算每一个数据点与所有数据点之间的距离来确定它们之间的相似性, 距离越小、相似度越高, 由此将距离最近的两个数据点或类别进行组合, 生成聚类树。由于数据点是表示生态系统服务簇类别标签变化轨迹的字符串, 因此采用最小编辑距离(Minimum Edit Distance, MED)测算各轨迹间的距离。MED作为字符串对比的算法, 被广泛应用于计算应用领域^[30], 是指两个字串之间, 由一个转换成另一个所需的最少编辑操作的次数, 距离越大, 说明两个字符串越是不同。

2 结果分析 2.1 生态系统服务簇类型

在最佳聚类情况下得到了6种生态系统服务簇：生态保育型、水源消耗型、水文调节控制型、服务枯竭型、文化支持潜力型、水文调节主导型, 每类生态系统服务簇的结构如图 3所示, 其中水源消耗型的水资源供给为负值。(1)生态保育型：这类生态系统服务簇中服务价值比例超过10%的生态系统服务有4种, 且均小于30%, 因此没有明显主导的生态系统服务, 各类生态系统服务价值在此生态系统服务簇中差距相对其他类型的生态系统服务簇较小。(2)水源消耗型：这类生态系统服务簇的土地利用类型以水田为主, 水资源供给价值为负, 但食物生产价值远高于其他类型生态系统服务簇。(3)水文调节控制型：这类生态系统服务簇中水文调节价值占据绝对优势(比重超过75%), 其他生态系统服务类型的价值均不足10%。(4)服务枯竭型：这类生态系统服务簇的土地利用类型全部为建设用地, 各类生态服务价值均为0。(5)文化支持潜力型：这类生态系统服务簇的文化服务与支持服务价值比重超过20%, 大于除生态保育型的其他类型, 是最有潜力发展为生态保育型的类别。(6)水文调节主导型：这类生态系统服务簇与水文调节关系密切, 水文调节作为簇的主导生态系统服务类型, 其价值占生态服务总价值的比重为50%—75%。

2.2 生态系统服务簇动态变化

图 4是生态系统服务簇在研究时段内的空间分布。可以发现, 虽然不同时期各类生态系统服务簇的空间位置具有一定的变化, 但总体上变化不大。其中, 生态保育型主要分布在深圳中西部和东南部大鹏半岛等森林覆盖率较高的地区；水源消耗型靠近生态保育型与文化支持潜力型, 周边生态环境良好, 为农业生产提供条件；水文调节控制型主要分布在深圳西部和西南地区, 与水文调节主导型和服务枯竭型交错分布；服务枯竭型主要分布在西南沿海地区, 这些地区为深圳的政治、经济中心, 人为活动最为活跃, 土地利用类型均为建设用地；文化支持潜力型分布区域与生态保育型邻近, 位于生态保育型的边缘或过度地带；水文调节主导型主要分布在西部沿海地区, 呈现明显的纵向条带分布。

从时间序列来看, 生态保育型的数量呈现波动减小的趋势, 但始终为6种生态系统服务簇中数量最多的；水源消耗型的数量呈现先减少后增加的趋势, 始终都为6种生态系统服务簇中数量最少的；水文调节控制型的数量呈现逐渐减少的趋势；服务枯竭型的数量变化显著, 呈现出逐渐增加的态势；文化支持潜力型的数量虽然存在年际间的波动, 但较为稳定并没有明显的变化；水文调节主导型的数量呈现先减少后增加的趋势。

可以发现, 除水源消耗型外, 其他5种生态系统服务簇在1980—2015年均始终存在。同时, 有54%的街道生态系统服务簇类型始终保持不变, 说明深圳大部分街道的生态系统服务结构较为稳定, 并未发生明显变化。发生变化的区域大多毗邻香港, 这些地区是深圳经济社会发展最快的地区, 城市土地的快速扩张可能是引起生态系统结构变化的主要原因。

2.3 生态系统服务簇时空演变主导轨迹

基于K-means聚类的结果, 本研究使用1—6的数字字符依次表示生态保育型、水源消耗型、水文调节控制型、服务枯竭型、文化支持潜力型、水文调节主导型这6类生态系统服务簇, 由此57个街道在5个时间节点的类型形成了26条不同的变化轨迹。利用HCA根据每两条轨迹间的MED对轨迹进行聚类, 将MED小于等于2的轨迹聚为一类, 最终得到8种生态系统服务簇时空演变轨迹类型(表 2)。

每种生态系统服务簇时空演变轨迹类型包含若干条不同的轨迹, 同类型的轨迹间具有相似的演变特征。其中, (1)持续水文调节控制型, 数量占深圳6种生态系统服务簇的14.04%, 这种轨迹各阶段的生态系统服务簇类型较为单一, 基本保持为水文调节控制型；(2)水文调节控制向服务枯竭转变型, 数量占深圳6种生态系统服务簇的3.51%, 这种轨迹具有明显的转变趋势, 在1980—2000年间, 生态系统服务簇保持为水文调节控制型, 随后向服务枯竭型转变；(3)转变服务枯竭型, 数量占深圳6种生态系统服务簇的7.02%, 与水文调节向服务枯竭转变型轨迹相似, 这种轨迹也具有明显的转变趋势, 在2000年由其他类型的生态系统服务簇转变为服务枯竭型, 不同的是这类轨迹的生态系统服务簇初始类型不是水文调节控制型；(4)持续服务枯竭型, 数量与持续水文调节型相同, 占深圳6种生态系统服务簇的14.04%, 这种轨迹各阶段生态系统服务簇类型较为稳定, 基本持续保持为服务枯竭型；(5)持续水文调节主导型, 数量占深圳6种生态系统服务簇的5.26%, 这种轨迹的各阶段生态系统服务簇类型基本保持为水文调节主导型；(6)持续生态保育型, 数量最多, 占深圳6种生态系统服务簇的36.84%, 这种轨迹的各阶段生态系统服务簇类型基本保持为生态保育型；(7)生态保育转变型, 数量占深圳6种生态系统服务簇的3.51%, 这种轨迹的生态系统服务簇类型的初始状态为生态保育型, 之后向其他类型转变；(8)持续文化潜力支持型, 数量占深圳6种生态系统服务簇的15.79%, 这种轨迹的各阶段生态系统服务簇类型持续保持为文化潜力支持型。

从轨迹类型上看, 以上8种轨迹类型中有5种属于持续型, 数量占深圳所有街道的86%, 也就是说深圳市大部分街道的生态系统服务簇在近40年内结构相对稳定, 簇内各类生态系统服务变化不大。此外, 还有3种轨迹类型属于转变型, 其中有两种是向服务枯竭型转变, 一种是生态保育型向其他类型转变。由于生态保育型生态系统服务簇中各类生态系统服务均有贡献且较为均衡, 对于地区生态发展来讲, 生态保育型是生态潜力最强的生态友好型生态系统服务簇, 而服务枯竭型的各类生态系统服务价值均为0, 对生态系统的贡献最少。因此, 以上3种生态系统服务簇转变型轨迹在某种程度上均代表了该街道生态环境质量的降级。可以推测, 在近40年的城市发展过程中, 深圳的生态系统服务簇整体上没有提高, 均是保持原来水平或降级, 说明城市在社会经济发展的同时牺牲了一定的生态环境, 并不利于城市生态的可持续维持。

具体来看, 8种演变轨迹类型的空间分布如图 5所示。生态系统服务簇演变轨迹为持续水文调节控制型的有8个, 分布在深圳西部和西南地区；水文调节控制向服务枯竭转变型有2个, 集聚分布在福田区南部；转变服务枯竭型有4个, 零星分布在深圳西南地区；持续服务枯竭型有8个, 呈团状分布在福田区东部和罗湖区西部, 另有一个位于蛇口街道；持续水文调节主导型有3个街道, 两个分布在深圳西部沿海地区, 与持续水文调节控制型相邻。可以发现, 以上生态系统服务簇演变轨迹类型的街道大多分布在西部及西南沿海地区。演变轨迹为持续生态保育型的生态系统服务簇有21个街道, 是数量最多的演变轨迹类型, 主要分布在深圳中西部和东南大鹏半岛两大区域；生态保育转变型有2个, 分别位于龙华区和龙岗区；持续文化支持潜力型有9个, 在深圳东部、西部均有分布, 其中, 位于深圳东部的生态系统服务簇呈现集聚分布的格局, 位于深圳西部地区的生态系统服务簇则为零散分布。和其他轨迹类型相比, 持续生态保育型、生态保育转变型、持续文化支持潜力型这三种演变轨迹类型的生态系统服务簇在空间分布上与其他轨迹类型的生态系统服务簇具有明显的地域差异, 并且占地面积广阔, 占据了深圳72%以上的土地面积。

3 讨论

1980—2015年深圳市不同土地利用类型的面积变化显著(表 3), 其中建设用地的面积明显增加, 由1980年的274 km²增长到2015年856 km², 是所有土地利用类型中面积增长最多, 也是变化最剧烈的；裸地和湿地的面积虽然存在年际间的波动, 但总体变化不大, 其他的土地利用类型都呈现逐年减少的趋势。从各类土地利用类型的转换关系上也可以看出(图 6), 绝大部分地区的土地利用类型并没有发生改变, 主要的变化在于建设用地的增加, 其中大部分新增的建设用地是由水田、阔叶林、旱地转换而来。这解释了8种轨迹类型中持续型轨迹居多的原因：大部分街道的土地利用类型变化较小, 生态系统的结构相对稳定, 因此这些地区的生态系统服务簇类型在长时间序列并没有明显的改变。生态系统服务簇类型发生变化的街道可能是因为建设用地的增加, 由此可以推测城市扩张可能是影响深圳生态系统结构的重要原因。

将生态系统服务簇时空演变轨迹的分布图与深圳市建设用地的变更图进行叠加, 可以判断城市扩张与生态系统服务簇演变的一致性, 如图 7所示。整体来看, 建设用地在深圳各地均呈现扩张的趋势, 这种现象在深圳西南及西部沿海地区尤为明显, 这些地区也是生态系统服务簇转变型轨迹的主要分布地, 中部及东部地区城市扩张的速度较慢, 这些地区大多为持续型轨迹。以新安街道和沙河街道(图 7)为例, 两个街道的生态系统演变轨迹分别为持续型和转变型。新安街道1980年建设用地面积为10.29 km², 到2015年建设用地面积为16.73 km², 虽然有所增加, 但新增建设用地比例不足18%。建设用地的扩张对地区生态系统服务结构的影响并不明显, 因此该街道的生态系统服务簇类型较稳定的维持为水文调节控制型。相对照的, 沙河街道1980年建设用地面积为9.13 km², 到2015年建设用地面积增加到14.56 km², 新增建设用地比例超过34%, 而水田、阔叶林大量减少。建设用地的快速扩张使该街道生态系统服务结构发生明显的变化, 各类生态系统服务价值均减小, 伴随着城市扩张, 该街道的生态系统服务簇类型由水文调节主导型依次变为水文调节控制型、服务枯竭型。这种演变轨迹与城市扩张地域分布的耦合性表明城市扩张与生态系统服务簇演变具有一致性。

4 结论

本文基于生态系统服务簇, 使用K-means算法对深圳市各街道生态系统服务的组合特征进行聚类, 得到生态保育型、水源消耗型、水文调节控制型、服务枯竭型、文化支持潜力型、水文调节主导型等6种类型生态系统服务簇, 在此基础上对区域的生态系统服务簇进行长时间序列变化分析, 得到持续水文调节控制型、水文调节控制向服务枯竭转变型、转变服务枯竭型、持续服务枯竭型、持续水文调节主导型、持续生态保育型、生态保育转变型、持续文化支持潜力型等8种生态系统服务簇时空演变主导轨迹类型。根据轨迹变化特征, 发现深圳大部分地区的生态系统服务结构较为稳定, 发生明显变化的地区主要集中在深圳西南部建设用地扩张最迅速的地区。研究进一步证实, 深圳市生态系统服务簇演变与城市扩张具有一致性。

然而, 本研究仍存在部分有待改进之处。第一, 在计算最小编辑距离时, 本研究视各种生态系统服务簇间相互转化的难易程度相同, 并没有细致区分, 下一步研究中可以根据转化程度的难易对距离进行加权计算。第二, 研究单元的选择会对空间数据分析结果产生不确定性。对于本研究, 虽然栅格尺度精度高, 但确定适合的栅格大小较为困难, 且不便于实际的管理；区级尺度利于政策的实施和管理, 但尺度过大造成地区间差异不明显。因此, 本文以空间数据的精度为代价, 开展基于街道为单元的生态系统服务簇划分。在今后的研究中, 建议对可塑性面积单元问题重点探讨。