应用MODIS数据监测神农架林区大小蠹虫害时空分布

（东北林业大学 工程技术学院，黑龙江 哈尔滨 150040）

森林在全球生态系统中一直扮演着至关重要的角色，及时了解森林的健康状况对于维持生态系统的稳定具有极其重要的意义。森林虫害是导致森林健康质量降低的一个非常重要的因素，因此多年来对森林虫害的研究一直受到许多专家学者的密切关注[1-2]。

目前遥感影像监测是森林健康及资源状况调查的一个重要手段。MODIS 数据由于具有较高的时间分辨率，对于监测短期的、大面积的森林虫害的发生趋势应用广泛，并且MODIS 植被指数产品归一化植被指数（NDVI）是检测植被生长状态、植被覆盖度和消除部分辐射误差的重要依据。例如汪航等[3]选用2014—2017年上半年的MODIS 数据重构NDVI 时间序列数据集，并且利用3 种滤波方法对新疆巴楚胡杨林的春尺蠖虫害进行监测，分析飞机防护区和未防护区NDVI 的时间变化，证明MODIS-NDVI 时间序列数据对春尺蠖虫害监测是可行的。苗静等[4]利用2008—2010年新南威尔士州的MODIS 数据，通过植被生长期前后同期的NDVI 差值图像判断2010年植被长势情况，结果表明监测到的受灾范围与受灾程度基本一致。Czurylowicz[5]以加拿大不列颠哥伦比亚省内虫害发病区域为基础，对叶面积指数LAI（leaf area index）和净初级生产力NPP（net primary productivity）以及净生态系统生产力NEP（net ecosystem productivity）做了研究，指出在受到虫害感染的地区第一年的LAI 变化与第二年的NPP 和NEP 之间存在明显相关性。Buma[6]以落基山脉作为研究区，研究受灾区域叶面积指数LAI与物候指标之间的相关性，表明LAI 在受灾区域内明显下降而物候指标并不会明显变化。郭仲伟等[7-8]利用2002—2012年加拿大不列颠哥伦比亚地区MODIS 数据，分析不同受灾程度地区的叶面积指数LAI、归一化植被指数NDVI 和增强型植被指数EVI 之间的相关性，从而总结不同受灾程度区域内三者之间的相关关系，为利用遥感数据识别虫害提供基础。综上所述，MODIS 数据在监测森林虫害方面已经得到了广泛的应用。

神农架华山松大小蠹是典型的松树蛀干害虫，4月底开始羽化成虫，至10月止，主要危害30年生以上的健康华山松，以幼虫寄生在华山松韧皮部内侧及木质部表层取食韧皮部，钻蛀性强，被大小蠹侵害过的华山松，最快60 d 可致死亡。据2014年统计，虫灾已蔓延到神农架林区的3 个乡镇、4 个国有林场，华山松大小蠹对神农架林区的生态环境造成了严重的威胁[9]。目前，有关识别该林区大小蠧虫害发生区域以及蔓延规律的研究大部分是采用外业调查的方式[9]，但是由于华山松大小蠹虫害侵染面积较大，神农架林区多山峰，这无疑给外业调查带来了巨大困难。随着遥感的发展，利用遥感手段来监测虫害可以有效地解决外业调查困难的问题[10]。例如马望等[11]利用Landsat 数据研究了神农架华山松大小蠹危害程度的划分方法，准确地划分出健康、轻度和重度区域。本研究利用MODIS 数据，通过分区统计法、趋势线分析模型和相关系数模型观测神农架林区华山松大小蠹虫害2008—2015年的受灾趋势，估测虫害发生区域并分析其与海拔之间的相关性，从而判定MODIS 数据在识别植被虫害发生区域及观察虫害时空分布的有效性，期望为森林虫害的预警和科学防治提供技术手段和科学支撑。

1 研究区域和数据处理

1.1 研究区概况

