基于小波变换的卧龙国家级自然保护区植被时空变化分析

朴英超，关燕宁，张春燕，郭　杉，阎保平

1 中国科学院遥感与数字地球研究所，北京　100101 2 中国科学院计算机网络信息中心，北京　100190 3 中国科学院大学，北京　100049



基于小波变换的卧龙国家级自然保护区植被时空变化分析

朴英超1,2,3，关燕宁1，*，张春燕1，郭杉1，阎保平2

1 中国科学院遥感与数字地球研究所，北京100101 2 中国科学院计算机网络信息中心，北京100190 3 中国科学院大学，北京100049

摘要:提出一种基于小波变换从长时间序列、大范围遥感数据中快速、自动化检测植被动态变化的方法。以MODIS 500m空间分辨率，16d合成的NDVI数据为数据源，对受到2008年5月12日汶川地震严重影响的卧龙国家级自然保护区内2003年至2012年的植被动态变化进行时空分析，为保护生态多样性及生态系统的稳定性提供依据。研究表明：1)地震后保护区内植被指数减少的面积大范围增加，且波动较震前更为明显，统计分析结果能够更为直观地反映地震及其次生灾害等极端现象对该地区植被的破坏程度；2)保护区内植被指数极值变化多发生在夏季或秋季，较低海拔地区极值变化多发生在夏季，而在高海拔地区则多发生在秋季；3)在大熊猫最适宜栖息的区域(2600—2800m)植被指数极值减少量大于0.4的范围大于增加量大于0.4的范围，反映出植被在震后的恢复状况并没达到理想的水平。同时发现在该海拔区域范围内植被指数减少的面积在春夏两季较大，表明在该时间段卧龙地区大熊猫最适宜生存区域的植被情况较为不稳定，需更为关注其动态，采取适当的保护措施。

关键词：植被变化；小波变换；时间序列分析；卧龙国家级自然保护区；汶川地震

Analysis of temporal and spatial changes in vegetation cover using wavelet

植被是反映生态环境变化的重要指标，也已经成为生态研究领域的一个热点[1- 3]。自然因素(降水，气温等)、人为因素(砍伐，城镇化等)或其他特异因素(火灾，地震等)都能够影响植被的生长状态[2,4]。卧龙国家级自然保护区是中国大熊猫的重要栖息地，该地区的植被是大熊猫赖以生存的根本保障，因此研究栖息地植被的动态变化能够为保护生物多样性及生态系统功能的完整性和稳定性提供依据。

归一化植被指数(NDVI)被广泛应用于植被变化的检测中，因其能较好的反映植被的生长状态及植被覆盖度，其时间变化曲线可以反映季节和人为活动的变化。利用NDVI进行植被动态变化检测的研究也在不断发展中,包括统计方法[5- 7]，频谱分析[8- 9]等等，但对长时间序列分析方法的研究还是很有限的[4,10- 12]。之前的研究包括对阈值的设定或置信区间的计算并不普遍适用且准确度不高，时间序列波动的拟合效果也并不是十分理想。时间序列的变化检测方法本身面临很多挑战，包括：分析结果通常并不直观，物候变化和突变无法明确区分以及在大范围的数据上进行变化检测通常要耗费大量的时间和成本。

小波变换的方法最初应用在信号和图像处理领域，并被认为是图像处理的重要方法之一[13]。随着其在地球物理领域的应用[14]，小波分析在时间序列数据和空间数据分析中得到广泛关注[15]。与其他将光谱转换到频域的方法(如：傅里叶变换)相比该方法更具优势，因为小波变换对数据的稳定性要求不高而且可以同时在时域和频域上进行分析[3]。目前小波变换在植被动态变化检测上的应用主要集中于发现特征样点变化或提取物候特征的关键因子方面[6-7,16]，而在大范围区域空间上进行植被动态变化检测的研究十分有限。

卧龙国家级自然保护区建立于1963年，主要保护西南高山林区自然生态系统及大熊猫等珍稀动物，是我国大熊猫的重要栖息地，大熊猫的数量约占全国总数的10%。2008年5月12日汶川大地震给该保护区带来严重的影响和破坏。主震及多起强余震引起的山体滑坡、崩塌以及泥石流等次生地质灾害120多处，造成森林植被损毁，保护区内大熊猫栖息地的动、植物生境遭到严重破坏，区域生态环境受到重大影响[17- 19]。

栖息地是野生大熊猫赖以生存的物质基础，研究和评估栖息地质量是有效保护和管理大熊猫的重要手段[20]，而其中地形、土壤、气候和植被是大熊猫生存的重要限制因素[21]。 本文提出一种从长时间、大范围遥感影像中快速、自动化地检测植被动态变化的方法，可以尝试发现地震、火灾和砍伐等极端自然现象或人类活动等其它特异变化因子驱动的植被变化。研究汶川地震及其次生灾害对卧龙国家级自然保护区大熊猫栖息地植被的影响，同时分析2003—2012年间卧龙国家级自然保护区植被变化的时空特性，以期为卧龙大熊猫栖息地的恢复与保护提供科学依据。

