Windthrow Detection in European Forests with Very High-Resolution Optical Data

With climate change, extreme storms are expected to occur more frequently. These storms can cause severe forest damage, provoking direct and indirect economic losses for forestry. To minimize economic losses, the windthrow areas need to be detected fast to prevent subsequent biotic damage, for example, related to beetle infestations. Remote sensing is an efficient tool with high potential to cost-efficiently map large storm affected regions. Storm Niklas hit South Germany in March 2015 and caused widespread forest cover loss. We present a two-step change detection approach applying commercial very high-resolution optical Earth Observation data to spot forest damage. First, an object-based bi-temporal change analysis is carried out to identify windthrow areas larger than 0.5 ha. For this purpose, a supervised Random Forest classifier is used, including a semi-automatic feature selection procedure; for image segmentation, the large-scale mean shift algorithm was chosen. Input features include spectral characteristics, texture, vegetation indices, layer combinations and spectral transformations. A hybrid-change detection approach at pixel-level subsequently identifies small groups of fallen trees, combining the most important features of the previous processing step with Spectral Angle Mapper and Multivariate Alteration Detection. The methodology was evaluated on two test sites in Bavaria with RapidEye data at 5 m pixel resolution. The results regarding windthrow areas larger than 0.5 ha were validated with reference data from field visits and acquired through orthophoto interpretation. For the two test sites, the novel object-based change detection approach identified over 90% of the windthrow areas (≥0.5 ha). The red edge channel was the most important for windthrow identification. Accuracy levels of the change detection at tree level could not be calculated, as it was not possible to collect field data for single trees, nor was it possible to perform an orthophoto validation. Nevertheless, the plausibility and applicability of the pixel-based approach is demonstrated on a second test site.


Introduction
In the last few decades, Europe was hit by a series of heavy storms, such as Vivian and Wiebke in 1990, Lothar in 1999, followed by Kyrill  forest damage in South Germany.It is expected that in the future the frequency of severe storms will further increase in Europe [1,2].These storms damaged areas were often correlated with subsequent insect outbreaks, mainly European spruce bark beetle Ips typographus (L.).In Europe, abiotic wind disturbances and biotic bark beetle infestations are the main factors for causing forest disturbances [2].
Between 1950 and 2000 the approximate average annual storm damage in Europe was 18.7 million m 3 of wood, with most of the storm damage occurring in Central Europa and the Alps [1].In the same period, the average annual wood volume damage by bark beetles was 2.9 million m 3 per year [1].Thus, the affected forest areas need to be localized quickly to assess the damage and to minimize the impact of biotic subsequent damage, as storm damaged trees are providing breeding materials for the insects [1,3].Efficient and cost-effective tools for a fast and reliable detection of the vast windthrow areas are still missing.Previous studies have shown the general suitability of passive optical Earth Observation (EO) sensors for forest disturbances detection [4][5][6][7][8][9][10][11].
These studies mainly apply pixel-based approaches to medium-resolution (MR) data (Table 1).Compared to MR data, high (HR) to very high resolution (VHR) images display much more detail of the Earth's surface, thus enabling the detection of fine-scale tree damage, which MR data is either unable or limited in its ability to detect.Despite these characteristics, VHR images are less often applied for evaluating forest cover change [12], probably related to the high costs involved with data acquisition and data storing for large areas.Only a few recent studies [7,10,13,14] use HR to VHR satellite imagery for forest disturbance monitoring.As more and more VHR satellites are launched into orbit, and costs for VHR data are continuously decreasing, this data can be a sufficient data source to detect mainly small scale disturbances in forests.a broad overview of established CD techniques.For a synopsis of forest disturbance CD studies, the reader is referred to Chehata et al. [12].
Generally, CD techniques can be differentiated in pixel-based analysis (PBA) and object-based image analysis (OBIA) [15,17,26].In OBIA, instead of analyzing single pixels, typically image segmentation is applied first.OBIA is particularly useful when carrying out a CD with VHR imagery, since VHR pixels are in the size of 1 m (measurement unit) and therefore are often smaller than the land surface object (i.e., tree or group of trees) that needs to be detected (mapping unit).Thus, the object size can be scaled to the size of the studied feature [27].A further advantage of OBIA is that per image-object, statistical, and geometrical measures can be computed in addition to spectral information, which can increase change detection accuracy [22,28].
In this paper, we aim to demonstrate that passive VHR EO satellites provide detailed information about the size and location of the storm loss areas.Large damaged areas are usually easily detected and cleared, smaller areas often not.The remaining dead trees can be very problematic for further bark beetle infestation.Thus, in our study we specially focus on detecting small-scale damages to provide fast information for foresters.
Therefore, the images are taken in a narrow period before and after the storm.As such, mainly forest changes caused by the storm are detectable, which is the focus of our interest.Further, we wanted to capture by applying Random Forest feature selection the spectral characteristics, texture, vegetation indices, and layer combinations that are most suitable for detecting changes caused by windthrow.
Furthermore, we wanted to explore the potential of the red edge spectral region for windthrow detection.In our study, we applied VHR RapidEye data comprising a red edge channel that was centered at 710 nm.Recent studies have shown the performance of red edge bands in facilitating important information on the state of vegetation and their benefit for vegetation monitoring [29][30][31].In 2015, Sentinel-2A satellite has been launched.The data is globally and freely available, and within the next couple of months, two identical satellites (Sentinel-2A and 2B) will be in orbit, acquiring imagery of land's spectral properties in up to 10 m spatial resolution every five days, thereby providing dense time series.Sentinel-2 multispectral imager is the first HR sensor to include three red edge bands.RapidEye's red edge band is close to Sentinel-2's band 5 [32], hence our findings can support future windthrow analysis with Sentinel-2 data.
As storms such as Niklas often cause different kinds of damage patterns, ranging from scattered tree damage to compact, homogeneous areas with up to several hectares, we aimed for the identification of two sizes of damaged forest areas: • windthrow areas ≥0.5 ha, and • groups consisting of only few fallen trees (both tree fall-gaps and freestanding groups).

