Assessment of Poplar Looper (Apocheima cinerarius Erschoff) Infestation on Euphrates (Populus euphratica) Using Time-Series MODIS NDVI Data Based on the Wavelet Transform and Discriminant Analysis

Abstract

Poplar looper (Apocheima cinerarius Erschoff) is a destructive insect infesting Euphrates or desert poplars (Populus euphratica) in **njiang, China. Since the late 1950s, it has been plaguing desert poplars in the Tarim Basin in **njiang and caused widespread damages. This paper presents an approach to the detection of poplar looper infestations on desert poplars and the assessment of the severity of the infestations using time-series MODIS NDVI data via the wavelet transform and discriminant analysis, using the middle and lower reaches of the Yerqiang River as a case study. We first applied the wavelet transform to the NDVI time series data in the period of 2009–2014 for the study area, which decomposed the data into a representation that shows detailed NDVI changes and trends as a function of time. This representation captures both intra- and inter-annual changes in the data, some of which characterise transient events. The decomposed components were then used to filter out details of the changes to create a smoothed NDVI time series that represent the phenology of healthy desert poplars. Next the subset of the original NDVI time series spanning the time period when the pest was active was extracted and added to the smoothed time series to generate a blended time series. The wavelet transform was applied again to decompose the blended time series to enhance and identify the changes in the data that may represent the signals of the pest infestations. Based on the amplitude of the enhanced pest infestation signals, a predictive model was developed via discriminant analysis to detect the pest infestation and assess its severity. The predictive model achieved a severity classification accuracy of 91.7% and 94.37% accuracy in detecting the time of the outbreak. The methodology presented in this paper provides a fast, precise, and practical method for monitoring pest outbreak in dense desert poplar forests, which can be used to support the surveillance and control of poplar looper infestations on desert poplars. It is of great significance to the conservation of the desert ecological environment.

1. Introduction

Euphrates or desert poplars (Populus euphratica) are medium-sized deciduous trees mainly found in fluvial and floodplain areas in arid and semi-arid regions, where the groundwater is close to the surface. Their unique physiological traits make them tolerant of saline and brackish water and resistant to droughts and sand storms. They form the near top-level natural arbor community in the desert environment and are a significant component of Tugai floodplain ecosystems in arid climates. Desert poplars play an important role in maintaining ecological balance, establishing windbreaks for the protection from wind and the encroachment of sands, fixing shifting sands, regulating the climate and improving the health of ecosystems [1,2]. However, poplar plantations in China have been constantly threatened by defoliating insects, particularly the poplar looper Apocheima cinerarius Erschoff (Lepidoptera: Geometridae). In recent years, there are widespread outbreaks of insect pests in desert poplar plantations in northwest China due to the weakened resistance of water-stressed poplars to pests and diseases resulted from excessive unregulated water diversion and climate change. Poplar looper is the most serious insect pest affecting desert poplars in the Tarim Basin, ** strategies to protect desert poplars.

Poplar looper is a geometrid moth, distributed in the arid climate zone in the East Asia monsoon region, including the entire area of northern China [5]. A lot of studies have been done on this moth species, including its habitat, physiology, life cycle, evolutionary history, geographical distribution, population dynamics, the nature of its attack and damage it caused, as well as forecasts and control of its outbreaks [6,7,8,9,10,11,12]. However, Euphrates are distributed in deserts, which are usually remote and poorly accessible by transport. It is difficult to carry out large-scale and in-depth studies on the pest infestation using traditional field survey techniques. So far, there has been a lack of information on the extent and damage level of the poplar looper infestation in Euphrates forests in China.

