Investigating the Correlation between Multisource Remote Sensing Data for Predicting Potential Spread of Ips typographus L. Spots in Healthy Trees

Abstract

In the last decade, thousands of hectares of forests have been lost in the Czech Republic, primarily related to European spruce bark beetle (Ips typographus L.), while more than 50% of the remaining Czech forests are in great danger, thus posing severe threats to the resilience, stability, and functionality of those forests. The role of remote sensing in monitoring dynamic structural changes caused by pests is essential to understand and sustainably manage these forests. This study hypothesized a possible correlation between tree health status and multisource time series remote sensing data using different processed layers to predict the potential spread of attack by European spruce bark beetle in healthy trees. For this purpose, we used WorldView-2, Pléiades 1B, and SPOT-6 images for the period of April to September from 2018 to 2020; unmanned aerial vehicle (UAV) imagery data were also collected for use as a reference data source. Our results revealed that spectral resolution is crucial for the early detection of infestation. We observed a significant difference in the reflectance of different health statuses, which can lead to the early detection of infestation as much as two years in advance. More specifically, several bands from two different satellites in 2018 perfectly predicted the health status classes from 2020. This method could be used to evaluate health status classes in the early stage of infestation over large forested areas, which would provide a better understanding of the current situation and information for decision making and planning for the future.

1. Introduction

During the past decade, rising mean annual temperatures have caused drought prolongation and shifts in growing seasons, which have led to an increased frequency of forest insect attacks and mass propagation of the European spruce bark beetle (Ips typographus L.) in the temperate forests of Central Europe [1,2]. For example, spruce (Picea spp.) and pine (Pinus spp.) stands currently account for 49.2% and 15.9% of forests, respectively, in the Czech Republic [3], and usually grow in homogenous forest stands affected by deer damage and root or stem diseases. Such forest stands are more vulnerable to abiotic factors, particularly strong winds [4]. European spruce bark beetles are capable of eruptive population growth, which usually leads to widespread tree mortality, particularly for Norway spruce (Picea abies (L.) Karst.) stands [5]. Several studies have also observed a connection between European spruce bark beetle outbreaks and other forms of forest disturbances, such as snow damage and windthrows [6,7,8,9]. However, the recent unprecedented bark beetle outbreak in Czechia indicates a shift to more drought-driven dynamics [10]. Unfortunately, the situation continues to escalate, thus posing serious threats for the health and future of forests in central Europe [11]. Despite the abundance of literature and research about the European spruce bark beetle, the driving forces behind such eruptive population dynamics are not fully understood. From a management perspective, attacked spruce forests are difficult to control. Forest managers require thorough knowledge of the: (i) impacted forest area, (ii) bark beetle spatial spread patterns, and (iii) factors influencing the risk of abiotic disturbances on forest stand structures (e.g., windthrow) [12].

In forest monitoring and mapping, extensive scale coverage of damaged forests by field observations is practically impossible [13]. However, remote sensing approaches can quantify the expansion and frequency of European spruce bark beetle outbreaks by monitoring and mapping the spatiotemporal changes over time [14,15]. In addition, drought stress affects the biochemical and biophysical properties of trees, and creates spectral signature variations that can be inferred from moisture content decreases and chlorophyll fluorescence changes [16,17]. Different satellite sensors have been used to identify disturbed forest areas because they offer data collection ability over broad geographic regions at multiple resolutions [18,19]. In addition, very-high-resolution satellites, such as Pléiades 1A/1B (PHRB) and WorldView-2 (WV2), can be utilized to detect disturbed areas (e.g., harvested areas), and estimate and extrapolate individual tree parameters with significant accuracy [20,21,22].

