EWMACD Algorithm in Early Detection of Defoliation Caused by Dendrolimus tabulaeformis Tsai et Liu

Abstract

The persistent increase in forest pest outbreaks requires timely detection methods to monitor the disaster precisely. However, early detection is challenging due to insufficient temporal observation and subtle tree changes. This article proposed a novel framework that collaborates multi-source remote sensing data and uses a change detection algorithm to archive early detection of infestation caused by Dendrolimus tabulaeformis Tsai et Liu (D. tabulaeformis) attacks. First, all available Sentinel-2 images with less than 20% cloud cover were utilized. During periods with long intervals (>16 days) between Sentinel-2 images, Landsat-8 images with less than 20% cloud cover were downscaled to a spatial resolution of 10 m using a deep learning algorithm to meet the requirement for a high temporal frequency of clear observations. Second, the spectral index differences between healthy and infested trees were examined to address the challenge of detecting subtle changes in pest attacks. The Enhanced Vegetation Index (EVI) was selected for early defoliation detection. On this basis, the EWMACD (Exponentially Weighted Moving Average Change Detection) algorithm, which is sensitive to subtle changes, was enhanced to improve the capability of detecting early D. tabulaeformis attacks. The assessment showed that the overall accuracy of the change detection (F1 score) reached 0.86 during the early stage and 0.88 during the late stage. The temporal accuracy (Precision) was 84.1% during the early stage. The accuracy significantly improved compared to using a single remote sensing data source. This study presents a new framework capable of monitoring early forest defoliation caused by D. tabulaeformis attacks and offering opportunities for predicting future outbreaks and implementing preventive measures.

1. Introduction

Chinese pine (Pinus tabulaeformis) is one of the endemic coniferous species in China. With increasing stand age and environmental changes, forest pests and pathogens, especially Chinese pine Caterpillar (Dendrolimus tabulaeformis Tsai et Liu, D. tabulaeformis) have become increasingly severe and have killed a significant number of pine trees [1]. This has resulted in substantial losses and persistent threats to forestry, impacting services, structure, carbon, and functionality of forest ecosystems [2,3,4]. D. tabulaeformis is a defoliating pest with a rapid rate of spread [5]. Widespread defoliation could potentially occur in just one month [6]. Therefore, early warning for D. tabulaeformis is crucial for timely intervention. Because of the widespread coverage of pine forests and the high spread area of this pest infestation, traditional forestry inventory methods, such as trapping of larvae, are difficult to conduct and the spatiotemporal observation scale of the results is insufficient [7].

Previously, monitoring forest pests typically used hyperspectral data or high-spatial-resolution multispectral data [2,8]. However, the limited time frequency and high cost make it usually not suitable for large-scale analysis. The changes in crowns caused by pest infestation are gradual and subtle. It could be best characterized through time series analysis [9]. Some scholars have applied time series data with coarse spatial resolution but high temporal resolution, such as MODIS [10,11,12]. This method sacrifices spatial resolution, making it difficult to delineate the area of infestation. Due to its free availability and long-term monitoring capabilities, Landsat data is frequently used for forest pest monitoring [13,14,15]. However, the influence of cloud cover can sometimes result in the actual revisit period of Landsat data exceeding 16 days. Missing one image could lead to missing the entire early stage of pest infestation (≈30 days). With the launch and free data access of the Sentinel-2A and -2B sensors, Sentinel-2 data is widely used for pest monitoring [16,17,18]. However, the insufficient data spatiotemporal resolution used in current research limits the spatiotemporal accuracy of pest monitoring. The location provided by studies based on coarse spatial resolution images is not accurate enough. The studies based on finer spatial resolution images have the potential to improve identification accuracy, such as the results from Daniele et al. [18], and the monitoring results are close to monthly temporal resolution. There is still room for improvement in the monitoring frequency, which is crucial for early pest infestation detection.

To address the insufficient spatiotemporal resolution of observations from a single data source, some scholars have advanced their research by integrating multi-source remote sensing data to obtain denser time series data. Shang et al. [19] near-real-time monitored forest disturbance with harmonized Landsat-7, -8 and Sentinel-2 data. Christopher et al. [20] utilized a Harmonized Landsat Sentinel-2 (HLS) dataset with a temporal resolution of 7 days to monitor forest changes in high-latitude regions. However, previous studies mostly resampled the finer-resolution image to match the coarser-resolution image when integrating data with different spatial resolutions. This approach sacrifices the image’s spatial resolution in exchange for higher temporal resolution [21]. Currently, deep learning-based super-resolution reconstruction is one of the most popular image fusion methods, enabling images to be collaborative at high spatial resolution. Zhang et al. [22] utilized deep learning to fuse Sentinel-2 and Landsat images to track small-scale tropical forest disturbances and mapped the disturbance with a spatial resolution of 10 m. It detected 11% to 21% more disturbance areas compared to using Sentinel-2 or Landsat-7/8 time series alone. Therefore, using deep learning algorithms to downscale low-resolution data and collaborate with high-spatial-resolution data can increase observation frequencies, thereby timely monitoring early pest infestations.

