生态系统服务是指自然生态系统以及物种所提供的能够满足和维持人类生活需要的条件和过程, 是人类直接或间接从生态系统获得的所有收益^[1]。由于生态系统服务的多样性、空间分布的不均衡性以及人类使用的选择性, 在人类活动及自然因素的共同作用下, 生态系统服务内部和不同生态系统服务之间呈现出了动态变化, 表现为此消彼长的权衡关系和相互促进的协同关系^[2]。充分认识生态系统服务间的权衡/协同关系是实现区域生态系统服务可持续供给的基础, 对区域生态系统服务优化具有重要意义^[3]。

喀斯特地貌是一种特殊的景观类型, 其在可溶性岩石上发育而成, 约占全球陆地总面积的15%^[4]。作为世界上最大的连续型喀斯特区, 中国西南喀斯特山区土壤厚度较薄(通常小于10 cm), 地下洞穴系统发育强烈, 地表径流严重不足, 加之强烈的人类开发活动, 造成了严重的生态环境问题(石漠化等), 威胁着区域生态系统服务功能的可持续供给, 确保区域生态系统服务安全迫在眉睫^[5]。目前, 将生态系统服务及其关系纳入西南喀斯特地区的研究较多, 主要集中于生态系统服务时空变化研究^[6]及其权衡/协同关系识别^[7]。在权衡/协同关系识别过程中, 学者们发现由于水土组合、光热条件、人类活动的差异, 喀斯特地区生态系统服务关系表现出了明显的空间异质性^[8]。于是, 学者们开始衡量其空间分异的程度, 比如利用变异系数^[9]和空间自相关指数^[10]等定量揭示其空间变异大小, 或者在不同梯度(环境)条件下探索生态系统服务关系的空间差异^[11]。在此基础上, 学者们还对喀斯特山区生态系统服务关系异质性产生的内在机理进行了剖析, 发现自然气候条件^[12]、生态环境因子^[13]、土地利用变化^[14]、生态修复工程^[15]等都会对其产生影响。只是, 这些研究多局限于一种或几种特定影响因素, 较少从自然条件-社会-经济-生态等多方面对生态系统服务关系的空间异质性进行解释。

贵州省作为中国西南典型的喀斯特地貌区, 生态本底条件较脆弱, 历史上由于自然与人为因素影响, 区域出现了严重的生态破坏和功能退化问题；为此, 近年来区域实施了一系列生态修复措施, 以增加地表植被覆盖, 改善区域生态环境。正如预期, 这些措施在增加区域植被覆盖^[16]、增加碳储量^[17]、减少土壤侵蚀^[18]等方面发挥了积极作用。同时, 也有研究指出生态修复措施的实施影响了地表产水, 增加了供给与调节服务间的冲突关系。只是这些研究多将贵州省作为一个整体研究, 忽略了其空间差异, 同时对生态系统服务关系变化的内在驱动机制也鲜有研究。基于此, 本文以2000—2015年为时间节点, 在对贵州省植被净初级生产力(Net primary productivity, NPP)、土壤保持(Soil conservation, SC)、产水(Water yield, WY)3类生态系统服务估算的基础上, 通过相关性分析、均方根误差法对生态系统服务关系变化及其分异特征进行分析, 并借助地理探测器探测各因子对区域生态系统服务关系变化的影响及因子间的交互作用关系, 旨在明晰生态系统服务关系变化的驱动机理, 以此为区域生态建设与工程布局提供科学参考。

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

贵州省地处中国云贵高原, 地理范围介于103°36′—109°35′E、24°37′—29°13′N(图 1)。全省包含89个行政区(县), 总面积约17.2万km², 约占全国总面积的1.8%。气候上属于亚热带湿润季风气候, 冬暖夏凉, 年均气温10—18℃, 年均降水量1000—1500 m^[19]。地质构造上主要由石灰岩、白云岩组成, 面积约为全省70%, 广泛发育着世界五大造型地貌之一的喀斯特地貌。区域由于自然本底条件的制约, 土层较薄, 土壤有机质含量较低, 基岩裸露, 生态环境敏感且脆弱。同时, 贵州省也是我国贫困问题较突出的省份之一, 人类活动频繁且不尽合理。在自然与人为因素的共同作用下, 区域出现了石漠化、水土流失等生态环境问题, 区域生态系统服务功能下降、供给与调节服务冲突不断, 制约着地区经济社会发展和生态协调。