神农架林区位于湖北省西北部，地理坐标为31°15′—31°57′N，109°56′—110°58′E，林区总面积为3 253 km2，林区内多山峰，地势起伏较大，平均海拔为1 800 m，最低点与最高点相对高度差达2 900 m 以上。神农架属于典型的亚热带季风气候，气温偏凉且多雨，林区内气温随海拔变化明显，林间湿度较大。神农架林区作为中国东南西北植被种类的过渡区，林区内动植物资源丰富，垂直自然带谱明显，自下而上主要分为常绿阔叶林、常绿落叶阔叶混交林、落叶阔叶林、针阔混交林和针叶林，其中华山松作为林区内的主要树种，分布面积广，为有害生物的侵害提供了有利条件（图1）。

1.2 数据来源及预处理

MODIS 植被指数产品来自NASA-Land Processes DAAC 数据中心的MOD13Q1（16 d合成产品，空间分辨率为250 m，行列号为h27v05），使用 NASA 官方提供的MODIS Projection Tool 工具提取 NDVI 数据，对所获得的数据进行裁剪便可获得研究区域的NDVI 栅格数据。数字高程模型（DEM）数据来自于空间地理数据云（http://www.gscloud.cn/），空间分辨率为30 m，主要用于研究该区域的地形特点。研究中利用 MODIS Projection Tool 工具对数据进行数据格式及投影的转换，统一采用Krasovsky Albers投影。

图1 研究区位置Fig.1 Location of the study area

2 研究方法

2.1 趋势线分析模型

趋势线分析模型能够模拟每个栅格的变化特征，根据变化特征可以反演出该栅格的植被生长趋势，从而反映研究区在不同研究时段内植被覆盖度变化的空间特征[12-14]。本研究主要利用趋势线分析模型分析神农架林区2008—2015年植被覆盖度年际变化与2014年年内变化情况，从而分析虫害发生区域的空间特征。其公式[14]如下：

式（1）中：n为监测累计年数；Valuei为第i个数据的值（比如第i年NDVI 值）；θslope为趋势线的斜率，若θslope＞0，说明NDVI 在n年间的变化趋势是增加的，反之则减少。

趋势线分析的显著性检验利用归一化植被指数NDVI 序列和时间序列（年份）的相关系数R来检验NDVI年际间变化的显著性，显著性仅代表趋势性变化可置信程度的高低，与变化快慢无关。其公式如下[15]：

通过计算出的相关系数R，查找对应的相关系数显著性检查表，显著性水平设置为0.05，将植被变化趋势分为变化显著（P＜0.05）和变化不显著（P＞0.05）。

2.2 相关系数模型

相关系数模型又称线性相关系数，是衡量变量之间线性相关程度的指标。样本相关系数用r表示，相关系数的取值范围为[-1,1]。本研究主要利用该模型计算林区内像元海拔与NDVI 变化值之间的相关性。利用SPSS 软件进行相关系数计算，可以得到海拔与NDVI 变化值之间的相关系数r以及对应相关程度的显著性水平P值。

式（3）中：Xi和Yi分别表示随机选取像元的NDVI 变化值和DEM 值两组观测数据；和表示随机选取像元的NDVI 变化的平均值和DEM 两组观测数据的平均值。

3 结果与分析

3.1 虫害发生区域空间分布格局

3.1.1 年际虫害发生区域空间分布格局

根据趋势线分析模型，利用2008—2015年的MODIS-NDVI年平均数据得到整个研究区8 a 间的NDVI 变化趋势（图2），趋势线斜率θslope在-0.073 6～0.037 3 之间，其中图2黄色区域趋势线斜率θslope≤0，表示NDVI 值在这8 a 间呈现下降趋势；绿色区域趋势线斜率θslope＞0，表示NDVI 值在这8 a 间呈现上升趋势。

