Estimating VAIA Windstorm Damaged Forest Area in Italy Using Time Series Sentinel-2 Imagery and Continuous Change Detection Algorithms

Abstract

Mapping forest disturbances is an essential component of forest monitoring systems both to support local decisions and for international reporting. Between the 28 and 29 October 2018, the VAIA storm hit the Northeast regions of Italy with wind gusts exceeding 200 km h⁻¹. The forests in these regions have been seriously damaged. Over 490 Municipalities in six administrative Regions in Northern Italy registered forest damages caused by VAIA, that destroyed or intensely damaged forest stands spread over an area of 67,000 km². The present work tested the use of two continuous change detection algorithms, i.e., the Bayesian estimator of abrupt change, seasonal change, and trend (BEAST) and the continuous change detection and classification (CCDC) to map and estimate forest windstorm damage area using a normalized burned ration (NBR) time series calculated on three years Sentinel-2 (S2) images collection (i.e., January 2017–October 2019). We analyzed the accuracy of the maps and the damaged forest area using a probability-based stratified estimation within 12 months after the storm with an independent validation dataset. The results showed that close to the storm (i.e., 1 to 6 months November 2018–March 2019) it is not possible to obtain accurate results independently of the algorithm used, while accurate results were observed between 7 and 12 months from the storm (i.e., May 2019–October 2019) in terms of Standard Error (SE), percentage SE (SE%), overall accuracy (OA), producer accuracy (PA), user accuracy (UA), and g_mean for both BEAST and CCDC (SE < 3725.3 ha, SE% < 9.69, OA > 89.7, PA and UA > 0.87, g_mean > 0.83).

1. Introduction

The effects of climate change and climate variability on forest ecosystems are evident around the world and further impacts are unavoidable, at least in the short and medium period [1]. Under global warming, in this century the frequency and severity of wind gust events is increasing [2,3], especially in Europe where since 1950 over 130 autumn and winter windstorms have caused notable damages to European forests [4,5,6]. In fact, autumn and winter windstorms are responsible for more than 50% of reported primary damages to European forests (in terms of growing stock volume), including all biotic and abiotic causes [4,6]. At the European level, most of the autumn/winter windstorm damage is reported in coniferous forests probably because their leaf-on canopies offer more resistance against wind compared to deciduous canopies which have shed their leaves in the winter [4,7,8].

Understanding the effects of wind gusts on forests is of great interest for environmental and economic issues [4,9,10] and for the assessment of greenhouse gas inventories [11,12]. In fact, after a windstorm, it is necessary to rapidly obtain estimates about the damaged forest area, to support forest operations, and future management strategies [13,14,15], and to report forest indicators in the context of international reporting obligations [2] such as FAO [16], the Ministerial conference on the protection of Forest Europe [17], and IPCC [11].

Remote sensing (RS) is considered the most efficient way to map forest disturbances over large area. Many RS data were tested to map and estimate forest windstorm disturbances as for example Airborne Laser Scanner data (ALS) [8], Sentinel-1 SAR data [13,18], high and medium spatial resolution multispectral satellite imagery [15,18,19,20]. However, some of these approaches present relevant limitations. For example, the one proposed by Chirici et al., (2018) [8] requires the acquisition of post-event ALS data, which is expensive for large area application (i.e., 0.8 €/ha). The one using SAR data proposed by Rüetschi et al. (2019) [13], is not able to determine the exact boundary of damages that is mandatory to estimate the damaged area and cannot be applied for winter storm and/or in absence of a high-resolution digital terrain model (i.e., 2 m resolution). The approaches proposed by Kislov & Korznikov (2020) [19] and Hamdi et al. (2019) [20] using high resolution imagery and deep learning seems to be promising in mapping the boundaries of windstorm disturbances; however, the cost of high resolution images over large areas is high.

Within the vast availability of different remotely sensed data, the freely available multispectral optical satellite images are still considered the most feasible and cost-effective option to derive information for extensive forest areas [21] since such data are freely available online at medium and fine resolution, e.g., Landsat (spatial resolution: 30 m, revised time: 16 days) and Sentinel-2 (S2) (spatial resolution: up to 10 m, revised time: 5 days). Dalponte (2020) [15] reported that the efficiency in detecting windstorm damages using optical RS data depends mainly on three main aspects: temporal resolution, spectral resolution, and spatial resolution.