Most studies of the detection of disturbed trees by pests are based on optical remote sensing sensors. One of the main advantages of optical satellite sensors is the increased number of spectral bands; the spectral properties of each band can be effectively used for forest pest disturbance monitoring and mapping. The red edge and green regions of the electromagnetic spectrum (400 to 700 nm) have proven to be very sensitive to changes in chlorophyll content [22]. Furthermore, WorldView-2 has a red edge band with superior sensitivity for early detection of damaged trees by bark beetles [23,24]. Immitzer and Atzberger [25] used satellite images acquired by WorldView-2 for early-stage detection of European spruce bark beetle attacks, with an overall accuracy of 70% for distinguishing green-attacked trees from healthy trees. Mullen [26] also used WorldView-2 for early-stage tree damage detection by montane pine beetles, with an overall accuracy of 75% to distinguish between green-attacked trees and nonattacked trees. In another study, Flichev [27] used WorldView-2 images to assess the damage caused by European spruce bark beetles in the Bistrishko Branishte UNESCO Man and Biosphere (MAB) reserve using several vegetation indices (VIs). His results showed that the plant stress index CSc (λ = 605/760 nm) performed better in detecting stressed trees than the anthocyanin reflectance index (ARI). Finally, Huo et al. [28] used Sentinel-1 and Sentinel-2 images (with Pléiades images as reference data) to detect forest areas attacked by bark beetles. To obtain the precise locations of the attack, they conducted visual interpretations from bitemporal data. Based on images from 2018, the attacked-tree crowns were visible as gray in R-G-B (red-blue-green) and as green in IR-R-B (infrared-red-blue).

Within the Czech Republic, there has been an effort to conduct a systematic forest health assessment. Since 2016, Sentinel-2 cloudless images have been used to produce leaf area index (LAI) maps based on annual intervals. Forest health status is then assessed based on LAI changes during the selected time interval [29,30]. To our knowledge, there have been no previous studies that used UAV images as reference data together with a variety of very-high-resolution satellite platforms and time series data for the detection of European spruce bark beetle dispersal patterns. Additionally, this is the first study that used tasseled cap transformation for the Pléiades 1B satellite platform. The main objectives of this study were to test:

(1)

The correlation between tree health status (standing and lying trees) and time series remote sensing data using different processed layers;

(2)

The best correlated layers between healthy and nonhealthy trees to investigate all those significant differences in, for example, the amount of brightness, greenness, and wetness, among other indicators that could lead to early detection of European spruce bark beetle in healthy trees.

This study aimed to create a blueprint for predicting potential dispersal patterns of European spruce bark beetle from the unhealthy (attacked) to healthy trees, hoping that it could lead to more efficient forest management strategies that provide increased protection for the healthy trees.

2. Materials and Methods

2.1. Study Area

The study was conducted on the Czech side of the Smrčina Mountain; the nature conservation area is located in the eastern portion of Šumava National Park (Figure 1a). The park is part of the Natura 2000 network, which aims to protect the most endangered species and habitats in the European Union. The peak of Smrčina Mountain (1333 a.s.l.) and the saddle to Hraničník Mountain are populated by natural spruce forests that are almost 250 years old, including Calamagrostio villosae-Piceetum and Athyrio alpestris-Piceetum [31], which, at the time of study, were part of a nonintervention regime. The broader area of Smrčina was subject to salvage logging after a serious bark beetle outbreak that ended approximately 10 years ago. During the Herwart windstorm (28–29 October 2017), the forest stands in the Smrčina territory were heavily damaged, and the result was several large windthrow areas. The downed spruce trees were quickly occupied by bark beetles (dominantly Ips typographus L.), allowing the European spruce bark beetle to increase its population density, which resulted in an outbreak. Due to the disturbance event, large, patched mosaics of standing dead trees, lying dead trees, and snags remain in the area (Figure 1b,c).

2.2. UAV Data Acquisition—Preprocessing and Sampling

We used the DJI Phantom 4 Pro quadcopter (Dá-Jiāng Innovations Science and Technology Co. Ltd., Shenzhen, China), which was equipped with an integrated R-G-B sensor of 20 Mpix resolution, a sensor size of 1 inch, aperture f/2.8–f/11, mechanical shutter, and maximum ISO 12,800. For this study, Phantom 4 Pro offered two main advantages; first, the mechanical shutter does not cause geometric distortions, and, second it provides image acquisition with less noise and greater details. Both were significant advantages because the UAV data, apart from being interpreted visually to provide reference data (Figure 2a,b), also georeferenced the available satellite images (see details in Section 2.3).

