Detecting Shoot Beetle Damage on Yunnan Pine Using Landsat Time-Series Data

Tomicus yunnanensis and Tomicus minor have caused serious shoot damage in Yunnan pine forests in the Yunnan Province of China. However, very few remote sensing studies have estimated the shoot damage ratio (SDR). Thus, we used multi-date Landsat satellite imagery to quantify SDRs and assess the possibility of using spectral indices to determine the beetle outbreak time and spread direction. A new threshold-based classification method was proposed to identify damage levels (i.e., healthy, slightly to moderately infested, and severely infested forests) using time series of moisture stress index (MSI). Permanent plots and temporal field inspection data were both used as references for training and evaluation. Results show that a single threshold value can produce a total classification accuracy of 86.38% (Kappa = 0.80). Furthermore, time series maps detailing damage level were reconstructed from 2004 to 2016. The shoot beetle outbreak year was estimated to be 2013. Another interesting finding is the movement path of the geometric center of severe damage, which is highly consistent with the wind direction. We conclude that the time series of shoot damage level maps can be produced by using continuous MSI images. This method is very useful to foresters for determining the outbreak time and spread direction.


Introduction
Yunnan pine (Pinus yunnanensis Franch) is a local afforestation tree species in southwest (SW) China.It covers around five million hectares of land surface, which accounts for 52% of the forest area in Yunnan Province [1][2][3].More than one-fourth of the territory of the Yunnan Province belongs to the Karst landform.Rock desertification areas have continually expanded in recent years.Because of these reasons, it is very important to preserve the current vegetation cover and to protect Yunnan pine forests [4][5][6][7].
However, due to global climate change, the growth of the Yunnan pine in Yunnan Province has been poor in recent decades [8][9][10][11][12].In 2009, widespread drought occurred in Yunnan Province, which continuously weakened the trees [13][14][15][16][17]. Serious defoliator and woodborer pest damage then followed.Major impacted cities include Kunming, Yuxi, Chuxiong, Qujing, and Dali.Massive death of the Yunnan pine led to huge losses in the forestry industry and reduced ecological function.Because the damage mechanism of defoliators is clear, their monitoring and control are relatively successful [18][19][20][21].Meanwhile, the attacking process of woodborers is still vague and complex.As a result, there are fewer successful case studies to monitor and control these pests.One of these studies focused on the mountain pine beetle (MPB), which attacks only the trunks of pine.Based on canopy color, the damage stages caused by MPB are classified into green, yellow, red, and gray stages [22].Therefore, MPB outbreak time could be reported by aerial detection surveys or field surveys from local agencies, such as the USDA (United States Department of Agriculture) Forest Service [23][24][25].Due to their clear Forests 2018, 9, 39 2 of 14 life cycle and detailed field data, the monitoring accuracy of MPB using remote sensing technology is acceptable [26][27][28][29][30][31][32].
In Yunnan Province, Tomicus spp.are considered to have been responsible for destroying 1.5 million hectares of Yunnan pine forests in SW China over the past 20 years [33][34][35][36].Therefore, it is urgent and crucial to develop an effective method to monitor Tomicus spp.Two members of the Tomicus genus, (1) Tomicus yunnanensis Kirkendall and Faccoli and (2) Tomicus minor Hartig, are the most destructive pests of Yunnan pine in Yunnan Province.To our knowledge, unlike MPB, the life cycle of Tomicus spp. is unclear.It can be roughly divided into two stages: the trunk-shoot stage and the shoot-trunk stage.In the trunk-shoot stage, beetles coming out of trunks attack the fresh shoots of Yunnan pines from May to November.The attacked shoot needles gradually change from green to yellow and red.Correspondingly, the spectral reflectance of the needles changes.In severely infested forests, the beetles may harm whole shoots of Yunnan pine and represent the direct cause of tree mortality.During December, the sexually mature beetles begin to migrate from the shoots to the trunks of weak trees.In the trunk, they mate and lay eggs.During May of the second year, a new trunk-shoot stage begins.New adults drill out of the trunk and move to the shoots for feeding [37][38][39][40][41][42].Because of this complex life cycle, the damage trend of T. yunnanensis and T. minor cannot be effectively mitigated despite many efforts to control its initiation and the spread of the beetles.As a result, it is urgent to gain additional knowledge on the attacking stages of these species.
Remote sensing (RS) is useful in detecting forest disturbances.Therefore, we tried to use RS methods to monitor the damage caused by these shoot beetles.Compared to single-date RS data, time series of Landsat images contain more information to detect pest disturbances, especially disturbances caused by bark beetles.For example, a Landsat time series change detection algorithm (LandTrendr) was used to draw spectral trajectories associated with pest disturbances of different times and severities, as well as to investigate the relationship between these trajectories and ground-based measurements of tree mortality and surface fuels in the Cascade Range of Oregon, USA [43].The researchers found that spectral changes were linked to tree basal area mortality (Adj.R 2 = 0.40) and coarse woody detritus (Adj.R 2 = 0.29).Hais et al. [44] used Landsat imagery to identify the differences in spectral responses between the two types of forest disturbances (i.e., bark beetle outbreak and clear-cuts) and their temporal dynamics in the central part of the Šumava Mountains, at the border between the Czech Republic and Germany in Central Europe.They found that DI (Disturbance Index), wetness, and brightness indices were connected with both forest disturbance types.In addition, the multi-temporal analysis could distinguish different stages of development.Meddens et al. [45] used both single-date (using maximum likelihood classification) and multi-date (using time series of spectral indices) classification methods of Landsat imagery to investigate the change of accuracy in mortality severity in a study monitoring coniferous forests across North America.Results show that the overall accuracy is 91.0% and 89.6% for the single-date image classification (using the tasseled cap transformation indices) and the multi-date analysis (using the Band5/Band4 anomaly), respectively.Therefore, multi-date RS data have the potential to enhance the monitoring accuracy for shoot beetles.
However, since most studies are focused on tree mortality, it is not sufficient to define the shoot damage process.To our knowledge, there has been no previous study using multi-date Landsat imagery to investigate tree shoot damage.As a result, the time series method should be re-evaluated for the monitoring of Yunnan pine.To address this gap, the main objectives of this study were (1) to investigate the capacity of multi-date Landsat satellite imagery to classify damage levels following T. yunnanensis and T. minor attack and (2) to explore the efficacy of spectral indices to detect the beetles' outbreak time and spread direction.