Usually to map forest disturbances common change detection (CD) methods with optical data consist of analyzing the pre- and post-disturbance images are used [15,19]. For example, Dalponte et al. (2020) [15] in a study area of 732 km² mapped forest windthrow caused by VAIA with a change vector analysis (CVA) bi-temporal approach between two dates as a vector within the features space, using two high resolution planet scope (PS) images and two medium resolution S2 images. Dalponte et al. (2020) [15] showed that the most accurate map of windstorm damages areas can be derived after 8 months from the storm event comparing two images acquired in summer time, one pre- and one post-event (S2 Kappa accuracy (KA) = 0.74, PS KA = 0.70), while less accurate results can be obtained using images acquired closer to the storm event e.g., after 14 days PS and after 45 days S2 (KA ≤ 0.29). Vaglio Laurin et al. (2020) [18], in an area of 20,000 km² damaged by VAIA windstorm found that S2 is more efficient to detect forest damages than SAR Sentinel-1 (S1). In fact, Vaglio Laurin et al. (2020) [18] with a bi-temporal change detection approch obtained an overall accuracy of 86% for S2 using an image acquired after 7 months from the storm and 68% for S1 using data acquired after 15–20 days after the storm. Also, Olmo et al. (2021) [22] using a bi-temporal approach between vegetation indices (vegetation index differences) calculated on S2 bands produced accurate results in detecting VAIA forest windthrow in Friuli Venezia Gilia. Olmo et al. [22] showed with a multitemporal analysis based on vegetation indices that, starting from April 2019, it was possible to observe a deviation of reflectance values for damaged and undamaged areas comparing to the years before the storm. Piragnolo et al. [23] showed that the multitemporal analysis of vegetation indices calculated on the S2 imagery are useful to predict severity classes of damaged areas using aggregational statistics of VIs as input to random forest machine learning algorithm.

Several authors have demonstrated that an analysis based on the multitemporal spectral signatures, derived by satellite time series (TS) imagery, can be more efficient [19,24,25,26,27,28,29] than a bi-temporal approch. In fact, in recent years, new methods based on the multitemporal spectral signatures [19,24,25,26,27,28,29] were developed to map forest disturbances and forest interannual changes using multispectral optical satellite TS (i.e., Landsat, MODIS, and in the last year S2 data) [24,27,28,29,30,31,32,33,34] thanks also, to the development of cloud computing platform such as Google Earth Engine (GEE). Cloud computing platform, in fact, gives access to a complete catalog of remotely sensing data, with the capability to process efficiently large multitemporal set of remotely sensed data [35].

Within the vast array of algorithms developed to map land use changes using TS optical satellite imagery, those that are able to detect interannual changes using trajectory analysis appears to be adequate to detect interannual changes of forest area [36,37]. Some of the most promising CD approaches are: the continuous change detection and classification (CCDC) [38]; the exponentially weighted moving average change detection (EWMACD) [39]; or the Bayesian estimator of abrupt change, seasonal change, and trend (BEAST) [40] methods. These algorithms were successfully applied to monitor decline processes [41], changes in vegetation growing [42], seasonal forest changes [43], and spatiotemporal forest dynamics [44]. In addition, these algorithms are able to split the time series into three adaptative components (i.e., trend, seasonal, and remainder) [42,44,45] for each investigated year. Many of these studies use dense Landsat or MODIS TS as input data [40,46], while, the use of S2 and CCDC algorithm to map forest areas damaged by windstorms is limited [37]. Many recent studies on VAIA windstorm underlined that S2 is adequate to detect windstorm since S2 program offers innovative features for forest remote sensing by combining high spatial resolution (i.e., 13 bands, from 0.443 to 2.190 µm with the visible (i.e., R, G, B,) and the near infrared bands at 10-m spatial resolution and four red-edge bands at 20-m spatial resolution), wide coverage and a quick revisit time (i.e., every 5 days after the launch of Sentinel-2B satellite in 2017) [15,18,22,23]. However, they considered only a small portion of the area hit by VAIA windstorm and no-ones provided an estimation of damaged area.

The aim of this work is to test the accuracy of continuous change detection algorithms in mapping forest damaged forest areas by VAIA windstorm over approximately 67,000 km², of which 144,651.1 ha are forested following the FAO forest definition. We tested two of the most used TS Continuous Change Detection algorithms (BEAST and CCDC) to detect VAIA windstorm damage forest areas using three years continually S2 normalized burn ratio (NBR) time series (i.e., January 2017–October 2019). The intent was to analyze how the accuracies of the maps and area estimation change within 12 months from the storm. The accuracies of the maps derived by the two algorithms were evaluate using a ground truth point dataset. Moreover, since area estimates is requested in the context of international reporting, such as greenhouse gas inventories, we calculated with probability-based stratified estimator the total damaged area and the standard error of the damage area estimates. To the best of our knowledge this contribution presents the first tentative of mapping the total forest area damaged by the VAIA storm.

