生态系统服务被广泛定义为直接或间接促进可持续人类福祉的生态特征、功能或过程, 并按照生态系统对人类产生的直接影响分为供给、调节、文化服务, 以及对维持生态系统的其他功能具有重要作用的支持服务^[1]。当前, 伴随着生态系统服务对人类发展的重要性被广泛认识, 探求城市群生态系统服务的演变规律与驱动因子关系逐渐成为学界的关注点之一。相关研究主要将生态系统服务驱动因素划分为自然因素和社会经济因素两方面。在自然环境方面, 纪建万、闫国振、苏常红等学者研究发现, 气候、土壤、高程等自然气候因素变化与生态系统服务价值演变存在明显相关^[2—4]。在社会经济方面, 徐煖银, 耿甜伟、璩路路等学者发现社会因素对生态系统服务价值有较大的影响, 并在空间上与生态系统服务具有高度相关性^[5—8]。由此可见, 自然环境因素和经济社会因素变化对生态系统服务有显著影响和密切联系。因此, 识别城市群生态系统服务的主要驱动因素, 可以为未来规划者在城市群发展过程中做出决策提供依据。

2018年11月中共中央、国务院印发《关于建立更加有效的区域协调发展新机制的意见》, 提出“以北京、天津为中心引领京津冀城市群发展, 带动环渤海地区协同发展”的目标。环京津地区在自然资源的提供、废气废物的处理等方面为京津发展提供支持, 是京津地区生态安全的绿色屏障和重要生态缓冲地带, 因此对京津及环京津六市进行生态系统服务评价和驱动因素分析, 可以为未来协调该地区生态保护与社会经济发展提供支撑, 并更好地促进京津冀一体化发展^[9—11]。

同时, 为避免当前广泛使用的“基于单位服务功能价格的方法”和“基于单位面积价值当量因子的方法”分别存在的忽略生态系统服务在时空尺度上的内在丰富信息, 不利于对生态系统保护的管理决策^[12—13]和可以转化为市场交易价格的生态系统服务指标较少, 限制了对区域关键生态系统服务的评估^[14—15]的问题。本次研究在系统梳理研究区域的发展脉络的基础上, 选取长时间跨度中的多个时间节点, 并以当前京津及环京津六市面临的生态挑战作为研究指标构建符合区域特征的生态系统服务评价体系, 生成生态系统服务综合指数。通过PCA-OLS-GWR空间统计分析, 较为科学地分析研究时期内, 生态系统服务演变的特征, 辨析其发展的动态驱动机制。通过分析实现主要研究目标如下：(1)分析2000—2018年间生态系统服务的时空演变特征；(2)构建自然环境因素和社会经济驱动因子指标体系, 阐释各项驱动因子对生态系统服务演变的驱动机制；(3)结合政策导向提出区域生态系统服务优化提升路径。

1 实验数据与方法 1.1 研究区域与研究数据

本次研究将北京、天津以及环京津六市(承德、唐山、张家口、保定、沧州、廊坊)8个城市划定为研究区域, 总面积16万km²。该区域地处东经113°40′—119°19′, 北纬37°29′—42°37′, 位于黄河下游以北, 北靠燕山山脉, 南朝华北平原, 西倚太行山, 东临渤海湾, 属于典型温带季风气候。地势由西北向东南逐渐倾斜, 降水量自东南向西北递减, 地貌类型和生态系统类型多样, 从西北向东南依次分为坝上高原区、燕山和太行山区、京津保中心区、海岸海域区以及低平原区五大区域(图 1)。2018年常住人口为7486.8万人, 2018年GDP总值高达68081.46亿元, 是中国北方地区经济发展、人口聚集的城市群, 在中国经济社会发展中具有重要的战略地位。

