生态系统服务(Ecosystem Services, ES)是指人类为满足自身需求从生态系统中所获取的直接或间接收益^[1-2], 主要包含供给服务、调节服务、支持服务和文化服务等4种服务类型^[3]。自1997年Constanza等率先估算全球尺度生态系统服务以来, 国内外学者陆续对全球^[4]、国家^[5]、区域^[6]、省市^[7]、山区^[8]、流域^[9]等不同空间尺度生态系统服务进行评估, 并测度冰川^[10]、森林^[11]、草地^[12]和湿地^[13]等类型生态系统服务。不同生态系统服务受到自然环境和人为干扰的综合影响呈现此消彼长的权衡关系和互相增益的协同关系^[14-16], 正确认知生态系统服务相互关系是制定多种生态系统服务可持续管理决策的前提, 有助于区域生态系统服务效益最优化, 从而实现经济发展与生态环境保护双赢^[17]。权衡全球或区域尺度上生态系统服务相互关系成为学界研究的热点^[18], 统计描述法、空间统计制图法、情景分析法和模型模拟法是权衡与协同常用分析方法^[19-20], 如王鹏涛运用二阶偏相关分析逐像元计算土壤保持、NPP和产水的相关系数, 研究表明3种服务皆以权衡关系为主导^[21]；钱彩云利用双变量空间自相关测度白龙江流域产水量、土壤保持、碳储量和粮食生产权衡程度, 高高集聚和低低集聚为协同关系、高低集聚和低高集聚为权衡关系, 识别空间尺度上生态系统服务关联特征^[22]。以往研究多集中于测度时间尺度上生态系统服务两两之间权衡关系并进行空间可视化, 生态系统服务整体间权衡程度研究相对薄弱^[23-24]。在此情形下, 亟需评估区域层面上生态系统整体权衡关系及空间特征, 对生态系统服务两两之间非线性权衡关系研究进行补充, 揭示生态系统服务权衡程度及空间差异, 为制定流域生态环境保护政策提供科学依据。因此, 本文以淇河流域为对象研究生态系统服务权衡及空间分异机制, 探讨小流域生态系统服务类型之间的权衡与协同规律。

自然和人文因素综合影响生态系统服务的形成及其权衡关系, 探究生态系统服务关系的外部驱动机制受到学界的广泛关注^[25], 如张宇硕通过层次分析法构建结构方程模型分析京津冀湿地生态系统服务的影响因素, 结果得出气候、土壤等内部驱动力对生态系统服务起促进作用, 外部驱动力抑制生态系统服务提升^[26]；郭慧基于生态系统服务价值的均衡因子和产量因子测度方法, 评估了1990—2015年北京市门头沟区生态承载力, 发现均衡因子和产量因子较通用因子更好体现区域生态承载力^[27]。上述研究大多通过分析外部驱动因素对各项生态系统服务的影响程度, 但对于生态系统服务关系对环境变量的响应机制亟需厘清^[28], 已有研究表明植被覆盖度^[29]、土地利用类型^[30]、地形^[31]、人类活动影响^[32]对生态系统服务之间相互关系的作用方向及影响程度存在差别, 研究手段多为统计学方法, 对外部驱动因素在空间尺度上影响生态系统服务动态权衡关系解释程度不够^[33]。地理探测器是一种新型的空间统计方法, 能有效识别同一地域单元上相邻要素之间空间关联程度^[34-35], 主要用于测度驱动因子影响程度及其交互作用, 已广泛应用于城市建设用地驱动^[36]、生态系统服务价值空间分异^[37]、粮食安全^[38]、疾病监测^[39]、贫困化差异机制^[40]等方面。流域是具有明显物理边界且综合性极强的地理单元, 运用地理探测器分析流域尺度生态系统服务权衡关系的外部驱动因素, 有助于识别生态系统服务权衡关系影响因素的地理影响及交互作用, 能够为流域的生态环境规划及决策提供科学依据和实证分析。

淇河流域位于中国第二级阶梯向第三级阶梯过渡区域, 地处太行山与华北平原的地理过渡地带, 具备丰富的水源涵养功能及供给服务能力^[41]。本文以淇河流域为研究对象, 利用CASA模型、InVEST模型测度2000—2015年流域内固碳、土壤保持、水源供给等生态系统服务, 刻画流域尺度生态系统服务时空变化特征, 运用空间叠置法度量不同时段生态系统服务空间权衡程度, 进一步利用相关性分析精准分析生态系统服务两两之间权衡关系, 结合地理探测器中因子探测和交互探测模块分析生态系统服务权衡关系外部驱动机制, 从而量化自然环境和人为扰动对生态系统服务权衡的影响程度, 为改善区域生态环境质量及提升人类福祉提供理论依据和实证分析。