2. Materials

2.1. Study Area and VAIA Windstorm

The study area coincides with the Italian area hits by the VAIA windstorm as reported by Chirici et al. (2018) [47] in the following Regions of Italy: Trentino Alto Adige, Veneto, Friuli Venezia Giulia, Lombardia, Piemonte, and Valle d’Aosta (Figure 1), for a total of 67,000 km² covered by forest for 144,651.1 ha following the FAO forest definition [48,49].

The VAIA storm hit the area in the night between 28 and 29 October 2018 with winds up to 200 km/h [45]. Local authorities on the basis of field surveys reported that most of the damages were in coniferous forests dominated by Spruce (Picea abies L.), Silver Fir (Abies alba Mill.), and Pinus (Pinus spp.), however some damages were reported also in broadleaves forests such as in beech stands (Fagus sylvatica L.) [48]. Preliminary field recognition and manual photointerpretation of satellite data showed that a total of 42.525 ha of forests were damaged by the storm [47,49]. The minimum mapping units used to delineate damaged area was 2000 m² and the areas were considered damaged if at least 50% trees fallen.

2.2. Sentinel-2 Time Series Preprocessing

We used surface reflectance Bottom-Of-Atmosphere S2 images available in GEE with cloud cover lower than 10% acquired between the 1 January 2017 and the 31 December 2019 resulting in a total of 1360 images of which 707 were acquired before VAIA storm and 653 after the storm (Figure 2). From Figure 2 it is possible to note that most of the images were acquired in late spring and summer (June to September) and from Figure 3 it is possible to see the monthly cloud composite of S2 imagery derived by the use of image collection after VAIA windstorm.

The pixels in the images of the collection covered by clouds and cloud shadows were masked using the algorithm Fmask implemented in GEE [50]. Moreover, to reduce the area of the analysis just to forest areas, according to the FAO definition, we excluded the pixels of each images outside the forests using the high resolution forest masks [49].

For each forest pixel two indices (i.e., the NBR and the normalized difference snow index (NDSI)) were calculated. NBR is calculated as the ratio of NIR (S2-Band 8) and the SWIR (S2-Band 12) bands

$$NBR= \left(NIR-SWIR\right)/\left(NIR+SWIR\right)$$

(1)

and the NDSI is calculated as a ratio between the Green (S2-Band 3) and the SWIR (S2-Band 11)

$$NDSI= \left(GREEN-SWIR\right)/\left(GREEN+SWIR\right)$$

(2)

The two indices were calculated because the NBR was used to detect damaged forest areas, and the NDSI to obtain information on snow cover to reduce noises [51]. In details, the image collection was used to create the S2 time series data. To reduce spikes due to snow covers in the NBR TS we used a despike approach [52] in correspondence of high values of NDSI. In details, when NBR spikes were observed in correspondence of high values of NDSI (i.e., NDSI > 0.1), we modeled the NBR value using a linear regression model with the two closer NBR values without snow (i.e., before and after). In fact, from Figure 3 it is possible to observe that NBR anomaly spikes (i.e., light blue line) are present when NDSI is higher than 0.1 (i.e., black points).

2.3. Training Dataset

To study the S2 time series and to calibrate the classification algorithms we created a training dataset of forest polygons covering both damaged and undamaged areas. Using a grid of 30 km × 30 km we selected a total of 100 polygons damaged by the VAIA storm from the European dataset of windstorm recently published by Forzieri et al. (2020) [4]. In details, we extracted at least one damaged polygon for each cell of the grid. Then in the same cell of the grid we manually photointerpreted a total of 100 undamaged forest polygons using high resolution satellite Planet scope images [53]. Overall, 45 ha of damaged (i.e., 1125 S2 pixels) and 47 ha (i.e., 1175 pixels) of undamaged areas were acquired for the training dataset.

2.4. Validation Dataset

To assess the accuracy of the classification methods we created a validation dataset. We started from an existing point sampling created for the Italian Inventory of Land Use (IUTI—Inventario dell’Uso delle Terre d’Italia), that is a permanent inventory done by ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale) [54]. IUTI is a permanent monitoring system to estimate the extent of six land use categories identified in the IPCC Good Practices Guidance for Land Use, Land Use Change and Forestry: 1—forest land; 2—cropland; 3—grassland; 4—wetland; 5—settlements; 6—other lands. IUTI is based on a systematic unaligned sampling design with a 0.5 km × 0.5 km grid cells. A complete description on IUTI can be found in Vangi et al. (2021) [55].