Because the subtle and gradual changes are difficult to detect, current research has developed change detection algorithms sensitive to pest infestations to detect early forest pest outbreaks. Huo et al. [23] proposed a novel vegetation index to detect early European spruce bark beetle infestations based on Sentinel-1 and -2 images. The method achieved identification accuracy ranging from 0.80 to 0.82 during the early stages of infestation. Vojtěch et al. (2021) [24] utilized the random forest algorithm to conduct classification between healthy and bark-beetle-infected trees based on Sentinel-2 images, achieving a 78% overall accuracy in the classification. Christopher et al. [20] employed the BEAST algorithm to monitor forest changes in high-latitude regions. The results indicate that this method can serve as an early warning signal for changes caused by pest infestation. However, most of the aforementioned methods still need improvement in terms of sensitivity to early pest infestations. In studies such as [13], over 50% of Landsat pixels are detected during the transition period between green (early stage) and red/gray stages after being infested, which greatly affects the response speed and accuracy. In studies such as [24,25,26], the accuracies of green-attack (early) stages detection are significantly lower than red or grey-attack stages, which could reach 90% [27,28]. These studies point to a large insufficient in early infestation detection. Identifying such subtle changes often requires more sensitive approaches. The EWMACD (Exponentially Weighted Moving Average Change Detection) algorithm is a change detection algorithm sensitive to subtle changes based on time series data. It has been used to identify minor disturbances in forests, such as selective logging [29,30]. Research has also confirmed the potential of using the EWMACD algorithm for detecting pest infestations [31]. Therefore, the EWMACD algorithm has the feasibility to monitor early D. tabulaeformis infestation.

The objective of this study is to achieve early warning of D. tabulaeformis infestations through collaborative multi-source remote sensing time series data and apply EWMA (Exponentially Weighted Moving Average) change detection algorithm. By collaborating Sentinel-2 and Landsat-8 data, higher temporal resolution time series could be obtained, enabling timely and spatially explicit monitoring of D. tabulaeformis infestation. On this basis, sensitive features of infestation are extracted, and the EWMA change detection algorithm, which is sensitive to subtle changes, is applied to achieve early warning and accurate detection.

2. Materials and Methods

2.1. Study Area

The study area, the mid-west part of Lingyuan City, is located in Liaoning Province, China (Figure 1). Lingyuan City has a total forest area of 136,000 hectares, of which 69,600 hectares are pine forests. The forest in the study area is mostly composed of Pinus tabulaeformis, covering a total area of 3335 hectares. The average temperature in March 2018 was 6.39 °C, and the average temperature in October 2018 was 16.58 °C. As a primary tree species for afforestation on barren hills, Pinus tabulaeformis forests are predominantly artificial pure forests [32,33]. D. tabulaeformis infestation broke out in 2018, with one generation occurring each year in the study area. D. tabulaeformis began climbing trees from mid-March to mid-April after overwintering, feeding on pine needles in small quantities, marking the early stage of infestation. In mid-June, they begin pupation, with adults emerging in early July to lay eggs, during which time they do not harm pine trees. Larvae start appearing from late July to mid-August, feeding extensively on pine needles, displaying the most severe damage, indicating the late stage of infestation. By early October, D. tabulaeformis started descending from trees to overwinter.

2.2. Data and Pinus Tabulaeformis Forest Mask

Based on the physiological habits of D. tabulaeformis, the study period is determined to be from 1 March 2018 to 31 October 2018. As described in Section 3.2.2, the change detection algorithm uses information from the target year and the reference year to detect forest pests. Therefore, all Sentinel-2A and -2B images with cloud cover less than 20% in the study area from 2017 to 2018 were collected, totaling 56 images. Additionally, during periods with long intervals (>16 days) between Sentinel-2 images, Landsat-8 images with cloud cover less than 20% in the study area were downscaled, totaling 16 images. The average revisit period during the study period was approximately 8.6 days. Figure 2 shows the temporal distribution of all the images used from 2017 to 2018.

Sentinel-2 Top of Atmosphere (TOA) images were obtained from the European Space Agency (ESA) and then atmospherically corrected using Sentinel Application Platform (SNAP), and converted into surface reflectance (SR) images. Landsat-8 SR images were downloaded from the Google Earth Engine (GEE) platform.

