经济的快速发展给资源与环境带来了影响, 生态环境问题日益显现。自然保护区的建立对我国自然资源及生态环境保护有着举足轻重的作用^[1]。国家公园是维持生态系统的原真性和完整性, 保护面积较广且具有国内代表性的生态环境系统, 在我国自然保护区体系中居主要地位^[2]。近年来, 我国国家公园在动植物生境保护方面取得实质性进展, 但不可否认的是历史及现存的资源采集、农业开垦、居民生活等人类扰动依旧给国家公园生态环境带来极大压力, 同时, 由于国家公园植被覆盖度高, 其地表生态环境对气候及环境变化也较为敏感^[3]。因此, 定量评价国家公园生态环境质量变化和气候影响规律显得尤为迫切, 对逐步恢复国家公园的原始自然生态环境及其高质量发展有着重要的现实意义。

迄今我国设立了10个国家公园试点区, 涉及12个省, 总面积达到22.3万km^2[4], 国家公园的生态环境质量相关问题受到学者广泛关注, 左婵等^[5]利用GLOPEM-CEVSA模型模拟植被潜在NPP, 分析其时空格局及人类活动对三江源国家公园生态环境质量的影响。刘晓娜等^[6]按照“因素识别-指标构建-单因子评估-综合评估”的思路, 评估色林错-普若岗日国家公园的生态环境脆弱性, 提出提升该地区生态环境质量的对策建议。贾岚等^[7]采用层次分析法构建武夷山国家公园生态环境质量评价指标体系, 系统分析其生态环境的变化动态并探究驱动机制。

作为现有的生态环境质量评价标准规范, 《生态环境状况评价技术规范(试行)》(HJ 192-2015)通过计算生态指数(EI)来评价区域生态环境质量^[8], 但该方法存在难以提取评价指标、不能很好反映生态环境质量空间变化等问题^[9]。近年来, 许多学者使用基于遥感的综合因素指数评价法来评价生态环境质量^[10—12], 其中结合了绿度、干度、热度、湿度的遥感生态指数(RSEI)^[13]在国家公园生态环境质量评价方面应用最为广泛^[14—15]。然而, 多数RSEI模型选用的绿度指标为归一化植被指数(NDVI), 虽然NDVI目前使用较广^[16—17], 但相比于NDVI, 增强型植被指数(EVI)改良了植被指数的算法及合成方法^[18—20], 弥补了NDVI易饱和的不足^{[19, 21—22]}。另外, 长时间序列、高精度的遥感数据集的获取及繁杂的数据预处理和RSEI综合指数计算工作通常费时费力^[23], 遥感计算云平台Google Earth Engine(GEE)可以便捷地利用高性能资源用以计算大型地理空间数据集, 在大空间范围长时序遥感应用研究中具有较大优势^[24—25]。

东北虎豹国家公园2021年正式成立。20世纪50年代, 由于森林过度采伐、农田开垦, 加之全球气候变暖的背景下, 致使生态环境加剧恶化导致东北虎豹及珍贵植物濒临灭绝。随着天然林保护工程实施、自然保护区建立, 东北虎豹栖息地生态环境逐步改善。目前, 对东北虎豹国家公园的研究多集中在保护地体系政策规划^[26—27]、野生动植物^[28—29]、生物多样性^[30—31]、景观生态^[32]等方面, 而园区生态环境质量问题鲜有人研究。因此, 本文基于GEE云平台, 利用MODIS影像数据, 选用EVI代替其绿度指标中的NDVI, 重新构建RSEI指数, 对2001—2020年植被生长季(4—10月)东北虎豹国家公园生态环境质量的时空变化进行动态评价, 并分析生态环境质量变化与气候因素中气温、降水及日照时数的响应关系, 旨在为东北虎豹国家公园的生态环境保护和可持续发展提供科学的参考依据。

1 研究区和数据来源 1.1 研究区概况

东北虎豹国家公园(以下简称“虎豹公园”)地处我国吉林、黑龙江两省交界的老爷岭南部区域, 地理坐标北纬42°31′06″—44°14′49″, 东经129°5′0″—130°18′48″, 总面积14926 km²(图 1)。虎豹公园所在区域是中俄朝三国交界地带, 园内分区实行差别化管控, 按照生态系统完整性、原真性和保护发展协调性三大原则将虎豹公园分为核心保护区和一般控制区, 核心保护区由“一大两小”共3个区域构成(以下分别简称①、②、③号保护区, ①号保护区即为东部核心保护区), 一般控制区又分为虎豹潜在栖息地和主要人口聚集区^[33]。

