2008年的“5.12”汶川大地震对区域生态系统造成极大的破坏, 灾区总面积约50万km^2[1-2], 受损生态系统的空间分布与地震烈度高度相关, 在9级烈度区域就导致65584 hm²大熊猫生境丧失, 97748 hm²森林生态系统受损^[3]。地震造成植被破坏、生物栖息地受损、水土流失加剧等问题^[2-6]。在震区生态恢复研究方面, 杨斌等^[7]对汶川县震后5年生态环境质量恢复评价研究的结果表明2013年汶川县55.84%的区域环境质量得到了明显改善；杨渺等^[8]对汶川地震极重灾区10个县市震后5年的生态服务功能评估结果表明区域生态服务功能得到了一定程度的恢复, 但未达到震前水平。震后10年来, 不同地震烈度区的生态恢复效果如何？生态系统服务功能和生态系统健康如何变化？这是对震区后期生态恢复策略制定具有重要意义的关键科学问题, 已有的研究对不同地震烈度级别区的生态恢复效果缺乏综合比较。

生态系统健康是描述相对于理想或目标状态的生态系统状况, 是保证生态系统功能正常发挥的前提^[9-10], 尤其是在强烈依赖自然资源的西南震区, 生态系统的健康程度直接影响震区生态安全和可持续发展。因此, 本文首先利用遥感数据识别震后恢复期生态系统类型、格局以及质量变化, 然后运用InVEST(Integrated Valuation of Environmental Services and Tradeoffs, InVEST)模型定量评估关键的生态系统服务功能格局和变化, 最后耦合灾区生态系统活力、组织力以及关键的生态系统服务, 综合定量评估不同地震烈度区的生态系统健康状况, 以期为现阶段制定震区生态恢复策略提供科学依据。

1 研究区

本研究范围为四川省的汶川地震极重灾区和重灾区, 包括39个县市区(图 1, 地震烈度数据来自中国地震局), 总面积9756025.56 hm²。该区域位于青藏高原向成都平原的过渡地带, 地貌类型多样, 海拔在311—6117 m之间, 分布多种植被类型, 岷山和邛崃山系是全球生物多样性保护热点与关键地区^[11-12], 该区域也是岷江、沱江、嘉陵江、涪江的发源地。汶川地震受损生态系统的空间分布与地震烈度高度相关, 主要分布在地震烈度Ⅹ度及以上区域^[3], 绵竹、汶川、安县和什邡等是受损最严重的区域。

2 数据与方法 2.1 数据来源

