武陵山区是我国集中连片特困地区之一, 中共中央、国务院及各级政府为武陵山区的发展从各方面提供政策支持。退耕还林工程作为我国“六大林业工程”之一, 同时也是我国实施西部大开发战略的重要举措^[1-2], 新一轮退耕还林工程自2014年启动以来, 武陵山区内各省退耕还林工程任务主要向片区倾斜, 坚持发展特色产业精准助推扶贫, 取得了很好的成效。

随着生态系统为人类福祉做出的贡献被普遍认可, 国内外大量关于生态系统服务的研究涌现, 部分学者从生态系统服务量化评估、时空演变、模拟预测以及权衡协同关系等多个角度进行了的研究^[3-7], 较为常用的生态系统服务价值量评估方法一般为实际市场法、替代市场法和虚拟市场法^[8]；目前生态系统服务在生态学领域已经有较为成熟的评估方法, 同时, 生态系统服务在景观规划领域也成为了研究热点, 有学者从农田、林地等不同景观尺度对生态系统服务进行研究, 在此基础上进行景观格局演变对生态系统服务价值影响的研究^[9-17]；此外, 还有学者针对景观破碎化与生态系统服务功能之间的联系提出了一个概念框架, 认为景观破碎化主要对生物多样性保护功能、涵养水源和森林游憩等功能产生正面或负面的影响^[18]。总体而言, 目前已有的研究主要侧重于生态系统服务本身的估算与演变, 即便是从景观角度进行分析, 也是对单一功能或单一景观进行研究, 针对生态系统服务与景观结构两者之间整体的相互作用机制的研究较少。因此本文以武陵山区退耕还林工程为例, 对其生态系统服务进行估算、景观格局演变进行分析, 在此基础上研究生态系统服务与景观格局耦合协调度, 一方面对片区实施退耕还林工程取得的生态效益情况进行总结, 另一方面揭示生态系统服务与景观格局在时间和空间上的动态关系, 填补生态系统服务功能和景观格局演变耦合协调研究的空白, 同时可为片区生态系统管理和土地利用资源合理有效配置提供决策辅助, 对人与自然和谐耦合发展具有重要意义, 对于全面实施长江经济带发展战略和加快推进生态文明建设具有时代意义。

1 研究区概况

武陵山区(图 1)是我国以少数民族聚集为主的集中连片特困地区之一^[19], 地处重庆、贵州、湖北、湖南四省交界地带, 包含64个县, 片区以武陵山脉为中心, 地形复杂, 山峦起伏, 小型盆地和丘陵遍布其中^[20], 片区总人口为2300多万人, 少数民族占总人口数量将近一半^[21], 受民族宗教、自然因素等影响, 片区长期处于相对封闭的状态, 经济发展落后。此外, 武陵山区位于我国重要的长江经济带, 同时是我国的重点生态功能区, 典型的亚热带动植物区, 生物多样性保护区, 长江流域水土保持、水源涵养以及生态安全屏障功能区, 是连接中原与西南的交通纽带, 中国区域经济的分水岭和西部大开发的最前沿^[22], 片区脆弱的生态环境和缓慢的经济发展对长江经济带整体的建设起到一定的阻碍作用, 基于此, 武陵山区经济发展以生态文明建设与扶贫开发相结合, 强调生态保护, 把修复生态环境放在首要位置。

2 研究方法 2.1 数据来源

(1) 由于退耕还林工程开始于1999年, 前期缺乏相关生态参数数据, 故本文选用2005年和2017年各省退耕还林工程资源数据进行生态系统服务功能估算, 该数据由国家林业和草原局退耕还林(草)工程管理中心提供(表 1), 退耕还林工程共有退耕地还林、宜林荒山荒地造林和封山育林3种植被恢复模式, 本文以3种模式进行数据展示和计算, 需要注意的是, 退耕还林工程并不是武陵山区唯一的生态工程, 本文以退耕还林工程为例, 进行生态系统服务功能与景观格局之间的耦合协调关系, 不涉及其他林业工程；

