Monitoring Bark Beetle Forest Damage in Central Europe. A Remote Sensing Approach Validated with Field Data

: Over the last decades, climate change has triggered an increase in the frequency of spruce bark beetle ( Ips typographus L.) in Central Europe. More than 50% of forests in the Czech Republic are seriously threatened by this pest, leading to high ecological and economic losses. The exponential increase of bark beetle infestation hinders the implementation of costly ﬁeld campaigns to prevent and mitigate its e ﬀ ects. Remote sensing may help to overcome such limitations as it provides frequent and spatially continuous data on vegetation condition. Using Sentinel-2 images as main input, two models have been developed to test the ability of this data source to map bark beetle damage and severity. All models were based on a change detection approach, and required the generation of previous forest mask and dominant species maps. The ﬁrst damage mapping model was developed for 2019 and 2020, and it was based on bi-temporal regressions in spruce areas to estimate forest vitality and bark beetle damage. A second model was developed for 2020 considering all forest area, but excluding clear-cuts and completely dead areas, in order to map only changes in stands dominated by alive trees. The three products were validated with in situ data. All the maps showed high accuracies (acc > 0.80). Accuracy was higher than 0.95 and F1-score was higher than 0.88 for areas with high severity, with omission errors under 0.09 in all cases. This conﬁrmed the ability of all the models to detect bark beetle attack at the last phases. Areas with no damage or low severity showed more complex results. The no damage category yielded greater commission errors and relative bias (CEs = 0.30–0.42, relB = 0.42–0.51). The similar results obtained for 2020 leaving out clear-cuts and dead trees proved that the proposed methods could be used to help forest managers ﬁght bark beetle pests. These biotic damage products based on Sentinel-2 can be set up for any location to derive regular forest vitality maps and inform of early damage.


Introduction
In the last three centuries, the forests in the Czech Republic have changed dramatically. Triggered by economic reasons, natural mixed forests were replaced by even-aged spruce or pine monocultures. The spruce currently accounts for 49.2% of stands and pine does for 15.9% [1]. These stands are less stable and more threatened by abiotic factors (e.g., drought, wind, snow, glazed frost) and, more indices have proven to be the most effective to detect bark beetle effects at the green attack phase. The correlation between satellite-derived vegetation indices and areas affected by bark beetle at early stages has been tested with Sentinel-2 and Landsat-8 [22] and Landsat-5 images [19]. These studies have proven that changes in photosynthetic activity, vegetation water content, and leaf area index (LAI) are correlated with bark beetle infestation at green attack phase, especially with Sentinel-2, yielding a user's accuracy of 67% for green attack mapping. Since 1986, a regular forest health assessment has been performed in the Czech Republic using the systematic network of the International Cooperative Programme on Assessment and Monitoring of Air Pollution Effects in Forests (ICP Forests). This network is based on core monitoring sites of 16 × 16 km and selected areas of 8 × 8 km, with a total of 306 areas. At altitudes from 150 m to 1100 m, approximately 11,000 trees are evaluated every year, representing 28 species of trees in different age classes. In addition to ground-based evaluation, Landsat-5 satellite images have been used since the mid-1980s, allowing a uniform systematic evaluation unbiased by a subjective factor in ground visual assessment [31]. Since 2016, satellite images from the Sentinel-2 satellites have been analysed at one-year intervals. Only cloudless images of the Czech Republic are used, selected during the phenological summer (i.e., June to August). Leaf area index (LAI) maps are derived from Sentinel-2 images, which have been validated with LAI ground surveys and defoliation values from ICP Forests data. The health status of the stands is assessed on the basis of LAI changes, which significantly correlate with defoliation in the selected time interval [32,33]. For more detailed mapping of bark beetle outbreaks at higher spatial resolution and more frequent updating, a monitoring system to detect salvage cutting and standing dead wood using Planet images has been used since 2018 in the Czech Republic. The following main layers are used to build this system: (1) Areas dominated by spruce and higher than 12 m. The forest tree map is created by the classification of Sentinel-2 images, performed on the basis of the spectral response of tree species in different phenological phases using training data collected during the National Forest Inventory (NFI) ground survey. The Planet system also uses a digital surface model (DSM) derived from aerial images taken in 2016 and 2017. The DSM is normalized using digital terrain model derived from airborne laser scanning to obtain a canopy height model. (2) Map of LAI decrease, which detects either the clear-cuts or dead wood. Normalized difference vegetation index (NDVI) images at 3 m spatial resolution and two additional categories are detected: standing dead forest and newly established clear-cuts. These categories are distinguished based on the Triangular Greenness Index (TGI) [34]. Different Czech institutions are also conducting research into the use of multispectral and hyperspectral cameras carried by unmanned aerial vehicles (UAVs) [25,35,36]. These systems can provide images at local scales almost anytime with very high spatial resolution, which allows the evaluation of health at the level of individual trees. The disadvantage of UAVs compared with high resolution satellites such as Sentinel-2 or Landsat is the cost of data acquisition and the complexity of data processing, although this disadvantage is gradually disappearing with the development of cloud computing. Another shortcoming is that most multispectral cameras are used in UAV-based studies, which allow the calculation of only a few vegetation indices (e.g., NDVI or Normalized Difference Red Edge Index-NDRE) [37]. Hyperspectral cameras are more expensive, and data processing is more complicated. In contrast, Sentinel-2 provides free, easy to process data with a short revisit time, thus allowing more regular estimates about vegetation status.
In this study, a multi-temporal regression-based change detection method is proposed to map areas affected by bark beetle with different degree of severity at 10 m spatial resolution. Moreover, a comparison is carried out between forest vitality estimates and field records of forest damage caused by Ips typographus L. infestation.
The aim of this study is to test a methodology for the automatic mapping of bark beetle-induced damage using Sentinel-2 images. The specific objectives of this work are (a) to evaluate the performance of the damage detection algorithm, (b) to study the feasibility of discriminating areas affected with different intensity or level of mortality, and (c) to assess the influence of possible clear-cuts not related to bark beetle in the final damage maps.