Since the focus of this study is on identifying forest pests, and the forest area in the study area is mostly pure pine forest, the Pinus tabulaeformis forest mask was created based on the Environmental Systems Research Institute (ESRI) land use/land cover (LULC) dataset during 2017 and 2018, along with field-collected Pinus tabulaeformis samples, and used the random forest classification algorithm.

2.3. Training and Validation Data

The training data used for feature selection: a field survey was conducted in October 2018. D. tabulaeformis infestations were identified through visual inspection of pine needle color and signs of feeding. A total of 5 infestation points and 5 healthy points were collected through field measurements, and their locations were recorded. The damage level of the infestation points is severe. Visual interpretation was also conducted based on available Very High Resolution (VHR) images from Google Earth Pro (GEP) and Sentinel-2 images obtained from 2017 to 2018. Combining field surveys and visual interpretation, a total of 200 training data points were selected, with 100 infestation points and 100 healthy points.

The validation data used for accuracy assessment: validation data were selected through visual interpretation based on available VHR images from GEP and Sentinel-2 images from 2017 to 2018. A total of 100 validation points were selected for both the early and late stages of infestation, with 50 infestation points and 50 healthy points for each stage.

3. Methodology

The overall workflow of this study is illustrated in Figure 3. Firstly, Sentinel-2 images were preprocessed to obtain surface reflectance data. Then, based on the ESRCNN data fusion algorithm, Landsat-8 images were fused with Sentinel-2 images to generate a dense time series with a spatial resolution of 10 m (see Section 3.1). Next, for each pixel, the Enhanced Vegetation Index (EVI), Normalized Burn Ratio (NBR), Red-Green Index (RGI), and SWIR2 were calculated. Then the Jeffries–Matusita (JM) distance for the indices in both healthy and infested points were computed (see Section 3.2.1). EVI was selected as the optimal input feature and EVI time series data were input into the Exponentially Weighted Moving Average Change Detection (EWMACD) algorithm. Infestation is detected when the EWMA value consistently exceeds the control limit (see Section 3.2.2). Finally, post-processing is applied to the detection results, followed by an accuracy assessment (see Section 3.3). The following sections will delve into the key processes in detail.

3.1. Fusion of Sentinel-2 and Landsat-8

Collaborating multi-source data could enhance the temporal resolution of observations, enabling timely observation. To harmonize the resolution disparities between images from different sources, data fusion algorithms could downscale coarse spatial resolution images to match fine spatial resolution images. Therefore, the ESRCNN algorithm proposed by Shao et al. [34] based on deep learning for image super-resolution reconstruction can automatically select the optimal features from one or multiple images for high spatial resolution reflectance prediction. Compared to previous algorithms such as STARFM and ATPRK [35,36], ESRCNN improves image quality and preserves the reflectance distribution of the original images better.

ESRCNN uses fine-resolution Sentinel-2 image as auxiliary data to downscale the coarse-resolution Landsat-8 image. Firstly, ESRCNN employs an adaptive fusion algorithm, instead of the simple resampling method, to downscale Sentinel-2 image bands 11–12 from 20 m to 10 m resolution. Secondly, ESRCNN utilizes Landsat-8 panchromatic bands and Sentinel-2 bands 2–4, 8, and 11–12 to downscale Landsat-8 bands 1–7. This fusion network is designed to accommodate flexible numbers of Sentinel-2 images as auxiliary datasets. In this article, for each downscaling of the Landsat-8 image, three auxiliary Sentinel-2 images are required.

Table 1. The dates of downscaled Landsat-8 images and auxiliary Sentinel-2 images.

Landsat	Sentinel-2
Date	Date1	Date2	Date3
24/1/2017	2/1/2017	12/1/2017	4/2/2017
18/2/2017	12/1/2017	4/2/2017	26/3/2017
13/3/2017	4/2/2017	26/3/2017	15/4/2017
29/3/2017	26/3/2017	15/4/2017	22/4/2017
9/5/2017	15/4/2017	22/4/2017	25/5/2017
1/6/2017	22/4/2017	25/5/2017	14/6/2017
26/6/2017	25/5/2017	14/6/2017	11/7/2017
29/8/2017	5/8/2017	17/9/2017	22/9/2017
9/5/2017	5/8/2017	17/9/2017	22/9/2017
24/11/2017	8/11/2017	3/11/2017	28/11/2017
20/1/2018	2/1/2018	12/1/2018	26/2/2018
2/5/2018	2/1/2018	12/1/2018	26/2/2018
21/2/2018	12/1/2018	26/2/2018	23/3/2018
16/3/2018	26/2/2018	23/3/2018	15/4/2018
1/4/2018	23/3/2018	15/4/2018	2/5/2018
12/5/2018	2/5/2018	7/5/2018	27/5/2018

