长三角城市群地区生态系统服务价值时空演变及驱动因素研究

2022-12-06马伟波赵立君张孝飞张龙江李海东

生态与农村环境学报 2022年11期

马伟波，杨帆，王楠，赵立君，谭琨，张孝飞①,张龙江，李海东

(1.生态环境部南京环境科学研究所，江苏南京 210042；2.陕西省生态环境厅，陕西西安 710004；3.华东师范大学地理信息科学教育部重点实验室，上海 200241)

生态系统服务代表了人类从生态系统中获得的利益[1-2]，通过其生态特征、功能、过程、产品直接或间接促进人类福祉[3]。城市是一个社会-经济-自然复合生态系统[4]，由于快速工业化和经济扩张，城市生态服务供需严重失衡[5]，潜在生态风险不断增加，城市生态福祉问题日益受到关注。城市化与生态系统服务之间存在复杂关系[6-8]，通过核算生态系统服务价值(ecosystem service value, ESV)，剖析ESV空间特征可提升对城市生态系统服务的认识[9-11]，但这对解决更深层次的城市生态问题仍然缺乏有效支撑[12-14]。因此，探索城市群ESV驱动效应对于更好理解城市生态环境问题的空间模式、过程和机制是有必要的[15-17]，也对于提升城市群尺度生态福祉和人居环境健康具有重要借鉴意义[18-19]。

ESV受到生态系统内部和外部的直接或间接影响。其中，自然因素来自生态系统内部，对ESV有最直接的影响。研究发现，气候通过调节地表水热条件进而影响植被生长[20]，生物生长和栖息地的空间特征演变对ESV也至关重要[21-23]，地形可以影响土壤保持力、供水能力和作物产量等多种生态系统服务[24]。同时，随着人类对地球生态环境的影响和破坏，人为因素变得越来越重要[15]。研究证实，土地利用变化反映了人类的资源足迹[25-26]，是人类影响生态系统最直接的表现形式，也是ESV变化的关键影响因素；人口密度、经济发展水平和城市化等因素在不同程度上对ESV也均有影响[15,23]。然而，现有研究对驱动因素的认识还不全面，对驱动关系的讨论并不深入。在驱动方法研究方面，Pearson相关系数能够确定变量间的相关性，但不能指明变量间的明确因果方向；传统回归分析方法不能厘清变量之间相互作用关系，以至于无法精确测度各变量对被解释变量的重要性。因此，目前的研究尚不能阐明ESV各影响因素之间的相互作用关系，各影响因素对ESV的直接或间接效应并没有在统一框架下进行定量测度，这种直接或间接效应的演变特征鲜有被深入探讨。

长江三角洲城市群(Yangtze River Delta Urban Agglomerations, YRDUA)是中国现代化发展和城市化水平最高的地区之一，然而密集的经济发展和快速城市扩张对区域生态福祉获益直接或间接地构成了威胁。基于此，以长江三角洲城市群地区为研究对象，基于修订当量价值系数和土地利用等数据分析长三角城市群地区ESV时空演化特征，并采用随机森林(random forest, RF)和结构方程模型(structural equation model, SEM)方法探究人为活动和自然变化对ESV的驱动特征及驱动路径，以期为城市群一体化生态保护修复和生态福祉提升提供参考。

1 材料与方法

1.1 研究区概况

长三角城市群是长三角一体化发展中心区，位于长江入海之前的冲积平原，根据《长江三角洲区域一体化发展规划纲要》[27]，包括：上海，江苏省的南京、无锡、常州、苏州、南通、盐城、扬州、镇江和泰州，浙江省的杭州、宁波、温州、嘉兴、湖州、绍兴、金华、舟山和台州，安徽省的合肥、芜湖、马鞍山、铜陵、安庆、滁州、池州和宣城27个中心城市(图1)，陆域土地利用面积为22.24万km2。

审图号: GS(2022)4739号

1.2 数据来源

1.2.1ESV核算数据