(2) 生态参数数据来源于武陵山区内及周边的10个森林生态站和4个退耕还林工程专项及辅助监测站, 以及多个辅助观测点和固定样地(见图 2)；

(3) 1995、2005、2015年3期土地利用数据从中国科学院资源环境科学数据中心(http://www.resdc.cn/)下载。

2.2 研究方法 2.2.1 生态系统服务功能评估

本研究选取保育土壤、林木养分固持、涵养水源、固碳释氧、净化大气环境5项主要生态服务功能, 参照中华人民共和国国家标准《森林生态系统服务功能评估规范》(GB/T 38582—2020)和《退耕还林工程生态效益监测与评估规范》(LY/T 2573—2016)和《森林生态系统生物多样性监测与评估规范》(LY/T 2241—2014)两项行业标准, 采用分布式测算方法对武陵山区退耕还林工程生态系统服务功能进行评估, 即对武陵山区每个县, 按退耕还林工程3种植被恢复模式的优势树种(组)进行计算, 相关计算公式详情请参照文献^[23-27]。

2.2.2 景观格局分析

景观格局分析基于1995、2005、2017年3期的景观类型分布图, 景观核密度分析结果利用ArcGIS的核密度分析工具获得^[28-29], 景观格局指数利用fragstats专业软件获得, 各指标含义参照文献^[30-31]。

2.2.3 耦合协调度模型

(1) 指标无量纲化

考虑生态系统服务和景观格局变量在水平和量纲上的影响, 采用极差标准公式, 对数据进行标准化, y由于指标具有正、负两种效应, 所以正效应指标选用公式：, 负效应指标选用公式：y′_ij=其中i为评价样本, j为指标, y_ijmax和y_ijmin分别表示系统第j项指标在所有样本中的最大值和最小值。

(2) 嫡权法确定权重模型

熵权法是利用各指标值的差异程度来计算各指标的权重。第j项指标占所有样本的比重P_ij：P_ij= , 之后, 求指标的熵, 最后计算评价指标的熵权。

(3) 耦合协调度

根据综合评价函数：计算生态系统服务和景观格局的综合评价指数U₁和U₂, 当U₁>U₂, 景观格局分布较差, 反正, 则生态系统服务较差。

耦合关联度函数为：, 此外还需计算综合调和指数, T=αU₁+βU₂, α和β分别代表生态系统服务和景观格局的权重, 由于二者同样重要, 所以α和β均取值0.5；最后计算耦合协调度。

根据已有研究, 当C =0时, 无耦合, 当D∈(0, 0.3]时, 为低度耦合, D∈(0.3, 0.8]时, 为中度耦合, D∈(0.8, 1)时, 为高度耦合。当D∈(0, 0.3]时, 为低度协调耦合, D∈(0.3, 0.5]时, 为中度协调耦合, D∈(0.5, 0.8]时, 为高度协调耦合, D∈(0.8, 1)时, 为极度协调耦合^[32-33]。

3 结果与分析 3.1 退耕还林工程生态系统服务功能物质量

武陵山区退耕还林工程生态系统服务功能物质量如表 3所示, 2005—2017年间武陵山区退耕还林工程各项生态服务功能物质量均增长程度均在15倍左右, 其中林木养分固持物质量增长程度最高为16.75倍, 其次为释氧和固碳功能, 增长程度均在16倍之上。整体来看, 武陵山区退耕还林工程2005—2017年间, 生态系统服务功能发挥稳定增长, 对于改善生态环境, 保护生态安全具有重要作用。

以武陵山区2005年和2017年涵养水源、固碳和滞尘功能为例, 分析片区退耕还林工程生态系统服务功能物质量空间分布特征。由图 3可知, 2005年三项功能整体空间分布格局基本遵循由中间向南北两个方向逐渐降低, 且片区南部整体大于北部的趋势；2017年, 片区内各县三项功能均有较大增长, 空间分布格局呈现生态系统服务功能物质量较高的地区逐渐由西部向东南方向蔓延, 分布较2005年更为均衡。

3.2 景观格局变化空间集聚特征分析

1995—2005年10年间, 武陵山区景观变化主要呈片状集中分布在西部地区的重庆市和铜仁地区(如图 4), 景观格局变化最大值为148.086。2005—2015年间, 武陵山区景观变化景观变化向东南部扩散, 主要集中在片区内重庆市、铜仁地区、娄底市、和邵阳市, 恩施土家族苗族自治州和湘西土家族苗族自治州次之, 景观格局变化最大值为38.2869, 与1995—2005年间对比景观变化不太明显, 景观格局变化空间上更为分散。从景观格局变化的主要区域来看, 20年间, 景观格局变化主要发生在坡度大于25°的区域, 这与退耕还林工程实施情况极为相似, 坡度大于25°耕地是退耕还林工程实施主要区域之一, 同时, 工程实施期间2005年前资源面积增长幅度较高, 景观格局变化最大值也较大, 反映了退耕还林工程对景观格局的影响。此外, 武陵山区20年间景观格局演变主要分布在林地、耕地和草地之间, 其中林地和耕地之间的转移率最高, 1995—2005年间, 三者相互转化占整个转移网络的86.91%, 2005—2015年间, 三者相互转化占整个转移网络的76.52%, 表明武陵山区景观格局变化主要受这3种景观类型影响。

3.3 景观格局指数变化分析

在景观尺度上进行景观格局指数的计算, 如表 4所示。

从景观水平指数上可以看出, 1995—2015年武陵山区的景观斑块数量和斑块密度先减少后增加, 增幅不太明显, 蔓延度指数逐年减少, 聚集指数先增加后减少, 表明景观中存在大量破碎的小板块, 连接性在下降, 原因一方面为城市扩张影响, 另一方面与退耕还林工程有关, 增加的林地和减少的耕地造成了景观的破碎化, 影响到景观格局中原有优势景观及其连通性, 但在类型水平上计算时, 2005—2015年间, 耕地和林地的蔓延度指数和聚集度指数均呈现增加趋势。香浓多样性指数和香浓均匀度指数整体呈增加趋势, 表明片区土地利用景观丰富度和复杂度不断升高, 景观类型逐渐趋向均匀分布^[34], 不同类型对景观整体格局的影响作用逐渐平衡^[35]。总的来看, 随着林业工程的实施, 武陵山区土地利用景观类型呈现出复杂性和离散性趋势, 这表明武陵山区当前还没有出现优势景观, 景观中存在大量的小斑块, 在下一步生态工程建设中, 更应注意景观格局的连接度以及生态网络的流通性, 对武陵山区的景观格局进行优化。

3.4 生态系统服务于景观格局耦合协调度分析

从上述结果可以发现, 武陵山区退耕还林工程生态系统服务和景观格局的时空演变存在一定的相似性, 为了进一步阐明二者之间的关系, 使用耦合协调模型对二者关系进行测度, 并揭示其协调程度。由于武陵山区退耕还林工程植被恢复模式包含封山育林, 即退耕还林工程开始实施前, 已有部分林分在发挥生态系统服务, 但是由于缺乏工程开始前的相关生态参数, 无法进行准确的计算, 故生态系统服务选取武陵山区内各市级行政区2005—2017年间的差值进行计算, 相应的, 景观格局指数选取2005—2015年间的差值进行耦合协调度分析, 生态服务功能选取涵养水源、固碳、滞尘和生物多样性保护功能, 景观格局选取NP、PD、CONTAG、SHDI、SHEI、AI, 结果如下表所示：

生态系统服务和景观格局综合评价指数, 在各个地市中存在差异, 多数城市为景观格局分布情况较好, 生态系统服务有待提升, 只有重庆和铜仁地区的景观格局分布情况有待进一步进行规划。每个市级行政区生态系统服务与景观格局的耦合程度均为中度耦合, 表明二者在武陵山区有较强的相关性, 生态互动关系较好, 但生态系统服务和景观格局分布还有提升和优化的空间。协调度都在高度协调之上, 且多数为极度协调, 表明武陵山区生态系统服务和景观格局10年间在发展程度上步调一致, 同时, 景观格局发展受生态系统服务影响较大。

由于退耕还林工程林分幼林龄、中林龄占比偏大, 所以各市退耕还林工程生态系统服务功能单位面积物质量与市区内所有森林生态系统服务功能单位面积物质量相比偏小。以第八次全国森林资源连续清查生态系统服务功能评估结果的单位面积物质量为依据, 合理放大生态系统服务所选指标, 以分析生态系统服务功能各指标对生态系统服务与景观格局的耦合协调度的影响。结果表明：涵养水源、固碳、滞尘和生物多样性保护功能中, 12个地市耦合协调度变化程度排名为固碳(-0.06)、涵养水源(0.05)、滞尘(0.02)、生物多样性保护功能(-0.01), 其中, 随各功能增加耦合协调度增加排名为：涵养水源(0.06)、滞尘(0.3)、固碳(0.02)、生物多样性保护(0.01)；各地市中, 以娄底市变化情况最佳, 随涵养水源和固碳功能增加, 生态系统服务与景观格局的耦合程度从中度协调均提升为极度协调, 娄底市在下一步退耕还林工程建设中, 应注重涵养水源和固碳功能的提升。

4 讨论

在生态文明建设已经成为既定国策的今天, 武陵山区作为我国重点生态功能区, 生态恢复和生态保护将直接关系到国家的生态安全^[36]。同时, 武陵山区作为集中连片特困地区之一, 是片区内各省份新一轮退耕还林工程的主要实施区域, 片区内64个县均实施了退耕还林工, 对片区生态文明建设起到了一定的成效。截至2017年底, 武陵山区退耕还林工程每年涵养水源35.09×10⁸m³, 相当于三峡水库蓄水深度达到15.64m时的库容量；固土4712.45×10⁴t, 是2014年长江(2.75×10⁸t)土壤侵蚀量的17.14%^[37], 水土保持生态服务价值占总价值的43%, 处于主导地位, 与退耕还林工程建设的初衷一致。生物多样性保护价值排在第三位, 占总价值量的19.48%, 与生态区位需求基本吻合。从2005—2017年12年间, 各项生态系统服务功能以每年约1.25倍的速度增长, 为武陵山区生态修复提供了保障。从现有研究中可以发现, 影响生态系统服务功能大小的最主要原因是资源面积的大小, 水热、地形以及经营等同样会影响生态系统服务功能^[23-26], 随着工程的进一步实施, 退耕还林工程成果的进一步巩固, 未来退耕还林工程在武陵山区将会发挥更高的生态系统服务功能, 为区域水土保持、生物多样性保护, 以及维护国家生态安全发挥巨大作用。

利用武陵山区生态系统服务功能空间分布图和景观格局核密度分析结果对比发现, 20年间景观格局空间变化情况与生态系统服务功能的空间变化基本一致, 均为由西部地区向东部、东南部、东北部呈发散状进行转变, 坡度较大的区域景观格局变化程度也较大, 同时, 景观格局变化区域主要集中在林地、耕地和草地之间。景观格局不仅是景观异质性的表现, 也是生态系统的生态过程的表现^[38], 景观格局通过影响生态过程, 从而对生态系统服务功能造成一定的影响^[39]。武陵山区20年间, 林地的斑块形状越来越规则, 呈现集中分布的趋势, 林地面积大幅增加, 片区林业工程实施情况较好, 对生态系统服务功能的提升具有一定的帮助；但重要的生态源地的水体景观越来越破碎, 作为片区内林业工程实施的重要保障, 在片区生态环境建设中应进一步加强水体景观建设, 保证生态源地的稳定, 增加水体景观的连接度, 形成水体网络^[29], 进而提高生态系统服务功能。

使用耦合协调度模型得到的结果, 证明武陵山区的退耕还林工程生态系统服务功能和景观格局存在较好的耦合关系和极强的关联关系。武陵山区整体处在生态系统服务提升和景观格局演变协调度较高的水平, 二者发展情况一致性较高, 相互影响较大, 多数市的景观格局趋向正向演变, 生态系统服务功能提升较为滞后, 重庆市和铜仁地区恰好相反, 生态系统服务功能提升幅度较高, 但景观格局有待优化；此外, 生态系统服务功能指标对耦合协调度的影响因子以固碳和涵养水源排在前列, 随后是滞尘和生物多样性保护功能。固碳排在首位主要与退耕还林工程植被恢复模式有关, 由表 1可知宜林荒山荒地造林和退耕地还林占比最大, 工程区林分多为幼龄林或中龄林^[40], 处于生长期、生长季的林分固碳功能较高^[41-42], 退耕还林工程固碳功能与工程实施情况关系密切, 随着新造林面积增加, 固碳功能也会增加, 但景观格局势必会受到影响, 尤其在新的景观格局形成的初期, 景观格局指数各项指标会出现一定的波动；涵养水源功能排名靠前与退耕还林工程实施的主要目的和武陵山区生态区位有关, 工程实施主要为缓解我国当前所面临的水土流失问题, 涵养水源对景观格局起到正面影响, 与工程建设的初衷一致, 达到了预期目的, 同时, 武陵山区是长江流域水土保持、水源涵养以及生态安全屏障功能区, 随涵养水源功能的增加, 除重庆外, 各地市耦合协调度均呈增加趋势, 片区的涵养水源功能还需进一步进行提升, 而重庆市由于其特殊的地形条件, 目前景观格局演变与涵养水源功能相比, 处于劣势地位, 景观分布情况较差, 涵养水源功能提高, 这种差异也将被放大, 景观格局有待优化。滞尘排名较为靠后主要是由于退耕还林工程一般远离城市, 大气颗粒物却离城市很近, 二者在空间距离上明显存在差异。生物多样性保护功能排名最后, 主要是由于武陵山区作为生物多样性保护的重点生态功能区, 生物多样性保护工作进行的较为到位, 也很难再发生较大的改变, 所以生物多样性保护对武陵山区退耕还林工程与景观格局耦合协调度影响不大。

退耕还林工程实施保持情况较好, 生态系统服务稳步提升, 但提升幅度相较于景观格局演变情况还有很大的提升空间, 主要是由于退耕还林工程以幼龄林和中龄林为主, 林分在很长的一段时期内多数生态系统服务功能都处在较低水平, 随着时间的推移以及经营等手段的实施, 未来生态系统服务功能将有更大的提升。同时, 景观格局在演变过程中, 分布越来越均匀、形状越来越规则, 景观丰富度和复杂度不断升高, 其中林地和草地受退耕还林工程影响, 呈现集中分布趋势, 向着优势景观发展, 但在武陵山区下一步生态建设时应注意增强各景观类型的连接度, 加强生态网络的建设, 增加区域内物质和能量的流动。

武陵山区实施退耕还林工程, 有效提高了生态系统服务功能和改善景观格局。在今后武陵山区发展建设过程中, 可以结合退耕还林工程, 进一步优化景观格局, 提高二者耦合协调度, 规划以保护生态环境和资源可持续发展为目的的景观格局, 对区域经济发展和长江流域经济带建设都将带来巨大收益。

5 结论

本文利用分布式测算方法对武陵山区退耕还林工程生态服务功能进行评估, 与景观格局演变情况进行分析, 主要得出以下结论：

(1) 2005—2017年间武陵山区退耕还林工程生态服务功能时空变化为：时间上, 各项生态服务功能物质量增长程度均在15倍左右, 空间上, 分布格局最佳区域由西部地区逐渐向东南扩大；

(2) 武陵山区景观格局的时空变化以林地、耕地、草地为主, 变化剧烈区域由西北地区逐渐向东南方向蔓延, 景观格局的连通性下降, 抗干扰能力提升, 耕地类型和林地类型的聚集程度和连通性增加, 有形成优势景观的趋势, 与武陵山区退耕还林工程生态服务功能时空演变存在一定的相似性；

(3) 武陵山区退耕还林工程成效明显, 与景观格局演变具有较高的耦合协调度, 以固碳和涵养水源功能对生态系统服务功能与景观格局演变耦合协调度影响最高, 片区内以水体为主的景观有待进一步优化。