1研究区域

卧龙国家级自然保护区建立于1963年，是我国最早建立的综合性国家级保护区之一。1980年，卧龙自然保护区加入联合国教科文组织“人与生物圈”保护区网，并与世界野生生物基金会合作建立中国保护大熊猫研究中心，主要保护西南高山林区自然生态系统及大熊猫等珍稀动物。2006年被列入世界自然遗产名录。该保护区的地理位置特殊，位于四川盆地西缘，邛崃山脉东南坡(东经102°52′─103°25′，北纬30°45′─31°25′)，成都平原向青藏高原过渡的高山峡谷地带，地势由西北向东南急剧递减。海拔高度从5795米减至1195米(ASTER GDEM2全球数字高程数据)相对高差达4600 米。区内主要有正河、西河和中河等河流，河谷呈“V”形状，落差较大。卧龙国家级自然保护区气候属典型的亚热带内陆山地气候，夏季温暖湿润，冬季寒冷干燥。从山谷到山顶形成了亚热带(2000m以下)、温带(2000—2600m)、寒温带(2600—3600m)、寒带(3600—4400m)、高寒带(4400—5000m)和极高山冰冻雪带(5000m以上)等气候垂直带。保护区内沙湾气象站(海拔1920m)年平均气温8.5℃，1月平均气温-0.9℃，7月平均气温17.1℃，年日照时数950 h，年降雨量890mm，相对湿度80.3%以上。

大熊猫对生境的植被群落有较广的适宜性。在卧龙，落叶阔叶林、落叶-针叶混交林及针叶林3种植被类型均为它的适宜生境[22- 23]。落叶阔叶林主要分布在海拔2000—2600m之间，卧龙的落叶阔叶林是常绿阔叶林或常绿-落叶阔叶林；针阔混交林主要分布在海拔2200—2700m之间；亚高山针叶林主要分布在2500—3200m之间，林下分布大熊猫可食竹类，是卧龙保护区大熊猫最重要的生境。

2基于小波变换的植被动态变化检测方法

2.1数据及处理2.1.1MODIS 数据

在连续或者大范围的数据分析中，中、高分辨率的遥感影像如TM，SPOT和雷达数据并不合适，因为这类传感器会生成大量的影像数据，进行变化检测分析所需计算时间长且成本较高，实现变化检测的算法上也会有相当的难度[24]。而相对来说，MODIS数据以其具有高时间分辨率的重要特征，在全球及大尺度区域植被覆盖变化研究中显示出其它数据所无法替代的作用。本研究选择2003年1月至2012年12月的MODIS植被指数产品(MOD13A1)作为主要数据源。MOD13A1数据是MODIS的3级产品，内容为归一化植被指数和增强型植被指数(NDVI/EVI)，该产品在1B数据的基础上，对由遥感器成像过程产生的边缘畸变进行校正，提供16d合成，每年23期的500m空间分辨率的3级正弦投影产品。本研究的研究区域内共包含9485个MODIS网格。

2.1.2DEM 数据

DEM数据选用ASTER GDEM2全球数字高程模型数据。ASTER GDEM2由美国国家宇航局(NASA)和日本经济产业省(METI)联合研制，其垂直精度为10—25m。数据经过云掩膜处理和残差异常处理，最终形成1°×1°的数据块。ASTER GDEM2数据以GEOTIFF格式分发，坐标系统为WGS84，数据空间分辨率约为30m。在进行空间统计时，由于ASTER GDEM2数据与MODIS NDVI数据的空间分辨率不同，利用ARC/GIS软件自动选取较低分辨率进行重采样，统计结果的空间分辨率与MODIS NDVI数据一致。

2.2基于小波变换检测植被动态变化

小波变换的方法是采用正交或复正交变换，应用滤波器对时间序列进行分析的一种新兴技术,其优于傅里叶分析之处就在于具有良好的局部化性质[6]。小波变换的方法可以用来分析植被年际变化的周期特征以及自然与人类对植被的极端影响，譬如：火灾、地震以及地震次生灾害带来的影响。该方法已经被成功地应用在众多领域，包括医学、地球物理和遥感图像处理等等[25]。与其它提取局部频谱信息的分析工具如窗口傅里叶变换相比，小波变换可以分析主频域比较广的数据，因为这类数据不适合预先定义尺度参数[26]。小波变换分为离散小波变换和连续小波变换，连续小波变换的连续性使变换结果与原始信号相比较时更加直接和明显[25]。