1.2 数据来源与处理

省域、区(县)等边界数据来源于国家测绘地理信息局；地貌类型根据张雪梅等^[9]的研究绘制而成；岩性数据来源于全国1∶250万地质图空间数据库；数字高程数据来源于中国科学院资源环境中心的SRTMDEM-V2产品, 空间分辨率为30 m, 网址为https://www.resdc.cn/, 坡度数据是根据高程在ArcGIS软件中生成；土壤数据采用的是基于世界土壤数据库(HWSD)的中国土壤数据集, 来源于兰州寒区旱区科学数据中心；土地利用类型数据来源于中国科学院资源环境科学数据中心, 空间分辨率为30 m；气象数据, 主要包括气温和降水, 来源于中国气象数据网, 采用ANUSPLIN软件的薄盘样条法对其进行插值, 获得了多年逐月时间序列的栅格图像, 并通过最大合成法获取了年序列的栅格数据, 空间分辨率为1 km；归一化植被指数(Normalized Difference Vegetation Index, NDVI)和植被类型数据来源于中国科学院资源环境科学数据中心, 网址为https://www.resdc.cn/, 其中植被数据为1∶100万的矢量图层, NDVI为空间分辨率为1 km的栅格图层；实际蒸散发数据来源于国家青藏高原科学数据中心(http://data.tpdc.ac.cn/zh-hans/), 空间分辨率为0.1°；灯光数据来源于美国海洋和大气管理局地球观测组(https://www.ngdc.noaa.gov/eog/viirs.html), 空间分辨率为500 m。统计数据中的人均地区生产总值(Gross Domestic Product, GDP)、人口密度数据来源于贵州省统计年鉴, 造林总面积来源于中国林业统计年鉴。对于统计数据的空间化, 人均GDP和人口密度是通过灯光数据插值所得, 即计算各栅格DN值在区(县)总DN值上的占比, 然后乘以相应的统计数据^[20]；参考陈金珂等^[21]研究, 通过NDVI计算得到贵州省的植被覆盖度, 进而得到区域裸岩率的空间分布图。按照裸岩率越大的区域实施生态修复的可能性越大, 将造林总面积落实于空间。其他部分有关数据来源于相关文献或者采用的模型默认值, 例如InVEST模型计算土壤保持服务时所需的降水侵蚀因子R、土壤可蚀性因子K等主要是参考相关文献, 而泥沙输移比等校准参数采用的是模型默认值。在获得相关数据后, 对其进行了矢栅转换, 并将数据的坐标系和分辨率分别统一为WGS 1984和1 km。

1.3 研究方法 1.3.1 生态系统服务估算

植被净初级生产力(NPP): 基于CASA光能利用效率模型, 利用相关数据估算得到贵州省植被净初级生产力, 其计算公式如下：

(1)

(2)

式中, APAR(x, t)代表t月像元x吸收的光合有效辐射(MJ/m²)；ε(x, t)对应t月像元x的实际光能利用效率(gC MJ^-1)；SOL(x, t)表示t月像元x的太阳总辐射(MJ/m²)；FPAR(x, t)表示植被层对入射光合有效辐射的吸收系数。具体计算过程主要是参考朱文泉等^[22]的研究。

土壤保持(Soil conservation, SC): 本文基于InVEST模型中的泥沙输送比例模块计算了潜在土壤侵蚀量与实际土壤侵蚀量, 通过两者差值估算得到贵州省土壤保持服务, 计算公式具体如下:

(3)

(4)

(5)

