DETECTING AND EVALUATING DISTURBANCE IN TEMPERATE RAINFOREST WITH SENTINEL-2, MACHINE LEARNING AND FOREST PARAMETERS

Earth observation via remote sensing imagery provides a fast way to define alteration levels. In this work 12 stands of AraucariaNothofagus forests were selected in southern Chile, which represented four alteration levels: (i) None (ii) Low (iii) Medium and (iv) High. The stands were surveyed measuring 379 field plots and Google Earth Engine was used to collect a composite of Sentinel-2 images over a one-year range, from June 2019 to June 2020. The following approaches were tested: (i) aggregating the normalized difference vegetation index (NDVI) of each image and selecting the 95 and 99 percentile values of NDVI for each pixel; (ii) creating a composite imagery with best pixels over one year timeline using NDVI as weighting factor and NDVI value band itself (NDVI) – this is similar to the 99 percentile in the previous point, but with maximum values of NDVI; (iii) aggregating the composite as in the previous approach, but using the full spectral information of Sentinel-2 and then random forest machine learning for classification over alteration areas with k-fold validation with k=5. Results show that the 95 and 99 percentile of NDVI values from approach (i) do not discriminate the four classes correctly. The maximum NDVI from approach can distinguish all four classes. It must be noted through that statistical significance does not necessarily imply a strong practical significance; medium and high alterations have very similar NDVI distributions. Random forest results provided an F-score for each class higher than 80% except for the “medium alteration” class. * Corresponding author


INTRODUCTION
The Araucaria-Nothofagus forests are located in the southern part of Chile. The Araucaria araucana [Mol.] K. Koch, in association with other Nothofagus spp. are constantly affected, either by natural disturbances such as natural wildfires, alterations of volcanism, snow avalanches, droughts, landslides or wind and geological instability (Heusser et al., 1988;Burns 1993;Puchi et al., 2021). As well as by anthropogenic disturbances, such as forest fires caused by human actions, harvesting of edible seeds, firewood, and timber extraction (in case of Nothofagus spp.), domestic livestock and introduction of exotic species e.g., Pinus contorta Dougl. Ex Loud. (Peña et al., 2008;González et al., 2011;Drake et al., 2012). Therefore, this high dynamism in these forest ecosystems creates the need to develop tools that can identify and quantify such disturbances and analyze these changes at a high spatial scale and in temporal scenarios. A general approach to identify and determine the degradation of a forest is through a standard forest inventory survey using dasometric parameters. The forest structure can be assessed by a dominant tree height (m), canopy cover (%) or either by a total aboveground tree biomass (Mg ha -1 ) at stand level (Coops et al., 2020), allowing a quantitative comparison between forests stands with different level of degradation. However, the conventional method through field measurements is expensive, time-consuming and with some restriction to measuring in remote areas with lack of accessibility, being a not suitable approach in a large spatial scale (Soenen et al., 2020). Remote sensing is an alternative to optimize resources and time to identify and quantify such changes by different natural or human forest disturbances (Schultz et al., 2016;Senf et al., 2017;Bullock et al., 2020). Nowadays, there are a wide range of satellite images with a high spatial resolution, which are free and open access data worldwide (Zhu et al., 2019). On the other hand, spectral bands and vegetational indices, which can be extracted from satellite images e.g., normalized difference vegetation index (NDVI), is a common approach used to assess disturbances on forest vegetation (Lambert et al., 2013;Hojas-Gascón et al., 2015;Lastovicka et al., 2020), the advantage of this vegetation index is the capacity to detect anomalies changes pattern on small scales, which can be related to several biotic or abiotic agents, being identified at pixel level in the vegetation by positive or negative changes (Nolè et al., 2018). The NDVI is highly related to the aboveground tree biomass (and to other forest parameters), but it could be conditioned on the season when the analysis is carried out. Zhu & Liu (2015) mentioned that NDVI in the fall season shows stronger correlation to aboveground tree biomass than in the peak season due the saturation issue in canopies with a high density. Therefore, seasonality is an important aspect to be considered when correlating NDVI and the aboveground tree biomass, especially because the peak season contains information that comes mostly from the top of the canopies, which is sensitive to the presence or absence of the leaves. Therefore, it is important to collect a high number of images that can represent the entire year, creating a single image in order to avoid saturation or lack of vegetation due seasonality and compare different areas contemporaneously.
The goals of this study can be summarized to two main points. First is the use of Sentinel-2 optical imagery data to compute NDVI using different percentiles (95, 99 and maximum NDVI value) in 12 Araucaria-Nothofagus forests stands in southern Chile divided in four different alteration levels (three stands for each alteration level) and identify which percentile can recognize all alteration levels. Secondly analyze if NDVI and forest parameters such as tree density (N ha -1 ), basal area (m 2 ha -1 ), stem volume (m 3 ha -1 ) and total aboveground tree biomass (Mg ha -1 ), can distinguish between different levels of alterations using ground truth values.