(1)

(2)

当原始信号为离散时间序列xn时，连续小波变换可以利用以下公式计算：

(3)

式中,s和n分别表示尺度参数和位移参数。公式(3)的结果是不同尺度下的小波系数集合。假定尺度和位移参数为s和n，则通过公式(3)可以得到所需尺度及位置上的小波系数，该系数能更好地反映原时间序列的周期性波动(图1)，过滤掉部分高频噪声信号，为下一步的植被变化检测奠定基础。

图1　NDVI时间序列在设定不同尺度参数进行小波变换后的小波系数曲线Fig.1　Original NDVI curve and coefficients curve after continuous wavelet transform with different scale factors尺度参数越小，变换后的小波系数能保留越多数据波动的细节信息，反之能更好地体现时间序列数据的整体波动趋势

本研究中，通过对2003—2012年共230期NDVI数据进行连续小波变换，选择合适的尺度参数，并结合时间序列异常点检测算法分析变换后的小波系数曲线，进一步分析统计数据发现保护区内植被指数的时空变化特性。

2.2.1连续小波变换

首先将样点的NDVI数据提取出来作为下一阶段的原始输入信号。在选择小波母函数时主要考虑以下几个因素：正交或非正交，实小波或复小波以及小波的波形。正交小波可以给出原始信号最紧凑的表示方式，但非周期的偏移会导致产生不同的小波频谱，而非正交小波光滑和连续变化的特性使其在计算上会产生冗余，但这在长时间序列的分析上是十分重要的。另外，基于复小波的变换能够保留振幅及相位信息，与实小波相比在捕捉变化时是更好的选择[26]。在常用的非正交复小波中本文选择了Morlet小波(图2)，其振幅光滑且连续，而且波形更能反映NDVI时间序列的波动特征[25]。利用小波变换的方法在处理时间序列时，在数据的两端可能会产生所谓的“边界效应”，为了避免边界数据变换结果发生严重偏差，本文采用了“对称性延拓”的方法，即将数据进行对称拓展，长度为n的数据拓展为长度为2n的时间序列。

尺度参数越小越能反映原始信号的细节信息，而反之能更好地提取波动的整体趋势(图1)。另外尺度参数越大，在信号的低频率部分具有高的频率分辨率，但在时域上时间分辨率较低；相反，尺度参数越小，在信号的高频率部分具有低的频率分辨率，但在时域上具有高的时间分辨率。因此选择合适的尺度参数尤为重要。

在连续小波变换中，尺度参数与频率的关系如下面的公式所示：

(4)

另外，在使用显微镜时一定要对镜头与玻片之间的观察距离有一个大致的概念，做到心中有数，不可盲目调焦。要根据物镜放大倍数与工作距离的关系进行正确的调节。这样既可以提高调焦的速度和准确性，又可避免物镜与标本片碰撞，污染和损坏镜头和标本，造成不必要的损失。

本文主要研究植被的年间动态变化，因此所适合的尺度参数需要使得到的小波系数曲线更能反映年度间的数据波动趋势。因此，将接近原信号采样频率1/23≈0.0435的尺度参数(16—19)进行比较(图3)。从图中可以发现，在这些尺度参数之间，当尺度参数为16时得到的小波系数曲线更为光滑，在时域上有较好的时间分辨率，当尺度变大时，系数曲线波动更频繁，无法更好地表现原信号的整体变化趋势。因此本文中选择尺度参数为16进行连续小波变换。

图2　Morlet小波函数Fig.2　Morlet Wavelet Funtion[-4, 4]区间中1000个点的Morlet小波相应的数值

尺度参数ScalefactorS伪频率Pseudo-frequencyFa尺度参数ScalefactorS伪频率Pseudo-frequencyFa10.813110.073920.406120.067730.271130.062540.203140.05850.163150.054260.135160.050870.116170.047880.102180.045190.0903190.0428100.0813200.0406

图3　尺度参数为16到19的小波变换结果Fig.3　Results from wavelet transform with scale factors from 16 to 20

2.2.2系数曲线分析

通过连续小波变换得到研究区域的相应小波系数，小波系数能够更加清楚的反映变化的趋势。而通过计算小波系数的峰值方差能够获得植被变化发生的幅度及时间点，系数中峰值方差越大说明该地区发生变化的概率越大，而峰值下降最大的时间即为变化发生的时间点。

在以上分析方法的基础上，由于植被指数本身会受到气候等因素的影响呈现季节性的变化特性，因此在小波变换方法的基础上结合季节和趋势叠加模型以及迭代检测结构异常点的方法[4]来进一步检验研究区域的植被变化特征。首先定义趋势-季节模型：

Yt=Tt+St+et,t=1,...,n

(5)