Study Area
The study site was located in a forest area owned by Mendel University in Brno and managed by the University Forest Enterprise (UFE), in Křtiny, Czech Republic (49.170 • N, 16.440 • E. Figure 2). This property covers around 10,000 ha of forest stands. Mixed forests are prevailing where the proportion is 38% of conifers (19% spruce, 8% pine, 8% larch, and others) and 62% of broadleaves (33% beech, 15% oak, 8% hornbeam, and others). Altitude ranges between 200 and 570 m. Mean annual rainfall is 610 mm and mean annual temperature is 7.5 C. UFE forest management is focused on close-to-nature methods, using minimum clear-cuts and a high share of natural regeneration. In the Czech Republic, salvage cutting has increased dramatically in the last three years. Figure 3 shows the total wood production in millions of cubic meters. Table 1 shows salvage cutting volume at the UFE forest.  A total of 32.58 million cubic meters of raw wood was harvested in the forests of the Czech Republic in 2019, which means a further increase of 6.89 million m 3 compared with the previous year. The salvage cutting of 30.94 million m 3 of wood contributed significantly to this volume. The share of salvage cutting in 2019 was 95%. Thus, the initial conditions for planned forest management continued to deteriorate. In terms of the composition of harvesting by tree species, the volume of coniferous wood harvested increased by 7.1 million m 3 compared with 2018 to a total of 31.31 million m 3 . The share of coniferous wood harvesting in total harvesting was approximately 96%. The proportion of deciduous and coniferous wood harvesting is due to salvage cutting, especially the wood infested by bark beetle (Figure 4).
In 2019, 20.7 million m 3 of harvested spruce wood infested by bark beetle was registered in the Czech Republic, which represents an increase compared with 2018 by more than 70%, when approximately 12 million m 3 was recorded (2017-5.34 million m 3 ) ( Figure 5). Harvested wood in 2019 was practically exclusively infested with European spruce bark beetle (Ips typographus L.), which is usually accompanied by six-toothed spruce bark beetle (Pityogenes chalcographus L.) and double-spined bark beetle (Ips duplicatus Sahlberg), especially in northern and central Moravia and Silesia, but locally often elsewhere.  The massive spruce infestation with bark beetles has been triggered by higher average temperatures and lower precipitation over the last five years compared with the long-term normal. Extreme drought prevailed, especially in 2015 and 2018 ( Figure 6).

