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

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

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

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

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

卧龙国家级自然保护区建立于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之间，林下分布大熊猫可食竹类，是卧龙保护区大熊猫最重要的生境。

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

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

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

小波变换的过程中，通过小波母函数{ψ(t)}的尺度和位移变化产生一组小波基函数{ψ_λ,τ(t)}^[16]。

点击浏览公式

式中，λ和τ分别表示小波的尺度参数和位移参数。这两个参数分别决定了小波基函数的宽度和位置，而且小波母函数必须满足Foufoula-Georgiou和 Kumar曾提出的相关条件^[27]。S(t)是在时域上的实值信号函数，关于ψ_λ,τ(t)的连续小波变换如下式所示：

点击浏览公式

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

点击浏览公式

式中,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数据进行连续小波变换，选择合适的尺度参数，并结合时间序列异常点检测算法分析变换后的小波系数曲线，进一步分析统计数据发现保护区内植被指数的时空变化特性。

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

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

图选项

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

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

点击浏览公式

式中,F_a为尺度参数为S时的伪频率，f为采样周期，F_c为所选小波母函数的中心频率。在不同尺度下通过连续小波变换得到的小波系数越大，说明这个位置上原信号含有F_a附近的频率成分比较多，从而能够反映在该位置原信号与ψ_λ(t)的相似度，小波系数越大则相似度越高。根据所选的小波母函数Morlet小波，计算尺度参数与频率的关系如表 1所示。

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

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

图选项

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

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

点击浏览公式

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

假设T_t是分段线性且有m个断点t₁^*,...,t_m^*，定义t₀^*=0，则趋势部分表示如下：

点击浏览公式

与构建趋势部分相似，假设季节部分有p个断点τ₁^#,...,τ_p^#，定义τ₀^#=0，τ_p+1^#=n,则季节部分表示如下：

点击浏览公式

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

点击浏览公式

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

利用基于局部多项式回归拟合(LOESS)的时间序列季节分解算法(STL)估计季节部分的初始值(所有季节子序列的均值)。^[28]

1)利用最小二乘基于迭代残差的移动和方法(OLS-MOSUM)^[29]检测趋势部分(Y_t-)是否存在突变点，如果存在突变点则表示为t₁^*,...,t_m^*；

2)基于M估计(M-estimation)^[30]计算公式(6)中趋势部分的系数α_i和β_i。设置趋势部分为,t_i-1^* < t≤t_i^*；

3)利用OLS-MOSUM检测季节部分(Y_t-)是否存在突变点，如果存在突变点则表示为τ₁^#,...,τ_p^#；

4)基于M-estimation计算公式(8)中季节部分的系数γ_j,k和θ_j,k。设置季节部分为。

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

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

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

图 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

图选项

图 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

图选项

选取地震前(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年植被指数极值减少面积发生了偏差。

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

图 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

图选项

植被变化在三维空间上(包括经度、纬度和高程)存在较大的空间异质性。在卧龙国家级自然保护区这样不足一个经、纬度的区域内，植被垂直高程的空间变化远大于水平方向的变化。因此选用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)，植被恢复状况并未达到最理想的水平。植被指数减少的面积在春夏两季较大，表明在该时间段卧龙地区大熊猫最适宜生存区域的植被情况较为不稳定，需更为关注其动态，采取适当的保护措施。

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

图选项

本研究主要是以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数据在这两年受云影响，质量存在问题，使得统计植被指数极值减少面积发生了偏差。因此，数据质量对植被指数极值也存在非常大的影响。

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