采用当量因子法核算ESV时，需结合研究区土地利用、社会经济、农作物产值和生态指标对模型参数进行修订。采用净初级生产力(net primary productivity, NPP)数据(MOD17A3H)和降水数据(https:∥crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04/)对长三角地区生态系统服务价值系数进行修正。粮食产值、播种面积主要来源于上海市、江苏省、浙江省和安徽省统计年鉴。长三角城市群陆域土地利用数据(https:∥www.globallandcover.com)显示，2000—2020年长三角城市群呈现出耕地、林地和湿地面积持续下降，而水体和人造地表面积逐步上升的土地利用格局变化趋势(表1)。

表1 长三角城市群土地利用数据

1.2.2ESV驱动数据

为揭示长三角城市群地区ESV驱动效应，选择人为活动和自然条件2个方面10项具体指标作为解释变量(表2)。人为活动驱动类型包括人口密度、夜间灯光、土地利用结构和PM2.5浓度，自然条件包括高程、坡度、归一化植被指数(normalized difference vegetation index, NDVI)、降水、温度和河网密度。夜间灯光被视为反映城市扩张和电力消耗影响的指标，是衡量人类经济活动的有效变量，其能较好地捕捉城市、城镇和工业场所亮度以及城市交通和住宅的扩展特征；研究[28]表明夜间灯光指标是社会经济发展较好的测度指标，笔者采用夜间灯光指标分析社会经济发展对ESV的驱动效应。同时，长三角城市群大部分处于我国东部平原地区，降水充沛，河网密集，传统的土地利用数据由于空间分辨率及地图综合，容易丢失河网密度信息。而河网密度在对土地利用类型进行补充的同时，其对湿地保护也至关重要，因此笔者将其纳入ESV的驱动分析。研究[15]发现，空气环境质量与ESV存在驱动关系，而长三角城市群在2000—2020年存在一定程度大气颗粒物污染，因此将PM2.5选为驱动因素。同时，由于快速城市化导致土地利用结构剧烈变化，生态用地是驱动ESV演变的重要影响因素[15]，因此，将土地利用中的耕地、林地和水域面积等生态用地纳入分析范围。人口密度、高程、坡度、降水和地表温度均为常规驱动因素。所有数据统一为WGS_1984_UTM投影系统，采用渔网方式严格对齐空间格网并提取驱动因素特征值。

表2 长三角城市群生态系统服务驱动数据

1.3 研究方法

1.3.1ESV核算与系数修订

采用当量因子法[1,29]对长三角城市群ESV进行估算，当量因子法适合大区域尺度的ESV核算，其计算过程为

(1)

式(1)中，VES为生态系统总服务价值，元·a-1；Ei,j为第i种生态系统类型第j类生态系统服务价值系数，元·hm-2·a-1；Ai为第i种生态系统类型面积，hm2。同时，采用地区农作物产值面积、NPP和降水数据修订价值系数[30]。无人力投资的自然生态系统输出价值为农田食物生产服务单价的1/7，因此修订标准当量因子公式为

(2)

式(2)中，E为单位面积粮食产量的经济价值，元·hm-2；T为研究区粮食总价值，元；X为研究区粮食播种面积，hm2。进一步依据条件因子法利用NPP和降水数据对部分价值系数进行修订，计算公式为

λ=B/B0，

(3)

Ej=λ·E0j。

(4)

式(3)～(4)中，λ为生态系统服务当量的区域修正系数；B0为全国平均NPP数值；B为长三角城市群逐像元NPP，g·m-2·a-1；Ej为第j类生态系统经过区域修正后的生态系统服务价值当量；E0j为第j类生态系统全国平均生态系统服务当量，元·a-1；j=1,2,…,8，依次为食物生产、原料生产、气体调节、气候调节、净化环境、维持养分循环、维持生物多样性和提供美学景观服务。在利用NPP对8种生态系统服务价值系数进行逐像元修订的基础上，再通过降水产品数据对水资源供给和水文调节服务价值当量继续进行修订。2000—2020年长三角地区农作物种植面积为1 047.24万hm2，对应总产值为3 465.93 亿元，则得到单位面积粮食产量的经济价值为4 727.96元·hm-2。通过统计得到研究区不同土地利用生态服务功能价值当量(表3)。

表3 长三角城市群单位面积生态服务价值当量

1.3.2随机森林

