基于多源遥感数据的青藏高原植被变化特征评价

2022-09-15赵耀华彭小清金浩东彭思佳

冰川冻土 2022年4期

赵耀华，彭小清，2，金浩东，杜冉，陈聪，彭思佳

（1.兰州大学资源环境学院，甘肃兰州 730000；2.兰州大学祁连山冻土生态环境野外科学观测研究站，甘肃兰州 730000）

0 引言

全球平均表面温度持续上升，青藏高原气候变暖更为显著。近几十年来，青藏高原升温速率约为全球同期升温速率的2倍，达到每10年升高0.3～0.4℃［1］。快速的增温变暖给作为我国重要生态屏障的青藏高原带来巨大的生态安全挑战［2］。一方面，植被是陆地生态系统的必要成分，是生态系统中最为活跃的因素，具有连接土壤、大气、水分与其他生态环境因子的桥梁作用，其变化影响着地球重要圈层的物质循环与能量流动。另一方面，生态系统的变化可以被敏感的植被所反映［3］。科学地评估青藏高原植被变化有助于维护高原地区生态平衡以及恢复高原生态系统，以应对全球气候变化带来的负面影响。在生态文明加速发展的背景下，《生态文明体制改革总体方案》的部署也明确指出青藏高原生态屏障修复的重要性［4］。正确认识植被的变化情况将促进区域经济与生态文明协调发展。

全球环境的变化正在迅速影响着陆地植被动态，但植被应对快速变化的环境的机制尚不明确。因此，植被生态对于气候变化的反馈与响应受到专家学者持续性的关注。青藏高原属于高寒生态系统，主要覆盖着高寒草原与高寒草甸，相比其他生态系统更为脆弱，也由此成为相关研究的重点［5］。近些年来，对于青藏高原的生态系统的研究主要围绕着植被物候、植被绿度、林线位置、生物多样性、生态系统生产力与碳汇等主题［2，5-6］。其中植被绿度及物候变化是生态系统监测的重要内容，利用卫星遥感的衍生指标，比如归一化植被指数（NDVI），评估青藏高原植被长期状况，揭示植被变化的研究日渐增多。Yu等［7］通过GIMMS NDVI对青藏高原植被生长季动态进行研究，发现在气温的影响下，高原春季物候在20世纪90年代中期发生了延迟。Chen等［8］虽然表示并未在野外检测中观察到相关现象，但整合了草原退化导致的变冷效应及冻土冻结提前对春季物候延迟现象的解释。Yi等［9］认为可能是青藏高原春季的较高气溶胶指数影响了NDVI的监测值，造成春季物候延迟的现象。Piao等［10］也使用GIMMS NDVI得到了青藏高原春季物候在1982—1999年以0.88 d·a-1的速率提前，随后春季物候出现了延迟，与Yu等［7］的研究结果一致。Zhang等［11］基于MODIS、SPOT与GIMMS NDVI数据集，揭示在1982—2011年间青藏高原春季物候以1.04 d·a-1的速度持续提前。并且，Zhang等［11］对三种NDVI数据集进行了不同时间尺度的年际变化分析，认为GIMMS NDVI存在质量问题。但Shen等［12］在随后利用SPOT与MODIS NDVI指出，尽管青藏高原近10年春季温度上升，但在区域尺度上植被的返青期并未显著提前。关于青藏高原春季物候的提前、静止或延迟，学界一直缺乏共识［13］，不同研究对于青藏高原植被评价的结论存在差异性，而这种不一致性很可能与研究所选取的卫星产品相关［11，14］。不同时间、空间分辨率的多源NDVI数据通常作为直接表征区域植被绿度的指标，用来探究植被变化，同时也是植被物候提取的数据来源，因此基于多源NDVI数据集对青藏高原地区植被变化特征进行评价并探讨不同NDVI数据之间的差异性尤为重要。但是，现有研究对于多源遥感数据在植被变化特征差异性方面的探讨缺少针对性：一方面，研究更多地探讨不同卫星衍生指标之间的差异，而非聚焦同种指标［14-16］，如归一化植被指数之间的差异；另一方面，即使基于归一化植被指数进行研究，所使用的NDVI数据集种类并不全面［11，17］，研究范围也往往并非针对青藏高原［18-21］，而光谱植被指数受到的区域限制因素众多，存在极大的空间异质性。所以，本研究采用多种不同时空尺度的归一化植被指数在青藏高原范围内，对植被状况及其长时间变化进行分析，探讨和对比不同数据之间的差异性。