式中，Yt时间为t时的观察数据，Tt为趋势部分，St为季节部分，et为剩余部分(包括受云或传感器自身影响产生的噪声等)。

(6)

(7)

式中，K为谐波级数，公式(7)中f为采样频率。为了更容易拟合季节部分，公式(8)将公式(7)转为多个线性谐波回归模型。

(8)

迭代的检测异常点算法如下：

上述4个步骤迭代进行，直到检测出的突变点不再改变。

利用该算法对卧龙国家级自然保护区的NDVI数据进行突变点检测，统计分析突变点发生在汶川地震之前及之后的区域面积，从而反映地震对卧龙国家级自然保护区内植被的破坏程度。

3结果分析

3.1地震对卧龙国家级自然保护区内植被的影响

图4　卧龙国家级自然保护区植被指数极值减少(a)或增加(b)时间点发生在地震前、后Fig.4　NDVI′s extreme values greatly decreased (a) or increased (b) within Wolong National Natual Reserve before and after the earthquake

连续小波变换后得到的小波系数曲线能够更为直观地反映植被指数的波动趋势，因此通过计算系数曲线中的极值变化进一步分析植被指数曲线中极值的变化趋势。小波变换的分析结果及空间统计结果表明，地震后的植被指数极值减少面积大范围增加(表2，图4)，且植被指数极值波动较震前更为明显(表3，表4，图5)。由于2008年5月12日汶川地震是大陆内部的浅源地震，波及范围较广，震后卧龙自然保护区出现大面积山体滑坡，地表下陷，导致熊猫栖息地植被遭到严重的损毁。

图5　卧龙国家级自然保护区植被指数极值不同程度减少(a)或增加(b)的范围Fig.5　The area of NDVI′s extreme values decreased (a) or increased (b) at different levels in Wolong National Natural Reserve

时间Time减少Decreased/km2增加Increased/km2震前Before379.5165567.3429震后After1656.52101468.6946

表3　卧龙国家级自然保护区内植被指数极值发生不同等级减少的时间点发生在地震前、后的面积

表4　卧龙国家级自然保护区植被指数极值发生不同等级增加的时间点发生在地震前、后的面积

3.2卧龙国家级自然保护区植被的时间动态变化

选取地震前(2008年5月2日)与地震后(2008年6月3日)两期MODIS植被指数产品，计算地震前、后卧龙国家级自然保护区植被指数大于0.4的斑块数；斑块平均周长、总周长；斑块平均面积、总面积。计算结果见表5。虽然随着季节进入初夏，植被生长使得植被指数大于0.4的范围略有增加，但地震前后植被指数大于0.4的斑块数从4个增至7个，斑块平均面积缩小近一半。斑块数面积比从0.0033增加至0.0055。说明：地震引起的地表变化使得卧龙国家级自然保护区大熊猫生境景观生态格局发生了变化，大熊猫主要栖息地地更加破碎，生态系统稳定性受到了严重干扰。同时，从另一角度印证了通过小波变换方法提取变化信息的结果，统计数据更为直观地反映了地震以及地震次生灾害等极端现象对该地区植被的破坏程度的可靠性。

通过按照季节时间研究小波变换后的系数曲线中极值的变化，结果表明：极值的增加与减少范围都主要发生在夏季，而春季和秋季是植被指数极值减少和增加的第二重要季节 (表6，表7)。通过对植被指数极值变化按照年份进行分析研究，结果表明：在地震发生的2008年，植被指数极值减少的面积在地震前后的10a当中，并不是范围最大的，只排在第三位，位于震前2007年、震后2012年之后(表8)，经过数据质量检查发现：2007年第33天与2012年第273天卧龙保护区NDVI数据受云影响，质量存在问题，使得2007、2012年植被指数极值减少面积发生了偏差。

表5卧龙国家级自然保护区地震前、后植被指数大于0.4的斑块数、斑块面积

Table 5The number and area of patches where NDVI values were higher than 0.4 before and after the earthquake within Wolong National Natural Reserve

统计数据StatisticalData地震前Before(2008-05-02)地震后After(2008-06-03)平均周长AveragePerimeter/km173.510294.8238平均面积AverageArea/km2306.6797160.4169总面积TotalArea/km21226.71861283.3352斑块数NumberofPatches47斑块数/面积NumberofPatches/Area0.00330.0055

表6　卧龙国家级自然保护区植被指数极值增加发生在不同海拔高度各季节的面积/km2

为了更加全面的分析小波变换之后的结果，结合上文所提到的趋势-季节模型以及时间序列分析算法进一步研究。地震及其次生地质灾害破坏了大面积保护区内的植被，因此我们选取检测出的突变点中对应系数曲线数值最小值的时间点即NDVI数据发生了急剧下降情况的时间。由于原始NDVI 数据在2007年和2012

