生态系统为人类提供支持、调节、供给、文化等多种服务。近百年, 剧烈的土地利用与土地覆盖变化改变了陆地生态系统的格局、过程、功能, 导致生态系统服务及其相互关系的极大变化。生态系统服务以多种复杂方式相互作用, 此消彼长的权衡、相互增益的协同是生态系统服务之间的典型关系^[1-4]。权衡通常出现在供给与调节服务之间^[5]或两种调节服务之间^[6-7]。探讨生态系统服务的时空动态及其相互关系, 分析其类型特征、响应速率、时空格局、区域差异和驱动机制等^[8-9], 是可持续性研究领域的紧迫问题, 可为土地利用决策和区域生态系统管理提供理论依据^[10-11], 通过增强协同效应、削弱权衡效应^[8]实现生态系统服务利益最大化具有重要意义。然而, 目前关于生态系统服务相互作用、相互联系的机理研究及尺度效应尚不明确, 仍不足以支持当前的土地利用决策, 生态系统服务权衡与协同关系的时空动态变化特征及驱动机制有待深入研究^{[4, 9, 12]}。

生态修复工程对生态系统服务相互关系的影响是近年来的研究焦点。覆盖范围最广、影响人数最多的退耕还林还草工程, 通过提供现金和粮食补贴在耕地上造林、种草, 以增加植被覆盖、控制水土流失、改善农民生计^[13-14]。正如预期, 工程在增加植被覆盖度与碳蓄积^[15-17]、控制土壤侵蚀^[18-19]、增加季节产水量^[11]、减少河流年输沙量^[20]、促进非农就业^[21-22]等方面发挥了重要作用。然而, 局部研究也得出工程对其他生态系统服务产生不利影响, 比如导致耕地面积减少、威胁粮食安全^[23-25]、加剧社会不平等^[26], 加剧黄土高原缺水并影响水沙关系^[27], 对西北土壤肥力和植被覆盖不利^[28], 导致四川的生物多样性损失和天然林损毁^[29]。因此, 退耕还林还草工程虽然提升了部分关键生态系统功能, 却造成了生产和生态用地之间的潜在冲突, 改变了生态系统服务之间的相互作用^[30], 降低了生态系统服务协同效应, 使其对生态系统服务及其相互关系的影响存在争议^[11]。

退耕还林还草增加了调节服务、减少了供给服务, 可能加剧粮食生产与生态保护之间的权衡^[23-25]。工程规划了退耕还林还草目标, 但是, 尚不清楚哪种方式能更好地增加生态系统服务的提供、同时最大限度地减少生态系统服务的权衡^[11]。迄今为止, 针对退耕还林还草影响生态系统服务之间关系的研究集中在黄土高原、云贵高原等区域, 然而较少关注其时空尺度变化^{[1, 31]}。已有研究表明, 退耕还林导致黄土高原植被覆盖与产流、土壤含水量的权衡^[23], 土壤含水量与林下植物多样性之间的权衡最高, 且沿降水梯度、随林龄增加而变化^[32]。碳固定、土壤保持和水文调节服务之间的关系不仅是权衡或协同, 还有约束效应, 其中降水供给是关键制约因素^[31]。最大限度地减少生态系统碳与水服务之间不必要的权衡, 被认为是促进黄土高原可持续发展的必要条件^[33]。因此, 考虑各类生态系统服务相关关系进而制定优化的退耕还林还草空间规划及相应方案, 以提高生态修复效果、减少负面影响^[11]。

可持续发展目标下, 我们面临着保障粮食供给的同时减弱其对生态环境负面影响的挑战。退耕还林还草一期工程始于1999年, 并于2014年启动新一轮工程, 林草植被恢复同时具有正、负效应, 对于权衡绿色生态空间和粮食生产空间的土地利用与管理优化带来更大的挑战^{[11, 25, 30]}。本研究拟解决以下问题：(1)退耕还林还草地块上生态系统调节或支持功能是否有所提升？(2)退耕还林还草工程对工程县域的粮食供给服务有何影响？(3)哪些区域的工程县域粮食供给损失最小、调节服务收益最大？通过量化生态系统水源涵养、土壤保持、防风固沙和粮食生产, 分析成对的生态系统服务之间的权衡与协同关系及其影响因素, 明确协同高、权衡低的退耕还林还草最佳目标县域, 探索哪些区域存在突出的权衡问题, 对于退耕还林还草工程的持续滚动实施、从空间角度提出平衡绿色和粮食的土地利用方案, 具有极为重要的科学意义。

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

