随着1997年Costanza等^[1]发表在Nature上的一篇关于核算全球生态系统服务价值的文章以及2001年联合国千年生态系统评估项目^[2]的开展, 生态系统服务得到了越来越多的关注。生态系统服务作为连接自然和人类社会的桥梁和纽带, 是指生态系统与生态过程所形成及所维持的人类赖以生存的自然环境条件与效用, 提供了人类赖以生存和发展的资源与环境基础^[3-4]。对生态系统服务的量化评估可以为生态安全、土地规划等部门提供科学定量的生态环境信息, 将有利于开展针对生态环境的管理和决策, 这也使得近年来通过模型模拟的方式评估一区域生态系统服务并分析它们的特征与相互关系成为了热点研究领域^[5-6]。如Li等^[6]以潮白河沿岸的城市带为研究区, 模拟了当地固碳、土壤保持、防风固沙、产水以及生境质量五项生态系统服务, 并分析了五项服务之间的相互关系。Wang等^[5]分析了京津冀地区产水、固碳、土壤侵蚀等服务的多年变化及权衡协同关系。

然而不同地区有着不同的环境及社会本底特征, 各项生态系统服务的供给量必定有所区别, 由此造成了生态系统服务在空间分布上的异质性^[7-8], 识别生态系统服务空间异质性的主要驱动因素成为了一项亟待解决的问题。目前对生态系统服务空间异质性驱动因子的研究多是通过建立生态系统服务与影响因子的相关关系、回归模型或是通过情景模拟的方式分析驱动因子在不同情境下对生态系统服务的影响, 如饶恩明等^[9]和张玲玲等^[10]用建立相关的方式分析了土壤保持、产水、生境质量等服务与人口密度、GDP、气候因子等的关系；侯文娟等^[11]通过地理加权回归分析了产流量与NDVI、景观破碎化指数等因子的相互关系；Bai等^[12]和Clerici等^[13]则通过设置不同的情景, 揭示气候变化和土地利用类型对美国肯塔基州以及安第斯山脉东部两个小流域生态系统服务时空格局的影响。上述研究虽然都在空间上分析了生态系统服务空间格局与他们影响因素间的关系, 但没有阐明具体哪些影响因素对该生态系统服务空间异质性的形成起到主要作用。地理探测器是用来分析数据空间异质性的有力工具^[14], 可用于揭示数据空间异质性背后的主要驱动因子, 并定量分析单个驱动因子及因子交互对因变量的解释力度。

探索生态系统服务空间异质性并开展其主要影响因子的识别对地区生态系统服务功能调控和优化具有重要意义。那曲市位于青藏高原的中心, 是一个典型的生态脆弱地区^[15]。同时那曲也是我国“两屏三带”^[16]生态安全战略格局中青藏高原生态安全屏障的核心区域, 在水源涵养、气候调节、生物多样性维持等方面都具有重要的生态意义。随着2018年5月那曲撤地设市, 那曲市各层级、各项社会经济职能也逐渐的步入正轨, 在未来经济建设过程中思考如何协调生态保护, 做到生态保护与民生改善同步提升意义重大。因此, 本文以那曲市为例, 详细介绍了如何运用模型评估其多年的生态系统服务并分析生态系统服务的空间格局和相互关系, 随后引入地理探测器, 选择气候、植被指数、土壤、地形、人类活动及土地利用等影响因子探究各生态系统服务空间异质性的主要影响因素, 试图找出影响各生态系统服务空间分布的主要影响因子, 研究结果可以为那曲市未来的生态系统功能优化、生态安全格局构建及生态文明建设等提供科学的决策依据。

1 研究区概况

那曲市(83°55′E-5°5′E, 29°55′N-36°30′N)位于青藏高原生态安全屏障的核心区域, 长江、怒江、澜沧江、拉萨河等大江大河皆发源于那曲或支流流经那曲^[17], 水资源极其丰富。那曲市湖泊湿地遍布且孕育着西藏最大面积的天然草场, 在涵养水土、调节气候的同时也为大量高原动植物物种提供了必要的生存场所, 那曲市提供的多项生态系统服务功能必将会为整个高原乃至全国的生态安全保障做出积极贡献。那曲市平均海拔约4500 m, 整体地势由西北向东南倾斜, 西北部地势高亢, 海拔高度平均在5500 m左右, 中西部地势开阔, 属于高原丘陵地形, 东部属河谷地带, 多高山峡谷且山势陡峭, 海拔高度在3500-4500 m之间。那曲市属亚寒带气候区, 基本特点是气温低, 空气稀薄, 大风天气多, 太阳辐射强, 同时其东南部地区易受夏季风影响, 降水与温度条件明显好于中西部地区, 那曲市的东南部地区也是藏北仅有的农作物产区, 人口相对稠密, 同时拥有少量的森林资源和灌木草场。那曲市的土地利用类型包括耕地、林地、草地、水域、建设用地、未利用地、永久性冰川积雪七大类, 其中覆盖面积最大的草地约占那曲市总面积的52%, 林地、耕地、建设用地的占比较低(图 1)。