表7　卧龙国家级自然保护区植被指数极值减少发生在不同高度海拔各季节的面积/km2

表8　卧龙国家级自然保护区植被指数极值减少、增加发生在各年份的面积/km2

图6　利用时间序列突变点检测算法分析小波系数曲线研究卧龙国家级自然保护区植被在震前和震后发生急剧退化情况的区域Fig.6　The regions in Wolong National Natural Reservewhere vegetation had a obvious degeneration before or after the earthquake using the break points detection approach from time series data

年受云量影响，质量较差，在检测算法中也发现大量突变发生在该时间段。因此在算法执行过程中对突变值发生在2007年或2012年时的像元进行调整，选取该像元其他突变点中对应最小值的时刻。图6表明震后卧龙国家级自然保护区内的植被发生退化的区域明显大于震前。表9中为突变点发生在各年份的范围，表明2008年地震发生时，植被指数急剧下降的面积增幅明显(比2006年增加25%以上)。地震后卧龙国家级自然保护区内山坡上的岩石受到地震波的影响破碎，石块堆积和浅薄植被大量脱落，进而由于降水和重力的影响下引发泥石流等地质灾害[31]，对卧龙自然保护区内的植被造成严重破坏。因此，表9中突变点发生的范围在2008年后呈逐渐增长的趋势，卧龙国家级自然保护区内的植被在震后的恢复情况并不理想。

3.3卧龙国家级自然保护区植被的空间动态变化

植被变化在三维空间上(包括经度、纬度和高程)存在较大的空间异质性。在卧龙国家级自然保护区这样不足一个经、纬度的区域内，植被垂直高程的空间变化远大于水平方向的变化。因此选用ASTER GDEM2全球数字高程模型数据与通过小波变换后的植被指数极值变化结果进行空间叠加分析统计，从而得出卧龙保护区植被指数极值变化的空间垂直分布。表5，表6，图7是卧龙保护区植被指数减少与增加极值大小发生在不同垂直高度的面积。根据周世强、陈立顶等人的研究[21,32], 亚高山针叶林下分布大熊猫可食竹类，是卧龙保护区大熊猫最重要的生境。加之与植被指数统计结果比较发现：大熊猫最适宜生存区域植被指数一般大于0.4。在2600—2800m高度段，植被指数极值减少量大于0.4的范围,占总面积的53%,大于植被指数极值增加量大于0.4的范围,占总面积的34%(图7)。地形垂直高度变化是大熊猫生存活动的最重要影响因子，随着垂直高度的增加，在极大程度上限制了大熊猫的活动范围和觅食能力。据野外调查发现，大熊猫通常在1400—3600m的范围内活动。不同海拔高度范围内出现的频率明显不同，大熊猫出现频率较高的地段，是大熊猫最适宜的高度，2600—2800m熊猫出现的频率最高，因此，可以认为这个高度段最适合于大熊猫的活动和生存[32]。而结果表明在这个最适宜大熊猫生存的区域较高覆盖率的植被减少大于植被增加(表6，表7，表10)，植被恢复状况并未达到最理想的水平。植被指数减少的面积在春夏两季较大，表明在该时间段卧龙地区大熊猫最适宜生存区域的植被情况较为不稳定，需更为关注其动态，采取适当的保护措施。

表9　卧龙国家级自然保护区植被指数时间序列突变点发生在各年份的面积/km2

表10　卧龙国家级自然保护区植被指数极值减少发生在不同海拔高度的面积/km2

图7　卧龙国家级自然保护区植被指数极值随垂直高度变化Fig.7　The extreme values changed at different elevations in Wolong National Natural Reserve

4结果与讨论

