云南省墨江县森林对虫害和干旱干扰的响应特征

曹雪佳，沈欣悦，沈 亲，黄华国

（北京林业大学 省部共建森林培育与保护教育部重点实验室，北京 100083）

植被生长受到气候变化的影响，温度和降水量是气候的2个最主要因子。探究植被变化与气候因子的关系是研究生态系统对全球气候变化响应研究的基础[1]。增强型植被指数EVI的动态变化可以体现植被状态对全球气候变化的响应程度[2]，通过植被指数EVI与气候因子的相互作用，探究全球及区域气候、生态系统的反馈响应和滞后效应是当下的研究热点。

云南素有四季不分明、干湿季明显的气候特征，冬春季为季节性干旱期，降水稀少且蒸发量大，极易爆发干旱[3]。近年来，云南干旱灾害层出不穷，其中又以2009—2010年西南干旱事件最为严重。墨江哈尼族自治县（以下简称墨江县）是云南省普洱市下辖县之一，旱灾爆发频率高、影响范围广并且持续时间长[4]。其中2009—2010年的西南干旱事件是该地区有气象记录以来最为严重的干旱事件[5]。干旱又进一步诱发了大量病虫害，前后两年的连旱会进一步促发松毛虫大暴发[6]。

为探明在干旱和虫害干扰下云南森林植被的状态及其变化机理，本文以墨江县为例，基于MODIS卫星的EVI连续遥感产品数据、多年气象数据和研究区虫害信息，探究2006—2015年的云南植被指数时间序列的变化规律及气候变化趋势，试图阐明云南植被覆盖在不同生长状态下（干旱、虫害等）与气象因子间的相关性以及植被生长主要限制因子，同时进一步研究植被对气象因子的滞后效应及滞后期，此研究对于林业上开展干旱与植被覆盖变化及植被对气候变化的响应的研究起到技术和理论上的支撑作用，同时有助于帮助决策部门采取积极有效的防旱抗旱措施，将干旱灾害对森林植被带来的损失降到最低。

1 研究区概况与数据

1.1 研究区概况

墨 江 县 地 处 22°51′～ 23°59′N、101°08′～102°04′E 之间，面积 5 312 km2（见图 1），低纬度高海拔。四季冷暖不分明，但有明显干湿两季：墨江县当年11月到次年4月属干季，晴多雨少，气候较为干燥；5—10月属湿季，晴少雨多。降水量空间差异较大，呈南多北少之势。其森林覆盖率高至59.5%，森林类型多样，主要包括思茅松、云南松、季节性河谷和亚热带常绿阔叶林[7]。其中，思茅松Pinus kesiyavar.langbianensis是墨江县林区主要树种之一。

思茅松毛虫Dendrolimus kikuchii，鳞翅目枯叶蛾科松毛虫属，以松针为食，是危害思茅松的主要森林虫害。由于思茅松一年抽发两梢，其中春秋各一梢。相应地，思茅松毛虫一年在林区有二个世代，第一个世代（夏代）大约有4—5个月的生活周期，第二代（越冬代）有长达8—9个月的周期。一代虫害爆发盛期在每年7—8月，二代虫害爆发盛期在每年3—4月。松毛虫害发生过的思茅松林多成片死亡，被暴食后的针叶多枯黄，与火烧痕迹相似，被称为“不冒烟的森林火灾”。

气象资料显示，2009—2010年研究区全面严重干旱，虽有持续中度干旱。另外当地森防站资料显示，2009、2010、2013和2015年研究地区松毛虫大面积爆发，伴随有少量小蠹虫危害。

图1 墨江县地理位置及土地利用Fig.1 Geographical location and land use of Mojiang

1.2 数据源及预处理

1.2.1 遥感数据

在全球区域和大尺度的植被覆盖动态研究中，MODIS/EVI产品连续遥感数据被大范围运用[8]。MODIS/EVI产品利用了其辐射计的优势，经蓝波段减少了烟雾以及子像元薄云残余的大气污染并耦合冠层背景信号，达到削弱大气影响，提高植被监测精度之目的，于此同时，订正后的地表反射率对与高植被覆盖区的反应更敏感，保持了植被状态的稳定性。相较于NDVI，EVI的时间序列季节性更明显，其特点主要表现在高植被覆盖敏感性、降低土壤背景影响、消除大气传输干扰等，可以更好的反映高植被覆盖区季节性变化特征，时间序列曲线比较平滑[9]。同时，墨江县土地利用主要是林地、农田、水体等，其他土地利用类型面积相对较小（见图1）。因此，本研究选取MODIS/EVI产品作为植被覆盖度时空变化监测数据。

