Optimizing Near Real-Time Detection of Deforestation on Tropical Rainforests Using Sentinel-1 Data

: Early Warning Systems (EWS) for near real-time detection of deforestation are a fundamental component of public policies focusing on the reduction in forest biomass loss and associated CO 2 emissions. Most of the operational EWS are based on optical data, which are severely limited by the cloud cover in tropical environments. Synthetic Aperture Radar (SAR) data can help to overcome this observational gap. SAR measurements, however, can be altered by atmospheric e ﬀ ects on and variations in surface moisture. Di ﬀ erent techniques of time series (TS) stabilization have been used to mitigate the instability of C-band SAR measurements. Here, we evaluate the performance of two di ﬀ erent approaches to SAR TS stabilization, harmonic deseasonalization and spatial stabilization, as well as two deforestation detection techniques, Adaptive Linear Thresholding (ALT) and maximum likelihood classiﬁcation (MLC). We set up a rigorous, Amazon-wide validation experiment using the Google Earth Engine platform to sample and process Sentinel-1A data of nearly 6000 locations in the whole Brazilian Amazonian basin, generating more than 8M processed samples. Half of those locations correspond to non-degraded forest areas, while the other half pertained to 2019 deforested areas. The detection results showed that the spatial stabilization algorithm improved the results of the MLC approach, reaching 94.36% global accuracy. The ALT detection algorithm performed better, reaching 95.91% global accuracy, regardless of the use of any stabilization method. The results of this experiment are being used to develop an operational EWS in the Brazilian Amazon.


Introduction
The latest world data on primary forest losses in tropical regions report a 3.8 Mha loss in 2019, where Brazil accounts for one third of this loss, reaching 1.36 MHa of deforested area [1]. The most recent official data from the Brazilian government reveal a significant recent rise in annual rates of deforestation, from 0.75 Mha·yr −1 (2018) to 1.01 MHa·yr −1 (2019), which is an increment of 34% in a single year [2].
In this rapid degrading scenario, Early Warning Systems (EWS), defined as a collection of algorithms and procedures able to identify tree loss or disturbance, on a periodic basis (monthly, weekly or even daily) [3], become an essential tool to initiatives of deforestation reduction. EWS have proved to be a crucial element in reinforcing public policies that have led to significantly reduce the deforestation rates in Brazil [4][5][6], Peru [7], and, more recently, in Indonesia [8]. Nevertheless, cloud cover constitutes a serious obstacle to EWS based exclusively on optical data in tropical forest environments. Table 1. Summary of the most important sources of errors and ambiguity in Synthetic Aperture Radar (SAR) measurements over forested surfaces found in the literature.

Effect
Band
Regarding C-band measurements, backscattering can suffer attenuations between −2.0 and −2.4 dB when crossing dense storm cells [22,23] and can increase from 1.0 to 1.5 dB due to the effect of canopy water content [17]. More recent work reports increases of 0.5 dB after medium to severe rain showers, both in temperate forest [26] and tropical forests [27]. Seasonal changes in water content of the canopy can lead to oscillations of 2.5 dB in ERS-C [21] and to 1.5 dB [26] in Sentinel-1 measurements over mature temperate forests. Slight increases in backscattering (+0.5 dB) due to rain interception were found in a basin-wide survey in the Amazon rainforest [27]. The same study reports a heavy attenuation of backscattering in the first hour after a heavy rain event. L-band measurements can equally be affected by meteorological and moisture variations. Although there is no account for L-band attenuation due to dense storm raincells, Koyama et al. reported a 1 to 2 dB increase in backscattering due to rain interception and to the increase in water content of the canopy [24].
These instabilities can be treated using different approaches. One straightforward approach consists of the mathematical treatment of the time series, using temporal pixel-wise moving windows [26] or harmonic functions [28]. Although this approach is powerful in an operational context, it can be misleading, as it will not account for the interannual variability of rain and temperature regimes. For instance, a harmonic fit trained in a humid year would probably overestimate backscatter values in the next, drier year, and would not deal properly with short-term variations of backscattering due to heavy-rain episodes. A spatio-temporal approach to stabilization was proposed in [29]. In that work, the authors use the regional backscatter P95 percentile value to define a mean forest response for every recorded date and then to normalize all the image pixels by this value. This approach has an added advantage: it will deal with sensor oscillations as those reported in [30]; however, it can be difficult to apply to large regions. Other mitigation techniques make use of global, high-frequency rainfall records to mask measurements made during or after heavy rain episodes [26,27]. Meanwhile, masking using globally available coarse resolution data, such as NASA's Integrated Multi-satellitE Retrievals for GPM (IMERG) 10 × 10 km hourly rainfall can mask too many valid measurements contained in a rainy cell [27]. The use of fine spatiotemporal resolution meteorological data, such as those made available by meteorological radars, and an attenuation model, such as the one proposed in [14], can help to overcome this problem [16]. Unfortunately, fine spatiotemporal resolution meteorological data are far from being widely available over entire tropical basins. Recently, a hybrid random forest model, using meteorological and forest structure data, has been successfully applied to stabilize data over a medium-sized region in Eastern Amazon [31]. Although the proposed method is very promising, it does not account for undocumented variations in S1 calibration, such as those reported in [26], and it can fail to model extremely heterogeneous regions such as the entire Amazon basin.