RF是一种集成学习方法[31]，能用于充分探测变量对解释变量的贡献程度[32]。训练一个RF的本质是训练多个分类与回归树(classification and regression trees, CART)。CART是一个二叉树模型，其核心是切割变量和切割点的选择[32]。在RF中，单棵CART首先遍历一部分变量和变量数据，然后根据切割后节点的不纯度，确定最佳切割变量和切割点，综合所有树的结果得到最终模型。选择RF回归方法，节点不纯度的计算公式为

(5)

同时，RF通过袋外误差(out-of-bag error, OOB)样本来估计变量的重要性，计算公式为

(6)

式(6)中，I(vari)为变量i的重要性；Eij1为根据CARTj中变量i出袋数据计算的误差；Eij2为根据CARTj中变量i出袋数据加上噪声干扰计算的误差；n为CART数量。10项驱动因素对ESV响应的解释程度采用RF方法进行建模，方差var解释程度为RF模型的解释能力，最大为100%；采用I(vari)测度各驱动因素对ESV的相对重要性。

1.3.3结构方程模型

SEM与回归分析的重要区别是在SEM中，某个响应变量(回归分析中也叫因变量)还可以作为其他响应变量的预测变量(自变量)。换言之，SEM包含了多个变量之间直接和间接的因果关系，并通过图示中的箭头(路径)方向展示所有的因果关系[33-35]。与传统的基于方差-协方差的SEM相比，分段SEM[33-35]具备以下3个优点：(1)将多个独立线性模型组合成一个单一的因果网络；(2)使用Shipley分离检验来检验模型中是否缺少任何路径；(3)使用Akaike信息准则(AIC)来比较嵌套模型。

采用分段SEM研究10项指标对长三角城市群ESV的驱动路径，同时也考虑以林地面积和NDVI为中间路径间接驱动ESV的情形。分别采用线性模型和拟合分组模型，对驱动因素的驱动路径进行标准化处理，最终采用Shipley分离检验评估分段SEM的整体拟合解释度。采用R 4.1.2软件piecewiseSEM工具实现分段SEM。参考前人研究[36-37]，采用Chi-Squared、Chi-Squared/P-value、Fisher′s C、Fisher′s C/P-value、AIC和R2作为模型检验指标，其中，Fisher′s C检验值低于显著水平(P<0.05)表明模型假设与数据不一致，应予以拒绝。

2 结果与分析

2.1 长三角城市群ESV时空特征

2.1.1时空变化

由表4可知，2000—2020年长三角城市群ESV总体呈现先下降后上升的变化趋势，由2000年的24 508.53亿元下降到2010年的24 285.57亿元，再上升至2020年的25 217.34亿元，2000—2020年ESV整体增加708.81亿元，增幅为2.89%。2000—2020年，研究区食物生产、原料生产、气体调节、气候调节、维持生物多样性、维持养分循环和美学景观类型生态服务价值均呈现下降趋势，而水资源供给呈减少损失趋势，水文调节呈上升趋势，这主要是由于耕地和林地面积下降而水域面积上升驱动所致(表1)。具体来看，2000—2010年，水文调节服务价值下降120.52亿元，而水资源供给服务价值上升75.81亿元。2010—2020年，水文调节和水资源供给服务价值增加1 096.2亿元，是整体ESV提升的主要贡献者。

表4 2000—2020年长三角城市群生态系统服务价值(ESV)变化情况

从空间变化(图2)来看，长三角城市群形成重点湖泊、湿地和水系地区ESV最高，南部丘陵地区ESV次之，北部农田地区ESV再次之，都市群建成区ESV最低的空间分布格局。由图2可知，2000—2020年合肥市巢湖北部、杭州市东部、盐城市北部、苏州市东部和北部、无锡市北部、常州市、扬州市、南京市、宁波市和金华市均出现ESV下降趋势。泰州市、扬州市北部和南通滨海湿地ESV有所增加；苏州市城区ESV呈集聚态势发展，离散化程度有所降低，ESV有所增强；太湖、巢湖、高邮湖、石臼湖、滆湖和长江干流水系以及滨海湿地等地区ESV呈稳定发展态势。

图2 长三角城市群ESV时空变化

2.1.2空间异质性

