生态系统服务表示生态系统结构、过程和功能通过直接或间接途径提供给人类的各种惠益^[1], 自20世纪90年代以来生态学、经济学和地理学学者针对不同时期、不同空间尺度、不同类型的生态系统服务开展了大量评估工作^[2-6], 为生态政策制定、生态系统决策管理提供了重要依据^[7-8]。研究方法上主要基于单位服务功能价格的方法和基于单位面积价值当量因子的方法^[9], 相较而言, 后者对中国生态系统服务价值评估具有更强的针对性^[10], 研究特征上呈现出跨学科属性的交叉融合、静态分析向动态时空研究深入拓展, 单一系统向多系统维度比较演进, 区域生态系统服务空间差异性和不同地理尺度的关联特性受到进一步重视^[11]。

提高人类福祉是可持续发展最终目的^[12], 也是生态系统服务研究的根本出发点^[13]。2005年千年生态系统评估报告明确提出生态系统服务对人类福祉贡献的概念框架, 这是福祉研究发展中的重要里程碑^[14], 至此, 生态系统服务研究如何延伸至人类福祉提升等应用指向是重点关注方面^[15], 二者相互关系成为可持续科学的核心议题^[12]。福祉由多种要素构成, 其含义与幸福、福利或利益相接近^[16], 分为主观和客观两个方面, 其中由联合国开发计划署基于能力评估框架提出的人类发展指数(HDI)影响最为广泛和深远, 作为人类福祉的有效衡量指标受到普遍认可^[17]。生态系统通过供给、调节、支持、文化四大功能为人类维持基本物质生存需求和实现高品质精神生活提供几乎所有的福祉要素^[18-19], 研究普遍发现生态系统服务与福祉间存在显著的相关性, 其中, 人类活动是其联系的中介^[20], 但二者之间并非呈现出简单线性关系, 随着服务类型^[13]、空间尺度^[21]、时间范围^[19]、群体差异^[22]等表现出不同形式的反馈。

关于生态系统服务与居民福祉关系展开的研究从定性向定量转变, 集中于相关关系、耦合协调、时空格局研判与分类、空间权衡/协同、影响与溢出效应等方面^[23-25], 并将生态系统服务价值评估纳入各级决策管理^[26]；研究手段上主要采用spearman等级相关系数、双变量空间自相关分析、结构方程模型、Logistic回归、灰色关联模型等^{[25, 27-29]}；地域上以全球、全国、省域等宏观尺度单元入手^{[14, 30-31]}, 或以不同地貌类型及流域为空间范围^{[21, 31-33]}, 基于耕地、草原、湿地等特定的生态服务类型, 从当地居民视角探讨生态系统的福祉效应较多^{[20, 25]}, 以县域为研究尺度将生态系统服务与居民福祉二者变化过程纳入分析视野的研究较少。山区作为我国主要地貌空间单元, 绿色生态基底好但生态脆弱敏感, 是国家经济发展重要的资源支撑地和生态系统服务供应地, 同时也是社会经济发展洼地^[34], 人与自然协调问题突出, 且内部生态、经济地域空间差异较大, 现有研究对山区不同生态服务类型的福祉影响方向和程度的时空变化具体刻画不足。

基于此, 本文以西部典型山区——川西为研究案例, 通过空间计量分析手段探讨2000年《长江上游、黄河上中游地区天然林资源保护工程实施方案》批准以来, 川西山区居民福祉在生态环境保护与生态价值提升过程中是否受到积极反馈效应, 不同生态系统服务类型的福祉影响有何差异, 其在空间上影响方向和强度有何不同。该研究在视角上将生态系统服务作为居民福祉的关键影响要素, 并从时空异质性角度出发把握生态系统服务对福祉影响的规律性特征；在研究手段上融入GWR计量分析手段, 打破了对研究区整体或平均的拟合方法, 有助于进一步把握影响结果的局域特征, 可为地区制定科学合理的生态系统保育政策提供一定依据, 助推西部山区生态系统服务和居民福祉的协调可持续发展。

1 研究区概况与数据来源