1.2 数据来源

用于构建RSEI所涉及的四个分量指标的MODIS产品数据来源于美国地质调查局(United States Geological Survey, USGS), 此数据集成于GEE云计算平台中, 时间范围为2001年1月1日—2020年12月31日。气象数据为2001—2020年中国1 km分辨率逐月平均气温数据集、中国1 km分辨率逐月降水量数据集以及2001—2015年中国1 km分辨率逐月日照数据集, 该数据来源于国家地球系统科学数据中心(http://www.geodata.cn/), 2016—2020年日照数据来自中国气象数据网(http://data.cma.cn/)的气象站点日数据, 利用ArcGIS 10.8软件对该数据进行空间插值，制作研究区1 km分辨率逐月日照数据。地表覆盖数据采用2020年GlobalLand30全球地表覆盖数据(http://www.globallandcover.com/), 数据共包括10个一级类型, 分别是：耕地、林地、草地、灌木地、湿地、水体、苔原、居民用地、裸地、冰川和永久积雪。数据详细信息见表 1。

2 研究方法 2.1 RSEI模型构建

RSEI模型包含绿度、湿度、热度、干度4个自然指标, 在此基础上, 选用EVI代替绿度指标中的NDVI, 弥补了NDVI易饱和的问题, 能够稳定地反应虎豹公园植被的生长状况, 改进后的RSEI模型方程为：

(1)

式中, EVI代表研究区绿度指标, 研究中采用MOD13A1中的EVI指数；WET代表研究区湿度指标, 它映了水体和土壤、植被的湿度, 计算方法见公式(2)^[34—35]；LST地表温度代表研究区热度指标, 研究中采用MOD11A2中的LST指数；NDBSI代表研究区干度指标, 由裸土指数(soil index, SI)^[36]和建筑指数(index-based built-up index, IBI)^[37]组合而成, 能够反映一个地区地表干燥程度, 计算方法见公式(3)—(5)^[38]；f是主成分分析方法(Principal Components Analysis, PCA), 是一种多变量统计方法, 它根据所用数据本身的性质、上述四种指标对各主分量的贡献度来自动客观地确定各指标权重, 从而避免人为主观因素在权重设定过程中的结果偏差^[34]。

(2)

(3)

(4)

(5)

式中, ρ₁—ρ₇为MOD09A1影像数据集中的波段1至7, 分别对应红波段、近红外1波段、蓝波段、绿波段、近红外2波段、短波红外1波段、短波红外2波段。

最后, 对上述4个自然指标去除量纲并将他们的数值映射到[0, 1]区间, 通过主成分分析法构建RSEI, 为便于比较, 同样对RSEI进行归一化处理至[0, 1]区间。生态环境质量分为五个级别, 对应值分别为：差(0—0.2)、较差(0.2—0.4)、中等(0.4—0.6)、良好(0.6—0.8)、优(0.8—1)^[39]。以上RSEI计算工作均在GEE平台完成。

2.2 变异系数

变异系数(Coefficient of Variation, CV)分析方法可以有效地揭示数据波动, 本文采用2001—2020年虎豹公园RSEI进行波动分析, 公式如下：

(6)

式中, σ_RSEI为RSEI数据的标准差, RSEI_mean为RSEI的平均值, CV可用于衡量虎豹公园整体生态环境的变化。CV值越高表明RSEI时间序列的波动性越大, 反之则越稳定。为了更直观地反映虎豹公园RSEI的变化情况, 将CV值细分为5个级别：低波动(CV≤0.1)、相对低波动(0.1＜CV≤0.2)、中等波动(0.2＜CV≤0.3)、相对高波动(0.3＜CV≤0.4)、高波动(CV＞0.4)^[40]。

2.3 趋势分析法

为分析虎豹公园2001—2020年RSEI的时空变化特征, 本文采用一元线性回归法逐像元地分析RSEI年际变化趋势, 计算公式如下^[41]：

(7)

式中, θ_Slope为一元线性回归方程的斜率, i为年份, RSEI_i为第i年的RSEI, n为年份总数。

2.4 Hurst指数

本文采用Hurst指数分析法, 分析虎豹公园2001—2020年RSEI的持续特征, 计算公式、原理参考谢平、江田汉、严恩萍等^[42—44]。Hurst指数取值包括三种形式: 当0.5＜H＜1时, 表明RSEI时间序列为一个持续型序列, 未来变化趋势与过去变化趋势一致, 越接近1, 其持续性越强；若H=0.5, 则说明RSEI时间序列为随机序列, 未来变化趋势未知；若0＜H＜0.5, 表明RSEI时间序列具有反持续性, H值越接近于0, 其反持续性越强。

2.5 相关性分析

利用Pearson相关系数法, 基于像元尺度, 对各像元植被生长季RSEI与气温、降水、日照因子进行相关性分析^[45]。

3 结果与分析 3.1 RSEI各指标主成分分析

在近20年研究时段内, 由于数据量较大, 本文选取2001年、2005年、2010年、2015年、2020年的主成分分析结果。从表 2可知, 与其他成分相比, 第一主成分(PC1)特征值贡献率在4项指标中均占到70%以上, PC1集合了各指标的绝大多数特征, 且PC1各指标载荷值相对均匀, 证明RSEI是一个可行的遥感生态指标, 基于PC1提取的指标信息来表征RSEI是合理的。

3.2 虎豹公园生态环境质量时空分布特征

根据2001—2020年虎豹公园整体RSEI多年均值结果显示, 各管控区生态环境质量空间分布表现出明显的分异性(图 2)。其中, 核心保护区方面, ①号保护区总体呈中间高南北低的分布态势, ②保护区北高南低, ③号保护区总体偏低。虎豹潜在栖息地高值区主要集中在中部偏南地区, 低值区集中在西部。主要人口聚集区生态环境质量相对最低。分析2001、2005、2010、2015、2020年的生态环境质量可以看出(图 3、表 3), 2005年前核心保护区及主要人口聚集区生态环境质量总体好于虎豹潜在栖息地。生态质量中等以下区域主要分布在虎豹潜在栖息地中西部。2010—2015年期间, 生态质量中等以下区域主要分布在①号核心保护区、主要人口聚集区大部及虎豹潜在栖息地东部。2020年园区生态环境质量等级总体达到良好。分析结果显示, 自园区西北向南又转向东部至①号核心区腹地可见一个较为固定的“弧形生态走廊”, 近20年间该走廊地带生态环境质量随时间发生波动变化。

不同土地覆被类型的多年RSEI均值也存在显著差异。虎豹公园的地表覆盖分为耕地、林地、草地、湿地、水体、居民用地和裸地七种类型, 前三种类型土地面积占全域的98.6%。其中耕地、林地和草地覆被面积依次顺序为：林地＞耕地＞草地, 而各地类RSEI值的排序依次为：耕地＞草地＞林地(表 4), 耕地具有相对较高的湿度和绿度值, 且受人类活动影响较大, 故RSEI值相对最高。

3.3 虎豹公园生态环境质量动态变化 3.3.1 多年生态环境质量变化分析

近20年虎豹公园RSEI变化范围为0.316—0.657, 园区历年RSEI在波动中总体呈微弱下降趋势, 其变化大致分为“上升-下降-上升”三阶段(图 4)。其中, 2001—2006年RSEI年均值波动上升, 变化趋势线斜率为0.03/a。2006—2015年RSEI年均值整体呈波动下降趋势, 由2006年的0.594降至2015年的0.359, 变化趋势线斜率为-0.013/a。2015—2020年园区RSEI年均值快速上升, 由2015年的0.359升至2020年的0.604, 变化趋势线斜率为0.034/a。随着2015年地方政府陆续停止虎豹公园内天然林商业性采伐, 园内重点国有林区进入全面保护的新阶段, 这与2015年以后生态环境质量的向好发展具有一定的关联性。

通过多时相的RSEI差值以5年为跨度分析虎豹公园生态环境质量变化情况(图 5), 结果表明：相较于5年前, 2005年生态环境质量变化不大, 2010年发生较大退化, 2015年后生态环境质量明显提高, 总体上2020年较2001年有11.6%的区域生态环境质量得到提升。其中, 2005年生态环境质量等级相较于2001年以不变为主, 降低的区域占全域面积的17.6%；相比于2005年, 2010年生态环境质量下降明显, 下降区域达到56.61%, 主要分布在①号保护区北部和南部、潜在栖息地东部和主要人口聚集区。各地类生态环境质量下降的比例依次为：耕地88.5%、草地79.6%、林地53.9%；2015年与2010相比总体环境质量继续下降, 有33.8%的区域生态环境质量进一步退化, 分布较为分散, 各地类生态退化区域各自占比分别为：草地40.3%、林地34.2%、耕地22.7%；2020年较2015年生态环境质量有较大程度提升, 除“弧形”走廊带呈降低情况外, 其他大部地区高达65.7%的区域生态环境质量有所提高。生态环境质量提升的地类占各自面积的比例顺序为：草地85.1%、耕地73.3%、林地64.6%。

3.3.2 生态环境质量稳定性与持续性分析

2001—2020年虎豹公园生态环境质量表现出较高的波动性, 稳定性较差。其中, 高波动区域面积占总面积的42.7%, 主要分布在①号保护区北部和南部及主要人口聚集区。相对低波动和中等波动区域分别占总面积的26.9%和20.2%, 主要分布在①号保护区中部、②号保护区以南及虎豹潜在栖息地中西部。相对高波动以及低波动区域仅占全域面积的9.9%和0.3%(图 6)。

虎豹公园大部地区未来生态环境变化态势具有较强持续性, 且生态环境发展趋势总体向好。近20年RSEI呈减小(Slope＜0)趋势的区域占园区总面积的58.63%, 生态环境恶化的地区主要分布在①号保护区南部、虎豹潜在栖息地东部和主要人口聚集区大部(图 6)。Hurst指数范围在0.25—0.71之间, 其中高于0.5的区域面积占64.03%(图 7)。将虎豹公园过去20年的RSEI变化趋势与Hurst指数值相叠加, 把未来RSEI变化趋势分为持续上升、持续下降、上升转下降、下降转上升四种类型(图 7)。这四种类型的面积占比分别为28.80%、33.07%、6.07%和32.06%。生态环境持续改善及下降趋势未来将转为上升的地区主要分布在虎豹潜在栖息地中部、①号、②号保护区大部。生态环境持续恶化的地区主要分布在虎豹潜在栖息地西部和东部、①号保护区南部局部和主要人口聚集区大部。

3.4 虎豹公园生态环境质量对气候的响应 3.4.1 气候因子变化趋势

2001—2020年虎豹公园4—10月平均气温、降雨量及日照时数均呈年际波动上升趋势, 且空间分布表现出明显的空间集聚特征(图 8、图 9)。其中, 平均气温在12.24—13.31℃之间, 多年均值为12.79℃, 趋势变化率为0.0025℃/a。降雨量在423.2—773.0 mm, 多年均值为603.8 mm, 最大和最小降雨量分别出现在2016年和2011年, 趋势变化率为6.9369 mm/a。日照时数在1340.27—1621.22 h之间, 多年均值为1451.88 h, 趋势变化率为2.0817 h/a。气温高值主要分布在①号保护区的东南部和南部、虎豹潜在栖息地的东部及主要人口聚集区, 低值区主要分布在①号保护区中部。年均降雨量从虎豹公园西部到东部呈递增趋势, 高值区主要分布在①号核心保护区大部。日照时数从虎豹公园西部到东南部呈递增趋势, 高值区主要分布在①号核心保护区东南部。

3.4.2 生态环境质量与气候因子的相关性

逐像元计算2001—2020年植被生长季RSEI与同期平均气温、降水量、日照时数之间的相关性(图 10), 结果表明：降水量对虎豹公园生态环境的影响程度大于气温与日照时数的影响程度, 且降水量以正相关为主。其中, 年均RSEI与年均气温及日照时数之间的相关性不大, 气温总体呈不显著的负相关, 而日照时数表现为不显著的正相关。但年均RSEI与年平均降水之间的相关关系较大, 相关系数变化范围介于-0.79—0.87之间, 呈正相关关系的区域面积占总面积的88.1%, 其中显著性正相关(P＜0.05)区域占全域的67.3%。主要集中分布在①号保护区中部及虎豹潜在栖息地的中部及主要人口聚集区大部。

4 讨论与结论 4.1 讨论

本文构建RSEI指数, 在做绿度指标的选取时, 利用2020年4—7月的多幅同期MOD13A1的NDVI及EVI影像, 抽取各月相同经纬度像元数据作为样本点, 用于定量统计不同时期样本点数据NDVI及EVI的分布情况, 结果发现：4—5月NDVI及EVI的分布趋势大致相同, 6月份NDVI及EVI相比5月总体呈增长趋势, 但NDVI在高于0.62以后分布曲线与5月部分重合, 而EVI仍在明显增长, 7月低值被覆盖样本NDVI和EVI分布相似, 而在高植被覆盖区样本NDVI曲线大多维持在6月份水平, 而EVI分布仍在增长, 这点在林地中尤为明显。这证明相较于EVI, NDVI在虎豹公园植被茂盛期的长势监测上出现了饱和现象, 故本文选用EVI代替原有绿度指标中的NDVI, 可以更真实的反映园区不同时期植被生长情况。同时, 对于主成分分析中的PC1, 绿度和湿度指标均高于零, 对RSEI有正向影响, 热度和干度均小于零, 为负向影响, 这个结论与前人多数研究结果相同^{[16, 38, 46]}。

近20年生态环境质量稳定性较差, 一是受人类活动和相关政策影响较大, 另一个主要原因是气候条件。2002年生态环境质量为近20年来最差, 当年生长季平均气温为近20年最低。与2005年相比2010年发生较大退化, 2015年总体环境仍继续下降。分析2005、2010、2015年虎豹公园植被需水关键期(5—8月)累积降水量分别为546.9 mm、427.9 mm、442.5 mm, 2010年较2005年同比下降29%, 而2015年累积降水量虽略高于2010年, 但2015年6月、7月平均降雨量较2010年同期分别下降45%和37%, 虽然后期降水量增加, 但难以弥补前期土壤水分亏缺带来的影响, 这可能是导致2010年虎豹公园生态环境质量较2005年退化明显且2015年继续下降的自然因素。2009年6月平均降雨量高达145.6 mm为历史同期最高, 相比历年同期平均降雨量增加65.9 mm, 也是2009年生态环境质量较高的主要原因之一。

根据《东北虎豹国家公园总体规划(2017—2025年)》, 2014—2016年黑龙江、吉林、长白山森工及汪清县地方国有林业相继全面停止天然林商业性采伐, 加之2017年东北虎豹国家公园体制试点的建立^[33], 标志着园区进入全面保护的新阶段, 这表明保护政策对2015年以后园区生态环境逐步恢复起到较为明显的促进作用。虎豹公园是亚洲东北部温带针阔混交林生态系统的中心地带, 混交林是野生东北虎豹偏好利用的生境类型^[47], 园区内植物资源丰富, 主要有东北红豆杉、红松、钻天柳等上百种植物物种。植物作为生产者是生态链中最为重要的一环, 生态环境的改善其能够为虎豹喜食的植食性动物提供更为丰富的食物和更佳的虎豹隐蔽环境等, 从而为各类野生动物长期定居、栖息和繁衍提供了保障^[48]。虽然虎豹公园大部地区环境质量呈向好趋势, 但局部地区仍存在生态环境质量持续下降的情况, 这些地区可能受人类活动的影响较大。随着近代以来人类生产生活等开发建设, 导致森林破碎化严重, 园区内落后的生产生活方式至今仍未发生改变^[32], 人类的过度开发、开采活动造成不同程度环境破坏难以短期内恢复。

本文利用RSEI分析生态环境质量, 由于生态系统复杂多样且具有许多特征, 该方法虽然总体上能反映生态环境质量状况, 但仍有一定片面性, 今后可以进一步扩大指标的数量更加精准反映区域生态系统质量。另外, 生态环境质量的驱动因素很多, 除自然因素的影响, 人类活动对生态环境仍起较大作用, 后续研究可以结合如夜间灯光、人口迁徙等表征人文经济的数据, 从而更加全面的分析区域生态环境影响因素。

4.2 结论

(1) 2001—2020年虎豹公园生态环境质量呈“上升-下降-上升”的微弱下降趋势, 2015年后生态环境质量逐步提升, 2020年达到良好水平。生态环境质量较好区域主要分布在①号保护区中部及虎豹潜在栖息地南部。不同土地覆盖类型的历年RSEI均值依次为耕地＞草地＞林地。

(2) 相较于2001年, 2005年生态环境质量变化不大, 2010年发生较大退化, 退化的区域达到56.61%, 主要分布在①号保护区北部和南部、潜在栖息地东部和主要人口聚集区。2015年总体环境仍继续下降, 有33.8%的区域生态环境质量进一步退化。2015年后生态环境质量明显提高, 总体上2020年较2001年有11.6%的区域生态环境质量得到提升。

(3) 虎豹公园近20年生态环境比较脆弱, 环境质量总体表现出较高的波动性, ①号保护区北部、南部及主要人口聚集区RSEI的波动性要明显高于虎豹潜在栖息地。未来环境质量变化趋势持续向好、持续变差、向好转变差、变差转向好的面积分别占整个园区面积的28.80%、33.07%、6.07%和32.06%, 表明生态环境趋势总体向好, 但局部地区仍存在持续恶化态势。

(4) 虎豹公园历年气温、降水和日照在波动中呈上升趋势, 生态环境质量与气温及日照的相关性不高。而降水与生态环境质量相关性较强, 主要呈正相关关系, 其中显著性正相关区域主要分布在①号保护区中部及虎豹潜在栖息地的中部及主要人口聚集区大部。