UAV image acquisition was conducted during the period of 2018 to 2020, and included 12 flights in total, with 4 flights per year at a height of 100 m above the ground with 80% frontal overlap and 70% side overlap. The UAV data were acquired from May to mid-June for 2018, and in June of 2019 and 2020. The following properties was used in preprocessing of UAV data: alignment: high, 40,000 key points, 4000 tie points, medium-quality dense cloud. The first alignment was refined using the optimize cameras method. The spatial resolution for all UAV images was 3 cm. R-G-B aerial images from all years were photogrammetrically processed using structure-from-motion (SfM) in Metashape (Agisoft LLC: St. Petersburg, Russia). To optimize the final 3D models and before the alignment process, we set the accuracy to high. The georeferencing of the point cloud was done exclusively by the onboard global positioning system (GPS) sensor, which geotagged each raster; the GPS signal had high accuracy given the mountaintop location of the study area.

The UAV orthomosaic images were used as reference data to create samples for the years 2018, 2019, and 2020. The samples were selected based on visual interpretation, and they were divided into four different feature classes including standing healthy, standing dead, lying attractive, and lying dead (Figure 3). The primary reason for using visual interpretation was that remotely sensed data are less costly and time consuming than traditional methods. All samples were non-fixed-size polygons (Figure 3). Although samples were randomly selected, we distributed them to cover the entire study area, and the number of selected samples in each class was not fixed (see details in

Table 1. Statistical summary table of the sample classes.

Year	Classes	Sample Frequency	Average Area (m2)	Sum of Area (m2)
2018	1 LA	30	218.62	6558.69
	1 LD	55	164.68	9057.24
	1 SD	18	185.67	3342.04
	1 SH	109	284.23	30,981.56
2019	1 LA	21	228.74	4803.56
	1 LD	66	171.21	11,299.85
	1 SD	19	168.38	3199.16
	1 SH	106	289.03	30,636.98
2020	1 LA	10	349.22	3492.23
	1 LD	89	185.08	16,471.75
	1 SD	45	350.71	15,781.75
	1 SH	68	208.73	14,193.81

¹ LA: lying attractive, LD: lying dead, SD: standing dead, SH: standing healthy.

). After the sample selection, we used the zonal statistics as table tool of ArcGIS Pro V2.7.2 (ESRI Inc., Redlands, CA, USA) to extract the spectral information in each class to evaluate the relationship between the spectral data and sample classes for each year.

2.3. Acquisition and Preprocessing of Satellite Imagery

We used seven satellite images, two by Pléiades 1B (April 2019 and September 2020), three by SPOT-6 (September 2018, July 2019, and July 2020), and two by WorldView-2 (August 2018 and September 2019); see

Table 2. Specifications of all satellite platforms and their sensors used for this study.

Satellite Platform	Imagery Products Resolution	Spectral Bands (nm)	Image Location Accuracy
Pléiades 1B	2 m multispectral—0.50 m PAN	Red: 600–720	With 1 GPSs is 1 m; without GPSs is 3 m (CE90)
		Green: 490–610
		Blue: 430–550
		NIR: 750–950
		PAN: 480–830
WorldView-2	1.64 m multispectral—0.41 m PAN	Coastal blue: 400–450	≤3 m (using a GPS receiver, a gyroscope, and a star tracker) without any GCPs
		Blue: 450–510
		Green: 510–580
		Yellow: 585–625
		Red: 630–690
		Red edge: 705–745
		NIR1: 770–895
		NIR2: 860–1040
		PAN: 450–800
SPOT-6	6 m multispectral—1.5 m PAN	Red: 620–690	35 m CE 90 without GCP within 30 viewing angle cone; 10 m CE90 for Ortho where Reference3D is available
		Blue: 450–520
		Green: 530–600
		NIR: 760–890
		PAN: 450–750

¹ GCPs: ground control points.

for a detailed description of the characteristics of each satellite platform. Unfortunately, no data were found for the study area in 2018 for Pléiades 1B and for 2020 from WorldView-2. To achieve accurate georeferencing for all the satellite images, we used the mosaicked georeferenced image provided by the UAV with an accuracy of 1 cm. In addition, an additional pan-sharpened frame of TripleSat from the year 2018 was used, with an error of less than 25 cm. For georeferencing of Pléiades 1B and WorldView-2, the total error was less than 1 m, while for SPOT-6 it was less than 3 m. Radiometric evaluation for all satellite images was performed visually using a histogram of values for each band. The entire process of satellite image preprocessing was conducted in ArcGIS Pro V2.7.2 (ESRI Inc. Redlands, CA, USA).

2.4. Processing of Spectral Satellite Data