Deforestation Detection with SAR
Deforestation is formally defined as the permanent change from a given forest area to another type of land cover [32]. This change may be of anthropic origin or not. According to the Food and Agriculture Organization of the United Nations (FAO), a forest is defined as an stretch of land larger than 0.5 hectares covered by trees higher than 5 m, whose canopy represents more than 10% of the surface of the area considered, and which is officially considered as forest [33]. Any event, natural or not, that permanently reduces the capacity of the forest cover to reach 10% of the canopy will be considered deforestation. According to [33], the removal of forest cover in the context of timber management cannot be considered as deforestation, as it is expected that the forest will recover the original level of coverage, either naturally or with the help of silviculture techniques. At a national level, the National Institute for Space Research (INPE), officially responsible for the monitoring and accountability of Brazilian Amazon deforestation, considers deforestation as an anthropogenic process that can last several months [34].
Regarding SAR response to the deforestation process, some studies [35][36][37][38] point out that there is an ambiguity in the sign and intensity of the SAR signal change immediately after forest degradation. This change would be positive or negative, depending on the local management practices, the type of SAR sensor, and the span of time between measurements. Following on from this, Watanabe et al. [38] report an increase of 1.2 dB in ALOS 2-PALSAR 2 HH backscattering values in the early stages of deforestation, followed by a decrease of 1.2 dB in HV backscattering values after the fallen trunks were removed from the deforested patch. Meanwhile, the results of a recent simulation [39] indicate a decrease in backscattering after clear-cutting or selective logging in temperate forests. A complete Remote Sens. 2020, 12, 3922 4 of 31 field study in Peru's Madre de Dios area [40] reported a significant decrease in L-Band HV backscatter in forest patches subject to logging and/or deforestation. Additionally, a clear decrease in both L-band and C-band cross-polarized backscattering values was found in Bolivian dry forests subject to deforestation [29]. This decrease would be short-termed in the case of the C-band measurements. More recently, Bouvet et al. [41] found a short-term (~5-6 months) decrease in deforested C-band VV values, similar to the patterns reported in an Indonesian deforested tropical forest [28]. Thus, we can conclude that, for C and L-bands, there is an overall decrease in SAR backscattering measurements following a deforestation event, which can be quite short-termed (4 to 6 months), in the case of the C-band time series.
More details about backscattering changes in deforestation can be found in Table 2, which summarizes referenced measurements of backscattering change in deforested tropical forests. Automated or semi-automated methods for detection of deforestation can be categorized in two general kinds: (1) those which perform supervised classification on different dates to find differences among the computed maps (classification approach) and (2) those which compare the original values of the images on different dates (with or without previous segmentation), and then define changes as a function of the magnitude of the difference among images. We will call this second approach SAR-CD (SAR Change Detection).
Regarding the classification approach, regional classification of detailed forest classes including degraded or logged areas using SAR backscattering information is still challenging. While good results discriminating high-level forest classes in temperate forest exist in the literature [43], regional studies over tropical forests have found little contrast among forest classes while using incoherent attributes of SAR times series of ALOS-PALSAR [40,44], JERS-1 [35,45], SIR-C [46], or Sentinel-1 [47,48] satellite data. In [39], a detailed model found C and L-band backscattering to be responsive to defoliation and clear-cuts, but hardly useful to diagnose changes in forest structure or selective logging.
In this context, polarimetric (PolSAR) and interferometric (InSAR) attributes have been used as an attempt to enrich backscattering information to optimize forest cover classification. Incoherent polarimetric target decomposition techniques for forest classification were thoughtfully examined in [49] using C-band full polarization data over tropical dry forest stands in central India.
Yamaguchi decomposition [50] demonstrated a good performance, overcoming optical methods, but a significant confusion between different forest classes was still reported. Cloude-Pottier H/A/α decomposition [51] had poor results, mixing open, degraded, and dense forests. Those results are consistent with those presented in [52] which, although found an improvement in full-pol ALOS/PALSAR Remote Sens. 2020, 12, 3922 5 of 31 classification using polarimetric decomposition, especially Cloude-Pottier, still reported very low accuracies (around 18%) while classifying intermediate stages of forest regeneration.
Interferometric methods are typically used to generate high-resolution topographic maps [18] or to determine surface deformation [53]. Meanwhile, in the mid-90s interferometric coherence was found to be a useful parameter to measure boreal forest characteristics [54,55]. Further research proved that interferometric coherence may be used over tropical forests to increase classification accuracy and deforestation detection [56][57][58][59]. Experimental, tower-based results point out that time series of P and L-bands may be suitable for coherence use for classification and change detection, while C-band interimage coherence will degrade quicker, being unusable after 5-10 days [60,61]. Some recent results, however, relate good results using C-band Sentinel-1 interferometric coherence [62,63].
Regarding backscattering SAR-CD techniques, early attempts used manual thresholding to detect changes in forest cover [64]. Expectation-Maximization techniques [65] have been successfully applied in automatic thresholding [66]. A robust, generalized version of the Kittler-Illingworth thresholding method [67] to classify non-filtered and filtered ratio images of C and X-band data led to good results in filtered data on an agricultural-forest landscape. Generalized Gaussian (GG) models for the change and no-change classes in the log-ratio image, following a moment-based estimation method, have been used with good results [68]. An operational and reliable method based on Bayesian probability applied to deseasonalized backscattering time-series analysis is proposed in [69]. Applications can be reviewed in [28,29], which report excellent values of user accuracy (96.8%) and producer accuracy (95.8%) using Sentinel-1 VH data, although non-forest classes characterization is somewhat unclear. These allegedly constitute today's most accurate methods and will form the core of RADD (Radar Alerts for Detecting Deforestation), a forthcoming EWS method based on SAR data [70]. Other SAR-CD amplitude approaches, such as the one recently proposed in [71], are been used to complement the Bayesian approach. In the approach in [71], VV and VH change information to form a Polarimetric Change Vector (PCV), which can be helpful to discriminate the type of detected changes. This method seems promising to identify degradation types such as those caused by fire and flooding.
Regarding polarimetric SAR-CD methods, PolSAR data can overcome the lack of intensity data sensitivity to changes, using Wishart distribution to model L-band full-polarimetric data and to detect changes in forest cover [72]. Wishart distributions became a standard way to model PolSAR data, and many studies still use it to detect changes [73]. Muro et al. [74] applied the polarimetric omnibus technique based on Wishart distribution modeling, to dual-pol C-band time series over European wetlands, obtaining quite robust results on qualitative change detection [75]. Entropy and Alpha angle (H-α) polarimetric decomposition methods, originally developed for full polarimetric data, have been adapted to dual-pol data [76]. Despite certain controversy [77], some recent publications used this method in SAR-CD using Sentinel-1 and ALOS-PALSAR dual-pol data, reporting promising results using Entropy [78] and Alpha angle [79]. Concerning interferometric techniques, a fast, maximum-likelihood probabilistic approach to interferometric Coherence Change Detection (CCD) has been proposed [80]. Change detection in tropical forest cover using interferometric coherent decay constant constitutes a promising approach [63,81]; however, there is no account for accuracy or detection rates yet. Classic machine learning methods, such as neural networks, have been used in SAR-CD as early as 1991 [82], but its use was of lesser importance until the emergence of deep learning (DL) methods in 2013 [83]. A quick search in the Web of Science service informs us that today DL/ML techniques represent around 40% of the research on SAR-CD (author work, non-published). Convolutional Neural Networks (CNN) are the most used DL frameworks. Experimental results [84][85][86] seem to overcome conventional SAR-CD techniques. Other types of deep learning architectures, such as contractive autoencoders (sCAE) [87,88] or Generative Adversarial Networks (GANs) [89], are being applied with success in experimental SAR-CD, but still no operational algorithms have been developed to date [90].

SAR Data
The data used in this research correspond to the Sentinel-1A (S1A) IW GRD dataset, made available by the European Space Agency (ESA) and pre-processed in the Google Earth Engine (GEE) platform [91]. The main characteristics of the input images are presented in Table 3. Table 3. Main characteristics of the SAR data during the study.
Orbit file correction: Updates orbit metadata with a restituted orbit file.

2.
GRD border noise removal: Removes low-intensity noise and invalid data on scene edges.

3.
Thermal noise removal: Removes additive noise in sub-swaths to help reduce discontinuities between sub-swaths for scenes in multi-swath acquisition modes.

4.
Radiometric calibration: Computes backscatter intensity using sensor calibration parameters in the GRD metadata.