displays the dates of downscaled Landsat-8 images and auxiliary Sentinel-2 images.

The data fusion results were evaluated by comparing the spatial distribution and spectral value between the fused images and the reference images (Sentinel-2 image on the same or near date). Root Mean Square Error (RMSE) and the Pearson correlation coefficient (Correlation) were used to assess the accuracy of reflectance values. These two metrics are widely employed for evaluating the performance of data fusion algorithms.

The RMSE is used to evaluate the spectral difference between the reference image and the fused image. RMSE is calculated as:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}(\mathrm{R},\mathrm{F})=\frac{1}{\mathrm{M}\mathrm{N}}\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{M}}\sum _{\mathrm{j}=1}^{\mathrm{N}}{\left[\mathrm{R}\left(\mathrm{i},\mathrm{j}\right)-\mathrm{F}(\mathrm{i},\mathrm{j})\right]}^2}$$

(1)

where R represents the reference image, F represents the fused image, M is the number of pixels on the row, and N is the number of pixels on the column.

Correlation is used to indicate the correlation between the reference image and the fused image, calculated as:

$$\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\left(\mathrm{R},\mathrm{F}\right)=\frac{\sum _{\mathrm{i}=1}^{\mathrm{M}}\sum _{\mathrm{j}=1}^{\mathrm{N}}\left[\mathrm{R}\left(\mathrm{i},\mathrm{j}\right)-\mathsf{\mu }\left(\mathrm{R}\right)\right]\left[\mathrm{F}\left(\mathrm{i},\mathrm{j}\right)-\mathsf{\mu }\left(\mathrm{F}\right)\right]}{\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{M}}\sum _{\mathrm{j}=1}^{\mathrm{N}}{\left[\mathrm{R}\left(\mathrm{i},\mathrm{j}\right)-\mathsf{\mu }\left(\mathrm{R}\right)\right]}^2\sum _{\mathrm{i}=1}^{\mathrm{M}}\sum _{\mathrm{j}=1}^{\mathrm{N}}{\left[\mathrm{F}\left(\mathrm{i},\mathrm{j}\right)-\mathsf{\mu }\left(\mathrm{F}\right)\right]}^2}}$$

(2)

where $\mathsf{\mu }$ refers to the mean value.

To validate the applicability of the ESRCNN algorithm for this study, it was necessary to have Sentinel-2 images available on dates close to those of the fusion data. Therefore, we selected Sentinel-2 images collected on 11 July, 16 July, and 31 July 2017 as auxiliary data to downscale the Landsat-8 image collected on 12 July 2017. The result is shown in Figure 4. Visually, the downscaled images exhibit finer spatial details across all bands compared to the original Landsat-8 images. Moreover, the bands closely resemble those of the neighboring Sentinel-2 images.

Figure 5 displays the accuracy assessment results of the downscaled image. Comparing the downscaled images with the original Landsat-8 images and the adjacent date Sentinel-2 images, the RMSE is generally less than 0.1, and the Correlation is generally greater than 0.85. These metrics indicate that the downscaled images are consistent with the original images in terms of spectral and spatial characteristics.

3.2. Defoliation Detection Algorithm

3.2.1. Feature Selection

Selecting features sensitive to D. tabulaeformis infestation is crucial for inputting into the change detection algorithm. D. tabulaeformis infestation results in pine defoliation, needle desiccation, and reduced chlorophyll content. Therefore, we referred to existing studies and selected several vegetation indices and bands (

Table 2. The candidate features to detect D. tabulaeformis infestation.

Index	Formulation	Indicated Change	Reference
EVI	$\frac{2.5\left(\mathrm{N}\mathrm{I}\mathrm{R} - \mathrm{R}\mathrm{e}\mathrm{d}\right)}{\left(\mathrm{N}\mathrm{I}\mathrm{R}+6\mathrm{R}\mathrm{e}\mathrm{d} - 7.5\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}\right) - 1}$	Needle structure	[12,37]
RGI	$\frac{\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$	Coloration	[8,27]
NBR	$\frac{\mathrm{N}\mathrm{I}\mathrm{R} - \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}{\mathrm{N}\mathrm{I}\mathrm{R} + \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}$	Moisture stress& needle structure	[38,39]
SWIR2		Moisture stress	[13,23]

) related to water stress, color, and needle structure. Using the training dataset, we calculated the Jeffreys–Matusita (JM) distance for each candidate feature in both the early and late stages of the damage.