以实施退耕还林还草一期工程的1897个县域作为工程区边界, 县域尺度生态系统变化反映气候变化背景下人类活动特别是生态修复工程的综合影响, 难以拆解出具体某个工程的作用。因此, 本研究基于中国科学院土地利用与土地覆被变化数据库, 提取县域内2000—2015年耕地转变为林地和草地的空间范围作为退耕还林还草工程地块, 统计分析地块上的水源涵养、土壤保持、防风固沙功能指标变化, 理论上可认为是退耕还林还草工程导致的地块变化。此外, 由于工程必然导致地块上粮食产量归零, 因此粮食供给服务以工程县域为研究单元, 分析工程对县域粮食供给服务的影响。工程县域按地理位置和地形地貌分为西南高山峡谷、川渝鄂湘山地丘陵、长江中下游低山丘陵、云贵高原、琼桂丘陵山地、新疆干旱荒漠、黄土丘陵沟壑、华北干旱半干旱、东北山地及沙地9个区域。

1.2 核心生态系统功能估算 1.2.1 粮食供给

工程县域的粮食供给以粮食产量统计数据为基础, 统计数据源于2000—2015年《中国县域统计年鉴》, 空间上根据归一化植被指数(NDVI)的时空分布特征将粮食产量分配至县域耕地地块得到粮食供给空间分布。首先, 基于2000—2015年1 km空间分辨率、16天时间分辨率的MODIS-NDVI数据, 通过格式转换、重投影、拼接、重采样和S-G滤波处理, 采用最大合成法得到连续时间序列的年最大NDVI数据, 然后叠加县域耕地分布得到耕地年最大NDVI空间分布, 最后依据NDVI值高低将县域粮食产量分配到每个栅格^[34]。

1.2.2 土壤保持

利用土壤保持量(SC, 单位：t/hm²)衡量生态系统减少降水导致土壤侵蚀的能力, 表征植被作用引起的水蚀减少量, 即实际地表覆盖条件下(M_sc)与极度退化状态下土壤水蚀量(M_sce)的差值。采用修正通用水土流失方程(Revised Universal Soil Loss Equation, RUSLE)估算地块尺度土壤水蚀模数：

(1)

(2)

式中, R是降雨侵蚀力因子, 采用基于日降雨量资料的半月降雨侵蚀力模型^[35]来估算；K为土壤可蚀性因子, 采用Nomo图法计算；L是坡长因子, S是坡度因子, 基于McCool等^[36]和Liu等^{[35, 37]}的方法计算；C是覆盖和管理因子, 依据蔡崇法等^[38]的方法计算；P为水土保持措施因子。

(3)

式中, OM为土壤有机质含量百分比(%), M_s为土壤颗粒级配参数, 即粉粒、极细砂与粘粒百分比之积, St为土壤结构系数, P_d为渗透等级, 土壤属性数据源于1∶100万中国土壤数据库。

(4)

(5)

(6)

(7)

式中, λ为坡长(m), θ为坡度, 计算L时把道路、河流、沟塘湖泊等地表要素作为径流的阻隔因素, 改进传统算法中通过相邻栅格间的坡向以及坡度变化率确定坡长终止点的方法, 避免坡长因子的高估。f为植被覆盖度, 利用NDVI数据和像元二分法计算得到。

1.2.3 防风固沙

利用防风固沙量(SF, 单位：t/hm²)衡量生态系统削弱风力、降低风蚀量的能力, 表征植被作用引起的风蚀减少量, 即实际地表覆盖条件下(M_sc)与极度退化状态下(M_sce)土壤风蚀量的差值。在充分考虑气候条件、植被状况、地表土壤的粗糙度、土壤可蚀性、土壤结皮的情况下, 利用修正土壤风蚀方程(Revised Wind Erosion Equation, RWEQ)估算地块尺度土壤风蚀模数^{[39, 40]}。

(8)

(9)

式中, Q_x表示x处的沙通量(kg/m), x表示地块长度, Q_max表示风力的最大输沙能力(kg/m), s表示关键地块长度(m)。

(10)