Storm Niklas and Study Area
Storm Niklas hit Bavaria, Germany, on 31 March 2015 and caused severe damage.Niklas reached peak gusts up to 192 km/h in the Bavarian Alps and up to 120 km/h in the lowlands, and was one of the strongest spring storms in the last 30 years [33].The core area affected by the storm was the region between the south of Munich and the Alps.
Niklas caused a heterogeneous damage pattern: mainly small groups and single fallen trees, as well as large-scale storm loss in a temperate forest.We performed our study on two test sites in Bavaria, Germany (see Figure 1), which presented both of the two damage patterns.The first site was selected for method development, while the second site was used for method validation.
The first study area Munich South is situated south of the city of Munich (48 • 0 N, 11 53% of the site is covered by temperate forest, with Norway spruce (Picea abies (L.) Karst) being the dominant tree species (Bavarian State Institute of Forestry: pers.communication).The forests have been hit by the European windstorms Vivian and Wiebke in 1990; since then, forest owners pushed for a conversion into hardwood rich mixed forests [34].
The second study site Landsberg (840 km 2 ) is located west of Munich, near the city of Landsberg (47 • 53 N, 10 • 58 E, 550-820 m a.s.l.).The region is characterized as alpine foothills with a hilly morainic topography.Approximately one third of the area is covered by forests.Norway spruce forests, pure as well as mixed with European beech (Fagus sylvatica L.) are the main forest types (Bavarian State Institute of Forestry: pers.communication).
Niklas caused a heterogeneous damage pattern: mainly small groups and single fallen trees, as well as large-scale storm loss in a temperate forest.We performed our study on two test sites in Bavaria, Germany (see Figure 1), which presented both of the two damage patterns.The first site was selected for method development, while the second site was used for method validation.
The first study area Munich South is situated south of the city of Munich (48°0′ N, 11°32′ E, 530 to 750 m above sea level (a.s.l.)), covering an area of approximately 720 km 2 .The zone is located within the Munich gravel plain and is characterized by very low terrain variation.Approximately 53% of the site is covered by temperate forest, with Norway spruce (Picea abies (L.) Karst) being the dominant tree species (Bavarian State Institute of Forestry: pers.communication).The forests have been hit by the European windstorms Vivian and Wiebke in 1990; since then, forest owners pushed for a conversion into hardwood rich mixed forests [34].
The second study site Landsberg (840 km 2 ) is located west of Munich, near the city of Landsberg (47°53′ N, 10°58′ E, 550-820 m a.s.l.).The region is characterized as alpine foothills with a hilly morainic topography.Approximately one third of the area is covered by forests.Norway spruce forests, pure as well as mixed with European beech (Fagus sylvatica L.) are the main forest types (Bavarian State Institute of Forestry: pers.communication).

RapidEye Data Set
For the windthrow detection, RapidEye scenes before and after the storm were acquired.RapidEye offers five spectral bands at 5 m spatial resolution: blue (B), green (G), red (R), red edge (RE), and near infrared (NIR).Data are delivered in tiles of 25 × 25 km 2 .As five identical satellites are in orbit, in principle a daily revisiting can be achieved (ignoring cloud cover).When satellites need to be pointed, however, substantial off-nadir views can occur.
For both test sites, the pre-storm images were acquired on 18 March 2015, 13 days before the storm.The post-storm image of Munich South was taken on 10 April 2015, 10 days after the storm.The post-storm image of Landsberg was acquired on 19 April 2015, 19 days after the storm.All images are nearly cloud free, but with some haze.The pre-storm images also have some snow patches.Vegetation changes in the image are expected to be mainly caused by the storm event, due to the short period between the acquisitions dates (23 days for Munich South and 32 days for Landsberg).Before mosaicking the single tiles (25 × 25 km 2 ), radiance was transformed into top of atmosphere (TOA) reflectance following BlackBridge [35].A normalization for differences in sun zenith angles in the TOA correction is expected to be sufficient for the detection of abrupt vegetation changes [36].No information about aerosol loads and water vapor contents were available to run a proper atmospheric correction.
Correctly co-registered images are crucial for change detection analysis.Therefore, all scenes were acquired as orthorectified (level 3A) products.No further georectification was necessary, since the images were overlaying accurately with orthophotos.

Reference Data
Field surveys were conducted by Bavarian State Institute of Forestry in the study regions in April 2015, a few days after the storm.Several damaged areas were thereby mapped using Garmin GPSMAP 64 devices.In addition, information was collected about windthrow size and type (e.g., single fallen tree, stand, or large-scale windfall), tree damage (e.g., broken branches, tree split, uprooted trees), and the main affected tree species.Table 3 summarize the main characteristics of these field surveys, which were used for a first validation of the developed method.
In addition, and for large-scale validation, another reference data set was created from orthophotos (see Table 4).For this purpose, false-color infrared orthophotos with 20 cm spatial resolution were used, which were acquired in June and July 2015.The orthophotos were visually interpreted and windfalls (≥0.5 ha) were manually digitized.The delineated windthrow areas were cross-checked with older pre-storm orthophotos and RapidEye scenes.It should be noted that the reference obtained from the orthophoto is potentially biased.The orthophotos were acquired three (Munich South) to four months (Landsberg) after storm Niklas.During this time, some windfall areas had been cleared and adjacent (not damaged) trees were removed for forest protection reasons (Bavarian State Institute of Forestry: pers.communication).Thus, some cleared areas are larger in size than the windthrow damaged areas.

Method/Approach
To map and identify windthrow areas, a two-step approach was applied.First, a bi-temporal object-based change detection was performed, applying a feature selection and classification with Random Forest (RF).Image segmentation was carried out with large-scale mean shift (LSMS) algorithm.In a second step, smaller windthrow areas were detected at pixel-level using a hybrid technique.The overall approach is outlined in Figure 2; the four main parts of the approach, described in Sections 2.4.1, 2.4.2, 2.5.1-2.5.4 and 2.6 were: It should be noted that the reference obtained from the orthophoto is potentially biased.The orthophotos were acquired three (Munich South) to four months (Landsberg) after storm Niklas.During this time, some windfall areas had been cleared and adjacent (not damaged) trees were removed for forest protection reasons (Bavarian State Institute of Forestry: pers.communication).Thus, some cleared areas are larger in size than the windthrow damaged areas.

Method/Approach
To map and identify windthrow areas, a two-step approach was applied.First, a bi-temporal objectbased change detection was performed, applying a feature selection and classification with Random Forest (RF).Image segmentation was carried out with large-scale mean shift (LSMS) algorithm.In a second step, smaller windthrow areas were detected at pixel-level using a hybrid technique.The overall approach is outlined in Figure 2; the four main parts of the approach, described in Sections 2.4.1 to 2.6, were: 1. calculation of input layers,   Flowchart of the developed windthrow analysis procedure.Input features for Random Forest are highlighted (orange).Abbreviations: TOA (Top Of Atmosphere correction), VI (Vegetation Indices), MAD (Multivariate Alteration Detection), SAM (Spectral Angel Mapper), TCC (Tasseled Cap Coefficients), DI (Disturbance Index), LSMS (Large-scale Mean Shift), NDVI (Normalized Difference Vegetation Index), NIR (Near Infrared), R (Red), and RE (Red Edge).
At pixel scale, 175 feature layers were generated (Table 5), belonging to three groups: (i) spectral input layers; (ii) transformation-based input layers; and (iii) textural input layers.Additional details are given in Appendix A. From the 175 input layers of groups (i) to (iii), seventeen statistical measures were calculated per image-object (Table 6), leading to 2975 object variables (175 × 17).For each object, five additional geometrical variables were calculated.These geometrical variables are listed in Table 6.This resulted in 2980 features per bi-temporal image-object.

Segmentation
For image segmentation, the large-scale mean shift (LSMS) algorithm was used [41,42].LSMS is provided within the Orfeo ToolBox (OTB) 5.0.0 [43].The mean shift procedure, a non-parametric clustering, was originally proposed by Fukunaga and Hostetler [44].For image segmentation, the preand post-storm images were trimmed.Extreme values (beneath 0.5 percentile and above 99.5 percentile of all pixel values) were removed and images were linearly rescaled to values between 0 and 255.
LSMS requires the setting of three parameters: spatial radius (hs), range radius (hr), and minimum region size (ms).hs is the spatial radius of the neighborhood, hr is defining the radius in the multispectral space, and ms the minimum size of generated image objects.If a region is smaller than the ms, the region is merged to the radiometrically closest neighboring object [43].
As we were interested in forest change, we tried to find the best combination of input layers to get objects well delineating changed forest areas.Segmentation tests were performed with both spectral bands and VIs, but we choose to perform the segmentation with spectral bands only.This was done so as not to distort the Random Forest feature selection result, which might otherwise rank the VI high, which was used for segmentation.LSMS was applied to the following layer combinations: The results were visually compared to the pre-and post-storm images as well as to the orthophotos.The best results were obtained by using one difference layer calculated from the mean of the difference layers R, RE, and NIR (Figure 3).Only this segmentation is considered in the remaining study.
Additional trial-and-error tests (not discussed here) were conducted to assess the optimal settings for hs, hr, and ms.For the test site Munich South, following Boukir et al. [45], we decided to set hs to half the standard deviation (SD) of all image digital numbers and for radiometric radius hr a quarter of the SD was used.Minimum object size (ms) was set to 10 pixels (250 m 2 ), representing a small group of trees (see Table 7).The results were visually compared to the pre-and post-storm images as well as to the orthophotos.The best results were obtained by using one difference layer calculated from the mean of the difference layers R, RE, and NIR (Figure 3).Only this segmentation is considered in the remaining study.
Additional trial-and-error tests (not discussed here) were conducted to assess the optimal settings for hs, hr, and ms.For the test site Munich South, following Boukir et al. [45], we decided to set hs to half the standard deviation (SD) of all image digital numbers and for radiometric radius hr a quarter of the SD was used.Minimum object size (ms) was set to 10 pixels (250 m 2 ), representing a small group of trees (see Table 7).For Landsberg, these values had to be adjusted to take into account the finer granularity of the forest cover.The best segmentation result was obtained by setting hs to a third of SD and hr to a sixth of the SD of the image digital numbers.ms has been set to five pixels (125 m 2 ), since the forest is more mixed (and more heterogeneous) compared to the Munich South area.
To reduce the amount of data in the subsequent analysis, only objects that were located within forests were kept.For this purpose, a forest mask was used [46].The mask was buffered with 10 m to possibly capture windfall areas on the edge of the forest mask.For the classification and feature selection we used the Random Forest (RF) classifier [47].This machine learning algorithm consists of an ensemble of decision trees and is currently widely used in remote sensing [25,27,[48][49][50].RF was chosen for classification as the algorithm can deal with few training data, multi-modal classes, and non-normal data distributions [47,51].Some further advantages of RF are the included bootstrapping which provides the relatively unbiased "out-of-bag" (OOB) classification results and the calculation of feature important measures such as Mean Decrease in Accuracy (MDA).This important ranking can be used for feature selection [29, 50,52,53].More details about the algorithm can be found in Pal [54], Hastie et al. [51], and Immitzer et al. [55].
In this study we used the RF implementation randomForest [56] in R 3.2.3[57].For the modelling we set the number of trees to be grown in the run (ntree) to 1000.The number of features used in each split (mtry) was set equal to the square root of the total number of input features [47,56].

Training Data
For the training of the supervised RF classifier, 25 reference polygons for class "windthrow" had been manually selected in the RapidEye imagery.In addition, 25 polygons representing the class "non-windthrow" were selected, such as intact forest, water bodies, fields, shadows, and urban areas.Examples are provided in Figure 4.The RF classifier is sensitive to the number of input data, but as we did not know the amount of windthrow areas, we chose to have equal proportions of training data in the two classes.

Feature Selection
To reduce the number of features in the final RF model, we started with a model including all input features and made a step-wise reduction by removing each time the least important variable based on the MDA values [50,53,58].MDA importance values were recalculated at each step.At the end, the model with the lowest number of input features which reached min 97.5% of the maximum OOB overall accuracy was used.

Classification Margin
Besides indicating the most voted class of each object (majority vote), the frequency of the most common and the second most common class was calculated.By subtracting these two values, a measure for classification reliability can be derived [29, 50,59].In case that the RF model is certain about the class assignment, most trees will be classified in the same class.Consequently, this will result in a high difference of the classification frequency compared to the second voted class.In the same logic, if two classes are more or less equally often assigned to a given object, the difference is low, indicating a small confidence.

Feature Selection
To reduce the number of features in the final RF model, we started with a model including all input features and made a step-wise reduction by removing each time the least important variable based on the MDA values [50,53,58].MDA importance values were recalculated at each step.At the end, the model with the lowest number of input features which reached min 97.5% of the maximum OOB overall accuracy was used.

Classification Margin
Besides indicating the most voted class of each object (majority vote), the frequency of the most common and the second most common class was calculated.By subtracting these two values, a measure for classification reliability can be derived [29,50,59].In case that the RF model is certain about the class assignment, most trees will be classified in the same class.Consequently, this will result in a high difference of the classification frequency compared to the second voted class.In the same logic, if two classes are more or less equally often assigned to a given object, the difference is low, indicating a small confidence.
After classification, for each as 'windthrow' classified image-object the classification security was calculated and grouped in five classes (see Table 8).The regrouping of results permits in a subsequent step to merge adjacent polygons belonging to the same security class.

Pixel-Based Change Detection for Identifying Smaller Groups of Fallen Trees
To identify changes smaller than the smallest segment defined by ms (Table 7)-less than 250 m 2 for Munich South and 125 m 2 for Landsberg-a pixel-based classification was performed.We intended to do so to provide the foresters an additional map at hand, which they could use to faster detect small-scale damage (e.g., some few individual or clusters of damaged trees).This is important because in the region we often had problems with the undetected damaged trees providing nesting material for subsequent bark beetle infestations.
For this, only the two most important predictor layers of the RF model were used as input data.Therefore, the layers with the highest number of counts (several statistical metrics) were chosen among those highest ranked in the RF feature selection.Correlated variables can influence accuracy importance measures (like MDA or mean decrease in Gini (MDG)), but are able to deliver reliable results; MDA is more stable when correlations among the predictor variables exist [55].
Different CD methods have been tested, which are fast and relatively easy applicable, such as difference images, SAM and MAD.The best results were achieved with a combination of SAM and MAD, yielding satisfactory results also in areas with haze and shadows.For this purpose, MAD and also SAM was calculated on each of the layer stacks of the two most important variables of the preand post-storm images.The two generated MAD layers and the generated SAM layer were then joined to form a layer stack of three change layers.

Validation
The quality of the object-based change detection was evaluated against the validation data set (see Table 5).Contingency matrices were calculated for the two binary variables, distinguishing between four possibilities [60,61]  The case "true negatives" did not apply in our analysis since only windthrow areas have been digitized in the validation data set.
From the contingency matrix, we calculated sensitivity (Equation ( 1)) [62].Sensitivity, also known as true positive rate, is defined as: This criterion indicates the percentage of segments that have been properly classified as windthrow area, subject to the totality of all by the algorithm recognized windfalls.
For the pixel-based approach, visual validation with the orthophotos was conducted, since it was not possible to map all single fallen trees in the test areas.

Classification Results
The results of the object-based windthrow detection for areas greater than 0.5 ha were satisfactory.An exemplary change detection result of the study site Munich South is shown in Figure 5.The contingency table for the study site Munich South is given in Table 9.  Sensitivity for the test area Munich South was 93%, indicating a high conformity between the classified windfalls and the validation data set.The detected windthrow areas also showed a good  Sensitivity for the test area Munich South was 93%, indicating a high conformity between the classified windfalls and the validation data set.The detected windthrow areas also showed a good surface sharpness.Of the 295 areas 71% were almost congruent (with size deviations between 1 and 3 m), 13% of the surfaces were smaller in size than the actual windfall extent, 15% were larger in size.Twenty-four areas were wrongly classified as windfall.These areas often belonged to security classes 1 or 2 (see Table 8), and were often adjacent to bright gravel pits or to fallow grounds.This indicates that misclassification was mainly caused by over-exposure from neighboring bright areas.
Similar results were obtained for the Landsberg area; validation results for areas larger than 0.5 ha are listed in Table 9.The sensitivity for this second area was 96%, thereby demonstrating a good match between the classified windthrow areas and the validation data.The detected windthrow areas had a good congruency, however, they were less good than for the Munich South study site.Approximately 51% had an overlapping extent, 47% of the detected areas were smaller than the validation areas, and 2% were larger.
The reason for underestimating the extent of the detected areas is probably caused by the acquisition date of the orthophotos.They were acquired four months after the storm event.Therefore, the actual windfall areas could not be identified perfectly.Instead, what was mapped corresponds to the forest state at the time of orthophotos acquisition.At that time, some of the damaged areas were already cleared.The cleared sites were often greater than the actual affected area, due to additional timber outtake for protective measures.This was also confirmed by visual inspection of RapidEye and Google Earth images [63], which were acquired closer to the storm event-these images reveal a better agreement.
In Landsberg, only one area was incorrectly classified as windfall, whereas four windthrow areas were not detected.These areas were mainly mixed forest stands with a relatively strong increase in vegetation cover coupled with the release of understory [64] and the onset of spring and during the 32 days between the storm event and image acquisition.

Selected Features
Not all features were equally important for the detection of windthrow areas using the RF model.The nine most important variables for the Munich South study area are listed in Table 10.Four of the selected variables were directly based on the Plant Senescence Reflectance Index (PSRI); a fifth one (RE reflectance) entered into the calculation of PSRI.Two of the remaining selected variables were associated with the NIR channel of the pre-storm image.The last variable was based on Haralick texture.The strong contribution of PSRI is interesting because the index had been initially developed to determine the onset of senescence and to monitor the vegetation's vitality [65].In our study, PRSI was found to be very important in detecting forest disturbances caused by windfalls.This demonstrates the suitability of the index and the RE channels for the detection of stressed and damaged vegetation.
Besides the PSRI, the NIR band of the post-storm image was selected twice among the important features.As PSRI is calculated with bands of the visible spectra (B, R, and RE) and does not consider the NIR region, it is sensitive to plant pigments.NIR reflectance provides information about the cell structure of the plant and is widely used in vegetation monitoring [66].The post-NIR reflectance indicated changes in leaf structure of the damaged trees.
Similar to the Munich South area, the PSRI was found to be one of the most important variables for the Landsberg area, and the variables based on RE reflectance appeared several times (Table 11).The RE channel of the pre-storm image was ranked highest.Three variables were based on the difference of Normalized Difference Red Edge Blue Index, also closely related to the RE channel.Only the variable difference of 2nd Modified Soil Adjusted Vegetation Index was based on NIR channel (see Table 11).No texture measure was observed to be amongst the most important variables.Together, Table 10 (Munich South) and Table 11 (Landsberg) show that mainly vegetation indices and spectral bands were important for windthrow detection.Specified Tasseled Cap Coefficients (sTCC) and DI features were not selected, indicating that these transformations perform better with sensors they originally have been developed for, such as, sensors with shortwave infrared (SWIR) bands (e.g., Landsat).Textural features are in principle suitable for object-based change detection [67], but played a minor role compared to the variables based on spectral information.

Detection of Smaller Groups of Fallen Trees Using Pixel-Based Change Detection
For the pixel-based CD, the key input layers of the RF feature selection, PSRI and NIR (see Table 10) were used.In this MAD-SAM combination, both the larger windfalls as well as smaller windfall patches appear in blue (Figure 6).To test the applicability of the pixel-based approach, the MAD-SAM combination was accordingly computed for the Landsberg study area and led to similar results as those observed for Munich South.
Visual comparison of the change images with orthophotos and validation data was conducted.Figures 6-8 illustrate the result of the PBA approach.The utilized layer combination also shows single fallen trees.Furthermore, the PBA is fast to compute.
Applying SAM and MAD yields better results in comparison to a simple layer stack of the selected features in the areas covered by haze.Since MAD [39] and SAM, given the same number of bands [68], can handle data acquired with different sensors, this approach could be suitable if there is no pre-and post-storm image available from the same sensor, or if a cloud free scene was acquired faster by a different sensor.However, this case is subject to further research.Most of the damage in Munich South occurred in the western and northern part of the study area.The damage happened often, but not only, in forests where canopy gaps already existed.These canopy gaps were present for different reasons, like forest thinning (in these areas often skid roads are visible) and forest conversion.The forest conversion is done to adapt the forest to the local conditions and to minimize windthrow damage as well as biotic damage (Bavarian State Institute of Forestry: personal communication), changing the forest type from spruce monoculture to mixed forest.Furthermore, the area was hit by winter storms Elon and Felix in January 2015, which minor damage compared to storm Niklas.Less forest damage occurred in the southeastern part of the study area (see Figure 10), which had less canopy gaps.In the Landsberg study area mostly spruce and mixed forests were affected.Overall, this region was also affected by windfall, but the individual windthrow patches showed smaller sizes as in Munich South (see Figure 9).In both study sites, mainly smaller areas, with a size up to 0.9 ha were blown down.Only a few areas were larger than 2.0 ha.Often these areas were closely located to each other, only being separated through a group of remaining trees.
The windthrow detection performed well for both flat conditions (e.g., Munich South) and in hilly terrain (e.g., Landsberg).In mountainous areas, however, the classification outcome is expected to be less accurate, as illumination conditions may vary sharply between adjacent sites.Difficulties are also expected regarding the detection of damaged leaning trees.In our test sites, the spectral differences between leaning trees and intact trees were only minor, prohibiting reliable results.
Both data sets featured the effects of snow, shadows, and haze.Despite these difficult conditions, the algorithms could identify practically all windthrow areas with good area sharpness.Shadow had minor influence in the pixel-based analysis.This demonstrates that the methods are robust enough for handling difficult and varying conditions during image acquisition as long as the area is not covered by clouds.10), which had less canopy gaps.In the Landsberg study area mostly spruce and mixed forests were affected.Overall, this region was also affected by windfall, but the individual windthrow patches showed smaller sizes as in Munich South (see Figure 9).In both study sites, mainly smaller areas, with a size up to 0.9 ha were blown down.Only a few areas were larger than 2.0 ha.Often these areas were closely located to each other, only being separated through a group of remaining trees.
The windthrow detection performed well for both flat conditions (e.g., Munich South) and in hilly terrain (e.g., Landsberg).In mountainous areas, however, the classification outcome is expected to be less accurate, as illumination conditions may vary sharply between adjacent sites.Difficulties are also expected regarding the detection of damaged leaning trees.In our test sites, the spectral differences between leaning trees and intact trees were only minor, prohibiting reliable results.Both data sets featured the effects of snow, shadows, and haze.Despite these difficult conditions, the algorithms could identify practically all windthrow areas with good area sharpness.Shadow had minor influence in the pixel-based analysis.This demonstrates that the methods are robust enough for handling difficult and varying conditions during image acquisition as long as the area is not covered by clouds.

Comparisson with Other Studies
Our findings correlate with other studies applying HR [12] and VHR satellite data [14] for windthrow change detection.Chehata et al. [12] implemented an object-based approach using HR Formosat-2 data on a 60 km 2 study site to detect windthrow.In Chehata et al.'s study, 87.8% of the disturbed areas were identified, with green band shown to be the most useful.It should be noted that Formosat-2 does not have a red edge band.The study performed by Elatawneh et al. [14] detected forest cover loss with OBCD after a windfall with accuracies ranging between 87% and 96% on a 130 km 2 area.RapidEye data were applied, and the red edge band was labeled as best performing band, as a result of a visual pre-study [69].However, our study included two larger test sites (720 km 2 and 840 km 2 ) compared to the studies by Chehata et al. [12] and Elatawneh et al. [14] and various predictor variables (texture, spectral, statistical, and geometrical) were used to test for windthrow detection with Random Forest.
NDVI is widely used for vegetation and disturbance monitoring.However, in our study, the NDVI was not ranked highly as a suitable predictor layer.This is in line with studies applying HR or VHR for both disturbance [12,14] and severity [7] detection after windthrow.This further confirms the findings of other studies [5,6,17] that used medium resolution data for forest disturbances monitoring.
Some uncertainties exist in our study.First, we used two RapidEye datasets with 23 and 32 days of difference, and we performed our analysis under the explicit assumption that the changes we discovered were mainly caused by windthrow, due to the narrow time window.As salvage logging is the common harvesting type in our study region, we cannot, however, exclude the possibility that some of the detected changes were caused by forest thinnings, which would present similar spectral changes [7].Secondly, we focused only on the detection of changed forest areas.The scope of our study did not consider the identification of several different types of land change, like in other studies [9,21].In the future, different types of disturbances should be included in the approach.Further, while the impact of changes in illumination and phenology between pre-and post-storm images in our study were minor, they could generally lead to misinterpretations.

Conclusions
In our study, we applied a two-step approach to identify windthrow areas in a temperate forest.To detect damaged areas greater than 0.5 ha, we developed an object-based method using Random Forest (RF) classification that demonstrated high accuracies, which were validated with independent reference data.The built-in feature selection identified similar predictor variables in both areas.After having identified the most important input layers, as derived from the object-based approach, we developed the subsequent pixel-based method to visualize windthrow areas on pixel level.With this, we intended to give the foresters an additional map to have at hand, which they could use to faster detect small-scale damage (e.g., just some few damaged trees or canopy gaps).This is important because in the past, problems occurred with undetected clusters of damaged trees which provided breeding material for bark beetles.
Our study proved the applicability and robustness of RF classifier for windthrow detection at object-and pixel-level and its benefit for ranking features according to their importance.As RF is easy to implement and to handle, we recommend its use for windthrow detection.
According to our results, vegetation indices proved suitable measures for the detection of windthrow.In particular, vegetation indices achieved better results compared to individual spectral bands or texture.Plant Senescence Reflectance Index (PSRI) and Normalized Difference Red Edge Blue Index (NDREDI), both applying the red edge and blue spectral bands, were identified as the most applicable indices.The use of sensor data providing information in red edge and blue spectral bands, and vegetation indices including those bands, is therefore potentially beneficial, but any definite suggestion would require further research under different test conditions.In future work, the use of high resolution Sentinel-2 data is potentially beneficial for windthrow detection.However, it needs to be tested whether small-scale changes can be captured with the Sentinel-2 data, which have a coarser resolution than RapidEye.Otherwise, Sentinel-2, especially considering its three red edge bands as well as its shortwave infrared bands, has the necessary spectral band setting to perform well in large-scale windthrow monitoring.Our developed method as well as existing forest disturbance monitoring methods based on Landsat data need to be further tested with Sentinel-2 data.

Figure 1 .
Figure 1.False color composites of RapidEye scenes (band combination: near infrared-red-green) acquired after storm Niklas covering the two test sites Munich South (west) and Landsberg (east) in Bavaria, Germany.

Figure 1 .
Figure 1.False color composites of RapidEye scenes (band combination: near infrared-red-green) acquired after storm Niklas covering the two test sites Munich South (east) and Landsberg (west) in Bavaria, Germany.

2 .
image segmentation, 3. RF feature selection and object-based change classification, and 4. development of pixel-based approach

•
stacked pre-and post-storm TOA reflectances (10 layers) • difference images of the TOA reflectances (5 layers) • difference images of the three spectral channels R, RE, and NIR (3 layers) • mean of the three difference images (R, RE, and NIR) (1 layer)


stacked pre-and post-storm TOA reflectances (10 layers)  difference images of the TOA reflectances (5 layers)  difference images of the three spectral channels R, RE, and NIR (3 layers)  mean of the three difference images (R, RE, and NIR) (1 layer)

Figure 3 .
Figure 3. Detail of Munich South scene showing input (a) and resulting image segmentation (b).The input layer (a) consists of the mean values of the three difference layers from red, red edge, and near infrared bands.Darker regions within the forested area indicate windthrow.

Figure 3 .
Figure 3. Detail of Munich South scene showing input (a) and resulting image segmentation (b).The input layer (a) consists of the mean values of the three difference layers from red, red edge, and near infrared bands.Darker regions within the forested area indicate windthrow.

Figure 4 .
Figure 4. Example reference polygons showing different windthrow areas (yellow polygons) as well as non-windthrow polygons (green polygons).Polygons are overlaid on a false color composite orthophoto (band combination: near infrared-red-green) acquired four months after the storm.

Figure 4 .
Figure 4. Example reference polygons showing different windthrow areas (yellow polygons) as well as non-windthrow polygons (green polygons).Polygons are overlaid on a false color composite orthophoto (band combination: near infrared-red-green) acquired four months after the storm.

ForestsFigure 5 .
Figure 5. Object-based change detection results for an exemplary area in Munich South.Polygons are overlaid on false color composite of post-storm scene (band combination: near infrared-red-green).(a) Reference windthrow areas greater than 0.5 ha; (b) detected windthrow areas greater than 0.25 ha, the minimum mapping unit of the segmentation, with classification margins.

Figure 5 .
Figure 5. Object-based change detection results for an exemplary area in Munich South.Polygons are overlaid on false color composite of post-storm scene (band combination: near infrared-red-green).(a) Reference windthrow areas greater than 0.5 ha; (b) detected windthrow areas greater than 0.25 ha, the minimum mapping unit of the segmentation, with classification margins.

Figure 6 .
Figure 6.Detail of false color composite (band combination: near infrared-red-green) pre-storm RapidEye scene (a); post-storm RapidEye scene (b) and pixel-based change detection result based on RapidEye data (c) for Munich South.Windthrow areas appear in blue.

Figure 7 .
Figure 7. Examples of detected windthrow areas depicted with false color composite orthophoto (band combination: near infrared-red-green) (a,c,e) and pixel-based result based on RapidEye data (b,d,f) for Munich South.Windthrow areas appear in blue.Scattered tree damage (a-d), the damage often occurred in recently thinned forest (a) where skid trails are visible.Small groups of affected trees could also be detected (e,f).Here small thrown stands, comprising eight (left) and three (right) trees, respectively, were detected.

Figure 7 .
Figure 7. Examples of detected windthrow areas depicted with false color composite orthophoto (band combination: near infrared-red-green) (a,c,e) and pixel-based result based on RapidEye data (b,d,f) for Munich South.Windthrow areas appear in blue.Scattered tree damage (a-d), the damage often occurred in recently thinned forest (a) where skid trails are visible.Small groups of affected trees could also be detected (e,f).Here small thrown stands, comprising eight (left) and three (right) trees, respectively, were detected.

Figure 10 .
Figure 10.Occurrence of windthrow areas ≥0.5 ha (green areas) in the Munich South study area.As background is the false color composite (band combination: near infrared-red-green) post-storm RapidEye scene depicted.

Figure 10 .
Figure 10.Occurrence of windthrow areas ≥0.5 ha (green areas) in the Munich South study area.As background is the false color composite (band combination: near infrared-red-green) post-storm RapidEye scene depicted.

Table 1 .
Selection of essential previous studies for forest disturbance detection.Most studies used medium resolution (deca-metric) satellite data and pixel-based approaches on temperate and boreal forest ecosystems.Abbreviations: TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), • 32 E, 530 to 750 m above sea level (a.s.l.)), covering an area of approximately 720 km2.The zone is located within the Munich gravel plain and is characterized by very low terrain variation.Approximately

