生态系统服务是指人类从生态系统中所获取到的所有惠益^[1-2], 同时具有多重特性, 如服务类型的多样性、服务关系的不均衡性、空间分布的差异性等, 但是在外力的作用下, 生态系统服务之间的关系发生转变, 权衡或者协同关系逐渐显现^[3]。所以, 明晰生态系统服务之间的关系和作用已经成为区域发展必不可少的步骤, 也是制定区域绿色经济发展和生态环境等一系列可持续发展政策的前提条件。

随着人们对生态环境重视程度的加深, 促使越来越多的学者对生态系统服务进行探究, 并且取得了丰厚的学术成果。在探究生态系统服务的方法上普遍采用空间制图^[4-6]、情景分析^[7-9]、玫瑰图^[10]、模型模拟^[11]等来研究生态系统服务的权衡协同关系的空间特征、尺度效应以及影响机制等。而在探究生态系统服务关系时发现, 供给和调节服务是相互抑制的, 处于权衡的关系^[12-13], 例如植被生长状况良好时, 能提供的调节服务水平高, 但是粮食生产的能力会降低, 导致供给服务的水平下降^[14-16]。由于生态系统服务关系存在区域差异性, 不同服务之间在不同的研究区域所表现的权衡和协同关系不同^[17-18]。Qiao等^[19]通过不同的空间尺度对太湖流域上生态系统服务间的关系进行探究, 研究发现全域尺度和县域尺度下生态系统服务间的关系均会有所改变; Turner等^[4]从不同的空间尺度对丹麦11种生态系统服务进行研究, 发现植被固碳、旅游休憩等文化服务在50km的空间范围内聚集分布, 而畜牧量、粮食产出和淡水产出等供给服务在150km的空间范围内聚集分布。上述的研究大多集中在行政区划尺度和整体尺度上探究生态系统服务, 同时对于权衡与协同关系分析多基于统计关系的数量分析, 来反映区域的整体差异性, 缺少区域内部时空差异的空间表达和系统的了解生态系统服务关系形成的内在机制以及对自然生态系统内部异质性的研究。

秦巴山区是中国地理格局的重要枢纽、大尺度东西向的生态廊道, 以复杂、多样、独特的自然环境孕育了丰富的天然动植物资源, 对其生态环境变化的研究具有重要意义。根据秦巴山区当地的区域特点和特殊的地理位置选取了五种生态系统服务(粮食生产、NPP、土壤保持、水资源供给、生境质量)探究彼此之间的关系, 从不同尺度分析生态系统服务的异质性和定量评估生态系统服务权衡关系, 从宏观和微观两方面表现出秦巴山区多维地带性对生态系统服务关系的影响, 运用分解和综合的视角凸显生态系统服务功能间的相互变化。该研究可以协助决策者在适当的空间尺度制定相应的管理策略, 同时为区域间的山水林田湖草生态保护提供基础性参考。

1 研究区概况及数据来源 1.1 尺度划定

全域尺度：秦巴山区是我国南北过渡地带, 特殊的地理位置使得该区域在地理上有明显的地域分异规律, 作为南北生态屏障和安全区具有很突出的研究价值; 综合分区尺度：综合自然地理分区主要结合自然地理环境各组成成分(气候、地貌、地形、土壤等)的相互关系, 进行自然区划研究, 探讨自然区划之间的差异性。本文根据秦巴山区的气候、地貌、土壤等差异性, 将区域划分为14个子区域, 研究生态系统服务关系在子区域的异同; 典型样区尺度：秦岭-大巴山是南北过渡带, 其气候条件、地貌特征在南北坡存在差异性, 因此选取了秦岭和大巴山的南北坡做为典型样区; 样点间隔尺度：本研究通过ARCGIS软件建立渔网, 按照一定的间距设置随机点, 分别设定500m、1km、3km、5km、10km、20km、30km、40km、55km、85km、100km、200km共12种空间尺度, 探讨每种间隔尺度下生态系统服务的权衡协同关系, 4种尺度划分如图 1所示。

1.2 数据来源及处理