式中, WF表示气象因子；EF表示土壤可蚀性成分因子；SCF表示土壤结皮因子；K′表示土壤糙度因子, 反映地表抗风蚀能力, 即地表对风速减弱作用及对风沙活动的影响, 其大小取决于地表粗糙元的性质, 通过文献参数整理得到不同地表类型的土壤糙度因子值；COG表示植被因子, 是枯萎植被、直立残茬和生长植被覆盖的土壤流失比率的乘积。

(11)

(12)

(13)

式中, WS₂为2 m处风速(m/s), WS_t为2 m处临界风速；N为风速观测次数；N_d为试验天数；ρ为空气密度(kg/m³), g为重力加速度(m/s²)；SW为无量纲的土壤湿度因子；SD为雪覆盖因子, 即计算时段内积雪覆盖深度大于25.4 mm的概率。S_a为土壤砂粒含量, S_i为土壤粉砂含量, C_l为土壤粘土含量, OM为有机质含量, CaCO₃为碳酸钙含量。

1.2.4 水源涵养

利用降水贮存量法计算森林、草地生态系统水源涵养量(WR, 单位：m³), 以生态系统水文调节效应衡量其涵养水分的能力, 即与裸地相比涵养水分的增量^[41]。

(14)

(15)

(16)

式中, A为森林、草地和湿地生态系统面积(hm²), J为产流降水量(mm), J₀为年降水量(mm), K为产流降水量占降水总量的比例, R_e为与裸地相比, 生态系统减少径流的效益系数, R₀为产流降水条件下裸地降水径流率, R_r为产流降水条件下生态系统降水径流率。通过搜集已发表文献的实测降雨产流临界值, 以临近国家气象台站实测日降水数据修正同时期的TRMM逐日3小时降水量数据, 累积单次降雨量大于降雨产流临界值的数值得到单点产流降雨量占降雨总量的比例, 再与多年平均河川径流系数建立线性关系, 得到K值空间分布, 同时利用坡面径流实测数据确定了生态系统与裸地相比减少径流的效益系数。

1.3 生态系统服务权衡与协同关系分析

利用最小二乘法分析2000—2015年退耕还林还草工程县域粮食供给以及工程地块上土壤保持、水源涵养、防风固沙的年际变化趋势。

(17)

式中, Slp为年际变化率, X_i为第i年的功能量, i=1, 2, …, n。年际变化率反映2000—2015年某像元的总体变化趋势, 斜率为正表示年际变化趋势增加, 反之则为降低。

生态系统服务权衡与协同关系研究多基于相关分析、聚类分析、极值法、情景分析、图形比较和模型模拟等方法^{[4, 42]}, 从定性分析、区域差异、尺度特征、情景模拟预测、动态变化及其应用等方面^[43], 开展时间权衡、空间权衡、可逆权衡及多种服务之间的权衡^[44]。通过空间叠置分析, 基于逐像元空间自相关分析方法, 探讨两种服务间相关系数的绝对值大小及正负方向, 判断两两生态系统服务权衡与协同关系。

首先, 计算两类生态系统服务之间的相关系数、一阶偏相关系数和二阶偏相关系数。

相关系数计算公式如下：

(18)

一阶偏相关系数计算公式如下：

(19)

式中, r_ij·h表示在变量i和j控制变量h后的一阶偏相关系数, r_ij、r_ih、r_jh分别为两个变量的相关系数。

二阶偏相关系数计算公式如下：

(20)

式中, r_ij·hm表示在变量i和j控制变量h和m后的二阶偏相关系数, r_ij·h、r_im·h、r_jm·h分别为两个变量的一阶偏相关系数。

如果二阶偏相关系数为正, 表明在排除另外两个服务的影响下, 这两个生态系统服务之间为协同关系；相关系数为负, 则为权衡关系；相关系数为零, 则表明无相关关系。

根据偏相关系数的零假设检验t检验方法, 判断两两之间相互关系的显著性, 公式如下：

(21)

式中, r是对应的偏相关系数, n是样本观测数, k是可控制变量数, n-k-2是自由度, 当t > t_0.05(n-k-2)时即P < 0.05, 拒绝原假设, 相关性结果显著, 当t > t_0.01(n-k-2)时即P < 0.01, 拒绝原假设, 相关性结果极显著。t_0.05(n-k-2)和t_0.01(n-k2)查阅t检验表确定。

2 结果与分析 2.1 退耕还林还草对核心生态系统功能的影响