利用ArcGIS 12.0 中的地图代数工具，在0.05显著性水平上对NDVI 变化趋势结果进行显著性检验，得到2008—2015年年际NDVI 显著变化图（图3）。根据检验结果将变化趋势分为3 个等级：显著减少（θslope＜0，P＜0.05）、变化不显著（P＞0.05）和显著增加（θslope＞0，P＜0.05）。黄色区域NDVI 显著减少，灰色区域NDVI 变化不显著，绿色区域NDVI 显著增加。

图2 2008—2015年年际NDVI 变化趋势Fig.2 Interannual NDVI trends from 2008 to 2015

图3 2008—2015年年际NDVI 显著变化Fig.3 Significant change of NDVI during from 2008 to 2015

随着时间的推移，林区内大部分区域的趋势线斜率θslope＞0，而有部分区域的趋势线斜率θslope＜0，NDVI 值呈现下降的趋势。由于处于同一林区内，NDVI 值下降几乎可以排除气候因子等其他影响因子的影响。因此，NDVI 值下降的主要原因可以判定为外来物种的入侵或人工干预。华山松大小蠹钻驻华山松韧皮部，侵入木质层，被侵染的华山松由于无法向上层输入养分，可在两周内导致华山松树冠上层叶尖变黄，一个月后上层针叶全部变为黄色，下层尖端也开始变为黄色，NDVI 值能够有效反映植被冠层生长状况。基于此，将NDVI 显著减少（θslope＜0，P＜0.05）的黄色区域看作虫害发生区域，结果如图3所示。华山松大小蠹虫害的实际调查数据从2008年开始有记载，虫害在自然保护区最早被发现，其后依次为红坪林场、温水林场、红花朵林场、木鱼林场和红坪镇[9]。从图3中可以看出，2008—2015年神农架林区华山松大小蠹虫害受害重度区域主要集中在自然保护区、木鱼林场、木鱼镇、红坪林场和红花朵林场，这与实际调查结果基本一致[9]，说明基于MODIS 数据和趋势线分析模型识别大面积虫害受灾区域是可行的。

3.1.2 年内虫害发生区域空间分布格局

查阅相关文献了解到2008—2015年间的神农架林区华山松大小蠹虫害在2014年最为严重，达到高峰期[10]。由于华山松大小蠹在4月底开始羽化成虫，到10月份止[12]，所以选择2014年第97—305 天的MODIS-NDVI 数据，来分析神农架林区华山松大小蠹虫害年内变化趋势。根据趋势线分析模型，得到2014年NDVI的变化趋势（图4）。同样，黄色区域趋势线斜率θslope≤0，表示NDVI值在2014年第97—305 天间呈现下降趋势；绿色区域趋势线斜率θslope＞0，表示NDVI 值在2014年第97—305 天间呈现上升趋势。

图4 2014年年内NDVI 变化趋势Fig.4 NDVI trends in 2014

同样，对2014年第97—305 天NDVI 的变化趋势进行显著性检验，得到2014年年内NDVI 显著变化图（图5）。其中，黄色区域NDVI 显著减少（θslope＜0，P＜0.05），灰色区域变化不显著（P＞0.05），绿色区域NDVI 显著增加（θslope＞0，P＜0.05）。从图5中可以看出，2014年虫害受灾区域主要集中在自然保护区、酒壶林场、红坪林场和神农架东北部地区。实地调查虫害受灾区域分别为自然保护区、酒壶林场、红坪林场，不包括神农架东北部地区。神农架东北部地区NDVI 的减小，可能和森林采伐作业有关[11]。

3.2 虫害随海拔变化的分布特征

3.2.1 不同海拔范围内虫害分布特征