本研究主要利用的数据包括遥感数据、高程数据、土地利用数据、气象数据、统计数据等多源数据集。DEM数据主要来源于地理空间数据云; 土地利用数据来源于中国科学院资源环境科学数据中心; 2005、2010、2015年陕西省、河南省、湖北省、甘肃省、四川省和重庆市的粮食产量、种植面积等统计数据, 来源于各省及市域的统计年鉴; 遥感数据主要是2005年、2010年、2015年250m空间分辨率MOD13Q1数据; 气象数据源为中国气象科学数据共享服务网。

2 研究方法 2.1 植被净初级生产力估算(NPP)

CASA模型最早是在1993年由Potter等^[20]提出的, 主要是从遥感影像中提取的NDVI数据, 并通过光能利用率的原理来估算全球生态系统的净初级生产力(NPP)。计算表达式如下：

(1)

式中, APAR(x, t)代表的是t月中像元x吸收的光合有效辐射(MJ/m²); ε(x, t)代表的是t月中像元x的实际光能利用率(gC/MJ); NPP(x, t)为像元x在t月份的净初级生产力(gC/m²)。

2.2 土壤保持估算

Wischmeier和Smith^[21]提出通用土壤流失方程(Universal Soil Loss Equation, USLE), 后经改进形成修正的通用土壤流失方程^[22] (Revised Universal SoilLoss Equation, RUSLE), 其计算方式如下：

(2)

式中, A为年均土壤侵蚀模数(t hm^-2 a^-1); R为降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1), 采用Wischmeier和Mannering ^[23]提出的公式估算所得; LS为坡长、坡度因子, 无量纲, 坡长因子采用Wischmeier和Mannering ^[23]提出的计算公式, 坡度因子采用Mc Cool等^[24]提出的计算公式; K为土壤可蚀性因子(t hm² h MJ hm^-2 mm^-1); C为植被覆盖与管理因子, 无量纲, 参照了蔡崇法等^[25]的研究成果; P为水土保持因子, 无量纲。

2.3 水资源供给估算

本文采用In VEST模型中的产水模块计算秦巴山区的产水量, 基本原理如下：

(3)

式中, Y(x)为栅格单元x的年产水量(mm); AET(x)为栅格单元x的年实际蒸发量(mm); P(x)为栅格单元x的年降水量(mm)。

(4)

(5)

(6)

式中, PET(x)为栅格单元x的潜在蒸散量(mm); ET₀(x)表示栅格单元x的参考作物蒸散(mm); K_c(x)为作物蒸散发系数; AWC(x)为土壤有效含水量(mm); W(x)为经验参数; Z为Zhang系数^[26]。

2.4 生境质量估算

InVEST软件的生境质量模块测算生境质量指数评价生境质量, 计算结果范围为：0—1, 值越大代表栖息地质量越高, 生境质量质量越高。计算公式为：

(7)

式中, Q_xj是土地利用类型j中栅格x的生境质量; D_xj是土地利用类型j栅格x的生境胁迫水平; k是缩放参数(常数)。

(8)

栅格y中胁迫因子r(r_y)对栅格中生境的胁迫作用为i_rxy, 结合相关文献^[27-29], 选取了研究区域的耕地、园地、公路、铁路、城镇村及工矿用地和其他用地作为其生境质量的威胁因子, 对其赋值, 如表 1所示。

(9)

式中, d_xy为栅格x与y之间的线性距离; d_rmax为胁迫因子最大有效威胁距离; W_r为胁迫因子的权重; β_x为栅格x的可达性水平; S_jr为土地利用类型j对胁迫因子的敏感性。

2.5 粮食产量估算

粮食生产与农田的面积呈现正向关系, 但是不同的种植物对粮食产量具有不同的影响。根据相关文献^[30], 统计分析建模是实现粮食产量空间化的一种常用的方法, 所以本文选用旱地和水田的面积作为自变量, 县级尺度的粮食总产数据作为因变量, 建立多元线性回归方程来实现粮食产量的空间化。

3 结果与分析 3.1 不同尺度下生态系统服务关系分析 3.1.1 全域尺度