2000—2015年, 我国退耕还林还草地块的土壤保持量与水源涵养量均上升, 而防风固沙量则下降(图 1)。过去16年, 地块土壤保持功能以增强为主, 38%工程县域无变化, 46%工程县域的地块土壤保持量增加, 其中14%增势显著, 主要分布在黄土高原、秦巴山地；约16%工程县域的地块土壤保持量下降, 主要在云南省和湖北省。49%工程县域的地块水源涵养量增加, 其中仅约3%增势显著；约10%工程县域的地块水源涵养量下降, 主要在新疆、湖南、湖北等省域。风蚀区515个工程县域中, 26%工程县域的地块防风固沙量下降, 其中约46%降势显著；仅约2%工程县域的地块防风固沙量增加。工程县域粮食供给以增加为主, 76%工程县域有所增加, 其中16%增势显著, 约19%工程县域则下降。

不同区域比较(图 2), 多年平均土壤保持量最高值在西南高山峡谷区(172.2 t hm^-2 a^-1), 最低值在华北干旱半干旱区(8.2 t hm^-2 a^-1), 其年际增加趋势在黄土丘陵沟壑区(2.7 t hm^-2 a^-1)极为显著, 在长江中下游低山丘陵区、琼桂丘陵山地区、华北干旱半干旱区、东北山地及沙地区增势缓慢, 在西南高山峡谷区、云贵高原区、新疆干旱荒漠区呈现递减趋势。多年平均水源涵养量最高值在琼桂丘陵山地区(22.7 m³ hm^-2 a^-1), 最低值在华北干旱半干旱区(2.0 m³ hm^-2 a^-1), 各区域均呈增加态势, 特别是琼桂丘陵山地区(0.34 m³ hm^-2 a^-1)、云贵高原区(0.29 m³ hm^-2 a^-1), 而新疆干旱荒漠区、华北干旱半干旱区的增势缓慢且波动较大。多年平均防风固沙量最高值在新疆干旱荒漠区(40.7 t hm^-2 a^-1), 最低值在黄土丘陵沟壑区(8.63 t hm^-2 a^-1), 各区域均呈波动减少态势, 特别是新疆干旱荒漠区(3.75 t hm^-2 a^-1)。

2.2 退耕还林还草地块生态系统服务权衡与协同关系

从两两相关关系来看(图 3), 约28%工程县域粮食供给与地块水源涵养之间呈现显著相关性, 其中86%工程县域为显著协同关系(55%极显著), 14%为显著空间权衡关系。约21%工程县域粮食供给与地块土壤保持之间呈现显著相关性, 其中87%工程县域为显著协同关系(50%极显著), 13%为显著空间权衡关系。风蚀区515个工程县域中, 约44%工程县域粮食供给与地块防风固沙之间呈现显著相关性, 其中87%工程县域为显著权衡关系(34%极显著), 13%为显著空间协同关系。通过空间叠置分析, 约31%工程县域为多重生态系统服务供给区, 其中43%(292个县域)为双重生态系统服务供给区, 53%(357个县域)为三重生态系统服务供给区, 约4%(31个县域)为四重生态系统服务供给区。

2.3 退耕还林还草影响生态系统服务权衡与协同关系的时空分异特征

不同区域生态系统供给与调节服务的相关关系差异明显(图 4)。退耕还林还草导致工程县域粮食供给与地块土壤保持的相关关系呈现4种变化趋势：1)西南高山峡谷区、川渝鄂湘山地丘陵区、云贵高原区、新疆干旱荒漠区由权衡转变为协同；2)长江中下游低山丘陵区、华北干旱半干旱区的权衡关系轻微地减弱；3)黄土丘陵沟壑区、东北山地及沙地区的权衡增强；4)琼桂丘陵山地区是相互增益的协同关系增强。从波动幅度和离散程度来看, 新疆干旱荒漠区、西南高山峡谷区、川渝鄂湘山地丘陵区、云贵高原区均大于其他区域。

退耕还林还草导致工程县域粮食供给与地块水源涵养的相关关系呈现5种变化趋势：1)西南高山峡谷区、川渝鄂湘山地丘陵区、琼桂丘陵山地区的权衡减弱并趋向协同；2)长江中下游低山丘陵区的极显著权衡减弱；3)云贵高原区的权衡增强；4)新疆干旱荒漠区、华北干旱半干旱区的极显著协同减弱；5)黄土丘陵沟壑区、东北山地及沙地区从协同变为权衡而近几年又趋于协同。

