森林的健康状况直接关系到全球生态安全和人类社会的可持续发展。随着自然资源过度开发，森林资源锐减，森林生态系统的服务功能大大减退，全球出现了一系列的生态危机，诸如生物多样性减少，水、旱、火灾、森林病虫害等等频发，导致了世界大多数国家的森林生态系统都在不同程度地退化。森林灾害监测为森林受损评价提供基础数据，为受害森林恢复的经营管理提供依据。

目前的森林状况监测多是基于现地实测的，效率较低，无法实现对大区域的快速连续监测。随着空间信息技术的不断进步，遥感技术在林业中的应用日益广泛，越来越多的研究者利用遥感手段定性或定量地获取森林参数(Anderson et al., 2006; Zhao et al., 2009; 李崇贵等，2005; 庞勇等，2008; 覃先林等，2005; 王海鹏等，2008)，以及进行森林健康监测(Vogelmann et al., 1989; 肖文发等，2001; 赵春燕等，2006; 朱振华等，2008)。MODIS中等分辨率成像光谱仪因其时间分辨率高(每1~2天观测地球表面1次)，可以连续观察地物的变化，因而多用于区域尺度的研究(Zhang et al., 2008; 方墨人等，2004; 宫攀等，2006; 王红说等，2009; 于信芳等，2006; 李月臣等，2005)，但在森林生长异常监测方面的应用尚不成熟(Eklundh et al., 2009)。

本文尝试采用高时间分辨率的MODIS数据对区域尺度的森林生长状况进行监测，力图找出生长季内的异常点，提高大范围森林生长异常监测的效率，以便森林决策者对异常情况做出快速判断并采取相应措施。

1 研究区概况及数据获取 1.1 研究区概况

本文选取黑龙江、吉林、辽宁(118° 30'—136°E，38°40'—54°N)作为主要的研究区域。东北林区地处我国大、小兴安岭和长白山脉，生态地位非常重要。土地总面积78.73万km²，人口约10 715.4万人，属大陆性季风气候，四季分明，雨热同期，日照丰富，积温较高，冬长夏暖，春秋季短。地形以山地、平原、河流为主。

东北三省的林业用地面积达3 466.46万hm²，占土地面积的43.92%，森林覆盖率达37.99%，远远高出全国森林覆盖率20.36%的水平; 该区还是我国最主要的天然林区，天然林面积2 392.63万hm²，占该区有林地面积的80.26%。主要森林类型为针叶林、阔叶林和针阔混交林，大兴安岭北部山地寒温带针叶林区，是我国唯一的大片落叶松林区。主要森林树种有兴安落叶松(Larix gmelinii)、红松(Pinus koraiensis)、樟子松(Pinus sylvestris var.mongolica)、白桦(Betula platyphylla)、黄檗(Phellodendron amurense)、水曲柳(Fraxinus mandschurica)、鱼鳞云杉(Picea jezoensis var.microsperma)、山杨(Populus davidiana)、蒙古栎(Quercus mongolica)、榆树(Ulmus pumila)等。

随着经济的发展，东北地区天然林大量被消耗，人工纯林面积不断增大，生物多样性降低，森林生态系统失衡，导致森林鼠害，落叶松毛虫、光肩星天牛、青杨虎天牛等病虫害大面积发生; 并且该地区的森林覆盖率很高，林下可燃物堆积厚，年降雨量小，导致5月、10月的火险等级很高; 加之北部寒潮时有入侵，天气多变，易发生雪灾、风倒等灾害。

1.2 数据获取 1.2.1 时间序列MODIS数据获取

本文拟通过遥感技术连续地观察林木生长状态，进而进行森林生长的异常监测，因此需获取时间序列的遥感影像作为数据源。考虑到研究区覆盖面积较大，且对时间分辨率要求较高，选择Terra-MODIS地表反射率数据(MOD09A1)作为数据源(其空间分辨率为500 m，时间分辨率为8天)，时间范围为2004—2006年。根据NASA对全球区域的划分，东北三省MODIS影像主要覆盖h25v03，h26v03，h26v04，h27v04四景，下载的原始影像为HDF-EOS(hierarchical data format-earth observing system)存储格式，投影方式为ISIN(integerized sinusoidal)投影。