Table 2
summarizes the main characteristics of the acquired RapidEye data.

Table 2 .
Main characteristics of the acquired RapidEye data used to detect windthrow areas caused by storm Niklas on 31 March 2015 in South Bavaria.

Table 3 .
Number of reference polygons obtained from field surveys per affected forest type.

Table 4 .
Number of reference polygons (≥0.5 ha) obtained from orthophoto interpretation.The size distribution of the affected windthrow areas is indicated.

Table 6 .
Computed statistical and geometrical variables calculated for each image object.The seventeen statistical variables were calculated per input layer (175) whereas only one set of (five) geometrical variables was derived.

Table 7 .
LSMS segmentation parameters for the two test sites Munich South and Landsberg.The two test sites were segmented using a single layer input image (average of red, red edge, and near infrared channels).

Table 8 .
Overview of the windfall securities classes.

Table 9 .
Contingency table of binary classifier for Munich South and Landsberg.

Table 9 .
Contingency table of binary classifier for Munich South and Landsberg.

Table 10 .
Random Forest model with the most important object variables for the Munich South study site.The last column indicates the value of the Mean Decrease in Accuracy (MDA).

Table 11 .
Most important object variables in the Random Forest model for the Landsberg study site.The last column indicates the value of the Mean Decrease in Accuracy (MDA).
21rests 2017, 8,2118 of 27 canopy gaps were present for different reasons, like forest thinning (in these areas often skid roads are visible) and forest conversion.The forest conversion is done to adapt the forest to the local conditions and to minimize windthrow damage as well as biotic damage (Bavarian State Institute of Forestry: personal communication), changing the forest type from spruce monoculture to mixed forest.Furthermore, the area was hit by winter storms Elon and Felix in January 2015, which caused minor damage compared to storm Niklas.Less forest damage occurred in the southeastern part of the study area (see Figure