式中, SEDRET_x和RKLS_x为栅格x的土壤保持能力和潜在土壤侵蚀量(t hm^-2 a^-1)；USLE_x是栅格x的实际土壤侵蚀量(t hm^-2 a^-1)；R_x是降水侵蚀因子(MJ mm hm^-2 h^-1 a^-1), 依据Wischmeier等^[23]提出的基于月均降水量和年均降水量的关系计算所得；LS_x为地形因子, 通过高程数据模拟而得；K_x为土壤可蚀性因子(t hm² h MJ^-1 hm^-2 mm^-1), 通过高敏等^[24]的研究和区域土壤类型数据获取；C_x为植被覆盖管理因子, P_x是水土保持因子, 通过韩会庆等^[25]的研究成果和区域土地利用类型而获得。

产水服务(Water yield, WY)：本文通过计算降水量与实际蒸散发量之间的差值, 用以量化区域产水服务。其计算公式如下:

(6)

式中, Y_x代表研究区栅格x的年产水量(mm/a)；AET_x是研究区栅格x的年实际蒸发量(mm/a)；P_x是栅格x的年降水量(mm/a)。

1.3.2 生态工程影响区识别

趋势分析法: 本文通过趋势分析计算了2000—2015年NDVI变化斜率(Slope), 用以明晰研究时段内NDVI变化趋势。其计算公式如下：

(7)

式中, Slope表示回归方程斜率；i为时间变量；n为研究时段长度即n=16；NDVI_i为像元在第i年的值；Slope＞0和Slope＜0分别表示NDVI随时间呈现出上升或下降趋势；Slope=0表示研究时段内NDVI无明显变化趋势。同时, 基于此得出贵州省2000—2015年NDVI的Slope平均值为0.0045, 用于后续研究。

偏相关分析: 气候影响着研究区NDVI变化。本文利用偏相关分析法计算得到NDVI与年均温、年降水总量间的偏相关系数, 用于检测气候与NDVI之间的相关程度。其计算公式如下：

(8)

(9)

式中, R_xy是2个变量间相关系数, x_i和y_i分别代表第i年的值, x和y代表平均值。R_{xy, z}是自变量z固定后因变量x与自变量y的偏相关系数。偏相关系数一般用t检验进行显著性检验, 其公式如下：

(10)

式中, t为显著性检测统计量, n(n=16)为样本量, m自变量。

基于以上研究方法, 本文参考相关研究^[26], 对生态工程显著影响区进行了识别, 该区域必须同时满足3个条件：(1) 2000—2015年NDVI水平必须得到显著提升；(2) NDVI与气候(年均温、降水总量)间正相关性不显著；(3) 年NDVI的Slope值大于区域平均值, 0.0045。据此条件得到生态工程显著影响区和非生态工程影响区, 结果如图 2所示。从图中可知, 研究区大部分区域NDVI的增长受到生态修复工程的显著影响, 面积达63.2%, 主要分布于研究区北部和西南部, 这与Tong等^[17]得到的研究结果相吻合。

1.3.3 权衡/协同关系评价

相关分析法：为揭示生态系统服务间的权衡/协同关系, 本文采用Pearson相关系数法计算生态系统服务间的相关性。在通过1%显著性检验的条件下, 相关性系数为正, 表明生态系统服务之间为协同关系, 反之则为权衡关系。

均方根误差法：均方根误差(Root Mean Square Error, RMSE)不仅可以表征生态系统服务间的协同/权衡关系, 也可以表达生态系统服务关系间同向变化率的不均衡性和同向间协同/冲突程度, 是量化两个或多个生态系统服务间权衡度的一个简单而有效的方法。其计算公式如下：

(11)

式中, RMSE为生态系统服务权衡强度, n为样本数量, ES_i为第i种标准化后的生态系统服务价值, ES为各类生态系统服务的期望值。

1.3.4 影响因子探测