Study Area
The study area (5248 hectares) is located in Pupeng town of Dali City in Yunnan Province (Figure 1).In this area, a major T. yunnanensis and T. minor outbreak was recorded in 2013.The two beetle species have caused extensive shoot damage ratios in the study area.The area undergoes two seasons: the dry season (i.e., November-May) and the wet season (i.e., June-October).The mean average annual precipitation is 783.7 mm, and the average annual temperature is 14.7 • C. The average monthly maximum and minimum temperatures are 27 and 7.6 • C [46], respectively.The elevation ranges from 1720 m (Yunlichang Valley) to 2745 m (Yingge Mountains).

Field Experiments
Forty permanent field plots (30 m × 30 m) were set up in 2016 and distributed as evenly as possible over the study area (see Figure 1c).The geographic coordinates for all plots were measured using differential GPS.At each plot, we sampled all trees with a diameter at breast height (DBH) greater than 4 cm.Three major parameters were measured including DBH, height (H), and the shoot damage ratio (SDR) of each tree.The plot SDR was the averaged shoot damage ratio of single trees (p i ).The average equally treats all trees without considering tree size with regards to weight.The plot size of 30 m × 30 m was to match the spatial resolution of Landsat images.

Study Area
The study area (5248 hectares) is located in Pupeng town of Dali City in Yunnan Province (Figure 1).In this area, a major T. yunnanensis and T. minor outbreak was recorded in 2013.The two beetle species have caused extensive shoot damage ratios in the study area.The area undergoes two seasons: the dry season (i.e., November-May) and the wet season (i.e., June-October).The mean average annual precipitation is 783.7 mm, and the average annual temperature is 14.7 °C.The average monthly maximum and minimum temperatures are 27 and 7.6 °C [46], respectively.The elevation ranges from 1720 m (Yunlichang Valley) to 2745 m (Yingge Mountains).

Field Experiments
Forty permanent field plots (30 m × 30 m) were set up in 2016 and distributed as evenly as possible over the study area (see Figure 1c).The geographic coordinates for all plots were measured using differential GPS.At each plot, we sampled all trees with a diameter at breast height (DBH) greater than 4 cm.Three major parameters were measured including DBH, height (H), and the shoot damage ratio (SDR) of each tree.The plot SDR was the averaged shoot damage ratio of single trees (pi).The average equally treats all trees without considering tree size with regards to weight.The plot size of 30 m × 30 m was to match the spatial resolution of Landsat images.