Remote sensing is a proven technology for pest monitoring and damage assessment over large areas. Many studies have used remotely sensed data with various spatial resolutions to detect and assess forest insect infestations, based on the analysis of spectral responses of trees to the biophysical phenomena such as water stress caused by insect attack [13,14,15,16,17,18,19]. Over the years, a number of spectral indices (such as the Normalised Difference Vegetation Index or NDVI, the Normalised Burn Ratio or NBR, the Moisture Stress Index or MSI, the Enhanced Wetness Difference Index or EWDI, and the Leave Area Index or LAI) have been used for the detection and differentiation of insects and diseases using multispectral and hyperspectral data with different spatial resolutions. According to a recent review by Senf et al. [20], medium and coarse spatial resolution remote sensing data and NDVI were most often used for map** insect induced broadleaved defoliation, and intra-annual time series analysis was frequently applied in the analysis of these data sets. As many insects are active for a short period when an outbreak is detectable, most studies used very dense time series of remote sensing data such as the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), and SPOT VEGETATION to achieve daily temporal resolution. MODIS data have both high spectral and temporal resolutions, and have been more frequently employed [21]. For example, they have been used to detect forest defoliation induced by the European gypsy moth in North America [22], the autumnal moth and winter moth in northern Sweden [23], the European pine sawfly in southeastern Norway [24], and by the bark beetle in North America [25]. These studies achieved various levels of success. For example, Olsson et al. [26] reported that with MODIS data only 50% of damaged Scots pine forest stands caused by pine sawfly in eastern Finland were detected and 22% of healthy stands were misclassified, while in northern Sweden 75% of damaged birch forest stands were detected and 19% of healthy birch forest stands were misclassified. However, these applications have shown potential of MODIS data for large area monitoring and early detection [27]. Our research used MODIS NDVI data to detect defoliation of poplar trees induced by poplar looper and assess the severity of the insect outbreak.

In recent years, for large-scale pest infestation monitoring, Landsat and Sentinel-2 data with a higher spatial resolution have been widely used [28,29,30,31]. Since poplar forests are most vulnerable to poplar looper attacks in less than two months [32], it requires monitoring data to be updated daily to quickly monitor the status of poplar forests damaged by the insect. Although the spatial resolution of MODIS data is coarser than Landsat and sentinel-2, its long timespan and high temporal resolution are more suitable for this monitoring task. With the temporal resolution of 16 days, Landsat images could miss the foliation peak starting of defoliation, reducing the information about the pest infestation. Though Sentinel-2 data have higher spatial and temporal resolutions, hence a higher monitoring accuracy than Landsat 8 [33,34], the data with a 5-day temporal resolution only started to become available in 2017, which limits their use in the analysis and modelling of the long-term trend. Although the spatial resolution of MODIS data is coarser than Landsat and Sentinel-2, its long timespan and high temporal resolution are more suitable for the monitoring of poplar looper attacks on poplar forests.

Time series fitting techniques have been used to analyse MODIS time series data to detect insect disturbance [20]. For example, TIMESAT, a computer program for time series analysis of satellite sensor data, fits smoothed functions (including Savitzky-Golay fitering, least-squares fitted asymmetric Gaussian functions and double logistic functions) to time series data and extracts seasonality parameters [35], which was applied in the studies by Eklundh et al. [24], Anees and Aryal [25] and Olsson et al. [23,26]. Jia [36], Huang [37] and Wang [38] used the S-G method to filter the MODIS NDVI time series and combined it with the data on the poplar phenology and poplar looper life history to extract the poplar looper hazard information in a poplar forest area using a change detection method. The wavelet transform is a signal processing tool for analysing both the local frequency and temporal behaviour of signals by decomposing a signal into frequency components and then studying each component with a resolution matched to its scale [39]. Wavelet analysis has the merit of decomposing the signal into components of different scales (or frequencies) thus facilitating the detection of subtle signals. An NDVI time series can be considered as a signal. The wavelet transform was used in this study to decompose the NDVI time series into components of different scales, facilitating the detection of trends and subtle signals. Discriminant analysis is a statistical technique for analysing a data set with individuals classified into particular groups, and then using the results to classify new individuals that are not included in the above data set [40].

In this paper, we present an innovative and effective method for detecting and monitoring the poplar looper infestation of dense desert poplar forests using MODIS time series data through the wavelet transform and discriminant analysis. The objectives of this study can be formulated as the following:

(1)

Proposing a method of decomposition of mixed signals from the NDVI time series to isolate and enhance the pest infestation signals; and

(2)

Building a predictive model for identifying severity levels and dates of pest outbreaks.

2. Materials and Methods

