基于时间序列分割算法的雅鲁藏布江流域NDVI(1985—2018)变化模式研究

王 塞,王思诗,樊风雷,2,*

1 华南师范大学地理科学学院, 广州 510631 2 西藏大学高原地表环境遥感监测联合实验室, 拉萨 850000 3 天津大学环境科学与工程技术学院, 天津 300350

植被是联结水分、大气和土壤的重要纽带,也是评价生物多样性和生态服务的重要指标,在全球碳平衡和气候变化中扮演重要角色[1-3]。雅鲁藏布江流域生态环境脆弱,对气候变化和人类活动的响应显著[4]。近30年来,流域所属的青藏高原年平均气温每10年上升0.7℃,增温趋势总体呈现由南向北逐渐降低,降水量受复杂地形影响在不同区域差异显著,流域气候暖干化趋势明显。由此,干旱程度加剧导致的植被稀疏化对流域脆弱的生态系统产生不利影响。此外,土地过度开垦、过度放牧和城市化进程等人类活动也加剧了流域植被退化[5-7]。因此,雅鲁藏布江流域的植被变化及其影响因素倍受学界关注。目前研究的方法与途径主要是,以遥感影像为数据源,通过对不同类型的植被覆盖指标与气温、降水等气候因子进行耦合测度,进而判别该流域内生态环境状况[6,8-10]。开展高原高寒植被物候期变化研究,讨论气候变化对植被生长的影响与作用规律[11]。然而,对雅鲁藏布江流域植被变化模式研究鲜有相关报道。

