Assessing Combinations of Landsat, Sentinel-2 and Sentinel-1 Time series for Detecting Bark Beetle Infestations

ABSTRACT

Bark beetle infestations are among the most substantial forest disturbance agents worldwide. Moreover, as a consequence of global climate change, they have increased in frequency and in the size and number of affected areas. Controlling bark beetle outbreaks requires consistent operational monitoring, as is possible using satellite data. However, while many satellite-based approaches have been developed, the full potential of dense, multi-sensor time series has yet to be fully explored. Here, for the first time, we used all available multispectral data from Landsat and Sentinel-2, Sentinel-1 SAR data, and combinations thereof to detect bark beetle infestations in the Bavarian Forest National Park. Based on a multi-year reference dataset of annual infested areas, we assessed the separability between healthy and infested forests for various vegetation indices calculated from the satellite data. We used two approaches to compute infestation probability time series from the different datasets: Bayesian conditional probabilities, based on the best-separating index from each satellite type, and random forest regression, based on all indices from each satellite type. Five different sensor configurations were tested for their detection capabilities: Landsat alone, Sentinel-1 alone, Sentinel-2 alone, Landsat and Sentinel-2 combined, and data from all satellite types combined. The best overall results in terms of spatial accuracy were achieved with Sentinel-2 (max. overall accuracy: 0.93). The detections of Sentinel-2 also were the closest to the onset of infestation estimated for each year. Sentinel-2 detected infested areas in larger contiguous patches with higher reliability compared to smaller patches. The results achieved with Landsat were somewhat inferior to those of Sentinel-2 (max. accuracy: 0.89). While yielding similar results, the combination of Landsat and Sentinel-2 did not provide any advantages over using Landsat or Sentinel-2 alone (max. accuracy: 0.87), while Sentinel-1 was unable to detect infested areas (max. accuracy: 0.62). The combined data of all three satellite types did not achieve satisfactory results either (max. accuracy: 0.67). Spatial accuracies were typically higher for Bayesian conditional probabilities than for random forest-derived probabilities, but the latter resulted in earlier detections. The approach presented herein provides a flexible disturbance detection pipeline well-suited for the monitoring of bark beetle outbreaks. Furthermore, it can also be applied to other disturbance types.

KEYWORDS:

• forest disturbance
• multispectral
• SAR
• time series
• Bayesian probabilities
• random forest regression

Introduction

Natural disturbances are key to the functioning of forest ecosystems. They include abiotic processes, such as wildfire, drought, and windthrow, as well as biotic processes, such as insect outbreaks, pathogens, and herbivory. By increasing ecosystem heterogeneity in terms of structure, composition, and function, disturbances are essential drivers of ecosystem renewal, succession, and biodiversity (Seidl et al. 2017; Thom and Seidl 2016). In recent decades, however, distinct changes in disturbance regimes have occurred in forests around the globe, with disturbances becoming more frequent and more severe (Seidl et al. 2017). In particular, so-called mega-disturbances, which can alter a forest’s characteristics beyond its recovery potential, are emerging (Millar and Stephenson 2015). Climate change has strongly contributed to shifting disturbance regimes and its impact can be expected to increase in the coming decades (Seidl et al. 2017). Among others, this includes disturbances caused by bark beetles.

Bark beetle outbreaks are among the most common disturbances in forest ecosystems (Hlásny et al. 2021; Schelhaas, Nabuurs, and Schuck 2003; Seidl et al. 2017). Under favorable conditions, they can induce major forest die-off events. For instance, outbreaks of the mountain pine beetle (Dendroctonus Ponderosae Hopkins) led to widespread forest mortality in Northwestern America (Raffa et al. 2008). In Europe, the combination of high growing stock of Norway spruce (Picea Abies L. H. Karst.) and a series of hot, dry summers have recently facilitated unprecedented outbreaks of the European spruce bark beetle (Ips typographus L.). This was the case in 2018, when the hot, dry summer severely weakened host trees (Schuldt et al. 2020) and supported massive outbreaks of I. typographus, resulting in large-scale forest die-off and thus substantial economic losses (Hlásny et al. 2021; Senf and Seidl 2021). In Germany alone, 16% of all spruce growing stock was affected between 2018 and 2020 (Bundesministerium für Ernährung und Landwirtschaft 2021). These large-scale events and their distinct consequences demonstrate the need for early intervention and thus for consistent, operational monitoring over large areas.

Remote sensing has become a key asset in the monitoring of forest ecosystems (Lausch et al. 2016, 2017; Lechner, Foody, and Boyd 2020). In addition to enabling area-based assessments of forest condition, it provides a complementary perspective to ground-based analyses. In particular, time series obtained with Earth observation satellites have enabled continuous, frequent, and consistent monitoring of forests at large scales. In recent years, satellite monitoring of forests has benefited from multiple technological advancements, including the availability of satellite data free of cost, an increase in available computational power and data storage capacities, and from new and improved techniques of data analysis, such as deep learning and dense time series analyses (Holzwarth et al. 2020; Hostert et al. 2015). These tools have been employed in multiple, large-scale, high-impact studies conducted over the last 10 years (Hansen et al. 2013; Hermosilla et al. 2019; Reiche et al. 2021; Thonfeld et al. 2022). Among forest disturbances, bark beetle infestations have been studied particularly frequently, especially in Europe and North America where large outbreaks have occurred in recent decades (Senf, Seidl, and Hostert 2017).

During a bark beetle infestation, various processes change the spectral traits of an attacked tree, enabling the infestation detection via remote sensing (Lausch et al. 2016). By feeding on the living tissue under the bark and inducing pathogens, bark beetle infestations disrupt water and nutrient transport and thus typically kill host trees (Hlásny et al. 2021; Krokene 2015). During the first phase of the infestation, the tree’s water and nutrient transport is disturbed but its crown remains visibly green. This so-called green attack phase is the most difficult to detect (Abdullah et al. 2018; Wulder et al. 2009) but detection during this phase is crucial in managed forests, as the trees have to be removed within a few weeks to prevent the spread of a new generation of beetles (Hlásny et al. 2021). In the subsequent red attack phase, tree crowns turn visibly yellow, red, and brown before the trees starts to lose their needles, while during the final, gray attack phase, only the tree trunk and branches remain. In these two latter phases the spectral difference compared to healthy trees is more distinct (Huo, Persson, and Lindberg 2021) and, accordingly, they are more easily detected.

Most studies assessing bark beetle infestations via satellite remote sensing have relied on multispectral data, especially from the Landsat satellite family: in addition to detections across Landsat’s visible and near-infrared bands (Assal, Sibold, and Reich 2014; Bryk, Kołodziej, and Pliszka 2021), its shortwave infrared (SWIR) bands enable early detection, as they are sensitive to leaf water content, the alteration of which is an early effect of a bark beetle infestation (Abdullah et al. 2018; Goodwin et al. 2008; Meigs, Kennedy, and Cohen 2011). Landsat’s thermal instruments have also been used (Abdullah et al. 2019; Hais and Kučera 2008). In recent years, Sentinel-2 complemented Landsat for mapping bark beetle infestations, with the additional advantage of a higher spatial and temporal resolution (Abdullah et al. 2019; Bárta, Lukeš, and Homolová 2021; Fernandez-Carrillo et al. 2020; Huo, Persson, and Lindberg 2021). Its higher spectral resolution – incorporating three red edge bands – has proven to be beneficial, too, but has seldomly been exploited (Abdullah et al. 2019; Dalponte et al. 2022). To our knowledge, bark beetle infestations have yet to be analyzed based on a combination of Landsat and Sentinel-2 time series. However, this would yield time series with a higher observation density and hence higher chances for an early detection, especially in cloudy areas.

Synthetic aperture radar (SAR) data have also been used to assess bark beetle infestations (Senf, Seidl, and Hostert 2017), although a recent review found only four relevant articles (Hollaus and Vreugdenhil 2019). Of these, Tanase et al. (2018) used L-band ALOS data to map both windthrow and bark beetle infestations, but the accuracies for the latter were lower than those for windthrows. Ortiz, Breidenbach, and Kändler (2013) combined TerraSAR-X data and RapidEye imagery and found that the combination of multispectral and SAR data was the most informative. Ranson et al. (2003) also noted the advantages of combining radar data from different sources but they pointed out that the accuracies were consistently lower than those of Landsat 7 imagery. Nonetheless, SAR imagery, with its cloud-independence, may be a valuable tool in the detection of bark beetle infestations, especially given the importance of quick responses during an ongoing infestation in managed forests. While Huo, Huo, Persson, and Lindberg (2021) found little benefit of including Sentinel-1 data in the monitoring of bark beetle infestations, other studies have shown the value of Sentinel-1’s high spatial and temporal resolution (Hollaus and Vreugdenhil 2019), such as in the detection of insect defoliation (Bae et al. 2022). The combination of multispectral and SAR time series would enable time series with a higher temporal density than currently possible with either method alone.