2.1. Study Area

The study area is located in **, land surface reflectance, atmospheric interference, and influence of other types of noise.

2.5. Methodology

NDVI time series are non-stationary and usually characterised by long term trends, seasonal variations, and localised abrupt changes resulting from disturbance events (such as insect outbreaks, droughts and sand storms) as shown in Figure 3. They allow for the analysis of the inter-annual and intra-annual changes in the vegetation canopy. Smoothing functions and curve-fitting algorithms have been used to analyse the variability of vegetation using MODIS time series [35]. Their limitations were discussed in Galford et al. [44]. In this study, we adopted the wavelet transform to analyse the MODIS NDVI time series in order to extract information about poplar looper infestation. Our approach involved using the wavelet transform to:

(1)

decompose the NDVI time series into approximation (low-frequency) components representing trends and detail (high-frequency) components representing localised changes or details at multiple temporal scales;

(2)

de-noise the NDVI time series using the decomposed time series components from 1) to remove all sorts of noises and create a new NDVI time series representing the intra- and inter-annual variations in the NDVI for healthy poplar tree canopy;

(3)

compare and match the field observed leaf phenology of desert poplars with the de-noised or smoothed time series derived from (2) to identify the overall trend, seasonality and periodicity of the NDVI of poplar tree canopy;

(4)

compare and match the field observed poplar looper infestation dynamics with the decomposed time series components derived from (1) and information from (3) to extract signals of defoliation (abrupt changes in the NDVI) caused by the pest infestation;

(5)

reconstruct the NDVI time series by blending the extracted insect infestation signals from (4) to the de-noised NDVI time series from (2), which preserves the signals of the insect infestation, but removes other noises;

(6)

decompose the reconstructed NDVI time series from (5) to enhance the signals of the pest infestation;

(7)

build a predictive model via discriminant analysis using the information generated from (6) to detect pest infestation and classify the level of severity of infestation;

(8)

assess the detection and severity classification accuracy.

Figure 4 shows the process.

2.5.1. Wavelet Transforms

Basically, the wavelet transform is to decompose a time series into different components to determine various modes of variability and how those modes change in time, which allows for a full representation of the localised and transient events occurring at different temporal scales [45,46,47,48,49,50,51]. It transforms a time series using a set of basis functions, called wavelets. There are two types of wavelet transforms: continuous wavelet transform (CWT) and discrete wavelet transform (DWT).

The DWT provides an efficient alternative. It can be considered as a sequence of pairs of high-pass and low-pass filters, known as a filter bank, applied to a time series. Low-pass filters remove details or high-frequency components, smoothing the data series and producing approximation coefficients, whereas high-pass filters enhance details, producing detail coefficients. The DWT transforms a time series by sub-band coding where the series is passed through filters with different cutoff frequencies at different scales, and typically the half-band filters down-sample the series by a factor of two at each level of decomposition, which is called dyadic sampling [52]. Mathematically, the basis function used in the DWT is defined as [53]:

$${\psi }_{j,k}(t)=\frac{1}{\sqrt{a_0^j}}\psi (\frac{t-a_0^j{b}_0^j}{a_0^j})$$

(1)

where a₀ and b₀ are constants, j is the decomposition level or scale, and k is the time translation factor. For dyadic sampling of DWT, a₀ = 2 and b₀ = 1; and Level j approximation coefficients are decomposed into Level j + 1 approximation coefficients and Level j + 1 detail coefficients with the same length. The information of the original time series is retained in wavelet coefficients.

The original time series can be reconstructed from the approximation and detail coefficients at each level by up-sampling, passing through high- and low-pass synthesis filters and adding them. This process is also called multi-resolution analysis (MRA). Martinez and Gilabert [48] offered a detailed discussion on MRA. NDVI time series are discrete in nature. Therefore, we applied MRA with the dyadic DWT to the analysis of poplar canopy NDVI time series.

The choices of the mother wavelet $\psi (t)$, the number of vanishing moments and the decomposition level are three key issues in applying a DWT to the NDVI time series. As a wavelet used for the DWT must meet the regularity condition, i.e.,