To develop a deeper understanding of the spatiotemporal patterns related to the European spruce bark beetle, several VIs (

Table 3. Reflectance indices for each satellite platform.

Satellite Platform	Used Indices	Formula	Citation
Pléiades 1B SPOT-6	Soil Adjusted Vegetation Index (SAVI)	[1.5 × (NIR − Red)]/(NIR + Red + 0.5)	Huete [38]
	Normalized Difference Vegetation Index (NDVI)	(NIR − Red/NIR + Red)	Rouse et al. [39]
	Simple Ratio (SR)	(NIR/Red)	Birth and McVey [40]
	Transformed Vegetation Index (TVI)	1 SQRT((NIR − Red)/(NIR + Red) + 0.5)	Deering et al. [41]
WorldView-2	Soil Adjusted Vegetation Index (SAVI)	[1.5 × (NIR − Red)]/(NIR + Red + 0.5)	Huete [38]
	Normalized Difference Vegetation Index (NDVI)	(NIR − Red/NIR + Red)	Rouse et al. [39]
	Simple Ratio (SR)	(NIR/Red)	Birth and McVey [40]
	Transformed Vegetation Index (TVI)	SQRT((NIR − Red)/(NIR + Red) + 0.5)	Deering et al. [41]
	Modified Chlorophyll Absorption Ratio Index (MCARI)	[(Red-edge − Red) − 0.2 × (Red-edge − Green)] × (Red-edge/Red)	Daughtry et al. [42]
	Plant Senescence Reflectance Index (PSRI)	(Red − Green)/(Red-edge)	Ren et al. [43]
	Normalized Difference Red-edge Index (NDRE)	(NIR − Red-edge)/(NIR + Red-edge)	Clarke et al. [44]

¹ SQRT: square root.

, see Appendix A), including various wavelength bands (Table 2), were applied to evaluate the status of the vegetation cover between the years 2018 and 2020. WorldView-2 contains a red-edge band; thus, we additionally applied the following vegetation indices (MCARI, NDRE, and PSRI; see details in Table 3, see also Appendix B). Moreover, we used principal component analysis (PCA; see Appendix C) [32] and tasseled cap transformation (Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8; see Appendix D) [33,34,35,36,37] to convert image bands into secondary channels and enhance the spectral differences between trees’ health statuses. For more efficient classification, tasseled cap was used to reduce the multiband dataset to three channels (brightness, greenness, and wetness). Tasseled cap transformation was conducted using Catalyst (PCI Geomatics, Markham, ON, Canada) for all the images.

2.5. Accuracy Assessment

SPSS Statistics V.24 (IBM Corporation, Armonk, NY, U.S.A), Statistica V.12 (Dell Inc., Round Rock, TX, USA), and MATLAB^® R2017b professional edition (MathWorks, MA, USA) were used for the statistical analyses and visualization of the collected data. Based on Kolmogorov–Smirnov test results, our data were not normally distributed, thus Levene’s test was used to assess the equality of variances for all the tested statistical variables for each processed layer. The following descriptive statistics were calculated for each layer within each sample polygon area: minimum and maximum values, range, mean, standard deviation, sum, median, and variance.

To determine the strength of the association between the categorical (independent) and numerical (dependent) variables, the Eta (η) coefficient was used (Equation (1)). Eta-squared (η²) was calculated to measure the coefficient of determination (Equation (2)):

$$\mathrm{Eta} \left(\eta \right)=\sqrt{\frac{S{S}_{effect}}{S{S}_{total}}}$$

(1)

$$\mathrm{Eta} \left({\eta }^2\right)=\frac{S{S}_{effect}}{S{S}_{total}}$$

(2)

where $S{S}_{effect}$ is the sum-of-squares between the groups and $S{S}_{total}$ is the total sum-of-squares.

Dunnett’s T3 test was used to compare the means between the variables (Equation (3)). Dunnett’s T3 test does not assume any equality of variance for analyzing the differences between the groups. A significance level of 0.05 was used to determine whether the above test was statistically significant.

$$D_{Dunnett}=t_{Dunnett} \sqrt{\frac{2M{S}_{S/A}}{n}}$$

(3)

where $t_{Dunnett}$ refers to the critical value in the Dunnett critical value table, $2M{S}_{S/A}$ is mean square within groups, and n is the sample size.

The feature selection tool was used in Statistica to calculate Cramer’s V, which was used to understand the strength of the relationships between the tested variables. Cramer’s V was calculated as:

$$V=\sqrt{\frac{{\phi }^2}{\min \left(k-1, r-1\right)}}=\sqrt{\frac{x^2/n}{\min \left(k-1, r-1\right)}}$$

(4)

where φ is the phi coefficient, $x^2$ is derived from Pearson’s chi-squared test, n is the total observations, k is the number of columns, and r is the number of rows.

3. Results

Spectral Analysis

(1). Eta-Squared (η²)

The results of Eta-squared (η²) showed that for the Pléiades images, although there was a medium-size effect for the blue and red bands in 2020, in general, there was a low-size effect for the Pléiades main bands (more details in Table S1). Our results showed a moderate effect for the tasseled cap transformations derived from Pléiades in 2019 and a low effect for Pléiades in 2020. In addition, a significant size effect (~30%) was detected for VIs, including NDVI, SAVI, SR, and TVI, derived from Pléiades in 2019. Based on Eta-squared (η²) results, SR in 2019 was the most important layer to predict the health status classes for 2018, 2019, and 2020, with size effects of 38%, 38%, and 34%, respectively. In general, Pléiades in 2019 performed better in predicting the health classes compared with the data from 2020 (

Table 4. Summary of the most important results from Eta-squared (η²), PHRB.

Satellite	Year of Imagery	Band	ETA-Squared (2018)	ETA-Squared (2019)	ETA-Squared (2020)
PHRB	2019	NDVI (mean)	0.34	0.33	0.28
		SAVI (mean)	0.34	0.33	0.28
		SR (mean)	0.38	0.38	0.34
		TVI (mean)	0.33	0.32	0.27
		Wetness (mean)	0.21	0.18	0.15
	2020	Greenness (mean)	0.12	0.14	0.04

).

Regarding WorldView-2, coastal blue, blue, green, and red bands had a very high impact (>50%) in detection of health classes. The bands derived from WorldView-2 in 2019 had slightly better performance than in the 2018 image. Greenness derived from the tasseled cap had a low impact on the health status, but contrary brightness (~26–38%) and wetness (32–51%) showed high-size effects for all years (2018–2020). VIs, such as NDVI, SR, TVI, and PSRI, had a moderate (WorldView-2 for 2018) and high (WorldView-2 for 2019) impact on predicting the health status classes for all years (2018–2020; more details in

Table 5. Summary of the most important results from Eta-squared (η²), WV2.

Satellite	Year of Imagery	Band	ETA-Squared (2018)	ETA-Squared (2019)	ETA-Squared (2020)
WV2	2018	Coastal Blue (mean)	0.37	0.36	0.36
		Blue (mean)	0.36	0.34	0.35
		Green (mean)	0.33	0.31	0.32
		Yellow (mean)	0.42	0.40	0.40
		Red (mean)	0.41	0.39	0.40
		Brightness (mean)	0.36	0.38	0.28
		Wetness (min)	0.30	0.30	0.36
	2019	Coastal Blue (mean)	0.52	0.49	0.58
		Blue (mean)	0.51	0.48	0.57
		Green (mean)	0.45	0.43	0.49
		Yellow (mean)	0.46	0.44	0.54
		Red (mean)	0.48	0.46	0.60
		PSRI	0.25	0.24	0.29
		Wetness (min)	0.34	0.33	0.51

and Table S1).

A low association impact was detected for PCA in all years (2018–2020) for Pléiades and WorldView-2 images, and for SPOT-6 in 2020. For SPOT-6 (2018 and 2019), the PCA results showed a moderate impact. In addition, R-G-B bands derived from SPOT-6 (2018 and 2019) had a significant impact on health status prediction, with an impact ranging between 34 and 56%. A moderate impact was detected in the case of NIR bands. In addition, the results showed that tasseled cap channels (from SPOT-6 in 2018 and 2019), such as brightness and wetness, had significant importance in detecting health classes (more details in Table S1).