对MODIS/EVI进行投影变化、最大值合成等数据预处理，获得墨江县2006—2015年空间分辨率1 km，时间分辨率1个月的EVI时空序列数据集。

1.2.2 土地利用数据

结合2010年11月初Landsat7影像，选择最大似然法监督分类法将研究区分为4种类别：林地、农田、水体和其他。综合野外调查数据和Google Earth图像，选取验证样本进行验证。分类结果精度是88.24%，其Kappa系数是0.84，分类可靠性较高。最后将土地利用分类图的林地分类结果和EVI图层进行叠加，最终获得EVI值。

1.2.3 气象数据

本文所用气象研究数据均来源于中国气象局气象数据中心（http://data.cma.cn/），主要包括云南省全境站点的逐月降水量和温度数据等。应用ArcGIS进行克吕格空间插值，获得墨江县的平均气象数据（见图2），降水量呈现明显的由北往南，由东向西增加的趋势。

1.2.4 虫害数据

图2 降雨量分布插值Fig.2 The precipitation interpolation

当地森防部门提供了全县宏观的虫害年份和每年防治面积。本文将研究时段分别划分为无明显虫害和干旱年（2007—2008年）、显著干旱年（2011—2012年）、虫害年（2013和2015年）和干旱和虫害同时爆发年份（2009—2010年）。

野外调查于2014年在危害典型的联珠镇进行。调查步骤如下：在松毛虫危害的林分中设置了25块样地，其中轻度受害9个、中度受害8个、重度受害4个，以及往年受害的但现在已经更新恢复的4个样地，样地具体坐标及受害分布状况如下图3。

2 研究方法

2.1 低谷深度指数DEVI分析

图3 墨江县联珠镇危害现状Fig.3 Location of the insect pest plot in Lianzhu

夏代松毛虫的爆发盛期是7—8月，松毛虫危害过后的思茅松，生长速度缓慢，危害更为严重时思茅松叶被完全吃光，在8—9月EVI会呈现低谷。基于此，本文提出了低谷深度指数（DEVI）构建虫害年份判别模型：

式中，DEVI表示EVI低谷深度指数，i代表EVI最高月份，EVIi、EVIi+1和EVIi+2分别代表第i、i+1和i+2个月份的EVI。

当DEVI显著大于0时，极有可能为虫害年。

2.2 EVI与气象因子的偏相关分析

对2个及以上具有相关性的变量研究分析，以此衡量变量之间的相关性及其密切度称为相关性分析。公式如下：

式中，rxy是变量x与y两者的相关系数。n代表样本数量，和分别代表变量的平均数。xi和yi分别代表第i个月的两个变量值。

为了消除其他变量影响的条件下研究单独两个变量间的相关程度，本文采用偏相关分析研究了植被与气象因子（温度和降雨）相关性。公式如下：

式中，rxy,z为变量z固定后y与x的偏相关系数。rxy、rxz、ryz分别为x和y变量、x和z变量、y和z变量的相关系数。

2.3 滞后性分析

对林地EVI和前3个月、2个月、1个月以及当月的气象因子相关系数大小进行比较，分析研究林地EVI对于气候变化的响应时滞，用最大偏相关系数相对的滞后时间来代表滞后期。

进一步，构建2008—2014年林地EVI和降水量、温度这二者的多元线性回归模型，包括同期及滞后模型，依据气象因子变化来定量预估植被生长状况，探究EVI改变的原因。选择同期降水量和温度作为同期模型中变化因子，选择对应滞后期前的降水量及温度作为滞后期模型中的因子。

3 结果与分析

3.1 墨江县降水量和温度动态变化趋势