由图3可知，长三角城市群在不同空间分辨率尺度下ESV空间Moran′sI指数均大于0.45,Z得分大于1.96(均通过1%显著性检验)，表明长三角城市群ESV有较强的空间聚集模式。具体来看，在1和5 km空间分辨率下，ESV的Moran′sI指数和Z得分测度接近，2000—2020年Moran′sI指数平均值分别为0.83和0.78，Z得分平均值分别为553.23和541.76，均高于其他空间分辨率，表明1和5 km是ESV空间异质性表现最为强烈的空间尺度。同时发现，在不同空间分辨率尺度下，2000、2010和2020年ESV空间异质性强度呈微弱下降趋势，这表明长三角城市群地区ESV的空间离散性特征随城市化剧烈演变而逐渐增强，生态空间格局呈破碎化和离散化。

图3 不同空间分辨率下长三角城市群ESV空间异质性

2.1.3冷热点分析

综合空间异质性表达能力和计算工作量，确定5 km分辨率作为研究区ESV驱动分析的空间单元。由图4可知，2000—2010年，安庆市南部、池州市北部、南通市与上海市接壤地区、苏州市与无锡和上海市接壤地区、杭州市与湖州和嘉兴市接壤地区以及舟山市出现空间冷点，表明这些地区ESV下降呈现出空间聚集效应；而常州、无锡和苏州市南部地区，南京、马鞍山、芜湖和宣城市接壤地区，铜陵、池州、安庆、泰州北部、滁州和合肥市接壤地区，宁波、台州和温州市滨海地区出现空间热点，表明这些地区ESV上升呈现出空间聚集效应。与2000—2010年相比，2010—2020年(图4)扬州、泰州和盐城市接壤地区，盐城和南通市滨海湿地ESV上升空间聚集效应呈扩大趋势；而宁波、台州和温州市滨海地区，滁州和合肥市接壤地区，南京、马鞍山、芜湖和宣城市接壤地区则由2000—2010年的空间热点转换为空间冷点，杭州市与湖州和嘉兴市接壤地区、南通市与上海市接壤地区由2000—2010年的空间冷点转换为空间热点，这表明2010—2020年这些地区ESV变化呈现剧烈波动。整体来看，2000—2020年扬州、泰州和盐城市接壤地区，盐城和南通市滨海湿地，常州、无锡、苏州和湖州市接壤地区，南京和马鞍山市接壤地区ESV上升热点效应非常显著，而上海、苏州、南通、无锡、镇江和南京6个城市以及合肥、芜湖和铜陵3个城市建成区ESV下降冷点效应非常显著，长三角城市群ESV上升和下降的空间聚集效应均非常显著。

图4 2000—2020年长三角城市群ESV变化冷热点

2.2 长三角城市群ESV驱动因素

2.2.1相关性分析

Pearson相关性分析(图5)显示，2000—2020年土地利用结构对ESV的相关性整体较高，水体面积与ESV的相关性稳定在0.8，林地面积则从0.5下降到0.3，耕地面积从-0.9变化到-0.8，PM2.5从-0.7 变化到-0.8。由于长三角城市群大部分地区是平原，因此地形因子与ESV相关性极低。河网密度与ESV的相关性在-0.5～-0.4之间波动。降水、人口密度、NDVI和夜间灯光与ESV的相关性不高。此外，耕地面积与林地和水域面积呈显著负相关，而耕地面积与河网密度和PM2.5呈显著正相关；温度与降水和NDVI呈显著正相关；夜间灯光与人口密度呈现一定程度正相关。

图5 2000—2020年长三角城市群ESV与驱动因素相关性

2.2.2重要性分析

表5显示，2000、2010年2020年RF模型的方差解释度分别为91.53%、91.83%和89.97%，均超过89%，这表明RF模型通过驱动因素可以实现对ESV较强的解释，因此可基于该模型对驱动因素变量重要性进行测度。

表5 长三角城市群ESV随机森林模型精度评价

由图6可知，水域面积、林地面积、耕地面积、PM2.5和河网密度是2000—2020年长三角城市群ESV重要驱动因素，且水域面积重要程度远远高于其他影响因素，平均值达到1 781。

图6 2000—2020年长三角城市群ESV与驱动因素相关性