For SPOT-6 (2018 and 2019), Vis, such as NDVI, SAVI, TVI, and SR, mainly had a moderate to high importance in detecting health status, while for SPOT-6 2020, there was no significant importance for these indices. However, there was a moderate impact for R-G-B bands and the brightness channel derived from SPOT-6 in 2020. In general, SPOT-6 for 2018 and 2019 had much better performance in distinguishing the health classes than 2020 (see

Table 6. Summary of the most important results from Eta-squared (η²), SPOT6.

Satellite	Year of Imagery	Band	ETA-Squared (2018)	ETA-Squared (2019)	ETA-Squared (2020)
SPOT6	2018	Blue (mean)	0.40	0.38	0.39
		Green (mean)	0.38	0.37	0.39
		Red (mean)	0.41	0.40	0.41
		Brightness (mean)	0.36	0.36	0.28
		Wetness (mean)	0.40	0.38	0.43
		SR (mean)	0.21	0.21	0.26
	2019	Blue (mean)	0.40	0.38	0.49
		Green (mean)	0.41	0.40	0.43
		Red (mean)	0.38	0.36	0.46
		Brightness (mean)	0.58	0.58	0.42
		Greenness(mean)	0.27	0.22	0.04
		Wetness (median)	0.27	0.21	0.12
	2020	Blue (min)	0.19	0.18	0.26
		Green (min)	0.15	0.14	0.19

and Table S1).

The R-G-B bands and most satellites were highly significant in the prediction of health classes, most likely because the associated wavelengths were highly correlated with photosynthetic activities. In addition, coastal blue derived from the WorldView-2 images was crucial for health status assessment because it is absorbed by chlorophyll-a in healthy plants, thus it is very effective when conducting vegetative analyses.

Between all the calculated descriptive statistics (mean, max, min, STD, etc.), mean values were determined to be more promising in the prediction of classes (further details in Table S1). Therefore, for further analyses, we only focused on the mean values. Figure 9 shows the differences in mean values for the calculated transformations and VIs from all satellites.

(2). Feature Selection

A feature selection algorithm was used to identify the most important layers in the modeling process for each dependent variable separately (Table S2). The feature selection approach was in agreement with the Eta-squared (η²) results. Some of the main influential layers for the classes in 2018 included: brightness SPOT-6 2019 (Cramer’s V = 0.5453), coastal blue WorldView-2 2019 (Cramer’s V = 0.4962), blue WorldView-2 2019 (Cramer’s V = 0.4855), red WorldView-2 2018 (Cramer’s V = 0.4665), green WorldView-2 2019 (Cramer’s V = 0.4531), and red WordView-2 2019 (Cramer’s V = 0.4522). For the classes in 2019: brightness SPOT-6 2019 (Cramer’s V = 0.5166), coastal blue WorldView-2 2019 (Cramer’s V = 0.4814), blue WorldView-2 2019 (Cramer’s V = 0.4709), brightness WorldView-2 2018 (Cramer’s V = 0.4522), and red WorldView-2 2019 (Cramer’s V = 0.4515). For the classes in 2020: red PHRB 2020 (Cramer’s V = 0.5655), red SPOT-6 2020 (Cramer’s V = 0.5524), blue PHRB 2020 (Cramer’s V = 0.5379), green PHRB 2020 (Cramer’s V = 0.5368), and coastal blue WorldView-2 2019 (Cramer’s V = 0.5063; further details in Table S2).

In general, WorldView-2 had better performance than the other two satellites for predicting the health status classes in 2018 and 2019. For 2020, there was not a significant difference between the satellites. The main bands (R-G-B and coastal blue) and tasseled cap channels (brightness and wetness) had the highest significance and were the most important layers in identifying the health status classes (see Table S2 for further details).

(3). Post Hoc Analysis

In most cases, the results of the multiple comparisons for all years showed that there was a significant difference in the mean values between the health status classes. In Figure 10, the results of Dunnett’s T3 are presented for some of the most important layers according to the feature selection algorithm (Table S3), while Figure 11 shows the bands that were helpful to differentiate all the classes from each other. Our results demonstrated that it was generally easier to distinguish the healthy standing class from the other classes; however, in some cases, it was difficult to distinguish the lying attractive (LA) from standing healthy (SH). For example, the yellow band from WorldView-2 in 2018 could distinguish all the classes except the standing healthy (SH) and lying attractive (LA). Similar to this were the red band from WorldView-2 in 2018, and the coastal blue, blue, green, yellow, red, brightness from WorldView-2 in 2019. In addition, in Figure 10, it is evident there were some challenges to identifying lying attractive (LA) from lying dead (LD) trees. We attributed this to the nature of sample collection that caused a mixture of values from the soil and leaves. Although prediction of all classes was not fulfilled in many of the bands, several bands could separate the classes successfully (Figure 11).