研究数据包括基础地理信息数据、土地利用数据、遥感(NNP、NDVI)数据、数字高程(DEM)数据、气象数据、土壤属性数据、社会经济统计数据7类。基础地理信息数据及行政区划、土地利用数据、温度、降雨量等气象数据和土壤属性数据均来源于中国科学院资源环境科学数据中心(www.resdc.cn)；遥感数据和高程数据来自国家地理信息空间数据云平台(http://www.gscloud.cn/)；人口等社会数据来源于2019年北京、天津及河北各市县的统计年鉴；基础交通设施POI矢量数据来自于世界数据地图(www.openstreetmap.org), 以上数据分辨率均为1km×1km。通过对比1998—2018年NDVI的平均值, 根据突变点并结合研究区域的发展历程筛选出具有代表性的4个时间点——2000年、2005年2015年和2018年, 划分3个时期进行研究。

1.2 研究方法

相关学者的研究表明, 研究区域近年来面临着植被减少、水土流失、水源匮乏、城市热岛、动物多样性等一系列问题^[16—19]。因此, 本研究根据“千年生态系统评估”框架并依据当前研究区域生态系统的特征和面临的生态压力, 从供给服务、调节服务、支持服务和文化服务4种服务类型出发, 综合净初级生产力(NPP)、水土保持、固碳释氧、水源涵养、热岛缓解、生物多样性保护、游憩利用7种指标表征区域生态系统服务整体水平(表 1)。通过ENVI遥感解译、地理信息叠加分析等技术进行生态系统服务评估。

1.2.1 净初级生产力(NPP)

净初级生产力(NPP)是绿色植物呼吸后所剩下的单位面积单位时间内所固定的能量或所生产的有机物质。公式如下：

式中, NPP单位为：g m^-2 a^-1, R_a表示自养呼吸的消耗量, GPP表示总初级生产量。

1.2.2 水土保持

水土保持是生态系统通过其结构与过程减少由于水蚀所导致的土壤侵蚀的作用。公式如下：

式中, S_pro为水土保持服务能力指数, NPP_mean为植被净初级生产力平均值, F_slo为坡度因子, K为土壤可蚀性因子。

1.2.3 固碳释氧

固碳释氧是生态系统通过绿色植物的光合作用进行氧气和二氧化碳的物质交换过程, 生态空间每生产1g的植物干物质, 可以固定1.63g的CO₂, 同时会释放出1.2g的O₂。公式如下：

式中, G_v为植被年固碳量, 单位为g/d, G₀为生态空间年释氧量, R_c为CO₂中碳的含量, 取值27.27%, A为生态空间面积(hm²), NPP为植被净初级生产力平均值。

1.2.4 水源涵养

水源涵养是生态系统通过其特有的结构与水相互作用, 对降水进行截留、渗透、蓄积, 并通过蒸散发, 实现对水流、水循环的调控。公式如下：

式中, WR为生态空间水源涵养服务能力指数, NPP_mean为植被净初级生产力平均值, F_sic为土壤渗透因子, F_pre为平均降水量因子, F_slo为坡度因子。

1.2.5 热岛缓解

热岛缓解是通过增加植被、水体及湿地等生态系统, 减轻由于人为原因导致的城市温度高于周边地区的热岛效应。公式如下：

式中, T_s为地表温度, 单位为℃/a, T_i和T_j为热通道i和j的亮度温度, A和B为系数, 由大气透过率和地表比辐射率等因子决定。

1.2.6 生物多样性维护

生物多样性维护是生态系统支撑服务的重要内容之一, 旨在维持物种、基因多样性中发挥的作用。通过该项指标的评价, 识别现状和未来不同区域承担生物多样性保护的能力强弱。公式如下：

式中, S_blo为生物多样性维护服务能力指数, NPP_mean为植被净初级生产力平均值, F_pre为平均降雨量, F_tem为平均气温, F_alt为海拔因子。

1.2.7 游憩利用

游憩娱乐为人类生活提供了休闲、娱乐与美学享受。通过对游憩资源、游憩设施、游憩区位三方面多个指标进行评价, 判断生态系统文化服务的能力水平^[12]。公式如下:

式中, S为游憩利用综合评价指数, W_i为游憩因子等级值, X_i为不同游憩因子权重值, n为评价因子的个数。

1.3 生态系统服务综合指数

为避免同一类型服务因指标因子多而造成数值过高的问题, 本次计算首先对上述7项服务指标进行最大值-最小值归一化, 然后按照生态系统服务4大分类等权叠加后再次进行归一化, 最终将四项服务指标累加获取生态系统服务综合指数, 表征区域生态系统服务的总体水平。具体公式如下：

式中, y_ij、x_ij为生态系统服务因子i的标准化数据和原始数据; x_max和x_min分别为状态指数j的原始数据的最大值和最小值。

CIES为生态系统服务综合指数; ES为第i类生态系统服务的标准化结果; n为4。

1.4 因子筛选与相关性分析 1.4.1 普通最小二乘法(Ordinary Least Squares, OLS)

OLS模型作为全局回归模型, 常用于解释单一因变量与多个自变量之间的驱动关联。本研究将其作为筛选影响城市群生态系统服务相关因子的评价途径, 计算生态系统服务价值与各项驱动因子基础回归关系, 并进行因子筛选。计算公式如下：

式中, β₀是常数项；β_i为回归系数；ε_i为随机误差项

1.4.2 地理加权回归(Geographical Weighted Regression, GWR)

在建立OLS模型的基础上, 通过GWR模型建立嵌入空间加权函数的线性回归模型, 用以探索研究对象在某一尺度下的空间变化关系及驱动因素。相比使用OLS模型进行回归分析, GWR模型能够反映参数在不同空间的非平稳性, 使变量间关系可以随空间位置变化, 同时考虑空间对象的局部效应, 其结果具有更高的准确性^[20—21]。其评价公式如下：

式中, β_{0(Ui, Vi)}是(U_i, V_i)空间位置的地理加权回归截距, β_{k(Ui, Vi)}是k项自变量(驱动因子)在(U_i, V_i)空间位置的加权回归系数, xk(U_i, V_i)为第k项自变量(驱动因子)在(U_i, V_i)空间位置的取值, ε_i为算法残差, XT是自变量(驱动因子)的转置, W(U_i, V_i)是距离权重矩阵, h为AIC准则的带宽, 为空间i位置与空间j位置的距离。

2 研究结果 2.1 生态系统服务时空演变特征

对研究区域生态系统服务指数进行具体分析发现, 4类服务指数在各个时期表现出的特征各异(图 2), 具体情况如下：1)供给服务指数在2000—2018年间总体上呈现从上升向平稳过渡的特征, 并在2015年达到研究时期最大值0.316。2)调节服务指数在2000—2018年间呈现波动逐渐增强的变化特征, 在2015年达到研究时期最高值0.333后并在2018年骤降至0.292。3)文化服务指数在2000—2018年间在总体上呈现先上升后下降末期数值高于初期的特征, 在2005年达到研究时期最高值0.057。4)支持服务指数在2000—2018年间变化幅度较大, 总体呈现先上升后下降再上升的特征。在2005年达到研究时期最高值0.057并在2015年降至研究时期最低值0.020, 整体波动较大。

