呼包鄂榆城市群生态系统服务价值驱动因素及其交互效应

乐荣武,李 巍, *,周思杨,宋南奇

1 北京师范大学环境学院环境模拟与污染控制国家重点联合实验室,北京 100875

2 交通运输部水运科学研究院,北京 100088

生态系统服务(ES)被广泛定义为直接或间接促进可持续人类福祉的生态特征、功能或过程[1—2]。生态系统服务价值(ESV)是人类直接或间接地从生态系统功能中获得的有形或无形的利益的客观体现, 是区域生态文明与可持续发展的重要表征[3]。ESV 的空间格局和演变受到自然和社会经济等多方面驱动因素的影响[4],而且各驱动因素通过直接或间接途径产生复杂的相互作用关系。因此,识别ESV多个驱动因素的交互效应,对于区域人类活动协同管控和系统提升区域生态系统服务功能具有重要意义[5]。

当前,对ESV的研究逐渐从ESV的评价和分析其时空变化转向到分析驱动机制。要摸清ESV的驱动机制,必须识别ESV的驱动因素并探究其交互作用关系[6—7]。探究ESV驱动因素常用的方法包括冗余分析[8—9]、相关性分析[10]、主成分分析[11—12]、逐步回归[13]和Logistic回归[14]等传统统计分析方法,这些方法仅能识别驱动因素的作用大小和方向,无法识别多个驱动因素之间的交互作用。有较多的研究使用地理探测器模型探究ESV的主导驱动因素及各驱动因素的交互作用[15—18],虽然地理探测器可以对ESV驱动因素进行非线性归因且能识别驱动因素之间的交互作用,但不能判断驱动因素的作用方向,也无法刻画ESV对各驱动因素的变化响应特征[19—20]。地理探测器对驱动因素的数据要求转化为类型数据,该处理方法涉及可变面元问题增加了研究结果的不确定性[21—22]。在交互路径识别方面,已有学者使用结构方程模型对驱动因素间的作用路径进行定量测度[23—27],但这类研究仍然较少,有待进一步探究。总之,当前的研究方法尚未系统全面分析ESV驱动因素的重要性、作用方向、交互作用和刻画ESV对驱动因素的变化响应特征。如何系统性揭示ESV变化的驱动机制,是当前亟需解决的问题。近年来,机器学习模型因其可以很好地表征驱动因素与ESV之间的非线性关系而被大量用于ESV的驱动分析[28—30]。然而,机器学习模型的大多数算法都被视为黑箱模型,并不能充分解释驱动作用机制。沙普利加性解释(SHAP)方法是近期开发的机器学习模型解释工具,不仅可以衡量驱动因素的贡献大小,还能识别驱动因素的影响方向及交互效应[31—32],但将其用于ESV驱动作用机制的研究尚未见到。因此,将SHAP方法用于探索ESV驱动因素的交互作用机制方面具有较大的应用潜力。

呼包鄂榆城市群作为我国典型的快速城市化地区,是黄河流域人口、生产力布局的主要载体。城市群地处农牧交错带,生态本底脆弱,资源开采、城市用地扩张等导致生态用地破碎化,生态系统极易受到破坏。本研究重点关注 ESV 变化与土地利用类型转移模式之间的关系,探索ESV驱动因素交互作用的潜在机制。具体研究内容为:(1)分析研究期间土地利用和ESV的动态演变;(2) 耦合随机森林模型与SHAP 方法分析ESV驱动因素的重要程度、变化响应特征及其交互效应;(3)采用偏最小二乘路径模型识别驱动因素的交互路径。本研究系统性探究了呼包鄂榆城市群ESV驱动因素的交互作用机制,以期为黄河流域生态保护和高质量发展提供科学参考。

1 研究区概况与数据来源

1.1 研究区概况

呼包鄂榆城市群是以内蒙古自治区呼和浩特市、包头市、鄂尔多斯市和陕西省榆林市组成的国家级城市群[33]。该地区国土面积约17.5万km2,位于我国中西部的干旱半干旱区,土地利用以草地为主,约占总面积的54%(图1)。2020年总人口约1197万人,平均海拔约1300 m,年均气温约8℃,年均降水约320 mm。该区域是我国重要的能源、煤化工基地、农牧产品加工业基地和稀土新材料产业基地,也是黄河流域经济发展的重要增长极。