The JM distance quantifies the separability of selected features between infested and healthy pine forests. It ranges from [0, 2], with higher values indicating greater separability. The JM distance is calculated as:

$$\mathrm{J}\mathrm{M}=2(1-{\mathrm{e}}^{-\mathrm{B}})$$

(3)

$$\mathrm{B}=\frac{1}{8}{({\mathrm{m}}_{\mathrm{i}}-{\mathrm{m}}_{\mathrm{j}})}^2\frac{2}{{{\mathsf{\sigma }}_{\mathrm{i}}}^2+{{\mathsf{\sigma }}_{\mathrm{j}}}^2}+\frac{1}{2}\ln \left(\frac{{{\mathsf{\sigma }}_{\mathrm{i}}}^2+{{\mathsf{\sigma }}_{\mathrm{j}}}^2}{2{\mathsf{\sigma }}_{\mathrm{i}}{\mathsf{\sigma }}_{\mathrm{j}}}\right)$$

(4)

where B represents the Bhattacharyya Distance, ${\mathrm{m}}_{\mathrm{i}}$, ${\mathrm{m}}_{\mathrm{j}}$ represent the mean values, while ${\mathsf{\sigma }}_{\mathrm{i}}$, ${\mathsf{\sigma }}_{\mathrm{j}}$ respectively represent the standard deviation of a specific feature within the two categories.

3.2.2. EWMACD Algorithm

EWMACD is an online time series change detection algorithm sensitive to subtle change. The flowchart of EWMACD is illustrated in Figure 6.

Specifically, the EWMACD algorithm first establishes the following harmonic model during training to eliminate the impact of phenological changes on pest monitoring:

$${\mathrm{y}}_{\mathrm{i}}={\mathrm{a}}_0+\sum _{\mathrm{j}=1}^{\mathrm{k}}\left({\mathrm{b}}_{\mathrm{j}}\cos \left(\frac{2\mathsf{\pi }\mathrm{j}}{\mathrm{D}}{\mathrm{t}}_{\mathrm{i}}\right)+{\mathrm{c}}_{\mathrm{j}}\sin \left(\frac{2\mathsf{\pi }\mathrm{j}}{\mathrm{D}}{\mathrm{t}}_{\mathrm{i}}\right)\right)+{\mathsf{\epsilon }}_{\mathrm{i}}$$

(5)

where ${\mathrm{y}}_{\mathrm{i}}$ is EVI value, ${\mathrm{a}}_0$ is the intercept coefficient, $\mathrm{k}$ is the order of harmonic function and set = 2, ${\mathrm{b}}_{\mathrm{j}}$ is the coefficient of cosine harmonic, ${\mathrm{c}}_{\mathrm{j}}$ is the coefficient of sine harmonic, ${\mathrm{t}}_{\mathrm{i}}$ is Julian date, and $\mathrm{i}$= [1, 2, …, N], D is the number of days within a year. ${\mathsf{\epsilon }}_{\mathrm{i}}$ represents the regression residual of Julian date. These coefficients are estimated utilizing the ordinary least squares (OLS) method.

Secondly, all elements with standardized residual values greater than the specified threshold were filtered out through the X-Bar value to remove interference from clouds and outliers and then recalculate the harmonic coefficients:

$$\widehat{\mathsf{\beta }}=({\mathrm{X}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}^{\prime }{\mathrm{X}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}{)}^{-1}{\mathrm{X}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}^{\prime }{\mathrm{t}}_{{\mathrm{n}}_{{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}^{\prime }}}$$

(6)

where ${\mathrm{X}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ is a matrix composed of ${\mathrm{n}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\times \left(1+{\mathrm{c}}_{\mathrm{i}}+{\mathrm{b}}_{\mathrm{j}}\right)$, each harmonic coefficient is adjusted according to the Julian date. After all elements which standardized residuals greater than the threshold are filtered out, $\widehat{\mathsf{\beta }}$ is recalculated. Then, the residuals ($\mathrm{r}\mathrm{e}\mathrm{s}$) were calculated as:

$$\mathrm{r}\mathrm{e}\mathrm{s}=\mathrm{t}-\mathrm{X}\widehat{\mathsf{\beta }}$$

(7)

where $\mathrm{X}$ is the matrix derived from computing the full $\mathrm{n}\times (1+{\mathrm{c}}_{\mathrm{j}}+{\mathrm{b}}_{\mathrm{j}})$. The training period variance ${\mathrm{s}}^2$ was calculated as:

$${\mathrm{s}}^2=\frac{{\mathrm{r}\mathrm{e}\mathrm{s}}_{{\mathrm{n}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}}{\mathrm{r}\mathrm{e}\mathrm{s}}_{{\mathrm{n}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}}^{\prime }}{{\mathrm{n}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}-1}$$

(8)

where ${\mathrm{r}\mathrm{e}\mathrm{s}}_{{\mathrm{n}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}}$ represent the residual on the training period.

EWMA values were calculated as:

$${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}}=\left(1-\mathsf{\lambda }\right){\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}-1}+{\mathsf{\lambda }\mathsf{\epsilon }}_{\mathrm{i}}$$