通过对研究区域在2000年、2005年、2015年、2018年生态系统服务综合指数进行综合量化分析, 发现研究区域生态系统服务综合指数呈现出先上升后下降、研究期末生态系统服务综合指数略高于研究初期的总体特征。

通过对4类生态系统服务指数在空间变化上进行研究后发现(图 3)：1)供给服务指数在研究期间内呈现以太行山燕山山脉为界东南、西北两侧分别以“增强-减弱-增强”和“减弱-增强-减弱”的间隔规律差异发展的总体特征。2)调节服务指数在京津保中心区和燕山区域内逐渐增强, 但在坝上高原和西部太行山区呈现相反演变特征, 而海岸海域区、低平原区以2015年为界, 先增强后降低。3)文化服务指数在研究期间, 以2005年为界, 燕山和太行山区、京津保中心区、低平原区呈现先增强后减弱变化趋势。海岸海域区则呈现相反的态势。坝上高原区在研究期间变化趋势不明显。4)支持服务指数, 除了在太行山区呈现先增强后减弱再增强的变化波动现象外, 整体变化不明显。

生态系统服务综合指数的分级方法是在地理信息系统软件中将4类服务做叠加处理, 得到生态系统服务综合指数空间分布。同时, 进一步采用自然间断点分类对其进行5级分类。通过对在研究期间内的生态系统服务综合指数进行评估得出(图 4)：研究期间内, 总体上燕山山脉及太行山山脉北部区域为研究区域生态系统服务综合指数的高值区域, 坝上区域为低值区。生态系统服务综合指数波动主要集中在太行山和燕山区, 燕山区域生态系统服务综合指数在2000—2015年间呈下降的趋势, 2015—2018年间呈上升的趋势, 而太行山区域在整个研究时期呈现下降趋势。

2.2 构建驱动模型