$${\displaystyle {\int }_R{\mathrm{t}}^{\kappa }}\psi (t)dt=0\hspace{1em}\hspace{1em}\kappa =1,\dots ,K$$

(2)

where K is the number or order of vanishing moments, we used the Harr, Daubechies and Meyer wavelets to conduct a number of experiments. The Daubechies wavelet produced the best outcome and was selected for the DWT in this study.

The vanishing moment is a criterion about how a wavelet function decays toward infinity. When the wavelet has K-order vanishing moments, the wavelet coefficients of all polynomials are zero, that is, the high-frequency component of the wavelet decomposition is 0 [54]. A Daubechies wavelet with K vanishing moments should have the minimum filter size of 2K. In order to quickly detect the sudden changes, a higher number of vanishing moments is required. We selected six vanishing moments after several experiments when applying the Daubechies wavelet to our NDVI time series.

The highest decomposition level can be calculated as $\left[lo{g}_{2}n\right]$ where $\left[•\right]$ denotes taking the integer part of the real value in the bracket and n is the number of observations or time points in the time series [39]. In our study, the time series have 138 time points, therefore, we can decompose the time series up to eight levels. As the decomposition level increases, more subsets of the time series and detailed information contained in the series at a larger temporal scale would appear.

More information could contribute to better performance, but it may reduce the computing efficiency. Therefore, choosing a suitable decomposition level in DWT is of paramount importance to its performance and efficiency. We conducted experiments with the decomposition levels from one to eight, and found that eight is the optimum number of decomposition levels for the poplar canopy NDVI time series data as it produced the smallest standard error, and reserved key information and removed significant amount of noise in the time series.

Essentially, the DWT generates separate time series or components from the original one. The original information will be distributed into these different components in the form of wavelet coefficients. Therefore, it allows for reconstruction. Component x_j(t) at the decomposition level j can be reconstructed by the following equation [53]:

$$x_j(t)={\displaystyle {\sum }_{\kappa }{W}_x(j,k)}\psi (\frac{t-a_0^j{b}_0^j}{a_0^j})$$

(3)

The wavelet transform was implemented in MATLAB 2015b.

2.5.2. De-Noising of the NDVI Time Series

There are many techniques available for filtering noises in time series data, including the mean-value iteration filters (MVI), Savitzky-Golay filters, Fourier analysis, Harmonic Analysis of Time Series (HANTS), and wavelet transform [55,56,57,58,59]. We created an NDVI time series that represents the phenology of healthy desert poplars by removing noise (including the insect disturbance) in the original NDVI time series via wavelet de-noising. Wavelet de-noising is to remove the low amplitude noise from a time series [53] It involves first decomposing the original NDVI time series, then selecting an appropriate threshold limit at each decomposition level and removing all wavelet coefficients that are smaller than the threshold limit by setting them to zero, and finally reconstructing the time series from thresholded wavelet coefficients using Equation (3). Figure 5. shows the noise-filtered NDVI time series compared with the original one shown in Figure 3, in which the noise has been almost entirely suppressed and can be considered as the normal growth curve of desert poplars free of disturbance caused by insects, forest fires, droughts, harvesting, clouds and snows.

2.5.3. Detection of the Pest Infestation Signals

As discussed in Section 2.2, the poplar forest disturbance by poplar looper occurred in March and April in the study area. According to our field observations, the sudden changes in the original NDVI time series in the two months were mainly caused by poplar looper defoliation. When the pest infestation occurred, the NDVI value dropped sharply. At the end of the outbreak, the second foliation occurred [60] and, as a result, the NDVI value increased again. In order to identify exactly when pest infestations happened and assess how severe they were, we replaced the NDVI values at the time points in March and April each year from 2009 to 2014 in the noise-filtered NDVI time series with their original values from Figure 3 and created a blended NDVI time series. Figure 6 shows the blended time series created based on the smoothed time series in Figure 5. The blended time series contains only noise made by the pest disturbance.

The blended NDVI time series was then decomposed into eight levels using the DWT to accentuate the abrupt changes representing the pest disturbance in the time series. By analysing the detail components, we detected the time points of the insect disturbance.