在Modis Reprojection Tool批处理模式下，对四景MODIS影像进行镶嵌，转投影方式为Lambert _Conformal_Conic投影，利用东北三省行政区域对影像裁切，得到研究区MODIS影像图(图 1)。

1.2.2 时间序列植被指数计算

植被指数是利用遥感卫星探测数据的线性或非线性组合而形成的能反映绿色植被生长状况和分布的特征指数(陈述彭，1990)，而时间序列的植被指数能更准确地反映各个阶段的植被生长状况。归一化植被指数(NDVI)的应用最广泛且表现稳定，但由于它在植被高覆盖区有易饱和现象，而增强型植被指数(EVI)改进了NDVI的某些缺陷，特别是大气噪声、土壤背景、饱和度等问题(Huete et al., 2002; 王正兴等，2003)，能更好地反映高植被覆盖区的植被生长状况，因而本文选用年时间序列的EVI描述林地的生长状况。

在IDL中对时间序列的MODIS影像利用式(1)(Huete et al., 2002)计算EVI，生成研究区时间序列的EVI影像图(图 2)。

(1)

式中: ρ_blue，ρ_red和ρ_NIR分别为MODIS影像的蓝、红、近红外波段，L(=1)为冠层背景调节系数，C₁(=6)，C₂(=7.5)为气溶胶散射修正系数，G(=2.5)为增益因子。

2 研究方法 2.1 遥感数据预处理

尽管8天最大值合成的MODIS地表反射率数据在一定程度上去除了云和大气的影响，但仍然可能受太阳高度角、观测角度、气溶胶等一些随机因素的干扰，使得时间序列数据的变化呈不规则状态(于信芳等，2006)，通常表现为EVI曲线短期起伏很大，这是不符合植被正常生长状况的，因而需要对合成的EVI数据进一步滤波，去除噪声干扰。

本文选用Savitzky-Golay滤波函数平滑时间序列EVI(Savitzky et al., 1964)。Savitzky-Golay滤波是一种最小二乘法卷积平滑技术，其主要原理为:对长度为n的窗口邻域内各点数据，用一元p次多项式拟合，此多项式的系数借助最小二乘准则来确定，卷积核可以单独保存为系数矩阵，它仅由滤波器维数、窗口大小和阶数确定，多项式在窗口中心点的值即为经卷积平滑处理后的值。滤波过程需要经过多次试验，改变窗口大小和阶数等参数，既要达到平滑去噪的效果，又不可过度平滑弱化了植被原始信息。经比较选定窗口长度为7的一元二次多项式拟合，从平滑结果看曲线基本保持原有形状，但上下剧烈起伏的状况得到了改善(图 3)，认为Savitzky-Golay滤波很大程度上去除了噪声干扰影响，可以将平滑后的EVI值进行后续研究。

2.2 技术路线

生长受到胁迫的森林，其光合有效率、光谱特性、叶绿素含量、叶面积指数等往往会有所反应，生长季的EVI值也能够有所体现，其值一般会低于正常年份; 另外年时间序列的曲线上也必然起伏较大，复杂度较正常年份高。本文拟从这2个角度出发，计算EVI生长季面积(表征植被的生长状况)和EVI曲线的信息熵(表征曲线的复杂度)，联合判定异常点，确定灾害的发生。由于随机干扰的存在，单独从任一方面很难准确判定灾害的发生，但是二者的交集一定程度上可以剔除因偶然因素导致的异常。技术路线如图 4所示。

2.2.1 计算年EVI生长季面积

植被指数的生长季面积可以反映植被年内的生长状况。东北三省的常绿林所占比例较小，且其年内植被指数值变化不大，因而本文仅考虑主要的森林类型———落叶林和针阔混交林，计算生长季面积:

(2)

根据式(2)，首先确定积分限，即绿色开始出现的时间和绿色减到最少的时间，进而计算这一时间段EVI的积分。可见生长季开始和结束时间的算法决定着生长季面积计算的精度。本文采用Zhang等(2003; 2006)经过现地验证过的方法确定积分上、下限:每个独立的生长周期都包含绿色开始出现、绿色增加到最大值、绿色开始减少、绿色减少到最小值4个关键时刻点，通过S形函数(Logistic函数)分段拟合植被年生长曲线:

(3)

式中: EVI(t)为t时刻的EVI值; a，b为拟合参数; c+d为EVI的最大值; d为初始背景值。

计算该生长曲线的曲率及曲率的一阶导数，即曲率的变化率K' [(式(4)]，得到4个极值点，这4个极值点就是对应的植被生长的4个转折期，其中第1，4点分别对应积分的下限和上限。

(4)

式中: z=e^a+bt。

受灾害条件的影响，不同年份同一像元植被生长季的起止时间[式(2)积分上下限]及长度会有所不同，在计算生长季面积和曲线熵值时，时间范围如何选取是个关键的问题。如果假设其中一年是正常年，并将该年的生长季起止时间作为这3年统一的计算区间，则会受假设正常年的影响，一定程度上弱化了受干扰像元的信息，导致该像元不能被正常检测出来; 但如果同一像元取各自年内的生长季起止时间作为计算区间，可能会将由于全年气候影响导致的生长季起止时间发生变化的像元误检测，同时会降低年际间的可比性。因此，为了3年内同一像元间更具有可比性，本文生长季的起止时间选择同一像元3年的最晚发芽和最早落叶的时间。

2.2.2 计算年EVI曲线熵值

Shannon(1948)提出了信息熵的概念，信息熵是信息论中用于度量信息量的一个概念。一个系统越是有序，信息熵就越低; 不确定性越大，信息熵就越高。信息熵可以用来描述一个广义集合的内部状态的复杂程度，二者成正比例关系。根据香农信息论，本文中将信息熵概念推广到用于表达曲线的复杂程度，见式(5)。

(5)

式中: f(i)表示第i点的EVI值; p(i)表示f(i)出现的概率。熵值计算的起止时间选择与生长季相同的起止时间。

曲线的熵值越大，表明复杂度越高; 反之，熵值越小，复杂度越低。因此可以通过比较年EVI曲线的熵值来确定曲线的变异程度，从而判断是否发生异常。

2.2.3 生长季面积和熵值联合监测异常点

假设2004年某林地受到灾害影响发生异常，那么一般其后的1~2年内该处植被生长状况会比受灾之前差，对应像元的年时间序列EVI曲线的复杂度会高于正常年，而该年的生长季面积一般会小于正常年(受干扰程度不同，表现有所不同)，即:熵值越大，生长季面积越小，该像元发生异常的可能性越大，因而可以根据3年间生长季面积、熵值的变化来联合监测异常点。

分别计算3年时间序列EVI数据的熵值和生长季面积，并求其变化量，由于每年的气候条件不同，导致生长状况有差异，因而需通过设定阈值方法，剔除一些因气候等因素影响而异常的点或变化量相对较小的点; 然后取熵值和生长季面积变化量的交集，从而确定每年发生异常的像元; 画出3年异常点的EVI曲线，进行比较验证。

3 结果与分析 3.1 年生长季面积

8天合成的年时间序列EVI影像，每个像素都对应着46个值，根据式(2)和式(3)对这46个离散值用Logistic生长方程拟合年生长曲线，并计算该生长曲线的曲率变化率，得到的4个极值点分别对应一年当中绿色开始出现、绿色增加到最大值、绿色开始减少、绿色减少到最小值的时间点。经过计算, 2004, 2005, 2006年林地绿色开始出现的平均时间分别是119，125和119天，绿色减少到最小值的平均时间分别是286，298和301天。2005年绿色开始出现时间较另外2年稍晚，初步推断可能是2004年森林受到了干扰，导致次年展叶时间推迟，2006年干扰程度轻的植被又恢复了正常; 而一定程度的干扰可能会促进林下植被的更新，增强其活力，从而导致平均生长季略有延长，但不排除2005年由于气候等原因造成的绿色开始出现时间较晚，因而还不能仅凭关键时间点的变化判断是否发生异常。