Sentinel-2 Satellite Images
In this study, we used Sentinel-2 images provided by Copernicus Open Access Hub to detect bark beetle damage. These images were acquired in format Level-2A (bottom of atmosphere reflectance). Sentinel-2 L2A images are generated by applying atmospheric and topographic correction algorithms to the Level-1C (top of atmosphere reflectance) images [40]. In order to focus the change detection in areas prone to be attacked by the spruce bark beetle, forest/non-forest and dominant species maps were produced for the study area (Figure 7). These first maps were also developed using Sentinel-2 images ( Table 2).  Different Sentinel-2 bands (Table 3), vegetation indices, and texture indices were used for each process. The whole process, summarised in Figure 7, is explained in the subsequent sections.

Ground Truth Data
Records of salvage cutting and records of clear-cuts were used to build the ground truth dataset, together with a forest stand map derived from forest management plan (valid for the period 2013-2022). The records were collected continuously by foresters for the needs of the forest enterprise and state forest administration, and contained the ID of forest stands, month and year of the cutting, volume, type of cutting, cause, and area. A detailed explanation of ground truth generation can be found in Section 3.4.

Methods
Two algorithms were developed to map bark beetle damage: biotic damage A (BDA) and biotic damage B (BDB), considering as damage any negative deviation from the standard development of the forest (i.e., clear decrease of forest health). Both algorithms work by detecting changes between two dates. The BDA algorithm included masking Sentinel-2 surface reflectance images to areas dominated by spruce in the starting date of the pair. Hence, any area with dead trees or clear-cuts might be detected as bark beetle damage. In the case of BDB, satellite images were masked to the whole forest area in the final date of the pair (not in the initial one). Thus, all dead trees and clear-cuts could not be classified as damaged areas by the BDB algorithm. The objective of the difference between BDA and BDB was to test the ability of the algorithms to detect damage in non-cut areas and with most of the trees still alive. This allowed us to assess the true performance of our approach to detect the effects caused by bark beetle, as BDB leaves out possible over-detection caused by clear-cuts. Three different layers were ultimately developed (Table 2, Figure 7 Step 4): (1) the BDA layer using 2019 imagery (BDA19), (2) the BDA layer using 2020 imagery (BDA20), and (3) the BDB layer using 2020 imagery (BDB20) to map only changes in stands dominated by alive trees. The algorithms for BDA19 and BDA20 were identical, the only difference being the dates of the input imagery (2018-2019 and 2019-2020, respectively). The common steps for producing all damage maps from Sentinel-2 images are summarised in the flowchart ( Figure 7) and explained in detail in the following subsections.

Forest/Non-Forest Classification
To derive any of the damage detection layers explained above (Figure 7, Step 4.), we first apply a classification algorithm to discard all the pixels belonging to non-forest cover (Figure 7, Step 2.1.) [41,42]. To represent phenological variations in the forest canopy, two Sentinel-2 images were used, one for winter and one for summer (Table 2). Green, red, NIR, and SWIR bands were selected. SWIR bands, originally at 20 m spatial resolution, were resampled to the resolution of the other bands, 10 m. Normalized difference vegetation index (NDVI, Equation (1)) [43] and modified soil adjusted vegetation index (MSAVI, Equation (2)) [44,45] were computed together with two texture indices: homogeneity (Equation (3)) and entropy (Equation (4)) [46]. A supervised classification method based on random forest algorithm [47] was used to classify pixels in forest and non-forest categories with the selected Sentinel-2 bands and the vegetation and texture indices as input. The output of the classification was a binary forest/non-forest map of the study area at 10 m spatial resolution. A validation protocol was specifically designed to carry out an independent validation of the forest mask [42], based on a stratified random sampling strategy. According to this validation, the forest mask had an accuracy of 0.97 for the study area in 2018. This process was repeated to update the forest mask to 2019 and 2020.
Remote Sens. 2020, 12, 3634 9 of 19 where n is the number of grey levels and P(i, j) defines the entries of the grey-level co-occurrence matrix.