Two additional aspects would strongly improve existing bark beetle infestation assessments, referring to the temporal domain of analysis: the usage of temporally dense time series and the development of near-real-time (NRT) monitoring of insect disturbances. Many studies of bark beetle infestations have relied on a very limited number of satellite observations (e.g. Bryk, Kołodziej, and Pliszka 2021; Meddens et al. 2013; Skakun, Wulder, and Franklin 2003). In addition, many previous multi-date approaches were based on time series with relatively low, e.g. yearly, temporal resolution (Meddens et al. 2013; Meigs, Kennedy, and Cohen 2011) using techniques such as LandTrendr (Kennedy, Yang, and Cohen 2010). Other studies relied solely on cloud-free or mostly cloud-free images, which limits assessments to regions and/or seasons with low cloud cover (Bárta, Lukeš, and Homolová 2021; Huo, Persson, and Lindberg 2021) and increases the difficulty of detection in areas and/or years with higher cloud cover. Senf, Seidl, and Hostert (2017) therefore advocated the use of dense time series.

Oeser et al. (2017) applied dense time series to map (among others) bark beetle infestations in Central Europe but performed a retrospective analysis of disturbance types instead of assessing infestation timing itself. Thonfeld et al. (2022) assessed tree canopy cover loss in Germany between January 2018 and April 2021 at a high temporal resolution (monthly intervals) but did not specifically target bark beetle infestations. In addition to dense time series, the application of NRT approaches where new scenes are processed immediately may be especially advantageous in the case of bark beetle infestation where quick management responses are crucial (Fahse and Heurich 2011; Hlásny et al. 2021). Such approaches have successfully been applied in similar scenarios where fast action is advised, e.g. tropical forest loss: in a series of studies, Reiche et al. (2018); Reiche et al. (2015), used a Bayesian approach to combine Landsat and Sentinel-1 imagery to quickly detect newly logged tropical forests. Combining the datasets improved the speed of detection and overcame problems in the Landsat imagery caused by cloud cover. However, to our knowledge, this approach has yet to be applied to the detection of bark-beetle-infested areas.

Here, we try to overcome these research gaps by continuously monitoring bark beetle infestations in the Bavarian Forest National Park (BFNP), a mountainous protected forest area in Germany (Figure 1),

Figure 1. The Bavarian Forest National Park. a) Yearly bark-beetle-infested areas based on aerial imagery. b) the location of the Bavarian Forest National Park and the adjacent Šumava National Park in Central Europe.

based on multi-sensor satellite time series with full temporal resolution. We evaluate the performance of Sentinel-1A/B, Sentinel-2A/B, and Landsat 7/8 data as well as combinations thereof in the detection of bark beetle infestations in spruce forests of the BFNP between 2016 and 2019. Specifically, we

1. investigate the suitability of a large variety of multispectral and SAR-based indices in the separation of infested from healthy forest areas,

2. compare Bayesian and random forest regression approaches in the computation of infestation probabilities and infestation detection,

3. assess the spatial accuracy of the respective infestation maps derived from different sensor configurations (Landsat, Sentinel-2, Landsat and Sentinel-2 combined, Sentinel-1, all sensors combined),

4. and estimate the detection timeliness as well as determine the feasibility of early detection of infested areas.

Materials and methods

Study area

The BFNP encompasses an area of 242.5 km² and is located in the country’s southeast, along the border with the Czech Republic, where it adjoins the Czech Šumava National Park (680.6 km²; Figure 1b). Together, the parks form the Bohemian Forest Ecosystem, the largest strictly protected forested area in Central Europe (Heurich et al. 2015). The BFNP is characterized by altitudinal gradients between 550 m and 1453 m a.s.l. that strongly shape local climate and consequently vegetation composition. Depending on the altitude, annual mean air temperatures range between 2.0 and 7.2°C; annual precipitation is between 830 and 2280 mm. These diverse conditions give rise to three key forest types in the BFNP: (1) montane spruce forests, the dominant vegetation community in the park’s highest areas (above 1150 m a.s.l) and largely consisting of Norway spruce; (2) (700 to 1150 m a.s.l.), montane mixed forests in intermediate slope areas, where Norway spruce is interspersed with European beech (Fagus sylvatica L.) and silver fir (Abies alba Mill.); (3) wetland forests in the park’s valley bottoms (below 800 m a.s.l.), where colder and wetter conditions limit beech growth, dominated by Norway spruce as well (Heurich et al. 2010).

Given the dominance of Norway spruce at nearly all elevations, the BFNP is highly susceptible to outbreaks of the European spruce bark beetle. As in many areas of Central Europe, this is at least partly due to past forestry practices, in which the preference for spruce resulted in large-scale, evenly-aged spruce forests (van der Knaap et al. 2020). In the BFNP’s core zone, outbreaks are unmanaged such that infested and dead trees remain standing. This core zone is surrounded by an intervention zone where outbreaks are actively managed through salvage logging (Heurich et al. 2010). In the 1980s, the decision was made to not intervene in ongoing outbreaks inside the core zone, such that vast areas have been affected. Between 1988 and 2020, forest dieback covered an area of 8910 ha (Figure 2), of which 6920 ha have been left standing. Distinct peaks of infestation occurred in 1999 and 2019, causing extensive damage (991 and 624 ha, respectively) during those years (Figure 2).

Figure 2. Time series of area infested by bark beetles in the BFNP. Columns represent yearly increase in infested areas, and lines/dots the cumulative infested area. The study period is marked by the grey box. Infested areas that have been left standing, which were used in this study, are colored in blue.

The high infestation rates and vast uncleared areas allow for continuous monitoring of infestations, which makes the BFNP highly suitable for exploring the utility of remote sensing in the monitoring of bark beetle outbreaks. In fact, the park is already an important testing and validation site for remote-sensing-based approaches (Latifi et al. 2021) and is recognized as a “forest remote sensing hotspot” (Holzwarth et al. 2020). As a result, numerous studies have mapped and monitored bark beetle outbreaks in the park using remote sensing data, including hyperspectral aerial (A. Lausch et al. 2013), very high spatial resolution (VHR) satellite (Latifi et al. 2018), and medium-resolution satellite (Abdullah et al. 2019, 2019; Apdullah et al. 2019; Oeser et al. 2017) data as well as SAR imagery (Tanase et al. 2018).

Data

Reference data

Since 1989, color-infrared aerial imagery has been acquired over the BFNP during the summer months (June – August) and used to create orthoimages. Infested patches larger than five trees are delineated manually from these orthoimages each year using a stereoscopic screen. Each mapped polygon also encompasses additional information, including the date of image acquisition and whether the stand was salvage logged after the infestation (Heurich et al. 2010). This dataset provides a unique long-term and spatially explicit record of bark beetle infestations (Figure 2). While it does not contain temporal information for each infested patch apart from the flight year, the summer flight dates have been chosen to ensure that all die-offs that occurred during the previous year are captured. Consequently, we assume that polygons mapped in a certain year represent infestations from the previous year and adjusted the reference data accordingly. Only infested, uncleared patches in the core zone were considered, to ensure proper assessments of bark beetle infestations and to avoid confusing infested with salvage logged areas.

Satellite data

We used all available data for the study area between 1 January 2015 and 30 April 2021 as acquired by the Landsat, Sentinel-1, and Sentinel-2 satellites. The Landsat imagery includes Landsat 7 Enhanced Thematic Mapper (ETM+) imagery, as well as Landsat 8 Operational Land Imager (OLI) data (Wulder et al. 2019). From the Landsat imagery, only Tier 1 data were considered in this study. Sentinel-2 imagery from the twin satellites Sentinel-2A and Sentinel-2B, operating since 2015 and 2017, respectively, were included (Drusch et al. 2012). To ensure the inclusion of all potential observations, all Landsat/Sentinel-2 imagery, regardless of the cloudiness, was considered. In addition, all available Ground Range Detected (GRD) data from Sentinel-1A and Sentinel-1B were used. While the short-wavelength C-band data of Sentinel-1 are not well-suited for forestry applications, as this wavelength does not adequately penetrate the forest canopy (Woodhouse 2006), Sentinel-1 images were nonetheless included due to their spatial and temporal resolution that is unmatched by any other SAR satellite. In total, we used 406 Landsat, 662 Sentinel-2, and 1007 Sentinel-1 images.

Methods

The procedure used this study is outlined in Figure 3 and can be subdivided into data processing, sampling and class separability assessment, probability computation and infestation detection, as well as evaluation. Besides the Data processing, most steps were carried out using the R and Python programming languages (R Core Team 2023; van Rossum, and Drake 2009). Processing of the multispectral data was carried out on an Ubuntu virtual machine on the CODE-DE platform of the German Aerospace Center with 64 GB of RAM and 8 vCPUs. For the analysis, we used a Microsoft Windows personal computer with 32 GB of RAM and an 8-core CPU.

Figure 3. Workflow of the methods used in this study. Gray boxes indicate the datasets, and white boxes the applied methods.

Data processing