本研究主要是以MODIS的2003—2012年植被指数NDVI产品和ASTER GDEM2数字高程数据作为输入数据，针对小波变换方法在大范围区域空间上进行植被动态变化检测的应用不足[33- 34]，利用该方法对卧龙国家级自然保护区植被变化信息进行研究。Lhermitte等[35]认为与其它将光谱转换到频域的方法(如：傅里叶变换)相比，小波变换分析更具优势，因为该方法在分解光谱数据方面更能反映时间序列数据的局部变化[36]，且可以同时在时域和频域上对数据进行分析而受到不稳定数据的影响较小[3]。小波变换在分析长时间序列的变化特性时，可以同时发现时间序列中的波动模式以及这些波动的时间特性[26]。通过选取适合的小波母函数并计算相应的尺度参数能够对时间序列进行整体趋势和局部变化的分析[25- 26]，而且经过变换后的小波系数曲线能够更为直观地体现原始数据的波动情况，便于进一步提取不同时间尺度下的变化信息。研究首先将10a，230期时间序列的植被指数NDVI通过小波变换转换到频率域，再分析植被在频率域表现出的不同变化特征，得到变化极值以及极值发生的时间点。统计结果表明地震对保护区内的植被造成大面积破坏，虽然震后局部区域有植被恢复的趋势, 但恢复幅度小于破坏幅度，而且在震后三年内地震引发的次级地质灾害频发，植被的恢复情况并不理想。为了进一步分析植被指数极值的季节和年际变化，结合趋势-季节模型和时间序列异常点检测算法计算小波系数曲线中变化异常点的时间分布，发现该地区植被指数极值变化主要发生在夏季。夏季通常是植被覆盖达到最高值的时间，同时也是植被变化最不稳定的季节。卧龙国家级自然保护区是中国大熊猫的重要栖息地之一，栖息地的植被对大熊猫的生存有着重大影响，因此本文针对保护区内大熊猫最适宜栖息的海拔区域范围分析其植被在地震前后的变化情况。结果表明在该海拔高度段(2600—2800m)，植被指数极值减少量大于0.4的范围大于植被指数极值增加量大于0.4的范围，表明植被恢复状况并未达到理想的水平，需要更为关注并采取相关保护措施。尽管本文中的方法能够提取地震前后植被的主要变化信息，并定位异常变化时间点，但结果中发现2008年并不是该保护区的植被指数极值变化范围最大的年份，震前的2007年及震后的2012年植被指数变化都大于地震发生的2008年。经检查发现NDVI数据在这两年受云影响，质量存在问题，使得统计植被指数极值减少面积发生了偏差。因此，数据质量对植被指数极值也存在非常大的影响。

本文中基于小波变换从长时间序列、大范围遥感数据中快速、自动化检测植被动态变化的方法对于卧龙国家级自然保护区生态环境变化，尤其是检测地震、火灾和砍伐等极端自然现象或人类活动等其它特异变化因子驱动的植被变化，具有明确的指示意义。进一步的研究工作将聚焦本研究的不足，集中在地震、火灾、砍伐等极端自然现象或人类活动等其它特异变化因子驱动的植被极端影响变化与气候影响变化的区分等工作上。同时在更大区域范围、利用更高空间分辨率的植被指数数据进行尝试，实现植被动态变化检测方法的大数据密集型计算也是一项重要的工作。

参考文献(References):

[1]Anyamba A, Eastman J R. Interannual variability of NDVI over Africa and its relation to El Nio/Southern Oscillation. International Journal of Remote Sensing, 1996, 17(13): 2533- 2548.

[2]Justice C O, Townshend J R G, Holben B N, Tucker C J. Analysis of the phenology of global vegetation using meteorological satellite data. International Journal of Remote Sensing, 1985, 6(8): 1271- 1318.

[3]Sakamoto T, Yokozawa M, Toritani H, Shibayama M, Ishitsuka N, Ohno H. A crop phenology detection method using time-series MODIS data. Remote Sensing of Environment, 2005, 96(3/4): 366- 374.

[4]Verbesselt J, Hyndman R, Zeileis A, Culvenor D. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sensing of Environment, 2010, 114(12): 2970- 2980.

[5]Hall-Beyer M. Comparison of single-year and multiyear NDVI time series principal components in cold temperate biomes. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(11): 2568- 2574.

[6]Li Z T, Kafatos M. Interannual variability of vegetation in the United States and its relation to El Nio/Southern Oscillation. Remote Sensing of Environment, 2000, 71(3): 239- 247.

[7]Sarkar S, Kafatos M. Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters. Remote Sensing of Environment, 2004, 90(2): 268- 280.

[8]Azzali S, Menenti M. Mapping vegetation-soil-climate complexes in southern Africa using temporal Fourier analysis of NOAA-AVHRR NDVI data. International Journal of Remote Sensing, 2000, 21(5): 973- 996.

[9]Lim C, Kafatos M. Frequency analysis of natural vegetation distribution using NDVI/AVHRR data from 1981 to 2000 for North America: Correlations with SOI. International Journal of Remote Sensing, 2002, 23(17): 3347- 3383.

[10]Lunetta R S, Knight J F, Ediriwickrema J, Lyon J G, Worthy L D. Land-cover change detection using multi-temporal MODIS NDVI data. Remote Sensing of Environment, 2006, 105(2): 142- 154.

[11]Gross D, Dubois G, Pekel J F, Mayaux P, Holmgren M, Prins H H T, Rondinini C, Boitani L. Monitoring land cover changes in African protected areas in the 21st century. Ecological Informatics, 2013, 14: 31- 37.

[12]Mao D H, Wang Z M, Luo L, Ren C Y. Integrating AVHRR and MODIS data to monitor NDVI changes and their relationships with climatic parameters in Northeast China. International Journal of Applied Earth Observation and Geoinformation, 2012, 18: 528- 536.