鉴于青藏高原植被在地-气能量交换、水文过程等的重要作用，而过去研究缺乏多源数据在青藏高原植被变化评估中的差异性和不确定性，为此本文围绕不同时空分辨率的多源卫星遥感数据在青藏高原植被变化评估的差异性如何的科学问题开展相关研究。基于MODIS（250 m、500 m、1 km、5.5 km）、GIMMS 8 km和SPOT 1 km NDVI数据集，开展不同时空尺度下青藏高原地区的植被变化特征评价。具体包括：使用最大值合成方法（MVC）探究青藏高原植被2000—2014年的空间分布；使用Theil-Sen趋势估计和Mann-Kendall趋势检验探究时间变化特征；最终评价不同时间空间分辨率的多源卫星遥感数据在青藏高原植被变化特征研究上的差异性。

1 研究区概况

青藏高原是继南北极后的“第三极”，平均海拔高度超过4 000 m，面积约为2.5×106km2，是世界上海拔最高、面积最大的高原［22］。青藏高原范围为73°29′56″～104°40′20″E，25°59′37″～39°49′33″N［23］，山脉纵横，地域广阔。青藏高原气温差异大，年平均气温为-10～15℃［24-25］［图1（a）］。自1901年以来，气温增加趋势最高可达0.15℃·（10a）-1以上［图1（b）］，高原西北部增温趋势更为显著。降水主要集中在夏季和秋季，而冬、春两季相对干燥［26］。降水趋势在整个高原尺度上表现出较小的季节和空间波动，但总体上略有增加［26］。年降水量由东南向西北变化，从高于1 000 mm减少至低于100 mm［27］。复杂的气候特征造就了高原丰富的植被类型分布。相对更为温暖湿润的高原东南部与东部边缘，主要分布有森林和灌木，向西植被逐渐发展为高寒草甸，干燥的高原中部与东北部区域覆盖有温带与高寒草原，荒漠则分布在高原西北部的干旱地区［26］。

图1 青藏高原1901—2018年青藏高原年平均气温（网格温度数据来自CRU）（a）及其变化趋势（b）Fig.1 Mean annual air temperature（grid temperature data is from CRU）of the Tibetan Plateau from 1901 to 2018（a）and its change trend（b）

2 数据与方法

2.1 数据来源

本文研究数据来自MODIS、SPOT-VGT、AVHRR卫星遥感源传感器，涉及包括250 m、500 m、1 km、5.5 km、8 km共5种空间分辨率，跨越了从10-2km2至近102km2共5个数量级的空间尺度（表1）。

表1 卫星遥感数据集介绍Table 1 Description of satellite remote sensing datasets

本文使用搭载于EOS/Terra卫星上的中等分辨率成像光谱仪（MODIS）反演的植被指数产品，包括MOD13Q1、MOD13A1、MOD13A2、MOD13C1产品，空间分辨率为250 m、500 m、1 km、5.5 km，其提供了全球的植被状况与具体时空变化的信息，可以用于观测陆地植被光合作用的变化。MOD13Q1、MOD13A1与MOD13A2数据集来源于Google Earth Engine平台，MOD13C1数据来源于美国宇航局（https：//ladsweb.modaps.eosdis.nasa.gov/）。GIMMS NDVI 3g数据（http：//ecocast.arc.nasa.gov/data/）由NOAA卫星搭载的AVHRR传感器获取，是目前在植被研究中使用最广，涵盖时间序列最长的NDVI产品。SPOT NDVI产品（http：//land.coperni⁃cus.eu/global/products）为GIOGL1_ATBD_ND⁃VI1km-V2.2数据集，数据由SPOT-VEGETATION提供（2014年后由PROBA-V提供），基于最大值合成法合成以10天为周期的最佳NDVI值，减少了噪声的影响。