指标因子选取  影响生态系统服务关系变化的因素复杂多样, 大致可以分为自然地理条件、社会经济发展、生态修复工程3方面。结合区域现状及数据的可获得性、代表性, 借鉴前人研究成果^[27], 自然地理条件方面选取坡度、高程、年均温、降水总量4个因子, 社会经济方面选取人均GDP和人口密度2个因子, 生态修复工程方面以2000年开始实施的退耕还林(还草)政策为例, 选取造林总面积这个因子, 形成贵州省生态系统服务关系变化的驱动因子集。

地理探测器  喀斯特地理现象具有明显的空间分层异质性, 而地理探测器是近年来探测空间分异性并揭示其驱动机理的新型统计学方法^[28], 由风险探测、因子探测、生态探测、交互探测4个模块组成。本文在明确贵州省生态系统服务间权衡关系空间分异特征后, 利用地理探测器模型中的生态探测比较驱动因子对权衡分布差异形成的影响是否存在显著性, 然后通过因子探测和交互探测定量解析各因子及因子相互作用对生态系统服务权衡的影响, 旨在明晰区域生态系统服务权衡空间分异性的驱动机理。

生态探测主要用于揭示自变量解释力相对重要性的差异, 以F统计量来衡量, 计算公式如下：

(12)

式中, N_x1和N_x2分别表示两个因子的样本数量；L₁和L₂表示自变量x1和x2的分层数。

因子探测的核心是若存在某个自变量x对某个因变量y有影响, 那么自变量的空间分布和因变量的空间分布应该具有趋同性。其计算公式如下：

(13)

式中, q为探测因子对因变量的探测力值, q∈[0, 1], q值越大自变量x对因变量y的解释力越强, 反之则越弱；N, σ²代表研究区的总样本量和方差；L表示探测分区数；N_i, σ_i²代表不同分区的样本量和方差。

交互探测主要是通过计算自变量x1与x2对于因变量y的q值及x1与x2交互后对y的q值, 然后比对单因子的q值与双因子交互作用后的q值, 以此判定交互作用的类型和方向, 其交互关系可以分为5类, 具体见表 1。

2 结果与分析 2.1 生态系统服务时空变化特征

通过与NPP产品、统计公报中的径流深、侵蚀模数数据进行对比以及统计分析发现, 本文生态系统服务的估算效果良好, 结果较为可信, 其供给空间分布及变化趋势如图 3。2000—2015年, 区域NPP均值从403.3 gC m^-2 a^-1增加到了481.2 gC m^-2 a^-1, 变化幅度为-579.5—1243.2 gC m^-2 a^-1；土壤保持量均值从422.1 t hm^-2 a^-1增加到了512.8 t hm^-2 a^-1, 变化幅度为-11580.8—14477.3 t hm^-2 a^-1；产水服务从373.1 mm/a增长到了397.1 mm/a, 变化幅度为-943.2—1020.7 mm/a。整体而言, 3类生态系统服务均值均呈增长趋势。同时, 各项生态系统服务在空间分布上差异明显：NPP和产水服务均呈现出了北部低南部高的分布格局；土壤保持服务高值区分布于贵州省北部和南部, 中部属于低值区。就变化量的空间分布而言, NPP和土壤保持服务在除研究区东南部和西南部的部分区域外均呈增长趋势；产水服务在研究区北部和西部呈现出了减少趋势, 东南部增加趋势最明显。

2.2 全域生态系统服务关系探讨

借助ArcGIS软件, 通过建立渔网对研究区3类生态系统服务间的相关性进行分析, 结果如表 2所示。2000—2015年, NPP、产水服务、土壤保持间相关性变化明显。其中, NPP与土壤保持, NPP与产水服务间存在明显的正相关性, 即协同关系, 只是相关性系数在减小, 说明协同性在降低。土壤保持与产水服务从2000年的正相关性变成了2015年的相关性不显著。整体而言, 研究区生态系统服务关系以协同为主, 但协同关系的正相关性系数在降低；同时在2015年生态系统服务间关系有向权衡转移的趋势, 应引起适当关注。