We applied extensive processing steps to generate analysis-ready data (ARD) from the multi-source datasets (Frantz 2019; Truckenbrodt et al. 2019). A consistent preprocessing scheme was applied to both multispectral satellite types to ensure proper inter-calibration and enable the joint use of Landsat and Sentinel-2 as a virtual constellation (Wulder et al. 2015). This was done using the Framework for Operational Correction for Environmental Monitoring (FORCE) (Frantz 2019), as the usage of the ready-to-use Harmonized Landsat Sentinel-2 Dataset (Claverie et al. 2018) was not possible in our case due to its relatively low resolution of 30 m, which is lower than Sentinel-2”s native resolution of 10 m. FORCE applies a comprehensive radiometric correction procedure that includes the correction of atmospheric and topographic effects and accounts for bidirectional reflectance distribution function (BRDF) effects (Frantz 2019; Frantz et al. 2016). FORCE computes quality flags for each pixel, including snow, water, insufficiently illuminated areas, and, especially, clouds and cloud shadows, using a modified version (Frantz et al. 2015) of the Fmask algorithm (Zhe and Woodcock 2012; Zhe, Wang, and Woodcock 2015). Considering parallax effects in the Sentinel-2 imagery, the detectability of clouds is increased, which compensates for Sentinel-2”s lack of a thermal band (Frantz et al. 2018). We used the 30-m Shuttle Radar Topography Mission digital elevation model (SRTM DEM) to account for terrain in the radiometric correction process, and Moderate Resolution Imaging Spectroradiometer (MODIS)-based water vapor estimation in the atmospheric correction of the Landsat imagery (Frantz, Stellmes, and Hostert 2019).

After radiometric correction, the spatial characteristics of Landsat and Sentinel-2 were aligned. Thereby, both types of data were projected into one coordinate system (ETRS 89, UTM Zone 32 N) and cut into tiles of 30 km × 30 km size, i.e. 3000 × 3000 10-m Sentinel-2 pixels or 1000 × 1000 30-m Landsat pixels (see supplementary materials). This the default setting of the FORCE software. Since our study area precisely fit into one of these tiles, all other data were omitted from further analysis. During this process, the 20-m Sentinel-2 data was adjusted to 10 m using the ImproPhe approach (Frantz et al. 2016). The previously generated quality flags were used to remove pixels with invalid values, e.g. due to cloud and snow cover, as well as saturated bands.

Next, several vegetation indices were computed from the multispectral data (Table 1) at a pixel size of 10 m, whereby the Landsat data was resampled using nearest neighbor interpolation. All of the indices have been used in vegetation monitoring (as indicated by the references in Table 1). They cover the whole range of Landsat and Sentinel-2 bands across the electromagnetic spectrum. Besides, all indices could be easily computed and exported through the FORCE software (Frantz 2021).

Table 1. Overview of the vegetation indices computed from the Landsat and Sentinel-2 data. ρ represents the reflectance in a given band.

Name	Sensor	Formula	Reference
Normalized Difference Vegetation Index (NDVI)	Landsat/ Sentinel-2	$NDVI=\frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR}+{\rho }_{red}}$	Rouse et al. (1973); Tucker (1979)
Normalized Difference Moisture Index (NDMI)	Landsat/ Sentinel-2	$NDMI=\frac{{\rho }_{NIR}-{\rho }_{SWIR1}}{{\rho }_{NIR}+{\rho }_{SWIR1}}$	Gao (1996)
Normalized Burn Ratio (NBR)	Landsat/ Sentinel-2	$NBR=\frac{{\rho }_{NIR}-{\rho }_{SWIR2}}{{\rho }_{NIR}+{\rho }_{SWIR2}}$	López García and Caselles (1991)
Enhanced Vegetation Index (EVI)	Landsat/ Sentinel-2	$EVI=2.5\ast \frac{{\rho }_{NIR}-{\rho }_{red}}{{\rho }_{NIR}+6\ast {\rho }_{red}-7.5\ast {\rho }_{blue}+1}$	Huete et al. (2002)
Tasseled Cap Brightness (TCB)	Landsat/ Sentinel-2	$\mathrm{T}\mathrm{C}\mathrm{B}=\hspace{0.17em}0.2043\ast {\rho }_{blue}+0.4158\ast {\rho }_{green}+\hspace{0.17em}0.5524\ast {\rho }_{red}+0.5741\ast {\rho }_{NIR}\text{break}+\hspace{0.17em}0.312\ast {\rho }_{SWIR1}+\hspace{0.17em}0.2303\ast {\rho }_{SWIR2}$	Crist (1985)
Tasseled Cap Greenness (TCG)	Landsat/ Sentinel-2	$\mathrm{T}CG=-0.1603\ast {\rho }_{blue}-0.2819\ast {\rho }_{green}-0.4934\ast {\rho }_{red}+0.7940\ast {\rho }_{NIR}-0.0002\ast {\rho }_{SWIR1}\mathrm{\setminus }\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}-0.1446\ast {\rho }_{SWIR2}$	Crist (1985)
Tasseled Cap Wetness (TCW)	Landsat/ Sentinel-2	$\mathrm{T}CW=0.2043\ast {\rho }_{blue}+0.4158\ast {\rho }_{green}+0.5524\ast {\rho }_{red}+0.5741\ast {\rho }_{NIR}+0.3124\ast {\rho }_{SWIR1}\mathrm{\setminus }\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}+0.2303\ast {\rho }_{SWIR2}$	Crist (1985)
Tasseled Cap Disturbance Index (TCD); unscaled	Landsat/ Sentinel-2	$TCD=TCB-\left(TCG+TCW\right)$	Frantz (2021); Healey et al. (2005)
Kernel NDVI (kNDVI)	Landsat/ Sentinel-2	$kNDVI=\frac{1-k}{1+k}$ with $k=\mathrm{e}\mathrm{x}\mathrm{p}\frac{-{\left({\rho }_{NIR}-{\rho }_{red}\right)}^2}{2\ast {\sigma }^2}$ with $\sigma =0.5\ast \left({\rho }_{NIR}+{\rho }_{red}\right)$	Camps-Valls et al. (2021); Frantz (2021)
Normalized Difference Tillage Index (NDTI)	Landsat/ Sentinel-2	$NDT=\frac{{\rho }_{SWIR1}-{\rho }_{SWIR2}}{{\rho }_{SWIR1}+{\rho }_{SWIR2}}$	van Deventer et al. (1997)
Normalized Difference Red Edge Index 1 (NDRE1)	Sentinel-2	$NDRE1=\frac{{\rho }_{rededge2}-{\rho }_{rededge1}}{{\rho }_{rededge2}+{\rho }_{rededge1}}$	Gitelson and Merzlyak (1998)
Normalized Difference Red Edge Index 2 (NDRE2)	Sentinel-2	$NDRE2=\frac{{\rho }_{rededge3}-{\rho }_{rededge1}}{{\rho }_{rededge3}+{\rho }_{rededge1}}$	Ju et al. (2010)
Red Edge Chlorophyll Index (CRE)	Sentinel-2	$CLRE=\frac{{\rho }_{rededge3}}{{\rho }_{rededge1}}-1$	Gitelson et al. (2003)
Normalized Difference Vegetation Index Red Edge 1 (NDVIre1)	Sentinel-2	$NDVIre1=\frac{{\rho }_{BNIR}-{\rho }_{rededge1}}{{\rho }_{BNIR}+{\rho }_{rededge1}}$	Gitelson and Merzlyak (1998)
Normalized Difference Vegetation Index Red Edge 2 (NDVIre2)	Sentinel-2	$NDVIre2=\frac{{\rho }_{BNIR}-{\rho }_{rededge2}}{{\rho }_{BNIR}+{\rho }_{rededge2}}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)
Normalized Difference Vegetation Index Red Edge 3 (NDVIre3)	Sentinel-2	$NDVIre3=\frac{{\rho }_{BNIR}-{\rho }_{rededge3}}{{\rho }_{BNIR}+{\rho }_{rededge3}}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)
Normalized Difference Vegetation Index Red Edge 1 narrow (NDVIre1n)	Sentinel-2	$NDVIre1n=\frac{{\rho }_{NIR}-{\rho }_{rededge1}}{{\rho }_{NIR}+{\rho }_{rededge1}}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)
Normalized Difference Vegetation Index Red Edge 2 narrow (NDVIre2n)	Sentinel-2	$NDVIre2n=\frac{{\rho }_{NIR}-{\rho }_{rededge2}}{{\rho }_{NIR}+{\rho }_{rededge2}}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)
Normalized Difference Vegetation Index Red Edge 3 narrow (NDVIre3n)	Sentinel-2	$NDVIre3n=\frac{{\rho }_{NIR}-{\rho }_{rededge3}}{{\rho }_{NIR}+{\rho }_{rededge3}}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)
Modified Simple Ratio Red Edge (MSRre)	Sentinel-2	$MSRre=\frac{\frac{{\rho }_{BNIR}}{{\rho }_{rededge1}}-1}{\sqrt{\frac{{\rho }_{BNIR}}{{\rho }_{rededge1}}}+1}$	Chen (1996)
Modified Simple Ratio Red Edge narrow (MSRren)	Sentinel-2	$MSRre=\frac{\frac{{\rho }_{NIR}}{{\rho }_{rededge1}}-1}{\sqrt{\frac{{\rho }_{NIR}}{{\rho }_{rededge1}}}+1}$	Fernández-Manso, Fernández-Manso, and Quintano (2016)