为了探索自然社会因素与生态系统服务之间的关系, 本研究在可获取的数据范围内, 结合研究区域实际情况^[22—24]和相关研究成果^[25—28], 最终选择涵盖气候、地形土壤、人类活动及土地利用5类共8项因子代表区域自然地理特征和社会经济发展的特征(表 2)。这些因子作为研究区域生态系统服务的预选驱动变量, 并确定最邻近状态的2018年进行驱动要素分析评估。为避免自变量出现明显的共线性问题, 导致因子对模型的解释力出现偏差, 研究以拟合优度大于0.5、P值边界0.05、冗余检验小于7.5作为筛选条件, 运用OLS模型筛选出5个驱动因子, 即降雨量、温度、土壤含沙量、人口密度和基础设施密度。在对比OLS模型与GWR模型的基础上(表 3), 将驱动因子导入GWR模型中进行地理空间加权回归分析。模拟结果在ArcGIS 10.6软件中进行空间可视化, 获得的GWR模型拟合结果显示图(图 5)。在拟合优度(R²)模拟结果上, 除保定南部极少部分区域外整体拟合优度均高于0.6, 平均值为0.655, 模型效率可信；在拟合t值模拟结果上, 整体区域模型t值范围均在-2.58—2.58内, 模型结构可信。

2.3 驱动结果分析

全局回归结果表明(表 4), 各因子的驱动强度排序如下：降雨量(5.244405)>土壤含沙(-2.185510)>人口密度(-1.009130)>温度(0.788722)>基础设施密度(-0.649611)。由此可见, 与社会经济因素相比, 自然环境因素对研究区域生态系统服务能力影响更大。GWR回归结果如图 6所示。在自然驱动因子中, 降雨量对生态系统服务的驱动影响呈现明显的正相关性, 驱动作用强度从承德西北部和张家口东南部山地向南呈现递减的趋势。温度对生态系统服务驱动影响呈现由北向南从正相关向负相关过渡的趋势, 正向驱动作用在保定及廊坊以北逐渐增强, 在承德则出现最为明显的正驱动作用, 而在以南地区呈现负相关。土壤含沙量与研究区域生态系统服务整体呈负相关性, 土壤含沙量的增加会给生态系统服务带来负面影响, 且从北向南逐渐增强。南部地区保定、沧州市较中北部城市呈现的负驱动作用最为明显。人口密度和基础设施密度2项社会驱动因素都对生态系统服务驱动呈现负相关。人口的增长会给生态系统服务能力带来负面影响, 且从东北向西南逐渐增强, 南部地区保定、沧州市较中北部城市呈现出更为明显的负驱动作用。基础设施密度对生态系统服务的负向驱动作用由太行山脉部向东南和东北方向逐渐增强, 以坝上高原区及海岸海域区负向驱动作用最为明显。

3 讨论 3.1 自然-社会因子对生态系统服务的驱动作用分析

人为因素和自然因素是生态系统变化的主要驱动力因素, 在长时间跨度上, 自然因素对生态系统服务的影响是巨大的, 而在短时间内, 生态系统服务则更易受人为因素影响^[29]。京津及环京津六市是京津冀城市群中人口最为稠密、土地利用方式复杂、经济生态社会胁迫效应显著的区域。在研究期间内, 该区域的生态系统服务呈现明显且明确的空间分异, 即在空间维度上, 生态系统服务的高值区域集中在燕山山脉区域, 而生态系统服务低值区域集中在坝上高原地区。本研究运用“OLS降维模型+GWR回归模型”探索了驱动因素对生态系统服务的空间异质性影响^[30—31], 使得研究结果更加全面且符合区域实际。GWR回归结果表明, 由于研究区域位于属于中国半干旱向半湿润过渡的地带北方农牧交错带, 水资源较为匮乏, 生态本底较差^[32]。因此, 降雨量对该区域生态系统服务的正向演化具有积极作用, 并对较为干旱的坝上高原区域表现出更强的驱动力。同时温度也是限制植物生长的主要因素^[33], 由于太行山区及坝上高原区气候寒冷, 而南部平原地区较为温暖, 所以温度对生态系统服务的正向驱动作用由南向北增强。其次, 土壤的颗粒组成及含沙量是影响土壤可蚀性差异的主要因素。土壤中其含沙量越大, 土壤可蚀性越高, 易发生侵蚀^[34]。而在研究区域内, 南部地区用地性质以农田为主, 植被单一加之人类不合理的生产活动, 导致土壤养分流失和土地退化以及水土流失等生态问题, 从而影响生态系统的稳定性^[35—37]。社会因子对生态系统服务的影响也较为明显, 人口密度和人类活动强度的增加对生态系统产生较强的干扰和破坏作用^[38—39]。人口增长会增加对周边自然资源需求和改变土地覆盖变化, 对周边生态环境造成不利影响^[31]。本研究区域北部山地城市人口密度较小, 基础设施建设量也相对较少, 对自然系统的破坏较小。而中部南部平原城市由于北京首都功能疏解导致的人口外溢和人口的高流动性, 使该区域人口聚集, 对生态系统服务的负驱动作用更为显著。由此可见, 人类活动多集中在易于建设的平原地区, 社会因素对平原的生态系统服务负向驱动力明显高于山区与高原地区。