2005—2015年生态系统服务间的相关系数进行计算, 并对计算结果进行t检验, 根据t分位数表将结果划分为6个等级^[31], 协同^{* *}(r＞0, 0.01＜P＜0.05)、协同^*( r＞0, 0.05＜P＜0.1)、协同( r＞0, 0.1＜P)、权衡(r＜0, 0.1＜P)、权衡^*(r＜0, 0.05＜P＜0.1)、权衡^{* *}(r＜0, 0.01＜P＜0.05), 分析不同生态系统服务间权衡关系的空间格局和相关关系(图 2)。图 2的饼状图填充大小代表相关系数值, 颜色深浅代表值的大小, 值越大颜色越深, 蓝色代表正相关, 粉色代表负相关。

水资源供给与生境质量以权衡关系为主, 南阳盆地、汉中盆地、大巴山北部的权衡关系较突出; 水资源供给与土壤保持在部分区域(60.88%)呈现协同关系, 但权衡关系在汉江上游尤为显著, 主要是人口密集、植被对水分的消耗导致水资源供给量的减少, 土壤保持量在增加, 导致在汉江上游两者呈现显著权衡关系; 水资源供给与NPP以协同关系为主, 面积为56.16%, NPP有很强的区域依赖性, 植被覆盖度的增加会提升NPP指数; 水资源供给与粮食产量以协同关系为主, 汉中盆地以南呈现权衡的关系; NPP与生境质量是协同关系, 面积为59.69%, 显著协同的区域分布在伏牛山区以及秦岭山脉西侧, 大巴山的东南以及东昆仑山区两者关系呈现显著权衡的关系; NPP与土壤保持为协同关系, 显著协同比例高达36.36%, 主要集中在昆仑山、岷山以及大雪山一侧; 粮食产量和NPP、水土保持、生境质量之间均是权衡关系, 面积分别为61.27%、57.62%、68.69%, 因为粮食种植与地形有关, 主要分布在平原和盆地, 同时区域的NPP、水土保持和生境质量值都较低; 土壤保持和生境质量以协同关系为主, 面积为55.25%, 在豫西汉中盆地和秦岭山脉西侧呈现显著权衡关系。

3.1.2 综合分区尺度

根据划分的综合分区, 得到子区域两两生态系统服务的相关关系(图 3)。水资源供给和NPP整体上为协同的关系, 但在子区域4、6、7、9、12、13、14中呈现权衡关系, 原因是子区域内植被覆盖率提升的同时增加了地表植被的蒸腾作用, 黄淮海平原和南阳盆地人口密集, 生活和工业用水增多, 促使子区域水资源供给量降低; 水资源供给和粮食生产整体上是协同关系, 而在子区域1、5、6、12呈现权衡关系, 原因是坡度小的地区适合粮食生产, 地貌类型以高山为主的区域不适合粮食生产, 同时气候类型属于亚热带气候, 水资源供给量丰富, 促使子区域两者关系是权衡, 同样坡度较小, 土壤侵蚀量小, 所以在子区域1、2、5、7、12、13的水资源供给和土壤保持呈现权衡关系; 水资源供给与生境质量总体以权衡关系为主, 仅有区域6为协同的关系, 相关系数0.25, 主要因为秦岭山脉北坡生境质量丰富, 而区域内水资源供给量较为丰富; NPP与生境质量在区域5呈现权衡关系, 主要因为海拔高, 高山气候类型, 促使生境质量单一, 而原始植被状况良好, 破坏程度低; NPP与粮食生产整体上呈权衡关系, 而在子区域2、3、4、5为协同关系, 主要原因是气候和地形综合影响, 使得区域粮食生产能力和植被覆盖度低; 粮食产量和土壤保持在子区域1、2、3、4、5呈现协同关系, 主要原因是坡度在10—15°之间, 实际和潜在侵蚀量较小, 促使土壤保持能力弱, 而气候条件也会影响粮食生产; 生境质量和土壤保持以协同关系为主, 但在区域9却呈现权衡关系。

3.1.3 典型样区尺度