⁄ 100%
(1) where p is the single tree damage ratio, ni is the number of damaged shoots, Ni is the number of total shoots, and SDR is the averaged plot damage ratio of ntree trees.Each SDR is between 0% and 100%.We grouped the SDRs into two groups: that of an infested tree (0% < SDR < 100%), and a "red-attack" tree or a "gray-attack" tree (SDR = 100%).Based on the classification standard (Standard of Forest Pests Occurrence and Disaster LY/T 1681-2006) and considering the local damage situation, we divided the forests into three damage levels: healthy forest (0% ≤ SDR < 10%), slightly to moderately infested forest (10 ≤ SDR ≤ 50%), and severely infested forest (50% < SDR) (Table 1).
Based on the field plots, the boundaries of the three different damage level areas (healthy, slightly to moderately infested, and severely infested) were delineated through field inspections, where p is the single tree damage ratio, n i is the number of damaged shoots, N i is the number of total shoots, and SDR is the averaged plot damage ratio of n tree trees.Each SDR is between 0% and 100%.We grouped the SDRs into two groups: that of an infested tree (0% < SDR < 100%), and a "red-attack" tree or a "gray-attack" tree (SDR = 100%).Based on the classification standard (Standard of Forest Pests Occurrence and Disaster LY/T 1681-2006) and considering the local damage situation, we divided the forests into three damage levels: healthy forest (0% ≤ SDR < 10%), slightly to moderately infested forest (10 ≤ SDR ≤ 50%), and severely infested forest (50% < SDR) (Table 1).
Based on the field plots, the boundaries of the three different damage level areas (healthy, slightly to moderately infested, and severely infested) were delineated through field inspections, surveyed Forests 2018, 9, 39 4 of 14 with GPS, and marked on an ortho-photo from 2016 (see Figure 1c).We then used these specified areas to assess the spectral response of the different SDRs from attacks by the beetles.
Table 1.Three classes of Tomicus spp.severity using shoot damage ratios.

Landsat Data Preparation
Seven archived Landsat Thematic Mapper (TM) 5, ten Landsat Enhanced Thematic Mapper (ETM+) 7, and five Landsat Operational Land Imager (OLI) images from 2004 to 2016 were downloaded from the USGS GLOVIS (United States Geological Survey Global Visualization Viewer) website [47] (Table 2).The chosen images were acquired at time periods close to the end of shoot damage by these beetles (i.e., in November/December) and with minimal cloud cover.If there was significant cloud cover and SLC (Scan Lines Corrector)-off areas in a scene, we attempted to find other images within that same year to fill cloudy areas and SLC-off areas in the priority image.Using visual inspection, we noted that Landsat data from path/row 130/42, 130/43, and 131/42 were terrain-corrected (processing level: L1T) which resulted in high geo-registration accuracy between the images.The satellite scenes were orthorectified into the coordinate system and a DEM (Digital Elevation Model) was created from the contour lines (1:25,000).The calibration and atmospheric correction model (FLAASH) were applied to the images to convert digital number (DN) values into values for sensor radiance and reflectance at the bottom of the atmosphere.

Spectral Indices
The spectral regions of Landsat sensors are shown in Table 3.We calculated four spectral indices to aid in image classification (Table 4): the red-green index (RGI), normalized difference vegetation index (NDVI), normalized difference moisture index (NDMI), and moisture stress index (MSI).These four spectral indices were used to detect SDR because they have been successfully used to identify pest disturbances in previous research [48][49][50][51][52][53].