2000—2020年，水域面积、林地面积、PM2.5和河网密度重要性均呈上升趋势，耕地影响因素重要性呈微弱下降趋势。高程、坡度、降水、人口密度、夜间灯光、温度和NDVI等影响因素重要性较低，且无明显变化趋势。

2.2.3驱动路径

由SEM分析结果(表6)可知，所有年份Chi-Square和Fisher′s C检验的P-value均大于0.05，表明因果关系假设通过检验，SEM模型可信度高。由图7可知，总体来看，SEM模型评价结果显示，水域面积、林地面积、耕地面积、河网密度、PM2.5和人口密度等驱动因素对2000、2010和2020年ESV以及2000—2020年ESV变化的解释程度分别为85%、84%、83%和72%，表明SEM模型精度和解释能力较高。

表6 长三角城市群SEM精度评价

由图7可知，2000年水域面积、林地面积和PM2.5对长三角城市群ESV有直接正向驱动效应〔标准化路径系数(β)分别为0.79、0.43和0.05〕，耕地面积、河网密度、夜间灯光和人口密度则为直接负向驱动效应(β分别为-0.21、-0.06、-0.03和-0.01)。PM2.5对林地面积呈显著直接正向驱动效应，β值为0.24；耕地面积、水域面积、河网密度、夜间灯光和人口密度对林地面积均呈负向驱动效应(β分别为-0.72、-0.27、-0.06、-0.03和-0.02)。对NDVI呈直接正向影响的因素包括温度、坡度和林地面积(β分别为0.41、0.1和0.03)，而人口密度和水域面积呈负向影响(β分别为-0.04和-0.03)。PM2.5、人口密度和夜间灯光同时也通过林地间接影响ESV。而高程、温度和坡度对长三角城市群ESV和林地面积的驱动效应并不显著。

与2000年相比，2010和2020年长三角城市群ESV驱动模型差异较小(图7)，水域面积和河网密度对于ESV的作用强度有所增强，β分别从0.79上升到0.81和从-0.06变化到-0.11；林地面积β从0.43下降到0.41，耕地面积从-0.21变化到-0.16；PM2.5和夜间灯光对ESV的直接负向驱动影响显著下降，而高程和NDVI对ESV的正向驱动效应显著上升。耕地面积和PM2.5对林地的负向驱动效应下降，而水域面积和河网密度对林地面积的驱动效应则增强。由图7可知，2000—2020年ESV变化的显著正向驱动因素为水域面积、林地面积和PM2.5(β分别为0.83、0.4和0.05)，负向驱动因素为河网密度，β值为-0.06；耕地面积、NDVI、人口密度和夜间灯光对ESV变化的驱动不显著。同时，人口密度和夜间灯光通过驱动林地面积间接影响ESV变化的效应不显著；而PM2.5驱动林地面积则呈显著负向效应。

黑色箭头表示预测因子对反应变量有显著影响，蓝色虚线表示无显著影响。箭头方向表示因果作用方向，箭头上的数字为标准化路径系数(β)。

3 讨论

笔者基于修订的当量价值系数和土地利用数据分析长三角城市群ESV时空演化特征，采用RF和SEM模型探究人为活动和自然条件共10项指标对ESV的驱动特征及驱动影响演化特征，在统一框架下定量测度ESV各影响因素之间的相互作用关系及其对ESV直接和间接驱动效应，尝试探究长三角城市群ESV演化模式和驱动机制。

3.1 ESV演化特征

长三角城市群地形结构影响了土地利用结构基本格局特征，进而对长三角城市群ESV空间分布格局产生主要影响。研究区东部平原河网密集，水域湿地生态保护和土地利用结构是主导ESV变化的主要影响因素。2000—2010年，长三角城市群ESV下降主要原因是耕地和水域湿地等生态用地面积迅速下降，这与同时期该地区高速经济扩张和城市化建设是相吻合的，在上海、苏州和无锡市接壤地区，南通和上海市接壤地区，杭州、嘉兴和湖州市接壤地区ESV变化冷点效应非常显著(图4)。2010—2020年，南京、镇江、马鞍山和芜湖市连片接壤地区，苏州和南通市接壤地区，台州和温州市滨海地区ESV变化由2000—2010年的热点变为冷点或冷点效应显著(图4)，这从侧面表明这些地区是长三角城市群演化的新兴高速发展区域。