在明确了生态系统服务间的相关性后, 本文采用均方根误差法对区域生态系统服务间权衡度进行了计算, 结果如图 4。2000—2015年, 土壤保持与产水服务间权衡度较大, 矛盾最突出。就权衡度的变化而言, NPP与土壤保持间的权衡度从2000年的0.159降到了0.138, 而NPP与产水服务, 土壤保持与产水服务间的权衡度由2000年的0.129和0.161分别增加到了2015年的0.165和0.188。可见, 产水服务与其它两类生态系统服务间权衡度均呈增加趋势, 水资源是制约区域生态系统服务平衡的重要因素。

2.3 生态系统服务关系的空间分异特征

本文从岩性、地貌、土地利用及生态工程影响出发, 计算了不同背景条件下生态系统服务间的权衡度, 用以剖析生态系统服务关系的空间分异特征, 结果如图 5所示。就不同岩性而言, 非岩溶地区NPP与土壤保持间权衡度高于岩溶地区, 可见岩溶区的生态修复促进了区域植被生长、土壤保持；同时, 白云岩碎屑岩互层上NPP与产水服务的权衡度较高, 石灰岩碎屑岩互层上土壤保持与产水服务间的权衡度较高。就连续性石灰岩和白云岩而言, 石灰岩上NPP与土壤保持、土壤保持与产水服务间的权衡度要高于白云岩, 这可能与石灰岩本身更易溶蚀有关。

就地貌而言, 岩溶盆地上NPP与土壤保持间的权衡度较高, 岩溶峡谷上土壤保持与产水服务间的权衡度较高, 岩溶槽谷上NPP与产水服务间的权衡度较高, 整体表现为平坦的地区土壤保持服务与NPP间冲突较明显。就权衡度的变化而言, 岩溶盆地NPP与土壤保持间的权衡度减小最明显, 而峰丛洼地产水服务与其它两类生态系统服务间的权衡度减小趋势最明显, 这可能与区域蓄水条件的改善有关。

就土地利用而言, 建设用地上NPP与土壤保持、NPP与产水服务间的权衡度略高于高植被覆盖度地类(耕地、林地和草地), 而高植被覆盖度地类(耕地和林地)上的土壤保持与产水服务间的权衡度更高, 表明随着植被覆盖度的提升土壤保持与产水服务间的冲突会增强, 这也侧面暗示出以退耕还林(还草)为代表的生态修复工程不能仅追求植被覆盖度的提高, 同时也要关注生态系统服务间冲突关系及其变化。

就生态工程影响而言, 生态工程显著影响区产水服务与土壤保持间的权衡度略高于非生态工程影响区, 这与上述关于高植被覆盖度地类上产水服务与土壤保持间冲突更明显的结果相同；同时, 生态工程显著影响区NPP与产水服务间的权衡度较小, NPP与土壤保持间权衡度的下降程度高于非生态工程影响区, 这说明生态修复工程在一定程度上缓解了NPP与土壤保持、产水服务间的矛盾。

2.4 生态系统服务权衡变化的驱动因子解析

本文首先根据生态探测对驱动因子对生态系统服务权衡关系空间分异影响的显著性进行了检验。由表中可知：共有13对因子通过了NPP与土壤保持间权衡分异的显著性检验, 主要集中在坡度、造林总面积、高程、人均GDP等因子(表 3)；也有13对因子通过了NPP与产水服务间权衡分异的显著性检验, 主要为高程、人均GDP、降水等因子(表 4)；共有17对因子通过了土壤保持与产水服务间权衡分异的显著性检验, 集中于高程、人均GDP、降水等因子(表 5)。