Species Identification
The next step (Figure 7, Step 2.2.) was to identify the areas dominated by spruce (i.e., those that are prone to be attacked by the Ips typographus L.), for which a dominant tree species classification was carried out [48]. The bi-temporal Sentinel-2 bands and indices used in previous steps were masked with the forest/non-forest map, selecting only the forest pixels. A random forest model was trained with polygons containing the main tree species present in the area, namely, Carpinus betulus L., Fagus sylvatica L., Larix decidua Mill., Picea abies (L.) Karst., Pinus Sylvestris L., and Quercus robur L. Training samples were obtained joining forest stands with information from the local forest management plan. For each of the species considered, stands with an abundance of 70% or higher were pre-selected. Finally, the smaller stands that actually contain the specific tree species were selected for training.
The species attacked by Ips typographus L. (i.e, Picea abies (L.) Karst) were classified with an accuracy of 0.91 in 2018 according to the cross-validation performed during the model training. This dominant species map was produced for 2018, 2019, and 2020, as it was made with the forest mask.

Anomaly Detection and Damage Detection
Once spruce areas were clearly delimited for each year, a new method (Figure 7, Step 3) was proposed to map the bark beetle damage suffered by the spruce forests during a year (i.e., between a pair of dates). The first step was to compute vegetation indices able to represent vegetation condition for the two dates (t 0 and t 1 for each year, in Table 2). Based on the literature review about the effects of bark beetle [49], the following vegetation indices were computed for each date: NDVI (Equation (1)), MSAVI (Equation (2)), normalized difference moisture index (NDMI, Equation (5)), and green leaf area index (LAI green , Equation (6)) [50,51].
NDVI and SAVI are direct indicators of plant greenness and photosynthetic activity, NDMI is an indicator of vegetation water content, and LAI green is an estimate of the proportion of green leaves per area. As defoliation is one of the main effects caused by bark beetles [15], LAI has been considered as one of the key variables in this study. Moreover, we proposed the use of Sentinel-2-derived LAI green , as it considers only the area of green leaves, not only being an indicator of defoliation, but also of leaf browning [50,51].
To derive the BDA19 and BDA20 products, vegetation indices were masked with the spruce polygons identified in the dominant tree species map. For each index, a bi-temporal ordinary least squares (OLS) regression was carried out in order to model the typical behaviour of spruce forest patches between t 0 and t 1 . The first index was used as the independent variable, and the second one as the dependent variable in each regression. The results were images of estimated NDVI, MSAVI, NDMI, and LAI green for t 1 based on the average trend of spruce forests in the studied period. An error image was computed for each index using the estimated and real values of the indices in t 1 . These images show the deviations from the expected condition of healthy spruce masses, as it was assumed that regression would be adjusted to healthy forests and disturbed areas would appear as deviations from that ideal status. Hence, these error images showed estimates of temporal decrease or increase of photosynthetic activity, vegetation water content, and number of green leaves.
In order to obtain a single image of vegetation vitality changes (Figure 7, Step 4.1.), all error images were standardised (i.e., pixel values expressing differences between real and estimated indices were converted to z-scores). Z-scores (Equation (7)) represent the number of standard deviations between a value and the population mean.
where µ is the mean of the population and σ is the standard deviation of the population. The mean of all the standardised images was computed to obtain a single image of average changes in vegetation vitality. Pixel values in this image may be interpreted as mean standard deviations from expected NDVI, MSAVI, NDMI, and LAI green of undisturbed spruce forest. Finally, these vitality changes were reclassified (Figure 7, Step 4.2.) to obtain a categorical map of bark beetle damage (Figure 7, final outputs, Step 4.3.), using the following rule based on z-scores: 0-1 no damage, 1-2 minor damage, 2-3 moderate damage, and >3 severe damage.

Ground Truth Data and Validation
To validate the damage products, a ground truth dataset was generated through a field campaign. An equalized stratified random sampling was used for data collection. As field works are costly, a limit of 200 points per product was established. To stratify the samples, 50 random points were generated per class and product, obtaining 200 points per product and 600 in total.
The collected records containing the ID of forest stands, month and year of the cutting, volume, type of cutting, cause, and area were joined to the forest stand map using the forest stand ID. Different classes of forest stands representing the degree and year of damage were created using the actual vitality of the stands based on the information of the records. This modified forest stand map was used to decide to which class each point belongs in terms of actual damage. This decision was made based on the occurrence of damage in the given forest stand, not on the occurrence on the exact position of the point, as the records were focused only at the level of the given forest stand. As dead trees with inactive bark beetle are no longer harvested in 2019 and 2020, in the case of doubt (no record of harvesting, but occurrence of damage detected), a field evaluation was also performed.
An accuracy assessment was carried out based on the analysis of the confusion matrix derived from the damage maps and ground truth points. The overall accuracy was computed for each product. In order to obtain performance metrics per class, each class was binarised against the rest, considering as positive the category of interest and negative the sum of all other classes. The following performance metrics were computed per product and damage category: accuracy (Equation (8)), precision (Equation (9)), recall (Equation (10)), F1-score (Equation (11)), omission error (Equation (12)), commission error (Equation (13)), and relative bias (Equation (14)) [42].