4. Discussion

Currently, there is a severe lack of knowledge about the role of the European spruce bark beetle and its population dynamics via remote sensing, which is proving increasingly useful for predicting forest health issues. This study evaluated the use of multisource time series (2018–2020) remote sensing data; using different processed layers to predict the potential spread of European spruce bark beetle in healthy trees, we investigated spectral correlations between healthy and nonhealthy trees. The location of the study area was ideal for such a study due to the continuous and severe disturbances caused by abiotic and biotic factors during the period of interest of the study.

Notably, the date of image acquisition could significantly affect the image information content for the investigation and analyses of spectral correlations between different datasets. For example, in the image captured by Pléiades (acquired 25 April 2019), due to partial presence of snow on the forest floor, which has different reflectance properties, Eta (η) and Eta-squared (η²) had weaker association with the health classes compared to WorldView-2 (acquired 4 September 2019), which, in terms of spatial resolution, is very similar to Pléiades. In addition, unbalanced sample sizes resulted in insufficient information for some classes (e.g., lying attractive (LA) trees), leading to poorer classification performance; this result was in accordance with other similar studies [45,46].

Differences in the wavelengths of the same bands between different sensors (e.g., Pléiades blue band: 430–550 nm, WorldView-2 blue band: 450–510 nm, and SPOT-6: blue band: 450–520 nm) may also help explain the weaker spectral correlations found between the datasets. However, WorldView-2 has higher spectral resolution and provides an additional band (coastal blue: 400–450 nm), which provided better overall performance. More specifically, the coastal blue band was found to be crucial for tree health status assessment because it is absorbed by chlorophyll-a in healthy plants and aids in conducting effective vegetative analysis. As displayed in Figure 11, a higher amount of chlorophyll-a,b in healthy trees resulted in a lower amount of reflectance in coastal blue and blue bands (higher absorption meant lower reflectance). As the density of chlorophyll decreased, the reflectance of these bands increased proportionally. Essentially, changes in the amount of chlorophyll result in changes in the amount of absorption of coastal blue, as well as other bands in the visible spectrum [47].

The results revealed that spectral resolution was crucial for the early detection of infestation. For example, Figure 11 presents three bands from two different satellites in 2018 that perfectly predicted the health status classes from 2020. During the early stages of infestation, the density of leaf pigments changed; thus, it was more likely to detect the health status using the wavelength absorbed by the chlorophyll-a,b and carotenoids. Therefore, coastal blue and blue bands had much higher performance in tree health status prediction. After the infestation, the cell structure of leaves changed; thus, health status was detected using the wavelengths interacting with these properties (cell structure) [48,49]. NIR absorption/reflectance changed based on the cell structure of the vegetation. Therefore, this band was helpful in distinguishing current health classes.

We determined that spatial resolution was not as crucial as spectral resolution. In this study, SPOT-6 performed as well as Pléiades and WorldView-2. One reason for this finding may be the nature of the sampling method. In this study, each sample polygon covered several pixels, which led to mixed spectral values for each polygon. It can be assumed that if the sample size were smaller and representative of just one pixel (e.g., 2 × 2 m), the results may have been significantly different. Hypothetically speaking, in that case, SPOT-6 would not perform as well as the rest of the satellites due to the lower spatial resolution.

Issues differentiating between the lying attractive (LA) and standing healthy (SH) may be attributable to the high amount of chlorophyll-a,b in the lying attractive (LA) trees. Those trees were the result of windthrow, which meant that they were still alive and potential candidates for future bark beetle attacks. Although the VIs were significant for predicting health classes, and most of the VIs could distinguish healthy or green from the dead trees, distinguishing between standing healthy (SH) and lying attractive (LA) trees was challenging. This was likely due to the similar chlorophyll-a,b absorption and reflectance properties of the cellular structure of the leaves, especially in the NIR band. Our results were consistent with the results of several other studies that emphasized the importance of VIs in predicting health status classes [27,45,46,50,51,52,53].

In this study, the difficulties in separating the lying dead (LD) and lying attractive (LA) trees were attributed. Mixed values of soil, branches, and leaves caused this uncertainty in the samples, which resulted in the model being unable to distinguish between these two classes. Based on the visual interpretation of the 2018 UAV images, a large part of the study area was covered by lying attractive (LA) trees related to the recent windthrow, which were then defenseless against the European spruce bark beetle and were potentially a very good food resource and habitat for their propagation. In 2019, the potential food resources decreased, but it seems that they were still sufficient for the European spruce bark beetle regeneration. That is why approximately only 3% of the standing healthy (SH) trees turned to standing dead (SD).