[13]Truchetet F, Laligant O. Wavelets in industrial applications: a review. Proceedings of the SPIE Proceedings Vol. 5607, Wavelet Applications in Industrial Processing II, 1, 2004.

[14]Goupillaud P, Grossmann A, Morlet J. Cycle-octave and related transforms in seismic signal analysis. Geoexploration, 1984, 23(1): 85- 102.

[15]Percival D B, Wang M, Overland J E. An introduction to wavelet analysis with applications to vegetation time series. Community Ecology, 2004, 5(1): 19- 30.

[16]Martínez B, Gilabert M A. Vegetation dynamics from NDVI time series analysis using the wavelet transform. Remote Sensing of Environment, 2009, 113(9): 1823- 1842.

[17]余明艳, 高素萍. 震后卧龙——蜂桶寨生态廊道受灾状况及灾后恢复对策与建议 // 第四届中国竹业学术大会论文集. 雅安: 中国林学会, 四川农业大学, 2008: 58- 63.

[18]王岑涅, 高素萍, 孙雪. 震后卧龙-蜂桶寨生态廊道 大熊猫主食竹选择与配置规划. 世界竹藤通讯, 2009, 7(1): 11- 16.

[19]鄢武先, 张小平, 邓东周, 魏宗华, 尤继勇, 杨晓军, 徐海斌. 卧龙自然保护区地震灾后植被恢复主要植物种选择研究. 四川林业科技, 2012, 33(3): 32- 36.

[20]周世强, 古晓东, 黄金燕. 不同山系大熊猫栖息地植物区系特征的定量分析. 四川林业科技, 2006, 27(1): 14- 18.

[21]周世强, 冯莉, 张亚辉, 张和民. 野生大熊猫地理分布格局的空间尺度分析. 四川林业科技, 2009, 30(5): 53- 57.

[22]欧阳志云, 李振新, 刘建国, 安力, 张和民, 谭迎春, 周世强. 卧龙自然保护区大熊猫生境恢复过程研究. 生态学报, 2002, 22(11): 1840- 1849.

[23]周世强, 黄金燕, 张亚辉, 李德生, 黄炎, 周小平, 王鹏彦, 张和民. 卧龙自然保护区大熊猫栖息地植物群落多样性V:不同竹林的物种多样性. 应用与环境生物学报, 2009, 15(2): 180- 187.

[24]Lu D, Mausel P, Brondízio E, Moran E. Change detection techniques. International Journal of Remote Sensing, 2004, 25(12): 2365- 2401.

[25]Cheng T, Rivard B, Sánchez-Azofeifa G A, Feng J, Calvo-Polanco M. Continuous wavelet analysis for the detection of green attack damage due to mountain pine beetle infestation. Remote Sensing of Environment, 2010, 114(4): 899- 910.

[26]Torrence C, Compo G P. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 1998, 79(1): 61- 78.

[27]Kumar P, Foufoula-Georgiou E. Wavelet analysis for geophysical applications. Reviews of Geophysics, 1997, 35(4): 385- 412.

[28]Cleveland R B, Cleveland W S, Mcrae J E, Terpenning I. STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 1990, 6(1): 3- 73.

[29]Zeileis A. A unified approach to structural change tests based on ML scores, F statistics, and OLS residuals. Econometric Reviews, 2005, 24(4): 445- 466.

[30]Venables W N, Ripley B D. Modern Applied Statistics with S. New York: Springer, 2002.

[31]刘明冲. 浅析卧龙自然保护区地质灾害的危害. 四川林勘设计, 2011, (4): 46- 48.

[32]陈利顶, 刘雪华, 傅伯杰. 卧龙自然保护区大熊猫生境破碎化研究. 生态学报, 1999, 19(3): 291- 297.

[33]Cai D L, Guan Y N, Guo S, Yan B P, Xing Z, Zhang C Y, Piao Y P, An X D, Kang L H. Rapid assessment of large scale vegetation change based on multi-temporal phenological analysis // Proceedings of the SPIE Proceedings Vol. 8002: MIPPR 2011: Multispectral Image Acquisition, Processing, and Analysis. Guilin, China: SPIE, 2011.

[34]Piao Y C, Yan B P, Guo S, Guan Y N, Li J H, Cai D L. Change detection of MODIS time series using a wavelet transform // 2012 International Conference on Systems and Informatics (ICSAI). Yantai: IEEE, 2012: 2093- 2097.

[35]Lhermitte S, Verbesselt J, Jonckheere I, Nackaerts K, Aardt J A N V, Verstraeten W W, Coppin P. Hierarchical image segmentation based on similarity of NDVI time series. Remote Sensing of Environment, 2008, 112(2): 506- 521.

