Assessment of Sentinel-2 Images, Support Vector Machines and Change Detection Algorithms for Bark Beetle Outbreaks Mapping in the Tatra Mountains

Abstract

Cambiophagous insects, fires and windthrow cause significant forest disturbances, generating ecological changes and economical losses. The bark beetle (Ips typographus L.), inhabiting coniferous forests and eliminating weakened trees, plays a key role in posing a threat to tree stands, which are dominated by Norway spruce (Picea abies) and covers a large part of mountain areas, as well as the lowlands of Northern, Central and Eastern Europe. Due to the dynamics of the phenomena taking place, the EU recommends constant monitoring of forests in terms of large-area disturbances and factors affecting tree stands’ susceptibility to destruction. The right tools for this are multispectral satellite images, which regularly and free of charge provide up-to-date information on changes in the environment. The aim of this study was to develop a method of identifying disturbances of spruce stands, including the identification of bark beetle outbreaks. Sentinel 2 images from 2015–2018 were used for this purpose; the reference data were high-resolution aerial images, satellite WorldView 2, as well as field verification data. Support Vector Machines (SVM) distinguished six classes: deciduous forests, coniferous forests, grasslands, rocks, snags (dieback of standing trees) and cuts/windthrow. Remote sensing vegetation indices, Multivariate Alteration Detection (MAD), Multivariate Alteration Detection/Maximum Autocorrelation Factor (MAD/MAF), iteratively re-weighted Multivariate Alteration Detection (iMAD) and trained SVM signatures from another year, stacked band rasters allowed us to identify: (1) no changes; (2) dieback of standing trees; (3) logging or falling down of trees. The overall accuracy of the SVM classification oscillated between 97–99%; it was observed that in 2015–2018, as a result of the windthrow and bark beetle outbreaks and the consequences of those natural disturbances (e.g., sanitary cuts), approximately 62.5 km² of coniferous stands (29%) died in the studied area of the Tatra Mountains.

1. Introduction

Due to the rapid growth of biomass and low requirements of a habitat, Norway spruce (Picea abies) has played a leading role in forestry in the European mountain areas since the 18th century; for this reason, it was often planted in places inconsistent with its habitat, i.e., where deciduous and mixed forests previously grew [1,2]. The flat and horizontally developed root system of Picea abies prefers a cool and humid climate, not tolerating soil dryness and overheating [1]. However, with the aging of the stand, Picea abies becomes susceptible to attacks by insects (mainly spruce bark beetle Ips typhographus) and fungi [3,4]. In addition, progressive climate changes, manifested by higher air temperatures and reduced rainfall, allows Ips typographus to start reproducing earlier, reaching more generations annually [5]. In the years of 1958–2001, when the climate warmed by 0.13 °C per decade, the growth of forests covered by bark beetle gradations grew on average by 5.31% per year. The most intense dieback of Picea abies trees occurred in the years with low rainfall in spring and summer and high average annual temperature [6]. All these factors indicate that much more frequent and severe outbreaks of cambiophagous insects can be found in spruce forests. Nevertheless, it should be emphasized that the insects are naturally occurring to eliminate weakening spruce trees from their stands [4], and the outbreaks are a natural process that inhibits the expansion of spruce beyond its optimal climatic range [3].

The spruce bark beetle (Ips typographus) is an insect of the Scolytinae subfamily in the Curculionidae family. Adult beetles are 4–6 mm long and 2 mm wide, and are brown–black in color. It can be found all over Europe, in Russia from Siberia to Kamchatka, as well as in China, Korea and Japan. Its distribution is correlated with the range of its main host plant—Norway spruce (Picea abies) and closely related species, because Norway spruce grows only in Europe, and Ips typographus can be found in Asia on other spruce species, e.g., P. obovata, P. orientalis, P. yezoensis [7]. A large number of bark beetles are able to effectively stop the transport of water and nutrients in the spruce, and also introduce fungi that would not otherwise be able to attack trees [3]. For several months (often until winter), infested spruce trees have green needles and show no clear signs of attack [8]; in the next season, most of the needles are lost and the tree is dead, the so-called snags [9].

A forest regeneration capacity is significant; an example can be the Bavarian Forest National Park, where a spruce bark beetle outbreak began in 1993, and in 2005 the inventory indicated the regeneration index oscillated around 4500 20-cm-high trees per hectare, with 89 percent of Picea abies and 7.9% of rowan (Sorbus aucuparia). This shows that montane zones reproduce well and their continuity is not endangered, while the dynamics of the spruce forests in the foothill zone is more complicated [10].