All Sentinel-1 data were processed using the Google Earth Engine cloud processing platform (Gorelick et al. 2017). A consistent processing scheme was applied to the data. The Sentinel-1 Ground Range Detected (GRD) data available on Google Earth Engine are preprocessed, including thermal noise removal, radiometric calibration, and terrain correction based on SRTM 30 m DEM (Google Developers 2022). Additional processing steps were applied according to the workflow described in Mullissa et al. (2021). As a first step, we applied a border noise correction to remove faulty backscatter values along the edges of individual scenes and a multitemporal Quegan speckle filter with a temporal window size of 10 images (Quegan and Yu 2001). Since the BFNP is situated in a mountainous area, radiometric terrain flattening (Hoekman and Reiche 2015; Vollrath, Mullissa, and Reiche 2020) was used to reduce the influence of terrain on the backscatter values and to ensure that ascending and descending Sentinel-1 orbits could be used in conjunction. Again, the SRTM 30 m DEM was applied during this process. The indices listed in Table 2 were computed and all images were exported at a spatial resolution of 10 m. Analogous to the multispectral data, we reprojected the exported Sentinel-1 data, subdivided it into tiles, and discarded data not matching the study area. The final result was a consistent data cube comprising all available Landsat, Sentinel-2, and Sentinel-1 data, quality-masked and at a spatial resolution of 10 m.

Table 2. Overview of the indices computed from the SAR imagery. σ represents SAR backscatter in dB and $\gamma$ linearly scaled SAR backscatter.

Index	Sensor	Formula	Reference
VV Backscatter (VV)	Sentinel-1	${\sigma }_{VV}^0$	Holtgrave et al. (2020)
VH Backscatter (VV)	Sentinel-1	${\sigma }_{VH}^0$	Holtgrave et al. (2020)
VV/VH Ratio (RA)	Sentinel-1	$VVVH={\gamma }_{VH}^{\hspace{0.17em}}/{\gamma }_{VV}^{\hspace{0.17em}}$	Holtgrave et al. (2020)
Radar Vegetation Index (RVI)	Sentinel-1	$RVI=\frac{4{\gamma }_{VH}^{\hspace{0.17em}}}{{\gamma }_{VV}^{\hspace{0.17em}}+{\gamma }_{VH}^{\hspace{0.17em}}}$	Charbonneau, Trudel, and Fernandes (2005); Holtgrave et al. (2020); Nasirzadehdizaji et al. (2019)

The computed index time series were then smoothed to remove seasonal/phenological effects as well as noise in each pixel’s time series, which is especially important when data from multiple sensors are combined, and thus obtain a time series with stronger contrasts between healthy and infested forest plots. Here, we applied Locally Estimated Scatterplot Smoothing (LOESS) (Cleveland 1979). LOESS is an implementation of local regression in which a new regression model is fit to each data point individually (Hastie, Tibshirani, and Friedman 2009). The algorithm utilizes a polynomial fit to the data and is hence similar to the Savitzky-Golay filter (Savitzky and Golay 1964). Here, we used the default options implemented in the respective function of the stats R package and thus applied a smoothing degree of 0.75.

Sampling and separability assessment

The identification of infested areas required reference data for both already infested areas and healthy forest areas. However, the available reference dataset was derived from a yearly deadwood inventory and thus covered only infested forest areas. We therefore adopted a sampling strategy that enabled data sampling from both infested and healthy areas. This strategy is based on the assumption that forests patches infested in a recent year and inventoried as such could not have been infested during the previous years, or they would have already been inventoried.

First, we selected 50% of all pixels marked as infested between 2016 and 2019, leaving the other 50% for the validation data set and using stratified random sampling to represent all infestation years equally in the training and validation data. This timeframe marked the period common to all sensors but excluded the year 2015, which was required for the data sampling of infestations occurring in 2016. We chose a 50/50 training/validation separation in order to ensure a sufficient number of pixels from the years 2016 and 2017 in the validation data, as the total infested area and hence the corresponding pixel count was relatively low in these years (Figure 2).

For the training data, we computed yearly median composites of the satellite data for each calculated multispectral/SAR index. Median composites were suitable because seasonal and noise effects had already been removed with LOESS smoothing as described above. Consequently, forest areas that were (still) healthy as well as forests that had (already) fully died-off were expected to show relatively low changes in values during the sampling year. Every composite together with reference data from the year after each composite’s target year was then sampled to collect training data covering healthy forests. For example, the composite for 2018 was sampled together with 2019 infestation data, as areas marked as infested in 2019 were still healthy and stable in 2018. Conversely, the 2018 composite was sampled based on the 2017 reference data, as areas infested in 2017 were assumed to have died-off in 2018. This sampling design is illustrated in Figure 4.

Figure 4. Sampling design as exemplified for one pixel representing an infested area in 2018. The red line represents the changes in NBR values for one example pixel infested in 2018 over time. The left (right) gray column indicates the period covering the sampling of values representing healthy (infested) forest. The infestation occurred during the period represented by the orange box.

To determine the indices that best separated healthy from infested areas, we compared the sampled values of both healthy and infested areas for every remote sensing index. This assessment was performed to select the most appropriate indices for the models of Bayesian approach described below. Specifically, we aimed at selecting the vegetation indices with the highest separability from a) the joint multispectral indices of Landsat and Sentinel-2, b) the indices based solely on Sentinel-2, and c) the Sentinel-1 SAR indices. Separability was estimated by computing the Jeffries-Matusita (JM) distance metric for every index, according to the procedure outlined in Reiche et al. (2015); Reiche et al. (2018). The JM distance has a finite range from 0 to 2, with 0 indicating complete inseparability, i.e. full overlap between the classes, and 2 representing complete separability, i.e. no overlap between the classes (Laliberte, Browning, and Rango 2012; Reiche et al. 2015). This distance metric was chosen because it has been suitable in analyzing class separability in a large variety of remote sensing studies (Dalponte et al. 2013; Laliberte, Browning, and Rango 2012; Padma and Sanjeevi 2014; Reiche et al. 2015, 2018). Besides, one of its properties is that it overemphasizes small distance values. This was especially important in our case, as the separability between healthy and infested was expected to be relatively low for the four SAR-based vegetation indices. With a metric over-emphasizing small differences, the differences between these indices should be relatively large, allowing us to select the appropriate SAR vegetation index with higher confidence (Kavzoglu and Mather 2000).

Computation of infestation probabilities and infestation detection

With our sampling strategy, we gathered values for healthy forests on one side and gray-attack forest areas on the other. In the next step, infestation probability time series were derived from the computed index time series using two different approaches: a Bayesian approach, based on studies by Reiche et al. (2018, 2015) and a random forest (RF) approach (Breiman 2001). This is based on the assumption that – once an infestation began – the spectral indices would start to resemble the gray-attack state and hence infestation probabilities should increase, even during the previous infestation stages. Of the 25-index set, the index of each satellite type that best separated infested from healthy plots according to the JM distance was selected. $v1$ denotes the best-performing index from Landsat and Sentinel-2, $v2$ the best-performing Sentinel-2-only index (Table 1), and $v3$ the best-performing Sentinel-1 index.

Univariate probability density functions based on a Gaussian fit were computed for $v1$, $v2$, and $v3$, respectively, after which the respective probability densities $p\left(v1|H\right),p(v2|H)$, and $p(v3|H)$, i.e. that an observation represented a healthy forest, plus the probability densities $p(v1|I),p(v2|I)$ and $p(v3|I)$, i.e. that an observation represented an infested forest, were computed for all observations (Reiche et al. 2015). Then, assuming equal priors for $H$ (healthy) and $I$ (infested), we derived the conditional probability of an infestation using Bayes’ theorem (McElreath 2020):

(1) $P\left(I|v{1}_t\right)=\frac{p\left(v{1}_t|I\right)}{p\left(v{1}_t|I\right)+p\left(v{1}_t|H\right)}fort\in T_{v1}$(1)

where $t$ represents the date in the time series of $v1$ (Reiche et al. 2015). The conditional probabilities for $v2$ and $v3$ were computed analogously. For the detection of infested plots, five different sensor configurations and their capabilities were considered:

1. Landsat, using $P\left(I|v{1}_t\right)$ with a subset of $t$.

2. Sentinel-2, using $P\left(I|v{2}_t\right)$; if the best performing Sentinel-2-based index was one shared by Landsat and Sentinel-2, $v2=v1$, and $v1$ would be used with a different subset of $t$ that was different from that used by the Landsat configuration.

3. Landsat and Sentinel-2 combined, assessing the potential of this “virtual constellation” of satellites (Wulder et al. 2015) and using $P\left(I|v{1}_t\right)$ for all $t$

4. Sentinel-1 only, using $P\left(I|v{3}_t\right)$.

5. All three sensors combined, using both $P\left(I|v{1}_t\right)$ and $P\left(I|v{3}_t\right)$, arranged by $t$.