2006—2015年，墨江县多年平均年降水量为1 431.3 mm，呈现明显的季节波动，最大降水量发生在7—8月份，达到600.8 mm，占到年降水量的41.7%。最小降水量发生在2月和12月，达到34.7 mm，占到年降水量的2.4%（见图4a）。多年年均温度为19.6 ℃，同降水量分布一样，呈现明显的季节波动。最高温度出现在6月份，为23.5 ℃，最低温度出现在1和12月，为14.2 ℃（见图4b）。雨热同季分明且冬季温度较高，为松毛虫幼虫生活提供了有利环境，使得幼虫安全越冬，这是促使虫灾爆发的诱因之一。

图4 2006—2015年墨江县降水和温度的月动态变化趋势Fig.4 Temporal variations of precipitation and air temperature in the Mojiang during 2006—2015

图5是近10年墨江县降水量和温度距平。2006—2007年，降水量维持在平均值附近，2009年较2008年降水量较低。2010—2012年间出现连续负距平，较均值水平19.3%以下。2013年以后，降水量均在平均值以上。由降水量趋势可以看出，从2009年，降水量出现下降趋势。云南在2009年遭遇50年一遇严重旱情，到2010年百年一遇大旱，2011—2012年，旱情持续加重。本文将2009—2012年定义为干旱年。2007—2008年，温度距平低于平均值，2010年以后，温度距平高于平均值，2011年除外，温度整体呈升高趋势。

3.2 墨江县EVI动态变化和DEVI指数

图5 2006—2015年墨江县降雨和温度距平Fig.5 Anomaly of precipitation and air temperature in Mojiang county from 2006 to 2015

EVI月动态变化与降水和温度变化趋势一致（见图6）。EVI多年均值为0.41；最大值出现在7—8月份，其中，2011年7月达到0.54；最小值出现在2月份，其中，2010年3月份只有0.27。

图6 2006—2015年间墨江县林地EVI月变化趋势Fig.6 Temporal variations of EVI in Mojiang city from 2006 to 2015

多年低谷深度指数见表1。可以看出，虫害对应年的EVI低谷深度指数（0.045±0.014）显著高于非虫害年（-0.003±0.011）（t检验，p=0.001 9），表明DEVI指数可以很好的辨识虫害爆发的年份。

表1 2006年—2015年对应DEVI指数Table 1 The variations of DEVI index during the period from 2006 to 2015

3.3 EVI和降雨及温度的相关性研究

在年变化范围中，林地年平均EVI和年平均降水量、温度之间相关系数及偏相关系数都很低同时没有通过显著性检验。这反映了降水量、温度对年EVI值影响不大。进一步探究，发现在年内月尺度上，温度和降水量变化的幅度明显大于年际变化。EVI和降水量的关系显著相关（如2008年，r=0.740）。可见降水量和温度对林地EVI发挥着至关重要的作用。温度升高可以延长生长季时间，增加光合作用时间，最终促进植被的生长[10]。降水量增加同时伴随温度增高可以增加生态系统生产力，促进陆地碳吸收[11]。反映在遥感数据上，降水量和温度的增加会促进EVI增长，反之EVI则会下降。因此，总体上年均EVI持平，也解释了年际变化的不显著性。

为分析虫害和干旱的影响，将研究期间时间段分别划分为正常年（2007—2008年）、干旱年（2011—2012年）、虫害年（2013和2015年）和干旱虫害年（2009—2010年），分别与降水量和温度进行偏相关分析（见表2）。结果表明，降水量对EVI的变化解释程度由高到低依次为正常年（r=0.740，p=0.001）、虫害年（r=0.602，p=0.002）、干旱年（r=0.479，p=0.002）、干旱虫害年（r=0.480，p=0.002）。在干旱和虫害胁迫下，植被生长受到干扰，对气候的敏感性会受到抑制。

温度对EVI的影响未通过0.05水平检验。严建武等认为干旱期间西南地区对干旱敏感的用地类型如旱地和草地对气象因子变化反应剧烈，相关性明显比森林更强[12]。本文未考虑用地类型，在进一步的研究中，分用地类型讨论，研究结果可能会更客观具体。

表2 不同环境条件下林地EVI与同期气象因子偏相关系数†Table 2 Correlation coefficient for forest EVI and meteorological factor under different environmental conditions

3.4 降雨和气温对EVI的滞后影响