Study area
The study area was located in the North-East part of the Araucania region between the provinces of Malleco and Cautin. The 12 forest stands were mainly distributed in three communes and the areas of the stands were relatively similar with a range of 5.4 and 7.9 ha (table 1), being a total area of sampling of 76.0 ha. The climate is defined as cool temperature with a humid regime and Mediterranean tendency. The temperature varies between 19.7 °C and -1.1 °C and the mean annual precipitation is 2,470 mm (AGRIMED, 2017). On the other hand, the soils are derived from volcanic ashes, being relatively deep soils (~ 100 cm depth), its texture varies from clay loam to sandy and its drainage is moderate (Luzio, 2010

Araucaria-Nothofagus forest stands
The 12 forest stands were selected in January 2021 through a fieldwork prospecting and interviews to local people, considering mostly anthropogenic factors such as forest fires, seed harvesting of A. araucana (edible fruit with a high economical value in the market), livestock, firewood extraction of Nothofagus spp., introduction of exotic species (e.g., P. contorta), land division and land use changes, which were observed in situ in the fieldwork prospection and reported by local people. These activities allowed an objective classification of the stands by level of alteration (table 1). Then, between February and March 2021, the 12 forest stands were measured applying a systematic sampling technique and using horizontal point sampling in 379 field plots. In case of the stem volume and total aboveground tree biomass, these parameters were obtained through allometric equations reported by Drake et al. (2003) and Kutchartt et al. (2021). The assignation and categorization of the forest stands to one of the specific levels of alteration was made by some field measurements to characterize the forest structure through some forest parameters, using as a main indicator the basal area (m 2 ha -1 ) and the canopy cover (%) for a quantitative distinction between stands. In addition to the field measurements, it was taking under consideration the sum of the alteration factors identified in situ as it is described below in each of the alteration level category:

Area
No (almost none) alteration: Two of these stands were located in a protected area, one was located at the Conguillio National Park and the other at the Malalcahuello National Reserve, which are managed by the Chilean forestry service. The other stand was in a private forest called La Fusta. These stands presented a mean basal area of 41.9 m 2 ha -1 and a canopy cover of 87.0%, and in the case of Conguillio and Malalcahuello stands, these have not shown major changes since 1998. We evidenced logging in A. araucana mainly in the past (by stumps).

Low alteration:
The three stands were located in the forest El Naranjo, which belongs to an indigenous community. The mean basal area reported between the three forest stands was 44.3 m 2 ha -1 and the canopy cover presented in this level of alteration was 58%. The main disturbances were related to the edible seeds extraction of A. araucana and livestock (cows). There was evidence of logging in the past.
Medium alteration: Also, as in case of low alteration, the three stands were located at El Naranjo. The mean basal area here was 29.7 m 2 ha -1 and the canopy cover presented in this level of alteration was 42%. The main disturbances observed here were edible seeds harvesting in case of A. araucana, frequent livestock and wood extraction in N. pumilio.
High alteration: The three stands were located at La Fusta forest, of which some areas of this forest have historically undergone various types of damage. The mean basal area was 18.1 m 2 ha -1 and the canopy cover presented here was 42%. The alterations observed here were wood extraction of A. araucana in the past for plywood, extraction of wood of N. pumilio and harvesting of edible seeds of A. araucana. From 2012 the stands presented occasional livestock and controlled seed extraction of A. araucana. These areas experienced a huge forest fire in the season 2001-2002, which burned 20,000 ha in total, and it was caused by lightning, from dry thunderstorms during summer (González & Veblen, 2007).

Methods
Google Earth Engine (GEE) was used to create a composite of Sentinel-2 cloud-free imagery (Gorelick et al., 2017). The process consisted in filtering only images from June 2019 to June 2020 and calculating the well-known normalized difference vegetation index (NDVI) for each image. This stack of about 60 images was processed to create a single image using three approaches: it was reduced to a single NDVI image using (i) 95 th and (ii) 99 th percentile values of the NDVI vector at each cell and (iii) NDVI values were used as quality indicators to composite an image with all bands. The rationale of the first two approaches in that this should keep the highest NDVI value that was recorded over the year for each cell, thus removing any disturbance from clouds and keeping the NDVI signal from any vegetation that is present in each cell. The third approach filters the values to keep for each cell and each band by using the NDVI values. In practice, for each band, out of the n-size vector where n equals to the numbers of images in the year and each vector element is an image, the i-th and j-th cell, where i and j are the row and column indices of the specific cell of the n-th image, the final i-th and j-th cell value will consist in the value from the image where that cell had the highest NDVI value over all n cells. Also, in this case the rationale is that where there is vegetation, NDVI will be used as a proxy for "best" representative cell.
The three final images (NDVI_p95, NDVI_p99 and composite, maximum NDVI value) are then used to assess their ability to predict, for each of the 12 stands, alteration values and two other specific forest parameters that are also related: basal area and aboveground tree biomass. The composite image was further tested in two ways: one by creating an overall NDVI map from the composite, and one by keeping all the bands. All the bands of the composite imagery were assessed to know how much they contributed to defining the alteration class, and also were processed with random forest (RF) to assess how well they can be used to predict alteration levels and which bands are most important for this task (Pirotti et al., 2014;Vaglio Laurin et al., 2016). The RF method (Breiman, 2001) used the "ranger" package in the R software environment, with default parameters. The number of variables to split at each node (mtry) is the square root of the number of variables, thus rounded to 3. The number of trees (ntrees) is 500. All assessments over group differences (alteration levels) are assessed for statistical significance using Wilcoxon U Mann Whitey non-parametric test.

RESULTS
Collecting the imagery over 12 months from June 2019 to June 2020 provided with a total of ~60 NDVI images from both Sentinel-2 vectors (Sentinel-2A and Sentinel-2B).
The sections below are divided with respect to forest parameters that are confronted with Sentinel-2 data: alteration levels, tree density, basal area, stem volume and aboveground tree biomass.

Sentinel-2 vs Alteration
The cell values with NDVI timeseries over the whole year aggregated using 95 th and 99 th percentile shows that alteration values can be discriminated between none and low alteration, and medium and high alteration (figure 3 top part). It can be noted that the 95 th percentile provides significant difference also between medium and high alteration, but with opposite difference with respect to what can be expected, e.g., the median NDVI aggregated with 95 th percentile is higher in high alteration then in medium. The difference is not significant anymore if 99 th percentile is used.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France Furthermore, the NDVI created from composing using NDVI as quality factor (figure 3 bottom plot) shows a NDVI value distribution a bit lower for the high alteration. It can be interpreted as NDVI higher values are better suited than 95 th and 99 th percentiles. We must remember that the NDVI composite using timeseries NDVI values basically means keeping the highest NDVI value that was found in each cell over the whole year. Figure 4 shows the results from using random forest over all the bands available in the Sentinel-2 composite. The variable importance plot reflects what can be seen in the top plot as being significantly related to alteration. The red (Band 4) and a specific NIR-red edge (Band 7) and the short-wave infrared bands provide the most weight in the prediction of the alteration level. Table 2 below summarizes the classification accuracies over the independent validation set at each of the five runs. The medium alteration provides the lower score, probably due to mixing with the high alteration class. The five-fold validation allows to show also how consistent the accuracy metrics are over each independent run.   Table 2. Statistics of accuracy metrics from the five-fold random forest training-validation runs.

Sentinel-2 vs Forest Parameters
Forest parameters considered are tree density (N ha -1 ), basal area (m 2 ha -1 ), stem volume (m 3 ha -1 ) and total aboveground tree biomass (Mg ha -1 ). Table 3 report the coefficient of correlation with the average and median NDVI percentile metrics extracted over each stand. The higher values are seen over the 99 th percentile, for basal area, and also for the other forest parameters, but with lower values.  Table 3. Coefficient of correlation (r) for forest parameters and Sentinel-2 information.
A further indication of results can be seen from the plot in figure 5 below, where each stand mean NDVI value is plotted against basal area and total aboveground tree biomass.  We can see two distinct clusters, the high and medium alteration and the low and none alteration groups. Adjusted r-squared shows a good relationship between the two variables. Stand 3 is the only inconsistent one, as it has a low value of fieldmeasured basal area and total aboveground tree biomass, but the NDVI values are not as low as a linear model would predict. This part should be further investigated. Figure 6 shows the NDVI maps for each stand, showing that stand 3 does seem to have a dense patch of high NDVI values in the center, but lower values in the borders, in particular in the northern part of the stand. The plots in each stand are evenly distributed and therefore representative of the area, so a possible explanation can be that the alteration impacted to lower part of the forest vertical structure. Figure 6. NDVI values from reducing using 99 th percentile.

Forest disturbances detected by NDVI
Disturbances and degradation in forest can be detected through field observations and by simple field measurements. However, to detect changes at a large scale, analysis with different approaches and remote sensing advances are needed. Our findings using satellite data were in the same line with Hojas-Gascón et al. (2015), who detected and quantified differences in forest degradation using Sentinel-2 data and computing NDVI values. In our results, at least in three out of four alteration levels (95 th percentile) were categorized (figure 3), which showed a clear negative trend when the NDVI decreases, the levels of alteration are higher. Hojas-Gascón et al. (2015) could discriminated forest classes only with more than 40% of cover differences, identifying similar difficulties when the classification of forest degradation tends to be finer with more classes of degradation. Surprisingly, the high alteration level showed a NDVI value slightly higher in comparison to medium alteration in 95 th percentile. On the other hand, the 99 th percentile showed similar trend than in 95 th , but with the only difference that the categories of medium and high alteration were not possible to distinguished in 99 th between each other, which is in concordance with the ground truth values obtained in the canopy cover, where both categories obtained a 42% of canopy cover during the field campaign. On the other hand, using the maximum NDVI value (NDVI_p100), all the categories of The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France forest alteration were discriminated, but although they were distinguished by means of a non-parametric test, this does not mean that there is a clear distinction between groups, depending mainly on the type of statistical test applied, observing strong association between the none and low alteration as one group and medium and high alteration as a second group.

Degradation level through spectral indices on a single year
Several authors used NDVI to demonstrate the decline of a forest according to different anthropogenic drivers or through environmental changes (e.g., drought, high temperatures, pests by insects or pathogens) over time (Lambert et al., 2013;Neigh et al., 2014;Recanatesi et al., 2018). However, our study was focused on the same year avoiding a multi-temporal analysis, which normally with the inclusion of a longer time series, it provides more detailed information regarding trends related to climate change disturbances and all types of natural disturbances that can influence in the dynamism of a forest in a long term (Hall Bayer, 2003). However, long time series analysis demand high volume of data and it is time-consuming. Therefore, our approach was focused on the verification of anthropogenic alterations using satellite data in a single year related to the field measurements and observations in situ. As expected, we identified different values in their forest parameters through stands surveyed at the same time (monthyear). On the other side, similar response was found with Sentinel-2 and NDVI values, which also took in consideration the seasonality (> 60 satellite images per year) to avoid that phenological changes could influence in the NDVI saturation (Reed, 2006). Our finding using satellite data have the potential to increase the scale of the study area, reaching a provincial or regional scale, using the field measurements only as ground control points.

Anthropogenic drivers on forest status
Our finding showed that forest under similar environmental conditions, but under different anthropogenic pressures can vary on spectral indices, which were corroborated by the ground truth values obtained in the field campaign between February and March 2021. A clear example is La Fusta forest, which considered similar environmental conditions, but its historical management and damaged, classified stands in the same forest from none to high altered, according to its current status. On the other side, the anthropogenic drivers can be categorized in different levels, considering that some could be easily identified through satellite images in a single year, e.g. forest fires (Morresi et al. 2022) or wind-throw (Piragnolo et al. 2021;Vaglio Laurin et al. 2021), but other anthropogenic drivers which are slowly creating damages are not easy to identify in short terms (single year) through satellite data, such as illegal logging, overgrazing, seed harvesting, among others. Therefore, these types of alterations need to be analyzed in long terms using long time series in order to distinguish them in a properly way. It can be noted that other technologies can support more detailed definition of disturbances. The most commonly used in forest environments are laser scanner sensors working at treelevel (Pirotti et al. 2017) or at plot level (Pirotti et al. 2014) to define also a vertical structure and assess disturbance from tree cutting and felling.

CONCLUSIONS
Sentinel-2 data can be a powerful tool to detect and discriminate a level of alteration into a forest at a large scale. However, in some situations the level of alteration cannot be clearly distinguished. When the alteration difference is not clearly defined or when categories have similar parameters (e.g., tree heights, species composition), their classification becomes very complicated. Therefore, the proposed method would be a valid alternative, only if the forest classes are well defined because the detailed information at this level can be obtained only through field measurements.
Future work needs to address the integration of multiple data sources from optical imagery to improve the temporal density of images. For example, using both Landsat and Sentinel-2 missions can help to provide a better temporal composite. Also testing other vegetational indices and higher spatial resolution through commercial images to see if forest degradation classes can be distinguished in more specific groups and spatialized to a higher detail.