本研究采用的数据及其来源为：(1)2010年、2015年地震重灾区生态系统类型(90 m)、归一化植被指数NDVI、植被覆盖度栅格数据(250 m), 来源于国家重点研发计划“西南生态安全格局形成机制及演变机理”, 将研究区生态系统类型一共划分为8个一级和18个二级类(表 1)；(2)DEM数据(90 m), 来自中国科学院地理空间数据云平台(http://www.gscloud.cn/)；(3)1 : 100万土壤类型数据, 来自中国土壤数据集(http://westdc.westgis.ac.cn)；(4)四川、甘肃、陕西、重庆等80个气象站点1990—2015年的日值降雨数据, 来自中国气象局气象中心, 通过预处理、Kriging空间插值得到研究区年均降水量、月均降水量空间分布栅格数据；(5)月参考蒸散发栅格数据(1 km), 来自CGIAR CSI dataset(http://www.cgiar-csi.org/data)；(7)四川省自然保护区数据, 来自全国及各地区主体功能区划。

2.2 生态系统健康评价方法

生态系统健康是描述相对于理想或目标状态的生态系统状况, 是保证生态系统功能正常发挥的前提, 与人类的可持续生计和可持续发展密切相关^{[9-10, 13]}。结构和功能的完整性、具有抵抗干扰和恢复能力、稳定性和可持续性是生态系统健康的特征^[14-15]。Rapport等^[10]提出从生态系统活力—组织性—弹性(Vigor-Organization-Resilience model, VOR)三个层面来评价生态系统健康得到了广泛应用, 但是区域尺度的生态系统健康评价需要基于生态系统的结构、功能过程来确定指标^[16-18]。本研究基于VOR模型, 耦合灾区关键的生态系统服务, 定量评估灾区震后恢复期的生态系统健康状况, 方法如下：

(1)

式中, EHI_i表示栅格i的生态系统健康指数, V_i、O_i、R_i、S_i分别表示栅格i的生态系统活力指数、组织力指数、恢复力指数、生态系统服务指数。

生态系统活力V指根据营养循环和生产力所提供的能量^[19], 一般用植被生产力表示^[20]。由于归一化植被指数(NDVI)与植被生产力高度正相关^[21-22], 本文选取研究区植被生长季(5—9月)的NDVI作为生态系统活力指标。

组织力O指生态系统物种组成结构及其物种间的相互关系, 反映生态系统结构的复杂性。本研究选取景观水平上的蔓延度指数(CONTAG)、Shannon多样性指数(SHDI)和斑块类型水平上的森林聚集度指数(COHESION)三个指标的综合值来衡量生态系统组织力。蔓延度指数指斑块类型在空间分布上的集聚程度, 该指标单位为%, 取值范围为(0, 100], 指标值越大, 表明斑块集聚程度越高。Shannon多样性指数在计算生态群落多样性时应用十分广泛, SHDI无量纲, 取值范围为SHDI≥0。聚集度指数度量的是相关斑块类型的自然联通度, 该指标没有单位, 取值范围[0, 100), 值越小表明斑块类型分布越分散。在FRAGSTATS 4.2中计算各县2010年和2015年的CONTAG、SHDI、COHESION指数, 分别标准化到[0, 1]之间, 当年的生态系统组织力指数由标准化处理之后的CONTAG、SHDI、COHESION均值得到。

恢复力R指生态系统在胁迫下维持其结构和功能的能力。常用生态系统弹性作为恢复力评价指标。一个健康的生态系统, 对外部的干扰具有抵抗力和弹性, 使得生态系统的组织结构和生产力在很长一段时间内保持相对稳定状态^[23]。生态系统系统弹性大小与地形、气候、土壤以及植被状况密切相关^[24], 在参照相关文献^{[18, 25-26]}的基础上确定研究区各生态系统的弹性系数(表 2)。

一个健康的生态系统应具备良好的生态系统服务供给能力^[27]。由于震区是岷山-横断山重要生态功能区和生物多样性保护关键地区, 本文选取土壤保持、水源涵养、生境质量等三类生态系统服务作为生态系统服务评估内容, 采用InVEST模型v3.4.4版本进行评估。土壤保持服务采用SDR(Sediment Delivery Ratio, SDR)模块计算, 主要参数包括土壤可蚀性因子、降雨侵蚀力因子、植被覆盖和经营管理因子和水土保持措施因子。水源涵养服务采用SWY(Seasonal Water Yield, SWY)模块, 主要参数包括月均降雨量与月均降雨次数、月参考蒸散发量、土壤水文分类、各植被类型的月作物系数Kc和径流曲线数CN值, 本文用结果数据中的基流量(Baseflow)来衡量研究区的水源涵养能力。生境质量采用HQ(Habitat Quality, HQ)模块, 主要参数包括各期胁迫因子图层(本文将水田、旱地、居住地、工矿交通用地作为胁迫因子)、胁迫强度、最大胁迫距离、各生境类型适宜度及其对胁迫因子的敏感度。各模块所需的部分参数结合研究区特征并参考相关文献^[28-38]确定。每个栅格的生态系统服务指数由标准化后的土壤保持、水源涵养、生境质量图层求均值得到。

由于各评价指标的量纲各异, 无法直接进行生态系统评价, 将各指标标准化处理到0—1之间, 最后基于栅格单元在ARCGIS10.3中计算每个栅格的生态系统健康指数。生态系统健康指数(EHI)评价结果在栅格图层上呈现0—1之间连续变化的值, 值越靠近1, 生态系统就越相对完整和健康, 将当年的EHI评价结果划分为[0, 0.2)、[0.2, 0.4)、[0.4, 0.6)、[0.6, 0.8)、[0.8, 1]等5个区间。

3 结果与分析 3.1 生态系统类型变化

震区生态系统类型分布差异明显, 森林、灌木、草地、冰川积雪、裸地主要分布在龙门山地震断裂带核心区域和西部高山地区, 农田和城镇主要分布在龙门山山前平原和与龙门山山脉接壤的浅丘地区。2010—2015年震区各生态系统变化如表 3所示。

从面积变化上看, 总体上森林、灌丛、草地、农田、裸地的面积有不同程度的减少, 森林面积减少2975.94 hm², 减幅为0.08%；农田面积减少31362.39 hm², 减幅为1.33%；湿地和城镇面积分别增加7605.9、31268.43 hm², 由于湿地和城镇基数较小, 增幅分别为7.72%、31.35%；其他类型由于基数较大所以减幅不大。其中, 28240.65 hm²的农田和1286.28 hm²的森林转移为城镇用地, 且主要发生在Ⅶ—Ⅷ级烈度区；然而在Ⅹ—Ⅺ级烈度区, 森林面积增加1309.77 hm², 说明地震最严重地区内森林面积的恢复比较明显。

3.2 生态系统健康变化

各指标要素计算结果见表 4和图 2、图 3和图 4。标准化后的各要素指标值和生态系统健康指数(EHI)结果见表 5, 震区EHI空间评价结果见图 2。将两期生态系统活力、EHI结果图层分别做差值分析来评估2010—2015年的变化格局, 差值按照[-1, -0.1)、[-0.1, -0.01)、[-0.01, 0.01)、[0.01, 0.1)、[0.1, 1]区间划分为明显下降、稍微下降、基本不变、稍微升高、明显升高等5个级别。

3.2.1 生态系统活力变化

从2010—2015年的变化来看(表 4), 灾区整体的生态系统活力均值从0.7194提高到0.7587, 增幅为5.46%；各级烈度区活力均值也都明显上升, 其中Ⅹ—Ⅺ级烈度区活力增幅最大, 分别为14.84%、12.47%。从空间格局来看(图 2), 西南-东北走向的林灌山区生态系统活力最高, 东南和中部平原地区次之。从分级面积来看, 活力呈稍微升高、明显升高的面积占比分别为55.48%、12.96%, 其中明显升高的区域分布在龙门山地震断裂带核心带, 集中在安州、绵竹、什邡、彭州、都江堰、崇州和江油的西部、大邑的西北部以及汶川的东南部；少数局部地区活力下降, 面积占比达17.66%, 主要分布在宝兴、芦山、汉源、石棉、汶川西北部和北部部分区域。总的来看, 虽然2010—2015年灾区各植被类型的面积在不同地震烈度区有不同程度的增减, 但是大部分区域的生态系统活力得到了提高。

3.2.2 生态系统组织力和恢复力变化

2010—2015年整体上生态系统组织力变化较小(表 4)。在Ⅶ—Ⅺ级烈度区, 景观蔓延度指数比稍有降低, Shannon多样性指数稍有升高, 表明灾后的生态恢复和开发建设过程在一定程度上提高了景观破碎度, 同时丰富了景观多样性。从空间格局来看(图 3), 西部高山林灌区组织力最高, 东部浅丘地区次之, 中部平原由于城镇化水平高、生态斑块破碎导致组织力最低, 其中涪城区的生态系统组织力有轻微下降。各级烈度区生态系统恢复力基本无变化(表 4), 空间格局上表现为西部高山林灌区恢复力最高(图 3), 东部浅丘和中部平原次之, 西部裸岩、冰川积雪地区最低。

3.2.3 生态系统服务变化

从水源涵养服务空间格局来看(图 4), 震区西南部的冰川积雪和裸岩区域水源涵养能力最高, 中部和东部平原次之, 西北高山林区水源涵养能力最低, 由SWY模块计算原理和震区自然状况推断这可能与西北高山林区降雨量偏少、森林的作物系数和蒸散量较大有关。从2010—2015年的变化来看(表 4), 震区的水源涵养能力稍有降低, 水源涵养总量减少了0.42亿m³, 减幅为0.46%；单位面积生态系统水源涵养能力均值减少了425.93 m³/km²；从不同地震烈度级别区来看, Ⅹ—Ⅺ级烈度区的单位面积生态系统水源涵养能力减幅分别为0.29%、0.67%, 尽管Ⅹ—Ⅺ级烈度区植被生态恢复效果较好, 但是随着NDVI的提高, 植被的作物系数增大进而导致植被作物需水量和蒸散量增加, 所以其水源涵养能力降低；Ⅷ—Ⅸ级烈度区单位面积生态系统的水源涵养能力稍有提高, 增幅分别为0.81%、0.40%；Ⅵ—Ⅶ级烈度区单位面积生态系统的水源涵养能力分别降低了1.43%、0.72%。

土壤保持服务明显表现为从东南向西北递增的空间格局(图 4), 中部和东部的城镇、耕地区域土壤保持能力最低, 西部和北部高山林灌区由于植被覆盖度较高, 因此土壤保持能力较高。2010—2015年, 震区土壤保持总量增加了0.9亿t(表 4), 增幅为1.26%, 单位面积生态系统土壤保持能力均值提高了9.25 t/hm², 主要受生态系统类型和NDVI的影响；同时, 各级烈度区单位面积生态系统土壤保持能力均有所提高, 其中Ⅺ级烈度区单位面积生态系统土壤保持能力最高, 约为震区平均水平的2倍, 其单位面积生态系统土壤保持能力相较2010年提高了1.52%, 这得益于Ⅺ级烈度区森林面积和植被覆盖度的提高。

在生境质量方面, 不考虑西部的冰川积雪和裸岩(生境适宜度为零), 震区生境质量明显表现为从东南向西北递增的空间格局(图 4), 中部和东部的城镇、耕地区域生境质量较低, 南部和北部高山林区由于森林、灌丛类型占主导并且距离建设用地、农田等胁迫因子距离较远, 因此生境质量较高。2010—2015年震区单位面积生态系统的生境质量指数无明显变化, 各级地震烈度区生境质量指数也无明显变化。

3.2.4 生态系统健康变化

2010—2015年震区整体的生态系统健康有小幅度提升, EHI均值从0.6761上升到0.6853(表 5), 增幅为1.36%。从不同地震烈度级别区来看, 2000年Ⅸ级烈度区EHI均值最大, 为0.7467, 而2015年Ⅹ级烈度区EHI均值最大, 为0.7613, 主要原因是Ⅹ级烈度区森林面积增加较多, 且生态系统活力提升最明显；同时, 各级烈度区EHI均值均有不同程度的提升, 其中Ⅹ—Ⅺ级烈度区生态EHI增幅比其他地区明显, 增幅分别为4.01%、3.48%。

从空间格局来看(图 2), 西部高山林草区EHI最高, 东南部浅丘和中部山前平原的农用地区次之, 中部的城镇和西部冰川积雪、裸岩的EHI最低。从分级面积来看, EHI稍微升高和明显升高的区域面积占比分别为37.32%、5.70%, 明显升高区域主要集中在绵竹、什邡、都江堰、崇州、大邑县的西北部等龙门山地震断裂带核心区域；42.23%的区域EHI基本不变；同时局部地区EHI下降, 面积占比达14.75%, 主要集中在宝兴、汉源、石棉、涪城、广汉、汶川的西北部和震区北部部分区域。分析其原因, 主要由于这些区域重建过程中森林和灌木的减少以及城镇和裸地的增加, 导致生态系统活力和生态服务功能降低。其中宝兴、汉源、石棉的森林和灌丛在2010—2015年期间共减少了623.70 hm², 城镇和裸地共增加了503.01 hm², EHI均值减幅为0.6%, 活力指数、组织力指数、生态系统服务指数均值减幅分别为2.12%、0.08%、0.05%, 应该重视这些区域未来建设过程中的生态保护工作。此外, 从表 5中可以比较明显的看出各级地震烈度区标准化后的组织力指数、恢复力指数、综合生态系统服务指数变化极其微弱, 而活力指数和EHI变化较大；同时, 从图 2可知, 生态系统活力明显上升和退化的区域分布格局与EHI的变化格局趋向一致, 在本文的生态系统健康评价框架下, 表明生态系统活力是影响震区生态系统健康状况的主要因素。

4 结论与建议

(1) 2010—2015年震区的森林、灌丛、草地、农田、裸地的面积有不同程度的减少, 湿地和城镇面积有所增加。减少的农田主要转移为城镇用地, 且主要发生在Ⅶ—Ⅷ级地震烈度区；森林面积共减少2975.94 hm², 但是在Ⅹ—Ⅺ级烈度区域内, 森林面积共增加1309.77 hm², 说明受灾最严重地区森林面积的恢复状况比较好。

(2) 震区生态系统活力均值从0.7194提高到0.7587, 增幅为5.46%；Ⅹ—Ⅺ级烈度区生态系统活力均值增幅最明显, 分别为14.84%、12.47%。就生态系统服务而言, 震区生境质量指数空间格局无明显变化；土壤保持总量增幅为1.26%；水源涵养量稍有降低, 减幅为0.46%。

(4) 震区EHI均值有小幅度提升, 从0.6761上升到0.6853, 增幅为1.36%；同时, Ⅹ—Ⅺ级烈度区EHI增幅比其他地区明显, 增幅分别为4.01%、3.48%, 表明受灾最严重的地区生态恢复比较显著。但是, 局部地区如宝兴、汉源、石棉、汶川的西北部和涪城、广汉、震区北部部分区域由于灾后城镇的扩张、森林灌丛和农田减少、生态系统活力降低导致其EHI下降。

未来应该持续关注汶川地震重灾区的生态恢复效应, 结合自然恢复和人工干预, 加快生态恢复进程；重视生态系统健康状况下降的局部地区的生态保护工作, 构建长期的震区生态恢复监测网络；同时, 在震区中部以农田和城镇为主的人类活动密集的平原地区, 重建过程中要注意保护自然生态空间和基本农田, 实施人工林草恢复时可通过合理搭配植被种类以及优化景观空间格局来提升关键生态服务功能和区域的生态系统健康。