2.5.4. Discriminant Analysis and Accuracy Assessment

Discriminant analysis is to determine class membership using one or several discriminant functions based on linear combinations of the predictor variables that provide the best discrimination between the classes [61]. A general form of a discriminant function can be written as

$$y=b_0+b_1{x}_1+b_2{x}_2+\cdots +b_n{x}_n$$

(4)

where y is the categorical dependent variable or grou** variable used for classifying into two or more groups, x_i (i = 1, 2, …, n) is the predictor variable, and b_i (i = 0, 1, 2, …, n) is the discriminant function coefficient. The coefficients are derived from a sample of cases for which class membership is known. Once the function is generated, it can be applied to new cases of unknown classes that have measurements for the predictor variables. In our case, we only have one discriminant function where the y value is a severity discriminant score, and one predictor variable x is the amplitude of the changes in the NDVI between two time points extracted from detail components derived from the blended NDVI time series. There are several discriminant rules for predicting class membership based on discriminant scores, including the maximum likelihood rule, Bayes discriminant rule and Fisher’s linear discriminant rule [62]. We used the Bayes discriminant rule. SPSS was used to conduct the discriminant analysis.

The predictive model built through discriminant analysis was used for severity classification of the pest disturbance. The prediction accuracy of the model was assessed using the Press Q statistic. The Press Q is a statistical test for the discriminatory power of the classification matrix resulted from discriminant analysis when compared with the chance model. It is expressed as:

$$Press Q=\frac{{\left[N-lL\right]}^2}{N(L-1)}$$

(5)

where N is the total sample size, l is the number of correctly classified cases, and L is the number of classes. The critical value of Press Q is 6.63 at the 0.01 significance level, and 3.84 at the 0.05 significance level (subject to a chi-square distribution with a degree of freedom of 1).

3. Results

3.1. Trend and Seasonality of Poplar Canopy NDVI

The poplar canopy NDVI time series at each sampling site were analysed to understand the trend and seasonality of poplar canopy NDVI. For example, the NDVI time series curve shown in Figure 3 indicates that at Sampling Site 2, the maximum NDVI is 0.53, the minimum NDVI is 0.03 and the average is 0.23.

The NDVI time series was decomposed into eight levels through the DWT. Figure 7 shows eight reconstructed detail components (D1 to D8) and the level eight approximation component (A8) derived from the blended time series shown in Figure 6. Lower detail levels have higher frequencies, which represent the rapidly changing components of the dataset, whereas the higher detail levels have lower frequencies, which represent the slowly changing component of the dataset. The approximation component (A8) represents the slowest changing component of the dataset. It depicts the long-term trend of the poplar canopy NDVI, which exhibits a stable upward trend from 2009 to 2012, and then a slight decreasing trend after 2012. It indicates that the climate condition was generally improved to favour the growth of desert poplars, leading to the rise of NDVI. According to Liu et al. [63] and Liu et al. [64], in the Yerqiang River Catchment, the climate has been clearly becoming warmer and more humid in recent years, which promotes the plant growth. The trend of the poplar canopy NDVI coincides with the trend of the climate change. Therefore, the trend component of the NDVI time series reflects the climate change in the region to some extent.

The noise-filtered NDVI time series as shown in Figure 5 presents the seasonality of the NDVI time series. The growing season of poplars is observed every year. The maximum NDVI values occurred in the period of mid-June to mid-July every year, while the minimum values appeared in every late December to early February in the following year. This temporal pattern corresponds to the phenological cycles of poplar trees in the middle and lower reaches of the Yerqiang River.

3.2. Poplar Forest Defoliation Caused by Poplar Loopers

The detail component D1 shown in Figure 7b represents 32-day periodicity (two time points) of the original time series. It contains detail features that clearly describe the sudden changes and peaks or spikes in the observed NDVI time series. These sudden changes had various peak amplitudes and occurred at different time points. Those with a large peak amplitude can be considered as “resonance” points at which a sudden change in the NDVI was driven by an external force. By analysing these resonance points, we found that the sudden changes in the NDVI at the twotime points 5, and 122 appeared corresponding to the seasons when poplar loopers hatched and fed on the buds and leaves of poplars, therefore may be caused by defoliation due to the pest. These time points areearly March of 2009 and early April of 2014 respectively.