(9)

where the ${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_1$ = 0. The smoothing parameter $\mathsf{\lambda }$ ranges from 0 to 1. $\mathsf{\lambda }$ close to 0 indicates that a greater weight of the historical period data, ${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}}$, mainly relies on ${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}-1}$. The $\mathsf{\lambda }$ is set to 0.15 (Figure 7).

The control limit of EWMA chart were computed as:

$${\mathrm{C}\mathrm{L}}_{\mathrm{i}}=0\pm \mathrm{L}\mathrm{s}\sqrt{\frac{\mathsf{\lambda }}{2-\mathsf{\lambda }}\left[1-{(1-\mathsf{\lambda })}^{2\mathrm{i}}\right]}$$

(10)

where $\mathrm{L}$ represents the constant control limit, $\mathrm{s}$ represents the standard deviation of all residuals during the training period.

Then, the ${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}\mathrm{C}\mathrm{D}}_{\mathrm{i}}$ was calculated as:

$${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}\mathrm{C}\mathrm{D}}_{\mathrm{i}}=\mathrm{S}\mathrm{i}\mathrm{g}\mathrm{n}\left({\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}}\right)·\mathrm{i}\mathrm{n}\mathrm{t}\left[\left|\frac{{\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}}}{{\mathrm{C}\mathrm{L}}_{\mathrm{i}}}\right|\right]$$

(11)

where $\mathrm{S}\mathrm{i}\mathrm{g}\mathrm{n}\left(\right)$ is used to obtain the plus-minus sign of ${\mathrm{E}\mathrm{W}\mathrm{M}\mathrm{A}}_{\mathrm{i}}$, $\mid \mid$ represents the absolute value, and $\mathrm{i}\mathrm{n}\mathrm{t}\left[\right]$ represents the integer value. The pest infestation was signaled when EWMACD value less than 0.

It is significant to define smoothing coefficients $\mathsf{\lambda }$ when using the EWMACD algorithm. $\mathsf{\lambda }$ represents the weight of historical period data and affects the calculation of control limits of the EWMA chart, which is directly related to the sensitivity of change detection. Observing the detection result using differences $\mathsf{\lambda }$ in Figure 7, pest infestation was difficult to detect when setting $\mathsf{\lambda }$ high, such as $\mathsf{\lambda }$ = 0.2 (Figure 7d). When $\mathsf{\lambda }$ was set low (e.g., $\mathsf{\lambda }$ = 0.1), a significant increase in misclassified pixels (Figure 7f). Therefore, we will initially set $\mathsf{\lambda }$ to 0.15 as the smoothing coefficient (Figure 7e).

3.3. Accuracy Assessment

To assess the spatial accuracy of early and late detection results of insect damage, the precision (Precision), recall (Recall), overall accuracy (OA), and F1-score were computed as accuracy evaluation metrics based on the validation dataset.

$$\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}}$$

(12)

$$\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}=\frac{\mathrm{T}\mathrm{P}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

(13)

$$\mathrm{O}\mathrm{A}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}}$$

(14)

$$\mathrm{F}1-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{2\ast \mathrm{T}\mathrm{P}}{2\ast \mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}+\mathrm{F}\mathrm{N}}$$

(15)

where $\mathrm{T}\mathrm{P}$ represents true positive, $\mathrm{F}\mathrm{P}$ represents false positive, $\mathrm{T}\mathrm{N}$ represents true negative, and $\mathrm{F}\mathrm{N}$ represents false negative.

4. Results

4.1. Feature Selection Results

The time series of each candidate feature, including EVI, NBR, RGI, and SWIR2 were calculated, as shown in Figure 8. The JM distance for the early and late stages of infestation was computed, as shown in

Table 3. JM distance of each feature for the early and late stages.

Index	JM Distance (Early Stage)	JM Distance (Late Stage)
EVI	0.83	1.02
RGI	0.83	0.98
NBR	0.80	1.09
SWIR	0.87	0.86

