梵净山长时序植被覆盖变化趋势与对比分析

余 婷,杨 萍,邓 兴,刘绥华,周 洋

(1. 贵州师范大学地理与环境科学学院,贵州 贵阳 550000; 2. 贵州省测绘产品质量监督检验站,贵州 贵阳 550000; 3. 贵阳市测绘院,贵州 贵阳 550000; 4. 青岛四洲重工设备有限公司,山东 青岛 266000)

植被是生态系统中不可替代的关键组成部分,对于建立和维护良好的生态环境起着重要作用[1-2]。植被的变化不仅反映了地球表面能量交换、碳循环和人类活动,还连接了土壤、大气和水之间的复杂关系[3-5]。植被覆盖度(fractional vegetation cover,FVC)是评估植被覆盖程度和植被生长状况的常用指标,它是指地表植被的投影面积与地表总面积的比率,通常与植被类型、生长状况和高度等因素有关[6-7]。

近年来,随着遥感技术的蓬勃发展,遥感数据以覆盖面积大、精度高和节省时间等优势,成为生态环境监测领域中不可或缺的资源[8-9]。但传统的计算机平台难以应对日益复杂的数据存储及计算需求。这时,GEE(Google Earth Engine)的发展引起了国内外的广泛关注,并为海量数据的快速处理、分析提供了前所未有的发展机遇[10]。已有学者基于GEE对不同地区的植被覆盖进行了相关研究,如文献[11]研究了民勤绿洲NDVI的时空变化特征,发现绿洲有向荒漠扩张的趋势。文献[12]采用Landsat数据对黄土高原地区进行植被覆盖度时空变化的影响因素研究,得出人类活动是植被动态变化的主要驱动因素。

梵净山自然保护区是世界上同纬度保护最完好的原始森林。1987年,被列为国际“人与生物圈”保护区网络的一员;2018年,又被列为世界自然遗产名录中的第53项和我国第13项世界自然遗产[13]。该保护区具有丰富的植被资源,但由于气候、人类活动等多种因素的影响,植被覆盖在过去的几十年中发生了明显的变化。文献[14]利用遥感影像和土地利用数据对梵净山植被覆盖度时空变化特征进行了研究。结果表明,相比于提名地,缓冲区的植被覆盖度变化更为敏感。然而,该研究并没有对这种变化的趋势进行分析。常见的趋势分析方法包括回归分析和Sen+Mann-Kendall法分析,但大部分研究通常只采用其中一种方法进行研究[15-17],很少有研究对这两种趋势分析方法进行对比。

基于此,本文利用GEE云平台和Landsat 地表反射率(surface reflectance,SR)数据集作为试验数据,采用像元二分模型法、逐像元一元线性回归分析法和Sen+Mann-Kendall分析法对梵净山自然保护区1990—2020年共30年植被覆盖度的时空演变和长时序的变化趋势进行研究,并将两种趋势分析方法进行对比。通过分析植被覆盖度的变化趋势,可以揭示保护区环境的变化特征,为保护区的生态管理提供科学依据,制定合理的保护措施,促进生态系统的恢复和保护。

1 研究区概况

梵净山世界自然遗产地、国家级自然保护区位于贵州省东北部的印江县、松桃县和江口县三县交界之处,突立于云贵高原东部向湘西低山丘陵过渡的斜坡地带,海拔高程在414～2572 m之间,高差大于2000 m,最高峰凤凰山海拔为2572 m。遗产地面积为402.75 km2,缓冲区面积为372.369 km2。根据《中国植被》分类单元,将梵净山的自然植被划分为针叶林、阔叶林、竹林、灌丛(灌草丛)、沼泽等5个植被型组。

2 数据和方法

2.1 数据源

本文以Landsat SR数据集作为试验数据,包括Landsat 5 TM Surface Reflectance Tier1(1990—2011年)数据集和Landsat 8 OLI Surface Reflectance Tier1 (2013—2020年)数据集,共计180景遥感影像,影像分辨率为30 m。由于2012年图像质量不佳,故本文使用30 年数据。数据均来源于美国地质调查局(USGS),并通过GEE平台获取,该数据集已经过大气校正、几何精度校正和去云处理等,可在平台直接进行调用,确保了数据的科学性。为消除大气和传感器角度的影响,减少运算误差,采用了最大值复合法合成了研究区域1990—2020年的年度NDVI数据集。

2.2 研究方法

2.2.1 像元二分模型