In the RF regression approach to the derivation of infestation probabilities (Breiman 2001), the algorithm is based on two previously developed methods: decision trees and bagging. The RF procedure fits a multitude of decision trees, each time incorporating only a subset of the variables and observations in the dataset. The decision trees are thus decorrelated and predictive accuracy is increased (James et al. 2013). RF has been widely employed in remote sensing, due to its high prediction accuracies and easy applicability to high-dimensional data (Belgiu and Drăguţ 2016). We employed this approach to also test the performance of a multivariate procedure, compared to the univariate Bayesian approach. We fit classification RF models and used the information whether a reference area was healthy or infested as the binary response variable. Continuous class probabilities were derived by averaging the binary class votes over all trees (instead of choosing a majority vote as in a classification scenario). The same five sensor configurations as in the Bayesian approach were considered:

1. Landsat, using the same RF as Landsat/Sentinel-2, but with a different subset of $t$ than Landsat/Sentinel-2 (see below).

2. Sentinel-2, using a RF fit to all 21 Sentinel-2 indices.

3. Landsat/Sentinel-2, using a RF fit to all 10 indices shared by Landsat and Sentinel-2, as well as all available observations of these two sensors.

4. Sentinel-1, using a RF fit to all 4 Sentinel-1 indices.

5. Combination, combining the derived probabilities of Landsat/Sentinel-2 and Sentinel-1, arranged by $t$.

In this study, we used the RF implementation of the randomForest R package. All models were fit with 1000 trees, otherwise retaining the standard parameters of the R package (Liaw and Wiener 2022). In all random forest models, we iteratively removed the variable that had the lowest variable importance in the model and assessed the corresponding out of bag (OOB) errors in the models at each step. These errors consistently were the lowest with the largest respective model, i.e. the one incorporating all parameters for a given data configuration. Thus, only the largest model per sensor configuration was used. In summary, ten models were fit, consisting of five different sensor configurations and two different approaches (Bayesian and RF).

Evaluation

Every sensor configuration was assessed individually, and the following steps were applied: all time series of infested plots between 2016 and 2019 according to the reference data that had not been used to train the models were transformed into probability time series. For the probability thresholds, which increased from 0.01 to 1 by 0.01, each pixel’s time series was classified as healthy or infested, according to the respective threshold. The date at which this detection first occurred was then recorded. To reduce the possibility of false detections due to data variability, each detection had to be confirmed by two additional consecutive observations above the respective threshold, similar to various other remote sensing time series detection algorithms (Hirschmugl et al. 2020; Decuyper et al. 2022).

To increase confidence in the results, detections before the probable beginning of the infestation for each year, as well as after September 30^th of the following year (the approximate end of the bark beetle season) were marked as non-detections. To determine the approximate beginning of the infestations for each year, we implemented a modified version of the PhenIps model (Baier, Pennerstorfer, and Schopf 2007) which simulates bark beetle phenology based on topoclimatic data. Specifically, we used the RITY model (Ogris et al. 2019), a modified version of PhenIps that only requires air temperature data. Since PhenIps (originally developed for Austria) better represents climatic conditions in the BFNP, since RITY was developed for Slovenia, the specific thresholds in RITY were adapted to these in PhenIps (Bárta et al. 2022). Daily mean and maximum air temperatures were provided by the nearby weather station of the BFNP. The infestation onsets for each year are listed in Table 3.

Table 3. Infestation Onsets for the Investigated Years, as derived by using PhenIps.

Year	Date of Infestation Onset
2016	7 May 2016 (day of year: 128)
2017	15 May 2017 (day of year: 135)
2018	21 April 2018 (day of year: 111)
2019	18 May 2019 (day of year: 138)

Since the infestation dynamics in the BFNP typically involve one major outbreak wave in spring (see supplementary materials), the estimated infestation onsets should represent a good approximation to the actual infestation dates for many infested patches in the BFNP. The latter cutoff date of September 30^th of the following year was chosen to allow infestations occurring at the end of a given year, which likely only facilitate spectral changes in the following year, to be detected as well. All detections between these dates (beginning of the infestation in a given year and September 30^th in the year after) were considered.

For every threshold/sensor configuration, the overall accuracy plus the commission and omission errors of the infested class were evaluated, according to the recommendations of Olofsson et al. (2014). To determine whether the size of an infested patch influenced its detectability, the recorded measures were evaluated separately according to the number of adjacent 10 × 10 m pixels per patch.

Due to the importance of early infestation detection, we performed a statistical comparison of the detection date by sensor combination and method. While the absolute accuracy of the detection date could not be determined, due to the lack of respective information in the reference data, this timeliness analysis was still able to reveal differences between the sensor configurations and thus an earlier response of one sensor configuration or algorithm over another to the beginning of an infestation. Lastly, we chose the best possible index-configuration-threshold combination in terms of spatial accuracy and produced a final, temporally explicit infestation map for the whole BFNP based on it.

Results

Computed vegetation indices

While none of the computed vegetation indices were able to fully separate healthy from infested forest plots, there were distinct differences among them (Figure 5). Seven indices achieved a JM distance > 1.5, and ten a JM distance < 1. The best possible separability was achieved by the Chlorophyll Red Edge Index (CRE) and the Normalized Difference Red Edge Index 1 (ND1), with a distance of 1.73. The best multispectral index derived from both Landsat and Sentinel-2 data was the Normalized Burn Ratio (NBR), with a JM distance of 1.48, closely followed by the Normalized Difference Moisture Index (NDM). The results achieved with the “original” Normalized Difference Vegetation Index (NDVI) and the recently introduced Kernel NDVI (Camps-Valls et al. 2021) were comparably good as well (JM distances of 1.35 and 1.37, respectively). Of the ten indices with a JM distance < 1, four were SAR-based, whereby the JM-distances calculated from the original satellite bands (VH and VV) were somewhat higher than those resulting from their transformations which, for one red edge index allowed almost no separability at all.

Figure 5. Jeffries-Matusita distance for all computed vegetation indices. For the abbreviations of the indices, see Tables 1 and 2.

The best-performing multispectral and SAR indices were used in further analyses, i.e. the NBR of Landsat and Sentinel-2 ($v1$), the CRE for Sentinel-2 only ($v2$), and the VV band of Sentinel-1 ($v3$). The distribution of sampled healthy vs. infested values for the three indices is shown in Figure 6. The overlap between the value distributions representing healthy and infested pixels was relatively low for CRE and NBR whereas for VV the two populations largely overlapped. Figure 6 also shows the Gaussian probability density functions computed from the samples. The results of the separability analysis are mirrored by the variable importances in the random forest models fit for the respective approach. While the ranks are not always identical, the overall patterns are well in line (see supplementary figure S3).

Figure 6. Value distributions between healthy and infested forests for the three indices selected for processing: a) CRE, b) NBR, c) VH. Dashed lines represent the probability density functions derived from the data for the Bayesian probability computations.

Spatial accuracy

Bayesian approach

Overall, the Bayesian approach (Figure 7) achieved better results than the RF regression in terms of spatial accuracy. Especially at higher probability thresholds, comparably good results were obtained with all three multispectral sensor configurations (Landsat, Sentinel-2, Landsat and Sentinel-2 combined), with overall accuracies > 0.8 or even 0.9. Sentinel-2 consistently outperformed all other sensor configurations in terms of overall accuracy, with a maximum of 0.93 at a probability threshold of 0.98. Commission and omission errors consistently were lowest for Sentinel-2 as well, with omission errors (Figure 7b) close to 0 for most thresholds, and only increasing at very high probability thresholds. Commission errors were high for very low probability thresholds, but quickly decreased with increasing threshold values. Landsat performed slightly better than Landsat/Sentinel-2, even though their accuracy metrics are consistently very similar. Both exhibit similar values compared to Sentinel-2, but consistently score somewhat worse.

Figure 7. Overall accuracy (a), omission error (b), omission error of the infested class (c) commission error for the infested class based on the Bayesian approach.

The results achieved with Sentinel-1 and the Combination configuration were clearly inferior. Commission errors (Figure 7c) were consistently very high for both configurations, even exhibiting values of 1 for thresholds < 0.3. Only for higher values, they decreased. At the same time, omission errors quickly increased for thresholds > 0.3, whereby they were much higher for the Sentinel-1 configuration compared to the Combination configuration, at higher thresholds.

Random forest

The overall accuracies and errors of the three multispectral configurations obtained with the RF approach differed from those obtained with the Bayesian approach (Figure 8). While the results achieved with Landsat and Landsat/Sentinel-2 were similar, those of Sentinel-2 were inferior to both, especially at low probability thresholds, with lower overall accuracies compared to the Bayesian approach and especially high commission errors of 1 for thresholds < 0.25. Results were generally better at higher thresholds, but the RF-based overall accuracies of Sentinel-2 were never as high as either those of Landsat and Landsat/Sentinel-2 or those obtained with the Bayesian approach. Maximum overall accuracies of Landsat/Sentinel-2 did not greatly differ from the accuracies of Landsat. Similar to the Bayesian approach, the high overall accuracies obtained using RF were driven by low omission errors, which only increased for higher thresholds for all three multispectral sensor configurations. This was also the case for Sentinel-1 and the combination of all types of data, but the high commission errors (Figure 8c) of these two configurations resulted in overall accuracies ≤ 0.68 for the Combination configuration and ≤ 0.54 for Sentinel-1.