From the 128,548 IUTI points falling in forests in the study area we extracted a subsample of 700 points on the basis of a stratified random sampling. Therefore, for each administrative Province we extracted a different number of points calculated on the basis of the hectares of damage reported by local authorities and on the basis of a squared grid of 2 km × 2 km [47] with a minimum of 10 points per Province. The 700 points of the validation set were manually photointerpreted using high resolution Planet scope imagery [53] and classified into damaged or undamaged by the VAIA storm. To classify the points in damaged and undamaged forests we used a square grid of 2000 m² around the points. The point was classified as damaged if at least 50% of the square was covered by fallen trees.

3. Methods

We tested two different continuous change detection algorithms to map forest areas damaged by VAIA using the one year’s despike S2 NBR time series described in Section 2.2. Both the procedures were implemented in GEE using the function of the two algorithms available on GEE platform.

3.1. Mapping Vaia Damages

3.1.1. Breaks for Additive Seasonal and Trend Iterative Algorithm

The despiked NBR time series described in Section 2.2 was used as input data to decompose it into three components: trend, seasonal, and remainder components using the BEAST iterative algorithm [42]. The decomposition model was set as:

$$Y_t=T_t+S_t+e_t\hspace{1em}\left(t=1,\dots ,d\right)$$

(3)

where $Y_t$ was the observed data at time t, $T_t, S_t, \mathrm{and} e_t$, represent the trend, seasonal, and remainder component, and $d$ represent the number of observations. The $T_t$ was fitted using a piecewise linear model with specific slopes ${\alpha }_i$ and ${\beta }_i$ intercepts on different s + 1 segments as:

$$T_t={\beta }_i+{\alpha }_{i}t({\tau }_{i-1}^{*}<t \le {\tau }_i^{*})$$

(4)

where $i=1, \dots , s$ and m represent the number of breakpoints in the trend component. The seasonal component $S_t$ can be fitted as the piecewise harmonic models on different p + 1 segments as:

$$S_t={\sum }_{k=1}^k\left[{\gamma }_{j;k}\sin \left(\frac{2\pi kt}{v}\right)+{\theta }_{j;k}\cos \left(\frac{2\pi kt}{v}\right)\right]({\tau }_{i-1}^{\#}<t \le {\tau }_i^{\#})$$

(5)

where $j=1, \dots , p$ and p is the number of breakpoints in the seasonal component. Variable k is the order of harmonic function, which was set as 3 to accurately detect forest changes (i.e., double and triple-changes within a year).

The frequency v was set to 73 for annual daily observations for the NBR time series. The iteration process is initialized by estimating the seasonal component $S_t$ using the seasonal-trend decomposition (STL) algorithm, with the residual-based moving sum (MOSUM) [56] used to test whether breakpoints τ* 1,..., τ* s occur in the adjusted trend component $Y_t-\widehat{S_t}$. We used the minimizing Bayesian information criterion (BIC) and the residual sum of squares (RSS) to determine the number of breaking points in the trend component, and then we used the regression on the basis of M-estimation to estimate the coefficient and intercept of each of the segment, using to estimate the trend of Equation (5). Using a similar MOSUM test we identified the breakpoints ${\tau }_1^{\#},\dots , \le {\tau }_p^{\#}$, in the de-trend seasonal component $Y_t-\widehat{T_t}$, using as before BIC and RSS. Coefficient ${\gamma }_{j;k}$ and ${\theta }_{j;k}$ of the seasonal component for each of the segment are estimated with M-estimation, and then the seasonal estimate $\widehat{S_t}$ is computed based on Equation (5). The process is iterative and was performed until the number and position of breaking points remained stable. That allow to detect the breakpoints and estimate the parameters of the Equations (4) and (5).

Using the decomposition time series iteratively, the seasonality component before VAIA windstorm (i.e., from 1 January 2018 and 27 October 2018) was extracted. Then after the windstorm (i.e., from 31 October 2018 to 31 December 2019) we checked if damaged areas had a breaking point in correspondence of the remainder component (i.e., time where a series have a deviation from seasonality). The BEAST parameters were set using the training dataset and then BEAST was applied to each pixels and its own time series to detect, if present, a breakpoint of the seasonal component after the storm using the equation describe before. Therefore, applying the BEAST every month after the storm, we obtained a classification in all forested area as ‘damaged’ or ‘undamaged’ for each pixel of the study area. At the end we obtained for each month after the storm maps of damaged and undamaged area that were used to assess the accuracy of forest windstorm detection (Section 3.2) and to assess the damaged forest areas using the probability-based stratified estimators (Section 3.3). We refer to Cai & Liu (2015), Zhao et al. (2019), and Mardian et al. (2021) [40,46,57] for a complete description of the BEAST algorithm.

3.1.2. Continuous Change Detection and Classification Algorithm

As reported by Zhu (2014) [39] the CCDC algorithm assumes that noise is ephemeral and land cover change is persistent and uses all available satellite images observations at each pixel to simultaneously map land cover and land cover changes modeling temporal-spectral features including seasonality, trends, and spectral variability. We applied the CCDC implemented in GEE using the NBR despike TS as done for BEST. The CCDC time series model is able to detect three types of changes—i.e., interannual changes, gradual interannual changes, and abrupt changes. The changes that occur for pixel values over time were modeled with NBR using the ordinary least squares (OLS) method. The difference between the predicted and modeled pixel value and the true pixel value is calculated. When the difference between the values is three times greater than the root mean square error (RMSE), that pixel is flagged as a possible forest area damaged by the windstorm. Potential windstorm damaged areas were then assessed for true changes based on the number of consecutive observations. If a pixel’s value is markedly different from the model results only once, this is likely an outlier. If the pixel’s value is markedly different from the model results for a given number of consecutive observations, the algorithm considers that pixel as a real change. In our case, the minimum consecutive anomaly observations parameter was set to 2. Therefore, the pixel was considered a possible forest damage when the difference between observed and predicted images exceeds a threshold two consecutive times. Then to classify land changes the CCDC algorithm uses the coefficients of time series models as the inputs for land cover classification. In our case the coefficients of time series models were calculated on the basis of the training sample (polygons describe in Section 2.3.). At the end the coefficients from the time series models and the RMSE from model estimation were used as input to the random forest classifier (RFC). For more details on the use of CCDC algorithm we refer to Zhu et al. (2014) and Arévalo et al. (2020) [39,58]. Applying the CCDC we obtained every month after VAIA a forest damage map that is used to assess the accuracy of forest windstorm detection (Section 3.2) and to estimate the damaged forest area using the probability-based stratified estimators (Section 3.3).

3.2. Accuracy Assesment of Forest Windstorm Damage Maps

The accuracy of the pixel-level maps obtained for each month after the storm (November 2018–October 2019) by the two algorithms BEAST and CCDC was estimated against the validation dataset based on a 2 × 2 confusion matrix (

Table 1. Confusion matrix and estimators.

Map Class	Reference Class		Total	Stratum Weight ${\mathrm{w}}_1$	${\hat{\mathrm{p}}}_1=\frac{{\mathrm{n}}_{\mathrm{j}1}}{{\mathrm{n}}_{\mathrm{j}•}}$	${\mathrm{w}}_{\mathrm{j}}·{\hat{\mathrm{p}}}_{\mathrm{j}}$	$\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_{\mathrm{j}}\right)$	${\mathrm{w}}_{\mathrm{j}}^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_{\mathrm{j}}\right)$
	Damaged	Undamaged
Damaged	TP	FP	n1• = TP + FP	${\mathrm{w}}_1$	${\hat{\mathrm{p}}}_1=\frac{\mathrm{TP}}{{\mathrm{n}}_{1•}}$	${\mathrm{w}}_1·{\hat{\mathrm{p}}}_1$	$\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_1\right)=\frac{{\hat{\mathrm{p}}}_1·\left(1-{\hat{\mathrm{p}}}_1\right)}{{\mathrm{n}}_{1•}}$	${\mathrm{w}}_1^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_1\right)$
Undamaged	FN	TN	n2 = FN + FP	${\mathrm{w}}_2$	${\hat{\mathrm{p}}}_3=\frac{\mathrm{FN}}{{\mathrm{n}}_{2•}}$	${\mathrm{w}}_2·{\hat{\mathrm{p}}}_2$	$\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_2\right)=\frac{{\hat{\mathrm{p}}}_2·\left(1-{\hat{\mathrm{p}}}_2\right)}{{\mathrm{n}}_{2•}}$	${\mathrm{w}}_2^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_2\right)$

). We used the same approach of Giannetti et al. [36] and for each confusion matrix, we calculated true positive (TP), true negative (TN), false positive (FP), and false negative (FN) values. The TP, TN, FP, and FN were then used to calculate different accuracy indices such as the overall accuracy (OA), the producer’s accuracy (PA), the user’s accuracy (UA), and g_mean.

$$\mathrm{OA}=\frac{\left(\mathrm{TP}+\mathrm{TN}\right)}{\left(\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}\right)}$$

(6)

$$\mathrm{PA}=\frac{\mathrm{TP}}{\left(\mathrm{TP}+\mathrm{FN}\right)}$$

(7)

$$\mathrm{UA}=\frac{\mathrm{TP}}{\left(\mathrm{TP}+\mathrm{FP}\right)}$$