Generally, the average temperature in March determined when the resonance points occurred. According to the meteorological records available in the study area from 2009 to 2014, the average temperature in early March in 2014 was relatively low as shown in Figure 8, which delayed the poplar’s leaf unfolding and the poplar looper’s hatching. The growing season of desert poplars usually began at the Mid-February and ended in late November. The pupae of poplar looper hatched in early March. Mid-March to April is the peak period of its outbreak. But the temperature in late February 2009 was higher than usual, which resulted in earlier hatching of the pupae of poplar looper, thus an earlier outbreak of the insect. By examining the original NDVI time series and the detail component D1 in Figure 7b, we found that the NDVI value of poplar canopy increased rapidly when the growth period started each year, and then decreased sharply in early March of 2009, and in early April of 2014. According to our field observations, these sharp drops in the NDVI were all a direct result of poplar looper caused defoliation. From the detail component D1, the peak amplitudes of the sudden changes in the NDVI at the resonance points were 0.0269 in early March of 2009 and 0.0473 in early April of 2014. Larger the drop in the NDVI, severer the infestation. The NDVI started to increase in April after the insect-induced defoliation due to the second foliation.

3.3. Prediction Model of Infestation Severity

The discriminant analysis produced four classification functions to classify the infestation severity of the insect infestation, which are

$$y\left(0\right)=-3.145+366.557x$$

(6)

$$y\left(1\right)=-10.734+845.094x$$

(7)

$$y\left(2\right)=-34.494+1590.440x$$

(8)

$$y\left(3\right)=-59.798+212.522x$$

(9)

where $x$ is the amplitude of the NDVI change from the detail component D1 in Figure 7b and $y\left(0\right)$, $y\left(1\right)$, $y\left(2\right)$ and $y\left(3\right)$ are the discriminant scores of severity classification of the insect infestation. Substituting the existing x value into Equations (6)–(9) for calculation, the one which gets the highest score will be classified into that category. Thirty-six sample cases of 2014 were used to build the discriminant functions.

Table 2. Class statistics.

Classification	Average	Standard Deviation	Valid N (list Status)
			Unweighted	Weighted
0	0.0096	0.0052	22	22.00
1	0.0221	0.0052	9	9.00
2	0.0416	0.0053	3	3.00
3	0.0553	0.0000	2	2.00
Total	0.0179	0.0141	36	36.00

lists the class statistics, which presents the distribution, amplitudes of change, and standard errors of observations/sample cases in each class. We used the default weight of 1 for each observation. The eigenvalue is 7.36 > 1, the Wilks’ Lambda is 0.12 < 0.5 and p-value is less than 0.001, hence the classes are statistically distinct.

Table 3. Classification matrix.

Observed	Predicted
	0	1	2	3
0	20	2	0	0
1	1	8	0	0
2	0	0	3	0
3	0	0	0	2

lists the classification results. 91.7% of the observations/cases were correctly classified. According to Equation (5), the Press Q statistic is 85.3 > 6.63, the critical value at the significant level of 0.01. Therefore, the model has a very good predictive power and the classification was significantly better than chance.

Excluding the 2014 sample plots that participated in the establishment of the discriminant model, the thirty-five D1 values of the seven sample plots from 2009 to 2013 were substituted into the Equations (6)–(9) for hazard prediction, and the results were compared with the field observations.

Table 4. Confusion matrix of the severity of insect infestation.

Observed	Predicted
	0	1	2	3
0	12	3	1
1	2	11
2		2	4
3	0	0	0	2

shows the classification or confusion matrix, which gives an overall accuracy of 77.14%. It means that the model can effectively predict the severity of the insect infestation.

3.4. Insect Outbreak Date Prediction