Multi-Date Classification Method
Our study employed a Landsat time series change detection algorithm, which is similar to the one used by Meddens et al. [45].A temporal anomaly is a straightforward anomaly approach to reflect the difference of a spectral index value between a given time and the multi-temporal, undisturbed mean.If a bias threshold is defined away from the temporal undisturbed mean, temporal anomalies of the spectral indices (RGI, NDVI, NDMI, and MSI) can be detected in any given year if the bias is greater than the threshold value.Figure 2 shows the example of using MSI time series data to detect the anomaly years, where the temporal undisturbed mean is near 0 and the threshold line is around 0.13.If the MSI of one year is greater than 0.13, this year is highly likely to be an anomalous year.
In this study, we randomly selected 70% of the pixels from the healthy forest, slightly to moderately infested forest, and severely infested forest (n = 420 per class) to use as training data to determine the optimal threshold values when the maximum overall classification accuracy has been reached (Figure 3).For example, in Figure 3b, overall accuracy decreased when the value was above or below this threshold.Values above this threshold identified more forests with slight to moderate infestation than forests with severe infestation; in values below this threshold, the opposite was true.The remaining 30% of the pixels were used for validation (n = 180 per class).Using different anomaly threshold values, we generated confusion matrices [54] for each spectral anomaly.For each spectral index anomaly, we identified the highest overall accuracy across this range of thresholds in separating different beetle infestation locations.In this study, we randomly selected 70% of the pixels from the healthy forest, slightly to moderately infested forest, and severely infested forest (n = 420 per class) to use as training data to determine the optimal threshold values when the maximum overall classification accuracy has been reached (Figure 3).For example, in Figure 3b, overall accuracy decreased when the value was above or below this threshold.Values above this threshold identified more forests with slight to moderate infestation than forests with severe infestation; in values below this threshold, the opposite was true.The remaining 30% of the pixels were used for validation (n = 180 per class).Using different anomaly threshold values, we generated confusion matrices [54] for each spectral anomaly.For each spectral index anomaly, we identified the highest overall accuracy across this range of thresholds in separating different beetle infestation locations.In this study, we randomly selected 70% of the pixels from the healthy forest, slightly to moderately infested forest, and severely infested forest (n = 420 per class) to use as training data to determine the optimal threshold values when the maximum overall classification accuracy has been reached (Figure 3).For example, in Figure 3b, overall accuracy decreased when the value was above or below this threshold.Values above this threshold identified more forests with slight to moderate infestation than forests with severe infestation; in values below this threshold, the opposite was true.The remaining 30% of the pixels were used for validation (n = 180 per class).Using different anomaly threshold values, we generated confusion matrices [54] for each spectral anomaly.For each spectral index anomaly, we identified the highest overall accuracy across this range of thresholds in separating different beetle infestation locations.

Outbreak Time and Direction Spread Estimation
Except for a general record on whether there are outbreaks of beetles in each year, there are no historical records on the spatial and temporal information of beetle outbreaks.Therefore, we used the classification method and optimal threshold values to classify images from 2004 to 2016 and reconstruct the progression of beetle spread.The outbreak time for each pixel is defined as the year where an abrupt change larger than the threshold value occurred.It should be noted that one pixel may have an outbreak more than once.
By estimating the geometrical centroid of all outbreak pixels, the outbreak source for each year can be identified.The pixels of severely infested forests were screened as the outbreak forests.When all the outbreak source points by year were plotted, the moving or spreading path of shoot beetles could Forests 2018, 9, 39 7 of 14 be extracted.Since there was no beetle outbreak before 2010, the outbreak sources from 2011 to 2016 were calculated in ESRI ArcGIS.Otherwise, in order to explore the vector that could explain the spread direction of beetles, we calculated the percentage for wind direction in the study area from 2011 to 2017 [55].

Field Survey
Of our 40 permanent field plots, 10% had a tree mortality (or gray-attack trees) of >25% in the canopy within a pixel, more than 40% had a tree mortality of <5%, and almost 50% had a tree mortality from 5% to 25% (Figure 4b).The frequency of red-attack trees within plots is shown in Figure 4c.Red-attack tree frequencies of 0-5%, 5-10%, and 10-15% account for 72.5%, 20%, and 7.5% of plots, respectively.According to field plots, we found that a tree mortality of >25% and a red-attack tree percentage of >10% did not occupy a large proportion of the plots.The results showed that using percentages of tree mortality and red-attack trees to reflect damage levels was not suitable for Tomicus spp.because it neglects shoot damage.On the other hand, healthy forest (0 ≤ SDR ≤ 10%), slightly to moderately infested forest (10 ≤ SDR ≤ 50%), and severely infested forest (50% < SDR) account for 25%, 35%, and 40% of plots, respectively (Figure 4a).Based on the classification standard (Standard of Forest Pests Occurrence and Disaster LY/T 1681-2006), we found that shoot damage ratios could clearly classify the damage levels in these field plots.The results showed that for Tomicus spp., SDR was better than tree mortality or red-attack trees to reflect the damage degrees of forests.
reconstruct the progression of beetle spread.The outbreak time for each pixel is defined as the year where an abrupt change larger than the threshold value occurred.It should be noted that one pixel may have an outbreak more than once.
By estimating the geometrical centroid of all outbreak pixels, the outbreak source for each year can be identified.The pixels of severely infested forests were screened as the outbreak forests.When all the outbreak source points by year were plotted, the moving or spreading path of shoot beetles could be extracted.Since there was no beetle outbreak before 2010, the outbreak sources from 2011 to 2016 were calculated in ESRI ArcGIS.Otherwise, in order to explore the vector that could explain the spread direction of beetles, we calculated the percentage for wind direction in the study area from 2011 to 2017 [55].