在此基础上, 对研究区生态系统服务权衡关系变化的单因子进行了识别, 结果如图 6。就NPP与土壤保持间的权衡关系而言, 生态修复工程的影响力大于其它因子, 主要表现为造林总面积的解释率最高, 达40.2%(q值=0.402), 其次是高程与人均GDP, 解释率分别为30.9%和26.2%, 最后是降水总量(12.4%)和年均温(11.6%)。就NPP与产水服务间的权衡关系而言, 自然因子的影响力较强, 降水总量的解释率为17.5%, 其次是人均GDP和高程, 解释率分别为8%和7.7%。就土壤保持与产水服务间的权衡关系而言, 自然和经济因子的影响较强, 其中, 高程的解释率最大, 达21.3%, 其次是人均GDP, 解释率为19%。对人口密度而言, 其对贵州省生态系统服务权衡关系变化的影响并不显著。整体来看, 两两生态系统服务间权衡变化的影响因子不尽相同, 其中, 造林总面积是影响NPP与土壤保持间权衡变化的主导因子, 降水总量是影响NPP与产水服务间权衡变化的主导因子, 高程与人均GDP是造成土壤保持与产水服务间权衡变化的主导因子。

在对主导因子进行识别后, 本文又对因子间的交互作用进行了探讨, 结果表明(图 7), 任意两个因子交互作用的解释力均大于单因子驱动作用, 交互作用后因子的影响强度明显提升。就NPP与土壤保持间权衡关系而言, 其受造林总面积与高程的交互影响最强, 影响力为0.562, 其次是造林总面积与人均GDP、造林总面积与降水总量, 影响力均在0.4以上；NPP与产水服务间权衡关系受降水总量与高程的交互影响最强, 影响力为0.293, 其次是降水总量与人均GDP、降水总量与造林总面积。土壤保持与产水服务间权衡关系受降水总量与高程的交互影响较明显, 影响力为0.394, 其次是年均温与高程、高程与人均GDP, 影响力均大于0.3。整体来看, NPP与土壤保持间受因子交互影响作用最强, 其次是土壤保持与产水服务, 最弱的是NPP与产水服务。

3 讨论 3.1 生态系统服务权衡/协同特征

喀斯特地区由于特殊的自然地理环境及不合理的人类活动, 生态环境问题频发, 生态系统矛盾突出。为此, 区域开展了一系列生态修复工作, 以提升生态系统服务的综合水平。已有研究表明, 生态修复工程实施过程中喀斯特地区供给服务得到了明显提升, 而某些调节服务却表现出了下降趋势^{[12, 29]}, 也就是说生态系统服务间存在权衡/协同现象。通过本文研究发现, NPP、土壤保持、产水服务间确实存在一定的权衡关系, 特别是产水与其它两类生态系统服务间冲突明显。而从生态系统服务的角度来说, 生态保护和修复工程的目标应该是减少权衡、促进协同。准确识别区域生态系统服务权衡/协同关系及变化规律具有重要价值。

喀斯特地区是一个复杂的二元空间异质体结构, 其地貌、水热、植被、土壤条件具有明显的空间非平稳性, 这在一定程度上影响了生态系统服务及其关系的空间表达^[9]。从本文研究来看, 不同背景条件下, NPP与土壤保持、产水服务间权衡关系呈现了一定的空间分异特征。在岩性上, 非岩溶地区NPP与土壤保持间的权衡度较大, 这一结论从侧面印证了岩溶区生态修复工程对区域土壤保持的重要意义。相较于白云岩, 石灰岩更容易被溶蚀, 所以其土壤保持与产水服务间的权衡度更大, 这与王欢等^[30]得出的石灰岩、碎屑岩土壤流失程度高于白云岩的结论相契合。岩溶盆地上土壤保持与NPP间的权衡度较高, 主要是由于其地形平坦、光照充足、人类活动频繁, 人口压力较大, 存在较严重的土地退化、植被破坏问题^[31]。生态修复工程显著影响区及高植被覆盖区域, 土壤保持与产水服务间的权衡度更大, 原因在于高植被覆盖可以增强土壤抗蚀能力、减少土壤流失, 同时植被根系也具有较大的吸水能力, 会间接减少地表产水^[32]。为此, 在未来生态修复过程中要改进喀斯特区域背景条件下的措施类型, 实现生态系统的精准修复。

3.2 生态系统服务权衡的驱动机制