2.2 研究方法

2.2.1 最大值合成法与标准差

使用最大值合成法分别将不同时空分辨率的MODIS、SPOT、GIMMS NDVI数据合成2000—2014年多年季节性NDVI［28］，以探究青藏高原多年季节性植被空间分布状况。为了保证不同数据源NDVI之间的可比性，反映不同数据来源NDVI的差异性，使用最近邻插值将所有数据集的空间分辨率转换为250 m［29］，计算6种不同NDVI数据集之间的标准差［14］。

2.2.2 Theil-Sen趋势估计和Mann-Kendall趋势检验

由于NDVI时间序列通常不满足数据独立且正态分布的参数趋势检验要求［15］，在探究青藏高原植被变化时，使用了两种非参数方法——Theil-Sen趋势估计和Mann-Kendall趋势检验，以探究NDVI序列变化的趋势［30］。

研究使用Theil-Sen斜率估计［31］作为NDVI随时间变化程度的衡量指标，更为准确地定量化描述植被变化程度［32］。Theil-Sen斜率估计是一种稳健的趋势统计方法［30］，虽然与线性最小二乘回归方法相关，但该方法基于非参数统计，通过中值计算斜率，与前者相比能够抵抗噪声的干扰，对数据中的离群值有着良好的规避能力［6，18］。

NDVI时间序列趋势的显著性使用非参数的Mann-Kendall趋势检验，Mann-Kendall趋势检验比较样本数据的相对大小并计算得到统计量S及其方差Var（S），然后计算检验统计量Z，来确定数据变化的趋势显著性［15，33］。

Theil-Sen趋势估计和Mann-Kendall趋势检验相结合能更加有效地探究NDVI的时间序列变化，这两种方法被广泛应用于地球科学的许多领域［34］。

3 结果与分析

3.1 青藏高原NDVI空间分布特征

通过对不同遥感数据集的青藏高原2000—2014年季节性NDVI最大值空间分布的分析，结果表明，青藏高原区域植被分布情况有着明显的空间异质性。各季节青藏高原NDVI空间分布均呈现东高西低的特征，NDVI由东南到西北逐渐降低，空间分布格局明显（图2）。青藏高原植被分布格局的形成可能与植被对于温度的敏感性相关，高原年均气温的分布同样是东南部相对较高。虽然植被生长与气温之间的关系十分复杂且具有显著的区域差异［35］，但对于高纬度地区与青藏高原来说，气温依然被广泛认为是植被生长的重要限制因素。

图2 青藏高原生长季NDVI的空间分布Fig.2 Spatial distribution of NDVI in the growing season on the Tibetan Plateau

研究区NDVI在0～1的范围内均有分布（图2）。不同季节间，NDVI的主要分布范围有所差异。青藏高原夏季NDVI、秋季NDVI与生长季NDVI均呈现双峰式集中分布，夏季NDVI与生长季NDVI在0～0.3与0.7～0.9处集中。秋季NDVI主要在0～0.3与0.6～0.8相对集中，但集中在0.6～0.8的像元占比远小于前者。春季NDVI整体较低，约58.06%的像元NDVI集中在0～0.2范围。

不同NDVI数据集之间的标准差由北向南逐渐增加（图3）。秋季的高原整体标准差最小。春季，6种NDVI数据集差异最为明显，较大的标准差出现在高原东南部及南部。在夏季与生长季，不同ND⁃VI数据集主要在高原南部出现较大的差异。SPOT NDVI分布更多的集中在0.1～0.4之间（图4）。SPOT NDVI在0.1～0.4之间像元占比最大，平均多于其他数据集4.11%。在0.4～0.7的NDVI分布区间内，GIMMS NDVI较其他数据平均多约1.83%。而MODIS数据普遍在0.7～1区间的NDVI像元占比更大。6种NDVI数据集中，250 m分辨率的MODIS NDVI在0.7以上分布占比最具优势。而随着空间分辨率的降低，MODIS NDVI数据在NDVI>0.7的面积逐渐减少。