川西山区地处川、甘、青、藏、滇五省结合部, 是青藏高原东南缘向长江中下游平原过渡地带。行政区范围涵盖阿坝藏族羌族自治州、凉山彝族自治州、甘孜藏族自治州、攀枝花市, 辖53个区县, 介于97°20′—104°25′ E, 26°02′—34°18′ N之间, 幅员面积30.03×10⁴km², 占四川省行政区域面积61.79%(图 1)。

川西山区是世界生物多样性热点区和我国生物多样性保护关键区, 也是我国三大林区之一青藏高原东部林区、五大草原之一川西北大草原的重要组成部分。作为我国青藏高原、长江上游水源涵养地和生态屏障示范区, 生态系统服务变化对长江中下游地区具有较为深远的影响, 该地区素有“蜀地之肾”之称, 境内林地面积约1300万hm², 占全省林地面积的54.2%, 森林覆盖率达43.36%。此外, 川西山区是我国南方农耕文化和游牧文化交汇带, 是汉族聚居区向多民族聚居区的过渡带^[35], 是我国“藏羌彝走廊”的核心区之一^[36], 区域地理位置偏安一隅, 生态环境脆弱, 农户生计主要来源于传统农牧经济, 最基本食物和能源等生活必需品在很大程度上依赖于自然生态系统服务供给, 表现出自然资本的强依赖性。近年来在气候变化、灾害过程、人类干预的外部环境扰动下, 流域生态—经济系统受到较大影响^[24], 生态环境保护与经济社会协调发展面临困境。

