孙倩倩,刘 超
(安徽理工大学 空间信息与测绘工程学院,安徽 淮南 232001)
近一个世纪以来,全球气候变暖给陆地生态系统带来了重大的影响,如温度升高、冰川融化、海平面上升等一系列问题[1-3]。植被是组成陆地生态系统的重要成分,是连接大气、土壤和水分的枢纽[4-5],在全球能量循环、水文循环及生物化学循环过程中扮演着十分重要的角色,其生长趋势受一个地区自然因素多样性的长期影响[6],因此研究植被对气候变化的响应对于深度了解生态系统动态至关重要。黄土高原的生态环境脆弱,物种丰富,气候特征时空差异明显,研究黄土高原地区连续植被动态以及过去三十年植被物候变化对气候变化的响应为其他地区的生态环境变化可提供丰富的理论基础[7]。
目前,越来越多的研究者对黄土高原的植被动态变化进行了深入研究,植被生长对气候变化的响应是一个复杂的交互作用的过程,并且不同植被对气候变化的响应时间各不相同[8-9]。植被生长指数和气候因子的长时间序列均是非平稳、非线性的,包含不同的频率信息[10],如月际、季节、年际、长期或短期的变化[11-12]。对于长时间序列数据集,如果同一特征出现在特定的频率间隔内(例如,每两年发生一次干旱),则该特定的频率间隔称为时间尺度[13]。这些不同的频率分量可以提供详细的多尺度信息,这些信息有助于在不同地理时间序列中产生趋势并影响其相关性[14]。
本文在非平稳、非线性的长时间序列的基础上,选择一种多尺度分解的方法探讨气候变化与植被动态的关系,分析在不同时间尺度下气候对植被产生的作用。多尺度分解方法在分析植被生长的年际动态和长期趋势方面是具有有效性的,并且已经有研究开始使用多尺度分解方法探索植被生长和气候变化的总趋势以及不同尺度下气候因子与植被生长的关系[15-16]。例如,Qi等[16]整合了集合经验模态分解 (Ensemble Empirical Mode Decomposition, EEMD) 和剩余趋势方法去研究在多时间尺度上中国丝绸之路经济带的植被生长与气候变化及人类活动的关系。Liu等[17]利用具有自适应噪声的完整集合经验模式分解算法和快速傅里叶变换(Fast Fourier Transformation,FFT)方法探讨了内蒙古地区1982-2015年植被的时空变化规律及其与温度和降水的关系。同时,频谱分析可以将时间序列中的不同时间尺度的分量分离出来,获得植被生长的不同周期[16,18]。
本文以黄土高原为研究区域,采用改进的具有自适应噪声的完整集合经验模式分解算法(Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,Improved CEEMDAN)[19]对1982-2013年植被指数和气候因素的长时间序列进行分解,再利用FFT获得每个分量的时间尺度,并在不同尺度下,研究黄土高原近三十年的植被生长和气候变化之间关系,为黄土高原植被保护、生态恢复及生态建设等提供必要科学依据。
黄土高原 (100°54′~114°33′E, 33°43′~ 41°16′N) 是中国四大高原之一, 是世界上最大的黄土沉积区,位于中国中部偏北,黄河中上游,南北长约750 km,东西长约1000 km,总面积约为640 000 km2;海拔高度800~3000 m。其地貌复杂、地质环境丰富,水土流失非常严重且生态环境脆弱[20]。该区域属于大陆季风气候,冬季寒冷、夏季严热,年均温度3.6~14.3℃;降水呈现明显的时空差异(梯度变化),年降水量150~750 mm;蒸散发量大。该区域的自然植被自东北至西南呈条带状变化,依次为暖温带夏绿阔叶林、森林草原、干草原及中温带荒漠草原。森林/草原的生长季节通常是从4-10月。如图1,为黄土高原所处位置及植被类型和气象站点分布情况。
图1 黄土高原所处位置以及植被类型和气象站点分布情况
本研究使用的是1982-2013年的第三代全球清单建模和制图研究NDVI (GIMMS NDVI3g) 数据集。这些数据来自NOAA/AVHRR传感器,空间分辨率为8 km,每半个月采集一次数据。GIMMS NDVI3g数据集已经过仔细校正,以最大限度地减少校准损失,轨道漂移或火山爆发等有害影响[10]。目前其已经被认为是可用于长期NDVI分析的最佳数据集之一,并且已被证明有助于准确调查植被活动变化的真实特征[21]。为了减低云和气溶胶散射对数据的影响,本文所用的月数据是利用最大值合成方法[22]所获得的。
NDVIi=Max(NDVIi1,NDVIi2)
(1)
式中,NDVIi是第i个月的NDVI值;NDVIi1是i个月的前15天的NDVI值;NDVIi2是i个月的后15天的NDVI值。为了避免冬季和早春极端气候对植被的影响,本文选择4-10月代表生长季,并使用生长季范围的数据进行研究,共224个月。
气候数据包括月平均气温和月累计降水数据,来自于中国气象数据网 (http: //data. cma. gov. cn/) 的地面气候资料数据集,选用黄土高原的52个气象站点数据。为了确保插值数据的准确性,更好地凸显气候因素和NDVI之间的相关性[14],本研究采用反距离加权法进行插值,得到1982-2013年空间分辨率为1/12°的气候栅格数据集。
土地覆盖数据采用的是MODIS 产品 (MCD12C1),空间分辨率为0.05°,有17种土地覆盖类型。为了获得与NDVI数据同样椭球参数、投影方式和空间分辨率的植被分类数据,统一数据集形式,需要进行投影和像元的重采样,采用最邻近法采样,将0.05°数据重采样为1/12°。本研究选择2001-2012年黄土高原地区的混交林、草地、稀疏植被作为研究对象,并利用平均像元的方法将2001-2012年间土地覆盖类型没有发生变化的像元提取出来,再获得不同植被类型的平均降雨和温度数据,进行多尺度分解和回归分析。
Hilbert-Huang变换 (HHT) 是Huang等[23]提出的一种信号处理方法,该方法包括两个部分:EMD和Hilbert 谱信号分析方法。EMD算法属于自适应的局部时频分析方法,主要用于分析源自非线性系统的非平稳信号,可以将复杂信号分解为多个内部模式函数(Intrinsic Mode Function,IMF)分量和残余项Residual,分别表示原始信号的不同时间尺度[23-24]。EMD的基本思想是:给定某一信号,先获得信号的极值点,通过插值获得信号的上下包络,得到上下包络线的均值,与均值的差得到分解的第一层信号,重复此过程,直到将信号f(t)分解成有限个基本模式分量imfj(t)和残余项rn(t)的组合,其表达式为
(2)
Improved CEEMDAN是由Colominas等[19]提出的一种新的分解算法,该方法是在EMD的基础上,进行改进的一种气候变化领域的先进时频分析方法[15-16,19],在平均信号第一模式加入一定大小的白噪声,使得信号在不同的尺度之间具有一定的连续性,从而产生新的极值点,可以在一定程度上消弱频率混叠现象。设NDVI和气候变量(降雨、温度)的长时间序列为x,其算法描述如下:
(1)在长时间序列x中加入一定大小的白噪声后形成的待分解信号如下:
xi=x+β0E1(wi)
(3)
式中,xi为加入噪声的待分解信号,i为添加白噪声次数,i=1, 2, …,n;β0为噪声的大小,wi为零均值单位方差白噪声;E1(wi) 为白噪声wi的第1个EMD分量。高斯白噪声的标准差范围为 0.01~0.4,本文中对于降雨和温度加入的白噪声为0.1,NDVI和加入的白噪声为0.3,噪声添加的次数为50。
(2)采用EMD方法对n个待分解信号进行分解,得到Improved CEEMDAN方法的第1个IMF分量,表示为:
IMF1=x-〈M(xi)〉
(4)
式中,IMF1为第1个IMF分量;〈·〉代表取平均值,M(·)为满足IMF筛选阈值的包络线的局部平均值。
(3)构造第二次待分解信号〈M(xi)〉+β1·E2(wi)=r1+β1·E2(wi),采用EMD方法对信号进行分解,得到Improved CEEMDAN的第2个IMF分量为:
IMF2=r1-〈M(r1+β1·E2(wi))〉
(5)
式中,r1为第一次分解获得的残余项,β1为噪声的大小,E2(wi) 为白噪声wi的第2个分量。
(4)依次类推,第k次待分解信号为rk-1+βk-1·Ek(wi),得到的Improved CEEMDAN的第k个IMF分量为:
IMFk=rk-1-〈M(rk-1+βk-1·Ek(wi))〉
(6)
式中,rk-1为第k-1次分解获得的残余项,βk-1为噪声的大小,Ek(wi) 为白噪声wi的第k个分量。
(5)重复步骤(4),直到残余项Residual不符合分解条件,则停止分解。
当Improved CEEMDAN分解结束时,再利用FFT对分解的每个分量做频谱分析,获得每个分量的周期。最后一个残差项作为长时间尺度,频率显示为1 Hz的作为年际尺度,再根据不同尺度进行分析。
本文以黄土高原地区的3种自然植被将NDVI、降雨和温度分别划分为三类数据,图2是黄土高原1982-2013年生长季(4-10月)混交林的NDVI、降雨与温度的长时间序列,横轴表示长时间序列的长度 1982-2013年,纵轴从上到下表示的分别是NDVI值,降雨量(mm)和温度(℃)。从图2可以看出,NDVI的生长除了后期有一点小的波动,整体规律是比较明显的,具有一定的周期性,降雨的变化是一种复杂而随机的,没有规律可循,而Temperature具有较好的周期性。因此,在使用Improved CEEMDAN方法时,应针对不同的数据类型添加不同的高斯白噪声,降雨和温度加入的白噪声为0.1,NDVI和加入的白噪声为0.3。
图2 黄土高原1982-2013年混交林的NDVI、降雨与温度的长时间序列展示
利用Improved CEEMDAN方法对NDVI、降雨和温度长时间序列分别进行分解,获得若干个IMF,再利用FFT对每个分量做频谱分析。图3表示的是采用Improved CEEMDAN 分解温度的结果,一共有7个IMF,分解后的分量再利用FFT做频谱分析,获得每个分量的周期信息。图4表示的是采用Improved CEEMDAN 分解降雨时间序列的结果,降雨的长期变化是一种复杂而随机的,没有明显的规律,这里降雨一共被分解成了7个IMF分量。图5表示的是采用Improved CEEMDAN 分解NDVI数据的结果,从图5中可以看出,NDVI数据一共被分解成了6个IMF分量,与原始NDVI数据相比(见图2),经多尺度分解后的NDVI时间序列更加的平滑,并且具有非常明确的趋势性,能够直观地反映出植被生长的更迭变化情况。不同的分量具有不同的频率,反应了不同的周期信息,本文根据FFT获得每个分量的周期信息,选取了两个尺度进行分析,将频率等于1 Hz的分量看作是年际尺度,最后一个残余项Residual当作时长时间变化尺度,分析在两个尺度下植被生长和气候变化(降雨、温度)之间的关系。这里所采用的示例数据为混交林的数据,并对其他两类数据也做相同的处理。
左图纵轴表示NDVI分解后的值的波动范围;右图纵轴表示振幅。
左图纵轴表示NDVI分解后的值的波动范围;右图纵轴表示振幅。
左图纵轴表示NDVI分解后的值的波动范围;右图纵轴表示振幅。
为了揭示黄土高原不同气候带下不同时间尺度植被生长与气候因素的关系,选取3种具有象征性的自然植被,用于研究不同气候因素对不同植被的影响,并采用Improved CEEMDAN方法对1982-2013年间生长季平均NDVI、温度和降水时间序列数据进行多尺度分解。通过在分解时针对不同数据加入不同大小的噪声来构建相关实验。此外,本文根据频率信息在IMF分量中提取出年际尺度和Residual,再将不同尺度下的NDVI分别与对应的温度和降雨进行回归,获得不同尺度下气候变化对植被生长的影响。
在年际和Residual两个尺度下,温度和降雨对植被的影响是不同的,并在显著性水P<0.01和P<0.05的检验下进行结果分析。由于降雨的复杂性和随机性,从表1可以看出,年际尺度下,降雨与NDVI具有较好的相关性,其中,降雨对草地的影响最大,混交林的影响最小。在长时间尺度(Residual)下,降雨对混交林和稀疏植被的生长为负影响,但是对草地的生长仍然是促进作用,这可能是因为草地是一年生植被[25]。年际尺度和长时间尺度下,温度对不同植被的生长都表现出了较好的相关性,其中最明显的是混交林(见表2)。从总体上来看,不同尺度下温度对不同植被的影响大于降水。
表1 在不同尺度下不同植被生长与降雨的相关性
表2 在不同尺度下不同植被生长与温度的相关性
本文主要研究不同尺度下,黄土高原地区不同气候对不同植被产生的影响,假设气温和降水对植被的影响是相互独立的,选取黄土高原地区3种自然植被,研究自然植被和气候因素的相关性。利用Improved CEEMDAN方法进行多尺度分解,分为两个尺度进行研究,分析不同尺度下植被生长与温度和降水之间的关系。由于Improved CEEMDAN算法具有自适应性,因此本文无法基于像元进行分解,而是对相同类型气候带和植被像元取平均[17]。不同的气候因子对相同的植被具有不同的影响,不同植被对相同气候因子的反应也不同[5]。研究结果表明在不同的尺度下,温度和降雨所呈现的结果不同[24],研究发现植被与气候变化之间的关系随时间尺度的变化而变化,并且Residual尺度下,温度与植被之间的总体相关性大于降水。这与本研究的结果一致,说明气温是影响植被生长变化的一个重要因子,长期温度升高可以促进植被的光合作用,进一步促进植被生长[16],并且韩辉邦等[27]在黑河流域和王海军等[28]在中国西北地区也证明了温度的变化对植被生长的影响更大,这也验证了本文的结论。因此在多尺度下分析植被和气候因素的关系,揭示了黄土高原地区不同尺度下植被生长与气候之间的联系,可以更好地探索植被与气候等因素的关系。同时,监测植被动态及其与气候变化的关系被认为对全球变化研究和生态环境保护非常重要,本研究的结果有助于更好地理解多尺度下植被动态与气候变化之间的关系,对植被生长进行更精确的监测,并为全球脆弱生态系统的植被保护和恢复提供一些科学参考。