本文整理了2005、2010、2015年秦岭和大巴山的南北坡的五种生态系统服务值(图 4), 研究出生态系统服务关系在特定年份发生转变。水资源供给与NPP在2005年大巴山的北坡、2010年秦岭南北坡呈现权衡关系; 水资源供给和粮食产量在2015年四个典型区表现为权衡关系; 水资源供给和土壤保持在2010年秦岭北坡表现为权衡关系; 水资源供给与生境质量在2015年大巴山南坡表现为协同关系; NPP与生境质量在2015年大巴山南坡表现为权衡关系; 粮食生产和土壤保持在2005、2010、2015年大巴山南坡均表现为协同关系, 2010年秦岭北坡、大巴山北坡、2015年大巴山北坡均表现为协同关系; 土壤保持和生境质量在2005、2010、2015年大巴山的南北坡均表现为权衡关系等, 通过研究典型区服务间关系得出生态系统服务在随着时间的推移而发生了改变。

3.1.4 样点间隔尺度

根据计算结果, 将样点尺度下生态系统服务关系大致分为4类(图 5)：①生态系统服务以协同关系为主, 并且协同效应随着尺度的增大而增强。比如：NPP与生境质量、水资源供给, 水资源供给与粮食生产、土壤保持; ②生态系统服务以协同关系为主, 并且协同效应随着尺度的增大而减小。比如：NPP与土壤保持; ③生态系统服务以权衡关系为主, 并且权衡效应随着尺度的增大而增加。比如：NPP与粮食生产、粮食生产与生境质量、水资源供给与生境质量; ④生态系统服务之间的关系随尺度的增加由权衡变为协同和由协同变为权衡。比如：粮食生产与土壤保持以及生境质量与土壤保持。

3.2 4种研究尺度下生态系统服务关系对比分析

四种尺度下部分区域成对生态系统服务之间的相关系数(图 6)。在子区域4和大巴山南坡的NPP与粮食生产、粮食生产与土壤保持相关系数分别为0.056、0.032、0.061、0.025, 并且都通过了P < 0.05显著性检验, 在全域尺度和样点间隔尺度下两者关系却为权衡关系; 子区域7、12的水资源供给与土壤保持相关系数分别为-0.051、-0.07, 并且都通过了P < 0.01显著性检验, 而全域尺度、典型样区尺度、样点间隔尺度下两对关系却为协同关系; 水资源供给与粮食生产在秦岭、大巴山南坡、子区域12的相关系数分别为-0.042、-0.093、-0.087, 并且都通过了P < 0.01显著性检验, 而全域尺度和样点间隔尺度下两对关系却为协同关系。

4种尺度下完整的生态系统服务权衡与协同关系(表 2), 在全域尺度下探究生态系统服务之间的权衡与协同关系, 能比较宏观的揭示服务之间的关系。随着尺度的不断缩小, 在中尺度范围内(综合分区尺度)发现生态系统服务之间的差异性, 由于尺度之间基础条件的差异性, 造成在中尺度下服务之间关系的转变。在典型区之间生态系统服务关系也有所区别, 由于南北坡的气候、土壤、地貌等差异, 造成南北坡生态系统服务关系的差异性, 并在特定的年份服务之间的关系发生转变。在小尺度下研究生态系统服务之间的权衡协同关系也会随着尺度的差异而逐渐发生变化, 服务之间的关系随着尺度的增大而逐渐呈现不同方式的变化。

4 讨论与讨论 4.1 讨论 4.1.1 尺度差异性的驱动因素

秦巴山区是我国南北过渡地带, 过渡性对东西方向高大山脉的阻隔和分异作用产生影响, 坡向和垂直分异使得南北方向的土壤和植被具有多重分异的结果。东西向延伸的秦巴山脉, 不仅对南北方向的地理分异规律产生影响, 而且还具有更明显的阴阳坡和侧翼的分异^[31], 造成地形、地貌、气候等多方面的差异性, 会对生态系统服务关系产生影响, 促使生态系统服务关系的差异性。秦巴山区土壤层从南向北依次是半淋溶土、初育土、淋溶土、人为土、铁铝土等过渡形式; 气候类型从南温带、北亚热带、中亚热带过渡; 地貌类型也从秦岭大起伏高中山、豫西汉中中山盆地、南阳盆地、大巴山大起伏高中山过渡。地形、土壤、气候的差异性主要是受到纬度地带分异、经度分异、垂直地带性、坡向分异等自然规律的影响, 促使了秦巴山区高度复杂性、多样性和异质性。同时在区域内存在非地带性作用, 例如山体效应^[31-32], 也会增加秦巴山区环境的复杂性和异质性。