研究数据主要包括川西山区县域土地利用数据以及社会经济统计数据。2000、2005、2010和2015年土地利用类型数据来自中科院资源环境科学数据中心土地利用数据库(http://www.resdc.cn), 该数据以Landsat—TM/ETM/8遥感影像为数据源, 采用人机交互式目视判读方式构建, 分辨率为30 m, 精度达到90%以上^{[5, 37]}。社会经济数据主要来自阿坝州、凉山州、甘孜州、攀枝花市相应年份统计年鉴, 其中, 预期寿命、文盲率和平均受教育年限数据主要来源于第五次和第六次四川省人口普查资料, 其余年份通过内插法及外推法计算得到。

2 研究方法 2.1 生态系统服务价值评估模型

本研究根据土地利用类型将生态系统划分为耕地、林地、草地、水域、未利用土地以及建设用地六大类, 借鉴千年生态系统评估方法, 将生态系统服务分为供给、调节、支持和文化四大服务子类。以2015年最新的“中国生态系统单位面积生态服务价值当量”作为依据^[9], 确定川西山区生态系统服务价值。根据已有研究^[21], 令当年平均粮食单产市场价值的1 /7表示为1个生态服务价值当量因子的经济价值量；查阅同时期川西山区平均粮食产量情况以及全国农产品成本收益资料汇编中粮食平均出售价格, 计算得到2000、2005、2010和2015年4年农田自然粮食产量的经济价值；考虑到货币购买力随时间发生变化, 以2000年为基期通过购买力指数进行平减(表 1), 模型表达如下：

(1)

式中, ESV_i为i县t时间生态系统服务价值, A_ijt为i县t时间第j种土地利用类型面积, EV_ij是i县j种土地利用类型单位面积生态系统服务价值当量, VS_t为t时间1个生态服务价值当量因子的经济价值量, PI_t为t时间购买力指数, m为53个县域单元, n为6种土地利用类型, t为时间。

2.2 居民福祉的界定与测度

研究中参考2010年以后人类发展指数的计算方法, 从健康长寿的生活、知识以及体面的生活水平3个维度平均成就度量居民福祉水平^[38], 计算公式如下：

(2)

(3)

(4)

(5)

式中, HDI表示居民福祉, H、E、I分别代表健康指数、教育指数和收入指数。H_e、L_e、Y_e、I_e分别为i县平均预期寿命, 文盲率、平均受教育年限, 按照购买力平价折算的实际人均GDP, max和min表示4个年份中该指标出现的最大值和最小值。

2.3 空间自相关检验

空间自相关检验用于检验某种现象在空间上是否存在集聚, 是构建地理加权回归模型的基础和保障。本文运用Moran′s I指数对自变量与因变量进行空间相关性检验, 公式如下：

(6)

式中, n为县域个数；y_i和y_j分别为县域i和县域j的居民福祉指数；y为全体样本居民福祉指数均值；W_ij为二进制邻接空间权重矩阵。Moran′s I的值一般在[-1, 1]之间, 大于0表示区域间存在正自相关性, 小于0表示存在负自相关性, 等于0表示空间单元属性呈随机分布。

2.4 生态系统服务价值的福祉驱动效应模型

地理加权回归(Geographically Weighted Regression, GWR)是在普通线性回归模型基础上的拓展, 作为一种局域空间分析方法, 其采用局部加权最小二乘法进行逐点估计, 通过获取局部参数判断自变量与因变量关系在空间尺度的变异, 可直观反应研究变量之间的空间依赖与分异规律^[39]。其模型设定形式为：

(7)

式中, (u_i, v_i)为第i县空间地理位置, β₀(u_i, v_i)是第i县回归常数项, β_k(u_i, v_i)为第i县第k个回归参数值, ε_i表示随机误差项。由于人口的变化对区域土地利用格局会产生重要影响, 生态系统服务水平据此产生变化, 生态系统服务价值人均量对衡量区域人地关系状态具有重要意义^[40], 因此本研究自变量x_ik分别采用第i县人均供给服务价值、人均调节服务价值、人均支持服务价值、人均文化服务价值、人均生态系统服务价值变量, 因变量y_i为第i县福祉水平。地理加权回归计算在GWR 4.0中进行, 通过ArcGIS 10.2进一步作出空间可视化表达。

3 结果与分析 3.1 居民福祉指数空间自相关性

因变量存在空间自相关性是进行地理加权回归模型应用分析的前提条件, 可为模型的构建奠定基础。对居民福祉指数进行全域空间自相关检验, 4个研究年份中Moran′s I指数值均>0, 且Z-Score均明显大于检验临界值2.58, 表明在0.01显著水平下具有统计意义(表 2), 充分表明川西山区居民福祉水平在地域空间上具有显著的正向空间集聚性, 即福祉水平较高(较低)的县域趋向于和类似水平的县域集聚, 形成空间依赖关系。但从时间变化来看, 福祉水平地域空间集聚态势逐步弱化, Moran′s I指数值由2000年0.360减小为2015年0.297。

为进一步说明53个县域单元居民福祉水平空间上的相互关联类型和局域空间分布状况, 对居民福祉指数进行局域空间自相关分析, 从图 2可知：(1)川西山区居民福祉呈现同质集聚特征, 其中“高-高”集聚型数量稳定, 占比维持在11.32%, “低-低”集聚型趋于减少, 占比由2000年11.32%降低到2015年7.54%；(2)“高-高”和“低-低”集聚具有南北分异特征。居民福祉高水平县在2000年主要集中在攀西地区, 2010年开始“热点”区逐渐向岷江上游转移并进一步扩大, 攀西高水平集聚范围收缩；相反, 受到与落后县域为邻的不利影响, 福祉低水平县2000年在川西西北部呈现连片集中分布格局, 形成西部边缘“洼地”, 随着国家倾斜性政策的支持和基础设施的大规模建设, 居民福祉水平整体跃升, “低-低”集聚从西部边缘向东南部转移收缩, 聚集于凉山州东北部, 该地区为四川决战脱贫攻坚的“硬骨头”和地质灾害高发区。

总体而言, 全域和局域空间自相关结果均证实研究区居民福祉呈现出较强的空间依赖性, 支持回归模型对空间效应的引入, 因此运用GWR模型对生态系统服务的居民福祉效应进行分析具有必要性和可行性。

3.2 GWR模型的适用性

为避免多重共线性, 以居民福祉指数为因变量, 分别以人均生态系统服务价值及其构成的供给服务价值、调节服务价值、支持服务价值、文化服务价值进行单变量GWR建模。通过调整R²和赤池检验(AICc)两个参数进行模型适用性评估, 统计结果显示(表 3), 各自变量因子的OLS模型调整后的R²在0.2附近, 而GWR模型对因变量的解释力达到50%—60%, 较OLS模型提高了近40%；AICc值作为模型拟合显著性的标识, 由OLS模型的-77.582—-85.495下降到GWR模型的-114.433—-105.549, 减少了33.85%—36.05%, 远大于3^[39], 表明GWR模型改进回归的拟合度十分显著, 即使增加了计算复杂性, GWR模型也优于OLS模型。同时, GWR模型统计检验相应概率P值小于给定的显著性水平0.05, 具有统计学意义, 意味着居民福祉与生态系统服务价值各因素之间的整体关系表现出明显的空间不稳定性和非均质性特征。

3.3 生态系统服务对居民福祉影响的时空分析

生态系统服务因子显著负向影响居民福祉。GWR回归系数均值显示生态系统服务因子对居民福祉具有显著负向影响(图 3), 即生态系统服务水平越高的县域居民福祉发展水平反而越低, 表现出“资源魔咒”^{[37, 41]}, 研究结果与刘家根^[13]、任婷婷^[42]等结论具有一致性, 这与地区农牧业生产方式和生态经济发展程度密切相关, 同时也与生态系统服务变化对福祉影响具有一定的时间滞后性密不可分, 此外生态系统服务具有外部性特征, 服务与收益在空间上常常不一致^{[13, 43]}；不同生态系统服务类型中, 2000和2015年文化服务是福祉提升的重要制约变量, 2005和2010年供给服务对福祉的影响程度更大, 调节服务在4个年份中对福祉的制约作用均最小, 这可能由于生态系统服务价值具有尺度效应, 调节服务主要体现在全球和国家尺度^[44]；与2000年相比, 2015年各生态系统服务因子的抑制作用趋于缓解, 其中文化服务、供给服务的制约性大幅减小, 但依然是关键影响变量。原因可能在于供给服务与当地农户收入、能源获取密切相关, 但限制开发的政策约束了居民对大部分自然资源的开发利用^[19], 且地区对消遣娱乐、精神收益等文化资源的挖掘还较为薄弱。

局域空间上生态系统服务的边际效应在川西山区存在方向异性, 呈现出非均质性特征。影响因子在53个县域单元中正负波动(图 3), 即在不同县域生态系统服务水平对居民福祉会产生正向促进或逆向抑制的不同作用。对大部分县域而言生态系统服务价值是居民福祉提升的抑制力因子, 但受负向影响的县域数量占比随时间变化呈缩小趋势(图 4), 从2000年90.57%降低至2018年60.38%；不同生态服务类型中(图 4), 支持服务和文化服务的负向影响县域数量更多, 2000年占92.45%, 但随时间变化其中52.83%的县域影响系数由负转正；从系数波动来看, 2000和2015年文化服务弹性系数值空间波动最大(图 3), 2005和2015年供给服务系数值绝对差异最大, 随时间变化, 除文化服务弹性系数波动值先缩小后扩大外, 其余生态系统服务因子系数正负波动均缩小, 波幅减少约60%, 意味着生态系统服务对居民福祉影响的空间差异趋于缩小。

生态系统服务因子的福祉影响具有东西空间分异和集聚特征, 随时间变化呈现区域稳固性。人均ESV及不同生态系统服务类型对福祉的边际效应表现出东部正向驱动, 西部负向阻滞的特点(图 4)。具体而言, 大于0的影响系数以大渡河、岷江和安宁河流域上游为主, 2000年主要分布在阿坝州中部和凉山州东北部, 地区居民福祉受到生态系统服务因子的增量效应, 随时间变化进一步向阿坝州全域和凉山州东北部团聚状扩大, 形成两片正向影响“V”型区域；小于0的影响系数连片分布在川西西部, 生态系统服务变量对居民福祉均为减量效应, 表现出高生态系统服务价值-低居民福祉的权衡关系。这主要是因为川西西部地域面积广阔, 林草地资源丰富, 但山高谷深, 农业发展空间小, 农户生计方式单一, 生产生活方式封闭, 传统产业占比过高, 高等级交通运输网尚未全域覆盖, 生态资源红利释放不足；东部生态经济区率先建设, 拥有若尔盖草原、九寨沟、黄龙、松潘、马湖等重要旅游景点, 地理区位上更靠近平原经济区, 独特山地旅游景观对游客吸引力较强, 从而提高居民收入, 生态资源红利实现释放。

居民福祉对生态系统服务因子敏感性具有从东南向西北递减的特征, 负向敏感性高于正向敏感性。正向高敏感区位于凉山州东北部甘洛县、美姑县, 低敏感区位于阿坝州西北部九寨沟县、若尔盖县、松潘县(图 4), 2000—2018年人均ESV和支持服务正向影响区敏感程度递增, 其余服务类型弹性系数波动增大, 意味着地区居民福祉水平受到生态系统服务的正向带动效应整体增强。负向影响高敏感区集中于攀西南部, 低敏感区位于西北部白玉县、得荣县、巴塘县, 这主要是由于攀西南部地处长江干热干旱河谷退化区、金沙江下游水土流失退化区以及川滇黔喀斯特石漠化退化区, 工矿产业和地区农牧业过度开发引致水土流失问题严峻^[45], 脆弱的生态环境对地区农业增收和康养旅游的发展形成不利影响。得益于金沙江流域山水林田湖草生态保护修复工程实施, 与2000年相比2015年生态系统服务因子对福祉的减量效应减缓。

供给服务与文化服务是川西山区居民福祉的主导生态因子。通过比较影响系数绝对值, 确定影响程度最高的生态系统服务指标, 识别出居民福祉的主导生态系统服务因子及其作用区域(图 5)。其中2000年和2015年分别有75.47%和83.02%的县域居民福祉受到文化服务的影响最为强烈, 2005年和2010年供给服务对福祉的作用面积最大, 占区域总面积的81.13%和83.02%。供给服务和文化服务随时间变化主控区域具有空间交替性, 结合图 4可知, 供给服务在2000年对川西东北部、2005年和2010年对川西西部居民福祉提升形成关键的抑制力, 2015年对川西东北部居民福祉改善提供了重要驱动力, 这主要是由于阿坝州逐步探索出生态产品价值实现路径, 打响地区生态产品品牌, 如汶川甜樱桃、茂县花椒、小金苹果等；结合图 4可知, 2000和2015年川西西部县域主要受控于文化服务的抑制作用, 2005和2010年文化服务对川西东北部居民福祉的提升带来的驱动作用最大, 受益于旅游文化产业发展, 旅游业对阿坝州经济带动性强, 旅游精品目标逐步明确。

4 结论

本研究以川西山区为例, 基于人类发展指数和生态系统服务价值当量评估法, 利用空间自相关检验和GWR空间回归模型对生态系统服务因子的福祉效应时空变化进行探测, 并识别出不同空间关键生态服务因子, 所取得的主要结论如下：

(1) 川西山区居民福祉水平在地域空间上具有明显的同质集聚特征, 空间依赖关系随时间变化逐步减弱。两极在南北分异明显, “高—高”集聚区汇聚在攀西南部并逐渐向岷江上游转移扩大, “低—低”集聚区分布在川西西部边缘并向东南部转移收缩。

(2) 生态系统服务因子显著负向影响居民福祉, 但阻滞效应随时间变化逐渐缓解。边际效应在局域空间具有不稳定性和非均质性特征, 生态系统服务因子系数值均在正负间波动, 系数空间差异随时间变化趋于缩小。

(3) 生态系统服务因子的福祉影响具有东西空间分异和集聚特征, 随时间变化呈现区域稳固性, 川西西部居民福祉受生态系统服务价值变量的减量效应明显, 相反, 阿坝州和凉山州东北部居民福祉受生态系统服务的增量效应显著。敏感程度具有从东南向西北递减的特征, 其中负向敏感性强于正向敏感性。

(4) 供给服务与文化服务是川西山区居民福祉的主导生态因子, 主控区域具有空间交替性。供给服务在2000年对川西东北部、2005年和2010年对川西西部居民福祉提升形成关键的抑制力, 2015年对川西东北部居民福祉改善提供了重要驱动力；2000和2015年川西西部县域主要受控于文化服务的抑制作用, 2005和2010年文化服务对川西东北部居民福祉的提升带来的驱动作用最大。

5 讨论

本研究从时空异质性角度出发, 理清了2000—2015年川西山区生态系统服务价值及其不同类型对居民福祉的驱动关系及其变化过程。川西山区是西部山区的典型代表, 作为国家级重点生态功能区、“中华水塔”的重要组成部分, 同时也是西部生态脆弱区、敏感区和贫困集中分布地, 准确把握生态系统服务的福祉效应差异, 明晰不同地域空间的关键驱动因子, 可推进山区生态系统服务功能的优化。已有研究表明不同生态系统服务类型与福祉的关联不同^[20-21], 对于贫困山区而言, 增强生态系统的文化、供给服务对增进农户福祉具有重要意义^[29]；对自然保护区而言, 调节服务福祉贡献最大；而对草原地区来说, 供给服务与牧民福祉关系最为紧密^[20]。但已有研究主要从农户感知视角出发, 不同农户因生产方式、教育水平、性别差异、收入构成不同会影响生态系统服务福祉贡献的判断, 除受益者特征外, 空间与尺度也是福祉传递的重要原因^[43], 相关研究根据生态系统服务与居民福祉的相对关系分为了不同耦合模式^{[24, 27]}, 但无法直接体现生态系统服务对居民福祉的地域驱动关系。本研究通过空间计量分析手段探讨了不同空间生态系统服务的福祉效应差异, 是对生态系统服务与福祉空间关系研究的重要补充, 有助于在中小尺度单元制定针对性的生态管理政策。

研究中采用较为成熟的生态系统服务价值当量因子法, 该方法以土地覆被数据为源, 可直观便捷的表征生态系统服务结构、过程、功能动态变化, 被广泛应用于区域ESV评估^{[4-5, 46]}。但该方法依赖于单位面积生态服务基础当量表, 我国地域辽阔, 生态系统具有复杂性、动态性和异质性特征, 不同时间不同地域同一生态系统提供的主要功能会有所差异^[47], 不同时空转换条件下当量因子系数的修正值得进一步探索, 也需通过实地调研印证所得结果与实际服务价值的吻合性。人类福祉是多维度概念, 鉴于数据有限性和代表性, 研究中采用应用最广泛的人类发展指数作为福祉的代理指标, 但目前的指标体系以人文因素为主, 福祉的多维特性要求在评价福祉过程中从客观、主观以及社会、经济、生态多视角考虑福祉构成, 使福祉评价内容和方法更具有科学性^[27]。生态系统服务的贡献是人类福祉体系的子集而非全部^[43], 居民福祉既依赖于生态系统服务, 同时也依赖于非生态因素^[48], 生态系统服务与福祉关系较为复杂, 服务到收益的传递具有非线性特征, 后续研究中将结合国内外研究进展, 融入政策环境、社会经济等因素, 对多因素共同作用下生态系统的福祉效应展开科学识别和传导机制的系统诊断。

厘清生态系统服务的福祉效应是山区人地关系调控的重要方面, 可为实现生态文明建设和可持续脱贫提供重要参考。生态系统服务无疑对人类福祉具有重要贡献^[43], 当前国家高度强调推动长江、黄河流域生态保护和高质量发展, 筑牢江河源区生态安全屏障, 其关键在于贯彻新发展理念, 将绿色发展融入转型跨越、脱贫奔康的全过程。一是实行最严格的生态环境保护制度。严控生态保护红线, 全面启动实施天然林、草原、湿地、长江上游干旱河谷等地区生态修复, 加强西南部石漠化水土流失脆弱地区生态治理, 统筹推进生态文明建设, 构建金沙江、大渡河、岷江、安宁河等生态廊道, 增强地区居民生态保护意识。二是加快生态与产业融合发展, 深入实施产业富民战略。推进攀西经济区转型发展和川西北生态示范区建设, 推进农牧业供给侧结构性改革, 坚持生态经济提质, 提升特色农牧业产业化水平, 挖掘文化基因, 推进民族地区特色旅游和康养业发展, 加强园区和生态展馆建设, 壮大绿色现代工业, 激活绿色资本转换动能。三是健全自然资源有偿使用制度和综合生态补偿机制。深化岷江等重点流域上下游横向生态保护补偿机制、资源有偿使用制度试点, 建立全省统一的长江黄河上游生态系统保护与修复投入长效机制。