Field Survey
Of our 40 permanent field plots, 10% had a tree mortality (or gray-attack trees) of >25% in the canopy within a pixel, more than 40% had a tree mortality of <5%, and almost 50% had a tree mortality from 5% to 25% (Figure 4b).The frequency of red-attack trees within plots is shown in Figure 4c.Red-attack tree frequencies of 0-5%, 5-10%, and 10-15% account for 72.5%, 20%, and 7.5% of plots, respectively.According to field plots, we found that a tree mortality of >25% and a red-attack tree percentage of >10% did not occupy a large proportion of the plots.The results showed that using percentages of tree mortality and red-attack trees to reflect damage levels was not suitable for Tomicus spp.because it neglects shoot damage.On the other hand, healthy forest (0 ≤ SDR ≤ 10%), slightly to moderately infested forest (10 ≤ SDR ≤ 50%), and severely infested forest (50% < SDR) account for 25%, 35%, and 40% of plots, respectively (Figure 4a).Based on the classification standard (Standard of Forest Pests Occurrence and Disaster LY/T 1681-2006), we found that shoot damage ratios could clearly classify the damage levels in these field plots.The results showed that for Tomicus spp., SDR was better than tree mortality or red-attack trees to reflect the damage degrees of forests.

Classification of Damage Degree
The spectral anomaly of MSI resulted in the highest classification performance among all the spectral anomalies analyzed (i.e., RGI, NDVI, NDMI, and MSI) (Table 5).The classification accuracy of the MSI threshold value is shown in Figure 3. Using the threshold value of 0 for MSI, we distinguished slightly to moderately infested forests from healthy forests.Using the threshold value of 0.13 for MSI, we distinguished slightly to moderately infested forests from severely infested forests.The overall classification accuracy for the 2016 Landsat classification is 86.38% (Kappa = 0.80) with