(8)

$${\mathrm{g}}_{\mathrm{mean}}=\frac{\mathrm{TP}}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)*\left(\mathrm{TP}+\mathrm{FN}\right)}}$$

(9)

We used g_mean index because as reported in Giannetti et al. (2020) a small g_mean is an indication of poor performance in the classification of the positive cases (TP = “damaged forest”), even if the negative cases (TN = “not damage forest”) are correctly classified.

To compare the results obtained by the two algorithms we used the Cohens’s Kappa (K) for binary classification [59] to measure the agreement between the maps of each month obtained by the two algorithms. To calculate the K a 2 × 2 confusion matrix between BEAST map and CCDC map for each month was used. K was calculated as:

$$\mathrm{K}=\frac{P_o-P_e}{1-P_e}$$

(10)

where $P_o$ is the proportion of observed agreement and is calculated as:

$$P_o=\left(\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}\right)$$

(11)

and $P_e$ is:

$$P_e=\frac{\left(\left(\mathrm{TP}+\mathrm{FN}\right)*\left(\mathrm{TP}+\mathrm{FP}\right)+\left(\mathrm{FP}+\mathrm{TN}\right)*\left(\mathrm{FN}+\mathrm{TN}\right)\right)}{{\left(\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}\right)}^2}$$

(12)

The K can take values from −1 to +1. Negative values indicate that the observed agreement is worse than what would be expected with a random classification and K values below 0.60 indicate a significant level of disagreement.

3.3. Probability-Based Stratified Estimators

Because the estimation of confidence intervals of the estimated area for a map class among two forest classes (i.e., undamaged forests and damaged forests) is reduced to a two-class problem, as reported by McRoberts et al. (2018) [60], we used the stratified estimators provided by Cochran (1977) [61] that can be readily used for this purpose. Based on stratified estimators, the estimate of the total proportion of the area in the damaged class can be derived from the confusion matrix using the validation dataset as reference (Table 1), and is given by

$$\hat{\mathrm{p}}={\displaystyle \sum }_{j=1}^2{w}_j·{\hat{\mathrm{p}}}_{\mathrm{j}}$$

(13)

where $w_j$ is the proportion of the map in each map classes (i.e., $w_1$ damaged and $w_2$ undamaged forest), while the variance of esteems of the total proportion of the area in damaged class is:

$$\widehat{\mathrm{Var}}\left(\hat{\mathrm{p}}\right)={\displaystyle \sum }_{\mathrm{j}=1}^3{\mathrm{w}}_{\mathrm{j}}^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_{\mathrm{j}}\right)$$

(14)

On the basis of the confusion matrix we can produce a formal estimation of the damaged area as:

$${\hat{\mathrm{A}}}_{\mathrm{damaged}}={\mathrm{A}}_{\mathrm{tot}} \hat{\mathrm{p}}$$

(15)

Moreover, based on $\hat{\mathrm{p}}$ and $\widehat{\mathrm{Var}}\left(\hat{\mathrm{p}}\right)$ it is possible to calculate a 95% confidence interval for damaged area estimation, that is:

$${\hat{\mathrm{A}}}_{\mathrm{damaged}}·\hat{\mathrm{p}}\pm 2·{\hat{\mathrm{A}}}_{\mathrm{damaged}}·\sqrt{\widehat{\mathrm{Var}}\left(\hat{\mathrm{p}}\right)}$$

(16)

where ${\mathrm{A}}_{\mathrm{damaged}}$ is the mapped damaged forest area.

In addition, we calculated the standard error (SE) and the percentage SE (SE%) of the area estimates as:

$$\mathrm{SE}\left({\hat{\mathrm{A}}}_{\mathrm{damage}}\right)={\mathrm{A}}_{\mathrm{tot}}\sqrt{(w_1{}^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_1\right)+w_2{}^2·\widehat{\mathrm{Var}}\left({\hat{\mathrm{p}}}_2\right)}$$

(17)

$$\mathrm{SE}\%=\frac{\mathrm{SE}\left({\hat{\mathrm{A}}}_{\mathrm{damage}}\right)}{{\hat{\mathrm{A}}}_{\mathrm{damage}}}x100$$

(18)

4. Results

The two algorithms BEAST and CCDC were used to decompose the NBR time series to map abrupt changes due to VAIA windstorm in forests areas. In

Table 2. Results of damaged area estimates (${\hat{\mathrm{A}}}_{\mathrm{damage}; }\mathrm{SE}\left({\hat{\mathrm{A}}}_{\mathrm{damage}}\right);\mathrm{SE}\%)$ and accuracy assessment (OA = Overall Accuracy, PA = Producer Accuracy of damaged class, UA = User accuracy of damaged class; g_mean) obtained by BEAST and CCDC.