acc. =
True positive + True negative sample size (n) (8) prec. = True positive True positive + False positive As an additional exercise, each point of the ground truth dataset was characterised with the vitality value of the corresponding year. The objective of this procedure was to study the statistical robustness of each class and detect possible overlaps or other inconsistences across the categories defined in the final maps.

Bark Beetle Damage Maps
The BDA product was developed for 2019 ( Figure 8) and 2020, while the BDB was produced only for 2020. Each product had two outputs: (1)

Validation Results
The three categorical maps were validated following the methods explained in Section 3.4. Confusion matrices (Figure 9) show that the algorithm performed as expected, as errors were low compared with the agreement between the classification maps and ground truth. Errors seemed to be slightly higher for the BDB20 product. Overall accuracies were higher than 0.80 for all products (accBDA19 = 0.90, accBDA20 = 0.85, accBDB20 = 0.81). Metrics per class (Table 4) confirmed the good performance of the algorithms. BDA products showed a similar performance for different years, with accuracies increasing with the severity of the attack. Accuracy and F1-score were higher for the severe damage class (acc19 = 0.98, F1-19 = 0.96, acc20 = 0.95, F1-20 = 0.88), thus confirming that most affected areas were easier to detect. Nevertheless, cut areas might have influenced this result as they were classified into the severe damage class. On the other hand, the no damage class yielded lower accuracy and F1-score (acc19 = 0.90, F1-19 = 0.83, acc20 = 0.86, F1-20 = 0.76). The BDA20 product showed slightly worse results for all classes (acc no-damage = 0.81, F1 no-damage = 0.70, acc severe-damage = 0.96, F1 severe-damage = 0.90). Error metrics showed similar behaviour, being higher for BDB20 and lower for BDA19. Commission errors ranged from 0.08 (BDA19, severe damage) to 0.26 (BDB20, minor and moderate damage) and were higher than omission errors in all damage classes. The negative relative bias in these categories confirmed that the algorithm tends to over-estimate damaged areas, although values are not greater than 0.22. Omission errors were in general very low, between 0.00 and 0.08 for all damage classes in the three maps. The high recall in these classes (0.93-1.00) confirmed that the algorithms are able to detect affected areas, missing a small proportion of positive cases. On the other hand, omission errors were significantly higher than commission errors in the no damage class, leading to high positive relative bias ranging from 0.46 to 0.56. This was another signal of over-estimation of damaged areas.
Box plots of forest vitality per class ( Figure 10) confirmed the clear correlation between the actual areas attacked by the bark beetle and the decrease of forest vitality, measured as a standardised mean of changes in NDVI, MSAVI, NDMI, and LAI green . The interquartile range (IQR) represents the data dispersion and is calculated with the difference between 75th and 25th percentiles. IQR is greater in no damage (IQR BDA19 = 1.49, IQR BDA20 = 2.22, and IQR BDB20 = 2.63) than in the rest of the categories for the three products analysed. The overlap between categories decreases with the increase of severity, with moderate and severe damage not overlapping any other class in the three products. Outliers are present in the no damage, minor damage, and moderate damage categories, with severe damage being the only class without outliers in the three maps. This confirms that areas with severe damage are easier to detect and identify than the other classes, as inferred from performance metrics ( Table 4).
these categories confirmed that the algorithm tends to over-estimate damaged areas, although values are not greater than 0.22. Omission errors were in general very low, between 0.00 and 0.08 for all damage classes in the three maps. The high recall in these classes (0.93-1.00) confirmed that the algorithms are able to detect affected areas, missing a small proportion of positive cases. On the other hand, omission errors were significantly higher than commission errors in the no damage class, leading to high positive relative bias ranging from 0.46 to 0.56. This was another signal of overestimation of damaged areas.
Box plots of forest vitality per class (Error! Reference source not found.) confirmed the clear correlation between the actual areas attacked by the bark beetle and the decrease of forest vitality, measured as a standardised mean of changes in NDVI, MSAVI, NDMI, and LAIgreen. The interquartile range (IQR) represents the data dispersion and is calculated with the difference between 75th and 25th percentiles. IQR is greater in no damage (IQRBDA19 = 1.49, IQRBDA20 = 2.22, and IQRBDB20 = 2.63) than in the rest of the categories for the three products analysed. The overlap between categories decreases with the increase of severity, with moderate and severe damage not overlapping any other class in the three products. Outliers are present in the no damage, minor damage, and moderate damage categories, with severe damage being the only class without outliers in the three maps. This confirms that areas with severe damage are easier to detect and identify than the other classes, as inferred from performance metrics (Error! Reference source not found.). Analysing the three products separately, some differences were noted. Overlap of no damage and minor damage between the 25th and 75th percentiles was clearer in the BDB20 product. In this product, there was also overlapping between no damage and moderate damage. However, the severe damage class was clearly discriminated from no damage and minor damage. This behaviour might suggest that a more precise discrimination could be achieved by reducing the number of classes.