如表3所显示，林地EVI对降水量响应最强的是前1个月，反映为EVI和前1个月的降水量的偏相关系数最大，达到0.522。且林地EVI对前2个月的温度响应最强，表现为EVI与前2个月的温度的偏相关系数最大，达到0.772。因此可说明林地EVI对于降水量及温度响应的滞后期分别是1个月和2个月，且林地EVI对温度变化的时滞性要更强，说明前1个月的降水量和前2个月的温度会影响EVI的变化趋势。Piao等对全国植被NDVI与温度的研究发现，温度对NDVI有3个月的滞后期[13]。Kaspersen等发现生长季NDVI与降水量有1个月滞后期[14]。Xu等对全国分区NDVI变化趋势的研究发现，西部、中部和北部，NDVI与前1个月的温度和降水量显著相关[15]；南部可能由于丰富的降水量，NDVI与降水量和温度有2个月的滞后期。因此，NDVI对气候的响应普遍存在滞后效应[16-17]。本文研究结果和前人结果类似，滞后期则略有不同。

表3 2006—2015年月EVI与同期及前1－3个月气象因子（温度和降水量）偏相关系数分析†Table 3 Partial correlation coefficient between monthly EVI and meteorological factors (precipitation and air temperature) during the period from 2006 to 2015 with 0－3 month lag effect

3.5 根据前期降雨和气温预报后期EVI

林地EVI和降水量及温度的多元线性回归研究见表4。研究结果中滞后模型选择了前1个月的降水量和前2个月的温度，可以发现滞后模型的拟合精度较高，该结果和陈燕丽等在喀斯特地区的研究结果类似[18]。

结合平均EVI数据，滞后预测模型可以用于预测下一个月的森林植被的衰退可能性，有利于灾害防控。

表4 2006—2015年EVI与气象因子（温度和降水量）多元线性回归模型†Table 4 Multivariate linear regression model between EVI and meteorological factors (precipitation and air temperature) during the period from 2006 to 2015

4 结论与讨论

本研究以墨江县为例，基于MODIS卫星的EVI连续遥感产品数据、多年气象数据和研究区虫害信息，探究2006—2015年的云南植被指数时间序列的变化规律及气候变化趋势，开展了植被指数对气象数据和虫害信息的敏感性分析，并提出了一个新的低谷深度指数DEVI和一个考虑滞后性的EVI预测模型，结论如下：

（1）EVI对温度和降水量的响应存在时间尺度效应。降水量、温度对年EVI值影响不大，在月尺度上，EVI与同期降水量存在显著正相关（p＜0.05），降水量和温度的增加会促进EVI增长，反之EVI则会下降，这与张培松等的研究结果相似，他研究发现在我国大多数地区内EVI与降水量和温度等气象因子呈一定正相关性[19]。张艺成研究发现2000—2013 年九万山地区月均温度与植被 NDVI 呈现显著的正相关关系[20]。

（2）EVI对温度和降水量的响应存在滞后性，滞后期分别为2个月和1个月，植被对温度的滞后响应高于对降雨的滞后响应。考虑降水量和温度的滞后模型可以有效预测下个月的EVI（r=0.84）。陈燕丽研究指出喀斯特地区增强型植被指数EVI对降水量、气温、水汽压等气象因子响应的滞后性非常明显，据测算滞后期大概为16天[18]。何月等研究得出EVI对降水量的最大响应存有明显的时滞性，滞后期大概为1个月左右。而大部分地区植被生长季节变化对气温和降水量具有明显的时滞效应，滞后期大概为50天[22]。本文研究结果和前人结果类似，滞后期则略有不同。

（3）在干旱和虫害胁迫下，植被生长受到干扰，对气候的敏感性会受到抑制。研究中温度对EVI的影响未通过0.05水平检验，这与此前学者对西南地区的研究结果相符合，即干旱期间西南地区对干旱敏感的用地类型如旱地和草地对气象因子变化反应剧烈，相关性明显比森林更强[12]。

（4）虫害对应年的EVI低谷深度指数DEVI（0.045±0.014）显著高于非虫害年（-0.003±0.011）（t检验，p=0.001 9），表明新的低谷深度指数DEVI可以很好的辨识虫害爆发的年份。

本文研究区域较小，在接下来的研究中需在更大尺度范围内进行分析。同时本文并没有把植被类型考虑进去，综合严建武等人的研究，干旱期间西南地区对干旱敏感的用地类型如旱地和草地对气象因子变化反应剧烈，相关性明显比森林更强[12]，在进一步的研究中，分用地类型讨论，研究结果可能会更客观具体。