Traditionally pest outbreak is monitored through on-site surveys, while the outbreak time is predicted using meteorological data, especially temperature data, combined with the pest life history [65]. In this study, the maximum amplitude of the detail component D1 from March to April each year indicates the pest outbreak, thus its corresponding date is the predicted outbreak date. By obtaining the pest outbreak dates from the detail component D1 from March to April each year from 2009–2014 at each sampling site, a confusion matrix was formed to evaluate the accuracy of the method for predicting pest outbreak dates, which is shown in

Table 5. Confusion matrix of pest outbreak date from 2009 to 2014.

Observed	Predicted
	No Insect Infestation	Early March	Late March	Early April	Late April
no insect infestation	35	0	0	0	3
early March	0	16	0	1	0
late March	0	0	0	0	0
early April	0	0	0	13	0
late April	0	0	0	0	3

. This confusion matrix gives an overall accuracy of 94.37%, which means 94.37% of the dates of insect outbreaks are correctly predicted. Errors mainly occurred in late April. Therefore, the detected dates of insect outbreaks are considered to match the observed dates well.

4. Discussion

4.1. Poplar Looper Pest Monitoring

Remote sensing allows non-contact and spatially continuous monitoring of forest pest infestation [66]. In this paper, we presented an approach to the detection and severity classification of the infestation of poplar looper on desert poplars using MODIS NDVI time series through the wavelet transform and discriminant analysis. Gärtner et al. [67] did a similar research. They combined RapidEye high-resolution and Landsat 8 medium resolution data to produce synthetic images using the enhanced spatial and temporal adaptive reflectance fusion model (ESRARFM) [67] to detect poplar forest disturbance caused by poplar looper in a riparian Tugai forest at the Arghan forest station in **njiang. Their purpose was to assess whether the synthetic images could achieve a better disturbance detection accuracy than that achieved by using RapidEye data alone. They used the relative NDVI difference to measure forest defoliation and recovery, and classify the severity of defoliation based on the ratio of defoliation and recovery [68]. Eight RapidEye images (5 m resolution) and ten Landsat 8 images were acquired in the 2003 growing season from the end of March to mid-September with large time windows. One limitation with the study is that the acquisition dates for the high resolution images could miss the foliation peak starting of defoliation, reducing the information about the disturbance. The second limitation is that the relative NDVI differences might not capture defoliation caused by insects. The third limitation is the high cost of high resolution image acquisition, which would prevent the use of the method for real-time monitoring of forest insect infestation and for studying insect dynamics over a long term. MODIS data can overcome these limitations and allow for analysing long time series to capture the natural fluctuations that are inherent to insect dynamics [21,22,23,25,26]. Huang [32], Wang [38] and Qiu [69] used MOD13Q1 and satellite HJ A/B data to study poplar looper infestation. The main limitation of their research is that they only studied whether an infestation had occurred using very short time series data and did not assess the severity of pest infestation and detect the time of its occurrence. Jia [36] and Huang [37] used the same method as that used in Huang [32] to obtain poplar looper infestation information in different years in our study area, simulated the spatial spread of insect infestation, and identified the factors affecting the poplar looper outbreaks. Liu [70] used the MODLT1D temperature product to study the response of the poplar looper in the desert poplar forests to surface temperature in this area. These studies including ours analysed and modelled the infestation of poplar looper in desert poplar forests using remote sensing data from different perspectives, and provided new tools for future large-scale remote sensing monitoring of poplar looper infestation.

4.2. Identification of Outbreaks of Pests

Remote sensing images can capture the changes in vegetation canopy caused by pest infestation [71]. In this study, we applied the wavelet transform to filter the noise and graft the original NDVI values from March to April. The D1 detail component of an NDVI time series obtained after decomposition and reconstruction can highlight the sudden changes in NDVI values. As shown in Figure 9, when the NDVI time series curve suddenly increases or decreases, the D1 curve also fluctuates strongly. Because the poplar looper infestation starts in late March and ends in middle or late April in the study area, the maximum value on the D1 curve from March to April each year can be used to determine whether the poplar looper infestation occurred. For example, substituting the maximum values of D1 in March and April of 2009, 2010 and 2012 respectively (reference lines A, B, C in Figure 9) into the Equations (6)–(9) to determine the severity level of pest infestation at Sampling Site 4 (see Figure 1), the results show that poplar looper infestation occurred at this site in early March 2009, 2010 and late April 2012, which correspond to our field observations.