秦巴山区具有多维地带性结构, 高程、坡度因素会对整个生态系统服务造成不同的影响, 如图 7。大面积的耕地都处在DEM<2000m和坡度>25°, 大量的的梯田和坡耕地存在, 适合粮食的种植, 但是区域内人口密度较大, 在一定程度上会破坏植被的生长环境, 促使生境质量和植被的碳储量下降。森林和草地覆盖度高的区域大多集中在DEM>4000m和坡度在15—25°高海拔山区, 草地和森林的抗土壤侵蚀力比较强, 大面积的阔叶林分布, 人工林面积小, 故有较大的生态效益, 所以NPP、土壤保持、生境质量等功能较强。草地的分布也会增加地表覆盖, 更进一步的增强了地表抗侵蚀能力, 同时草地具有防风、固沙、保土、调节气候、净化空气、涵养水源等生态功能, 促使草地覆被的区域生态效益也非常突出。

近年来, 退耕还林还草、生态效益补偿、生境质量主体功能区的建设等国家政策的提出, 促使秦巴山区生态结构发生转变。相较于2005年, 2015年秦巴山区的耕地面积减少了10.04%, 林地面积增加了5.03%, 草地增加了5.24%, 城市用地增加了25.54%, 水域增加了4.13%, 未利用地减少了9.54%。土地利用形式的转变在一定程度上影响景观格局, 从而对生态系统生态系统服务产生影响。经研究发现仅有农作物生产服务与人类活动之间的关系是正相关, 凸显了人类活动对生态系统服务的影响程度^[33-34]。随着经济的发展, 大量的非建设用地转换为建设用地, 一方面导致了土壤保持能力的下降, 植被生长状态被破坏, 降低了NPP值, 另一方面也会降低地表的蒸散作用, 导致水资源供给量的增加, 从而促使一部分地区水资源供给和上述几种生态系统服务呈现权衡的关系, 这些结果和已有的研究成果相类似^[35-36]。粮食生产与土壤保持、NPP、生境质量之间是权衡关系。当为了提升粮食的产量, 会对其他的土地类型进行开垦和修复, 虽然提升了粮食产量, 但相对应的林地、草地的面积减少, 造成生境质量和NPP服务大幅度减少, 引发一系列的环境问题, 所以在增加粮食产量的同时也要考虑生态问题, 促使生态系统良性发展。

4.1.2 尺度方法的不统一性

本文为了凸显4种尺度下生态系统服务间权衡与协同关系的差异性, 采用了3种不同的方法来探究不同尺度下生态系统服务关系。在不考虑区域内部结构差异性的情况下, 采用相关分析法探究秦巴山区生态系统服务间关系的整体趋势, 凸显整体性; 基于R语言探究子区域(中、小尺度)生态系统服务的关联性; 随机样点法定量描述样点间隔尺度生态系统服务之间的关系; 样点尺度的选择具有随机性, 根据所设的空间间隔来进行相应的计算, 可以对比不同间隔尺度下成对服务关系, 来体现出小范围间区域特征的差异性。但是文中没有考虑在相同方法下生态系统服务关系究竟是如何变化的, 用相同的方法来衡量更具有说服力, 所以怎么样使用同一种方法在不同的时空尺度下都能清晰明了的展现出生态系统服务之间的关系转换是今后努力的方向。

4.2 结论

本文利用2005、2010和2015年数据, 探究秦巴山区大、中、小尺度下NPP、土壤保持、生境质量、水资源供给、粮食生产5种生态系统服务之间的权衡协同关系, 研究发现生态系统服务彼此关系在不同的研究尺度下存在差异性。结果表明：

(1) 全域尺度、综合分区尺度、典型样区尺度下生态系统服务之间的关系相同, 但在综合分区尺度下存在特定分区中生态系统生态系统服务之间的转变, 与全域间的生态系统服务关系出现差异; 典型样区尺度下与全域尺度在年际间发生变化, 在特定的年份服务之间的关系相互转变。

(2) 样点间隔尺度下协同与权衡关系在不同的间隔尺度中存在差异性, 协同、权衡效应随着尺度的增加而增强; 协同效应随着尺度的增加而减少; 权衡与协同关系随着尺度的增大相互发生转变。