然后对绿色开始出现至绿色减少到最小值的时间段积分，得到年EVI的生长季面积(图 5)。

从生长季面积图上可以明显地区分出林区及非林区，由于在生长季森林的植被指数值明显高于农田、河流以及建设用地等非林区，因此累积的生长季面积值大，图上表现为亮度值明显高于非林区。

东北三省境内有大、小兴安岭、长白山脉及千山山脉，林区纵贯南北。总体来看，小兴安岭至长白山区，森林类型是以红松为主的针阔混交林，其生长季面积值相对较高; 而黑龙江省北部大兴安岭地区是以兴安落叶松为主的针叶林区，生长季面积值相对较低; 辽东山区以次生林为主，森林类型多为油松(Pinus tabulaeformis)和落叶阔叶混交林，但由于该区域内的千山山脉是长白山支脉，部分植被与长白山脉又有相似之处，因而其生长季面积略低于小兴安岭林区但高于大兴安岭林区。

3.2 年EVI曲线熵值

根据式(4)逐年计算EVI曲线生长季的熵值并制图(图 6)。

曲线熵值是曲线复杂度的表达，森林植被季节性生长经历增加、生长盛期、凋落过程，因而森林植被的曲线复杂度要高于其他地物。但总体上看，林区与非林区生长季熵值不如生长季面积的差别明显，这主要是由于非林区地物EVI值的波动会增大曲线的熵值，但对总的EVI面积影响并不大。

比较3年每一点的熵值大小，可以直观地看到一些点在某一年熵值明显不同于其他2年，熵值越大说明EVI曲线越复杂，受灾害干扰的可能性也越大，但具体定量的分析还需要结合生长季面积来联合判定是否受到了干扰。

3.3 异常点提取

由于气温、降水等气候因素影响，正常生长的植被年际间也会有变化，只是其变化量会比受到干扰时的小，因而需设定相应阈值去除变化量较小的点，认为变化量大的点才是可能发生异常的点。

根据式(6)计算2004—2006年间生长季面积(area)和熵值(entropy)的变化量，参考该地区关于火灾(杜永胜等，2007)的历史记载(其他森林灾害事件的历史数据暂时缺乏)，进行阈值的设定，其原则是既保证过火面积的最大化提取，又保留其他的异常信息。经过反复试验，最终确定2004—2006年熵值增加范围[0.018 7，0.08]，面积减少范围[1000，2 900]。再取生长季熵值增加和面积减少部分的交集，即同时满足生长季EVI累积值减少且年EVI曲线复杂度高的点为最终提取的异常点，并制图(图 7)。

(6)

3.4 异常点曲线分析

如图 8，9所示为正常点与异常点3年EVI曲线的对比。可以看出正常点的EVI曲线整体波动不大，生长季长度、峰值差异并不显著，但异常点曲线可以明显看出年际间峰值有很大的波动，说明可以通过本文所述方法提取出受干扰的像元，并做进一步分析。

随机抽取异常点画出3年的EVI曲线(图 10)，以及对应生长季面积和熵值，比较其曲线形状、峰值，并验证本文方法所求熵值和生长季面积是否与曲线形状相符。