Inventory of dynamically changing stands is a demanding, labor-intensive and costly task due to the irregular network of forest roads. Additionally, standard grids of sampling plots (from 200 × 200 m to 1000 × 1000 m) and the subjectivity of the observed changes reduce monitoring quality [11,12]. This is the reason behind the enormous popularity of remote sensing data, which plays a primary role in inventory species composition [13,14,15]; vegetation parameters, e.g., Leaf Area Index (LAI), Leaf Chlorophyll Content (LCC) [16]; plant water stress [17], leaf defoliation [18] or dieback of trees as a result of insect attack [19]. An example of such research is the analysis of Landsat images to assess mountain pine beetle (Dendroctonus ponderosae) invasions on various pine species in the Rocky Mountain area [20]. The course of the invasion was similar to that of Ips typographus in the stands of Picea abies. The Maximum Likelihood classification of Landsat images and Tasseled Cap transform derivatives (Brightness, Wetness, Greeness) and methods have classified the affected pines with a producer accuracy (PA) of 88.3% and 97.3% user accuracy (UA). The overall accuracy of this classification was 91% and the Kappa coefficient oscillated around 0.90. Meddens et al. [21] also showed that the Sentinel-2 MSI was optimal for assessing the effectiveness of the classification of attacked lodgepole pine (Pinus contorta) and ponderosa pine (P. ponderosa) based on a multi-period analysis of changes, as the classification of attacked pines (over 70% of damaged trees on the studied plot) achieved 92.9% user accuracy and 87.4% produce accuracy, overall accuracy—89.6% and Kappa coefficient—0.86. Considering mixed pixels in which the classes have at least 50% of the pixel share, user accuracy was 81.3%, producer accuracy: 85.3%, overall accuracy: 83.2% and the Kappa coefficient—0.78 [21]. Similar studies were carried out in the Tatra Mountains to detect spruce bark beetle outbreaks [22]; based on the images of Landsat-5 and Landsat-7 (2006–2011), the researchers calculated the vegetation indexes: Normalized Difference Vegetation Index (NDVI), Moisture Stress Index (MSI), Normalized Difference Moisture Index (NDMI), Vegetation Condition Index (VCI) and Tasseled Cap (TC) indicators. It was shown that at the maximum points of gradation, the classification accuracy reached 94% (overall accuracy, and Kappa coefficient—0.70), confirming the usefulness of the MSI [22].

In 2018, Abdullah et al. [8], based on Sentinel-2 and Landsat-8 images, showed that affected tree needles have less chlorophyll and water, a similar amount of nitrogen and more dry matter than those from the undisturbed trees. Then, based on satellite images, chlorophyll indices (e.g., Normalized Green–Red Difference Index) confirmed the changes allowing to identify bark beetle outbreaks with 67% accuracy on Sentinel-2 images [8].

Sentinel-2 images were processed with change detection methods Multivariate Alteration Detection (MAD) and Maximum Autocorrelation Factor (MAF) and with Support Vector Machines (SVM) classification, allowing us to identify bark beetle outbreak in Białowieża Forest (NE Poland) with an overall accuracy of 97% (Kappa coefficient: 0.93; producer accuracy: 92.5%; and user accuracy: 92.0% for attacked trees [23]). Eight World-View-2 spectral bands presenting mountain forests of Austrian Styria and Random Forest classifiers identified bark beetle damages of Norway spruce in different conditions: healthy trees, or even still green, but attacked trees (green attack) and snags. The best results achieved snags (producer accuracy scored 99%, and user accuracy—100%), undisturbed forest and the attacked trees oscillated around 60–70% [24].

Remote-sensing detection of freshly/green-attacked spruce trees by the bark beetle insects is a challenge, but it allows us to react at the beginning of a new outbreak by setting pheromone traps and sanitation, protecting the rest of the area [3,8].

The current paper focuses on methodical and application objectives; in the first case, it is an assessment of the multitemporal Sentinel-2 images (2015–2018) for the identification of coniferous and deciduous forests, grasslands, rocks, snags and cuts or windthrows. For this purpose, the Support Vector Machines (SVM) algorithm was used as an optimal classifier for stand mapping. The studies based on accuracy metrics of Sentinel-2 bands, spectral signatures of individual land cover forms derived from different years and remote-sensing indices; the classification accuracy was compared between changes taking place in coniferous forests, focusing on areas with dieback of trees; cuts or windthrows; or no changes in individual years. The condition of coniferous stands was carried out on the basis of the analysis of the values of spectral features (raster stacking), trained spectral signatures, remote-sensing indices, as well as correlation analyses: Multivariate Alteration Detection (MAD), Multivariate Alteration Detection/Maximum Autocorrelation Factor (MAD/MAF), and iteratively re-weighted Multivariate Alteration Detection (iMAD). The application objectives are to assess the dynamics of changes taking place within coniferous stands, including the outbreak of the spruce bark beetle in coniferous forests for the years 2015–2018 for the whole Tatra Mountains (Polish and Slovak parts). These changes were determined by statistics and cartographic presentations. An important element of the study is reference data based on high-resolution image data (airborne orthophotomaps and multispectral WorldView 2 images), which allowed us to identify analyzed targets on detailed maps, which were field verified by employees of the Tatra National Park [25]. Having such detailed reference maps for the whole area of the Polish Park, the current study allowed us to evaluate the proposed methods. The achievements are not only maps and statistical data about the dynamics of the environmental changes, but above all an assessment of different methods, constituting a signpost for people planning to conduct long-term monitoring of mountain forest areas.

2. Research Area and Targets