退耕还林还草导致工程县域粮食供给与地块防风固沙的相关关系呈现3种变化趋势：1)新疆干旱荒漠区、东北山地及沙地区从权衡转变为协同；2)黄土丘陵沟壑区从协同转变为权衡；3)华北干旱半干旱区的相关关系在协同与权衡之间波动。

3 讨论

退耕还林还草工程增加了对于生态系统调节功能提高至关重要的林地、草地, 使得退耕还林还草地块上粮食供给服务呈现不同程度减少^{[26, 45]}, 也使得这些地块上供给与调节服务之间呈现极显著权衡关系。然而, 由于实施工程的耕地地块多属于水土流失、沙化、盐碱化、石漠化严重并且粮食产量低而不稳的陡坡耕地^[46], 因此, 从县域尺度粮食供给变化来看, 约80%工程县域粮食供给增加, 80%以上县域的粮食供给与地块水源涵养、土壤保持服务之间表现为显著协同关系, 说明大部分工程县域能够提供多重生态系统服务。

不同于水源涵养、土壤保持的上升趋势及其与供给服务之间的协同关系, 风蚀区防风固沙下降的工程县域数量远高于上升的县域, 约38%工程县域粮食供给与地块防风固沙之间呈现显著权衡关系, 仅6%为显著协同关系。分析原因发现, 近16年风场强度减弱特别是沙尘天气易发生的春季风场强度减弱, 以及对北方防风固沙服务最重要的草地和沙地植被恢复, 是北方土壤风蚀量下降的重要原因^[46]。通过与极度退化土地风蚀量对比得到的防风固沙量, 其大小取决于不同生态系统类型土壤风蚀量变化的相对大小。虽然, 春季植被盖度变化与防风固沙量呈显著正相关, 然而, 中度及以上植被覆盖增加的区域呈现植被覆盖变化与防风固沙能力的正相关关系, 中度以下植被覆盖度变化程度与防风固沙能力则呈现负相关^[48-49]。因此, 植被覆盖好转并不意味着防风固沙功能提升。

当然, 北方风蚀区退耕还林还草地块防风固沙功能下降进一步反映了干旱半干旱区开展不同生态修复措施的生态效应存在争议, 退耕还林、退耕还草、封育/封禁等都会产生差异效果^[28], 生态修复需遵循“宜草则草、宜荒则荒、宜林则林”原则, 以及多重服务权衡与协同关系评价的重要性。本文利用逐像元空间自相关分析来说明两两生态系统服务之间的关系, 然而, 多个生态系统服务权衡的检测仍然是有限的, 很少考虑联合权衡, 一个最优的、多目标的生态系统管理必须最大限度地提高协同、尽量减少权衡。因此, 还需要综合考虑气候调节、木材和纤维、基因和生物资源、病虫害调控乃至文化服务等, 在空间上对多个服务进行综合评价及权衡, 为区域生态修复与管理提供决策支持^[11]。

对于表现为权衡关系的区域, 长江中下游低山丘陵区、云贵高原区部分工程县域退耕还林还草地块的调节和供给服务均有不同程度的削弱, 此类区域成为水源涵养、防风固沙、土壤保持、粮食供给的低值中心, 同时也加大了粮食供给与生态需求的矛盾, 此为造成权衡程度较高的原因之一, 也由此说明退耕还林还草工程的区域适宜性和可持续性是工程持续开展必须考虑的重要因素。而对于表现为协同关系的区域, 即使近年已经取得了多重生态成效, 未来仍需加强生态修复与管理。

4 结论

退耕还林还草工程实施近十几年, 工程县域粮食供给与退耕还林还草地块的土壤保持、水源涵养、防风固沙之间权衡与协同关系的时空变化。工程县域地块土壤保持量与水源涵养量均上升, 而防风固沙量则下降。土壤保持年际变化在黄土丘陵沟壑区、川渝鄂湘山地丘陵区等区域呈现明显地波动增加态势, 而在西南高山峡谷区、云贵高原区、新疆干旱荒漠区等区域呈现递减趋势。水源涵养服务在各区域均呈现增加态势, 而防风固沙服务在各区域均呈现减少态势。不同区域退耕还林还草导致粮食供给与土壤保持、水源涵养、防风固沙相关关系呈现差异地变化趋势。工程县域粮食供给与地块防风固沙之间存在着显著的负相关性, 粮食供给与水源涵养之间、粮食供给与土壤保持之间存在着显著的正相关性。