通过分析神农架华山松大小蠹虫害年际和年内发生区域分布格局，发现大面积爆发的区域都为海拔较高的区域，例如神农架自然保护区，海拔在2 500 m 以上的山峰有20 多座，其中最高的神农顶高达3 105.4 m。因此将研究区划分为4 个不同海拔等级，分别为DEM≤1 600 m、1 600 m＜DEM≤2 000 m、2 000 m＜DEM≤2 400 m、DEM＞2 400 m（图6），以此来研究虫害发生区域与海拔高度之间的相关关系。利用ArcGIS 中分区统计方法，统计2008—2015年不同海拔范围内虫害发生面积，统计结果见表1。从表1中可以看出，当DEM＜1 600 m 的区域内虫害发生面积占总面积的18.59%，而在海拔位于2 000 m＜DEM＜2 400 m 的区域内虫害发生面积占总面积的比值达到41.52%，特别是在海拔高于2 400 m 区域内虫害发生面积占总面积一半以上。由此可见，随着海拔的增高，虫害发生面积占该区域总面积的比例逐渐增大，这也说明了在海拔较高的区域虫害爆发的范围增大，海拔较高的区域应该为防治虫害的重要区域。

图5 2014年年内NDVI 显著变化Fig.5 Significant change of NDVI in 2014

图6 海拔分区Fig.6 Distribution of different altitudes

表1 不同海拔范围内虫害发生面积统计Table 1 Statistics on the area of Dendroctonus armandi infestation at different altitudes

3.2.2 虫害区域与海拔的相关性

通过对不同海拔范围内虫害发生区域分布特征的分析，2008—2015年神农架大小蠹虫害发生区域主要集中在海拔2 000 m 以上的地区，表明虫害发生区域与海拔之间存在相关性。为进一步探究虫害发生与海拔之间的相关性，在研究区内虫害发生区域随机均匀地选取了43 个观测点，利用ArcGIS 工具获得这些观测点的NDVI 值和DEM 值，利用SPSS 工具进行相关性分析，得到虫害发生区域内NDVI 值与海拔之间存在较为显著的负相关关系（相关系数r=-0.491，显著性水平P=0.001）。这里相关程度不是很高的主要原因在于海拔只是影响虫害众多因素中较为重要的一个，虫害的发生还受其他一些环境因素的影响，例如温度、降水、坡度、坡向和风向等。

3.3 不同海拔范围内虫害蔓延趋势分析

通过对2008—2015年8 a 总体虫害发生区域的地形分布特征分析，了解到在海拔较高的区域虫害发生的几率较大。同样由于华山松大小蠹钻驻华山松韧皮部侵入木质层后导致华山松树冠上层叶尖变黄，NDVI 值能够有效反映植被冠层生长状况。不同的海拔范围内华山松的生长期不同，在统计不同海拔区域内NDVI 均值的变化过程中发现在华山松的生长期内存在NDVI 均值下降的现象，通过调查发现华山松大小蠧羽化成虫和扬飞的时间受温度的影响，这也就导致了不同海拔范围内大小蠧虫害爆发的时间不同。结合不同海拔范围内NDVI 均值的变化趋势以及大小蠧的生活习性分析可得到不同海拔范围内虫害蔓延趋势。因此利用2014年第97—305 天的MODIS-NDVI数据，同样利用ArcGIS 分区统计工具获得不同海拔等级内NDVI 的均值，统计结果见表2。根据不同海拔区域内NDVI 均值的变化以及大小蠧的生活习性，观察虫害的蔓延规律。

表2 不同海拔区域内NDVI 均值Table 2 NDVI mean values of forest areas at different altitude