Discussion
The results showed that the biotic damage algorithm worked consistently in different years. However, accuracy decreased by 5% between BDA19 and BDA20, while in these years, the proportion of spruce detected as damaged was very similar (7.2% and 7.4% of the total spruce area, respectively). Hence, although the results for both years were similar regarding performance metrics, a longer time series analysis would be needed to study the robustness of the method across time and quantify typical inter-annual variations. In the case of BDB20 (i.e., similar to BDA20, but considering areas dominated by alive trees of any species at the end of 2020, not only spruce areas), the lower accuracy (accHD20 = 0.81) may be explained by two main factors. Firstly, this product considers only alive trees, and minor and moderate changes have greater weight than severe changes compared with the other products. Hence, the overall metrics tend to be lower in those products where complex intermediate classes are more frequent than extremes that are easier to identify. Secondly, all species are considered within the BDB product, and not only spruce, subtle changes in areas dominated by Analysing the three products separately, some differences were noted. Overlap of no damage and minor damage between the 25th and 75th percentiles was clearer in the BDB20 product. In this product, there was also overlapping between no damage and moderate damage. However, the severe damage class was clearly discriminated from no damage and minor damage. This behaviour might suggest that a more precise discrimination could be achieved by reducing the number of classes.

Discussion
The results showed that the biotic damage algorithm worked consistently in different years. However, accuracy decreased by 5% between BDA19 and BDA20, while in these years, the proportion of spruce detected as damaged was very similar (7.2% and 7.4% of the total spruce area, respectively). Hence, although the results for both years were similar regarding performance metrics, a longer time series analysis would be needed to study the robustness of the method across time and quantify typical inter-annual variations. In the case of BDB20 (i.e., similar to BDA20, but considering areas dominated by alive trees of any species at the end of 2020, not only spruce areas), the lower accuracy (accHD20 = 0.81) may be explained by two main factors. Firstly, this product considers only alive trees, and minor and moderate changes have greater weight than severe changes compared with the other products. Hence, the overall metrics tend to be lower in those products where complex intermediate classes are more frequent than extremes that are easier to identify. Secondly, all species are considered within the BDB product, and not only spruce, subtle changes in areas dominated by deciduous species, which are also more sensitive to climate variations (e.g., drought), might have led to a greater confusion between the no damage and damage classes, as shown in the confusion matrices ( Figure 9). However, the accuracy of the BDB map was high and, more importantly, it confirmed that the high accuracy obtained for the first two products was not due to the detection of cut areas, thus damage caused by bark beetle was clearly detected in most cases.
The three models presented higher commission than omission and negative relative bias in all damage classes, suggesting a certain degree of over-estimation of damaged areas, as confirmed by higher recall and lower precision. Assuming that a precise balance between commission and omission errors was hard to achieve, we prioritised the commission of damaged areas, as the aim of the products is to help forest managers to optimise the control of the pest. Although this led to high positive bias in the no damage class (i.e., there was a relatively high omission), accuracy and F1 were still high for this category. For BDB20, the no damage class presented the worst results with a recall of 0.58 and an omission error of 0.42. This confirmed that further improvements should be made to avoid the confusion between bark beetle damage and minor changes in forest vitality caused by other factors across the different forest types present in the study area. Different solutions might be applied in order to balance the errors of the no damage category, such as adding more spectral indices or change the temporal windows. Time series could allow carrying out the damage assessment using a combination of shorter periods, thus minimising the possible influence of slow changes in vegetation cover.
The results clearly showed that the minor damage class was often confused with no damage, as the limit between small regular changes in vegetation condition and small, but anomalous changes is hard to define. This suggests that the models would improve significantly if the no damage category was removed. Nevertheless, we decided to leave this class because it might be directly related to green-attack detection [22], which would be essential to design an early detection system [52].
The approach proposed in this article may improve previous attempts to detect bark beetle infestation with Sentinel-2 data. The ability of Sentinel-2 for this purpose has been tested by different authors with variable results. Abdullah et al. [22] reported a user's accuracy of 67% (i.e., agreement between pixels detected as infestation and ground truth) and a producer's accuracy of 71% for green attack detection, while our user's accuracies were not lower than 74% and our producer's accuracies not lower than 95% (CE = 0.26 and OE = 0.05 in the worst model, BDB20) for the minor damage category, which, as stated before, might correspond to green attack. Moreover, our approach was validated with field data, while Abdullah et al. [22] derived ground truth information from the visual analysis of aerial photography, which might have added a certain level of uncertainty into the accuracy assessment. Both studies coincide in the affirmation that changes in photosynthetic activity, green leaves, and humidity triggered by bark beetle can be detected from early infestation stages using multispectral Sentinel-2 data. The relationship between LAI maps derived from Sentinel-2 and defoliation caused by bark beetle infestation has also been confirmed by Barka et al. [33] (r 2 = 0.58). An advantage of the methods proposed in this study compared with the previously mentioned ones is that it does not require training data to detect bark beetle infestation, as the algorithms used unsupervised methods based on a statistical change detection approach. Zimmermann and Hoffmann [23] used a change detection approach to detect areas affected by Ips typographus in Germany and Switzerland. Although they obtained low commission errors for the positive class (CE = 0.03-0.12 versus CE = 0.08-0.26 of our algorithm), omission errors were substantially higher than in our study (OE = 0.48-0.60 versus OE = 0.00-0.08), leading to a high bias in the final results.
Regarding very high resolution, Worldview-2 showed potential to discriminate green attack from healthy trees [53]. However, the spectral differences were so small that the accuracy of random forest classification and logistic regression did not exceed 70%. RapidEye data have been used to map green attack, nearly equivalent to our minor severity class, reaching a kappa coefficient of 0.51 [54]. In comparison with these studies, our performance metrics yielded higher values, with an accuracy of 93% for the worst model. Nevertheless, it must be noted that the spatial resolutions of multispectral Worldview-2 (1.84 m at nadir) and RapidEye (6.5 m at nadir) allow the detection of smaller groups or even individual trees when compared with Sentinel-2. Although bark beetle damage mapping has been carried out using UAV very high resolution imagery with promising results [24,[27][28][29], they rely on planned flights, which usually require high costs and cannot be repeated frequently. The value of using Sentinel-2 data is the free cost of imagery at relatively high spatial resolution (10 m), but mainly the revisit time [15] of 2-3 days at mid-latitudes. Relatively high performance metrics have also been reported using SAR data, with accuracies of 0.65-0.88 for L-band data [30] and kappa of 0.23 for green attack detection using X-band [55]. However, SAR-based detections were highly dependent on local environmental conditions [30]. The combination of very high resolution multispectral RapidEye and TerraSAR-X data yielded a kappa coefficient of 0.74, considerably improving the results based solely on optical and SAR data [55]. The results shown in this study are likely to improve if SAR data (e.g., Sentinel-1) are included in the models.
Using remote sensing products, we also need to consider that tree damage stages do not always correspond to the colour changes in the crown. In extreme cases, it occurs that the tree bark is totally peeled, the bark beetle has completed its development and flew out, while the crown is still green. This fact complicates the damage assessment with remote sensing optically derived products.
The next point of discussion focused on the local expertise of how the products developed in this study overcome some of the limitations of other datasets currently provided at the national level, as is the case of Czech Republic. Czech Hydrometeorological Institute provides a map of the sums of effective temperatures above 7.5 • C, knowing that the population dynamics of bark beetle are significantly affected by the weather, especially air temperatures, which determine the swarming and the number of generations. As the bark beetle swarms when temperatures above 7.5 • C sum up to 540 • C [56], this temperature map is intended to help forest managers in planning and prevention activities by estimating swarming dates. Even when the temperatures map is somewhat useful for managing spruce forests, there are other factors affecting the ecology of bark beetle that are not considered within this map. Moreover, it does not have the intention of being used for damage assessment purposes.
Next, also based on local expert knowledge, we assessed the monitoring system for the detection of salvage cutting and standing dead wood using Planet images established in the Czech Republic since 2018 (see Section 1) in comparison with the products produced and validated here. The Planet-based dataset is updated four times a year [57]. Although this product is useful to control dead trees and subsequent clear-cuts caused by bark beetle, it does not provide any information about the status of trees in intermediate infestation phases (i.e., it is equivalent to a binary classification detecting only dead trees). Contrarily, our bark beetle damage maps provide insights on the level of affection of different areas, even those dominated by trees that are still alive, as proven with BDB20. This suggests that the proposed method is able to improve Planet-based very high resolution maps by providing estimates of areas under red and green attack phase. Moreover, the costs of maps based on Sentinel-2 is lower because of the free availability of imagery. Nevertheless, the methods proposed here are limited to mapping relatively large areas, as the spatial resolution of 10 m does not allow detecting single trees affected by bark beetles. Hence, it is necessary that several trees within the pixel present enough spectral change to be detected as anomalous. The same limitation was found by other authors [23], still making these kind of methods useful for intermediate landscape scales.
Further research is still to be made on the field of bark beetle attack forecasting models, which are being highly demanded by forest managers [52]. Although the methods proposed here could be used for short-term forecasting because of their ability to detect areas with minor damage that might be in the first infestation stages, there is no evidence of a clear correlation between these areas and green attack phase. Moreover, future models should include climate data as input in order to perform mid-term forecasts allowing to mitigate the negative effects of this pest before its complete development.

Conclusions
In this study, a method was proposed to map bark beetle damage from satellite imagery. Using Sentinel-2 images as input, two multi-temporal regression models were built to detect and map the severity of bark beetle outbreaks on spruce forests in the Czech Republic at 10 m spatial resolution. The first model (BDA) was applied to map the damage occurred in years 2019 and 2020. The second model (BDB) was applied to map the bark beetle damage in 2020, in order to map only changes in stands dominated by alive trees. Both products were validated using a ground truth dataset generated through a field campaign and forest management plans. Different performance metrics were computed to assess the quality and errors of the three maps produced. Finally, forest vitality used to develop the biotic damage layers was compared against ground truth to study the overlaps between classes.
All products showed good performance, with accuracies higher than 0.80. The severe damage class yielded the best performance metrics in all products (acc > 0.95, F1 > 0.88), while the no damage yielded the worse metrics (acc > 0.81, F1 > 0.70). Commission errors were higher than omission errors in all positive damage classes, leading to a high relative bias in the no damage class (relB = 0.42-0.46). BDA products showed slightly better results than BDB (accBDA19 = 0.90, accBDA20 = 0.85, accBDB20 = 0.81) because of the easier identification of areas with clear-cuts or dominated by dead trees. However, the performance metrics yielded by BDB proved that the algorithm was able to identify areas affected with low severity unaffected by dead trees or clear-cuts, suggesting that some areas were correctly mapped at red and green attack phases. Changes in forest vitality grouped by ground truth classes confirmed that pixels were easier to classify in the correct class for the severe and moderate damage classes. Contrarily, the existing overlap of vitality values between classes of no damage and minor damage highlighted the difficulty of clearly discriminating infested areas with subtle signals of decay.
Comparing the proposed methods and outputs with the datasets currently used in the Czech Republic for bark beetle damage mapping, the products presented several advantages. Firstly, the cost is relatively low as it is based on freely available Sentinel-2 images. Secondly, the biotic damage maps provide information about damage intensity, suggesting that it might be used not only for damage assessment, but also to help forest managers in planning their prevention and mitigation activities. Biotic damage products can be set up for any location to monitor the forest vitality to derive regular maps as needed, for example, every month. Nevertheless, the presented methods are not valid to identify individual affected trees given their spatial resolution. Future research should be carried out to confirm the detection at the green attack phase and complement the existing studies with forecasting products.