5.
Terrain correction (orthorectification): Converts data from ground range geometry, which does not take terrain into account, to normalized backscatter coefficient using the SRTM 30 m DEM.
Our area of interest (AOI) covers the entire Amazon biome within the Brazilian borders. The S1A data collection over the AOI stabilized between October and November 2016. Thus, the time interval that can be potentially used ranges between November 2016 and December 2019, totaling about 3 years of continuous acquisition, and approximately 30 yr −1 samples per location.
Most of the AOI can be considered as being flat or, at most, hilly, with few mountainous areas. A simple approach to angle correction was applied, where the original backscattering was corrected following the expression [93]: where σ 0 represents the normalized backscatter coefficient, θ i the local incidence angle (LIA), and γ 0 , known as gamma naught, is the backscattering coefficient normalized by the incidence angle. LIA was obtained using the procedure discussed in [94]. Mountainous areas should follow a more complex procedure, such as the one recently proposed in [95].

Sampling Spaces
To test the proposed stabilization and detection techniques, we randomly sampled Sentinel-1A measurements in two different, basin-wide sampling spaces. The first sampling space corresponds to invariant forest regions. To define these areas with emerged, preserved primary forests, we excluded the forest areas affected by clear-cut deforestation and/or degradation from the Amazonian Biome region, as per the INPE up-to-date data available in the Terrabrasilis [96] website (http://terrabrasilis. Remote Sens. 2020, 12, 3922 7 of 31 dpi.inpe.br/). The resulting areas were then eroded using a 5000 m buffer to retain forest kernel areas. Finally, we retained all the pixels that had an NDVI median value higher than 0.85 during 2019, as measured in cloud-free observation of the Landsat-8 OLI sensor. As a result of those operations, performed using a GEE script, a total area of 18,630 Mha was delimited as the invariant forest sample region. The second sampling space corresponds to the deforested area sampling region-i.e., areas that have suffered a deforestation event. To define this space, we used the Mapbiomas Alerta confirmed and refined 2019 warning polygons (MBA2019). The MBA2019 collection (http://alerta.mapbiomas.org/) is the result of a sophisticated semi-automated processing chain that uses very high resolution daily optical images to confirm and refine warning polygons issued from different alert systems. Every MBA2019 warning polygon has two attributes, "before_dt" and an "after_dt", which gives the date of the last cloud-free optical image before the deforestation and the first cloud-free optical image after the deforestation. An example of a selected MBA2019 polygon is shown in Figure 1. For this research, we built a subset from the MBA2019 collection, retaining the polygons issued from the DETER EWS system [97] larger than 1 ha, which had a validation time (lapse between the "before deforestation" and "after deforestation" optical image) lower than 4 months and a 3-month median NDVI greater than 0.82 before deforestation and lower than 0.69 after deforestation, as measured in a cloud-free observation of Landsat-8 OLI sensor. This selection is intended to avoid the effects of vegetation regrowth on the S1 signal. After selecting, a total of 2658 polygons were retained for analysis, encompassing all the Brazilian Amazon regions.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 33 (http://terrabrasilis.dpi.inpe.br/). The resulting areas were then eroded using a 5000 m buffer to retain forest kernel areas. Finally, we retained all the pixels that had an NDVI median value higher than 0.85 during 2019, as measured in cloud-free observation of the Landsat-8 OLI sensor. As a result of those operations, performed using a GEE script, a total area of 18,630 Mha was delimited as the invariant forest sample region. The second sampling space corresponds to the deforested area sampling region-i.e., areas that have suffered a deforestation event. To define this space, we used the Mapbiomas Alerta confirmed and refined 2019 warning polygons (MBA2019). The MBA2019 collection (http://alerta.mapbiomas.org/) is the result of a sophisticated semi-automated processing chain that uses very high resolution daily optical images to confirm and refine warning polygons issued from different alert systems. Every MBA2019 warning polygon has two attributes, "before_dt" and an "after_dt", which gives the date of the last cloud-free optical image before the deforestation and the first cloud-free optical image after the deforestation. An example of a selected MBA2019 polygon is shown in Figure 1. For this research, we built a subset from the MBA2019 collection, retaining the polygons issued from the DETER EWS system [97] larger than 1 ha, which had a validation time (lapse between the "before deforestation" and "after deforestation" optical image) lower than 4 months and a 3-month median NDVI greater than 0.82 before deforestation and lower than 0.69 after deforestation, as measured in a cloud-free observation of Landsat-8 OLI sensor. This selection is intended to avoid the effects of vegetation regrowth on the S1 signal. After selecting, a total of 2658 polygons were retained for analysis, encompassing all the Brazilian Amazon regions. After defining the sampling spaces, a simple operation performed using the stratified sample routine of the GEE API allowed us to randomly sample 3000 locations in the invariant sample space and 3000 locations in the deforested area sample space. These locations are shown in Figure 2. After defining the sampling spaces, a simple operation performed using the stratified sample routine of the GEE API allowed us to randomly sample 3000 locations in the invariant sample space and 3000 locations in the deforested area sample space. These locations are shown in Figure 2.

Methods
The extraction and processing sequence followed for every sample location is illustrated in Figure 3. To handle the number of images and processing needed to extract, filter, and stabilize in a time-effective fashion, we used the cloud-computing processing platform Google Earth Engine (GEE), which allowed us to perform a complex global analysis of remote sensing data using a simple programming interface. In this case, we used the GEE Python API. To avoid the effect of the incidence angle in mountainous areas, we removed the observations that had local observation angles of more than 45° or less than 30°.

Time Series Stabilization
To stabilize SAR data time series, we used two different methods: First, we applied a simple pixel-wise no-trend harmonic fitting approach, like the one used in [28]. Additionally, we extended the spatial normalization technique used in [29] by improving its applicability to an operational

Methods
The extraction and processing sequence followed for every sample location is illustrated in Figure 3. To handle the number of images and processing needed to extract, filter, and stabilize in a time-effective fashion, we used the cloud-computing processing platform Google Earth Engine (GEE), which allowed us to perform a complex global analysis of remote sensing data using a simple programming interface. In this case, we used the GEE Python API.

Methods
The extraction and processing sequence followed for every sample location is illustrated in Figure 3. To handle the number of images and processing needed to extract, filter, and stabilize in a time-effective fashion, we used the cloud-computing processing platform Google Earth Engine (GEE), which allowed us to perform a complex global analysis of remote sensing data using a simple programming interface. In this case, we used the GEE Python API. To avoid the effect of the incidence angle in mountainous areas, we removed the observations that had local observation angles of more than 45° or less than 30°.

Time Series Stabilization
To stabilize SAR data time series, we used two different methods: First, we applied a simple pixel-wise no-trend harmonic fitting approach, like the one used in [28]. Additionally, we extended the spatial normalization technique used in [29] by improving its applicability to an operational To avoid the effect of the incidence angle in mountainous areas, we removed the observations that had local observation angles of more than 45 • or less than 30 • .

Time Series Stabilization
To stabilize SAR data time series, we used two different methods: First, we applied a simple pixel-wise no-trend harmonic fitting approach, like the one used in [28]. Additionally, we extended the spatial normalization technique used in [29] by improving its applicability to an operational Remote Sens. 2020, 12, 3922 9 of 31 environment. This proposed stabilization algorithm is built upon the premise that the variations of C-band backscattering over tropical dense forest areas are mainly due to variations in canopy moisture content, variations in the amount of intercepted rain, the occurrence of storm cells and variations in sensor calibration [26]. These variations can be slow pathed (alternation between dry and wet season) or quite fast (intense rain shower during the acquisition of the image data). Spatially, the variations have a global (calibration), regional (overall wetness of the canopy), or local (rain interception after a heavy shower) scale. Figure 4 shows an example of backscatter signal instabilities observed over an unaltered stretch of evergreen forest in Eastern Amazonia. In this example, one can see how the moisture content in the canopy directly influences the backscatter intensity in the long-term, following dry and wet season cycles, and in the short-term, due to heavy showers before the data acquisition.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 33 environment. This proposed stabilization algorithm is built upon the premise that the variations of C-band backscattering over tropical dense forest areas are mainly due to variations in canopy moisture content, variations in the amount of intercepted rain, the occurrence of storm cells and variations in sensor calibration [26]. These variations can be slow pathed (alternation between dry and wet season) or quite fast (intense rain shower during the acquisition of the image data). Spatially, the variations have a global (calibration), regional (overall wetness of the canopy), or local (rain interception after a heavy shower) scale. Figure 4 shows an example of backscatter signal instabilities observed over an unaltered stretch of evergreen forest in Eastern Amazonia. In this example, one can see how the moisture content in the canopy directly influences the backscatter intensity in the longterm, following dry and wet season cycles, and in the short-term, due to heavy showers before the data acquisition. In this regard, we may correct the variations of the SAR signal by computing a correction coefficient (CC) for every pixel of every image of a lengthy SAR image collection, taking the overall variations of the forested neighborhood of the pixel as a reference. Specifically, CC is computed as: where is the forest spatial mean of the pixel's neighborhood, and is the temporal mean of this same value, computed for every image of the time series.
"Forest spatial mean" denotes a spatial mean that solely considers the forest-covered areas of the pixel surroundings. The detailed steps in stabilization are as follows: 1. Masking of all the available images of the interest area, leaving only the forest pixels unmasked.
The forest mask is computed from previous knowledge of the deforestation history of the area and then applied to all the sensor images of the interest area. 2. For every pixel of each image, the mean forest backscattering value is computed as the forest spatial mean of a 5 km radius neighborhood. This radius value was fixed considering the general spacing of the deforestation patches of the colonization roads of the Brazilian Amazon (the wellknown "fishbones"). 3. For every image, the correction coefficient is computed as the ratio between the forest spatial mean to the temporal mean of the same forest mean computed along the entire time series. 4. The final, stabilized backscatter value is computed by dividing the actual backscattering value by the correction coefficient. In this regard, we may correct the variations of the SAR signal by computing a correction coefficient (CC) for every pixel of every image of a lengthy SAR image collection, taking the overall variations of the forested neighborhood of the pixel as a reference. Specifically, CC is computed as: where FS is the forest spatial mean of the pixel's neighborhood, and T FS is the temporal mean of this same value, computed for every image of the time series. "Forest spatial mean" denotes a spatial mean that solely considers the forest-covered areas of the pixel surroundings. The detailed steps in stabilization are as follows: 1.
Masking of all the available images of the interest area, leaving only the forest pixels unmasked. The forest mask is computed from previous knowledge of the deforestation history of the area and then applied to all the sensor images of the interest area.

2.
For every pixel of each image, the mean forest backscattering value is computed as the forest spatial mean of a 5 km radius neighborhood. This radius value was fixed considering the general spacing of the deforestation patches of the colonization roads of the Brazilian Amazon (the well-known "fishbones").

3.
For every image, the correction coefficient is computed as the ratio between the forest spatial mean to the temporal mean of the same forest mean computed along the entire time series.

4.
The final, stabilized backscatter value is computed by dividing the actual backscattering value by the correction coefficient.
This procedure allows the proposed stabilization to be applied to any area threatened by deforestation, with two caveats: (1) the user must have a reliable, up-to-date forest/non-forest map valid at the start of the detection period, and (2) the proposed procedure demands a great deal of computational resources. Its operational implementation should use a cloud-based solution, such as GEE.
To initially evaluate the effects of both stabilization algorithms, we computed the statistical parameters of the stabilized invariant forest time series and compared them with the original data. We used the standard deviation and coefficient of variation (CV) measures to perform this analysis. We then evaluated the stabilization algorithm using 1-and 2 year-long series.

Time Series Filtering
Speckle noise is a phenomenon common to all SAR measurement. In order to reduce speckle, we applied the spatio-temporal filtering procedure proposed in [98]. It has been applied with success over Sentinel-1A series [29,41,99], and it can be considered nowadays as a fundamental part of the processing chain of dense temporal SAR data. Following [21], the kth filtered image at position (x,y) is given by: where I k (x, y) is the backscattering value of image k at the position x,y and I k (x, y) is the speckle-free estimated value of the image k at the position x,y. Additionally, a spatial filter was applied to every temporally-filtered image, as recommended by [98]. Previous works [100] point to the Refined Lee filter as a good compromise between noise-removal performance, edge preservation, and operational usability, reaching a performance similar to advanced filters based on complex computations. Thus, the filtering step is composed by the Quegan and Yu temporal filtering with a 5 × 5 box-car internal filter and 7 × 7 Refined Lee as the spatial filter.

Deforestation Detection
After stabilization and filtering, two different deforestation detection techniques were used, corresponding to the two main approaches to deforestation detection: Maximum Likelihood Classification and Adaptive Thresholding.
The pixel-wise maximum likelihood classification (MLC) approach used a Bayesian classification framework to classify S1 samples as belonging to forest or deforested areas. This approach has been successfully applied over SAR time series [69], and is based on the classic Bayesian conditional probability equation. Applied to our context, it is expressed by the following equation: where P(ω i |x t ) is the posterior probability of the sample x t belong to the state ω i , P(ω i ) is the prior probability of the state ω i , and p(x t ω i ) is the likelihood of x t given the state ω i . The simplest Bayesian classification framework then builds a discriminant function g i (x t ) = P(ω i |x t ) for every possible state and then classify the state of x t based on the discriminant with the maximum value. Other monotonic discriminant functions can be used. Following [101], we used: Remote Sens. 2020, 12, 3922 11 of 31 In our case, we consider only two possible states: forest (F) and deforested (DF). Following (4) we can define a two-class discriminant function (called a dichotomizer) as: Substituting (4) in (6), we obtain the final expression of our discriminator: where p(x t |DF) is the likelihood of x f being deforested, given by the pdf of the deforestation state; p(x t |F) the likelihood of x f being forest, given by the probability density function (pdf) of the forest state; P(DF) the prior probability of the deforestation state; and P(F) the prior probability of the forest state. Finally, we classify x t as being a deforested location if g(x t ) > 0. As the ln P(DF)/P(F) constant (that we call the detection threshold) is initially unknown (although we know that generally, P(F) is much greater than P(DF)), we used the value that maximizes the accuracy of the detection results, following a heuristic approach. In this manner, we ran all the detection algorithms using a range of thresholds going from 0 (deforested and forested equally probable) to 10 (forest much more probable). A value of 3, for instance, means that the a priori probability of finding a forest pixel is e 3 ≈ 20 times higher than the probability of finding a pixel that was deforested during the detection period.
The maximum-likelihood classification of the F/DF time series requires the determination of the probability density function of the backscattering values, for the F and DF classes. While the determination of the parameters of the speckle-filtered forest series was trivial once we had a wealth of data points of invariant forest locations, the deforestation samples are somewhat trickier to study as we do not know exactly the time of the deforestation. Additionally, even if we know that backscattering values normally drop after a deforestation event, we do not know exactly the time that every deforested location would take to recover the previous backscattering values. Thus, to separate the forest and non-forest samples of the deforestation time series, we used a 2-component finite mixture statistical modeling technique [102], using a six-month sampling of the training TS, starting at some point after the deforestation event (this date is given by the after_dt parameter associated with every warning polygon), and assuming the normality of the subjacent distributions. Then, we studied the relationship between the modeled F and DF distribution parameters, looking for a simple and meaningful relationship that might be applied to other locations where the time or occurrence of deforestation is unknown.
The statistical model to be applied to the F/DF classis during the maximum-likelihood is defined by studying the statistical distribution of the extracted time series. We used the Shapiro-Wilk and Kolmogorov-Smirnoff tests to determine the statistical distribution models more suitable for the forest-invariant extracted samples, before and after the stabilization and filtering processes. This characterization was made on a pixel-wise basis, taking all the available samples at the forest locations, building the corresponding γ 0 time series (TS) and determining the Shapiro-Wilk and Kolmogorov-Smirnoff test p-values for normal and gamma distributions, respectively. We assumed the deforested samples to follow a similar statistical distribution model.
The second, more straight-forward detection approach uses an adaptive linear threshold (ALT) based on the relationship between the mean and the minimum value of the forest TS. The premise of this method is that the SAR signal variability over tropical forest tends to be constant, and thus, any sample that falls below an interval of backscattering values should be considered as being non-forest and flagged as a deforested sample. Following this premise, the forest threshold is defined as: where n is the number of collected forest locations, γ 0 Fi is the mean value of the forest TS of the ith location, and min γ 0 F i is the minimum value of the forest TS of the ith location. Following (8), we defined the forest threshold Th F using (a) a set of training samples, extracted from the invariant forest locations, and (b) a set of training samples extracted from the deforested locations samples predating the deforestation event (soon-to-be deforested series). Then, all the TS were thresholded using the computed Th F value for both cases. In order to optimize the detection performance, we tested the detection performance of a set of thresholds around the Th F value.

Detection Validation
As the final purpose of the proposed SAR data processing techniques to improve near real-time deforestation detection at the pixel-level, we applied a validation framework based on the statistical evaluation of the detection results at randomly sampled 20-m sized pixels. We decided to operate at the 20-m pixel, as this is the resolution of the original S1A images in range-azimuth geometry after multi-looking.
Once defined, the sampling spaces for invariant forest and deforested classes, and the images processed following the steps detailed until here, we proceeded to extract data series at random locations using the GEE python API.
For every random location, we extracted the time series corresponding to: Both VH and VV polarizations were used and evaluated during this phase. We also used two different time TS-lengths: 1 and 2 years. For every sampled location, we performed a complete detection processing using the two methods described above, looking for deforestation that occurred in a 4-month interval after the "after deforestation" date, in the case of the deforested locations, and after a random date, in the case of the unaltered forest locations. If deforestation is detected, the day of the first deforestation warning will be recorded, and the sample will be classified as "deforested". If no deforestation is detected, the sample will be classified as "no deforested". This procedure will be repeated for every sample, using different detection thresholds.
Once the detection process was completed, we computed the confusion matrix for every combination of Polarization, Treatment, and Threshold. The confusion matrix elements were combined to compute global accuracy (ACC), Specificity (TNR), and Sensitivity (TPR). They are given by: Speci f ity (TNR) = TN TN + FP (10) where: The Specificity, also called Selectivity or true negative rate (TNR), is the most important parameter to consider in an EWS, as it measures the ability of the algorithm to avoid false warnings (commission errors). The Sensitivity, also called Recall or true positive rate (TPR), complements the Specificity by informing the completeness of the algorithm detection (the lack of omission errors). The accuracy parameter is used as a measure of the global performance of the algorithms, but we use it with precaution, as it tends to give the same weight to omissions and commission errors, being our operational context much more sensitive to commission errors. Therefore, the different methods and treatments proposed here were compared using Accuracy, but also their relationship between Specificity and Sensitivity, by using Receiver-Operator Characteristic (ROC) curves and by comparing the TPR value for a given, optimal, TNR.

Code Availability
The above-described acquisition, processing, and validation of the Sentinel-1 data were performed in three different environments: • Earth Engine Javascript Code Editor: Definition of the sampling locations. All the used scripts and the resulting data are available in https://github.com/jdoblas/DETER_ SAR_RS.

Results
Following the proposed methodology, we extracted 24 different time-series of incidence angle-corrected intensity S1A backscattering (which we name as γ 0 ), in every one of the 6000 randomly selected locations. The number of series per location corresponds to the combinations of two different collecting periods (1 and 2 years) in two different polarizations (VH and VV), undergoing 6 different types of treatment, labeled as follows: • orig-original Sentinel-1 backscattering samples. • origf -original TS filtered. • harmon-original TS stabilized using harmonic fitting. • harmonf -original TS stabilized using harmonic fitting and then filtered. • spatial-original TS stabilized using spatial stabilization. • spatialf -original TS stabilized using spatial stabilization and then filtered Figure 5 shows two examples of each set of time series. As some of the randomly selected locations did not have an adequate Sentinel-1A covering, the initial extracting procedures resulted in 2978 deforested locations and 2983 invariant forest locations with complete (n ≥ 20) TS, rather than the expected 3000 locations for every class. The missing locations are located near the Brazilian southern border, where some marginal regions are being collected by the Sentinel-1B satellite, without any covering of the Sentinel-1A. At the initial stages of this research, we observed a substantial difference in the backscattering patterns of the satellites, and thus we decided to focus the remainder of this work on the S1A measurements.
After the Local Incidence Angle filtering procedure, we retained 2565 deforested locations and 2480 invariant forest locations, which represents 86.13% and 83.13% of the original deforested and forested points, respectively. The total number of collected SAR measurements is 8,726,898, divided into different polarizations, data collecting period, and treatments.
• origf-original TS filtered. • harmon-original TS stabilized using harmonic fitting. • harmonf-original TS stabilized using harmonic fitting and then filtered. • spatial-original TS stabilized using spatial stabilization. • spatialf-original TS stabilized using spatial stabilization and then filtered Figure 5 shows two examples of each set of time series. As some of the randomly selected locations did not have an adequate Sentinel-1A covering, the initial extracting procedures resulted in 2978 deforested locations and 2983 invariant forest locations The main characteristics of the times series are summarized in Tables S1 and S2 (Supplementary Materials), for the forest and deforested locations, and for every treatment and polarization.
The mean value of the forest γ 0 was determined as being −12.4 and −6.25 dB for VH and VV polarizations, respectively. The unfiltered series have a substantive variability in both polarizations of 1.75 dB (1σ), which is reduced to~0.57 dB when applied to the proposed speckle filtering. The mean backscattering values of the deforested series are consistently lower than the forest series, reaching −13.1 and −7.6 dB for VH and VV polarizations.
Regarding stabilization results, all the tested stability measures showed a statistically significant difference between the non-stabilized series and the stabilized ones. The stabilized TS show a diminution of the variability, as measured by the CV, of 3%-5% (for harmonic stabilization) and~10% for the spatial stabilization (Table S3 and Figure S1, Supplementary Materials). Other measures of variation, such as mean absolute deviation (MAD), standard deviation (SD), and interquantile ratio (IQR) yielded very similar results.
Before enterprise pixel-wise modeling, we tested the forest time-series for normality. The results of the normality Shapiro-Wilk testing showed that~80% of the filtered TS are normal (p = 0.05) (see Figure S2 and Table S4, Supplementary Materials). The non-filtered series were successfully fitted using a Gamma distribution (see Table S5, Supplementary Materials).

Detection of Deforestation
As detailed previously, the deforestation detection procedure followed two different approaches. The first approach, based on maximum-likelihood classification, uses the deforestation TS to build a 2-component mixture model of normal distributions, which should represent the statistical properties of the forested (F) and deforested (DF) samples ( Figure 6).

Detection of Deforestation
As detailed previously, the deforestation detection procedure followed two different approaches. The first approach, based on maximum-likelihood classification, uses the deforestation TS to build a 2-component mixture model of normal distributions, which should represent the statistical properties of the forested (F) and deforested (DF) samples ( Figure 6). To verify the separability of the obtained distributions, we computed the Jeffries-Matsushita (JM) distance [104] between the F and DF components, for every TS. All treatments presented a good separability with mean JM distances greater than 1.19 (see Figure S3, Supplementary Materials).
After deducing the parameters of the statistical distribution of the F/DF series in a set of training samples, we studied the relationship between them. As can be seen in Figure 7, the means of the deforested class have a good linear relationship with the means of the forest class, with R 2 values ranging from 0.59 to 0.64.
Meanwhile, the variance values of the deforested distributions did not correlate with the values of the variance, or the mean, of the forested samples, with R 2 values below 0.2. Thus, we used the mean value of the variance in every treatment to characterize deforestation variance.
Thus, the statistical distribution of the deforested class was modeled, in every location, as a normal distribution with a mean linearly related to the (known) mean of the forest class, and a constant variance, which was deducted from the training dataset.
Once we established the statistical parameters of the F/DF distributions, we were ready to start the maximum-likelihood detection process using the maximum likelihood approach. Initially, we computed the values of p(F) and p(DF) of every sample of the validation dataset during a four-month detection period. Then, the likelihood of the sample belonging to the deforested class was computed by dividing p(DF) by p(F) ( Figure S4, Supplementary Materials).  The final results of the detection procedure using the Maximum Likelihood approach are represented in Table 4, Figures 8 and 9.   We also considered the best results given an optimum level of TNR. Based on field experience, we fixed this level as 99.5%, which is equivalent to a rate of false positives of 0.5%. Table 5 shows the results of the MLC detection procedure given this level of TNR.  The second deforestation detection technique used-i.e., Adaptive linear thresholding (ALT)is based on the study of the variability of the values of the forest time series. For every forest TS, we computed the difference between its mean and its p1 (1% percentile) value. We also computed this relationship for the deforested values, selecting the samples measured before the deforestation event (as specified by the corresponding warning polygon metadata) in a 50%-split training subset.
The results of this analysis show that the distance between the forest mean and the minimum value of the forest TS is nearly constant (see Figure 10). Similarly, we may use this finding as a We also considered the best results given an optimum level of TNR. Based on field experience, we fixed this level as 99.5%, which is equivalent to a rate of false positives of 0.5%. Table 5 shows the results of the MLC detection procedure given this level of TNR.   The second deforestation detection technique used-i.e., Adaptive linear thresholding (ALT)-is based on the study of the variability of the values of the forest time series. For every forest TS, we computed the difference between its mean and its p1 (1% percentile) value. We also computed this relationship for the deforested values, selecting the samples measured before the deforestation event (as specified by the corresponding warning polygon metadata) in a 50%-split training subset.
The results of this analysis show that the distance between the forest mean and the minimum value of the forest TS is nearly constant (see Figure 10). Similarly, we may use this finding as a straight-forward approach to detect deforestation by using a range of values around the distance to threshold and flagging deforestation-suspected values of the time-series (see how this can be applied to some TS in Figure 11).
After learning with the training dataset about the relationship between the forests' backscattering means and minima, we applied this relationship to the validation dataset, flagging as deforested the samples that fell below the expected minimum, following a simple procedure to find the optimum thresholding level:

1.
Forest_mean = mean value of the validation TS samples before the deforestation date.

2.
Treatment_distance_mean = mean of the distances between the mean and the p1 value of the training dataset TS, for every treatment.

3.
Treatment_distance_sd = standard deviation of the distances between the mean and the p1 value of the training dataset TS, for every treatment.
Flag TS samples < Threshold.
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 33 straight-forward approach to detect deforestation by using a range of values around the distance to threshold and flagging deforestation-suspected values of the time-series (see how this can be applied to some TS in Figure 11). Figure 10. Relationship between the mean value of the forest VH backscattering and the minimum (p1) value of the time series in forest pixels (A). As an alternative analysis, chart (B) shows the relationship between the mean and the difference between the mean and the p1 value. After learning with the training dataset about the relationship between the forests' backscattering means and minima, we applied this relationship to the validation dataset, flagging as deforested the samples that fell below the expected minimum, following a simple procedure to find the optimum thresholding level: 1. Forest_mean = mean value of the validation TS samples before the deforestation date. 2. Treatment_distance_mean = mean of the distances between the mean and the p1 value of the training dataset TS, for every treatment. 3. Treatment_distance_sd = standard deviation of the distances between the mean and the p1 value of the training dataset TS, for every treatment. 4. Central threshold = Forest_mean-Treatment_distance_mean. 5. For thresholding_factor = −5 to 5: a. Threshold = Central_threshold-(Treatment_distance_sd*thresholding_factor). b. Flag TS samples < Threshold.
We ran this procedure for all the training datasets' TS and computed the TS that were correctly flagged and those that were not. Then, we computed the confusion matrix and accuracy measures (see Table 6, Figures 12 and 13). We ran this procedure for all the training datasets' TS and computed the TS that were correctly flagged and those that were not. Then, we computed the confusion matrix and accuracy measures (see Table 6, Figures 12 and 13).    We also considered the best results of the ALT method given an optimum level of TNR. Based on field experience, we fixed this level as 99.5%, which is equivalent to a rate of false positives of 0.5% (Table 7). We also considered the best results of the ALT method given an optimum level of TNR. Based on field experience, we fixed this level as 99.5%, which is equivalent to a rate of false positives of 0.5% (Table 7). The results of the validation phase showed that both methods yielded satisfactory results. However, the ALT method performed better than the MLC method, especially regarding the TPR measure of accuracy in a fixed false alarm rate context. This index reached 89.61% using ALT methodology, overcoming MLC's best result (78.92%).

Algorithm Scalation to Support an EWS
The preliminary results of the location-based algorithms led us to a further step, where we developed a set of algorithms to apply the proposed methodologies to arbitrary areas using the google earth engine (GEE) platform. Thus, we built a three-step procedure to create warning polygons associated with the backscattering anomalies, as detected by the ALT and MLC methods, in user-defined areas. The steps taken were: (1) processing of the image collections available for the area and period specified by the user (including a 2-year training period); (2) a pixel-based detection routine, which produces a multi-band image with the detection thresholding; and (3) vectorization of the warning raster to create a set of warning polygons.
This operational implementation is detailed in Figure 14.
A fully-operational, completely automated EWS based on the proposed methodologies was coded in the GEE Python API in the context of INPE's DETER project. The so-called DETERSAR system has been up and running since 21 June 2020, and it has already issued 25,107 warning alerts over five critical areas in the Brazilian Amazon Basin (as of 22 September 2020).
During the development of the EWS version of the algorithms proposed here, we were forced to reevaluate the algorithm's parametrization to fit the needs of the final users of the EWS. In this case, the user of the system will be willing to sacrifice accuracy in exchange for a very low false-positive rate (FPR). As a matter of fact, a false deforestation warning might misguide ground enforcement teams, and provoke an unnecessary waste of time and resources. This low-FPR aspect linked to the EWS is emphasized by [105].
To illustrate this point, we consider the case of the ALT method proposed here, in which the global accuracy optimum threshold (1.5) produces several false positives, linked to its FPR (1.9%). Even if this number is very low, a EWS user might want to choose a stricter threshold, in order to avoid those false positives. Visual inspection ( Figure 15) suggests a threshold value of 2.5. By shifting the threshold in this way, we lost global accuracy-from 95.9% to 94.4%-but we lowered the FPR from 1.9% to 0.5%, which dramatically reduced the number of false warnings of a potential SAR-based EWS system. google earth engine (GEE) platform. Thus, we built a three-step procedure to create warning polygons associated with the backscattering anomalies, as detected by the ALT and MLC methods, in user-defined areas. The steps taken were: (1) processing of the image collections available for the area and period specified by the user (including a 2-year training period); (2) a pixel-based detection routine, which produces a multi-band image with the detection thresholding; and (3) vectorization of the warning raster to create a set of warning polygons.
This operational implementation is detailed in Figure 14.  A fully-operational, completely automated EWS based on the proposed methodologies was coded in the GEE Python API in the context of INPE's DETER project. The so-called DETERSAR system has been up and running since 21 June 2020, and it has already issued 25,107 warning alerts over five critical areas in the Brazilian Amazon Basin (as of 22 September 2020).
During the development of the EWS version of the algorithms proposed here, we were forced to reevaluate the algorithm's parametrization to fit the needs of the final users of the EWS. In this case, the user of the system will be willing to sacrifice accuracy in exchange for a very low falsepositive rate (FPR). As a matter of fact, a false deforestation warning might misguide ground enforcement teams, and provoke an unnecessary waste of time and resources. This low-FPR aspect linked to the EWS is emphasized by [105].
To illustrate this point, we consider the case of the ALT method proposed here, in which the global accuracy optimum threshold (1.5) produces several false positives, linked to its FPR (1.9%). Even if this number is very low, a EWS user might want to choose a stricter threshold, in order to avoid those false positives. Visual inspection ( Figure 15) suggests a threshold value of 2.5. By shifting the threshold in this way, we lost global accuracy-from 95.9% to 94.4%-but we lowered the FPR from 1.9% to 0.5%, which dramatically reduced the number of false warnings of a potential SARbased EWS system.

Discussion
Initially, the vast number of time series collected for the present work allowed the examination of the sign of backscattering change related to deforestation events. The mean backscattering values of the deforested series are consistently lower than the forest ones, reaching −13.1 and −7.6 dB for VH and VV polarizations, respectively, against −12.4 and −6.25 dB for invariant forest areas. Thus, we can conclude that deforestation causes a decrease in C-band SAR backscattering. Meanwhile, a closer look at the extracted time series shows a more complex reality. The typical temporal profile of a clear-cut area in the Amazon tropical forest found during this work encompasses different stages, which we can tentatively summarize as:

1.
Forest canopy with little or no damage.

2.
Initial deforestation, where some tree trunks and standing trees remain. 3.
Fire removes vegetation and promotes double bounce (remaining standing trees with bare soil) returns to the sensor.

4.
Pasture grows with increasing precipitation and remnants of burned vegetation 5.
Pasture reaches a decimetric height and becomes hardly distinguishable from the original forest cover for SAR-C incoherent signal.
Up to the submission of this article, no research has been published focusing on those short-term patterns, but based the patterns of the extracted time series, we consider that clear-cut deforestation will cause a~6-month cycle of low backscattering in the C-SAR measurements as a working hypothesis. The length of this period will depend on a variety of factors, such as the type of deforestation, the season of the year, and the use of the deforested area.
Regarding variability, the extracted series variability (σ ≈ 1.75 dB in both polarizations) is highly consistent with the models proposed by Benninga et al. [26], which predict a radiometric uncertainty of 2.02 dB for a 20 × 20 m measurement of VH and 1.77 dB for VV backscattering measures.
The results of the stabilization procedures indicate that the stabilized series are more stable than the original ones, the spatial stabilization being more performant than the harmonic stabilization.
The detection algorithm results, though very satisfactory, did not register substantial differences between the un-stabilized series and the stabilized ones. It is true that spatial stabilization offers the best results in MLC workflow, but the difference in accuracy (94.36% to 93.93%) was not as prominent as expected. The advantages of spatial stabilization are more evident when using higher thresholds, such as those used in operational EWS. In this case, the differences in accuracy of the MLC method can improve from 84.16% to 89.05% using spatial stabilization.
Visual inspection of the stabilized time series (available on the online repository) reveals that some of the contrast and spikes related to deforestation are reduced by the stabilization processes. This is especially true for the dry season, where the forest canopy reduces its moisture content, and thus its backscattering intensity. Stabilization algorithms then tend to raise values of the TS, diminishing the contrast of eventual deforestation events. Our work proves that the loss of accuracy linked to this effect can overcome the benefits of stabilization.
These results can be seen as complementary to those presented by Reiche et al. [28] which first applied harmonic fitting to stabilize SAR TS in Bolivian dry forests, which are known to have a strong seasonality. Our results show that this correction may be unnecessary for evergreen forests. We still think that stabilization can be of great value in transitional forests, such as those of the Amazonia-Cerrado ecotone.
False deforestation alerts due to high-density storm clouds, relatively frequent in the change-ofseason months, are still a problem to be resolved. Several anomalies linked to this phenomenon were detected in time series while verifying false-alerts events, and in actual SAR images ( Figure 16). This kind of anomalies can disrupt the operation of a fully-automated EWS, as they can create vast deforestation polygons in untouched forest areas. Future attempts to control and stabilize Figure 16. Anomalies in Sentinel-1 backscattering measurements due to heavy rainstorms can produce false deforestation alerts. In the figure, a deforestation warning, measuring 359 ha. was related to a heavy-rain episode registered by Global Precipitation Measurement (GPM) records (left panel, the vertical dashed line shows the date of the image). Other peaks in rain profiles did not provoke anomalies. This is due to the difference in the resolution of the SAR data (20 m) and the GPM data (~5 km).
As an overall analysis of the detection results, it can be said that the obtained global accuracy and the false positive and negative rates were high, given the relative simplicity of the proposed methods and the continental extension of the experiment. An EWS user interested in reducing the false positive rate (FPR) to a minimum (e.g., 0.5%) will still have a true positive rate higher than 79%, in the case of the MLC method, and 90%, in the case of the ALT method. That means that using these principles we can build a very reliable EWS system with a low rate of omissions. The duration of the training period seems to improve the ALT method results. A 3-year period should slightly improve these results.
Another key aspect of the proposed TS-based pixel-wise detection methodologies is the ability to follow the dynamics of deforestation, as for every deforestation warning, it is possible to see the day of deforestation of every pixel ( Figure 17).

Conclusions
In this work, we investigated the performance of two different stabilization methods using two different deforestation detection algorithms in an evergreen tropical forest context. The results indicate that spatial stabilization is useful and can improve detection results if used within a maximum-likelihood detection approach. Meanwhile, stabilization did not seem to improve the detection results while using an Adaptive Linear Thresholding approach. Other peaks in rain profiles did not provoke anomalies. This is due to the difference in the resolution of the SAR data (20 m) and the GPM data (~5 km).
As an overall analysis of the detection results, it can be said that the obtained global accuracy and the false positive and negative rates were high, given the relative simplicity of the proposed methods and the continental extension of the experiment. An EWS user interested in reducing the false positive rate (FPR) to a minimum (e.g., 0.5%) will still have a true positive rate higher than 79%, in the case of the MLC method, and 90%, in the case of the ALT method. That means that using these principles we can build a very reliable EWS system with a low rate of omissions. The duration of the training period seems to improve the ALT method results. A 3-year period should slightly improve these results.
Another key aspect of the proposed TS-based pixel-wise detection methodologies is the ability to follow the dynamics of deforestation, as for every deforestation warning, it is possible to see the day of deforestation of every pixel ( Figure 17).
Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 33 Figure 16. Anomalies in Sentinel-1 backscattering measurements due to heavy rainstorms can produce false deforestation alerts. In the figure, a deforestation warning, measuring 359 ha. was related to a heavy-rain episode registered by Global Precipitation Measurement (GPM) records (left panel, the vertical dashed line shows the date of the image). Other peaks in rain profiles did not provoke anomalies. This is due to the difference in the resolution of the SAR data (20 m) and the GPM data (~5 km).
As an overall analysis of the detection results, it can be said that the obtained global accuracy and the false positive and negative rates were high, given the relative simplicity of the proposed methods and the continental extension of the experiment. An EWS user interested in reducing the false positive rate (FPR) to a minimum (e.g., 0.5%) will still have a true positive rate higher than 79%, in the case of the MLC method, and 90%, in the case of the ALT method. That means that using these principles we can build a very reliable EWS system with a low rate of omissions. The duration of the training period seems to improve the ALT method results. A 3-year period should slightly improve these results.
Another key aspect of the proposed TS-based pixel-wise detection methodologies is the ability to follow the dynamics of deforestation, as for every deforestation warning, it is possible to see the day of deforestation of every pixel ( Figure 17).

Conclusions
In this work, we investigated the performance of two different stabilization methods using two different deforestation detection algorithms in an evergreen tropical forest context. The results indicate that spatial stabilization is useful and can improve detection results if used within a maximum-likelihood detection approach. Meanwhile, stabilization did not seem to improve the detection results while using an Adaptive Linear Thresholding approach.

Conclusions
In this work, we investigated the performance of two different stabilization methods using two different deforestation detection algorithms in an evergreen tropical forest context. The results indicate that spatial stabilization is useful and can improve detection results if used within a maximum-likelihood detection approach. Meanwhile, stabilization did not seem to improve the detection results while using an Adaptive Linear Thresholding approach.
Both deforestation detection methods offered good results, and its implementation as the base of an operation Early Warning System seems very promising, especially if combined with the use of thresholds that can spatially change as a function of the risk of deforestation.
The validation of our results used the DETER polygons refined by the MapBiomas Alerta initiative. These polygons constitute an extraordinary source of training and validating samples, and should be used in future developments, as the initiative continues to generate warnings polygons.
Global application of the proposed methodologies should be straightforward, by adapting the provided GEE scripts. The used mask layers, which are specific to the Brazilian Amazon, can be substituted by global datasets such as Hansen's Global Forest Change data [106].
Further research should focus on the use of polarimetric and interferometric attributes of change detection, meteorological stabilization using high-resolution meteorological data, integration with other sensors, and detection in transitions to savanna regions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/23/3922/s1, Table S1: Characteristics of the collected forest Time Series, Table S2: Characteristics of the collected deforestation Time Series, Figure S1: Box-plot of the coefficient of variation of the forest speckle-filtered time series, Table S1: Mean reduction of the variability of the time series, Figure S2: Results of the Shapiro-Wilk normality test applied to the treated and non-treated forested TS, Table S4: Results of the Shapiro-Wilk test for normality of Forest time series, Table S5: Results of the Kolmogorov-Smirnoff test for gamma-distribution fitting of unfiltered Forest time series, Figure S3: Jeffries-Matsushita distance of the F/DF classes as computed by the finite mixture algorithm, Figure S4: Density distribution of the values of the deforestation log-likelihood of the validation samples.