3.2 生态系统服务时空演变分析

本研究通过总结研究区域当前面临的关键挑战, 选取确定了净初级生产力(NPP)、水土保持、固碳释氧、水源涵养、热岛缓解、生物多样性保护、游憩利用7种评价指标, 生成该区域生态系统服务综合指数, 与相关学者对该区域的生态系统服务时空分布的研究结果相似^[28], 但由于在指标选取方面存在差异, 生态系统服务的空间分布也存在一定的差异, 同时相比与运用生态系统服务价值时空分布, 该实验结果可以更加准确地反应时空尺度上的内在丰富信息。通过对2018年研究区域生态系统服务驱动因子的分析, 结合研究区2000—2018年生态系统服务指标时空演变的特征, 对近20年的生态系统服务的驱动因子作用机制进行分析：2000—2015年间, 退耕还林等生态治理工程、京津风沙源治理工程等生态治理工程在研究区域内不断推进。这些生态保护项目注重地表覆盖的绿化, 植被覆盖的明显提升, 能够降低土壤侵蚀, 促进土壤保持服务、固碳等服务能力的提升^[40—41]。2015年, 中共中央政治局召开会议, 审议通过《京津冀协同发展规划纲要》, 为京津冀地区注入了新的发展活力。在这4年间, 中国政府开始大规模推进京津冀城市群一体化发展战略, 开展一系列的基础设施建设活动, 改变了区域土地利用的方式, 造成了绿地面积的减少和生态空间的破碎化。这无疑加大对生态资源的消耗, 对当地的降雨、温度等自然环境因子产生了影响, 造成生态系统服务能力的退化^[42—43], 也使得研究区域生态系统服务总体则呈现下降的趋势。因此在未来该区域发展的过程中, 应将与人口密度、基础设施建设相关的政策作为必要的辅助行动, 来维护区域的生态系统健康。

3.3 生态资源与空间管控优化策略

社会与生态是相辅相成的, 二者只有相互协调达到良好的耦合状态, 才能实现研究区域生态保护和社会发展的双重目标。因此, 未来应该依据分区规划, 因地制宜的采取有针对性的生态资源与空间管控优化策略, 并积极推动生态补偿体系建设, 以实现区域协调发展^[44—45]。在京津保中心区, 通过疏解非必要的首都功能, 控制人口聚集、城市扩张来减弱对周边生态系统的影响, 并适当增加生态用地比重。在燕山和太行山区内应继续开展一系列的生态修复工程。在保护生物多样性的前提下, 可以适当促进文旅产业发展, 增加当地收入。在坝上高原区, 受到降雨等影响较大, 且生物多样性低, 应通过增加植被覆盖、保护河流, 增强该区域的水源涵养和地表蓄水能力, 提升区域生物多样性^[46]。在低平原区应关注土壤沙化问题, 合理开发水资源, 保护生态系统免受退化^[47]。

4 结论

由于城市化进程的加快, 人类活动极大地改变了城市群的生态系统服务能力。本研究选择以2000—2018年近20年的长时间跨度, 探索生态系统对研究区域的功能反馈效益。并探讨了生态-社会驱动因子与生态系统服务的关系, 揭示其在空间尺度上对生态系统服务的影响程度与作用机制。结果表明：1)2000—2018年, 研究区域的生态系统服务综合指数呈现出先上升后下降, 在2015年达到最高值。2)2000—2018年生态系统服务综合指数呈现明显的空间分异特征, 以燕山和太行山为界, 燕山和太行山区域呈现先减弱后提升的趋势, 西北区域则整体呈现减弱趋势, 东南区域从2015年起逐渐降低。3)自然-社会因子对生态系统服务的驱动作用从强到弱排列依次是：降雨量>土壤含沙>人口>温度>基础设施4)。在驱动因子中, 分析结果表明：人口密度、土壤含沙量、基础设施密度呈较为明显的负相关驱动特征, 降雨量呈较为明显的正相关驱动特征。温度呈南部负相关, 北部正相关。同时, 本文还依据实验结果提出了相应的优化策略, 为研究区域实现生态系统可持续管控提供有价值的政策指引。