After two years of exponential procreation (each beetle: 25³; meaning that one bark beetle could procreate at least 15,000 beetles per year), the number of beetles in this area increased to some millions, and, because the available food resources were insufficient as food supply, this meant most of the lying attractive (LA) trees had transitioned to lying dead (LD) trees in 2019; thus, by 2020, the beetles infested the healthy trees. In 2020, approximately 28% of the standing healthy (SH) trees turned to standing dead (SD; Table 1). It can be assumed that the number of the attacked trees was higher, but ~28% could not resist the massive outbreak, which may have occurred due to the health condition of these trees. The results showed a significant difference between these trees (~28%) and the remaining standing healthy (SH) trees, even two years before the European bark beetle attack. To summarize, we surmised that these trees had some health status issues that could not survive the bark beetle attack; the findings of WorldView-2 coastal blue and blue bands from 2018, and the SPOT-6 blue band in 2018 (Figure 11), detected photosynthetic differences between the standing healthy (SH) and standing dead (SD) trees in 2020. It should also be noted that both classes were standing healthy (SH) in 2018. Alternatively, another possible scenario is that these trees in the study area were already under attack (the ~28%) and could be seen in the reflectance, but not phenotypically (after two years, they dried out).

Furthermore, this study showed that a single spectral band could distinguish among some or all of the health status classes. Conclusively, a single band (coastal blue, brightness, wetness, etc.) can be used to successfully detect early European spruce bark beetle spreading. In a study by Liu and Liu [54], multidate images were used to detect dead Robinia pseudoacacia L. plantations using the tasseled cap transformation algorithm. Their results, which correspond with our results, showed that tasseled cap can be efficiently used to monitor forest health status. The findings of this paper also were in line with another study conducted by Erener [55], who also highlighted the importance of tasseled cap for monitoring and health status assessment in a case study of the Seyitömer Lignite Enterprise in Kütahya, Turkey.

We expect that soon, most forest managers will be able to operate all types of UAV platforms that will allow for rapid data collection and real-time data streaming over large forested areas. Such qualitative data could be then used in combination with advanced satellite data to support innovative methods for early detection of European spruce bark beetle spreads into healthy trees.

5. Conclusions

Early detection of bark beetle spreading in healthy trees is of great importance for the effectiveness of forest management strategies and policies. Our study showed that annual index changes from very high and high spatial and spectral image resolutions can effectively be used to detect European spruce bark beetle spreading trends to healthy trees in advance. This new information is of practical value for monitoring purposes, because it can predict a priori European spruce bark beetle outbreaks in healthy trees over large forest areas. A single spectral band can distinguish some or all of the health status classes from each other. Conclusively, a single band can be used successfully for the early detection of European spruce bark beetle spreads. We determined that spectral resolution was more important than spatial resolution for detecting health status classes. Therefore, identifying physiological stress in earlier stages or detecting smaller spectral differences could be more feasible, with higher reliability and accuracy, using hyperspectral sensors in wavelengths that multispectral sensors cannot cover. Different stages of an infestation can change the absorption/reflectance in different wavelengths. Changes in chlorophyll activities result in different absorption/reflectance between the range of 400 to 700 nm. In addition, for trees that are dehydrated or afflicted with a disease, spongy layers within the leaf surface become deteriorated; thus, those trees absorb more NIR than healthy trees do.

Unfortunately, the cloudiness of the Šumava region due to its specific microclimate, made it difficult to obtain clear images over short time intervals. Satellite images with higher temporal resolution (e.g., Sentinels) can potentially increase image capturing during cloud-free days. However, the main disadvantage of these satellites is the lower spatial resolution that would not be suitable for our study due to the size of our samples. Limitations in sampling from certain classes with multispectral sensors (e.g., lying attractive (LA) from 2019 and 2020 and standing dead (SD) from 2018 and 2019) accounted for some of the uncertainties associated with our proposed approach. We believe that thermal sensors could improve the classification of different tree health statuses. Finally, we assume that the composition of several bands using parametric or nonparametric algorithms can lead to early detection of tree health status with high accuracy.