图3 6种NDVI标准差的空间分布Fig.3 Spatial distribution of standard deviation of six NDVI：spring（a），summer（b），autumn（c）and growing season（d）

图4 青藏高原各季节NDVI的频率分布Fig.4 Frequency distribution of NDVI in each season on the Tibetan Plateau

3.2 青藏高原NDVI时间变化特征

3.2.1 区域平均NDVI变化趋势

通过对6种NDVI数据集在青藏高原范围内进行区域平均计算和Theil-Sen趋势估计分析（图5）。结果表明，在青藏高原的区域尺度上，SPOT NDVI与MODIS 250 m、500 m、1 km NDVI数据集除秋季外，均表现显著的绿化趋势（表2），MODIS 5.5 km与GIMMS NDVI数据集分别在春季、夏季、生长季呈现显著绿化。

表2 青藏高原平均NDVI变化趋势［单位：（10a）-1］Table 2 Trends in average NDVI changes on the Tibetan Plateau［unit：（10a）-1］

图5 青藏高原平均NDVI变化Fig.5 Average NDVI changes on the Tibetan Plateau

对于不同数据源而言，NDVI的显著变化趋势相差较大，从0.0092（10a）-1到0.0184（10a）-1。即使季节相同，不同NDVI数据集之间的显著绿化趋势差异最大可以达到0.0090（10a）-1（生长季）。SPOT NDVI在四个季节均表现显著的绿化趋势，且平均绿化趋势达到0.0162（10a）-1，高于其他数据平均变化趋势0.0085（10a）-1，绿化显著且迅速。GIMMS NDVI仅在2000—2014年间的生长季变化分析中具有0.0092（10a）-1的显著绿化趋势。具有最高空间分辨率的MODIS 250 m NDVI数据集往往表现出最小的显著变化趋势值。

就季节而言，除秋季外其他季节均具有显著的变化趋势。各个NDVI数据集在15年间春季的显著变化趋势介于0.0098～0.0184（10a）-1之间。多数数据集夏季显著绿化趋势明显高于其他季节，夏季平均显著绿化趋势达到0.0140（10a）-1，而青藏高原夏季的NDVI往往深刻影响整个生长季的植被变化。

3.2.2 青藏高原多年植被变化趋势分析

青藏高原植被的长期变化存在空间差异，并非完全遵循在区域平均尺度上分析得到的规律均匀发展［35］，受到复杂的生物地球物理与生物地球化学过程影响，而在气候变化的背景下，这些循环过程变得尤为复杂。此外，植被受到气温、降水、海拔、经纬度等环境变量的影响，植被分布状况具有明显的空间异质性，因此对于植被变化的评估十分困难。但其是探究陆地生态系统对气候变化响应的关键［36］。对于青藏高原这一巨大的地理单元，在探究区域植被的具体变化趋势时，应考虑广阔的高原因自然条件的不同而导致的差异性与复杂性。因此，进一步基于年际季节性最大NDVI数据逐像元进行Theil-Sen趋势估计并结合Mann-Kendall趋势检验，对青藏高原植被随时间变化趋势及其相应空间格局进行探究。

在2000—2014年15年间，青藏高原大部分区域的NDVI呈现上升趋势，通过Theil-Sen趋势估计确定的植被变化斜率在平均74.04%的范围呈现正值，反映高原植被逐渐绿化，其中29.29%的像元表现为显著的绿化趋势（图6～9）。但是绝大部分区域NDVI变化趋势较小，集中在0～0.02（10a）-1。绿化趋势较大的区域基本分布于高原的东部，具体的分布格局因数据源与研究季节而异。

图6 青藏高原春季NDVI变化趋势的空间分布Fig.6 Spatial distribution of trends in NDVI changes in spring on the Tibetan Plateau