Algorithm	Month from the Storm	${\widehat{\mathbf{A}}}_{\mathbf{d}\mathbf{a}\mathbf{m}\mathbf{a}\mathbf{g}\mathbf{e}}$ [ha]	$\mathbf{S}\mathbf{E}\left({\widehat{\mathbf{A}}}_{\mathbf{d}\mathbf{a}\mathbf{m}\mathbf{a}\mathbf{g}\mathbf{e}}\right)$ [ha]	SE%	OA	PA	UA	${\mathbf{g}}_{\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{n}}$
BEAST	1	11,119	26,345.6	236.93	48.2	0.15	0.18	0.17
	2	12,369	23,246.7	187.93	51.4	0.21	0.2	0.2
	3	20,377	22,354.1	109.7	60.2	0.37	0.46	0.42
	4	22,614	199,922.9	88	62.1	0.39	0.46	0.42
	5	27,527	14,802.1	53.77	67.2	0.46	0.49	0.47
	6	36,766	13,149.1	35.76	75.2	0.57	0.7	0.64
	7	38,416	3725.3	9.69	89.7	0.87	0.87	0.83
	8	38,819	405.4	1.04	97.1	0.95	0.95	0.95
	9	40,018	402.1	1	97.8	0.95	0.97	0.96
	10	39,931	346.5	0.87	98	0.95	0.97	0.97
	11	40,126	346.5	0.86	98.1	0.96	0.98	0.97
	12	39,954	238.2	0.6	98.4	0.97	0.98	0.97
CCDC	1	10,203	28,631	280.6	43.5	0.09	0.18	0.09
	2	11,388	28,099	246.7	44.3	0.1	0.1	0.1
	3	13,160	26,560	201.8	46.4	0.12	0.12	0.12
	4	14,268	21,110	148	52	0.18	0.16	0.17
	5	25,349	12,041	47.5	69	0.49	0.49	0.47
	6	32,355	10,295	31.8	75.5	0.59	0.63	0.61
	7	39,204	3254.52	8.3	91.1	0.82	0.9	0.86
	8	38,632	405.4	1.04	97	0.95	0.94	0.95
	9	40,008	402.1	1	97.8	0.95	0.97	0.96
	10	39,929	346.5	0.87	98	0.95	0.97	0.97
	11	40,116	346.5	0.86	98.1	95.8	98.1	0.97
	12	39,951	238.2	0.6	98.4	96.7	98.1	0.97

, we present for each algorithm and for each month after the storm (i.e., 1–12: November 2018–October 2019) the estimated damaged forest area (${\hat{\mathrm{A}}}_{\mathrm{damaged}}$), the SE, SE%, and the accuracy of the maps calculated using the validation dataset. The most accurate results in terms of SE, SE%, OA, PA, UA, and g_mean for both algorithms were achieved between 7 and 12 months after the storm (i.e., May 2019–October 2019) (Table 2 and Figure 4). In fact, from Figure 4, it is possible to see that NBR trajectory of damaged area shifts from the normal seasonality trajectory between April 2019 to October 2019 when the NBR of damaged areas remain always lower than 0.45, while for undamaged areas the trajectories are equal with the ones of the previous years with the NBR that remain always higher than 0.4. The agreement between the maps obtained by the two algorithms increases during the period (May 2019–October 2019) (K > 0.8), while in autumn–winter–early spring (November 2018–April 2019) the two algorithms map differently damaged forest areas (K < 0.35). In Figure 5 we report the maps of VAIA damaged forest area obtained after 8 months (June 2019) with consistent results from both algorithms.

5. Discussion

In this study, we showed that forest area damaged by the VAIA windstorm can be mapped with both the algorithms we tested BEAST, and CCDC. The analyses we performed showed that S2 NBR TS and CDC algorithm (i.e., BEAST and CCDC) are able to map forest areas damaged by windthrows confirming the results, in terms of PA, UA, OA, obtained in Malawi and Austria by Puhm et al. [37], using S2 data a near real-time CDC that combines structural time series model and Kalman filter.