The Tatra Mountains are the highest part of the Carpathian range; it covers an area of 785 km², 22.3% in Poland (with the highest peak at 2499 m a.s.l.) and 77.7% in Slovakia (highest point, 2655 m a.s.l.). Natural vegetation shows specific patterns of altitudinal belts: lower and higher montane, subalpine and alpine. In the case of this analysis the research target was located in lower and higher montane zones (770–1600 m a.s.l.). In the upper montane belt natural Picea abies forest stands are dominating, sometimes of primeval character. The Tatra climate is also showing diversification in the vertical profile with high amounts of annual precipitation reaching 1850 mm at 1987 m a.s.l. This area is protected by bilateral national parks Polish (Tatrzański Park Narodowy—TPN) and Slovak (Tatranský národný park—TANAP). UNESCO Tatra Transboundary Biosphere Reserve was created in 1992, covering the whole Tatras and large adhering areas (Figure 1).

Picea abies is a coniferous tree species belonging to the Pinaceae family. It can reach a height of 50–60 m and 200–300 years of age (very rarely more). The range of Picea abies covers a large area of Northern Europe, from sea level till 2400 m a.s.l. in the Alps, where it grows as a krumholtz tree. It is shade- and cold-tolerant, so may persist where deciduous trees cannot withstand frost. Due to its flat root system, Picea abies do not tolerate droughts; therefore, they do not grow in central Poland, where precipitation is low. The most important habitat for this tree in the Tatras is the upper montane zone stretching at an altitude of approx. 1250–1550 m a.s.l., where it naturally forms extensive pure stands.

In the nineteenth century, most of the mixed forests in the lower montane zone of the Tatras (below 1200 m a.s.l.) were cut down for needs of mining and metallurgy, and herding caused another part of destruction. In the deforested areas, spruce was often planted as a fast-growing tree, a valuable resource. Spruce, as a pioneering species, also spontaneously enters open areas in lower zones, where higher air and soil temperatures prevail, which allow the bark beetle to limit the occurrence of spruce [7]; for this reason, increasing temperatures in the mountain environment leads to the dieback of mountain spruce forests [26].

Over the past two hundred years, the Tatra Mountains have experienced huge windthrow many times [27]. In 2004, there was heavy windthrow on the Slovak side, and a strong increase in the number of trees infested by the spruce bark beetle has been observed since 2005 [17,22,28]. At the end of 2013, heavy windthrow occurred in the Polish side, mostly in the western part of the National Park.

3. Materials and Methods

The basis of this research was the Sentinel-2 satellite images (13 spectral bands; four ten-meter, six twenty-meter, and three sixty-meter bands were used for atmospheric corrections only). Images were acquired on 3 October 2015, 30 September 2016, 2 October 2017 and 15 October 2018. The reference materials were provided by the Tatra National Park (TPN):

• WorldView-2 satellite image with a resolution of 0.5–2 m acquired on 7 September 2015;
• airborne orthophotomap (0.25 m; 2 October 2017);
• shapefiles of damaged forest stands by bark beetles (update for 2017), sanitary cuttings carried out in the park (updated for 2017), boundaries of protection zones, forest areas, windthrow and a forest map of the TPN area.

The first step was the atmospheric and topographic correction of the Sentinel-2 images; for this purpose, the Sen2Cor software (ESA) was used. Some areas were manually masked out due to clouds and their shadows. Twenty-meter bands were resampled to ten-meter resolution based on the nearest neighbor method (Figure 2).

The preparation of reference polygons based on a selection of six land cover classes: coniferous and deciduous forests, snags, sanitary cuts or overturned by windthrow, grasslands, and rocks (

Table 1. Number of reference polygons and pixels by a class for the 2015 Sentinel-2 image classification.

Class	Number of Reference Polygons	Number of Reference Pixels
Coniferous forests	39	9892
Deciduous forests	26	3114
Grasslands	60	9291
Rocks	40	9443
Snags	40	2347
Cuts or windthrow	32	5440

). In the first stage, the analyses were conducted on:

• a high-resolution WorldView-2 image;
• a vector layer updated for 2015: snags, sanitary cuts and windthrow and other forest parameters of TPN (all vector data were field verified).

Selected polygons were used to classify Sentinel-2 images from 2015. In the second stage of the work, reference polygons were created for the state for 2017:

• an actual airborne orthophotomap (acquired on 2 October 2017 with an accuracy of 0.25 m);
• an updated vector layer of snags, sanitary cuts and windthrow, which allowed us to assess changes between 2015 and 2017 (

  Table 2. Number of reference polygons and pixels between 2015 and 2017 used for change detection.

  Class	Number of Reference Polygons	Number of Reference Pixels
  No changes	40	21,301
  Dieback of trees	48	1260
  Cuts or windthrow	52	1865

  ).

It allowed us to prepare reference polygons to analyze:

• no changes between 2015 and 2017 (an undisturbed coniferous forest in 2015 remains undisturbed in 2017);
• dieback of trees—undisturbed coniferous forest in 2015, and dying in 2017;
• deforestation (sanitary cuts) or windthrow, undisturbed coniferous forest in 2015, and nonexistent in 2017 (this class may also include gaps in the stand, e.g., due to old age).

The spectral patterns acquired in 2015 and 2017, allowed to conduct analyses for 2016 and 2018 as one of the research scenarios.

3.1. Classification