不同数据集对探究青藏高原NDVI时间变化趋势所造成的差异难以忽视。数据的变化趋势值标准差说明了6种NDVI普遍在青藏高原南部及东南部差异较大（图10），在春季差异最为明显，且广泛分布在高原东南部。在其他季节中，不同数据集的变化趋势值标准差相对较小，趋势斜率的差异主要分布在青藏高原南部边界。与区域平均NDVI长时间序列的变化趋势结果一致（表2），SPOT NDVI反映的植被绿化趋势相对于其他数据更加显著，91.18%的像元在生长季中呈现绿化，29.01%的像元呈现显著绿化，分别超出MODIS 250 m、500 m、1 km NDVI、MODIS 5.5 km NDVI与GIMMS ND⁃VI 7.16%、5.53%、4.65%、2.87%与15.89%。显著绿化趋势达到0.02（10a）-1以上有23.86%，超出趋势分布具有较高一致性的MODIS数据11.25%，超出绿化像元占比最少的GIMMS数据16.76%（图11）。其次，春季SPOT NDVI呈现显著增长的像元占比最大，达到33.86%，大范围分布在青藏高原东南部；GIMMS NDVI呈现显著增长的像元则在春季占比最小，而显著下降的像元达到了7.88%，占比最大，分布在高原南部边界。这可能是春季相比其他季节标准差更大的原因。在其他季节中，GIMMS NDVI也具有呈现植被绿化趋势较小，褐化像元占比明显高于其他5种数据集的特点。GIMMS NDVI显著褐化的像元占比平均达到5.83%。MODIS NDVI随着分辨率的提高，显著绿化趋势的像元占比与绿化趋势却往往逐渐降低。总体上，逐像元的NDVI长时间变化趋势与青藏高原区域平均的多年植被变化相一致，GIMMS NDVI表现高原植被绿化范围与绿化趋势最小。绿化像元最多且趋势最大的SPOT NDVI，在区域平均NDVI变化斜率也最大。

图10 6种NDVI变化趋势标准差的空间分布Fig.10 Spatial distribution of standard deviation of six trends in NDVI changes：spring（a），summer（b），autumn（c）and growing season（d）

图11 青藏高原各季节NDVI变化趋势的频率分布Fig.11 Frequency distribution of trends in NDVI changes in each season on the Tibetan Plateau

4 讨论

MODIS NDVI与SPOT NDVI数据集是目前评估长时间序列植被变化的可靠数据来源，有着广泛的一致性。在青藏高原区域尺度NDVI时间变化分析与像元尺度NDVI时间变化分析中，MODIS 250 m、MODIS 500 m、MODIS 1 km与SPOT NDVI数据集显示青藏高原植被整体呈现绿化趋势，MODIS 5.5 km NDVI也在春、夏两季呈现显著的增长趋势。这与青藏高原正在变得温暖的认知相符合［1］，气温的上升往往被认为有利于植被的生长［6］。气候变暖主导了青藏高原的绿化趋势，可以通过促进植物光合作用及延长生长季来增加寒冷地区的植被绿度，对植被变化产生积极的影响［37-38］。但另一方面，在青藏高原的干旱与半干旱区域，高寒草甸与高寒草原对水分条件十分敏感，变暖引起的蒸散发增加可能致使植被受到干旱胁迫［39-40］，从而导致植被退化。冬季的增温过强可能使得植物休眠延迟，从而推迟春季物候，缩短生长季［7］，导致植被绿度降低。这能够解释在青藏高原的部分地区，NDVI呈现下降趋势的现象。

图7 青藏高原夏季NDVI变化趋势的空间分布Fig.7 Spatial distribution of trends in NDVI changes in summer on the Tibetan Plateau

