生态系统服务是指生态系统所形成和维持的人类赖以生存的自然环境条件与效用^[1], 是人类直接或间接地从生态系统中获取的所有收益^[2]。生态系统服务类型多样, 包括维持生态系统动态平衡的养分循环、土壤形成、气候调节、净化水源等支持服务和调节服务, 维持人类生存与发展的粮食、淡水、纤维等供给服务, 以及愉悦人类精神的美学、教育、宗教等文化服务^[3-4]。2005年, 联合国发布的千年生态系统评估(MA)报告指出, 其评估的24项生态服务中约有15项正处于退化或不可持续利用状态, 直接威胁着区域乃至全球的生态安全^[3], 此后生态系统服务逐渐成为国际生态学和相关学科研究的前沿和热点^[4-5]。由于生态系统服务的种类多样性、空间分布异质性, 以及人们对生态系统服务的需求偏好, 导致生态系统服务之间形成了相互交织、复杂的非线性关系^[6-7], 已有的研究通常将生态系统服务之间的关系抽象概括为此消彼长的权衡、相互增益的协同及无关三种^[8-9]。近年来, 探究不同生态系统服务之间的权衡与协同关系已成为生态系统服务研究的核心议题之一^[10-12]。从研究内容上看, 学者们主要集中在权衡与协同的理论基础、类型特征、驱动因素、动态变化及尺度效应^[13-14]；研究方法上则以统计学方法、空间分析方法、情景模拟方法及服务流动性分析方法为主^[14-17]；从研究尺度上来看, 集中在流域和城市群等较大空间尺度研究^[18-20]。针对生态脆弱的北方农牧交错带的研究主要集中在诸如土壤保持、碳固持及土壤风蚀等单一生态系统服务评估或气候、土地利用变化对不同生态系统服务的影响分析等领域^[20-23], 且研究区域多是在内蒙古赤峰市^[24]、鄂尔多斯市^[20]、呼包鄂榆城市群^[18]及整个北方农牧交错带^[25]等较大的空间尺度, 县域等较小空间尺度多种生态系统服务间复杂关系的研究较少。生态系统服务的尺度关联使得同一对服务在不同研究尺度上的权衡关系存在很大的差异^[13]。北方农牧交错带是从半干旱区向干旱区过渡的广阔区域, 区域内降水、温度和生物群落存在显著差异^[26], 大空间尺度上生态系统服务间的权衡/协同关系, 可能并不能代表小空间尺度上存在相同的关系。政府基于较大空间尺度上的生态系统服务间关系进行县域等较小行政单元的生态系统管理, 可能无法有效提高生态系统服务的供给能力^[27]。因此, 为了更好地理解生态系统服务间的复杂关系, 推进北方农牧交错带生态系统的管理实践工作, 亟需开展县域、乡镇等较小空间尺度下的多种生态系统服务间权衡/协同关系的对比分析研究。

地处干旱半干旱地区的北方农牧交错带是我国东部农区与西部牧区相连接的生态过渡带, 同时也是我国一个重要的生态脆弱区, 人类活动方式的多样性及特殊的地理位置导致该地区生态系统经常受到干扰和破坏^{[26, 28]}。内蒙古多伦县是北方农牧交错带的一个典型区域, 南部为上都河冲积区, 土壤肥沃, 是全县重要的粮食产区, 中部则是农牧混合地带, 北部以牧业为主。该县生态系统类型多样, 为人类提供了多种生态系统服务, 但随着气候变化和人类活动的增加, 多伦县生态系统服务供需矛盾日益突出, 产生了一系列的生态环境问题, 如草地退化、土地荒漠化及旱灾增加等^[29]。本文以多伦县为研究区, 基于2000年、2015年的遥感影像、气象、归一化植被指数(normalized difference vegetation index, NDVI)数据及社会统计资料等数据, 利用多种生物物理模型评估食物供给、碳固持、土壤保持、防风固沙及产水量5种关键生态系统服务的时空变化；分析多伦县5种服务热点区分布特征, 揭示该县不同区域的服务供给能力的强弱；此外采用相关分析法从县域、乡镇及不同地类3个尺度上分析生态系统服务权衡/协同关系。本研究着重于县域、乡镇等基本的行政管理单元, 研究结果旨在服务于管理实践, 为实现区域发展和生态保护提供科学支撑。