Figure 8. Overall accuracy (a), omission error (b), omission error of the infested class (c) commission error for the infested class based on the random forest approach.

In general, lower accuracies were achieved with the RF approach than with the Bayesian approach. The best results were those of the Bayesian Sentinel-2 sensor configuration based on the CRE and with a probability threshold of 0.98, which had the highest overall accuracies and lowest omission errors. The overall accuracy for this configuration-approach-threshold combination was 0.931, with an omission error for the infested class of 0.087 (i.e. a producer’s accuracy/sensitivity of 0.913) and a commission error for the infested class of 0.053 (i.e. a user’s accuracy/specificity of 0.947). This equals a precision of 0.962, a recall of 0.909, and thus an F1 score of 0.935. The respective Cohen’s kappa was 0.879. A detailed confusion matrix for this configuration-approach-threshold configuration are given in Table 4, expressed in terms of the proportion of area (Olofsson et al. 2014). The optimal probability thresholds and respective overall accuracies for every sensor configuration are listed in Table 5.

Table 4. Confusion Matrix for the Classification using the Sentinel-2 configuration with the Bayesian approach based on the CRE and a probability threshold of 0.98.

		Reference
		Infested	Healthy	Total
Prediction	Infested	43.93%	1.71%	45.64%
	Healthy	4.44%	49.92%	54.36%
	Total	48.37%	51.63%	100%

Table 5. Optimal thresholds and overall accuracies (OA) for every sensor configuration and approach.

Sensor configuration	Bayesian approach	RF approach
Landsat	Index: NBR; Threshold: .93, OA: .89	Threshold: .88, OA: .86
Sentinel-2	Index: CRE; Threshold: .98, OA: .93	Threshold: .90, OA: .78
Landsat/Sentinel-2	Index: NBR; Threshold: .96, OA: .87	Threshold: .96, OA: .87
Sentinel-1	Index: VH; Threshold: .67, OA: .62	Threshold: .93, OA: .58
Combination	Indices: NBR, VH; Threshold: .77, OA: .67	Threshold: .94, OA: .67

Influence of patch size

As shown in Figure 9 for Sentinel-2 (Bayesian approach), the spatial accuracies differed depending on the size of an infested forest patch, whereby the corresponding number of adjacent pixels within 10 × 10 m represented the patch size. The patch size classes in Figure 9 were chosen to ensure approximately equal pixel counts in every class. The overall patterns were similar between all patch sizes and resembled the general metrics from Figure 7. However, for pixels originating from smaller patches the overall accuracies were consistently lower and the omission/commission errors consistently higher than those from larger patches. For example, pixels from a patch with a size > 100 pixels had maximum overall accuracies > 0.95 (Figure 9a), with comparably low omission and commission errors. By contrast, pixels from patches with ≤ 5 pixels had lower maximum overall accuracies (<0.85) and correspondingly higher commission errors. The omission errors (Figure 9c) also differed but these differences were less pronounced.

Figure 9. Patch-size-specific spatial accuracies for the Sentinel-2 configuration based on the Bayesian approach: a) overall accuracy, b) omission error for the infested class, c) commission error for the infested class. The gray dashed lines represent the accuracy/errors regardless of patch size, as shown in Figure 7.

Timeliness of the detection

The detection dates obtained with the three sensors and their combinations were distinctly different (Figure 10). In each case, the probability threshold with the highest spatial accuracy was applied (Table 3). To ensure comparability between the included years, all detection dates were standardized to the onset of the bark beetle season of the of the respective year (Table 3). With the exception of Sentinel-1, two distinct peaks in the detection date were observed using either approach. The first peak developed 50–200 days after the onset of the bark beetle season and the second peak roughly one year later. Earlier detection dates were consistently achieved using the RF approach. For the majority of the infested areas, the earliest detections, achieved with Sentinel-2 and the RF approach, were within the first 200 days after the infestation. Only a minority of pixels marked as infested fell between the two peaks. There were especially few detections during the winter, although this varied among the configurations, with fewer detections for Sentinel-2 than for Landsat/Sentinel-2.

Figure 10. Distribution of values regarding the timeliness of detection for all sensor configurations based on the Bayesian approach (a) and the random forest approach (b). All dates are standardized to days since the beginning of the bark beetle infestation in the respective year.

The pattern of two distinct peaks encompassing a large share of the infestations was consistent not only between most sensor configurations (Figure 10) but also between the surveyed years. Figure 11 shows the detection periods for the Sentinel-2 RF configuration for the years 2016–2019. While there were some differences between years, e.g. later detections in 2017 and 2018, the general pattern of two peaks was consistent. Most infestations were detected during the first 200 days after the onset and between 300 and 400 days, i.e. roughly one year after.

Figure 11. Distribution of values regarding the timeliness of detection for the years 2016–2019, as determined using the Sentinel-2 random forest configuration.

Likewise, the pattern of two distinct peaks was consistent among infested patches of different sizes (Figure 12; using the same classes as in Figure 9), without major differences between the surveyed patch sizes. In contrast to the spatial accuracies depicted in Figure 9, detections in larger infested patches were not distinctly earlier than those in small patches. It should be noted that the detection dates changed if different probability thresholds were applied. Figure 13 provides an example for the Sentinel-2 RF configuration: at thresholds <0.9, infestations were typically detected earlier, whereby detections per date stabilized at around 50 to 75 days after the infestation onset.

Figure 12. Distribution of values regarding the timeliness of detection in patches of different sizes for the years 2016–2019, as determined using Sentinel-2 random forest configuratio

Figure 13. Distribution of values regarding the timeliness of detection according to three different probability thresholds (0.5, 0.75, 0.9), as determined using the Sentinel-2 random forest configuration.

Final infestation map

The final map of detected infestations in the BFNP was derived from the best index-configuration-threshold combination in terms of spatial accuracy, and thus based on Sentinel-2 based on the CRE and the Bayesian approach, with a probability threshold of 0.98 (Figure 14). As Sentinel-2 was not launched until mid-2015, the map represents only the infestations occurring since then, with the exception of the dark green areas seen along the eastern margin of the park. These correspond to high elevations dominated by mountainous spruce forests, where infestations began between the mid-1990s and mid-2000s, but they were still detectable in the early years following the launch of Sentinel-2.

Figure 14. Final infestation map for the Bavarian Forest National Park. Colors represent the detection date based on the Sentinel-2 configuration using the CRE and the Bayesian approach with a probability threshold of 0.98. The black box indicates the subset illustrated in Figure 15.

The subset of this map shown in Figure 15 demonstrates that the infested areas detected by Sentinel-2 were in good agreement with the reference data. Typically, the areas detected by Sentinel-2 were somewhat larger than those manually mapped from the aerial imagery. Some infested areas were missed, as indicated by the red boxes in the northwestern and southeastern parts of the map (Figure 15).

Discussion

Distinct differences in the spectral separability of healthy vs. infested forest plots could be discerned based on the JM distances of the tested indices. While for seven indices the JM distance was > 1.5, indicating only slight overlap, for 10 indices it was < 1. There were also differences between the five tested sensor configurations. Based on the Bayesian approach, which yielded a higher spatial accuracy, Sentinel-2 consistently outperformed the other four sensor configurations. Landsat and the combination of Landsat and Sentinel-2 achieved satisfactory results but failed to provide any advantages over Sentinel-2 alone. The two configurations that included SAR data performed worse than the three multispectral configurations.

Figure 15. Subset of the infestation map for an area in the northwest region of the BFNP. Polygons represent the disturbance patches in the reference data, labeled according to the infestation year in each one. The three red rectangles delineate areas that were not successfully detected. The satellite-based assessment agreed well with the ground data.

Influence of indices

All indices with the greatest ability to separate infested from healthy plots incorporated at least one Sentinel-2 red edge band. While red edge satellite imagery has rarely been applied in bark beetle infestation detections, its general suitability was previously documented. For example, the studies by Abdullah et al. (2018) and by Adamczyk and Osberger (2015) state that based on spectrometer data, the red edge region of the electromagnetic spectrum is relatively sensitive to bark-beetle-induced changes. Particularly, the red edge region is sensitive to losses in chlorophyll content (Abdullah et al. 2018) which in turn is an early biophysical change after a bark beetle infestation (Wulder et al. 2006), such that red edge indices should be used in assessments of bark beetle disturbance, and especially in early detections. This is also documented by a previous study utilizing multiple Sentinel-2-based indices (Dalponte et al. 2022), where the CRE, the best performing index in this study, achieved good results as well.