. Throughout the entire study period, EVI was the earliest and most consistent to display spectral differences between healthy and diseased trees. Also, EVI continuously maintained a higher JM distance. As the aim of this study is to achieve early detection of insect damage, considering the comprehensive comparison of time series indices and JM distances, EVI was ultimately chosen as the input feature for detecting pest infestation.

4.2. Feasibility Analysis of Using EWMACD Algorithm for Pest Detection

Figure 9 shows the time series of EVI and EWMA values for a typical infestation pixel to evaluate the feasibility of the EWMACD algorithm for pest detection. In mid-March, the D. tabulaeformis began to return to the pine trees and slightly nibble on the needles, indicating the early-stage onset of pest infestation. The EWMACD algorithm detected the pest occurrence starting from 23 March 2018. In mid-June, the D. tabulaeformis started cocooning and laying eggs after pupating, during which they did not harm the pine trees. The EWMA value remained unchanged from 16 June to 31 July. In late July, the larvae began to voraciously eat the needles, leading to a significant decrease in the EWMA value starting from 10 August. In early October, the D. tabulaeformis began to descend from the trees to hibernate, and the EWMA value gradually recovered from 14 October onwards, stabilizing after 24 October. In summary, the EWMACD algorithm effectively detects pine caterpillar infestation.

4.3. Assessment of Infestation Detection Accuracy

4.3.1. Overall Detection Results of EWMACD

The output of the EWMACD algorithm is a multi-band raster, with each band representing the detection results for the corresponding date. Based on the spatial clustering and temporal persistence of the pest infestation, we applied post-processing to the detection results. We aggregated all raster values for the early detection period (1 March 2018, to 15 April 2018) to obtain the early detection results, and aggregated all raster values for the entire study period (1 March 2018, to 31 October 2018) to obtain the late detection results, while also applying spatial filtering. The early and late detection results of the EWMACD algorithm for the entire study area are shown in Figure 10. The early damaged area is approximately 868 hectares, while the late damaged area is approximately 1183 hectares. Compared to the early stage, the late-stage outbreak shows a trend of spreading from the center to the surrounding areas, as adult D. tabulaeformis migrate and lay eggs in the surrounding areas after maturation, leading to the hatching of the next generation of larvae, which feed on pine needles.

4.3.2. Assessment of Accuracy in Spatial Domain

The accuracy of the EWMACD algorithm was assessed using validation data.

Table 4. Spatial accuracy assessment of EWMACD results in detecting early and late stage defoliation.

	Early Stage		Late Stage
	Multi-Source Image	Sentinel-2 Image	Multi-Source Image	Sentinel-2 Image
Precision	0.98	0.6	0.89	0.67
Recall	0.77	0.49	0.87	0.94
OA	0.86	0.59	0.87	0.74
F1 Score	0.86	0.54	0.88	0.78

presents the spatial accuracy of the early- and late-stage pest infestation detection results. The F1 score for the early-stage and late-stage detection results of EWMACD were 0.86 and 0.88. The precision of the late-stage results slightly decreased compared to the early-stage results, indicating a slight increase in misclassification errors, but the other metrics showed significant improvement, demonstrating overall good detection performance.

One of the main sources of commission errors is shadows, which can result from differences in image acquisition times and cloud cover, as shown in Figure 11. Despite using Shewhart control charts to remove outliers from the images, the impact of minor shadows may not be sufficient to be identified and removed by the Shewhart control chart. Additionally, due to the small smoothing coefficient used, the algorithm is sensitive to minor changes, leading to some shadows being misclassified as false changes. Figure 11 shows an area of misclassification by EWMACD, where three consecutive images starting from 2 May 2018 (Figure 11c–e) were affected by thin clouds and shadows, resulting in healthy pine forests being erroneously labeled as infested pine forests in panel (Figure 11f).

The main source of omission errors is the sparse distribution of pine forests within individual pixels. Because EVI changes are primarily caused by changes in pine needle structure, sparse pine forest pixels have lower overall leaf litter rates and EVI changes are not significant, leading to undetected changes by the EWMACD algorithm. As shown in Figure 12, in this area with low pine density, compared to Figure 12a, significant tree discoloration can be observed in the GEP image from 27 September 2018 (Figure 12b) and the Sentinel-2 image from 4 October 2018 (Figure 12c). However, both the early (Figure 12a) and late (Figure 12c) detection results show higher omission errors.

4.3.3. Assessment of Accuracy in Temporal Domain