Classification of Damage Degree
The spectral anomaly of MSI resulted in the highest classification performance among all the spectral anomalies analyzed (i.e., RGI, NDVI, NDMI, and MSI) (Table 5).The classification accuracy of the MSI threshold value is shown in Figure 3. Using the threshold value of 0 for MSI, we distinguished slightly to moderately infested forests from healthy forests.Using the threshold value of 0.13 for MSI, we distinguished slightly to moderately infested forests from severely infested forests.The overall classification accuracy for the 2016 Landsat classification is 86.38% (Kappa = 0.80) with accuracies across all classes (Table 6).The healthy forest class has the highest classification accuracies (e.g., producer's accuracy was 96.7%) while the greatest confusion existed between the slightly to moderately infested forest and the severely infested forest.The omission error for the severely infested forest class multi-date classification was 13.3% and the commission error was 17.5%.The omission error for the slightly to moderately infested forest class multi-date classification was 23.9% and the commission error was 18.0%.As a result, the classification using MSI was used in further analyses.

Outbreak Time and Spread Path Estimation
The frequency distribution of severely infested pixels within 30 m resolution of pixels generated from MSI threshold values is shown in Figure 5.The number of severely infested pixels show an abrupt jump from 2012 to 2013, which indicates a beetle outbreak in 2013.Based on the records of the local forestry station, the Yunnan pine suffered from serious shoot beetle damage in 2013.The outbreak area increased and reached its peak in 2014, then subsequently declined.
Forests 2018, 9, 39 8 of 14 accuracies across all classes (Table 6).The healthy forest class has the highest classification accuracies (e.g., producer's accuracy was 96.7%) while the greatest confusion existed between the slightly to moderately infested forest and the severely infested forest.The omission error for the severely infested forest class multi-date classification was 13.3% and the commission error was 17.5%.The omission error for the slightly to moderately infested forest class multi-date classification was 23.9% and the commission error was 18.0%.As a result, the classification using MSI was used in further analyses.

Outbreak Time and Spread Path Estimation
The frequency distribution of severely infested pixels within 30 m resolution of pixels generated from MSI threshold values is shown in Figure 5.The number of severely infested pixels show an abrupt jump from 2012 to 2013, which indicates a beetle outbreak in 2013.Based on the records of the local forestry station, the Yunnan pine suffered from serious shoot beetle damage in 2013.The outbreak area increased and reached its peak in 2014, then subsequently declined.Beetle damage maps based on MSI threshold values from 2011 to 2016 are shown in Figure 6a.The movement of the outbreak center is shown in Figure 6b.We found that from 2011 to 2016, the centers of severely infested forests have a significant trend going from a southwest to northeast direction.Furthermore, the southwest wind direction is the dominant wind direction in Pupeng Town (Figure 7), and it seems that wind direction may affect the moving direction of the beetles.
Forests 2018, 9, 39 9 of 14 Forests 2018, 9, 39 9 of 14 centers of severely infested forests have a significant trend going from a southwest to northeast direction.Furthermore, the southwest wind direction is the dominant wind direction in Pupeng Town (Figure 7), and it seems that wind direction may affect the moving direction of the beetles.

Discussion
Classification using MSI produced the highest accuracy for the multi-date approach.MSI is sensitive to water content and, thus, drought.As such, it is highly related to the weakness or vulnerability of the Yunnan pine.This can explain the good performance of MSI.However, the accuracy varies with SDR.High classification accuracies were obtained for high SDRs within pixels, however, there was low accuracy for lower SDRs.Our findings indicate the usefulness of multi-date Landsat imagery for detecting severe beetle disturbance (i.e., SDR > 50% within a pixel).However, there is still difficulty in accurately differentiating slightly to moderately infested forests (i.e., 10% ≤ SDR ≤ 50% within a pixel).
SDR is a good indicator for detecting forest disturbances caused by the beetles.Meddens et al. [45] found that Landsat is useful for detecting severe tree mortality (i.e., death of more than approximately 25% of trees in the canopy within a pixel) but cannot detect slowly progressing disturbances within a pixel (i.e., as fewer than 25% red trees at any one time) with high accuracy.In our 40 permanent field plots, tree mortality (gray-attack trees) of more than 25% in the canopy within a pixel account for 10% of plots, and red-attack trees of >10% within a pixel also account for less than 10% of plots.Thus, if we use gray-attack or red-attack tree percentages to reflect the damage degrees of forests, the detection accuracies would be low.According to the research of Meddens et al. [45] and the special damage characteristics of the beetles, we chose the SDR to define the damage degrees instead of tree mortality.
Our methods use a time window from November to December to detect beetle disturbance.This is different from previous studies which use the middle of the growing season (July or August) as a Forests 2018, 9, 39 9 of 14 centers of severely infested forests have a significant trend going from a southwest to northeast direction.Furthermore, the southwest wind direction is the dominant wind direction in Pupeng Town (Figure 7), and it seems that wind direction may affect the moving direction of the beetles.

Discussion
Classification using MSI produced the highest accuracy for the multi-date approach.MSI is sensitive to water content and, thus, drought.As such, it is highly related to the weakness or vulnerability of the Yunnan pine.This can explain the good performance of MSI.However, the accuracy varies with SDR.High classification accuracies were obtained for high SDRs within pixels, however, there was low accuracy for lower SDRs.Our findings indicate the usefulness of multi-date Landsat imagery for detecting severe beetle disturbance (i.e., SDR > 50% within a pixel).However, there is still difficulty in accurately differentiating slightly to moderately infested forests (i.e., 10% ≤ SDR ≤ 50% within a pixel).
SDR is a good indicator for detecting forest disturbances caused by the beetles.Meddens et al. [45] found that Landsat is useful for detecting severe tree mortality (i.e., death of more than approximately 25% of trees in the canopy within a pixel) but cannot detect slowly progressing disturbances within a pixel (i.e., as fewer than 25% red trees at any one time) with high accuracy.In our 40 permanent field plots, tree mortality (gray-attack trees) of more than 25% in the canopy within a pixel account for 10% of plots, and red-attack trees of >10% within a pixel also account for less than 10% of plots.Thus, if we use gray-attack or red-attack tree percentages to reflect the damage degrees of forests, the detection accuracies would be low.According to the research of Meddens et al. [45] and the special damage characteristics of the beetles, we chose the SDR to define the damage degrees instead of tree mortality.
Our methods use a time window from November to December to detect beetle disturbance.This is different from previous studies which use the middle of the growing season (July or August) as a

Discussion
Classification using MSI produced the highest accuracy for the multi-date approach.MSI is sensitive to water content and, thus, drought.As such, it is highly related to the weakness or vulnerability of the Yunnan pine.This can explain the good performance of MSI.However, the accuracy varies with SDR.High classification accuracies were obtained for high SDRs within pixels, however, there was low accuracy for lower SDRs.Our findings indicate the usefulness of multi-date Landsat imagery for detecting severe beetle disturbance (i.e., SDR > 50% within a pixel).However, there is still difficulty in accurately differentiating slightly to moderately infested forests (i.e., 10% ≤ SDR ≤ 50% within a pixel).
SDR is a good indicator for detecting forest disturbances caused by the beetles.Meddens et al. [45] found that Landsat is useful for detecting severe tree mortality (i.e., death of more than approximately 25% of trees in the canopy within a pixel) but cannot detect slowly progressing disturbances within a pixel (i.e., as fewer than 25% red trees at any one time) with high accuracy.In our 40 permanent field plots, tree mortality (gray-attack trees) of more than 25% in the canopy within a pixel account for 10% of plots, and red-attack trees of >10% within a pixel also account for less than 10% of plots.Thus, if we use gray-attack or red-attack tree percentages to reflect the damage degrees of forests, the detection accuracies would be low.According to the research of Meddens et al. [45] and the special damage characteristics of the beetles, we chose the SDR to define the damage degrees instead of tree mortality.
Our methods use a time window from November to December to detect beetle disturbance.This is different from previous studies which use the middle of the growing season (July or August) as a time window [45,56,57].A major reason for this is (1) the two stages of the life cycle for shoot beetles.They attack shoots from May to November, thus, the number of damaged shoots enters a rising trend in this time periods.The amount of shoot damage would reach its maximum peak sometime in November to December.After that, these beetles begin to transfer from the shoot to the trunks of trees.The color of the damage shoots would turn from green to red in November-December.This largely influences spectral reflectance, making it easier for RS to detect.Another major reason is (2) climate.
Our study area has a clear dry and wet season.The rainy season is from June to October, which is within the growing season.The significant cloud cover during July and August is not good for optical remote sensing.
Lacking historical monitoring data, it was difficult to study the previous and current state of the forests.In this study, we utilized multi-date Landsat imagery and successfully predicted the outbreak year (i.e., 2013) and reconstructed the locations of severely infested forests from 2011 to 2016.Although we could not truly evaluate the accuracy of beetle-disturbed regions before 2016, our study could detect the beetle disturbance to some degree, especially in severely infested forests.This method will be useful for similar regions without historical records.
Another interesting finding is that shoot beetles moved from the southwest to northeast direction from 2011 to 2016.Based on previous research, MPBs attack surrounding injured or diseased trees stressed by moisture limitation or disturbances during the period in that spot growth or areas within 2 km [58][59][60].Additionally, MPBs are transported long distances by the wind [61].Because T. yunnanensis and T. minor transfer between shoots during the shoot-feeding time period, their time exposed to air would be longer than MPBs.Ye et al. [39] found that the spread of T. piniperda was also greatly influenced by the wind.The number of beetles downwind of the bole mass of the attacked tree were greater than that upwind.This shows that wind is a vector for spread.In our study, the movement of the center of the severely infested forests from 2011 to 2016 was more than 2 km.Combining our findings with previous research, we speculate wind to be a vector for long-distance dispersal of Tomicus spp.
The rising and decreasing trends of severely infested forests from 2011 to 2016 appears strange.The major reasons for this include climate and human factors.Since 2009, widespread drought (i.e., the average annual precipitation from 2010 to 2014 was 810 mm) in Yunnan Province weakened trees, making it suitable for beetles to attack trunks.However, the precipitation rose after 2014 (i.e., the average annual precipitation from 2015 to 2016 was 1085 mm), which increased the resistance of Yunnan pine to beetles.This could partially explain the declining trend of pixels of severely infested forest from 2014 to 2016.Harvesting is the main human activity affecting this trend.During the investigation time, we found that villagers harvested infested trees to use as fire-wood, especially the red-attack trees which contain numerous adult beetles.This behavior reduced beetle density and the area of severely infested forests.However, further research needs to determine the other specific drivers of beetle spread besides wind.

Conclusions
Multi-data derived from time series data such as Landsat provide rich information for detecting forest disturbances such as T. yunnanensis and T. minor outbreaks in Yunnan Province.From our study, we can conclude two major points: (1) The SDR was a better measure than tree mortality for monitoring shoot beetle damage.
(2) Reconstructed SDR can be used to predict outbreak time and beetle moving patterns.
These results are important for the forest industry to implement forestry operations that could reduce the economic and ecological impacts of beetle attacks.Therefore, we discussed the relationship between the change trends of the severely infested forests and the precipitation from 2011 to 2016.These suggest that the beetle outbreak may be connected to the drought that occurred in the area during 2009.In future research, we will study the response of the beetle outbreaks to climate change.

Figure 1 .
Figure 1.(a) Yunnan Province in China; (b) The 2016 Landsat image with the study area in the right corner of the image; (c) Study area location.

Figure 1 .
Figure 1.(a) Yunnan Province in China; (b) The 2016 Landsat image with the study area in the right corner of the image; (c) Study area location.

Figure 2 .
Figure 2. Examples of temporal pixel trajectories representing MSI anomalies from the healthy forests, slightly to moderately infested forests, and severely infested forests.Dash straight line represents the threshold value (MSI anomalies = 0.13).

Figure 3 .
Figure 3. Overall accuracy, healthy forest class accuracy, slightly to moderately infested forest class accuracy, and severely infested forest class accuracy from confusion matrices with different MSI threshold values used to (a) separate the healthy class from the slightly to moderately infested class and to (b) separate the slightly to moderately infested forest class from the severely infested class.

Figure 2 .
Figure 2. Examples of temporal pixel trajectories representing MSI anomalies from the healthy forests, slightly to moderately infested forests, and severely infested forests.Dash straight line represents the threshold value (MSI anomalies = 0.13).

Figure 2 .
Figure 2. Examples of temporal pixel trajectories representing MSI anomalies from the healthy forests, slightly to moderately infested forests, and severely infested forests.Dash straight line represents the threshold value (MSI anomalies = 0.13).

Figure 3 .
Figure 3. Overall accuracy, healthy forest class accuracy, slightly to moderately infested forest class accuracy, and severely infested forest class accuracy from confusion matrices with different MSI threshold values used to (a) separate the healthy class from the slightly to moderately infested class and to (b) separate the slightly to moderately infested forest class from the severely infested class.

Figure 3 .
Figure 3. Overall accuracy, healthy forest class accuracy, slightly to moderately infested forest class accuracy, and severely infested forest class accuracy from confusion matrices with different MSI threshold values used to (a) separate the healthy class from the slightly to moderately infested class and to (b) separate the slightly to moderately infested forest class from the severely infested class.

Figure 4 .
Figure 4. (a) The frequency of shoot damage ratios within plots; (b) The frequency of gray-attack trees within plots; (c) The frequency of red-attack trees within plots.

Figure 4 .
Figure 4. (a) The frequency of shoot damage ratios within plots; (b) The frequency of gray-attack trees within plots; (c) The frequency of red-attack trees within plots.

Figure 5 .
Figure 5. Frequency distribution of severely infested pixels (30 m resolution) generated from MSI threshold values.Beetle damage maps based on MSI threshold values from 2011 to 2016 are shown in Figure 6a.The movement of the outbreak center is shown in Figure 6b.We found that from 2011 to 2016, the

Figure 5 .
Figure 5. Frequency distribution of severely infested pixels (30 m resolution) generated from MSI threshold values.

Figure 6 .
Figure 6.(a) The damage locations from 2011 to 2016 were speculated using the MSI threshold values (green corresponds to healthy forests, yellow corresponds to slightly to moderately infested forests, and red corresponds to severely infested forests); (b) The median center from 2011 to 2016 was calculated using ARCGIS.

Figure 7 .
Figure 7.The percentage of wind directions from 1 January 2011 to 31 May 2017.

Figure 6 .
Figure 6.(a) The damage locations from 2011 to 2016 were speculated using the MSI threshold values (green corresponds to healthy forests, yellow corresponds to slightly to moderately infested forests, and red corresponds to severely infested forests); (b) The median center from 2011 to 2016 was calculated using ARCGIS.

Figure 6 .
Figure 6.(a) The damage locations from 2011 to 2016 were speculated using the MSI threshold values (green corresponds to healthy forests, yellow corresponds to slightly to moderately infested forests, and red corresponds to severely infested forests); (b) The median center from 2011 to 2016 was calculated using ARCGIS.

Figure 7 .
Figure 7.The percentage of wind directions from 1 January 2011 to 31 May 2017.

Figure 7 .
Figure 7.The percentage of wind directions from 1 January 2011 to 31 May 2017.

Table 2 .
Dates, sensors, and cloud covers of Landsat images from 2004 to 2016 used in this study.
a Priority image for that year.Other images in same year were used to fill cloudy areas and SLC (Scan Lines Corrector)-off areas in the priority image.TM5: The matic Mapper 5; ETM+: Enhanced Thematic Mapper; OLI: Operational Land Imager.

Table 4 .
Landsat-derived spectral indices considered in the classification.

Table 5 .
Overall accuracy for different spectral anomalies used in the 2016 evaluation of Landsat multi-date classifications.

Table 6 .
Confusion matrix of the 2016 multi-data image classification using MSI indices and pixel locations from aggregated classified aerial imagery with >70% class proportions.

Table 5 .
Overall accuracy for different spectral anomalies used in the 2016 evaluation of Landsat multi-date classifications.

Table 6 .
Confusion matrix of the 2016 multi-data image classification using MSI indices and pixel locations from aggregated classified aerial imagery with >70% class proportions.