2 数据来源及研究方法 2.1 数据来源

本研究使用的数据主要包括：土地利用数据、气象数据、土壤数据、高程数据、归一化植被指数(NDVI)、行政区划数据、道路数据、人口密度数据等。其中2000、2010、2018年3期30 m分辨率的土地利用数据、土壤类型数据以及研究区县级行政区划数据均来自于中国科学院资源环境科学数据中心(http://www.resdc.cn/)。气象数据中的气温和降水使用了国家气象局CN05.1格点化观测数据集^[18], 日照时数数据来自中国气象科学数据共享服务网(http://data.cma.cn/)。土壤质地、土壤有机质含量、土壤容重等数据通过联合国粮农组织(FAO)世界土壤数据库以及国家青藏高原科学数据中心(http://data.tpdc.ac.cn)获得^[19]。DEM获取自地理空间数据云(www.gscloud.cn/), 并通过ArcGIS 10.3软件提取研究区的坡度。NDVI使用了美国地质调查局(USGS)陆地过程分布式数据档案中心(http://lpdaac.usgs.gov/)的MOD13Q1数据。研究区的道路数据来自于国家基础地理信息中心1：25万全国基础地理数据库(https://www.webmap.cn/)。人口密度数据来自于WorldPop的人口密度数据集(https://www.worldpop.org/)。香农多样性指数(SHDI)借助Fragstats 4.2软件计算得到。上述数据在预处理完成后均采用阿尔伯斯等积圆锥投影(Krasovsky_1940_Albers)并重采样为1 km分辨率用于后续分析。由于获取到的各类数据的年份集中在2000-2018年, 且那曲市自然环境与社会发展在此期间均表现出一定变化, 故本文选择2000年、2010年、2018年3个年份的数据进行后续研究。

2.2 研究方法 2.2.1 生态系统服务估算

(1) 产水服务

本研究使用InVEST模型的产水量模块进行了那曲市产水服务的模拟运算。产水模块运用研究区的降水量、蒸散量、土壤根系深度、植物可利用含水量、土地利用等数据, 基于Budyko水热耦合平衡假设^[20]计算研究区每个栅格单元的年总产水量。公式如下：

(1)

(2)

式中, AET_i表示栅格单元i上的年蒸散量(mm)；P_i表示栅格i上的年降水量(mm)；Y_i为栅格i上全年的产水深度(mm)；PET_i为栅格i上的潜在蒸散量(mm)；w_i为经验参数。

(3)

(4)

(5)

式中, K_c(l_i)表示栅格i上的蒸散系数, 数值的确定参考了模型手册^[21]及文献[22]；ET_0i为栅格i的参考作物蒸散, 通过Modified-Hargreaves法^[23]计算得到；root.depth为根系深度参考了模型手册^[21]及文献[24-25], Z为经验常数；PAWC_i为植物可利用含水量^[26], SAND%、SILT%、CLAY%、OM%分别表示土壤砂粒、粉粒、黏粒、有机质的含量。

(2) 土壤保持服务

InVEST土壤保持模块通过气象数据、土地利用数据、土壤数据以及DEM等数据计算研究区各栅格单元的潜在土壤侵蚀量和实际土壤侵蚀量, 两者之差则为研究区的土壤保持量, 本文将采用此方法模拟那曲市各栅格单元上的土壤保持量。公式为：

(6)

(7)

(8)

式中, RKLS_i表示栅格i上裸地的土壤侵蚀量(t/hm²)；USLE_i为栅格单元i在植被覆盖管理因子和水土保持措施因子作用下的土壤侵蚀量(t/hm²)；SD_i为栅格单元i上的实际土壤保持量(t/hm²)；R_i是降雨侵蚀力因子(MJ mm hm^-2 h^-1), 采用Wischmeier的月尺度公式^[27]计算；K_i是土壤可蚀性因子(t hm² h hm^-2 MJ^-1 mm^-1), 采用EPIC模型^[28]估算, 并通过修正公式^[29]进行本地化矫正；LS_i为坡长坡度因子, 参考InVEST模型手册与已有研究^{[9, 21]}, 模型中由DEM数据计算得到；C_i为植被覆盖管理因子, 参考了蔡崇法等^[30-31]的研究, 通过NDVI计算植被覆盖度最终由植被覆盖度计算得到各栅格单元的C_i值；P_i为水土保持措施因子, 参考书籍^[32]与前人研究^{[31, 33]}获得。

(3) 生境质量

InVEST生境质量模块以土地利用数据和威胁源数据为基础, 计算各栅格单元的生境质量得分, 得分在0-1之间, 数值越大表示生境质量越好^[21]。本文结合研究区实际情况与相关研究^[34-36], 在选择耕地、建设用地、未利用地、国道、铁路为威胁源数据的同时构建了生境质量参数表, 最终运行InVEST生境质量模块得到了那曲市各栅格单元的生境质量得分。生境质量计算公式如下：

(9)

(10)

式中, Q_ij为栅格i的生境质量；H_j为土地利用类型j的生境适宜度；D_ij为栅格i的栖息地退化程度；k为半饱和常数；Z为默认参数2.5；R为威胁因子的数量；P_r是r威胁因子的总栅格数；W_r是威胁因子r的权重；r_p是威胁强度；i_rip表示r_p对i的威胁程度；β_i表示i的可达性；S_jr土地利用类型j对r威胁因子的敏感性。

(4) 固碳服务

本文运用CASA模型^[37]模拟那曲市各栅格单元上植被对碳的吸收, 即固碳服务。CASA模型是一个典型的光能利用率模型, 用于模拟植被的净初级生产力(NPP), 该模型模拟植被NPP所需的参数相对较少, 同时部分参量可直接通过遥感手段获得, 使之成为目前应用广泛且不断发展的NPP模型之一^[38]。运用CASA模型进行NPP计算所需的数据包括土地利用数据、气象数据、NDVI数据等, 计算公式如下^[39]：

(11)

式中, NPP_i为栅格单元i处的植被净初级生产力(gC/m²)；R_s表示地表太阳辐射(MJ/m²), 参照童成立等^[40]的研究, 基于日照时数数据算得；0.5为经验常数, 与R_s相乘近似表示植物可利用的地表太阳辐射；FPAR为植物转化光合有效辐射的比例；f₁和f₂为低温和高温对植物光能利用率的影响；W为水分限制因子；e_max植被的最大光能利用率, 数值参考朱文泉等的研究成果^[39]。

2.2.2 生态系统服务综合热点区域

对研究区生态系统服务综合热点区域的识别可以为当地生态保护措施的制定与决策提供定性的依据。本研究在栅格尺度上将各项生态系统服务供给量超过各自均值的区域视为该项生态系统服务供给的热点区域, 通过叠加研究区各对应年份的产水量、土壤保持、生境质量、固碳4项服务的热点区域, 最终得到研究区各对应年份的生态系统服务综合热点区域的分布。根据区域上单项生态系统服务热点的数量, 本研究将热点区域划分为极重要区域、高度重要区域、中度重要区域、一般重要区域、非重要区域五类(表 1)。

2.2.3 生态系统服务间的相互关系

本研究基于各生态系统服务的模拟值进行了生态系统服务间逐栅格点的相关分析并检验其显著性水平, 探究那曲市各项生态系统服务间的相互关系。依据spearman相关分析, 得到4种生态系统服务之间的相关性。公式为：

(12)

式中, R_ab为两种生态系统服务的相关系数；a_i、b_i为生态系统服务a、b的第i个栅格值；a、b为生态系统服务a、b的均值；n为样本数。

相关性结果为正值时, 表明两种生态系统服务之间具有协同效应, 即两种生态系统服务在同一时间段具有同样的上升或降低趋势, 一种服务的增加会对另一种服务产生一定的促进和增幅作用, 结果为负值时, 表明两种生态系统服务具有权衡关系, 即某生态系统服务的增加会造成另一种生态服务供给量的下降。

2.2.4 生态系统服务空间分异的驱动因子识别

地理现象具有空间异质性, 地理探测器是揭示生态系统服务空间异质性背后主要驱动因子的有力工具, 它的核心假设为：若某自变量会对因变量造成一定影响, 则它们在空间上的分布会存在一致性^[42]。地理探测器可量化各驱动因子对地理数据空间异质性的解释度, 结果用q统计量表示, q值的含义为自变量解释了q×100%的因变量。为了探究各单项生态系统服务受哪些影响因素的驱动, 本文以研究区各单项生态系统服务为因变量, 以影响那曲市生态系统服务的各因素为自变量, 选用地理探测器中的“因子探测器”和“交互探测器”分析各驱动因子及因子组合对研究区生态系统服务空间异质性的解释程度。因子探测公式如下^[42]：

(13)

(14)

(15)

式中, i表示自变量或因变量的分层；N_i表示i层的单元数量, N表示全域的单元数量；σ_i²表示层i内因变量的方差, σ²表示全域内因变量的方差；SSW为不同层方差之和；SST为全区总方差。

在生态系统服务驱动因子的选择方面, 以那曲市生态环境与社会经济背景为特征, 最终考虑了气候、植被指数、土壤、地形、人类活动及土地利用六类共九个指标。气候要素中选择了年降水、年均温；植被指数要素中选择了归一化植被指数(NDVI)；土壤要素中选择土壤类型为主要驱动因子；地形要素中选取了高程和坡度因子；人类活动要素则主要考虑了人口密度；土地利用方面的驱动因子选择了土地利用类型及表示土地利用类型破碎化的香农多样性指数(SHDI)。

交互探测用于评估不同驱动因子组合在一起对因变量的影响程度, 两因子的交互作用共有5种类型^[42], 若两因子交互作用q值小于任意单因子q值, 为非线性减弱；两因子交互作用q值在两个单因子q值之间, 为单因子非线性减弱；两因子交互作用q值大于任意单因子q值, 为双因子增强；两因子交互作用q值等于两个单因子q值之和, 为独立关系；两因子交互作用q值大于两个单因子q值之和, 为非线性增强。

3 结果与分析 3.1 那曲市生态系统服务随时间变化趋势

2000-2018年那曲市产水量(图 2)总体呈明显的下降趋势, 其平均产水深度由174.45 mm下降至113.19 mm, 下降了35.1%。各县区多年间的变化趋势分为两种类型, 中部与东部县区多表现为先降后增, 西北部县区多为逐年下降。那曲市各县区2018年的平均产水深度普遍小于2000年的水平, 尼玛县、双湖县、申扎县的平均产水深度分别降低了76.9%、72.5%和68.5%。

那曲市土壤保持服务(图 2)在2000-2018年总体上表现为先降低后升高的趋势, 单位面积的土壤保持量由2000年的42.34 t/hm²下降至2010年的36.19 t/hm²后, 随后升高至2018年的40.32 t/hm², 2018年单位面积土壤保持量相较于2000年降低了4.8%, 总体变化较小。各县区多年间的变化趋势与产水服务类似。2018年巴青县、色尼区、聂荣县、索县单位面积的土壤保持量高于2000年, 其余县区的土壤保持量最终均低于2000年。

那曲市固碳服务(图 2)在2000-2018年总体呈上升趋势, 单位面积固碳量由2000年的219.07 gC/m²增长至2018年的231.40 gC/m², 增长了5.6%。那曲市各县区2000-2018年的单位面积固碳量除申扎县和班戈县略微下降以外其余县区均表现为增长, 嘉黎县、双湖县、安多县分别增长了11.0%、10.7%和9.8%。

2000-2018年那曲市生境质量(图 2)趋于稳定, 2000年、2010年和2018年那曲市生境质量平均得分分别为0.364117、0.364131和0.365114, 2018年得分相较于2000年提高了0.3%。那曲市各县区的生境质量得分在2000年至2018年间变化很小。

3.2 那曲市生态系统服务的空间异质性 3.2.1 单项生态系统服务的空间异质性

那曲市各项生态系统服务的空间分布格局如图 3所示。各单项生态系统服务的分布都呈现出明显的空间分异特征, 同时各项服务的分布趋势又存在一致性, 即东南部地区的各单项生态系统服务普遍高于西北部地区。

那曲市全市多年平均产水深度达133.79 mm, 产水服务的空间分异明显。那曲市产水服务的分布呈现出东南高西北低的特点, 与该地区降雨量的分布特征基本吻合。产水量高值区主要分布在那曲市东部水热条件较好的嘉黎县、索县、巴青县、比如县, 该地区多高山峡谷且夏季受印度洋季风的影响降水丰沛, 四县多年平均产水深度分别达到492.63 mm、486.11 mm、475.71 mm和470.66 mm, 产水服务显著高于西部各县。产水量低值区主要分布在广大的西部和北部地区, 该区地势高亢, 常年受西风带影响, 降水稀少。尼玛县、双湖县多年平均产水深度仅为38.76 mm和45.23 mm。

那曲市单位面积土壤保持量的多年均值为39.62 t/hm², 土壤保持服务的空间分布格局与产水服务类似。土壤保持服务的高值区主要分布在东南部县区的高山峡谷区域内, 东南部区域丰沛的降水会使该区域受降雨侵蚀力的作用高于西北部, 同时山区坡度变化大的区域会存蓄大量的流失土壤, 东南部山地的山腰和山脚发育有那曲仅有的森林, 保土能力强。索县、嘉黎县、比如县、巴青县多年单位面积土壤保持量分别为295.55 t/hm²、222.85 t/hm²、197.04 t/hm²、143.74 t/hm², 显著高于全市平均水平。而那曲市西部及北部的双湖县、安多县、尼玛县多年单位面积土壤保持量仅为6.02 t/hm²、10.49 t/hm²、20.28 t/hm², 明显低于全市平均水平。

那曲市固碳服务的空间分布也具有明显的分异性, 全市单位面积固碳量多年均值为224.40 gC/m²。单位面积固碳量高的区域分布在那曲市中部和东部地区, 该地区气候条件优于西部与北部, 同时该地区分布着大面积的高寒灌丛草甸与高寒草甸, 西部与北部地区则主要为高寒草原。东部的比如县、嘉黎县、索县的部分高山峡谷中还蕴机着以冷杉、柏木、高山栎等为主的森林。索县、聂荣县、色尼区、巴青县单位面积固碳量多年均值分别达到489.16 gC/m²、478.12 gC/m²、436.44 gC/m²、411.75 gC/m²。低值区双湖县、尼玛县单位面积固碳量多年均值为135.98 gC/m²和174.39 gC/m²。

那曲市生境质量得分的多年均值为0.36。生境质量高值区主要分布在湖泊、河流及植被覆盖状况较好的区域, 主要分布在那曲市的中部与东部地区。索县、班戈县、安多县、比如县的多年得分均值分别达到0.54、0.48、0.45、0.44。低值区即尼玛县、双湖县的多年得分均值仅为0.31和0.32。

3.2.2 生态系统服务热点区域的的空间异质性

如图 4所示, 那曲市生态系统服务综合热点区域的分布在2000-2018年间基本一致, 各项生态系统服务高值区的重叠比例即高度重要区域与极重要区域的占比较低。其中, 极重要区域在历年来的占比保持在5.3%-5.4%之间, 该区主要分布在那曲市中部与东部的高山峡谷与河谷地带, 植被与水热条件明显优于其他地区。高度重要区域在历年来的占比处于12.3%-13.8%之间, 分布在那曲市中部与东部极重要区域的外围。中度重要区域和一般重要区域的占比在多年间分别保持在16.2%-19.1%与32.2%-34.1%之间。非重要区域即没有一项生态系统服务供给量超越其均值的区域占比处于28.6%-32.1%之间, 该区域主要分布于那曲市西部与北部的尼玛县、双湖县及安多县部分地区, 该区普遍海拔较高且干旱, 温度与降水条件明显劣于那曲中部与东部区域, 环境条件极为严苛, 那曲北部区域也被划归在著名的“羌塘无人区”之中。

3.2.3 生态系统服务间的相互关系

在对那曲市四项生态系统服务进行估算之后发现各生态系统服务的时空分布格局具有一定的一致性, 各年份生态系统服务间的权衡与协同关系如表 2所示, 除固碳服务和生境质量服务在2000年与2010年没通过显著性检验之外, 其余年份各项生态系统服务之间的相关性均通过了0.01显著性水平的检验。2000、2010、2018三年间各生态系统服务间的权衡与协同关系变化不大, 稳定在一定数值。其中产水服务、土壤保持服务与其余所有生态系统服务均表现为协同关系。产水量高的区域植被普遍长势较好, 保土能力强, 研究结果表明土壤保持服务和产水量服务之间有很强的协同关系, 三个年份的相关系数保持在0.616-0.683之间。水和土壤是植被生长发育的基础, 各年份中固碳服务与产水量服务、固碳服务与土壤保持服务、生境质量与产水量服务间也表现出较强的协同关系。土壤保持与生境质量服务之间虽然存在协同关系, 但相关系数较低, 土壤保持服务和生境质量之间的相关系数多年间仅维持在0.11左右。

3.3 那曲市生态系统服务空间分异的驱动因子识别 3.3.1 生态系统服务空间异质性的单因子归因

那曲市各项生态系统服务均表现出了明显的空间分异性, 而这种分异性多是在自然因素和人类活动等的共同作用下形成的。本研究借助地理探测器的“因子探测”和“交互探测”模块得到了那曲市不同年份各驱动因子对不同生态系统服务空间异质性的贡献率(图 5)。由图 5可知, 各驱动因子对那曲市生态系统服务的解释力度在2000年、2010年、2018年间有较强的一致性, 故本研究后续以2018年数据为例进行那曲市生态系统服务影响因素的归因分析。

那曲市产水服务的空间分异受到多因子的共同作用, 各因子按照q统计量的大小排序分别为降水(0.8785)>土壤类型(0.5910)>NDVI(0.5101)>温度(0.2465)>坡度(0.1931)>土地利用类型(0.1828)>人口密度(0.1460)>高程(0.1212)>SHDI(0.0132)。因子探测结果表明, 降水对那曲市产水服务空间分异的解释力度最大, q值达到了0.8785, 那曲市产水服务的空间分布也大致与降水的空间分布一直, 表现出东南高、西北低的特点, 年降水量处于749.89-880.93 mm之间的区域所带来的产水量最高(图 6), 平均产水量可达483.96 mm。同时, 土壤类型、NDVI等自然因素对那曲市产水量服务的影响也极为显著, 那曲东部发育着大面积的草毡土和黑毡土而中西部和北部则多为寒钙土, 草毡土和黑毡土相较于寒钙土有更强的肥力, 且土层更厚, 植被发育更好。草毡土和黑毡土覆盖区的平均产水深度分别为314.92 mm、473.35 mm, 平均产水深度最低的寒钙土仅为33.92 mm。那曲市东部地区由于水热状况优越, 植被长势好使得NDVI也明显高于中西部与北部区域。NDVI处于0.73-0.99之间的区域产水量最大, 平均产水深度为406.99 mm。那曲市人类活动因素对产水服务的影响有限, q值处在一个较低的水平。

各影响因子对那曲市土壤保持服务空间分异的解释力度如图 5所示, 按照q统计量排序分别为坡度(0.5954)>土壤类型(0.2813)>降水(0.2634)>NDVI(0.1585)>高程(0.0748)>温度(0.0720)>土地利用类型(0.0653)>人口密度(0.0395)>SHDI(0.0059)。坡度是对那曲市土壤保持服务影响最为显著的因子, q统计量达到0.5954, 土壤保持服务功能较强的区域主要集中在那曲市东南部的高山峡谷与河谷等坡度变化较大的区域(图 3), 随着坡度的增大, 土壤保持量也呈现出一个增长的趋势(图 6)。土壤类型、降水以及NDVI对土壤保持服务的空间分布也起到了较为重要的作用。不同土壤类型中, 黑毡土覆盖区域的土壤保持能力最强, 单位面积土壤保持量为356.22 t/hm²。土壤保持服务高值区与降水、NDVI的高值区也基本吻合, 降水量处在749.89-880.93 mm之间、NDVI处于0.73-0.99之间的区域土壤保持服务也比较强, 单位面积土壤保持量分别为213.27 t/hm²、215.99 t/hm²。人类活动及土地利用因素对那曲市土壤保持服务空间分布的影响较小(图 5)。

各驱动因子对那曲市固碳服务空间异质性的解释力度比较高(图 5), 按照q统计量排序分别为NDVI(0.6999)>土壤类型(0.4599)>降水(0.3992)>人口密度(0.2441)>温度(0.2288)>高程(0.1586)>坡度(0.0558)>土地利用类型(0.0390)>SHDI(0.0234)。分析结果表明, NDVI对那曲市固碳服务空间分布的贡献率最高, q统计量达到了0.6999。NDVI是指归一化植被指数, 其高值区普遍为植被覆盖水平较高的区域, 故固碳量高, NDVI处于0.73-0.99之间的区域固碳量最大, 平均固碳量达到了585.36 gC/m²。土壤类型因子方面, 那曲市固碳服务高的区域发育有大面积的黑毡土、沼泽土、草毡土, 其中黑毡土上发育着那曲市仅有的林地, 沼泽土、草毡土上发育有大面积的优质牧草。那曲市黑毡土覆盖区域单位面积的固碳量最高(图 6), 可达到530.42 gC/m²。那曲市降水量高的地区固碳服务也比较强, 研究区降雨量处于642.67-749.89 mm之间的地区(图 6), 其单位面积的固碳量最大, 达到了477.01 gC/m²。

由图 5可知, 各因子对那曲市生境质量空间分异的解释度按照q统计量排序分别为土地利用类型(0.3356)>NDVI(0.1757)>土壤类型(0.1350)>人口密度(0.1002)>SHDI(0.0374)>降水(0.0233)>高程(0.0168)>温度(0.0083)>地形(0.0014)。土地利用类型是对研究区生境质量空间分布贡献率最大的影响因子, 其次为NDVI和土壤类型。那曲市生境质量的空间异质性明显, 主要区域即湖泊、河流及林草地覆盖区域的生境质量要明显优于其他区域, 生境质量得分较高, 而人口最密集、生产建设活动剧烈的建设用地区域生境质量低。NDVI值越高的地区生境质量的平均得分也越高, NDVI值处于0.73-0.99之间的区域, 其生境质量平均得分达到了0.54。除水域外, 那曲市生境质量最高的区域发育的土壤主要为黑毡土(图 6), 黑毡土覆盖区域的生境质量平均得分为0.62。土地利用类型、NDVI、土壤类型等因子的解释力度明显高于高程和地形因子, 地形与海拔对那曲市生境质量的影响有限。

3.3.2 生态系统服务空间异质性影响因子的交互作用探测

交互探测的结果用于表明任意两因子的组合是否会增强或减弱它们各自对生态系统服务空间分布的解释力度。各生态系统服务的因子交互探测结果显示(表 3), 任意两个因子的交互作用均大于单一因子对各生态系统服务空间分布的解释力度, 因子间两两交互的类型主要表现为“非线性增强”和“交互增强”, 这样的交互作用可以表明那曲市各项生态系统服务的空间格局是由多个因子共同作用的结果。

由表 3可知, 在产水服务中降水量与其他各因子的交互作用对那曲市产水服务空间格局的解释力度最大, 任意两因子组合的q值均大于单一降水因子的q值(0.8785), 其中降水∩土地利用类型(降水因子与土地利用类型因子的交互作用)的q值更是达到了0.9551, 说明在同一降水水平下, 土地利用类型的不同会对产水服务的空间格局造成很大的影响。NDVI∩其他因子以及土壤类型∩其他因子对产水服务空间格局的影响力也比较强, q值处在0.5101-0.8984之间。

对土壤保持服务空间格局的交互探测结果表明(表 3), 地形因子中的坡度与其他因子的交互作用对土壤保持服务的空间格局影响最大, 其中坡度∩NDVI、土壤类型等自然因子对土壤保持服务解释力度达到了0.7以上, 而坡度∩人口密度的q值仅比单一坡度的q值高了0.0082, 也可说明人为因素对那曲市土壤保持服务空间格局的影响程度有限。降水、土壤类型与其他因子交互对土壤保持服务空间格局的解释力度也比较大, q值处在0.2716-0.7248之间。

由表 3可知, 土地利用类型与其他因子交互作用对那曲市生境质量空间格局的解释力度最大, 其次是NDVI∩其他因子以及土壤类型∩其他因子。土地利用类型∩NDVI以及土地利用类型∩土壤类型的q值分别为0.3895和0.3783, 其次为土地利用类型∩人口密度和土地利用类型∩SHDI。

在固碳服务的因子交互作用中(表 3), NDVI与其他因子的交互作用对那曲市固碳服务空间格局的影响程度最大, q值均在0.7以上, 其次为土壤类型∩其他因子以及降水量∩其他因子。单因子的探测结果显示土地利用类型的q值要远小于土壤类型和降水, 然而交互作用中, NDVI∩土地利用类型的q值高于NDVI∩降水与NDVI∩土壤类型, 达到了0.7915, 表明在主导因素NDVI一致的情况下, 土地利用类型的变化成为影响那曲市固碳服务空间格局的重要因素。

4 讨论与结论 4.1 讨论

量化区域生态系统服务的空间异质性并揭示异质性形成的主要驱动因子可帮助决策者有针对性的施行生态优化和调控措施。然而近些年的相关研究多是以回归和相关为手段, 没能明确指出哪些驱动因子及因子组合会对生态系统服务空间异质性造成主要影响, 而地理探测器在这一方面拥有独特优势。本文尝试在模型模拟生态系统服务时空格局的基础上, 用地理探测器探究那曲市各项生态系统服务空间异质性的主要驱动因子, 较好的揭示了那曲市各项生态系统服务的主控因子和因子组合, 可为类似研究提供参考。本研究的模型模拟结果中那曲市各生态系统服务的值与临近研究区已有的研究结果类似^{[36, 44-47]}。如, 潘韬等^[45]的研究表明三江源2000年的平均产水深度为247.1 mm；杨蕾^[44]在三江源的研究表明三江源多年平均土壤保持量为72.03 t/hm²；陈童尧等^[46]在祁连山保护区的研究显示祁连山区2015年的单位面积土壤保持量为212.63 t/hm², 该区域环境本底特征与那曲东南部县区类似；邹德富^[47]的研究表明青藏高原2003-2010年的NPP均值为298 gC/m²。同时, 上述研究中各项生态系统服务的空间分布趋势也与本研究的模拟结果存在连续性, 故模型模拟结果可信, 能够用于生态系统服务空间异质性归因分析研究。

生态系统服务的空间异质性是在众多影响因子的共同作用下形成的。降水、土壤类型、NDVI这类自然因子对那曲市各项生态系统服务的空间异质性都起到了重要作用, 该结论与类似区域已有的研究结果一致^{[43, 48-49]}。在这些主要影响因子的共同作用下, 那曲市各项生态系统服务的空间格局存在相似性。除了上述不易受人为因素直接影响的自然因子外, 土地利用类型可以在很大程度上解释生境质量的空间格局, 从交互探测的结果来看, 土地利用类型分别对那曲市固碳服务和产水服务空间异质性起到重要影响, Hopping等^[50]和Tang等^[51]的研究也表明青藏高原土地利用及土地利用政策变化会对地区生态系统服务产生重大影响。在国家对口援藏政策^[52]的持续支持下, 近年来那曲市城镇聚集效应明显提升, 城镇基础设施建设也在逐步完善, 伴随而来的将是大量的其他用地类型向建设用地的转变。2018年那曲市建设用地相较于2000年增长了近95%, 一段时间内那曲市城市、农村居民点及工矿用地等建设用地面积还将持续增长, 而这势必会对当地的生态系统服务产生影响。对此, 建议那曲市应当：1)在未来城市建设与经济发展过程中继续将生态保护摆在优先位置, 继续推进生态安全屏障建设；2)严格落实生态保护红线, 对土地利用进行合理规划, 合理统筹生态用地、城镇建设用地、工矿用地、农业用地等功能空间；3)提升对全市各区域天然林、草地、湿地的保护力度, 同时对已受建设及工矿活动影响的生态脆弱区采取划区轮牧、退牧育草、大面积封禁管护等措施进行恢复重建。