1 研究区域与数据 1.1 研究区概况

多伦县(41°46′—42°36′ N, 115°51′—116°54′ E)位于内蒙古自治区锡林郭勒盟东南部, 地处阴山北麓的东段, 浑善达克沙地的东南端, 辖3镇2乡, 全县总面积3.95×10³ km² (图 1)。该地区地势总体由西南向东北逐渐低缓, 呈盆地状, 地貌主要有低山、丘陵、河谷洼地、山前倾斜平原及高平台、堆积型沙丘等5种类型, 海拔范围是1149—1796 m；土壤类型主要是栗钙土、草甸土和风沙土, 成土母质多以风积、冲积沙为主；该区域地处森林、草原与沙地过渡地带, 主要分布有典型草原植被、草甸草原植被、沙地植被、沼泽植被等植被类型；气候属于中温带半干旱向半湿润过渡区的典型大陆性季风气候, 冬季严寒漫长, 无明显夏季, 风力作用强, 光照充足, 水热同季；年平均气温3.5℃, 年平均降水量400 mm左右, 属于锡林郭勒盟境内的丰水带。

1.2 数据来源

本研究使用到的基础数据包括土地利用/覆被数据、气象观测数据、DEM数据、归一化植被指数(NDVI)数据、土壤数据及社会经济统计数据6类：(1)2000年、2015年土地利用/覆被数据, 来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/Default.aspx), 空间分辨率为30 m；(2)多伦县地区及周边区域200 km范围内的23个气象站点的逐日降水量、平均气温、平均风速及日照时数等气象观测数据来源于中国气象数据网(http://data.cma.cn/), 本文利用1998—2002年、2013—2017年两个时期的气象数据平均值进行插值以克服单一年份气象极端值的影响；(3)DEM数据来源于地理空间数据云(http://www.gscloud.cn/), 空间分辨率为30 m；(4)2000年、2015年MODIS-NDVI数据来源于美国地质调查局(http://ladsweb.nascom.nasa.gov/data/), 空间分辨率为250 m；(5)土壤数据来源于寒区旱区科学数据中心(http://westdc.westgis.ac.cn/), 基于世界土壤数据库(HWSD)的中国土壤数据集(v1.1), 包含土壤类型、土壤质地、土壤有机碳含量和粒径比等数据, 空间分辨率为1 km；(6)社会经济统计数据来源于《内蒙古统计年鉴2001》和《内蒙古统计年鉴2016》。将所有的空间数据投影统一到UTM投影系下, 空间分辨率统一为100 m。

2 研究方法 2.1 生态系统服务的估算 2.1.1 食物供给

多伦县地位于农牧交错地带, 食物供给是农田生态系统和草地生态系统提供的一项重要的供给服务, 对该县人类生存和社会发展起到了至关重要的作用。本研究选用武文欢等^[20]提出的公式估算单位面积上的食物供给量, 其具体过程是将相应年份统计年鉴中谷物、油料和蔬菜等农作物产量及肉类、奶类产量按照正值的NDVI分别分配给耕地和草地栅格。计算公式如下：

(1)

式中, P_i为第i个栅格所分配的粮食产量或肉类与奶类产量(t hm^-2a^-1)；P_sum为农作物总产量或肉类总产量(t hm^-2a^-1)；NDVI数据为夏季(6—8月)MODIS-NDVI产品数据的最大值合成数据^[30], NDVI_i表示第i个栅格的NDVI；NDVI_sum为研究区耕地或草地的NDVI之和。

2.1.2 碳固持

碳固持服务指的是自然植被将大气中的CO₂捕获并封存的过程, 能有效减缓全球气候变暖趋势。碳固持的估算选取植物净初级生产力(net primary productivity, NPP)作为指标进行表征, NPP是地表碳循环的重要部分, 不仅反映了植被的生产能力, 也是判定区域生态系统碳汇和碳源的重要因子。本研究采用CASA模型^[31-32]估算研究区NPP, 公式如下：

(2)

(3)

(4)

(5)

(6)

式中, NPP(x, t)是像元x在t月份的净初级生产力(gC/m²)；APAR(x, t)是指像元x在t月份吸收的光合有效辐射(MJ/m²), ε(x, t)为像元x在t月的实际光能利用率(gC/MJ)；SOL(x, t)是像元x在t月份的太阳总辐射量(MJ/m²)；FPAR(x, t)是植被层对入射光合有效辐射的吸收比例(无单位)；T_ε1(x, t)和T_ε2(x, t)分别为低温和高温对植被光能利用率的影响系数；W(x, t)表示水分胁迫影响系数(无量纲)^[32]；SR(x, t)为像元x在t月份比值植被系数, SR_max与植被类型有关, 取值范围在4.14—6.17之间, SR_min取值为1.05^[33]；ε_max是植被在理想状态下的最大光能利用效率(gC/MJ), 取值参考朱文泉等^[33]提出的参数；最后得出年度植被净初级生产力。

2.1.3 产水量

生态系统水供给服务与人类生活息息相关, 在多个方面影响着人类福祉, 水资源是该区域农牧业发展的重要限制因素之一。本研究使用InVEST模型中Water Yield模块^[34]估算产水量, 计算公式如下：

(7)

(8)

(9)

(10)

(11)

式中, Y_xj为栅格单元x中土地利用类型j的年产水量(mm)；AET_xj为栅格单元x中土地利用类型j的实际蒸散发量(mm)；P_x为栅格单元x的年均降水量(mm)；PET_xj为栅格单元x中土地利用类型j的年潜在蒸散量(mm)；w_x为修正植被年可利用水量与降水量的比值(无量纲)；Z代表Zhang系数^[35], 用于表征降水的季节性分配特征；AWC_x为栅格单元x的可利用含水量；D_s为最大土壤深度(mm), 从第二次全国土壤调查数据中提取获得；D_r为根系深度(mm)^[36]；PAWC_x为栅格单元x的植被可利用水(mm), 根据土壤质地计算^[37]；sand、silt和clay分别是砂粒、粉粒和黏粒的含量(%), OM为有机质含量(%)。

2.1.4 土壤保持

人类活动干扰及全球气候变化等多种因素会使植被退化, 加剧水土流失, 过量的土壤侵蚀会导致土壤肥力下降, 农牧业生产力下降。本研究采用RUSLE方程估算多伦县的土壤保持量, 公式如下：

(12)

(13)

(14)

(15)

式中, SC表示年土壤保持能力(t hm^-2a^-1)；R表示降雨侵蚀因子(MJ mm hm^-2h^-1a^-1), 采用Wischmeier和Smith^[38]提出的月尺度公式计算, 其中X_i为第i个月的降水量(mm), X为年降水量(mm)；K为土壤可蚀性因子(t hMJ^-1mm^-1), 本文采用Williams和Arnold^[39]提出的利用土壤有机碳和颗粒组成数据进行计算的方法, sand、silt和clay分别是砂粒、粉粒和黏粒的含量(%), OM为有机质含量(%)；C表示植被覆盖与管理因子, 无量纲, 通过蔡崇法等^[40]提出的公式计算；fvc为植被覆盖度, 利用每年植被生长季(5—9月)的NDVI均值计算；P为水土保持措施因子, 无量纲, 由土地利用类型决定^[41]；LS表示坡长坡度因子, 无量纲, 本文采用van Remortel等^[42]编写的AML代码从DEM中提取LS因子。

2.1.5 防风固沙

位于浑善达克沙地东南端的多伦县, 土壤风蚀经常发生, 致使土壤粗化、土壤肥力和土地生产力下降, 对农牧业等生产活动产生严重影响, 并威胁着周围城市的生态环境。基于GIS技术和遥感数据, RWEQ模型已经应用于区域尺度土壤风蚀量的估算, 本文使用RWEQ模型^[43-44]估算防风固沙服务, 具体公式如下：

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

式中, SL_pv为防风固沙量(kg/m²)；SL_p为无植被覆盖的潜在风蚀量(kg/m²)；SL_v为有植被覆盖的实际风蚀量(kg/m²)；z为距离上风向的长度(m)^[45]；S_p为潜在关键地块长度(m)；S_v为实际关键地块长度(m)；Q_pmax为潜在风力的最大输沙能力(kg/m²)；Q_vmax为实际风力的最大输沙能力(kg/m²)；WF为气候因子(kg/m²)^[44]；EF和SCF分别表示土壤可蚀性因子(无量纲)和土壤结皮因子(无量纲)^[44]；COG为综合植被因子(无量纲)^[44]；基于野外实验, 采用滚轴式链条法测定锡林郭勒草原区地表粗糙度因子K′, 沙地为0.96, 草地为0.69^[44]；F表示防风固沙保有率^[46], 消除气候波动影响的防风固沙保有率可较好地反映植被对防风固沙的贡献作用。

2.2 生态系统服务热点区识别与权衡/协同分析

在生态系统服务评估的基础上, 识别多伦县的服务热点区, 有助于了解不同区域服务供给能力强弱情况。本文将研究区各个生态系统服务超过各自平均值的栅格定义为该类服务的热点区, 借助ArcGIS 10.5软件将这5种生态系统服务的热点区进行叠加分析。若某一栅格单元内5种服务的值均未超过各自所对应的平均值, 则该类栅格单元被定义为非热点区；如果只有1种服务的值超过平均值, 则该类栅格单元被定义为1类热点区；以此类推, 分别定义2类、3类、4类和5类热点区, 并在ArcGIS 10.5软件中得到多重生态系统服务热点区的空间分布格局^[20]。

本文基于SPSS 20.0软件, 以2000年、2015年5种生态系统服务的评估值为基础数据, 采用相关分析法对县域、乡镇以及地类尺度上的生态系统服务间权衡/协同关系进行量化。相关分析是生态系统服务间相互关系研究中常用的统计分析方法之一。当两两生态系统服务之间的相关系数为负, 且通过0.05水平的显著性检验, 则认为两种生态系统服务之间具有权衡关系；反之, 相关系数为正值且通过显著性检验, 则认为两种生态系统服务之间具有协同关系^{[8, 18, 30, 47]}。

3 结果 3.1 生态系统服务时空变化

2000年和2015年5种生态系统服务空间分布以及变化量分布如图 2所示。总体来看, 除防风固沙量外, 多伦县其余4种服务呈增加趋势。2000—2015年食物供给增加两倍, 年平均值由50.66 t/hm²提高到157.67 t/hm²；碳固持年平均值在两个年份分别为367.95 gC/m²、372.36 gC/m², 仅增长0.1%, 总体增幅并不明显；产水量的年均值在两个年份分别为110.25 mm、113.5 mm, 总体增量较小, 约为0.3%；土壤保持量年平均值由400.9 t/hm²提高到485.8 t/hm², 增幅为21%；潜在风蚀量由2000年的11.5 kg/m²减少到2015年的1.7 kg/m², 减幅为85%, 实际风蚀量由2000年的5.6 kg/m²减少到2015年的0.6 kg/m², 减幅为89%, 防风固沙量由5.9 kg/m²减少到1.1 kg/m², 减幅为81%, 两个年份的防风固沙保有率分别是51%、61.2%, 提高了10.2%。

食物供给空间格局在2000—2015年之间基本保持不变, 高值区主要分布在耕地, 低值区主要分布在草地；从变化量来看, 增量最高的主要分布在耕地。碳固持的空间格局变化明显：2000年呈现出“四周高、中部低”的特征, 2015年高值区主要分布在东部和北部, 西南部高值区范围有所减少, 中部高值区范围有所增加；两个年份高值区的主要土地利用类型是林地、耕地；从变化量来看, 增加区域主要分布在多伦县的中部和北部。2000年、2015年产水量的空间格局基本稳定, 大致呈现出：中部高、南部和北部低的分布态势；高值区的土地利用类型主要是沙地；从变化量来看, 自西北向东南方向增量的数值呈现出增加的趋势。两个年份的土壤保持空间格局基本稳定, 表现为从西南向东北递减的特征；2015年的土壤保持量相比于2000年, 高值区范围稍有增加, 增值区主要分布在西南部和东部地区, 西北部主要为土壤保持的减少区, 变化量的空间分布格局同碳固持基本相似。防风固沙量的空间格局基本保持不变, 高值区主要分布在土壤类型多为风沙土的北部, 低值区分布在南部栗钙土地区和中部草甸土地区；从变化量上看, 多伦县大部分为减少区域, 其中增值区分布在由未利用地转化为耕地和林地的区域。

3.2 热点区识别

多伦县多重生态系统服务热点区分布如图 3所示。2000年和2015年, 非热点区的面积百分比分别为11%、9%, 且主要的土地利用类型为水域和草地。1类热点区由26%增加到36%, 该类热点区的主要土地利用类型为草地；该县北部1类热点区主要以防风固沙为主, 中部以产水量为主, 南部以食物供给为主。2类热点区在两个年份的面积百分比分别为35%、40.7%, 面积增加5.7%, 主要的土地利用类型为草地；北部以碳固持和防风固沙服务为主, 中部以产水量与防风固沙为主, 南部多为食物供给与碳固持。3类热点区在两个年份的面积百分比分别为23.5%、13.2%, 面积减少10.3%, 以食物供给、碳固持以及防风固沙为主要的生态系统服务, 该类热点区的土地利用类型以草地为主。4类热点区的面积百分比由5.11%减少到0.7%, 主要分布在耕地, 以食物供给、碳固持、产水量与防风固沙服务为主。5类热点区的面积0.1%, 主要土地利用类型为林地。整体来看, 多伦县3类热点区面积明显减少, 1类热点区面积明显增加。

3.3 不同地类生态系统服务供给能力对比分析

2000年、2015年多伦县耕地和草地的面积之和均超过研究区总面积的80%, 所以本文重点关注耕地和草地。按照各服务的均值和总量进行统计, 分析地类提供生态系统服务的能力强弱, 结果如图 4所示。

在单位面积上, 食物供给能力耕地大于草地；碳固持供给能力草地大于耕地, 二者的单位面积总供给量约为研究区的42%, 其余土地利用类型中供给量较高的为林地；土壤保持供给能力耕地大于草地, 二者的单位面积总供给量约为研究区的38%, 其余土地利用类型中供给量较高的为林地；防风固沙量草地大于耕地, 二者单位面积总供给量达到研究区的70%以上, 其余土地利用类型中防风固沙量较高的为林地；产水量供给能力耕地大于草地, 但二者的总单位面积供给量只达到研究区的22%左右, 其余土地利用类型中供给量较高的为林地。从县域供给总量上看, 食物供给能力耕地大于草地；碳固持供给能力草地最强, 草地和耕地的供给总量约为研究区的85%；2000年耕地面积虽然小于草地面积, 但土壤保持能力大于草地, 2015年土壤保持能力草地大于耕地, 二者的总供给量在两个年份均达到研究区的80%以上；两个年份草地的防风固沙量均大于耕地, 二者总的防风固沙量均达到研究区的83%以上；产水量供给能力草地大于耕地, 2000年、2015年二者的供给量分别为研究区的58%、66%。

3.4 生态系统服务的权衡/协同关系 3.4.1 县域生态系统服务权衡/协同关系

在县域尺度上5种服务之间的相关系数均通过了0.01显著性水平的检验。相关分析表明(图 5), 2000年、2015年食物供给、碳固持、土壤保持两两服务间以及防风固沙与碳固持、产水量之间均为协同关系, 土壤保持与产水量、防风固沙之间以及碳固持与产水量之间表现为权衡关系；食物供给与产水量由权衡转变为协同关系, 食物供给与防风固沙由协同转变为权衡关系。目前来看, 在县域尺度上管理生态系统, 当同时提升食物供给、碳固持及土壤保持3种服务供给能力时, 可能会导致产水量以及防风固沙能力的下降。

3.4.2 乡镇与耕地、草地中生态系统服务权衡/协同关系

在乡镇尺度上, 除多伦诺尔镇外, 4个乡镇中各服务间相关系数均通过0.01显著性水平检验, 多伦诺尔镇中2000年食物供给与防风固沙之间未通过显著性检验, 碳固持与防风固沙之间通过0.05显著性水平检验。相关分析表明, 食物供给、土壤保持、碳固持三者在各乡镇中均为协同关系, 并与县域尺度相同；而产水量与碳固持之间以权衡关系为主, 与县域尺度基本一致；产水量与防风固沙、土壤保持、食物供给间的关系以及防风固沙与食物供给、碳固持、土壤保持之间的关系随着乡镇的不同而发生变化；蔡木山乡与多伦诺尔镇中各服务间的权衡/协同关系与县域的相似度更高；在各乡镇中提高碳固持能力的同时, 会降低产水量的供给能力, 可能会提高土壤保持、防风固沙及食物供给能力。

耕地和草地中各服务间的关系均通过0.01显著性水平检验, 相关分析表明, 两种地类中5种服务间的权衡/协同关系差异明显：耕地中各服务之间的关系以权衡为主, 草地中以协同为主；草地中各服务间的关系比较稳定, 而耕地中有4对服务的权衡/协同关系发生变化；耕地和草地中存在明显差异的服务为食物供给与土壤保持、碳固持与土壤保持以及防风固沙与产水量, 这3对服务在耕地中为权衡关系, 在草地中为协同关系；与耕地相比, 草地中各服务之间的权衡/协同关系与县域的相似度更高；提高耕地食物供给、碳固持服务的同时, 可能会降低土壤保持、产水量的供给能力, 同时提高防风固沙能力；提高草地食物供给、碳固持服务供给能力的同时, 可能会同时提高防风固沙、土壤保持和产水量供给能力。

4 讨论 4.1 多伦县生态系统服务的时空变化分析

2015年食物供给相比于2000年增加两倍, 增幅非常明显, 主要原因是2000年多伦县开始实施退耕还林/草政策, 大面积的退耕使得粮食产量急剧减少；2003—2015年多伦县逐步利用废弃耕地, 保护和利用后备土地资源, 粮食产量逐渐增加^[48]。碳固持、土壤保持与产水量3种服务的总量分别增加1%、21%和2.5%, 这种变化在一定程度上受到降水量的影响：碳固持受到降水量的制约, 体现在CASA模型中的水分胁迫系数；降水量通过RUSLE方程中的降雨侵蚀R因子影响土壤保持功能；产水量为降水量与实际蒸散发量之差^[49]；此外碳固持与土壤保持还受到植被覆盖度的影响, 2000—2015年研究区的降水量增加5%, 植被覆盖度增加17%, 相应的3种服务也有所增加。2000—2015年研究区平均风速降低、降水量及植被覆盖度增加, 使该地区的潜在风蚀量和实际风蚀量都大幅下降；虽然防风固沙量减少, 但防风固沙保有率提高, 风蚀状况已有所好转, 这一结果与巩国丽等^[46]评估内蒙古典型草原区防风固沙能力的结果相似。

土壤保持供给量的高值区集中在研究区的南部, 与坡度值高的地区空间分布基本吻合, 表明该地区土壤保持也与坡度密切相关, 因此, 研究区应结合“京津风沙源治理”等生态恢复工程, 加强对坡地植被的保护。产水量除受降水影响外, 还与蒸散及植被覆盖度有关, 高值区多分布在中部的林地, 这可能是因为, 林地虽覆盖度较高, 蒸发较快, 但是林地雨水截留量大, 径流产生少, 所以产水量也相对较高^[50]；而北部和南部地区有大面积的草地和耕地, 植被覆盖度相对较高, 蒸散量较大, 易产生径流, 导致产水量较小, 这与张雪峰等^[50]在锡林郭勒草原区的研究结果相似。多伦县防风固沙量的分布趋势与风蚀量的分布趋势基本一致, 风蚀量大的区域主要分布在研究区北部, 靠近浑善达克沙地南端, 这部分区域土壤多为风沙土、易蚀性大, 潜在风蚀量很大^[46]。2000年多伦县的西南部耕地地区5种服务的部分热点具有空间异质性并表现为多重服务热点区, 而2015年南部耕地的食物供给量大幅度增加, 多重服务热点区明显减少, 这和Li和Wang^[12]以位于北方农牧交错带的银川为研究区的研究结果相符, 该研究认为农田不利于提高农牧交错区的防风固沙、土壤保持及产水量服务的供给能力。另外, 本研究发现研究区内林地的产水量、碳固持及土壤保持服务单位面积供给能力均强于草地与耕地, 适当增加林地的面积, 有助于提高多伦县的生态系统服务供给能力。

4.2 多伦县生态系统服务权衡/协同关系

多伦县碳固持与土壤保持之间为协同关系, 这是由于植被覆盖度的增加提高了碳固持能力, 而且可以减少雨水对土壤的侵蚀, 从而加强了土壤保持能力^[20]。食物供给与碳固持、土壤保持之间为协同关系, 这与之前研究^[51]发现的食物供给服务与调节服务之间为权衡关系有所不同, 原因可能是食物供给的评估方法不同, 因此, 由于生态系统服务的评估方法不同, 使得同一对服务之间的关系无法进行对比参照^[52]。本研究表明研究区的防风固沙与碳固持、产水量之间为协同关系：Du等^[53]在黄河流域的研究发现防风固沙量与植被覆盖度呈正相关关系, 多伦县植被覆盖度大的地区具有较高的碳固持功能, 同时良好的植被覆盖提高了防风固沙量, 所以防风固沙与碳固持之间为协同关系；此外, Du等^[53]认为土壤湿度大的区域防风固沙量也大, 本研究区中土壤湿度大的地区产水量的供给能力较强, 减少了实际风蚀量, 同时增加了防风固沙量, 所以防风固沙与产水量之间为协同关系。本研究中土壤保持与防风固沙服务之间为权衡关系, 主要是因为多伦县土壤保持高值区分布在地形起伏度大的南部, 而风蚀潜力较大, 即防风固沙量大的区域分布在地形平缓的北部, Jiang等^[54]以黄土高原为研究区的结果也表明土壤保持量较高的区域对应的地形起伏度较大, 而风蚀潜力较大的区域地形起伏度较小, 二者具有明显的空间异质性。产水量与土壤保持、碳固持服务之间主要为权衡关系, 这与Li等^[22]以内蒙古中部农牧交错带为研究区对这3种服务权衡/协同关系的研究结果相同, 增加植被覆盖度在提高碳固持及土壤保持能力的同时, 会增强地表蒸散量, 从而降低产水量^[50]。目前来看, 多伦县最重要的是产水量与土壤保持、碳固持以及食物供给之间的权衡关系, 在草地中产水量仅与土壤保持服务具有权衡关系, 因此, 该县可以通过适当增加草地面积来调节产水量与其他服务间的关系。

4.3 生态系统服务权衡/协同关系的尺度效应

通过对比2015年乡镇与县域尺度上10对服务间的权衡/协同关系发现, 乡镇中仅有4对服务(食物供给与碳固持、食物供给与土壤保持、食物供给与防风固沙、碳固持与产水量)之间的关系与县域相同, 其余6对服务间的关系会随着乡镇的不同而变化；与此同时, 在耕地、草地两种土地利用类型中仅有3对服务(食物供给与碳固持、碳固持与防风固沙、土壤保持与防风固沙)之间的关系与县域保持一致, 其余7对服务间的关系在两种地类中有明显的差异, 这说明县域尺度上的权衡/协同关系无法代表乡镇、地类尺度存在同样的关系；其他学者在中国北方呼包鄂榆地区^[18]和宁夏回族自治区^[21]的研究也发现, 服务间关系具有明显的空间尺度依赖特征。土地利用类型是影响生态系统服务变化的关键因素, 进而改变生态系统服务间的权衡/协同关系^[12]；本研究发现多伦县气候因素(气温、降水、太阳辐射等)空间差异较小, 而土地利用的空间异质性较大, 且研究区耕地中服务间的关系以权衡为主, 草地中服务间的关系以协同为主；例如耕地和草地面积比例为1.2的大北沟镇, 各服务间的权衡/协同关系与耕地保持一致, 耕地和草地面积比例为0.3的滦源镇, 服务间的关系与草地保持一致；因此, 土地利用类型的差异可能是乡镇之间、乡镇与县域之间服务的权衡/协同关系明显不同的原因之一。此外, 各乡镇的土壤类型、海拔高度也有较大差异：本研究区中栗钙土上食物供给、碳固持和土壤保持服务的能力较强, 草甸土上碳固持和产水量的能力较强, 风沙土上防风固沙的能力较强；海拔方面, Liu等^[55]发现服务间的关系会随海拔的不同而发生变化, 在本研究中仅发现防风固沙和产水量服务间的权衡关系会随海拔升高而加强, 未发现其余服务间的关系随海拔变化而发生明显变化的现象。土地利用类型、海拔和土壤类型的组合可能是权衡/协同关系具有空间差异性的原因之一, 这将是下一步研究中的关键问题。

北方农牧交错带位于东部季风区向西北干旱区的过渡带, 气温、降水等气候条件的空间异质性明显, 同时人类活动方式的不稳定性导致土地利用类型多样且易变, 选择适当的空间尺度对该区域进行生态系统管理, 对我国北方生态安全及促进经济、社会和生态的协调发展具有较强的现实意义。县域作为我国基本的管理单元, 本研究以多伦县为例, 从县域、乡镇以及不同土地利用类型的角度出发, 得出该县5种服务间的权衡/协同关系具有空间依赖性, 研究结果可为政府针对农/牧区制定生态保护政策提供依据。多伦县在实施退耕还林还草政策后, 林地面积不断增加, 逐渐成为该区域重要的土地利用类型, 本文仅考虑了耕地和草地上服务间的权衡/协同关系, 今后了解林地中服务间的关系也较为重要。此外, 在服务间关系的研究方法方面, Hao等^[25]在分析生态系统服务间相互关系时指出, 生态系统服务二维散点图中的点可能并不是围绕某一条线分布, 而是趋向于分布在某一条线(约束线)下, 约束线代表了一种服务在另外一种服务作用下的分布范围或者能够达到的潜在最大值；当约束线出现拐点(阈值点)时, 在点的两侧两种服务间的相互关系或相互关系的影响因素将会发生变化；相关性分析仅表征生态系统服务权衡/协同的数值关系, 无法揭示生态系统服务间相互作用的阈值特征。因此, 结合生态系统服务机理深入研究权衡与协同的阈值效应, 这可能是下一步工作的重要方面。

5 结论

通过对多伦县生态系统服务时空变化及权衡协同关系的研究, 本文主要得出以下结论：(1)2000—2015年多伦县防风固沙量减少81%, 其余4种生态系统服务呈现出增加趋势, 其中食物供给增加两倍, 土壤保持增加21%, 碳固持和产水量分别增加0.1%、0.3%；(2)该县1类以及2类热点区面积增加, 3类、4类以及非热点区面积减少；面积比例超过50%的草地在县域提供5种服务的能力强于耕地；(3)2000年、2015年多伦县主要的权衡关系存在于产水量与土壤保持、产水量与碳固持以及产水量与食物供给之间；草地中产水量与碳固持、食物供给之间主要为协同关系, 可以通过适当增加草地面积来调节产水量与其他服务间的关系；(4)由于土地利用类型的不同, 县域、乡镇、草地以及耕地上5种服务的权衡/协同关系存在明显差异：乡镇尺度上仅有4对服务间关系与县域保持一致；地类尺度上仅有3对服务间关系与县域保持一致。