多源NDVI数据集在青藏高原植被的空间分布与长时间序列变化中都存在着难以忽视的差异。尽管MODIS NDVI与SPOT NDVI在青藏高原区域范围内出现了较为清晰一致的变化模式，但具体的绿化幅度及绿化范围因分析的数据集而异。在探究2000—2014年的植被变化时，GIMMS呈现褐化趋势的像元占比显然多于不同分辨率的MODIS NDVI与SPOT NDVI，区域尺度上的GIMMS NDVI的变化斜率也明显小于二者。前人研究指出，青藏高原地区GIMMS NDVI要小于SPOT NDVI与MODIS NDVI，而GIMMS NDVI与SPOT NDVI、MODIS NDVI的差异大于SPOT NDVI与MODIS NDVI的差异［19］。SPOT NDVI与GIMMS NDVI数据在反映青藏高原1998—2006年植被变化时，呈现了较大的异质性，SPOT NDVI所反映的绿化范围与绿化趋势明显大于GIMMS NDVI［11］。此外，许多关于青藏高原植被动态的研究也得出了基于GIMMS NDVI所确定的春季物候变化趋势与基于MODIS NDVI、SPOT NDVI、MODIS增强型植被指数（EVI）或日光诱导叶绿素荧光（SIF）等其他遥感数据集所确定春季物候变化趋势不一致的结论［1，12，14，41-42］。部分学者认为这种物候变化趋势的不一致可能与GIMMS NDVI的降低有关［11，41］，并由此推断在青藏高原区域使用GIMMS NDVI探究植被变化会带来研究结果的不确定性。

图8 青藏高原秋季NDVI变化趋势的空间分布Fig.8 Spatial distribution of trends in NDVI changes in autumn on the Tibetan Plateau

不同数据集对植被变化分析的差异性可能由许多因素导致，例如空间分辨率、时间分辨率、传感器设计、校正方法、投影系统、地理定位误差之间的差异［15，17-18，21］。首先，从8 km到250 m的空间分辨率，随着观测尺度的变化，尺度依赖性和尺度效应的影响随之体现［19，43］。其次，卫星平台及其所搭载的传感器对NDVI产品有着重要影响。GIMMS NDVI、SPOT NDVI、MODIS NDVI数据集均来自于多传感器系统［21］，但相比AVHRR，其他二者的传感器拥有更稳定的轨道与改进植被监测的光谱配置［20］。MODIS传感器与SPOT-VGT传感器在红外与近红外的波段设计更为接近［19］，这可能也是其获取的NDVI更为接近的原因之一。来自于AVHRR传感器的GIMMS NDVI是目前时间序列最长的NDVI产品，但AVHRR系列传感器最初并非为植被研究设计［18］，缺少机载校准、光谱规格不连续［44］且光谱结构也难以进行更为精准的大气校正［18］。但大气成分的散射与吸收是NDVI监测的重要误差来源。有研究指出GIMMS NDVI与其他卫星遥感数据的差异可能受青藏高原春季气溶胶浓度增加的影响［11］。

5 结论

研究结果表明，SPOT NDVI所反映的植被绿化呈现显著且迅速的特征，27.44%的像元呈现显著绿化趋势，分别超出MODIS NDVI与GIMMS ND⁃VI 4.10%、15.89%；在生长季，SPOT显著绿化趋势达到0.0182（10a）-1，高于其他显著变化的数据0.0078～0.0090（10a）-1。MODIS NDVI随着分辨率的提高，显著绿化的像元占比与绿化趋势却逐渐降低；MODIS数据之间显著绿化的像元占比相差不足2.80%。GIMMS NDVI显著褐化的像元占比平均达到5.83%，超出其他数据集3.37%～5.51%。在春季，GIMMS NDVI表现显著褐化像元占比最大，达到7.88%，与显著绿化的像元占比相当。GIMMS的区域平均NDVI具有最小的显著绿化趋势，为0.0092（10a）-1，并且在春、夏两季呈现不显著的褐化趋势。基于GIMMS NDVI探究青藏高原植被变化特征时，特别是对春季物候研究，可能会造成结果具有较大的不确定性。而SPOT与MO⁃DIS NDVI体现了较高的一致性，可以互为补充，探究高原植被变化。