人类活动本身会对生态系统服务的空间分布带来很大影响^[53], 但截至目前人口密度因子对各生态系统服务空间异质性的解释力度并不高, Li等^[54]研究同样表明在1990-2010年间青藏高原整体的生态环境受人类活动的影响较低。那曲市幅员辽阔但人口稀少, 虽然近年来那曲市城镇集聚效应在逐步提升但城镇化整体仍处在起步阶段, 建设用地占比较低。同时, 随着近年来陆续施行的限牧禁牧、退牧还草、生态补奖等^[55-56]政策, 放牧活动对草地的破坏也得到了一定的限制。这使得在那曲全域尺度上人类活动的规模和强度仍相对较小, 对生态系统服务供给的影响尚处在较低水平。那曲大多数县区刚刚实现脱贫摘帽, 但大量居民尤其是牧民的生活水平仍然处在贫困线附近。今后, 那曲的经济发展及城镇化建设力度仍会增强, 生态保护和经济发展、民生改善之间的矛盾可能会更加突出。对此, 建议那曲市应当：1)在带领藏区群众走向全面小康的同时, 仍需加大力度落实生态扶贫及生态补偿政策, 在生态系统服务重要区域内按照草蓄平衡原则, 实施划区轮牧、阶段性禁牧和季节性休牧, 在高海拔生态脆弱区积极推进易地生态搬迁, 在提高民生的同时将人类活动对生态环境的影响降至最低程度, 努力做到生态环境与民生改善的协调发展；2)增强群众的环保意识, 对全市居民及外地游客广泛宣传生态环境相关的法律法规以及环保知识, 呼吁群众爱护自己赖以生存的生态环境；3)加强对各区域城镇及农村环境的综合治理, 实现城镇及广大农村区域生活垃圾的收运处置, 减少各类生活垃圾流通到环境中造成的环境污染；4)打造高质量生态经济产业, 推动生态旅游业等特色产业的发展, 在保护生态环境的同时获取经济利益, 提升附近居民生活水平。