计算植被覆盖度的主要遥感手段包括经验模型法、植被指数法和像元分解模型法。其中,像元分解模型法中的像元二分模型具有反演精度高、数据参数易获取、操作简单可靠等优点[18]。该方法的计算公式为

(1)

式中,NDVIsoil和NDVIveg分别表示无植被覆盖像元和纯植被覆盖像元的NDVI值,其理论值分别接近0和1。本文中,选用5%置信度确定NDVI取值区间。

2.2.2 一元线性回归趋势分析

一元线性回归趋势法通过对变量和时间之间的线性关系进行回归分析,来推断变量随时间变化的趋势[19]。公式为

(2)

式中,slope为植被覆盖度变化的斜率;n为时间序列长度(文中n=30);FVCi为第i年的FVC。当斜率为正时,表明研究区FVC在该时间序列下呈现增长趋势,反之则呈现减少趋势。

2.2.3 Sen+Mann-Kendall趋势分析

Theil-Sen Median方法又被称为Sen斜率估计法,是一种用于确定趋势的高效且稳健的非参数统计方法。Sen斜率估计用于计算趋势值时,通常与Mann-Kendall非参数检验结合使用[20]。Sen斜率计算公式为

(3)

式中,j、i为时间序数;Xj、Xi分别为第j和第i时间序列数据。当斜率β>0时表示FVC呈现上升趋势;β<0时表示FVC呈现下降趋势。Mann-Kendall统计检验公式如下。定义标准化检验统计量Z为

(4)

其中

(5)

(6)

(7)

式中,n为序列中数据的个数;Xj、Xi分别为第j和第i时间序列数据。

3 结果分析

3.1 FVC时间变化特征

根据梵净山植被覆盖度的实际情况,本文将梵净山的FVC平均值划分为5个等级:低等级FVC(<0.3)、中低等级FVC([0.3,0.45))、中等级FVC([0.45,0.6))、中高等级FVC([0.6,0.75])、高等级FVC(>0.75)。

如图1(a)所示,近30年来,梵净山植被覆盖度呈先下降后上升的“U”形增长趋势,增长率为3.46%,年均FVC为0.58,总体处于中等级水平,年际变化呈明显增加态势。从图1(b)可知,在1990年,高等级FVC面积占比最大,达24.24%;而低等级和中低等级FVC则在1996年占比最大,达13.40%;但在2002年之后,该等级比例减少到2.49%,这一阶段植被覆盖得到了显著改善,主要原因是2002年退耕还林还草工程在全国范围内全面启动,从而成为植被覆盖变化的一个重要转折点。然而,在2008年研究区域内的低等级FVC面积占比达3.06%,与1990年相比,增加了2.62%。这可能是人类活动导致的结果,在2018年7月2日的世界遗产大会上,梵净山自然保护区被列入世界自然遗产名录。随后,旅游业在这段时间内快速发展,而旅游业的发展需要大量的土地资源,这就导致原有的草地、森林等土地类型被转变为旅游景区或度假村用地,从而植被覆盖度急剧下降。

图1 1990—2020年梵净山平均FVC年际变化及面积占比

3.2 FVC空间演变特征

1990—2020年梵净山多年平均FVC空间分布如图2所示。总体而言,FVC呈现出中间高、四周低的空间分布特征。中等级及中高等级的FVC面积占比较多,FVC年平均值大于0.6的区域占整个区域的93.55%;低等级FVC面积占比最少,为0.66%;高等级FVC主要分布在梵净山自然保护区的山顶地区,具体为新金顶(海拔2342 m)、老金顶(2494 m)的西侧;低等级及中低等级FVC主要出现在梵净山自然保护区东侧及东南侧的沿山公路带上,以及北侧和西侧的村庄,均属于缓冲区内,这与文献[14]的研究结论一致。这可能是由于梵净山地形高差较大,山顶区域的环境条件相对较好,有利于植被的生长,而山脚下和缓冲区的土地利用、人类活动较为密集,对植被的生长有不利影响。

图2 1990—2020年梵净山多年平均FVC空间分布