Support Vector Machine (SVM) was used as the classifier based on the EnMap-Box 3.3 (QGIS 3.4) and R Studio software; all kernel functions (linear, polynomial, sigmoid and Radial Basis Function) were tested, but for the final classification, RBF function was applied. The sklearn library was applied and a tuning of SVM parameters allowed us to select the optimal parameters: gamma coefficient 0.01 was selected from the tested values: 0.001, 0.01, 0.1, 1, 10, 100 and 1000, and penalties C 1000 (from the checked: 0.001, 0.01, 0.1, 1, 10, 100, 1000). The reference polygons were sampled 10 times for training and validation in the 75/25 ratio. In the final version of the map, each pixel was assigned to the most frequently occurring class in subsequent iterations, the so-called dominant. The verification polygons were also rasterized and were used for verification and assessment of classification accuracy. The final accuracy of the classifications was measured by overall accuracy, Kappa coefficient, user accuracy (UA) and producer accuracy (PA), as well as the F1-score and error matrices.

After obtaining the final map containing 6 classes (coniferous forests, deciduous forests, grasslands, rocks, sanitary cuts and windthrow and snags) for 2015, all classes, except undisturbed coniferous forest, created a mask removing non-coniferous forest from further investigations, because the main idea was to follow the changes on the references coniferous forest patterns (initial moment of the bark beetle outbreak). On this basis, all satellite images were analyzed only on unmasked areas.

3.2. Multi-Temporal Analyses

The Sentinel-2 image from 2 October 2017, after pre-processing and masking undisturbed coniferous forest area from 2015, was compared with the image from 2015. Six methods were used to determine the changes in the state of coniferous forests from 2015–2017:

• remote-sensing indices (

  Table 3. Remote sensing indices used for forest analysis. Explanation: blue: Sentinel-2 band no. 2, green: 3, red: 4, Red Edge (RE1): 5, Red Edge (RE2): 6, REDEDGE: 7, NIR: 8, SWIR1: 11, SWIR2: 12.

  Index	Formula	Reference
  Atmospherically Resistant Vegetation Index (ARVI)	$\frac{NIR-\left(2RED-BLUE\right)}{NIR+\left(2RED-BLUE\right)}$	[29]
  Blue Normalized Difference Vegetation Index (BNDVI)	$\frac{NIR-BLUE}{NIR+BLUE}$	[30]
  Green Normalized Difference Vegetation Index (GNDVI)	$\frac{NIR-GREEN}{NIR+GREEN}$	[31]
  Enhanced Vegetation Index (EVI)	$2.5\ast \left(\frac{NIR-RED}{NIR+6RED-7.5BLUE+1}\right)$	[32]
  Enhanced Vegetation Index 2 (EVI 2)	$2.4\ast \left(\frac{NIR-RED}{NIR+RED+1}\right)$	[33]
  Green Soil Adjusted Vegetation Index (GSAVI)	$1.725\ast \left(\frac{NIR-GREEN}{NIR+GREEN+0.725}\right)$	[34]
  Normalized Difference Vegetation Index (NDVI)	$\frac{NIR-RED}{NIR+RED}$	[35]
  Normalized Green-Red Difference Index (NGRDI)	$\frac{GREEN-RED}{GREEN+RED}$	[36]
  Soil Adjusted Vegetation Index (SAVI)	$1.428\ast \left(\frac{NIR-RED}{NIR+RED+0.428}\right)$	[37]
  Simple Ratio / Short Wave Infrared (SR/SWIR)	$\frac{SWIR1}{SWIR2}$	[8]
  Green Leaf Index (GLI)	$2\ast \left(\frac{GREEN-RED-BLUE}{GREEN+RED+BLUE}\right)$	[38]
  Chlorophyll Vegetation Index (CVI)	$\frac{NIR\ast RED}{GREEN\ast GREEN}$	[39]
  Leaf Chlorophyll Index (LCI)	$\frac{NIR-REDEDGE}{NIR+RED}$	[40]
  Modified Chlorophyll Absorption in Reflectance Index (MCARI)	$(RE1-RED-0.2\ast$$\left(RE1-GREEN\right))\ast \left(\frac{RE1}{RED}\right)$	[41]
  Modified Chlorophyll Absorption in Reflectance Index 2 (MCARI 2)	$1.2\ast \left(2.5\ast \left(NIR-RED\right)-1.3\ast \left(NIR-GREEN\right)\right)$	[42,43]
  Normalized Difference Red-Edge 2 (NDRE 2)	$\frac{NIR-RE2}{NIR+RE2}$	[8]
  Transformed Chlorophyll Absorption Radio Index (TCARI)	$3\ast (RE1-RED-0.2\ast$$\left(RE1-GREEN\right)\frac{RE1}{RED})$	[42]
  Plant Pigment Ratio (PPR)	$\frac{GREEN-BLUE}{GREEN+BLUE}$	[44]
  Disease Water Stress Index (DSWI)	$\frac{NIR+GREEN}{SWIR1-RED}$	[45]
  Moisture Stress Index (MSI)	$\frac{SWIR1}{NIR}$	[46]
  Moisture Stress Index 2 (MSI 2)	$\frac{SWIR2}{NIR}$	[47]
  Normalized Difference Infrared Index (NDII)	$\frac{NIR-SWIR1}{NIR+SWIR1}$	[48]
  Normalized Difference Infrared Index 2 (NDII 2) Normalized Burn Ratio (NBR)	$\frac{NIR-SWIR2}{NIR+SWIR2}$	[49]

  );