在进一步的研究中还可以构建EVI与气象因子空间网格拟合模型，计算出精细格点尺度上的气象条件再分配结果，利于更准确地分析出植被与气候的相关性。同时，植被覆盖还会受到土壤类型、岩性等其他因素的影响，这在模型的优化中还有待于深入研究。

[1]韩春光.新疆植被指数与气象因子的响应关系研究[D].乌鲁木齐:新疆大学, 2009.

[2]彭代亮,黄敬峰,王秀珍.基于MODIS-EVI区域植被季节变化与气象因子的关系[J].应用生态学报, 2007,18(5):983-989.

[3]张冬冬,鲁 帆,严登华,等.云南省干旱时空演变规律及季节连旱的概率特征分析[J].应用基础与工程科学学报,2014(4): 705-717.

[4]彭贵芬,刘 瑜,张一平.云南干旱的气候特征及变化趋势研究[J].灾害学,2009,24(4):40-44.

[5]彭贵芬,刘盈曦.2009—2010年云南特大干旱动态风险分析[J].气象科技进展,2012,2(4):50-52.

[6]王学锋,周德丽,童 青,等.思茅林区松毛虫害发生与气象因子的关系[J].西部林业科学, 2003(3):77-82.

[7]郑天水.墨江县森林资源特点及林业发展对策[J].林业调查规划,2000,25(1):1-4.

[8]Galvão L S, Santos J R D, Roberts D A,et al.On intra-annual EVI variability in the dry season of tropical forest: A case study with MODIS and hyperspectral data[J].Remote Sensing of Environment, 2011, 115(9):2350-2359.

[9]李文梅,覃志豪,李文娟,等.MODIS NDVI与MODIS EVI的比较分析[J].遥感信息,2010(6):73-78.

[10]Slayback D A, Pinzon J E, Los S O,et al.Northern hemisphere photosynthetic trends 1982-99.Global Change Biol 9:1-15[J].Global Change Biology, 2002, 9(1):1-15.

[11]Nemani R R, Keeling C D, Hashimoto H,et al.Climate-driven increases in global terrestrial net primary production from 1982 to 1999.[J].Science, 2003, 300(5625):1560.

[12]严建武,陈报章,房世峰,等.植被指数对旱灾的响应研究——以中国西南地区2009年—2010年特大干旱为例[J].遥感学报,2012,16(4):720-737.

[13]Piao S, Fang J, Zhou L,et al.Interannual variations of monthly and seasonal normalized difference vegetation index (NDVI) in China from 1982 to 1999[J].Journal of Geophysical Research Atmospheres, 2003, 108(D14):ACL 1-1.

[14]Kaspersen P S, Fensholt R, Huber S.A Spatiotemporal Analysis of Climatic Drivers for Observed Changes in Sahelian Vegetation Productivity (1982–2007)[J].International Journal of Geophysics, 2011,(2011-12-27), 2011, 2011(1).

[15]Xu G, Zhang H, Chen B,et al.Changes in Vegetation Growth Dynamics and Relations with Climate over China’s Landmass from 1982 to 2011[J].Remote Sensing, 2014, 6(4):3263-3283.

[16]Nezlin N P, Kostianoy A G, Li B L.Inter-annual variability and interaction of remote-sensed vegetation index and atmospheric precipitation in the Aral Sea region.[J].Journal of Arid Environments, 2005, 62(4):677-700.

[17]J.Wang, P.M.Rich, K.P.Price.Temporal responses of NDVI to precipitation and temperature in the central Great Plains, USA[J].International Journal of Remote Sensing, 2003, 24(11):2345-2364.

[18]陈燕丽,莫伟华,罗永明,等.基于气候信息的喀斯特地区植被EVI模拟[J].农业工程学报,2015, 31(9):187-194.

[19]张培松.气象因子与NDVI的相关性分析[D].重庆:西南大学, 2007.

[20]张艺成,王光军,朱玉林.九万山近14年植被变化与气候因子的相关关系[J].中南林业科技大学学报, 2016,36(7):81-88.

[21]何 月,樊高峰,张小伟,等.浙江省植被NDVI动态及其对气候的响应特征[Z].厦门:第28届中国气象学会年会,2011.