[36]李弼程, 罗建书. 小波分析及其应用. 北京: 电子工业出版社, 2003.

transform method in Wolong Natural Reserve

PIAO Yingchao1,2,3, GUAN Yanning1,*， ZHANG Chunyan1, GUO Shan1, YAN Baoping2

1ChineseAcademyofSciences,InstituteofRemoteSensingandDigitalEarth,Beijing100101,China2ChineseAcademyofSciences,ComputerNetworkInformationCenter,Beijing100190,China3UniversityofChineseAcademyofSciences,Beijing100049,China

Abstract：On May 12, 2008, the Wenchuan earthquake struck the Sichuan province, triggering strong continuous aftershocks in its surrounding regions, which led to secondary geological disasters, such as landslides, collapses, and mudslides, in more than 120 regions in Wolong National Natural Reserve. This bioreserve consists of broad-leaved deciduous forests, deciduous coniferous mixed forest, and coniferous forests, and is considered an important habitat of giant pandas in China. The population of pandas in this area accounts for approximately 10% of the total number of pandas in the entire country. However, vegetation in this bioreserve was adversely affected by the 2008 earthquake. With the development of earth observation techniques over the last decade, long-term observations of natural vegetation have detected large-scale land-cover changes. Among the various long-term satellite data, MODIS data have significant potential in detecting vegetation changes from regional to global scale because of their high temporal resolution. In this research, we proposed a method based on wavelet transform to extract information regarding vegetation changes using normalized difference vegetation index (NDVI) time series data in Wolong National Natural Reserve. In particular, MODIS 16-daily products recorded from 2003 to 2012 with 500-m resolution were selected for investigation, and the effects of the earthquake and secondary geological disasters on vegetation in Wolong National Natural Reserve were quantified. We obtained numerous relevant findings. (1) The proportion of the areas where the NDVI greatly decreased after the earthquake is significantly larger than before. Moreover, by using the trend-seasonal model and breakpoint detection method from the time series data, we discovered that the NDVI curves fluctuated more obviously after the earthquake. (2) By performing a seasonal analysis on the extreme values of NDVI time series, we determined that the increased extreme values mostly occurred in summer because the vegetation cover during this season was the highest and also the most unstable one in a year. The second most important season was autumn. Compared with the increased values, the decreased NDVI extreme values often emerged in summer and spring. The results specified that the area with decreased NDVI data was the third highest when the Wenchuan earthquake transpired. Nonetheless, the statistical results obtained from comparing the inter-annual difference were biased because of the identified moderate cloud effect on NDVI data in 2007 and 2012. (3) The Wenchuan earthquake that occurred on May 12, 2008 was an intra-continental shallow-focus earthquake, which extensively affected the environment and areas surrounding Sichuan province. In this research, the data were also analyzed at land elevations where giant pandas are most comfortable in (2600—2800 m) to quantify the influence of the earthquake and secondary geological disasters on the habitat of the giant pandas. Accordingly, we realized that the extent of decreased NDVI extreme values remarkably exceeded the instances of increased NDVI. Vegetation recovery in Wolong National Natural Reserve did not achieve an expected level. Based on the above findings, the method introduced in this paper can be used to detect significant changes in long-time NDVI data; however, it cannot identify whether the causative factor is a disastrous natural phenomenon (e.g., earthquake, fire, etc.) or a human activity (e.g., deforestation, urbanization, etc.). In future, we intend to apply our method in a wider scope when data-intensive calculation technique is required.

Key Words:vegetation change; wavelet transform; time series analysis; Wolong National Natural Reserve; Wenchuan earthquake

基金项目:国家自然科学基金国际(地区)合作与交流项目中美软件合作研究项目(61361126011); 中国科学院信息化建设项目(XXH12504-1-12); 中国科学院遥感与数字地球研究所135计划项目突破二第二课题(Y3SG0500CX); 中国科学院信息化专项重点数据库项目(XXH12504-1-06); 中国科学院计算机网络信息中心一三五规划重点培育方向专项(CNIC_PY_1408); 中国科学院计算机网络信息中心一三五规划重点培育方向专项(CNIC_PY_1409)

收稿日期:2014- 07- 24; 网络出版日期:2015- 08- 26

*通讯作者

Corresponding author.E-mail: guan@irsa.ac.cn

DOI:10.5846/stxb201407241504

朴英超，关燕宁，张春燕，郭杉，阎保平.基于小波变换的卧龙国家级自然保护区植被时空变化分析.生态学报,2016,36(9):2656- 2668.

Piao Y C, Guan Y N, Zhang C Y, Guo S, Yan B P.Analysis of temporal and spatial changes in vegetation cover using wavelet transform method in Wolong Natural Reserve.Acta Ecologica Sinica,2016,36(9):2656- 2668.