1 研究区概况

淇河流域(35°32′—36°15′N, 113°15′—114°23′E)地处华北平原西南部和太行山南段之间, 源于山西省陵川县方脑岭, 经河南省汇入海河支流卫河, 流域面积2227 km², 呈现西高东低分布态势(图 1)。研究区主要气候类型为暖温带半湿润大陆性季风气候, 年均降水量574 mm, 年均气温11.9℃, 夏季降水集中；复杂的地形使淇河流域生态环境脆弱, 多样的地貌类型使流域生态系统服务呈现复杂的变化。

2 数据来源与研究方法 2.1 数据来源与处理

2000、2005和2010年淇河流域土地覆被/利用数据均来自中国地球系统数据共享平台-黄河中下游科学数据中心(http://www.geodata.cn/), 2015年土地利用数据在参考研究区历年土地利用类型图的基础上对LANDSAT5多波段遥感影像(来源于地理空间数据云, http://www.gscloud.cn/)进行人机交互目视解译并进行实地勘察, 其Kappa系数为86%；研究区内9个气象站点的气温和降水数据来于河南省气象局及中国气象科学数据共享网(http://cdc.cma.gov.cn)；土壤数据来源于国家地球系统数据共享平台-黄河下游科学数据中心(http://henu.geodata.cn)；植被NDVI数据来源于美国国家航空航天局(NASA) 16 d合成的空间分辨率为250 m MOD13Q1产品, 利用最大合成法提取每年7月中NDVI的最大值, 从而合成逐年NDVI数据；DEM数据来源于地理空间数据云平台(http://www.gscloud.cn), 分辨率为30m；土地利用强度和人类活动影响通过文献计算得到^[42-43]。

2.2 生态系统服务评估方法

本文利用CASA模型和InVEST模型等生态系统服务评估模型, 定量估算淇河流域固碳、土壤保持与水源供给服务时空变化特征。固碳利用CASA(Carnegie Ames-Stanford Approach)生态系统碳循环过程模型估算生态系统净初级生产力(Net Primary Productivity, NPP), 参数确定及模型计算方法, 具体见参考文献^[44]；土壤保持服务利用InVEST模型中SDR(Sediment Delivery Ratio Model)模块进行计算, 结合土地利用/覆被数据、DEM、土壤数据和降水等多种数据对土壤流失方程进行修正, 模型中所需的降雨侵蚀力因子(R)是利用Kriging插值得到并运用Wischmeier的月尺度经验模型进行计算得到^[45], 土壤可蚀性因子(K)是通过Williams等提出的EPIC模型对流域内1 ∶ 100万土壤数据集进行计算, 植被覆盖因子(C)参考蔡崇法的计算方法得到^[46], 水土保持措施因子(P)参考以往的研究并结合研究区土地利用情形及农业生产活动确定^[47-48]；水源供给服务根据InVEST模型中产水量模块计算得到, 该模块是根据水量平衡原理并通过Budyko曲线^[49]和年降雨量计算得到, 相关参数通过文献^[50-51]进行修正。

2.3 生态系统服务权衡与协同量化方法

基于Matlab平台对固碳、土壤保持和水源供给服务逐像元进行相关性分析, 量化固碳与土壤保持、固碳与水源供给及土壤保持和水源供给服务间的权衡协同程度并进行可视化表达。进一步利用t检验法对生态系统服务权衡强度进行显著性分析, 若P < 0.01, 则权衡/协同度通过了极显著性检验, 若0.01 < P < 0.05, 表明通过了显著性检验^[52]。为探究研究时段内不同时期生态系统服务权衡协同程度, 利用ArcGIS 10.3软件中的地图代数工具, 将生态系统服务标准化处理后并运用自然断裂点法进行分级, 按照低、中、高三类并进行重分类, 依次对应1、2、3, 参考叠置规则对标准化的各项生态系统服务值进行叠置^[53], 计算公式如下：

(1)

式中：ES_std为生态系统服务标准化的值(介于0—1), ES_pixel、ES_max、ES_min分别为某类型生态系统服务任意栅格值、最大值和最小值, 按照自然断裂点法将服务标准化的值分类为高、中、低；CODE为三位数叠置代码, 碳为固碳、水为水源供给、土为土壤保持, 代码序列依次为固碳、水源供给和土壤保持, 值介于111和333之间, 分别对应服务类型的供给能力, 即低、中、高。参考张静静等^[14]研究将一种生态系统服务能力高而其它生态系统服务能力较低视为强权衡, 两种生态系统服务供给能力高而另一种生态系统服务能力低定义为弱权衡；高协同表现为各项生态系统服务都处于较高水平、如供给能力组合中333、332等, 低协同是各项生态系统服务都处于较低水平, 如111、221等, 具体分级标准见表 1：

2.4 地理探测器

地理探测器能有效分析同一地域单元上两个要素之间的地理关联性以及不同区域变量的差异性, 可以揭示区域尺度上要素之间的因果性^[54-55]。本文利用地理探测器中因子探测和交互探测两个模块, 其中环境变量通过自然断裂点法进行分类, 以此探测各影响因素对生态系统服务权衡关系的影响及交互作用影响程度, 公式如下：

(3)

式中：q为探测因子A的探测力值, q∈[0, 1], q值越大, 表明因素A对于研究区生态系统服务权衡/协同影响度越高；N和N_h分别为研究区样本数；δ_h²为因素A在样本h内的离散方差；L为研究区各因素的类型。交互探测类型可以划分为双因素增强、非线性减弱、非线性增强、独立和单因子非线性减弱等五类, 识别自然环境和人为扰动交互作用后对生态系统服务关系综合影响程度^[55]。

3 结果分析 3.1 淇河流域生态系统服务时空变化特征

淇河流域固碳、土壤保持与水源供给空间分布及变化趋势(图 2)。近15年来研究区平均固碳、土壤保持及水源供给量分别为375.55 gC m^-2 a^-1、396.72 t hm^-2 a^-1和67.26 mm/a, 且各项生态系统服务空间分布上存在差别：固碳含量表现为西北高、东南低, 高值区集中分布淇河上游地区, 低值区分布于下游区域；土壤保持量高值区集聚态势不显著, 流域中上游、中下游分布较广；水源供给量高值区主要分布于流域四至及中下游区域, 低值区分布范围较广。2000—2015年间固碳含量以18.20 gC m^-2 a^-1的速率递减, 超过90%的区域呈减少态势, 淇河上游及下游地区减幅最大；土壤保持量递减速率为1.2 t hm^-2 a^-1, 减少幅度主要集中于[-8, -0], 流域内土壤保持含量均以一定速率递减；水源供给服务递减率为0.16mm/a, 95%的区域产水量减少, 流域中下游区域产水量增加较为明显。

3.2 生态系统服务权衡与协同关系定量分析 3.2.1 生态系统服务权衡分析

由图 3和表 3可知, 2000—2015年淇河流域生态系统服务关系以低协同和强权衡关系为主, 其中低协同主要分布于流域上游和下游, 固碳、土壤保持和水源供给等生态系统服务整体较低；研究区中下游生态系统服务关系以强权衡为主导, 表示一类生态系统服务供给能力较高而其它类型较低；从各时期生态系统服务权衡与协同关系像元个数及所占比例来看, 强权衡仅在2000年像元个数占比低于50%, 其它时期均高于50%, 表明淇河流域生态系统服务以强权衡为主导；在弱权衡方面, 2005年像元数量最多(1519个)、2015年次之(813个)、2000年最少(451个)；高协同是生态系统服务都处于较高水平, 生态系统服务关系达到最理想的状态, 2015年占比最多(10.74%)、2005年次之(10.63%)、2000年最少(7.43%)；低协同在2000年所占比例最高(50.64%)、2010年次之(37.20%)、2005年最低(31.69%), 反映固碳、土壤保持和水源供给服务都处于较低水平, 也是生态系统服务关系中最不理想的状态。研究时段内, 低协同分布的面积减少了16.21%, 强权衡、高协同和弱权衡分别增加了11.85%、3.31%和1.05%, 表明流域内生态系统服务协同程度加强, 生态环境呈好转趋势。

3.2.2 生态系统服务两两权衡分析

基于像元尺度计算2000—2015年固碳、土壤保持和水源供给相关性并进行空间制图, 图 4表示固碳与水源供给、固碳与土壤保持及水源供给和土壤保持服务之间的权衡协同程度及显著性检验分析。固碳与水源供给服务关系以权衡为主导, 比重为65.2%, 且通过显著性水平检验仅占3.04%, 表明两者生态系统服务权衡程度处于较低状态、显著性水平不高, 时间尺度上固碳服务的增加/减少对水源供给服务的减弱/增强的影响较为平缓, 协同关系比重为34.8%, 且通过显著性检验为2.63%；固碳与土壤保持权衡、协同占比分别为2.63%、97.37%, 主导服务类型为协同关系, 其中协同度大于0.8的像元占比8.08%、0.5—0.8占比64.76%, 即研究时段内固碳与土壤保持服务变化方向一致、协同程度高, 通过显著性水平的像元比重为10.56%, 其中约有8.54%像元表现极显著相关、1.68%像元为显著相关; 水源供给和土壤保持服务关系表明, 权衡关系比重为90.06%, 通过显著性检验的像元数为5969, 占比18.41%, 反映淇河流域水源供给服务和土壤保持服务空间尺度上表现为此消彼长的权衡关系, 即像元尺度上水源供给服务能力的增加会减损土壤保持服务供给量。

3.3 淇河流域生态系统服务权衡驱动力探测 3.3.1 因子探测

运用地理探测器中因子探测模块对生态系统服务权衡关系外部驱动因子进行分析, 揭示自然环境因子和人类活动影响空间尺度上对权衡关系的影响机制(表 4), 发现生态系统服务两两权衡关系对不同环境变量响应程度存在差别：在固碳与土壤保持权衡方面, 植被NDVI影响最强(0.306)、气温次之(0.247)、土地利用强度最弱(0.002)；固碳与水源供给方面, 植被NDVI影响强度远大于其它环境变量, 证明固碳与水源供给权衡关系主控因子为植被NDVI；水源供给与土壤保持服务外部驱动机制表明, 降水量因子解释力最强、为0.594, 气温次之(0.574), 土地利用程度影响能力最弱(0.004), 发现两类生态系统服务关系受自然环境因子影响显著高于人类活动及植被NDVI, 主要原因为地形及气候条件决定了水源供给及土壤保持服务供给强度, 显著影响生态系统服务相互关系。

3.3.2 交互探测结果分析

生态系统服务权衡的影响过程受到外部因素多重影响, 各影响因素彼此存在交互关系。由交互探测结果(表 5)可以看出, 两两因素交互探测后对各生态系统服务权衡关系因子解释力显著增大, 从不同交互类型来看, 主要表现为双因子增强和非线性增强两种类型：在固碳与土壤保持服务中, 植被NDVI与气温、降水和高程交互后因子解释力分别达到0.4以上, 植被NDVI、气温、降水和高程相互作用后对生态系统服务关系解释程度增强, 表明生态系统服务关系受到多种因素共同作用, 而非单一因子起决定性作用；在水源供给与土壤保持服务中, 气温、降水和高程与其它环境因素交互后解释力明显增强, 反映气温、降水和高程在水源供给服务与土壤保持服务关系中处于核心位置, 解释程度远高于其它环境变量。

4 结论与讨论 4.1 讨论 4.1.1 生态系统服务权衡与协同

目前, 生态系统服务权衡研究多集中于生态系统服务两两之间相互关系的时空权衡分析及空间表达上, 缺乏对生态系统服务整体权衡机制的深入剖析。本文定量评估了淇河流域固碳、土壤保持和水源供给服务, 生态系统服务整体关系以强权衡和低协同为主导, 其中低协同是生态系统服务关系都处于较低水平, 也是生态系统服务最不理想的状态, 迫切需要政府部门制定有关政策进行调控, 促使低协同区域逐渐过渡至高协同。运用空间叠置法对于研究区域尺度上生态系统服务整体权衡关系有积极意义, 但仍需进一步研究生态系统服务关系演变趋势及局部特征。相关性结果表明, 固碳与水源供给、水源供给与土壤保持以权衡为主导, 固碳与土壤保持以协同关系为主, 且空间格局上存在差异, 其原因为约束阈值对生态系统服务关系产生一定的影响^[56], 如固碳较高的区域植被覆盖水平较高, 抑制了土壤侵蚀作用, 这时固碳与土壤保持服务呈协同关系, 而植被覆盖度较高的区域降水量较为充沛, 土壤侵蚀程度加剧, 固碳与土壤保持服务表现为权衡关系^[21]。另外, 生态系统服务权衡关系存在尺度特征, 潘竞虎发现县域尺度上各项生态系统服务表现协同关系, 而区域尺度上土壤保持与水源供给、粮食生产为权衡关系^[57], 未来应加强流域生态系统服务权衡关系尺度特征与区域差异的研究, 构建栅格-乡镇-县域-流域多级评价单元, 明晰各尺度上生态系统服务能力及权衡关系。本研究集中评估固碳、土壤保持和水源供给等调节服务权衡与协同关系, 后期应开展其它类型生态系统服务的估测并权衡其相互关系, 从而建立流域生态系统服务簇并刻画各项服务权衡关系, 为流域国土空间规划及优化提供参考。

4.1.2 生态系统服务权衡驱动机制探讨

本文运用遥感、GIS和地理探测器技术, 系统分析了淇河流域生态系统服务权衡关系及其空间分异的驱动机制, 地理探测器有助于揭示淇河流域生态系统服务权衡程度空间分异的主要外部驱动因子及其交互作用, 较以往的统计学方法, 优势在于能够定量描述各环境变量对生态系统服务权衡关系的空间影响, 对于研究生态系统服务关系的驱动机制是一种有益的尝试, 为流域生态系统服务管制及生态安全格局的建立提供依据。综合来看, 自然环境因素对生态系统服务关系影响强度显著高于人类活动影响, 植被NDVI、气温、降水和高程是影响淇河流域生态系统服务权衡关系的主导因子, 且交互作用后因子解释力明显增大, 表明流域内生态系统服务关系受到多种因素共同影响, 而非单一因素起决定性作用。在CASA和InVEST模型参数输入过程中需要输入诸如植被NDVI、地形等要素, 一定程度上增加了自然要素的影响程度, 此外, 自然环境决定生态系统服务供给强度并形成生态系统服务权衡与协同关系, 这也可能对研究结果产生一定的影响。有研究表明, 土地利用/覆被变化是影响生态系统服务权衡与协同关系的主导因子^[58-59], 本文通过计算土地利用强度和人类活动影响表征人为扰动程度, 结果显示两者对生态系统服务关系影响程度较低, 主要是生态用地提供生态系统服务, 但土地利用强度及人为影响较弱, 其次, 数据分辨率较低也可能是影响地理探测结果的因素之一。因此, 今后研究的重点应包含：(1)提高土地利用/覆被及DEM数据的分辨率, 准确评估流域尺度生态系统服务能力；(2)选取自然环境-社会经济-人类活动等多种影响因素, 探究其对生态系统服务关系的影响机理；(3)基于土地利用/覆被变化, 进一步测度生态系统服务权衡与协同关系的变化机制；(4)关注极端气候变化对生态系统服务及其权衡关系的影响, 建立极端气候与生态系统服务时空匹配关系, 揭示全球变化背景下生态系统服务关系变化情形, 为区域生态环境建设、构建和谐的人地关系提供理论依据。

4.2 结论

(1) 2000—2015年, 淇河流域多年平均固碳量为375.55 gC m^-2 a^-1, 呈西北高、东南低分布趋势, 多年递减速率为18.20 gC m^-2 a^-1；多年平均土壤保持量为396.72 t hm² a^-1, 多年递减速率为1.2 t hm^-2 a^-1, 空间差异不显著；多年平均水源供给量为67.26 mm/a, 高值区主要分布于流域四至及中下游区域, 低值区分布广泛, 以0.16 mm/a的速率递减。

(2) 在生态系统服务权衡与协同关系中, 以强权衡和低协同关系为主, 其中强权衡主要分布于流域中下游、低协同主要分布于研究区下游及中上游：研究时段内低协同分布的面积减少了16.21%, 强权衡、高协同和弱权衡分别增加了11.85%、3.31%和1.05%。生态系统服务两两之间关系表明, 固碳与水源供给、水源供给与土壤保持服务以权衡关系为主, 水源供给与土壤保持服务以协同关系为主, 空间特征及显著性水平表现为异质性。

(3) 地理探测器结果表明：固碳与土壤保持方面, 植被NDVI影响最强(0.306)、气温次之(0.241)、土地利用程度最弱(0.002)；固碳与水源供给方面, 植被NDVI因子解释力显著高于其它因素, 为0.381；水源供给与土壤保持方面, 降水、气温和高程因子解释力均高于0.5；交互探测后各影响因子解释能力显著增强, 主要为双因子增强和非线性增强两种类型。