图 10中列举了2个异常点3年的EVI变化曲线，可以看出其中1年或2年明显出现了异常，2004年曲线形状符合正常的植被生长曲线: EVI逐渐增加，达到最大值后趋于平缓，持续一段时间后又开始下降; 2005年曲线形状基本正常，但峰值有所下降; 而2006年曲线明显异于前2年:峰值低，并且在生长旺季出现了双峰现象; 在2005年生长旺盛也出现了双峰现象，而2006年峰值变化更加明显，在生长阶段就出现了下降趋势，然后又缓慢上升达到峰值，但峰值阶段持续时间较短，再进行熵值和生长季面积的年际间比较: 2006年由于双峰的出现导致曲线复杂度高，曲线熵值明显大于前2年，另外峰值的降低以及较短的峰值持续时间，导致了生长季面积达到3年中最低值。可见，文中所述方法提取的生长季面积及曲线熵值可以反映曲线的趋势，作为异常点提取的指标是可行的。但具体的干扰因子尚不明确，如有的EVI出现了双峰，峰值部分趋势为上升-下降-上升，原因有可能是生长季早期发生了病虫害，生长季中后期逐渐恢复，EVI又呈上升趋势; 有的EVI达到最大值的时间较晚，峰值较正常年小，并且持续时间较短后就开始下降，这种情况分析可能是前一年发生了异常，导致下一年生长季增长缓慢等，因此具体的受干扰类型还需进一步的验证。

3.5 火点验证

为了明确发生异常的原因，以火点为例，选取2004年10月中旬黑龙江省黑河地区发生的林火灾害事件，对本文提取异常点进行验证。用于验证的火点数据来源于基于MODIS提取的火情监测指数NDTI(normalized difference thermal index)。NDTI的计算需要MODIS专门为火情监测设置的4 μm和11μm通道辐射数据(MOD021KM)以及定标数据(MOD03)。该指数经过现地验证(赵文化等，2008)，检测精度较高，并且计算速度要比MODIS推荐火情算法快2倍以上，可以快速实现火情提取。

图 11为截取局部的NDTI提取火点与本文提取异常点的比较。总体上看，2种方法提取的一致性相对较高，但由于本文提取的异常点中除火点外，还包含由于其他异常原因(如病虫害)导致的异常，因此异常点数稍多于NDTI的提取数，观察多提取点的EVI曲线，发现有些点与火点的曲线形状非常相近，但也有部分点相异较大。分析其原因: 1)由于NDTI的提取是基于单一时相，有些火点在当天未能检测但在时间序列中得以反映; 2)火场附近林地由于燃烧导致生境改变，生物多样性降低，从而易发生病虫害等灾害，但具体异常原因仍有待于结合其他数据辅助验证。

4 结论与讨论

本文通过对东北地区2004—2006年3年间的MODIS EVI数据的分析，结果表明时间序列的MODIS EVI数据可以用来进行森林生长异常制图。高时间分辨率的植被指数对森林生长状态敏感，甚至可以捕捉到细微的异常变化，一定程度上弥补了单一时相不能很好反映全年植被生长变化的不足。

针对时间序列数据的特点，文中将植被生长状况(生长季EVI面积)和时间序列EVI特征曲线(生长季熵值)的形态变化相结合进行异常分析，而不是局限于常规的基于分类的变化检测方法。由于这2个指标都与植被的生长状况密切相关，并且本文精确求取了生长季的时间范围，保证了指标提取的精度，因而一旦发生了灾害，指标上一定会有所反映，并结合适合的算法可以将异常点检测出来。

通过对东北三省林区的试验，本文提出的森林生长异常监测方法是基本可行的，克服了依赖现地调查的不足，特别适用于区域尺度，检测到的森林生长异常区域的EVI特征曲线也对异常的原因有一定的揭示; 观察异常点时间序列的EVI曲线发现:曲线的熵值明显大于正常年，生长季峰值低，导致生长季面积小于正常年，并且在生长旺季会出现峰值突然持续下降或出现双峰等异常现象。目前由于只对火点进行了初步的验证，虽然数据进行了滤波处理，但可能仍然会受到云的干扰导致误检测，因此下一步工作将利用MOD09A1数据的QA产品，去除有云像元，结合东北林区的森林灾害事件对本方法进行进一步的验证和改进，并分析曲线规律，找出各种异常的原因。