The best-performing multispectral indices were the NBR and the NDM, both of which incorporate the NIR and one of the two SWIR bands shared by Landsat and Sentinel-2. Their good performance is well in line with previous research showing that changes in needle water content are closely associated with changes in SWIR reflectance, which both NDMI and NBR use (Foster et al. 2017; Gao 1996). A loss of needle water is among the early responses of a tree to a bark beetle infestation (Abdullah et al. 2019; Wulder et al. 2006). Hence, the NDMI has been routinely employed in satellite-based detections of bark beetle infestations (Fernandez-Carrillo et al. 2020; Goodwin et al. 2008; König et al. 2020; Meddens et al. 2013; Ye et al. 2021). While the NBR was initially developed for the detection of burned areas (Chuvieco et al. 2019; López García and Caselles 1991) and is widely used in that setting, it has also frequently been applied to bark beetle detections (Meigs, Kennedy, and Cohen 2011; Ye et al. 2021).

Comparison of suitable sensors and sensor configurations

Despite their spectral similarities, the better performance of Sentinel-2 than Landsat is plausible, because Sentinel-2 offers a higher spatial resolution and is thus better able to capture small-scale dynamics, resulting in higher spatial accuracies. The maximum overall accuracies achieved by Sentinel-2 were very high compared to other studies with similar objectives. For example, in the study of Huo, Persson, and Lindberg (2021) the overall accuracy was 0.91, although only during late infestation stages (0.82 for earlier stages). In the study of Bárta, Lukeš, and Homolová (2021) the overall accuracy was 0.78, but only green-attack phases, which are more difficult to detect, were examined. Fernandez-Carrillo et al. (2020) reported accuracies of 0.95, but only for high-severity areas. Similar to our study, Dalponte et al. (2022) achieved maximum overall accuracies of 0.91 (even for early attack stages), but with individually delineated tree crowns based on lidar data as an auxiliary data set. The maximum overall accuracies of this study are the highest so far achieved by studies in the BFNP, where Abdullah et al. (2019) achieved a producer’s accuracy of 0.67 for Sentinel-2 and 0.36 for Landsat, whereby green attack was specifically targeted. The usage of L-Band SAR for this task (Tanase et al. 2018) achieved maximum overall accuracies of 0.88. An overall accuracy of 0.7 was achieved based on synthetic RapidEye data in an area in the northern part of the BFNP (Latifi et al. 2018).

Sentinel-2 also had a better temporal performance than Landsat, due to its higher temporal resolution. However, the performance differences between Sentinel-2 and the Landsat/Sentinel-2 combination are more difficult to explain. Although the use of FORCE ensured consistent preprocessing and the harmonization of Landsat and Sentinel-2 data, differences may have persisted, such as due to limitations of the applied radiative transfer models (Doxani et al. 2018). Another explanation is the mixed pixel effects arising from Landsat’s lower spatial resolution, which persisted even after Landsat was resampled to a 10 m pixel size. The most likely explanation though is probably that the NBR was used in the Landsat/Sentinel-2 configuration. This index was less able than the CRE of Sentinel-2 to separate healthy from infested areas, which is conductive to the hypothesis that a high spectral sensitivity to the physical processes induced by a bark beetle infestation is more important than a high temporal resolution in obtaining precise and rapid detections. Namely, based on our results, we found out that the usage of Sentinel-2 alone with a lower observation density in the time series but the greatest physical sensitivity achieved better results than the combination of Landsat and Sentinel-2 with more dense time series. This even holds true in a cloud-prone area like the BFNP and matches other studies reporting similar spatial accuracies in the detection of infested areas but utilized time series with lower temporal densities (Bárta, Lukeš, and Homolová 2021; Huo, Persson, and Lindberg 2021). In addition, this hypothesis is supported by Figure 11: even though data availability for Sentinel-2 was lower in the 2016–17 period compared to 2018–19, the overall detection patterns are still similar. As long as basic time series density per season is guaranteed, detections are not occurring later. Fewer observations are still suitable to capture the infestation-induced changes in canopy reflectance.

Besides, the usage of the NBR vs. the CRE also explains the high similarity between the Landsat and Landsat/Sentinel-2 configurations, both relying on this index. The better performance of Sentinel-2 in larger, continuous forest patches than in patches represented by only a few pixels can be attributed to the lower share of mixed pixels along the edges of each patch, as previously reported (Fernandez-Carrillo et al. 2020; Zimmermann and Hoffmann 2020).

The relatively poor performance of the SAR data from Sentinel-1 can be explained by the low separability between healthy and infested plots based on the calculated indices, as also reported by Huo, Persson, and Lindberg (2021). This can be attributed to the relatively low penetration depth of Sentinel-1’s C-band-based imaging (Reiche et al. 2018; Woodhouse 2006). Thus, the similarity in the infestation probabilities for pixels representing healthy and infested plots resulted in only small changes in the infestation probability after an infestation had occurred. Hence, at low thresholds the detections based on Sentinel-1 data were removed if they occurred before the beginning of the bark beetle season whereas at high thresholds the infestations were not detected because the probabilities remained low. In the combination time series, the consistently medium probabilities of Sentinel-1 added high variability to the time series if combined with multispectral probabilities that increased after an infestation. This reduced the probability of a detection that has to be confirmed by two subsequent observations.

The overall poor performances of the two sensor configurations that included Sentinel-1 data were in agreement with the study by Huo, Persson, and Lindberg (2021) in which the benefit of including Sentinel-1 data in their bark beetle infestation detection approach was minimal. While this configuration had the highest observation density, it appears that, as discussed above, sensor suitability is more important than a high time series density in early detection, as long as a minimum time series density is available to capture bark beetle-induced changes, i.e. to ensure that recent changes in canopy reflectance (or backscatter) are documented without the need for additional observations within each infestation stage. The limited suitability of the combined sensor configuration can also be traced to two different sensor concepts (optical vs. radar) and thus to differences in the sensitivities of the different spectral traits: while optical data is more sensitive to changes in leaf properties, SAR is generally more sensitive to canopy structure (Reiche et al. 2018; Woodhouse 2006). Combining these different approaches will therefore yield inferior results despite overall accuracies that are higher than those obtained with Sentinel-1-alone, as in the latter case the detections are driven by optical data and are not disrupted by less well-suited SAR data. This explains the difference between our study and that by Reiche et al. (2018), in which Sentinel-1 was used to capture logging, i.e. an activity resulting in more distinct structural change than that imposed by bark beetle infestations (Akbari and Solberg 2022; Reiche et al. 2015, 2018). In that study, both Sentinel-1 and the longer-wavelength ALOS-PALSAR data were able to capture the changes, but the use of Sentinel-1 data in the detection of bark beetle-infested areas, where tree trunks typically remain standing for years after tree die-off, is not possible. In summary, combining Landsat and/or Sentinel-2 with Sentinel-1 at the data level to map forest disturbances that are not stand-replacing or otherwise result in only small structural changes (e.g. bark beetle infestations), appears to be not possible using the methods applied here. This finding is backed up by the study by Balling et al. (2021) who assessed whether tropical forest fires – resulting in low structural changes as well – can be more easily captured with the combination of Landsat, Sentinel-2, and Sentinel-1, but also found that a data-level combination did not suit their analysis. This distinction may also explain while in our study, the simple thresholding of Bayesian and RF-derived probabilities provided better results than employing the iterative Bayesian updating of probabilities applied by the authors (Reiche et al. 2015, 2018) available through the bayts R package (Reiche 2019), which we also tested (see supplementary materials). This procedure appears to be better suited to their research problem where both the SAR and multispectral probabilities changed relatively quickly. In addition, the results of the bayts package for the single-sensor time series are likely inferior to our results because the iterative Bayesian updating added additional complicacy to the detection procedure, leading to more detections being filtered out by our evaluation procedure.

The insufficient applicability of SAR data found in this study contrasts with the results of Tanase et al. (2018) who found that infested areas are detectable with L-Band data after needle loss due to the higher penetration depth and lower saturation of their L-band data (Woodhouse 2006). An additional limitation of SAR that led to the low detectability of the infested plots might have arisen from the combined use of ascending and descending orbits. The radiometric terrain flattening according to Vollrath, Mullissa, and Reiche (2020) may not have fully accounted for the different viewing geometries between these two orbits, especially in steep terrain. In addition, the processing pipeline in Google Earth Engine relies on SRTM digital elevation data. Digital surface data (incorporating canopy height) would likely yield better results for side-looking SAR. In addition, other SAR-based metrics, e.g. interferometric coherence, may be tested for the detection of bark beetle-infested areas as well and may result in a higher sensitivity to small structural changes (e.g. a higher coherence as trees without needless move less during wind). In addition, more polarimetric data may be tested as well as other evaluation algorithms, e.g. based on shadow effects (Bouvet et al. 2018).

Possibilities for the early detection of infestations