• Sentinel-2 spectral signatures used from 2015 SVM classification to analyze class changes in 2017;
• Multivariate Alteration Detection (MAD);
• Multivariate Alteration Detection/Maximum Autocorrelation Factor (MAD/MAF);
• Iteratively Re-weighted Multivariate Alteration Detection (iMAD);
• raster stacking method of combining the 2015 and 2017 periods into one 20-band raster; 10 spectral bands acquired in 2015 and an additional 10 bands from 2017. The following spectral bands were used: 2 (Blue), 3 (Green), 4 (Red), 5 (Red Edge 1), 6 (Red Edge 2), 7 (Red Edge 3), 8 (NIR), 8A (narrow NIR), 11 (SWIR), 12 (SWIR); twenty-meter bands were resampled using the nearest neighbor method to a ten-meter pixel size.

Based on the mask of polygons (defined as healthy coniferous forests in 2015), 23 remote sensing indices were calculated (Table 3). The indices describe photosynthetically active energy accumulation, water content, biomass and general quantity and vigor of green plants [8].

All these indices were saved as separate raster images, which were combined into a single file through the layer stacking function. This resulted in a reference dataset for healthy coniferous forests. It allowed us to distinguish differences in the values for tree cuttings and dying trees. All calculations were made in the R package using the raster library [50,51,52]. In the next step, the significance of remote sensing indices in the classifier learning process was checked; for this purpose, the Recursive Feature Elimination (RFE) method was used, based on the subtraction of successive indicators from the full dataset and data ranking (from the most to the least informative).

Multivariate Alteration Detection (MAD) allowed the analysis of multispectral data of the same area from 2015–2017 [53,54]. Due to the propagation of differences between individual spectral bands in different periods, a transformation of the input bands was applied, giving more weight to those bands, which had more significant changes; in addition, a canonical correlation was used (identifying linear coefficients for each band of both images); it made it possible to create canonical weights of successive degrees, being orthogonal to the first degree up to the k-degree canonical weights. The analyses were performed in Orfeo ToolBox, generating a new image consisting of 10 bands (corresponding to the number of input Sentinel-2 bands from the 2015 and 2017 images, sixty-meter atmospheric bands were removed). In the next step, the Maximum Autocorrelation Factor (MAF) [55] method was used to rank the MAD bands from the most to the least autocorrelated. It consists of finding a linear transformation (as in the case of PCA and MAD) of the input bands to a band space in which the correlation of pixels with its surroundings is the greatest one.

Simultaneously with the MAD/MAF methods, the Iteratively Re-weighted Multivariate Alteration Detection was tested (iMAD) [56]. The IMad.py script written in Python provided by Dr. Mort Canty on GitHub was used for this. In the iMAD algorithm, the MAD transformation was performed in the first step. Then, more weights were given to those pixels that did not change significantly and less to those that showed a change. Then, the MAD was recalculated, but considering the weights given to the pixels, i.e., those pixels that showed a large change in the first step, were taken into account less for the correlation computation than those that had a small change. As a result, from two Sentintel-2 images for 2015 and 2017, one raster file with 10 spectral bands was obtained, but ranked from the highest probability of changes between images to the smallest.

In the R Studio, the layer was combined using the stack function from the raster library coniferous forests (raster stacking) from Sentinel-2 image from 2015 and 2017. In this way, a raster file with 20 bands was classified.

3.3. Map of Bark Beetle Outbreaks and Assessment of Classification Quality