随着遥感技术的发展,中分辨率遥感数据已经成为在国家和大陆尺度上进行环境监测的重要数据源[12]。基于中分辨率长时间序列遥感卫星影像是对大范围区域进行长期监测,分析其土地利用变化趋势,并对植被变化模式进行判别,已成为评价区域植被变化的重要手段[13-14]。归一化植被指数(Normalized Difference Vegetation Index,NDVI)在植被变化研究中得到广泛应用。NDVI反映植被冠层在近红外与红光波段的吸收和散射作用,与植被的光合作用与生理生化反应密切相关,因此既能监测植被对全球气候变化的响应,又能识别城市扩张对植被的影响[15]。中分辨率长时间序列NDVI数据集成以上优势,已逐渐成为研究植被变化重要的数据源。NDVI时间序列变化监测算法按原理可分为阈值、差异、分割、轨迹分类、统计和回归六类[16]。Landsat数据源的时间序列分割算法(LandTrendr)在森林植被覆盖变化检测与生长模式判别等研究方面取得了丰硕的成果[13,17-18]。Google Earth Engine(GEE)是处理和存储地理大数据的云计算平台,已在作物生长监测、土地覆盖变化和灾害预警等领域成功应用[19]。GEE平台实现的LandTrendr算法(https://emapr.github.io/LT-GEE)已经成为大区域长时间序列变化监测的有利工具[20]。

综上,本研究以GEE平台共享的Landsat TM/ETM+/OLI影像为数据源,利用LandTrendr算法提取雅鲁藏布江流域NDVI时空变化特征。在系统分析1985—2018年雅鲁藏布江流域NDVI变化特征基础上,按NDVI不同的变化强度和持续时间进行分类,进而对流域内NDVI变化模式进行判别,并进一步探查不同变化模式产生的可能原因,以期为雅鲁藏布江流域生态环境保护提供决策支持。

1 资料与方法

1.1 研究区概况

雅鲁藏布江流域位于北纬27°50′—31°15′,东经81°57′—97°5′,发源于喜马拉雅山北麓的杰玛央宗冰川,流经里孜、派乡和巴昔卡,最后流入印度,流域总面积25.8万km2(图1)。流域内自然环境复杂,海拔范围在149—7172m之间,多年平均气温4.7—8.3℃,年降水量251.7—580.0mm[6]。复杂的自然环境造就了多样的植被覆盖类型,上游段植被覆盖以高山草甸和高山草原为主,中游段主要为灌木草原,下游段以乔木和次生植被为主要类型[21]。雅鲁藏布江流域因其独特的自然环境,造就了该区域多样的植被覆盖类型,因其脆弱和敏感的自然环境而备受研究者关注。

图1 研究区示意Fig.1 Map of study area

1.2 数据与处理

本研究中使用的NDVI时间序列数据由Google Earth Engine地理云计算平台提供的1985—2018年Landsat TM/ETM+/OLI影像生成时间序列堆栈计算生成,其余数据均来源于中国科学院资源环境科学数据中心[22]。其中Landsat 5 TM影像共1737景,Landsat 7 ETM+影像共1877景,Landsat 8 OLI影像共1077景,并对数据堆栈中数据质量较差区域(云,云阴影和雪)进行掩膜。由于Landsat 7 ETM+相比Landsat 8 OLI数据具有较低的辐射分辨率,本研究利用反射率归一化最小化传感器之间的差异[23]。为获取以年为周期的NDVI时间序列数据,将研究区内同一年份所有NDVI数据进行年最大值合成,以削弱大气气溶胶和太阳高度角等因素的影响。

1.3 研究方法

1.3.1NDVI时间序列分割算法

时间序列分割算法(LandTrendr)对以年为间隔的NDVI时间序列进行分割、拟合和平滑,获取单个像元NDVI值在整个研究时间段内的变化特征[17](图2)。该分割算法主要分为如下4个步骤：

图2 LandTrendr分割算法示意Fig.2 Conceptual diagram of LandTrendr segmentation algorithm

(1)数据预处理。去除NDVI时间序列轨迹异常值,设定阈值对轨迹中的突变值进行筛选；去除超过NDVI正常变化趋势的异常值,并对缺省值进行插值填补。为削弱异常值对变化趋势的影响,增加信噪比,本研究采用3×3滑动窗口均值进行异常值剔除。

(2)时间序列分割。识别轨迹变化剧烈位置作为潜在分割点,并对潜在分割点进行判断。对比潜在分割点NDVI值,剔除由云、云阴影和雪等因素造成的伪分割点,并通过角度阈值剔除小于阈值的分割点。

(3)分割点拟合。对比第一和第二分割点间连线和拟合直线与原始NDVI轨迹间的均方误差(Mean Square Error,MSE),选取误差最小的拟合方法作为初始连接线,并在保证曲线连续的情况下按照MSE最小原则对NDVI时间序列进行拟合。

(4)拟合模型简化与遴选。通过迭代最大分割时间序列,利用定义的恢复率和置信度对模型进行筛选。考虑到植被覆盖恢复期应长于1年,剔除恢复率过高的分割段,使拟合模型具有更好的分割效果

本研究基于GEE平台实现的LandTrendr算法对NDVI时间序列进行分割与拟合[24]。分割算法需要设置适用于研究区的分割和拟合参数(表1),通过多次参数调整以适应该研究区的NDVI变化特征。

表1 时间序列分割模型参数Table 1 Time series segmentation model parameters

1.3.2NDVI变化概念模型

全球气候变化和人类活动对雅鲁藏布江流域脆弱的生态环境产生显著影响。温度和水分在整个生长季内对植被生长产生显著影响,水分控制作用是影响青藏高原植被生长的最主要因素[25],城市化过程加剧和过度放牧等人类活动促使该区域NDVI下降,一定程度的气温和降水量增加以及退耕还林还草等生态保护工程都对NDVI上升产生促进作用[26]。基于前人成果,本研究构建雅鲁藏布江流域NDVI变化概念模型(图3),即,全球暖湿化和生物质累积使年时间分辨率NDVI曲线随植被生长缓慢上升,人类活动加剧和自然环境恶化导致NDVI曲线下降,生态保护工程或气候暖湿化等自然环境改善又使NDVI曲线呈现上升趋势,生态保护工程效果减弱或者保护区域再次被破坏都将导致NDVI曲线的再次下降。影响NDVI曲线趋势的因素包括植被物候特征、人类活动和遥感影像采集的时间和质量,因此本研究将NDVI曲线波动最剧烈年份作为变化时间点,基于持续时间和变化强度对NDVI变化模式进行判别。

图3 NDVI变化概念模型Fig.3 Conceptual model for NDVI change

2 结果

2.1 干扰与恢复时空分布特征

2.1.1空间分布

利用时间序列分割算法对雅鲁藏布江流域1985—2018年NDVI时间序列数据进行变化特征提取,得出NDVI干扰与恢复空间分布情况(图4)。结果显示,自流域上游至下游发生干扰年份逐渐推迟的趋势明显,其空间异质性与地形特征存在明显的相关性。流域上游(贡嘎县以上)的广大区域NDVI干扰集中发生于2000年左右。流域中上游的贡嘎县、昂仁县和日喀则市NDVI干扰集中发生在1986—1990年间。中游当雄县等海拔较高区域NDVI干扰主要发生于在2000年之前,而拉萨河谷等海拔较低区域NDVI干扰主要发生在2010年后。流域下游工布江达县、米林县和林芝县海拔较低区域NDVI干扰年份发生在2000年左右,而海拔较高区域干扰多发生在2010年之后。自流域上游至下游恢复发生年份具有明显提前趋势,其与海拔、地形特征关联度较高。流域中上游的广大区域NDVI恢复年份发生在2010年之后。流域下游的米林县和林芝县NDVI恢复发生年份存在明显的海拔差异,高海拔区域恢复年份发生在2010年之后,低海拔区域恢复发生在2000年左右,但念青唐古拉山及唐古拉山以南的部分高海拔区域恢复发生在1990年之前,不同于以上规律。

图4 干扰与恢复空间格局Fig.4 Spatial pattern of disturbance and recovery

2.1.2时间分布

NDVI干扰与恢复面积随时间变化特征能表征植被覆盖变化活动的强弱,按NDVI时序变化趋势,流域内18.2万km2划分为干扰类别,24.3万km2划分为恢复类别(图5)。流域不同时间段植被覆盖变化强弱不同。1986—1990年,干扰和恢复面积皆呈先上升后下降趋势,干扰区域占流域总面积20.97%,恢复区域占总面积43.49%,恢复面积大于干扰二者相差22.52%,流域植被覆盖总体呈恢复态势。1991—2000年,干扰和恢复面积也呈现先上升后下降趋势,分别占总面积30.54%和27.67%,干扰面积略大于恢复,二者相差2.87%,研究区植被覆盖总体呈下降趋势。2001—2017年期间,干扰和恢复都呈现波动变化趋势,变化比例分别占流域总面积19.03%和23.02%,恢复面积略大于干扰,二者相差3.99%,流域植被覆盖度总体呈现恢复态势。综观1986—2017年间,研究区NDVI干扰和恢复植被变化面积总体呈先升后降趋势,恢复区域面积大于干扰区域面积,二者相差23.64%。

图5 1985—2018年干扰与恢复面积年际变化 Fig.5 Annual variation of disturbance and recovery area from 1985 to 2018

2.2 干扰和恢复持续时间与强度

2.2.1持续时间

1985—2018年间雅鲁藏布江流域NDVI干扰和恢复持续时间具有明显的时间序列非平稳性,干扰和恢复持续时间皆呈现先下降后上升的趋势,分别集中在1—10年和28—32年,其中二者下降持续期(1—10年)干扰略小于恢复,而上升持续期(28—32年)干扰明显大于恢复(图6)。NDVI干扰持续时间85%少于10年,干扰平均持续时间4.96年,其中49.08%的干扰持续时间少于等于1年。流域内NDVI干扰持续时间超过20年的区域主要集中在拉萨城区、墨竹工卡北部和仲巴县以南的部分区域(图7)。NDVI恢复持续时间50%少于6年,平均恢复时间12.55年。对比地形数据可见,NDVI恢复时间超过20年的区域主要分布在河谷等海拔相对降低的区域(图7)。

图6 干扰与恢复面积随持续时间和强度变化Fig.6 The duration and magnitude of disturbance and recovery area

2.2.2强度

NDVI干扰和恢复强度主要集中在0—0.75之间,恢复强度大于干扰强度(图6)。干扰活动剧烈区域主要集中在拉萨城区、林芝及其北部区域和里加(图7)。NDVI干扰强度95%都集中在0—0.42之间,平均下降0.14。NDVI干扰强度接近2时表明由植被转为水体。NDVI恢复强度95%集中在0—0.4之间,平均上升0.13,流域上游和中游区域NDVI上升强度较低,NDVI上升剧烈区域主要集中在里加和林芝以东区域(图7)。

图7 干扰与恢复持续时间和强度空间格局Fig.7 Spatial pattern of duration and magnitude of disturbance and recovery

2.3 干扰与恢复变化模式

1985—2018年NDVI时间序列为判别NDVI变化模式提供基础,将1987、2000和2016三年NDVI拟合数据进行假彩色合成,进而判别NDVI的不同变化模式。流域内NDVI变化主要分为上升和下降两种趋势,NDVI变化模式空间分布如图8所示,NDVI上升在整个流域内都有分布,而NDVI下降主要集中在河流谷地等海拔较低的区域,共划分为12种典型变化模式(图8),其中具有明显变化模式的面积占流域总面积87.23%。

流域内NDVI总体呈现上升趋势,其中上升模式E和下降模式E所占面积比例最高(图8)。下降模式E具有NDVI在1987—2000年间逐渐下降,2000—2016年NDVI快速下降的变化模式,主要集中在拉萨河和米林宽谷；下降模式F次之,表征NDVI在研究时间段内降低减速,主要分布在南木林和墨竹工卡等区域。上升模式E为流域内主要的NDVI增加变化模式,NDVI在1987—2000年快速上升,2000—2016年缓慢上升,主要集中在拉萨市北部和林芝县以东海拔较高的山区；先上升后下降的上升模式A主要分布在林芝和日喀则以东区域。

图8 雅鲁藏布江流域NDVI变化模式Fig.8 Different patterns of NDVI change in the Yarlung Zangbo River Basin

3 讨论与结论

3.1 讨论

本研究选择1985—2018年经过地形校正的Landsat TM/ETM+/OLI Level1数据,通过掩膜剔除云、云阴影和雪的区域,计算时间序列数据集的NDVI值,将一年内所有的NDVI数据进行最大值合成,相比较月或更短时间间隔的NDVI时间序列数据,该方法既能突出植被覆盖年际变化特征,又能削弱云、大气和太阳高度角等因素的干扰[13]。时间序列分割算法(LandTrendr)对NDVI时间序列进行分割和拟合,提取NDVI变化特征,广泛应用于植被覆盖变化及其变化原因分析[27]。相比于基于统计特征的时间序列变化检测方法,时间序列分割方法能提取子序列的时间变化特征,进而分析时间序列的变化模式[28]。

雅鲁藏布江流域植被对全球气候变化的响应显著,气候变暖和降水量增加导致植被覆盖增加,NDVI总体呈上升趋势,但受种间关系和降水时空分布的不均匀性,以及人类活动的影响,不同区域具有不同的变化趋势,这与已有学者研究结论相一致[29-30]。对比图4与图8, 念青唐古拉山及唐古拉山以南的部分高海拔区域恢复在1990年前的区域NDVI变化模式不显著,考虑到NDVI变化的复杂性,尚需进一步研究。

20世纪80和90年代,受全球气候变化的影响,流域上游和中游地区高原高寒植被物候期延长趋势明显,净初级生产力与得到显著提高[31],与研究结果上游和中游NDVI恢复持续时间较短且恢复持续强度较低相对应。上游和中游部分海拔较高区域NDVI干扰活动剧烈,可能与高海拔地区降水和长波辐射等因素对NDVI持续上升的限制因素有关[2]。下游森林覆盖区受采伐和保护双重作用,致使当地植被覆盖快速变化,如20世纪90年代前的森林大规模采伐、1998年后天然林资源保护工程实施[32],森林破坏面积下降,恢复面积缓慢增加,流域内恢复面积大于干扰破坏面积,NDVI干扰和恢复都经历了剧烈变化,但NDVI总体呈现增加趋[9],与本研究中恢复和干扰在2000年前变化剧烈,而2001年后变化活动逐渐减弱的结果相吻合。影响流域内植被恢复强度受不同植被覆盖类型存在差异,流域上至中下游存在明显的恢复强度差异。本文研究结果表明,干扰和恢复持续时间都出现了下降后上升的趋势,干扰和恢复持续时间主要集中在10年内的较短时期,分布集中在流域的中上游,这与近20年来青藏高原暖湿化,草地生态系统恢复期较短相对应,下游区域因森林经历大规模采伐,森林覆盖需较长的时间才能恢复。

NDVI增加模式以先加速上升后减速为主,主要位于雅鲁藏布江流域下游植被覆盖较好的区域,植被覆盖增加受气候和人类活动的共同影响,NDVI增加速率降低可能与降水的空间分布与树种组成结构变化有关[30]。随着生态保护工程的启动30多年[33],老年林所占的比例有所增加NDVI增加趋势逐渐减弱。NDVI下降速率逐渐增加主要集中在人类活动较多的宽谷地区,放牧和城市化过程加剧导致的植被破坏,NDVI下降模式在人口密集的海拔较低区域出现空间聚集[34],这与已有的研究结论相一致。考虑到影响雅鲁藏布江流域植被变化的因素众多,不同影响因素在不同的区域时空分布差异巨大,因此在分析NDVI变化模式特征及其时空分异时,还应分别对不同影响因素间相互作用模式进行耦合研究。

3.2 结论

本研究基于1985—2018年Landsat数据合成的年最大NDVI时间序列为数据源,通过Google Earth Engine平台实现的LandTrendr时间序列分割方法对NDVI时空变化特征进行提取,进而描述研究区内NDVI变化模式,主要结论如下：

(1)雅鲁藏布江流域植被时空变化显著,总体呈现上升趋势,仅在局部区域出现下降现象,自上游至下游NDVI变化强度逐渐增加。1985—1990年恢复面积大于干扰,NDVI变化剧烈；1990—2000年干扰和恢复变化面积差异不大；2000—2018恢复面积略大于干扰,NDVI变化活动逐渐减弱。

(2)1985—2018年间流域内NDVI干扰95%集中在0—0.42之间,平均下降0.14,平均干扰时间4.96年。NDVI恢复95%集中在0—0.4之间,平均上升0.13,平均恢复时间12.55年。NDVI干扰和恢复在强度上相近,但平均持续时间和变化面积差异显著,因此在近30年内NDVI总体呈现上升趋势,但上升趋势逐渐减弱。

(3)雅鲁藏布江流域NDVI变化模式以上升为主,流域内NDVI主要的干扰变化模式为先缓慢下降后快速下降；NDVI恢复模式为先快速上升后缓慢上升。