图1 研究区土地利用类型分布

1.2 数据来源

研究使用的主要数据和来源包括:①土地利用类型数据、行政边界矢量数据、气温、降水、潜在蒸散发、国内生产总值(GDP)格网和人口密度格网数据来源于中国科学院资源与环境科学数据中心(https://www.resdc.cn/),空间分辨率均为1 km,将土地利用类型重分类分为6个一级类, 分别为耕地、林地、草地、水域、建设用地、未利用地;②高程数据从地理空间数据云(https://www.gscloud/)获取,空间分辨率为90 m,在ArcGIS中使用分析工具计算坡度;③植被净初级生产力(NPP)数据从美国国家航空航天局(https://www.nasa.gov/)获取, 空间分辨率为500 m;④土壤属性数据从国家地球系统科学数据共享服务平台(http://www.geodata.cn/)获取, 空间分辨率为1 km,土壤保持量数据通过通用土壤流失方程(ULSE)进行计算;⑤其他社会经济数据来自《全国农产品成本收益资料汇编》《呼和浩特市统计年鉴》《包头市统计年鉴》《鄂尔多斯市统计年鉴》《榆林市统计年鉴》等。将所有空间数据统一到同一坐标系,并重采样至1 km。本研究的ESV时空分析部分均在1 km的栅格尺度上进行,综合权衡数据分析精度和模型计算效率,以5 km×5 km 的格网为基本单元,将研究区划分为6817个格网,在每个格网内统计所有变量的均值,用于本研究的驱动因素分析部分。

2 研究方法

2.1 生态系统服务价值核算

基于谢高地等的研究[34],1个标准单位ESV当量因子是指1 hm2当年平均产量的农田自然粮食产量经济价值的1/7。粮食产量价值主要依据稻谷、小麦和玉米进行计算,本文以研究区2015年的物价水平为基准,计算出研究区主要农作物平均单产价值,最终得到研究区1个标准当量经修正后的价值量为1463.69元/hm2(表1)。

表1 不同土地利用类型单位面积价值系数/(元/hm2)

由于研究区不同年份、不同区域的生态系统的基本情况是变化的,ESV也相应发生动态变化。本文参考谢高地等的研究[34],认为生态系统食物生产、原材料生产、气体调节、气候调节、净化环境、维持养分循环、生物多样性和美学景观功能与生物量在总体上呈正相关,水资源供给和水文调节与降水变化相关,而土壤保持与降水、地形坡度、土壤性质和植被盖度密切相关。故选取NPP、降水量和土壤保持量三项因子对当量进行动态调节,以此构建生态服务时空动态价值当量表,计算公式为:

(1)

式中,Fnij指某种生态系统在第i年第j地区第n类生态服务功能的单位面积价值当量因子;Fn指该类生态系统的第n种生态服务价值当量因子;n1表示与NPP相关的服务功能;n2表示与降水相关的服务功能;n3指土壤保持服务功能,Pij指NPP时空调节系数,Rij指降水时空调节系数,Sij指土壤保持时空调节系数。

研究区ESV计算公式为:

(2)

(3)

式中,c为第c种生态系统服务功能;Ec为第c种生态系统服务功能价值;Fnij表示某种生态系统在第i年第j地区第n类生态服务功能的单位面积价值当量因子;D为多年的1个标准当量因子的生态系统服务平均价值量(元/hm2),此处为1463.69元/hm2;Aij为第i年第j地区的面积。

2.2 生态系统服务价值变化趋势分析

运用一元线性回归方程拟合1990—2020 年ESV的时空演变特征,IESV表征每个栅格ESV线性变化斜率值,若IESV>0,则表明ESV有改善趋势; 反之则表示存在退化趋势。采用F检验ESV变化趋势的显著性。结合IESV及F检验结果,ESV变化趋势可分为5类:大幅改善(IESV>0,P<0.01),一般改善(IESV>0,P<0.05),变化不显著(P≥0.05),一般退化(IESV<0,P<0.05)和大幅退化(IESV<0,P<0.01)[35]。

2.3 随机森林模型与SHAP方法

本研究使用机器学习中的随机森林模型分析ESV的驱动因素[36],耦合SHAP方法对随机森林模型进行解释。随机森林相对于传统统计模型在进行预测时往往有更好的精度,但是同时也失去了统计模型的可解释性,所以随机森林通常被认为是黑箱模型。针对基于树的机器学习模型,2017年Lundberg和Lee提出了SHAP方法用来解释各种黑箱模型[37]。与以往研究中多采用线性拟合相关的方法或地理探测器模型量化ESV的驱动因素且无法排除驱动因素间的相互干扰相比,耦合随机森林模型与SHAP 方法不仅可以识别ESV和驱动因素的非线性关系,还能够分离出每个驱动因素对ESV的独立影响以及各个因素间的交互影响。

将GDP密度、人口密度、降水、气温、潜在蒸散发、高程、坡度、土壤容重、黏土含量9种空间变量(图2)以及6种土地类型面积比例共15个驱动因素作为随机森林模型的输入变量,将ESV作为模型的输出变量,模型样本量为6817个。

图2 驱动因素的空间分布

2.4 偏最小二乘路径模型

偏最小二乘路径模型(PLSPM)是一种分析多变量间复杂因果关系的综合分析模型,属于结构方程模型的一种方法。该模型不仅能解决变量间存在的共线性问题,对变量的分布状态无要求,还可以计算不同变量对响应变量的直接效应和间接效应[38]。模型中路径系数反映了变量之间关系的方向和强度, 而预测变量和响变量之间相乘的路径系数则显示出间接影响的强度。拟合度指数(GOF)用于评估模型的预测性能, 数值越大, 模型预测效果越好。本研究使用偏最小二乘路径模型识别气象、地形、土壤、土地利用和社会经济因素对ESV的交互影响作用路径。

3 结果与分析

3.1 土地利用类型演变

在1990—2020年间,耕地面积整体呈现下降趋势,在1995年面积达到峰值后逐渐下降,损失的耕地主要转化为草地(33.95%)和林地(5.12%)。林地面积在1995年达到最低值后逐渐上升,林地增加的主要原因是来自草地(35.03%)和耕地(19.21%)的转化。草地面积变化剧烈,草地与耕地和未利用地存在较大比例地互相转化;水域呈现先减后增的趋势;建设用地持续增加,主要来自于草地(40.11%)和耕地(28.89%)的转化(图3)。

图3 1990—2020年研究区土地利用转移变化

3.2 ESV的时空演变

3.2.1ESV的时间变化

1990—2020年研究区ESV呈先下降后上升趋势,其中在1990—1995年ESV下降的最多,在此期间林地面积和水域面积减少幅度均为最大,分别减少了7.64%、14.14%,损失的林地主要转变成耕地和草地,损失的水域主要转变成草地。在2000—2005年ESV上升的最多,在此期间林地面积增幅最大,达到了11.89%,增长的林地主要来源于耕地和草地。ESV在研究期内总体上增加了62.28亿元,主要在于由草地和耕地转化而来的林地面积增加了21.18%。(表2和图4)。

表2 1990—2020年研究区各地类ESV/亿元

3.2.2ESV的空间变化

在空间格局上,研究区ESV空间分异明显(图5)。研究区ESV总体呈现东高西低的分布格局,ESV高值分布以河流水系为中心,ESV较高值集中分布在呼和浩特市大青山保护区的林地,低值以未利用地类型为主要分布区。该分布趋势与土地利用分布基本吻合,研究区西北部干旱少雨,人类活动较少,土地利用以未利用地为主,植被覆盖率低,ESV较低。而东南部地势相对平坦,海拔较低降水较多,植被覆盖率高,ESV较高;河流成为高值中心主要与水域的单位面积ESV高有关。

图5 1990—2020年研究区ESV密度分布

总体来看,在1990—2020年间73.52%研究区的ESV变化不显著,ESV的退化区面积占比(9.76%)高于改善区(6.38%)(图6)。大幅退化和一般退化面积占比分别为1.62%、8.14%,大幅改善和一般改善面积占比分别为1.47%、4.91%。退化区主要集中分布在包头市北部和零星分布在鄂尔多斯市大部分地区,改善区主要分布在呼和浩特市和榆林市的东南部地区。虽然ESV的退化面积大于改善面积,但由于改善区主要位于土地类型由未利用地、草地向林地、水域转化的区域,单位面积林地和水域的ESV更高,因此整体上研究区的ESV呈增加趋势。

本文进一步通过热点分析揭示ESV的空间集聚特征及其演变规律。研究区ESV 冷热点总体上呈“西冷东热”的空间分布格局(图7),研究区约84%地区的ESV空间集聚特征不显著,ESV高值集聚和低值集聚所占区域面积比例分别约为10%、6%。ESV 热点区集中分布在研究区东北部呼和浩特市北部的武川县境内,该区域林地分布广泛,ESV 次热点区围绕热点区四周分布。 ESV 冷点区主要分布在研究区西北部的鄂尔多斯市杭锦旗境内,ESV 次冷点区聚集分布在包头市的最北部境内,该区域分布较广的是未利用地。从各时期阶段来看,热点区和冷点区面积呈波动下降趋势,次热点区和次冷点区面积呈波动上升趋势,表明研究区ESV 高值和低值聚类均逐渐弱化。

图7 1990—2020年研究区ESV冷热点的空间分布

3.3 ESV驱动因素分析

3.3.1驱动因素的重要性识别

耦合随机森林模型和 SHAP方法分析研究区ESV的驱动因素作用方向及各驱动因素的相对重要性(图8)。ESV的驱动因素重要性从大到小依次为水域比例、未利用地比例、林地比例、降水、坡度、草地比例、土壤容重、人口密度、潜在蒸散发、GDP密度、气温、黏土含量、高程、建设用地比例和耕地比例,表明土地类型中水域比例是影响ESV最重要的驱动因素,其次是未利用地比例,耕地比例在研究区的重要性最低。从驱动因素的作用方向来看,水域比例、林地比例、降水、坡度、草地、土壤容重、黏土含量和高程对研究区ESV有促进作用,未利用地比例、人口密度、潜在蒸散发、GDP密度、气温、建设用地比例和耕地比例对研究区ESV有抑制作用。从驱动因素的维度来看,土地利用类型的贡献度为61.24%,成为最重要的驱动力;其次是地形和气象的贡献度分别为17.59%和17.05%,土壤和社会经济的贡献度分别为2.39%和1.73%。

图8 驱动因素重要性排序图和模型SHAP摘要图

3.3.2ESV对驱动因素的响应特征

通常,ESV驱动因素间存在不同强度的相关性,这会干扰单一驱动因素对ESV的影响分析,SHAP 方法可以较好地排除其他因素的干扰,剥离出ESV随单一驱动因素的变化趋势(图9)。

图9 ESV对驱动因素的变化响应特征

总体来看,ESV对驱动因素的响应呈现出非线性变化特征。在社会经济因素中,当GDP密度和人口密度较低时,GDP和人口的SHAP值大于0,显示正贡献,表现出对ESV的促进作用;随着经济的发展和人口的增长开始对ESV产生抑制作用。在气象因素中,降水和气温对ESV的影响均存在明显的阈值效应,降水在低于和高于250 mm时分别对ESV起抑制和促进作用,气温在低于和高于7.5℃时分别对ESV起促进和抑制作用;在潜在蒸散发低于1000 mm时,潜在蒸散发对ESV的促进作用随潜在蒸散发的增加持续减弱,当潜在蒸散发高于1000 mm时,响应曲线在SHAP为0值的附近维持较为平稳,表明对ESV的影响较小。高程在低于和高于1300 m时分别对ESV起抑制和促进作用,随着高程的增加促进作用不再增强,坡度对ESV始终起促进作用。土壤容重和黏土含量在低值时对ESV无显著作用,在高值时具有促进作用。耕地比例较低时对ESV无显著影响,较高时开始呈现抑制作用。未利用地比例在低于和高于25%时对ESV分别起促进和抑制作用,建设用地比例对ESV始终起抑制作用。林地和水域比例始终对ESV起促进作用,草地比例在低于和高于50%时分别起抑制和促进作用。

3.3.3驱动因素的交互效应

ESV的变化受驱动因素的综合影响,不同驱动因素之间还存在相互影响,SHAP 方法可捕捉交互性最强的成对的驱动因素对ESV的相互作用效果(图10)。

图10 驱动因素的交互效应

GDP密度和人口密度与其他驱动因素无明显交互作用,这可能是由于研究区的社会经济因素对ESV影响的贡献率较低。在以水域比例为交互项中,在水域比例较高的地区,降水与水域比例的交互作用对ESV呈现出由抑制作用减弱到促进作用增强的变化趋势,而气温、潜在蒸散发与水域比例的交互作用对ESV均呈现出由促进作用减弱到抑制作用增强的变化趋势。在水域比例较低的地区,降水与水域比例的交互作用对ESV呈现出由抑制作用减弱到促进作用增强的变化趋势,但交互作用程度相对较低;气温、潜在蒸散发与水域比例的交互作用均不明显,未利用地比例对ESV具有更强的抑制作用。在降水量较多的地区,高程与降水无明显交互作用;但在降水量较少的地区,高程与降水的交互作用对ESV呈现出由抑制作用减弱到促进作用增强后又趋于稳定的变化趋势。在以林地比例为交互项中,在林地比例较高的地区,坡度、土壤容重和黏土含量的增加均对ESV呈现抑制减弱到促进增强的变化趋势,耕地比例的增加对ESV具有抑制作用;在林地比例较低的地区,坡度、土壤容重、黏土含量和耕地比例与林地比例无明显交互作用。林地比例在坡度大的地区较在坡度小的地区对ESV有更强的促进作用。在以耕地比例为交互项中,在耕地比例较大(较小)的地区,草地比例的增加对ESV的抑制(促进)作用减弱(增强),建设用地比例在耕地比例较大的地区对ESV具有更强的抑制作用。水域比例在降水较多的地区对ESV具有更强的促进作用。

随机森林模型耦合SHAP方法只能衡量两个驱动因素间的交互作用大小、方向及变化趋势,多种驱动因素间的交互作用路径如何需要进一步探讨。将驱动因素分为社会经济(GDP、人口)、气象(降水、气温、潜在蒸散发)、地形(高程、坡度)、土壤(黏土含量、土壤容重)和土地利用类型共5个潜变量,采用偏最小二乘路径模型探究各驱动因素潜变量间的交互路径(图11),模型拟合度指数(GOF)为0.64,表明模拟结果符合精度要求。

图11 驱动因素的交互作用路径和影响程度

土地利用类型对ESV的影响系数为0.696(图11),再次表明土地利用对ESV的影响起着决定作用。社会经济受到地形因素的制约影响,社会经济对土地利用和ESV的影响系数分别为0.222和-0.186,表明社会经济对ESV的直接影响是负效应,社会经济主要通过直接影响土地利用进而间接影响ESV。气象对ESV的直接作用较弱,主要是通过气象-土壤-土地利用或气象-土地利用作用路径影响 ESV。地形对ESV的影响路径有多条,地形对社会经济、气象、土壤和土地利用均有直接影响,其中对土壤的直接影响力最强。土壤受到地形、社会经济和气象的影响,土壤对ESV的影响主要是通过影响土地利用产生的间接影响。除土地利用以外,其他因素对ESV的间接影响均大于直接影响。

4 讨论

在ESV的变化响应特征方面,以前的研究尚未有ESV对驱动因素变化响应特征的探讨;在交互作用分析方面,当前大量的研究采用地理探测器方法,但该方法未能识别交互作用的方向性;耦合随机森林模型和SHAP的方法能识别驱动因素在不同值域内的交互作用大小、方向和非线性响应特征;在交互路径分析方面,当前研究主要采用结构方程模型识别交互路径,本文选择结构方程模型中的偏最小二乘路径模型,该方法对样本量和样本分布要求低[23,38],模拟结果精度较好,表现出了较好的适用性。综上,本研究耦合随机森林模型和SHAP方法并结合偏最小二乘路径模型,可以较为系统全面地分析ESV的驱动因素及其交互效应。

本研究表明土地利用的变化强烈影响 ESV,这与之前的研究一致[23,39]。ESV在2000年前后分别呈下降和上升的趋势,ESV的增加主要来自于耕地和草地向林地转换的结果。林地面积在2000年之后逐渐增加,这得益于当地“退耕还林还草”政策和“三北”防护林工程的实施。尤其是在ESV大幅改善的区域,例如,位于榆林北部的毛乌素沙地,由于近二十年的持续造林治沙活动,ESV得到了大幅提升。人类社会经济活动对生态系统服务的需求是其影响ESV变化的根本原因[24],本研究中GDP和人口的增长对ESV的影响由促进作用转变为抑制作用,表明较低程度的人类活动对生态环境具有改善作用,但在高收入地区大量人口的高需求和有限的供应水平将导致生态系统服务稀缺,从而对ESV产生抑制作用。偏最小二乘路径模型进一步揭示了社会经济因素对ESV的直接和间接影响是相反的,社会经济通过气象、土壤和土地利用对 ESV 产生间接的积极影响(图11),表明经济发展可能会在自然变化过程中改变生态系统服务演化的方向[40]。以往的研究表明,城市扩张和经济发展会直接加剧生态系统服务供需失衡[41]。此外,更多生态环境政策的实施也会间接促进了生态系统服务功能的恢复和提升[42]。因此,结合经济发展和自然条件的空间差异,优化土地利用配置,对提高ESV具有重要意义。考虑到水域是最重要的驱动因素,草地对研究区的ESV贡献最大。未来的生态恢复实践应优先考虑湿地和湖泊生态系统提供的生态价值和效益,通过跨省联动、多方补水、关闭周围煤矿等保护措施,加强红碱淖湿地国家自然保护区的生态保护;在沿黄河湿地建设沿黄生态廊道,严格保护基本林地和草原。

本研究发现降水和气温对ESV的影响具有阈值效应,在年降水量低于250 mm或年均气温高于7.5℃的地区,降水或气温成为该地区ESV提升的限制因子,即该地区降水和气温对ESV具有负贡献,随着降水的增加或气温的降低这种负贡献逐渐减小直至变为正贡献。在整个研究区内,降水的增加或气温的降低始终对ESV具有提升作用,原因在于研究区位于干旱半干旱区,热量充足而降水量小,有限的降水还未被植被充分吸收利用就已经蒸散发了。在交互分析中,水域比例是参与最多的交互项,主要在于ESV当量表中水域的值最大,水域比例成为最重要的驱动因素,这也导致降水和气温在水域比例较高的地区对ESV具有更大的影响程度。

本研究仍存在一定的局限性。首先,基于当量因子法的ESV的评价严重依赖于土地利用数据,在一定程度上可能夸大了土地利用类型对ESV 的影响,后期应尝试采用其他ESV 评价方法进一步验证不同驱动因素的重要性。其次,本研究ESV的驱动因素分析是基于分辨率为5 km的格网尺度,但是不同空间尺度上的驱动影响和交互关系尚不清楚,尺度效应增加了研究结论的不确定性,后续应该在不同空间尺度探究ESV的驱动因素。此外,本研究尚未考虑政策因素对ESV的影响,因此未来的工作应该将政策因素充分考虑进来。

5 结论

本研究核算了1990—2020年呼包鄂榆城市群的ESV并分析其时空变化特征,耦合随机森林模型和SHAP方法探究了ESV驱动因素的贡献程度、作用方向、变化响应特征及其交互效应,进一步采用偏最小二乘路径模型识别了交互路径。结论如下:(1)草地对研究区的ESV贡献最大,其次是水域;研究区的ESV在1990—2020年呈波动上升趋势,总体增加了62.28亿元,ESV的增长主要源自于单位面积ESV较低的草地和耕地向单位面积ESV较高的林地的转化。(2)ESV对驱动因素呈非线性变化响应特征;土地利用类型是最重要的驱动因素,其贡献度达到61.24%,地形和气象的贡献度分别为17.59%和17.05%,土壤和社会经济的贡献度较低,分别为2.39%和1.73%。(3)水域比例是最重要的交互项,其次是林地比例;不同因素间交互作用在因素处于不同范围内表现出不同的交互效应。(4)土地利用类型直接作用于ESV,其他因素主要通过直接影响土地利用类型进而间接影响ESV,且对ESV的间接影响均大于直接影响。