本文从自然气候条件-社会经济-生态政策3方面选取指标, 综合剖析贵州省生态系统服务权衡关系变化的内在机制。结果显示, 造林总面积是影响NPP与土壤保持间权衡关系变化的最主要因素, 这反映了生态修复工程在区域植被恢复、土壤修复方面的重要作用^[33]。降水总量是影响NPP与产水服务间权衡关系的最主要因素, 这与学者们得到的降水是影响植被生长和地表产水最重要的要素的结论相吻合^[13]。高程是影响土壤保持与产水服务间权衡关系变化的最重要因素, 主要是因为高程会影响喀斯特地区地表径流和土壤流失的强度^[34]。人口密度对区域生态系统服务关系变化影响不明显, 原因在于贵州省社会经济发展相对落后, 与快速城市化地区相比, 人口聚集程度相对较低, 包括气候、生态政策在内的其它因素可能会带来更大影响^[12], 本研究中气候条件、生态修复政策及地形对生态系统服务权衡变化的主导作用也证明了这一点。同时, 各因子对区域生态系统服务关系变化的影响并非完全独立, 而是会产生一定的交互作用, 并且交互作用起到了决定性的作用^[35]。未来在对影响因子进行解析时应更加关注因子间交互作用对生态系统服务关系的影响, 并尽可能的将因子影响落实于空间, 以更好地服务生态管理与决策。

3.3 不确定性分析与展望

本文通过相关模型和方法对3类生态系统服务进行估算, 对于模型中的常量采用的是默认值或参考相关研究, 虽然生态系统服务整体变化趋势与Peng等^[12]的研究结果相一致, 但对模型参数阈值等问题未进行深入探讨。在对造林总面积进行空间化时, 由于可参考性文献较少, 本文以地表裸露率为依据对其进行空间化, 虽然结果符合区域实际, 但内在科学依据需进一步梳理。同时, 喀斯特生态系统具有多种服务, 而本文仅考虑了3类, 未来可以将多种生态系统服务纳入研究。

生态和地理现象的空间异质性违反了经典统计学中的独立和相同分布假设^[36], 而基于空间分层异质性的地理探测器可以定量检测驱动因素及其相互作用对生态系统服务关系变化的影响, 这在一定程度上又有效地减少了信息冗余和共线性, 为本文提供了重要的方法支撑。同时, 时间和空间是生态系统服务及其关系研究的重要视角, 未来可以考虑多时间节点或更长时间序列, 系统性开展喀斯特地区生态系统服务关系的空间分异性及驱动机制研究, 以获得更客观的结果。

4 结论

基于对2000—2015年贵州省NPP、土壤保持、产水服务的模拟, 本文对不同背景条件下3类生态系统服务关系及其分异特征进行了剖析, 并通过地理探测器对其权衡关系变化的内在驱动机制进行了探讨。结论如下：

(1) 研究时段内贵州省3类生态系统服务均值均呈增长趋势, 表明生态系统服务供给能力的提升。生态系统服务关系以协同为主, 但协同关系在减弱, 并且在2015年表现出了一定程度的权衡趋势。

(2) 贵州省生态系统服务关系具有明显的分异特征。在岩性上, 由于石灰岩更容易被溶蚀, 因此相对于白云岩, 其土壤保持与产水服务间面临着更大的冲突；随着不同地貌类型的变化, 地势平坦地区即岩溶盆地上土壤保持与NPP间权衡度更大, 而陡坡地区NPP与产水服务、土壤保持与产水服务间的权衡度较大；由于植被显著的保土保水能力, 随着植被覆盖度的提高, NPP与产水服务间的冲突更明显。

(3) 研究区生态系统服务权衡关系受多种因素影响, 其中, 造林总面积对NPP与土壤保持间权衡关系的解释率最大, 达40.2%；降水总量对NPP与产水服务间权衡关系的解释率较大, 为17.5%；高程对土壤保持与产水服务间权衡关系解释率达21.3%。同时, 因子之间的交互作用均高于单因子的影响, 对生态系统服务权衡关系变化的影响更明显。