3.2 ESV驱动路径

2000—2020年各驱动因素的重要性(图6)以及驱动因素对ESV标准化路径系数显著性及大小的差异(图7)体现了各驱动因素在ESV及影响因素的演化过程中的分异特征，这种特征对理解ESV演化过程具有重要启示。在直接驱动路径方面，与LUO等[14]采用岭回归解释的城市群驱动因素结果相似，土地利用结构对研究区ESV驱动力非常显著，但笔者发现水域面积是ESV最主要的驱动因素，这与WANG等[37]的发现相似。笔者进一步阐明了水域、林地和耕地面积对ESV的直接影响以及水域和耕地面积通过林地面积对ESV的间接影响，剖析了完整驱动过程(图7)。同时，笔者选择河网密度和PM2.5作为长三角城市群ESV驱动因素是结合当地生态环境和社会经济条件的一种探索。由图7可知，河网密度对ESV的驱动影响非常显著，且标准化路径系数绝对值随时间呈上升趋势。同时，2000—2020年，长三角地区PM2.5呈下降趋势，其对ESV的驱动显著下降，对林地面积的驱动强度也呈下降趋势。

在间接驱动路径方面，NDVI和林地面积同时作为中间变量直接或间接作用于ESV，其中，其他驱动因素通过林地面积对ESV的间接作用则更为显著；同时，水域面积、耕地面积、河网密度、PM2.5、人口密度和夜间灯光对林地面积的驱动效应均非常显著，但对林地面积变化的解释程度则不是很高，为47%～55%。驱动因素直接和间接标准化路径系数的分异特征也体现出生态系统复杂的作用关系。笔者研究中，除了耕地面积通过林地面积间接作用于ESV的间接标准化路径系数小于耕地面积直接作用于ESV的标准化路径系数外，其他驱动路径均为间接驱动路径的显著性水平和驱动强度大于直接驱动路径，如PM2.5和河网密度作用于林地面积的显著性水平和标准化路径系数大于直接作用于ESV。

笔者研究对NDVI的解释程度较弱，低于YANG等[38]发现的江苏省人类活动变化对NDVI的解释强度为56%的结果。但笔者研究发现，在气候变化方面，温度对NDVI的标准化路径系数β为0.36～0.41，这可能是由于长三角城市群南部地区多为丘陵地貌，植被变化受到光热温度条件影响较大。

3.3 局限性

尽管在SEM驱动建模的ESV终端，SEM模型通过了假设检验且拟合精度较高，但其对建模路径中的NDVI和林地面积的解释能力较弱，这表明对于NDIV和林地面积仍然存在其他被忽视的驱动要素。同时，基于当量因子核算的ESV评价方法在一定程度上夸大了土地利用结构对ESV的驱动影响，后期应尝试采用其他ESV核算方法进一步验证驱动因素的重要性和驱动关系。在驱动指标选择中，建设用地扩张导致生态系统服务价值减少，也是重要的驱动因素，而且其与城市扩张代表性指标夜间灯光的相关关系也应该纳入分析；目前采用当量法测度的ESV，忽视了更加细分尺上度生态系统健康的内容，如原生森林生态系统单位面积ESV与人工林生态系统存在较大差异，这些可以在后期研究中通过试验分析完善。尽管在2.1节从空间角度阐述了ESV的空间特征，但2.2节分析是从统计学习角度开展分析讨论，并不能从空间角度深入分析，后期应注重从时空演变和地区政策实施效果角度深入阐述ESV演变时空效应。

同时，在区域尺度下不同生态功能区划主导生态功能存在差异，主导生态功能对生态系统提供的供给、调节、支持和文化等服务价值存在差异，则其对ESV演化的驱动关系网络必然存在分异特征，这是笔者研究未能考虑到的方面；而且，景观格局、交通网络、空气污染物(除PM2.5外)和土壤类型等均对ESV产生驱动作用[15,25,39]，以及受污染的水体、耕地和林地对ESV驱动的差异影响，在城市发展过程中城市-郊区-农村生态空间梯度演变特征对ESV演变也存在复杂影响，这些可以在后期研究中进一步探索。