Our results confirm the findings from Dalponte et al. (2020), and Vaglio Laurin et al. (2020) [15,18], which worked in two small study cases within the area hit by the VAIA windstorm. Our results confirm S2 imagery is adequate to map damaged forest area. The most accurate results can be obtained in spring–summer (i.e., after 7 months after the storm), independently of the CDC algorithm used. In fact, in our case it was not possible to produce an accurate map 1–6 months after the storm (November 2018–April 2019) (Table 2). The results we obtained after 7 months are in line with the ones obtained with Landsat bi-temporal approach in France by Haidu et al. (2019) [62] (OA = 86%) in the detection of windthrows in Voges mountains using “dark object” approach and the Disturbance Index, and by Baumann et al. (2014) [63] in European Russia and US using support vector machine classifier (i.e., OA ∼75%). However, after 8 months the accuracy we obtained was higher than the ones reported by the previous mentioned studies.

Analyzing the seasonality and the remainder components of damaged areas NBR TS (Figure 3) we observed a persistent deviation of NBR trajectory from the seasonality between May 2019 and October 2019 (Figure 3). On the contrary, these differences are not visible before April 2019 (Figure 3), this can be explained considering mainly four elements. Firstly, in coniferous forests fallen trees remained green on the ground for a couple of months after the storm [8,15], so differences in spectral trajectories with a pixel size of 10–20 m cannot be immediately detected, while in broadleaves forests (i.e., mainly beech) differences in photosynthetic activities in winter between fallen trees and not damaged trees cannot be detected because trees are leafless [8]. Secondly, snow cover in Alpine regions introduced noise in spectral trajectories, also applying a despike approach in correspondence of NDSI high values. In fact, the snow cover was persistent between autumn and early spring with different timing among different years, introducing noises in spectral trajectory data [64]. Thirdly, optical images acquired in mountains regions [15,65], especially during the winter period, are affected by shadows due to steep slopes which introduce large noises in NBR spectral trajectories that limit the accuracy of BEAST and CCDC methods. In fact, analyzing the omission and commission errors between November 2018 and March 2019 we found that 75% of errors are located in an area with steep terrains, while correctly classified areas are concentrated in the highlands or wide valleys where slope shadows are less present [47]. Fourth, as found by Puhm et al. [37], the accuracy of windstorm damaged areas increased with the incising of time series density (i.e., number of available imagery).

The two algorithms use two different strategies to decompose the NBR TS (i.e., number of continues persistent deviation observation from seasonality to detect changes and classification) and this can be the cause of the disagreement between the results obtained by BEAST and CCDC for the early months after the storm (K < 0.3), while since May 2019 the two algorithms produce consistent results (K > 0.8) with the vast majority of disagreements (98%) located on the borders of the windthrown areas. This is an expected result since pixel located on the border of damaged areas have mixed spectral behaviors between damaged and undamaged pixels [66].

We demonstrated that the map produced by two algorithms can be used to estimate the forest damaged area with a probability-based stratified estimator [67], obtaining an accurate estimate after 7 months from the storm (SE% < 10%) and even more accurate after 8 months (SE% < 2%). Moreover, we considered 7–8 months adequate timing to support decision regarding actions to remove wood materials, and to plan adequate mitigation strategies to limit the expansion of bark beetles outbreaks [68]. Using BEAST and CCDC implemented in GEE we were able for the first time to map all the damaged forest area interested by the VAIA windstorm.

The estimates done by local authorities immediately after the storm [4,47] resulted within the confidence intervals of the estimates we obtained with the two algorithms. However, the VAIA forest damaged area we estimated 7 months after the storm (i.e., from May to October 2019) is slightly smaller than that one reported by local authorities [4]. In fact, we found differences between 4109 ha and 2571 ha for BEAST and between 3321 ha and 2574 ha for CCDC [47]. Such differences can be explained considering that the spatial borders of windthrow patches obtained by satellite data can be substantially different from the perimeters shape mapped in the field [15,19]. In the future, the VAIA map we obtained (i.e., the most accurate one, at October 2019) can be used to support forest strategies regarding the management of instable trees placed on the borders of windthrown areas and future regeneration strategies.

6. Conclusions

Following our results, the forest area affected by the VAIA windstorm is 39,951 ha ± 450 ha. NBR S2 TS is useful to map windthrow damage within 12 months from the storm using CDC algorithms. However, both BEAST and CCDC produced accurate results in detecting damaged area after 7 months of the storm as the ones provided by less computationally intensive methods. In fact, in Alpine regions, the limitation in early mapping of forest disturbances, also if we applied CDC algorithm, are mainly due to slope shadows and snow covers that introduce noises in the TS analysis.

The use of the probability based estimator allows to assess the uncertainty of the estimates of forest damaged area providing a SE% that is mandatory to report the data in the context of international frameworks.

Future research efforts should focus on: (i) testing the procedure with different windthrows; and (ii) evaluating the use of global ecosystem dynamics investigation (GEDI) high resolution laser ranging data to quantify the growing stock volume fallen.