本文所采用的研究方法中存在有一定的不确定性。InVEST模型和CASA模型凭借其操作简单、参数输入相对较少、采用空间数据作为模型的输入与输出并且模型的输出结果可直接应用于分析生态系统服务的时空分布格局的优势, 获得了越来越多科研人员的使用和掌握。本研究运用InVEST模型和CASA模型进行生态系统服务的模拟评估, 在各模块的生物物理参数的计算与选择过程中均结合研究区特点参考了已有的计算方法和相关文献资料, 但在参数的确定中仍然存在着不确定性, 影响模拟精度。原因在于各参数(如产水模块中的植被根系深度、土壤保持模块中的水土保持措施因子、生境质量模块中各地类对威胁源的敏感性等)在基于研究区特点进行本地化修订时还需要针对研究区特点结合大量的野外观测调研及实验分析, 从而验证参数并确定参数应用的合理性, 进而增强模型的模拟精度。本文使用地理探测器分析生态系统服务空间异质性的主要驱动因子, 探测结果可以定量化的表现各驱动因子和因子组合对生态系统服务空间异质性的解释力度, 但此方法使用过程中在对连续的自变量因子进行离散化处理时没有明确的分类标准, 离散化方式的不同会对最终的探测结果产生一定影响, 日后应当针对此方面加强分析, 进一步提升地理探测器方法的可靠性。在生态系统服务空间异质性的评估过程中, 本研究目前仅在全域尺度上分析了那曲市各生态系统服务空间异质性的主要影响因素。然而在不同的空间尺度上同一生态系统服务的主要影响因子可能会存在差异, 今后的研究中应当着力在更小的空间尺度上进行生态系统服务空间异质性的归因分析, 为不同区域制定区别化的生态保护政策提供依据。