从图7中可以看出，NDVI 随时间变化总体呈现先增加后减小的趋势，这一趋势符合1年内第97—305 天植被生长规律。根据图6分析，在海拔低于1 600 m 的区域内NDVI 基本在平稳的变化，无异常点，说明该区域内植被生长状况良好。在1 600 m＜DEM＜2 000 m 区域内，第177—193天NDVI 值骤减，代表该区域在这个时间遭受虫害侵染并具有一定的表现，大小蠧属于钻蛀害虫，在侵染华山松两周后树冠上层针叶才开始脱水变黄[16-17]，由此推测该海拔范围内大小蠹在6月中上旬开始入侵，这与张子一等[17]研究的华山松大小蠹在海拔1 700 m 样地内6月中旬开始入侵10月上旬结束不谋而合。同样，分析海拔高于2 000 m低于2 400 m 的区域内NDVI 值变化曲线，可以发现在第193—209 天NDVI 值骤减，推测大小蠹在7月中上旬开始入侵。在海拔高于2 400 m 的区域第193—209 天NDVI 值有一个较小幅度得减少，第209—225 天减小幅度增大，推测在海拔高于2 400 m 的区域内，大小蠹在7月下旬或8月初开始入侵。这些与前人研究的在海拔2 100 m 的样地内7月中旬开始入侵至9月下旬结束，海拔2 500 m的样地内7月下旬开始入侵9月下旬结束的结果基本保持一致。由于研究区内的温度和降水每年都存在一定的差异，所以虫害每年发生的时间和区域也都会随着发生变化，但总体上与张子一等[17]研究的大小蠧虫害的生活习性保持一致。由此，可以得到不同海拔区域内虫害发生的时间不一样，这与华山松大小蠹的繁殖与羽化受温度影响具有直接关系。

图7 不同海拔区域内NDVI 均值变化趋势Fig.7 Trends of NDVI mean values at different altitudes

4 结论与讨论

4.1 结 论

利用2008—2015年的MODIS-NDVI 数据，对神农架林区华山松大小蠧虫害进行监测。主要从年际、年内2 个角度监测了虫害的发生区域，然后分析虫害发生区域与海拔之间的相关关系，最后分析了1年内不同海拔区域内虫害爆发的时间。

1）利用MODIS 数据监测华山松大小蠧虫害分布区域和虫害发生区域内的受灾程度与实际调查数据基本保持一致。在2008—2015年8 a 间，神农架林区华山松大小蠹虫害受害区域主要集中在自然保护区、木鱼林场、木鱼镇、红坪林场和红花朵林场，这与实际调查结果基本一致。根据2014年第97—305 天间NDVI 的变化值得到，2014年虫害受灾区域主要集中在自然保护区、酒壶林场、红坪林场，也与实际调查数据保持一致。可见，利用MODIS 数据识别大区域虫害受灾区域是可行的。

2）不同海拔高度范围内，大小蠹虫害受灾区域面积与总面积之比随海拔高度的增加而增大，大小蠹虫害受灾区域内的NDVI 值与海拔之间存在较为显著的负相关关系，相关系数r=-0.491，显著性水平P=0.001。

3）由于不同海拔高度存在的温度差异造成华山松大小蠹羽化成虫与扬飞的时间不同，最终导致虫害在不同海拔范围内爆发的时间不同，且每年由于气候的原因虫害爆发的时间存在着差异。

4.2 讨 论

本研究基于大小蠧侵害华山松后冠层树叶变黄的特点，以MODIS-NDVI 数据为基础，利用趋势线分析模型、相关系数法、分区统计方法来观察2008—2015年神农架大小蠧虫害发生区域的空间分布特征，并重点讨论了虫害发生区域与海拔之间的相关关系，结合大小蠧的生活习性分析了大小蠧虫害爆发时间与海拔之间的相关关系。

本研究局限于地形对大小蠧虫害的影响，没有分析气候对虫害爆发的影响，下一步应结合气候因素分析病虫害的爆发是否与气候条件相关，从而有利于全面掌握影响大小蠧虫害爆发的原因，期望能根据气候预测值以及地形的调查值对未来病虫害的爆发时间作出较为准确的预测。

由于MODIS 数据空间分辨率较低，对于神农架这样范围较大的林区研究虫害分布的大致位置具有一定的价值，但对于观察范围较小或研究大小蠧虫害侵害树木更加准确的范围时仅仅使用MODIS 数据是远远不够的。因此在接下来的研究中，也应该选择分辨率较高的遥感影像，挖掘更多的遥感图像光谱指数特征，从而更加准确地定位虫害发生的区域以及分析影响虫害爆发的各种因素，为病虫害的防治工作提供更加科学的依据。