We first compared the temporal accuracy of pest detection using synergistic multi-source data and using Sentinel-2 data alone (Table 4). Apart from the low false negative error in the late-stage results detected using only the Sentinel-2 time series, multi-source time series data significantly improved the detection accuracy in both the early and late-stage pest detection compared to using Sentinel-2 data alone. The lower false negative error in the late-stage detection results using only Sentinel-2 data was significantly better than its false positive error, which may be due to the uneven distribution and longer revisit period of the images. Images were used for model fitting during the study period in 2017, approximately 18 days, which resulted in the model not effectively distinguishing between seasonal, trend, and pest-induced changes during the fitting process. It also caused vegetation phenological changes in 2018 to be identified as pest induced.

In the assessment of temporal domain detection accuracy, the entire early-stage detection results from 1 March to 15 April 2018 are considered accurate detection. As shown in

Table 5. Temporal accuracy assessment results.

	Late = 0	Late ≤ 3	Late > 3	Total
Detection	37	2	5	44
Proportion	84.1%	4.5%	11.4%	100%

, the accuracy of EWMACD timely monitoring pest infestation was 84.1%. The EWMACD algorithm shows a timely response to pest infestation, mostly signaling in the early stage of infestation.

Although the EWMACD algorithm has demonstrated good timely response performance, solving the problem of delayed response still poses challenges. The area shown in Figure 13 for the pest infestation started on 23 March 2018, and it was not until 26 June that a large area was detected, with a delay of more than three-time steps between them. The main source of delay detection is insufficient model fitting, resulting in pest infestation not being marked as changes.

5. Discussion

This study synergy multisource remote sensing imagery aimed to enhance the temporal frequency of observations, select sensitive features of D. tabulaeformis infestation, and improve the EWMACD algorithm for early warning and accurate detection of D. tabulaeformis infestation. However, this method also has its limitations:

Firstly, due to the approximately 40-day time span of early infestation, this study downscaled Landsat-8 images and synergized them with Sentinel-2 images to improve the temporal resolution of observations (≤10 days). However, geographical constraints may still limit achieving the required temporal resolution, especially in cloudy or rainy areas, where even synergizing multiple sensor data may not be sufficient to observe early infestations. Sparse imagery may also fail to meet the requirements of the EWMACD algorithm, affecting its feasibility. Future work could focus on improving the selection of data sources. For example, previous studies have shown that SAR data can reflect forest structure and spectral characteristics [40], and infestations can lead to significant changes in radar signals [41]. Since SAR observations are not affected by cloud cover or rain, using SAR data or synergizing SAR data with optical time series data could significantly increase observation frequency and enhance the algorithm’s capability to detect forest pests.

Additionally, the EWMACD algorithm uses single-feature time series data as input, but many studies have shown that combining suitable bands or indices can better represent spectral changes of land cover, and synergizing multiple features as algorithm inputs could potentially improve change detection accuracy [42]. Future work could focus on improving feature inputs. Some scholars have introduced the MEWMA chart based on the EWMA control chart [43]. The MEWMA control chart offers a multi-feature approach and has shown good performance in quality control, offering potential for future research.

Lastly, many studies have demonstrated the correlation between environmental factors and forest pest outbreaks [44,45,46]. However, this aspect was not fully explored in this study. In addition to detecting pests, future research could incorporate environmental factors such as climate, soil, and vegetation, thereby predicting infestation and improving the early monitoring accuracy of pest infestation.

6. Conclusions

Timely monitoring of forest pests and accurately mapping their spatial distribution are crucial for sustainable forest management. In this study, we combined multi-scale remote sensing images to successfully achieve early detection and mapping of D. tabulaeformis disasters. As a basis for further detection, we first downscaled Landsat-8 images to a 10 m spatial resolution using deep learning algorithm ESRCNN instead of simple resampling, thereby improving the spatial accuracy of pest detection. We collaborated Sentinel-2 and downscaled Landsat-8 images to obtain time series at a higher temporal resolution for timely observation. Building upon this, we compared the performance of different indices in identifying early- and late-stage pest infestation, and EVI was selected as the sensitive feature. Finally, we improved the EWMACD algorithm to enhance its sensitivity to early D. tabulaeformis infestation. The assessment showed that the overall spatial accuracy (F1 score) reached 0.86 during the early stage and 0.88 during the late stage. The temporal accuracy (Precision) was 84.1% during the early stage. The spatial accuracy significantly improved compared to using a single remote sensing data source (The F1 score in early-stage detection has improved by about 0.3, and late stage by about 0.1). In summary, these results demonstrate the effectiveness of our proposed method in early warning and monitoring of D. tabulaeformis disasters. With frequent outbreaks and expanding the spatial scope of D. tabulaeformis disasters in recent years, our method can universally monitor forest pests and provide early warnings, thereby controlling the outbreak of D. tabulaeformis disasters.