4.2 结论

本文以模型模拟的方式评估了那曲市2000-2018年4项生态系统服务, 分析了它们的时空格局及相互关系, 随后借助地理探测器揭示了生态系统服务空间异质性的主要影响因子及因子组合, 得到如下结论：1)2000-2018年间, 那曲市产水服务总体呈下降趋势, 土壤保持服务表现为先降低后增高, 生境质量服务趋于稳定, 固碳服务呈现上升趋势。那曲市西北部各县区生态系统服务年际变化趋势相近, 东南部各县区生态系统服务的变化也存在一致性。2)2000-2018年间产水、土壤保持、固碳服务的高值区主要集中在那曲市东部与中部水热与植被条件好的区域, 尤其是高山峡谷与河谷地区, 生境质量的高值区分布在湖泊河流以及植被覆盖率高的地区。那曲市生态系统服务极重要区域与高度重要区域即各生态系统服务高值重叠区域则主要位于那曲市东南部高山峡谷与河谷区及其临近地区, 为那曲市整体生态环境最优的区域。3)2000-2018年间那曲市生态系统服务除固碳服务和生境质量服务之间不存在明显关系以外, 其余组合均表现为相互促进的协同关系。4)降水、土壤类型、NDVI等自然因子是对产水服务和固碳服务空间异质性解释力最强的驱动因子, 坡度、土壤类型、降水等因子对土壤保持服务空间异质性的解释力度强, 生境质量服务的空间异质性主要受土地利用类型、NDVI、土壤类型等的影响, 各生态系统服务空间异质性的主要驱动因子存在相似性。此外, 在交互作用中土地利用类型的重要性凸显, 土地利用类型∩NDVI对固碳服务和生境质量有最大的解释力度, 土地利用∩降水可在最大程度上解释产水服务的空间分异。截止目前代表人类活动的人口密度因子对各生态系统服务空间异质性的解释力度尚处在一个较低的水平。