1990—2020年梵净山各等级FVC空间分布变化如图3所示,1990—1996年,高、中高等级FVC面积占比减少了16.92%,主要分布在梵净山遗产地范围的中部海拔较高的地方;低、中低等级FVC开始大范围分布,从自然保护区东侧的缓冲区向核心区蔓延,1996年中低等级FVC面积占比达到11.32%。1996—2002年,中低等级以下的 FVC主要分布在梵净山自然保护区东南方的低平地区,该地区为江口县北部,从Google Earth影像上看,该区域是一些居民点及修建的盘山公路。这个时期在研究区北边出现了呈条带状分布的低等级FVC,与实际植被分布情况无关,这是由影像质量所导致的。2002—2008年,低等级FVC在梵净山核心区也有分布,且主要分布在几大景点高山之间;低等级FVC多分布于西侧的缓冲区域内;东南侧的低等级FVC面积占比相较于前十几年有所增加。2014—2020年,植被覆盖情况整体开始改善,尤其是在梵净山自然保护区的缓冲区域,低等级和中低等级FVC面积占比急剧下降,进而转变为中高等级FVC。

3.3 FVC变化趋势分析

3.3.1 一元线性回归法趋势分析

植被覆盖的改善面积远大于退化面积,分别为96.18%和3.82%(如图4、表1所示)。植被总体上呈改善趋势,在空间上,显著改善的区域主要位于梵净山北部、西部。显著退化和轻微退化分别占0.46%和3.36%,主要分布在自然保护区的缓冲区东南部、西北部及核心区中部凤凰山附近。这表明,旅游活动对植被退化有着直接的影响。

图4 一元线性回归法结果

3.3.2 Sen+Mann-Kendall趋势分析

Sen+Mann-Kendall趋势分析结果与一元线性回归趋势分析法结果基本一致。1990—2020年,植被改善区域(94.80%)明显大于退化区域(5.20%),见表2。如图5所示,明显改善区域主要分布在梵净山自然保护区缓冲区的北部、西部及东北部,这与梵净山近年来实施的生态移民搬迁保护措施密切相关。退化区域主要分布在核心区的西部、景点山峰周围以及缓冲区东南部和西北部。其原因在于这些区域分布着环山公路、零星居民点及旅游设施等。

表2 Sen+Mann-Kendall法趋势分析结果

图5 Sen+Mann-Kendall法结果

3.3.3 一元线性回归分析法与Sen+Mann-Kendall法的差异性

研究发现,两种趋势分析方法得出的梵净山植被覆盖度变化趋势基本相似,均表明在30年长时间序列中植被改善面积远大于植被退化面积。与一元线性回归法相比,使用Sen+Mann-Kendall法监测到的植被呈退化趋势的面积有所增加。此外Sen+Mann-Kendall法在准确模拟轻微退化和明显改善的区域方面优于一元线性回归法,对明显变化结果的敏感性提高了2.18%。因此,Sen+Mann-Kendall法的分析结果更精确,能更好地监测植被覆盖的变化,也为保护生态环境提供了重要的科学依据。

4 结 论

本文利用GEE云平台,基于梵净山自然保护区的Landsat SR数据集,研究了1990—2020年梵净山自然保护区FVC的时空分布和变化趋势。主要结论如下:

(1)近30年来,梵净山自然保护区的植被覆盖度呈先下降后上升的“U”形增长趋势,增长率为3.46%,年均FVC为0.58,年际变化呈明显增加态势。中低、低等级植被覆盖度面积占比随着时间推移显著下降,中等级及以上等级的植被覆盖度则呈现出显著上升的趋势,植被覆盖情况得到了明显改善。

(2)梵净山自然保护区植被覆盖度的空间分布呈现出中间高、四周低的分布特点,即遗产地内大多呈现中等级及以上等级FVC,较低等级和低等级FVC则出现在缓冲区域内。从变化趋势上来看,1990—2020年梵净山自然保护区植被覆盖情况明显好转,植被改善面积远大于植被退化面积。明显改善区域主要分布在梵净山自然保护区缓冲区的北部、西部及东北部。退化区域主要分布在核心区的西部、景点山峰周围及缓冲区的东南部和西北部。

(3)两种趋势分析方法得出的梵净山自然保护区植被覆盖度变化趋势结果基本相似。相比之下,Sen+Mann-Kendall趋势分析法在准确模拟轻微退化和明显改善的区域方面优于一元线性回归法。此外,Sen+Mann-Kendall 趋势分析法对明显变化结果的敏感性提高了2.18%。

本文仅使用了Landsat SR数据集,并未使用多源数据集对结果进行验证。由于不同影像数据集分辨率和质量不同,对于植被覆盖度变化的结果可能存在差异。因此,有必要进一步分析不同数据集之间的差异性,以评估它们的一致性和不确定性,从而获得更加可靠和稳健的结果。