As a result of the 10-fold randomized sampling of training and validation polygons, 10 maps of changes in the coniferous forests were obtained for all methods, then for each pixel from all maps, the modes were calculated, and on this basis, the maps of changes were prepared. This analysis was performed in the R program based on the e1071 library [57] and the one versus one method. These maps were generalized, i.e., eliminating less than 30 pixels for the whole Tatras. In the following step, a confusion matrix was created, presenting the number of pixels from the given reference classes and classified into individual classes on the diagonal of the matrix (reference classes were presented in columns, and the rows were classified according to the algorithm). An error matrix was calculated for each of the applied methods (remote sensing indices, Sentinel-2 spectral signatures, MAD, MAD/MAF, iMAD and raster stacking of the 20-bands file. In R Studio, a coniferous forest layer from a Sentinel-2 image from 2015 and 2017 was combined using the stack function from the raster library. Sentinel-2 satellite images from 3 October 2015 and 2 October 2017 after pre-processing and using healthy coniferous forest masks from 2015 and 2017 were merged into a 20-band raster file which was subjected to classification. Additionally, to evaluate the usefulness of the methods, the following metrics were calculated: overall accuracy of the classification (overall accuracy, OA, Aov), Kappa coefficient (κ), user accuracy (UA) and producer accuracy (PA) for each classification class and F-score.

4. Results

The result of the SVM RBF classification of Sentinel 2 images is the land-cover map for 2015 (3 October 2015; Figure 3). It allowed us to obtain a spatial distribution of coniferous forests, deciduous forests, grasslands, rocks, snags and sanitary cuts or windthrow. The coniferous forests layer served as a reference mask for the analyzes from the following years (2016–2018) and the assessment of the pace and extent of bark beetle outbreaks in the Tatra Mountains. The producer and user accuracies for individual classes is high (within the range of 90−100%;

Table 4. Error matrix for the Tatra land cover classification based on the Sentinel-2 image (3 October 2015); overall accuracy 97%, Kappa coefficient 0.96. Explanation: UA—user accuracy; PA—producer accuracy.

	Reference Data						UA (%)
	Snags	Coniferous Forest	Deciduous Forest	Grasslands	Rocks	Cuts or Windthrow
Snags	5410	1	1	29	80	474	90.2
Coniferous forests	2	24,958	52	1	97	2	99.4
Deciduous forests	0	113	6872	181	0	37	95.4
Grasslands	12	0	191	20,080	18	490	96.6
Rocks	23	0	0	2	22,054	46	99.7
Cuts or windthrow	211	2	1	597	131	11,798	92.6
PA (%)	95.6	99.5	96.6	96.1	98.5	91.8

). The error matrix presents that the places in the forest where cuts or windthrow occurred were also classified as snags or grasslands. High user accuracy for coniferous forests—99.4%—means that the created coniferous forest mask contains a high-accuracy reference layer for change detection.

Analyzing the dominant land cover forms for 2015, which becomes the reference year of observations in the analyzed period (2015−2018), it is worth noting that the undisturbed coniferous forest covers 42.27% of the area, and the snags—9.21% (

Table 5. Share of individual forms of land cover of the Tatra Mountains (both national Parks: TPN, TANAP), source: Sentinel-2 3 October 2015.

Class	Area (km2)	Percentage (%)
Snags	48.65	9.21
Coniferous forests (c forest)	223.39	42.27
Deciduous forests	66.54	12.59
Grasslands	120.93	22.88
Rocks	1.89	0.36
Cuts or windthrow	67.08	12.69

).

Analyzing changes that took place in coniferous forests between 2015 and 2017, it should be noted that all used methods allowed us to obtain very high and comparable results (

Table 6. Comparison of the accuracy of different methods of detecting changes between 2015 and 2017.

Method	OA (%)	Kappa Coefficient	Snags F1-Score (%)	Coniferous Forests F1-Score (%)	Cuts or Windthrow F1-Score (%)
Raster stacking	99.1	0.97	96.9	99.8	95.2
Remote sensing indices	99.0	0.97	96.4	99.9	93.8
iMAD	99.0	0.97	96.3	99.9	93.6
MAD	98.5	0.96	95.2	99.4	93.8
MAD/MAF	98.5	0.96	95.1	99.4	93.7
Spectral features obtained from Sentinel-2 2015 and SVM classified on Sentinel-2 2017	93.8	0.77	70.4	99.8	58.8

); Figure 4 presents differences of algorithm outcomes, where large-scale changes can be easily visible, which allowed us to show the influence of algorithms on the identification of the analyzed phenomena. Nevertheless, the visual interpretation of the results confirms that the outcomes based on the raster stacking were characterized by a large number of individual pixels (salt and pepper effect), while this unfavorable phenomenon was not observed in the case of remote-sensing indices; hence, this method was used to present the obtained results.

The most informative remote-sensing indices were selected by the Recursive Feature Elimination method. The informational character of individual indices depended on the analyzed object (Figure 5), while the use of several indices allowed us to differentiate analyzed classes with an accuracy level over 90% (SAVI, which almost unambiguously allows us to distinguish undisturbed forests from cutting and dieback of trees (F1-score is 99.60%)). SAVI is based on red and infrared spectral bands (with 10 m spatial resolution). Adding additionally MCARI2, which was the second most informative in the ranking, which also contains information about the green band, significantly improved the identification of snags (F1-score was 93.19%) from the cuts and windthrow area class (F1 = 87.49%). The analyses showed the usefulness of the set of 15 remote sensing indices: SAVI, MCARI 2, NDII 2 (NBR), LCI, PPR, GSAVI, BNDVI, NDII, MSI, NGRDI, SR/SWIR, NDRE 2, NDVI, EVI and CVI (Figure 5), adding the information contained in the following indices: MCARI, TCARI, EVI 2, GNDVI, GLI, MSI 2, ARVI and DSWI did not improve the obtained results. These indicators allowed us to distinguish the state of the coniferous forest with a very high accuracy (Table 6), as well as to precisely identify the damage in 2017 (Figure 6, Table 6).

In 2018, the total forest degradation compared to 2015 exceeded around 30% for both Parks (Figure 7 and Figure 8), while the share of cuts/windthrow areas in TANAP is higher by two–three percentage points (

Table 7. Snags and cuts/windthrow class statistics of the Polish (TPN) and Slovak (TANAP) Tatra National Parks. Data presented in hectares. Acquired on the SVM classification based on the 15 remote sensing indices.

Class	TPN (ha)	TANAP (ha)	Sum (ha)
Undisturbed coniferous forest	6014.46	13,546.43	19,560.89
Snags	726.47	822.00	1548.47
Cuts or windthrow	236.18	951.64	1187.82

).

5. Discussion

The achieved results allow us to assess the usefulness of Sentinel-2 satellite images and SVM-supervised classification for spruce bark beetle outbreak with high accuracies; in our case the worst results oscillated around 94% (overall accuracy). Different classifications of dying trees as a result of cambiophagous insects were analyzed in many parts of the world (

Table 8. Classification of the gradation of cambiophagous insect’s assessment of results. Explanation: OA—Overall Accuracy; PA —Producer Accuracy; UA—User Accuracy.

Reference	Sensor	Research Target	Classification Method	OA (%)	Snags PA (%)	Snags UA (%)
[23]	Sentinel-2	Ips tyhographus, snags	SVM	97	91.71	90.72
[21]	Landsat TM, ETM+	Dendroctonus ponderosae, snags	ML	91	88.3	93.7
[24]	WorldView-2	Ips tyhographus, green phase of the bark beetle attack and snags	RF	76.2	100	100
[22]	Landsat TM, ETM+	Ips typographus, snags.	MSI, VCI, NDMI–threshold	94
[59]	Landsat TM, ETM+, OLI	Tomicus yunnanensis and Tomicus minor	MSI	86.38
[62]	WorldView-2	Ips tyhographus, snags, mountain pine.	SVM	86.6	86.75	100
[58]	Sentinel-2	Dendroctonus valens: three classes: 0%, 0–15%, 15–50% of dead trees	SVM	59.5	-	-
[58]	Sentinel-2	Dendroctonus valens: two classes: 0%, 15–50%	SVM	81	-	-
[58]	GaoFen-2	as above, additionally: alive, attacked and dead tree	SVM	77.7	-	-
the current study	Sentinel-2	Ips tyhographus, snags, cuts and windthrow	SVM	97.2	95.6	90.2

), e.g., China [58,59], Europe [13,24,60,61], North America [21]; the authors based on different spatial and spectral resolution acquiring satellite or airborne images. The most similar analysis to this work, both in terms of geography and theused sensor (Sentinel-2) and methodology (SVM + MAD), is the work of G. Mikusiński and his team [23], who analyzed spruce bark beetle outbreaks and clearing cuts in the Białowieża Forest (NE Poland). A single-image Sentinel-2 classification for 2015 and a classification for 2018 change detection were performed; an overall accuracy oscillated at 97% (in our work—97.19%), and for the class of snags, the producer accuracy achieved 91.71%, and user accuracy—90.72% (we received 95.6% and 90.2%, respectively, which is a similar result to Mikusiński’s work [23]. In the multitemporal analysis (2015–2018), Mikusiński’s team achieved 97% (overall accuracy) and 92.2% (F1-score) for snags and 95.1% for the slice class [23]. In this paper, it is, respectively, 99.1% of overall accuracy, 95.7% (F1-score) for snags and 93.6% of F1-score for cuts/windthrow. The higher overall accuracy of this work resulted from correctly classified healthy forest, while the lower F1 score for cuts/clearing class, this is due to the fact that the Białowieża Forest is located in the lowland, where heavy equipment (e.g., harvester) can be used, which allows us to have homogenous large cuts/clearing areas. Moreover, the use of heavy equipment could have degraded the soil more, which could also have contributed to larger spectral differences between classes.

The team of A. Meddens [21] analyzed four classes: healthy forest, dead pines (over 70% in a 30-meter pixel), non-forest vegetation and masked areas (clouds, cloud shadows, topographic shadows) using the Maximum Likelihood, achieving 91% overall accuracy and 88.3% producer accuracy for snags and 93.7% user accuracy. The results are comparable to our work or slightly better; it may be a result of a lack of clear cuts, which were sometimes mixed up with the snags in our case.

M. Immitzer and C. Atzberger [24] applied the Random Forest and the WorldView-2 satellite images to classify individual trees during the outbreak of the spruce bark beetle in Austria. The authors focused on undisturbed forest areas, attacked green trees and dead trees. For the dead tree class, they achieved 100% for producer and user accuracies. This shows that the higher spatial resolution allows distinguishing snags and healthy trees. P. Rębiś [62] analyzed the bark beetle outbreaks in the Tatra Mountains; he used WorldView-2 images and Support Vector Machines to distinguish undisturbed forest, dwarf mountain pine (Pinus mugo), bark beetle nests and non-forest vegetation. The results oscillated at around 86.94% overall accuracy, 86.75% producer accuracy for snags and 100% user accuracy. The lower overall accuracy than in this study is due to the mixing of spectrally similar classes of dwarf mountain pine and healthy coniferous forests. This means that there were both green and dead trees in this class, which caused the mixing of the bark beetle nest classes with the healthy forest.

Z. Zhan’s team [58] analyzed a possibility of detecting dead Chinese pines among undisturbed forests in various shares of the species: 0% (undisturbed forest), 0–15% and 15–50%. For this purpose, authors used images of Sentinel-2 and Gaofen-2 (1–4 m pixel size); based on the CART (Classification and Regression Tree, which is a method of decision trees for the classification or regression of data), Random Forest and SVM methods, the authors obtained 59.5% overall accuracy for the three classes on Sentinel-2, but eliminating the class of 0–15%, overall accuracy increased to 81%. This shows moderate effectiveness in investigating snags in mixed pixels. The use of high-resolution Gaofen-2 images resulted in a 77.7% accuracy in distinguishing the three states of the attacked trees.

A similar level of multitemporal Sentinel-2 images and SVM RBF classification accuracies was observed for areas of Catalonia and the European Lowlands (near Warsaw, Poland), where coniferous forests were classified as the best (the median of the F1 score oscillated around 95%), deciduous (82–97%) and mixed forests (82–92%) [63]. Moreover, the values of the vegetation indices of stands in good condition were similar to those obtained for alpine grasslands [64], which were classified on a similar level of overall accuracy, reaching 84% for 22 vegetation communities based on airborne APEX hyperspectral images [65].

The second part of our study was focused on change detection between the reference year (2015) and the following years (2016–2018). Our results indicated very high accuracy between observed changes and reference patterns (oscillating between 93–99% of overall accuracy). Similar analyses were performed by the mentioned team of G. Mikusiński with the MAD/MAF algorithm [23]. Similar work was led by M. Havašová [22] on the bark beetle outbreak in the Tatra Mountains. Using the Landsat-5 and Landsat-7 satellite images and remote sensing indices, the threshold for the decline in indices between years, which would mean the occurrence of snags, was determined. For the years in which the gradation was the strongest, i.e., 2007 and 2011, 94% of the total accuracy was achieved, while this method fails in the remaining years. The aforementioned work by A. Meddens [21] investigated the deviation of remote sensing indices over the years. The highest overall accuracy was obtained for Tasseled Cap indices—91%. Landsat-5 and Landsat-7 images were analyzed, which presented a significant improvement in the overall accuracy of the multitemporal analyses based on the Sentinel-2 images due to higher spatial resolution (99.1% achieved).

Similar results on the same Landsat sensors were achieved by the team of L. Yu [59] in their study of Pinus yunnanensis pine and its defoliation by the lesser pineapple and the Tomicus yunannensis pineapple. Using the MSI remote-sensing index and grading based on the thresholding of the inter-year difference, they achieved 86.38% overall accuracy.

6. Conclusions

Among protective measures applied in both Tatra National Parks, the cutting and removal of infected spruces is quite common. Unfortunately, it is not easy to find infected trees from the ground survey, so it is usually performed when trees are dead, long after a new generation of beetles leave the tree trunk. This is totally ineffective. If it would be possible to distinguish infected trees in the early stage based on remote-sensing methods, protective measures could have been implemented sooner, and be more effective. So far, satellite images have not been implemented into protective measures in Tatra National Parks.

Sentinel-2 images, due to their spectral, radiometric and spatial resolutions, are a helpful source of information for the monitoring of large areas, including national parks, in assessing the course of the outbreak of cambiophagous insects. Most analyses are based on the spectral properties of research targets reflected in high resolution bands (10–20 m), which allows the identification of single or clumps of trees. Based on field polygons, it is possible to identify the first stages of bark beetle attacks (green attacks). Presented in the paper, accuracies (overall user and producer and F1-score) above 90% were obtained for snags and cuts/clearing patterns. Sentinel-2 satellite data is available every five days, so it is possible to monitor even mountain forests, which are clouded, but a high revisit time allows the masking of unexpected phenomena.

The conducted analyses identified optimal parameters of the SVM classifier (RBF, C = 1000, gamma = −0.01) for conifers, mixed forests, deadwood, clearings/windbreak and non-forest vegetation. Sentinel-2 images and the proposed methods allowed us to obtain very good results (overall accuracy 97%, and Kappa coefficient 0.96), which confirms the legitimacy of the proposed solutions; at the same time, the error matrix draws attention to misleading objects during classification, especially cuttings or fallen trees that at times had a spectral reflection very similar to insect-attacked trees.

It was noted that the most informative indices are SAVI, MCARI2, NDII 2 (NBR), LCI, GSAVI, BNDVI, NDII, MSI, NGRDI, SR/SWIR, NDVI, EVI and CVI; a further increase in the number of indices does not increase the accuracy of classified forests, cuttings and bark beetle outbreaks.

All analyzed transformations (MAD, iMAD, MAD/MAF) offer very high levels of accuracies, but the best results were obtained for iMAD (0.5–1% higher the others, but still lower than the raster stacking method). Subsequent change detection algorithms (MAD, iMAD, MAD/MAF, remote-sensing indices and raster stacking) did not bring significant changes in the results. The differences in the overall accuracy as well as the producer and user accuracies of each class were up by two percentage points, because all these algorithms are algebraic operations on spectral bands that are transformed by non-parametric algorithms, e.g., Support Vector Machines. A real qualitative change could already be brought by an improvement in the spectral or spatial resolution, as shown by the discussion of the results or the non-pixel but object-oriented approach, in which not only information from a given pixel is analyzed, but also the context in which it is located.

An interesting issue not addressed in this study was also the attempt to classify newly attacked spruce trees (green attacks). The high classification results in this work, as well as the promising results of H. Abdullah’s team [8], may indicate this direction of work on the course of the gradation as the most forward-looking.

The analysis showed that between 2015 and 2018, in the entire analyzed area of the Tatra Mountains, 30% of coniferous forests were damaged, of which approximately 5% were cuts and windthrow areas and 25% snags, probably mainly due to the bark beetle.