The outbreak time determined by the method and the outbreak time observed through field surveys match with an accuracy of 94.37%, which shows that the methodology is able to reliably determine the outbreak time.

4.3. Wavelet Analysis

Wavelet analysis has been applied in the field of radio [72] (such as noise reduction for acoustic or video signals), meteorological [73,74,75] and hydrological information analysis [76,77]., image enhancement [78] and image fusion [79]. Hyperspectral images have a very high spectral resolution with hundreds of bands forming an approximately continuous spectral curve [80]. Some scholars have tried to use continuous wavelet analysis for monitoring wheat yellow rust and powdery mildew for winter wheat based on hyperspectral images [81,82,83]. They expanded the spectral data of each pixel in a single hyperspectral image to form a spectral curve, which is similar to an acoustic audio signal, and then applied wavelet analysis to detect or extract specific information about the outbreaks of the fungal diseases in the crops.

In this study, we modelled the numerical curve of each image element of the long time series of NDVI images as a special signal curve based on the characteristics of the periodic and regular changes of the vegetation canopy NDVI, which contains information about the seasonality and trends of the poplar tree growth, as well as localized abrupt changes caused by disturbance events such as droughts, fires, pests, diseases, and leave and timber harvesting or caused by clouds and snows; decomposed NDVI time series data over multiple levels via the wavelet transform; and constructed a smoothed NDVI time series to represent the poplar tree growth under normal conditions. We then added the March and April NDVI data to the smoothed NDVI time series to extract the high-frequency components. The level 1 detail components of the blended time series provide the best indicators of the pest infestation. The result shows that we have successfully detected the pest information. This is a new attempt to apply the wavelet analysis method.

4.4. Limitations

Although the research focused on the poplar looper infestation on desert poplars, the proposed approach can be used to detect and assess forest insect disturbances in general. However, our approach has its limitations.

As we blended the smoothed NDVI time series with the subset of the original time series containing the pest infestation signals, other noises such as those created by clouds and dust storms in the same period remained in the blended time series, which may reduce the detection and severity classification accuracy. It is necessary to identify and remove those noises not caused by pest induced forest disturbances by comparing the original time series and the detail components of the blended time series.

In addition, MODIS NDVI data have a moderate temporal resolution of 16 days and a coarse spatial resolution of 250 m. The presented approach works well in dense poplar forests. But both temporal and spatial resolutions may not be sufficient for accurately determine when and where pest infestations occurred in the areas where poplars are sparsely distributed. Further research will be conducted to fuse MODIS NDVI time series data with other remote sensing data with higher temporal and spatial resolutions to improve the accuracy.

Poplar looper has one generation a year in the Tarim Basin. Its outbreak occurs in March and April. Our approach can be applied directly to an area with the similar natural environment and pest infestation dynamics. Tests are required to verify whether it is applicable to other regions where pests may live for more than one generation, such as the areas south to the Yangtze River in China where pests may have three generations a year and the areas north to the Yangtze River where pests may have two generations a year.

Detection and monitoring of poplar looper infestations on desert poplars by means of remote sensing is still in the stage of exploration and development. Further studies will be conducted to improve the proposed approach with the aim to more accurately detect the poplar looper infestation signals from remote sensing data and enable early detections of poplar forest pest-induced disturbance, which can provide timely information and early warning for effective poplar forest pest prevention and control. Without the disturbance of pests, the growth of desert poplars will be getting better, thereby improving the ecological conditions of its habitat in the Taklimakan Desert.

5. Conclusions

In this paper, we presented an approach to the detection and severity classification of the infestation of poplar looper on desert poplars using MODIS NDVI time series through the wavelet transform and discriminant analysis. Our research proves the potential of wavelet analysis technology based on MODIS NDVI time series in monitoring pests and diseases in dense desert poplar forests. Our method achieved 91.7% accuracy in detecting the poplar looper outbreaks in the middle and lower reaches of the Yerqiang River. Therefore, in the future, we can extend the application of the approach to monitoring the defoliation of trees caused by insect pests or diseases.