Regardless of the sensor configuration, both the spatial accuracies and the timeliness of detection should be interpreted jointly, because they will differ depending on the applied probability threshold. Low probability thresholds typically result in much earlier detections, but if a detection occurs too early, i.e. before the onset of the infestation in the reference year, it will not be considered (see Methods). This effect may be due, e.g., to the small changes in the Sentinel-1-derived probability time series, resulting in relatively high commission errors for low thresholds. High thresholds – although they are less likely to be met even in the case of an infestation – reduce commission errors because there will be fewer detections before the infestation onset. This will hold true up until a certain threshold specific to each sensor configuration, after which the overall accuracy will sharply decline because omission errors distinctly increase at high thresholds (Figures 7, 8). A scenario not taken into account in this study is the possibility of earlier stressors that change a tree’s spectral traits even before an infestation occurs, such as drought, storm damage, or snow- and ice-related crown damage (Hlásny et al. 2021) which may result in detections before the estimated infestation onset. Hence, they will be filtered out by the proposed procedure. In the absence of information in the reference data on these stressors, we cannot account for their influence in this study. This would have to be assessed with additional data.

While there is an empirically determined optimal probability threshold for all sensor configurations, there is a tradeoff, with low probability thresholds leading to early detections, many of which will be false-positives and thus decrease the overall accuracy. Low thresholds will include fewer infested pixels than higher thresholds but the infestations will be detected earlier. Higher probability thresholds result in higher overall accuracies but at the cost of later detections. If the procedure proposed herein were to be used for operational monitoring or in a NRT scenario, users could select a threshold that fits their needs: rapid, NRT detection with relatively high errors (which may lead to additional costs if larger areas are surveyed that are in fact not infested) or detections at later stages of an infestation and thus with lower error rates, such as may be preferred for yearly infestation/damage inventories.

Although a detailed temporal assessment could not be performed based on the available reference data and we could only assess the temporal offset from the estimated beginning of the infestation in each year, infestations were mostly detected several months after this date. Thus, while our approach is appropriate for relatively quick, automated inventories of areas infested by bark beetles, such as conducted by government agencies or protected area administrations, it likely lacks green-attack detection capabilities and cannot be applied for practical forest management. The generational development period of the European spruce bark beetle is between 7 and 11 weeks (Fahse and Heurich 2011), with the new generation able to infest previously healthy trees. Thus, to limit bark beetle spread, infested trees have to be removed within this 7–11-week period. This implies the detection of infested trees within a few initial weeks, which cannot be accomplished with our approach, as most detections occur only after this period.

Comparison of the Bayesian and random forest methods

Generally, both the Bayesian and the RF approaches provided reliable results in detecting bark beetle infestations. However, the performance of the Bayesian approach was better in terms of the achieved spatial accuracies, perhaps because of its better adaptability to time series analyses (Reiche et al. 2015; Zhao et al. 2019) whereas, due to the internal bootstrap process, RF approaches are hampered by the arbitrary sub-selection of samples of different time steps (Belgiu and Drăguţ 2016; Breiman 2001). In addition, the probabilities based on the Bayesian approach were very close to zero for most of the included time series but after the infestation they typically increased to values close to one. By contrast, the RF-based probabilities were typically higher from the beginning and the respective increase after the infestation was less pronounced. In addition, some RF-based time series never reached the same maximum probability values. Hence, at higher probability thresholds, the variability in the time series may have led to infestations not being detected as they would not have been confirmed by subsequent observations. All of this may be caused by correlations between the different indices in the training dataset which we observed. However, as described in the methods, for each sensor configuration, we used the random forest model with the most predictors as it achieved the lowest OOB errors.

Limitations and future requirements

Although the methods applied in this study resulted in very accurate detections of infestations and, more importantly, of infestation dates, several methodological adjustments might enhance the spatial and temporal accuracies. First, reference data with more detailed temporal information, including precise infestation dates for single trees or small stands, should be used in further studies. While temporal differences among the five different tested sensor configurations could be identified, we were not able to compare the detection dates computed here with any reference dates, such that several assumptions had to be made about the infestation dates in the reference data, e.g. no infestations before the modele onset of the infestation according to the PhenIps model and that identified patches in the reference data were infested the year before aerial imagery was acquired. Whether late detections stem from late (undocumented) infestation dates or the inability of the applied methods to detect early infestations is unclear. Similarly, the nature of the two distinct peaks obtained with most configurations (Figure 10) could not be resolved. The second peak may have represented infestations that occurred late in the year and thus remained undetected, due to the low observation densities during the winter, but which were then detected early in the next year, when observation densities were higher. Thus, field data with high temporal resolution remains a prerequisite for fully valid temporal assessments of bark beetle infestations via remote sensing.

Future studies should investigate the possibility of using the procedure proposed here without smoothing the time series. Applying the LOESS smoothing was necessary in our case and improved the results (see supplementary materials). We hypothesize that this is due to multiple reasons. First, artifacts in the time series prevail due to steep and rugged terrain and an imperfect co-registration between Landsat and Sentinel-2. Second, there was a mismatch of the reference polygons and the satellite imagery because the polygons were derived with a stereo workstation. This mismatch is difficult to account for, since it is not systematic. Third, the combined usage of ascending and descending orbit Sentinel-1 data required some smoothing and should itself be omitted in the future (see above). However, this smoothing may have removed some (early) signs of bark beetle infestation. Hence, removing it may improve not only the overall detectability, but particularly the early detection of infestations.

Our approach grants users some control over the detection timing. While it does target the detection during the green attack stage, it does not specifically relate to a specific attack stage, in contrast to other, similar approaches (Abdullah et al. 2019; Bárta, Lukeš, and Homolová 2021). Instead, the infestation probabilities based on healthy and gray attack spruce increase overtime and allow users the selection of an appropriate probability threshold. Hence, it can be employed for the continuous monitoring of bark beetle infestations and, as it is not aimed at a particular attack phase nor based on a change in a specific spectral trait (Lausch et al. 2016), in the detection and monitoring of other forest disturbances as well, if sufficient training data are available. If the focus is on the rapid detection of bark beetle infestations, other approaches that would also work in NRT may be more appropriate as well and should be tested, e.g. in which a baseline model of each pixel’s phenology or intra-annual variability is fit and then used as the basis for comparisons with newly acquired reflectance/index values (König et al. 2020; Löw and Koukal 2020; Mathieu et al. 2022). A possible improvement without completely altering our approach would be the use of minimum instead of median composites as the reference for infested areas; in this case the distance to index values representing healthy forests would likely be greater.

Our study showed that Sentinel-1’s C-band data are not well-suited to detect bark beetle infestations (although they have been successfully used to detect logging activity; Reiche et al. 2018), in contrast to SAR data with longer wavelengths (Tanase et al. 2018). However, there are currently no active SAR satellites offering the same continuous revisits as Sentinel-1, which means that SAR data cannot provide sensitive inputs for operational/NRT monitoring of bark beetle infestations: either the data are not suitable (Sentinel-1) or too few images are available to guarantee a basic time series density able to capture the different infestation stages (e.g. ALOS PALSAR; Reiche et al. 2018). New SAR satellites with L-band capabilities and consistent, frequent revisits, such as NISAR, would benefit forestry applications, especially studies of bark beetle infestation in areas with a high cloud cover probability.

Conclusions

In this study, we used full time series of multi-sensor satellite data, obtained from Landsat, Sentinel-2 (multispectral), and Sentinel-1 (SAR) to map bark beetle infestations in a mountainous forest in Central Europe and evaluated their respective performance. We applied a consistent evaluation procedure which made the comparison of the suitability of these different sensors possible. Additionally, as a new technique in this context, the combination of Landsat and Sentinel-2 as well as all three sensors (i.e. Landsat, Sentinel-1, and Sentinel-2) was tested.

Our results showed that Sentinel-2 is best equipped for mapping bark beetle infestations, based on its high spectral sensitivity to bark-beetle-induced changes in spectral traits and especially due to its red edge bands, which have not yet commonly been used for this task. Sentinel-2 also had a distinct advantage over Landsat, which when combined with Sentinel-2 did not provide any benefits, evidenced by the increased time series variability and the lower spectral resolution. Sentinel-2's high spectral resolution, frequent revisits, and relatively high spatial resolution make it an invaluable tool for monitoring bark beetle infestations. The red edge capabilities of Sentinel-2 should be further explored, such as by testing other infestation detection algorithms based on red edge vegetation indices. Other types of satellites with red edge capabilities, notably, the newer generation of PlanetScope, which also includes a yellow spectral band that may be useful for infestation detections, e.g. if combined with the red edge band (Planet, Inc 2022), will provide valuable information as well. Sentinel-1's high spatiotemporal resolution did not compensate for its limited suitability for forestry applications. Satellites offering similar resolution and cloud cover-independence but which are equipped with sensors working with longer wavelengths would greatly improve operational bark beetle infestation detections.

The high overall accuracies and low commission errors of Sentinel-2, especially in larger infestation patches, further recommend its use. In addition, our threshold-based classification allows users to select both a probability threshold based on the particular application and the optimal tradeoff between low commission errors and early detection. Hence, our approach can be operationally employed in NRT scenarios and in the inventory of bark beetle infestations and other forest disturbances.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author, S.K., upon reasonable request.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2226515

Additional information

Funding

This research was funded by the German Federal Ministry of Transport and Digital Infrastructure through the FirSt 2.0 project under grant number 19F2127D.