Use of the SAR Shadowing Effect for Deforestation Detection with Sentinel-1 Time Series

To detect deforestation using Earth Observation (EO) data, widely used methods are based on the detection of temporal changes in the EO measurements within the deforested patches. In this paper, we introduce a new indicator of deforestation obtained from synthetic aperture radar (SAR) images, which relies on a geometric artifact that appears when deforestation happens, in the form of a shadow at the border of the deforested patch. The conditions for the appearance of these shadows are analyzed, as well as the methods that can be employed to exploit them to detect deforestation. The approach involves two steps: (1) detection of new shadows; (2) reconstruction of the deforested patch around the shadows. The launch of Sentinel-1 in 2014 has opened up opportunities for a potential exploitation of this approach in large-scale applications. A deforestation detection method based on this approach was tested in a 600,000 ha site in Peru. A detection rate of more than 95% is obtained for samples larger than 0.4 ha, and the method was found to perform better than the optical-based UMD-GLAD Forest Alert dataset both in terms of spatial and temporal detection. Further work needed to exploit this approach at operational levels is discussed.


Introduction
At present, forest destruction is a major source of carbon emissions and a primary cause of biodiversity loss [1].Globally, deforestation and forest degradation are estimated to account for approximately 20% of anthropogenic CO 2 emissions [2], though the rate of net forest loss has been reported to have halved from 7.3 Mha•year −1 in the 1990s to 3.3 Mha•year −1 between 2010 and 2015 [3].However, estimates of the carbon emissions that result from forest disturbances still contain considerable uncertainty [4].The call to reduce uncertainties in estimating changes in forest cover is also driven by the reporting needs outlined in the Reducing Emissions from Deforestation and forest Degradation (REDD+) program.
Satellite imagery is the primary tool for providing information on newly deforested areas in vast and sometimes inaccessible forests [5], with most monitoring approaches relying predominantly on optical remote sensing.In particular, stimulated by the opening of the Landsat archive in combination with the ability to download fully pre-processed images, efforts in recent years shifted towards operational and large-scale deforestation monitoring systems based on Landsat time series, at annual scales [6,7] or even with near-real time (NRT) capabilities [8].Hansen et al. (2016) [8] demonstrated the potential for and constraints of operational Landsat based deforestation alerts for the humid tropics.A major limitation for optical-based NRT applications is the presence of haze in the dry season (caused by fire) and, more importantly, of clouds in the wet season [7][8][9][10].Some regions suffer from pervasive cloud cover throughout the entire year, and even more than one year [8,10,11].Incorporating Sentinel-2 (S2) data would increase data richness and improve the detection of change.Using synthetic aperture radar (SAR) is another option.
As a matter of fact, SAR is one of the most promising remote sensing tools for the NRT mapping of forest disturbances in the tropics, thanks to its ability to operate in all weather conditions at any time of day or night.However, the use of SAR data in forest disturbances mapping has not been well developed, so far, compared to optical data, and thus, its operational application has not yet been realized.This can be explained by the fact that many users are less able to process, interpret, and analyze SAR data, whereas optical data benefit from their relative ease of processing and interpretation, and the continuity of medium-resolution (10-30 m) observations since the 1970s [12].Fragmented and inconsistent data archives are a problem common to many SAR missions, as data have traditionally been collected for local or regional studies with no systematic acquisition planning, resulting in scarce spatial and/or temporal coverages.The L-band systems ALOS PALSAR (2006)(2007)(2008)(2009)(2010)(2011) and ALOS-2 PALSAR-2 (since 2014) launched by the Japan Aerospace Exploration Agency (JAXA) represent a significant exception, as their systematic acquisition strategies have been designed in such a way that global mosaics can be produced annually, resolving the spatial coverage issue.The release of the global JERS, PALSAR, and PALSAR-2 mosaics at 25 m resolution has fostered a large number of new studies for forest monitoring.For example, Shimada et al. (2014) [13] have used PALSAR mosaics to produce the first SAR-based annual (2007)(2008)(2009)(2010) global maps of forest and non-forest cover, from which some maps of forest losses and gain were generated based on thresholds.To the best of our knowledge, the largest area, i.e., at the subcontinental scale, where forest disturbances and regrowth have been assessed using SAR intensity changes, was achieved by Mermoz and Le Toan (2016) [14].Ryan et al. (2012) [15] has adopted a different approach in a small area in Central Mozambique, which consists in first estimating forest aboveground biomass (AGB) at two periods, and then estimating changes in AGB by subtracting the two estimates.This method is relevant in the context of REDD+ measurement, reporting, and verification (MRV), because both disturbance areas and biomass losses are estimated, but suffers from error propagation as errors in both maps may be summed.All these approaches based on ALOS PALSAR are, however, constrained by the small number of available observations (1 observation per year in the case of the mosaics, and 1 observation every 42 days at best with the original data) and are, therefore, mostly limited to bi-temporal analyses.
Due to their shorter wavelength, which limits the penetration into the canopy, C-band SAR data are less suited for forest disturbances assessment than L-band (as evidenced in [16]), and have been used less than L-band data.Actually, monitoring forests using C-band was hampered, so far, by the lack of availability of high temporal and spatial resolution time series data.A new era started with the launch of the Sentinel 1A (April 2014) and 1B (April 2016) satellites developed by the European Space Agency (ESA), providing a large and unprecedented amount of free data for the operational needs of the Copernicus program.The dense time series of the Sentinel-1 (S1) constellation offer a unique opportunity to systematically monitor forests at a repeat cycle of 6 to 12 days, depending on the data type and the region in the world.In addition, the continuity of Sentinel data is guaranteed up to 2030, and the next generation of Sentinel is planned beyond 2030, allowing long-term environmental monitoring.With such temporal series of C-band data, new methods for forest disturbances monitoring are to be developed.
One of the first papers dealing with the use of S1 for forest disturbances detection has been published by Lohberger et al. (2018) [17].However, the methodology developed in [17], dedicated to fire detection, is mostly based on bi-temporal S1 datasets.Reiche et al. (2018) [18] developed original methods based on a combination of several sensors, such as S1 and PALSAR-2, together with optical data.However, these methods suffer from a crucial flaw: assuming that disturbances are necessarily characterized by a decrease in C-band backscatter within the disturbed area, which is not always the case.
Radar backscatter is indeed affected by factors related to forest biomass, structure, and ground conditions.Past studies agreed that the observed backscatter at C-band is a combination of the ground backscatter attenuated by the canopy layer and the backscatter from the canopy, which includes simple and multiple scattering, and the vegetation-ground interaction [19,20].The backscatter from vegetation canopy is affected by vegetation 3D structure and water content (related to biomass) [21].The ground backscatter at C-band is affected by soil moisture, surface roughness, and terrain topography [22].The pixel-based detection of deforestation using C-band SAR backscatter time-series can therefore be difficult, as deforestation is not necessarily characterized by a drastic change in backscatter within the disturbed area.Deforested areas may exhibit the same backscatter values as intact forests, for example, after rainfall events.Management practices may also induce misdetections: deforested areas that are cleaned or burnt are easily detectable by their lower backscatter, whereas deforested patches with large branches remaining on the ground may show similar, if not higher, backscatter values than intact forests, because of a branches-ground double-bounce scattering mechanism.Note that the double-bounce also increases with soil moisture.In Brazil, some large trees, such as the Brazil nut trees (Bertholletia excelsa) are protected by law, and are not cut in disturbed areas, which prevents SAR backscatter to strongly decrease, and hampers the detection of deforestation.
In this work, we introduce a new, more reliable indicator of deforestation based on Sentinel-1 time series, and we demonstrate how this indicator can be used for NRT deforestation mapping over a test site in Peru.The choice of this study area is motivated by the fact that the Peruvian Amazon has become an important deforestation hotspot in the last years [23], and by the availability of Sentinel-1 acquisitions in both ascending and descending orbits, which is an asset for the considered method, as will be explained in Section 2. The paper is organized as follows: Section 2 gives general information on the study area, provides information on the data used in this study, both SAR data and reference data, and describes the rationale and the development of the deforestation detection method.Section 3 presents the resulting estimates of deforestation and a validation assessment based on reference samples and on a comparison with the optical-based University of Maryland (UMD) Global Land Analysis and Discovery (GLAD) Forest Alert dataset.In Section 4, the further works needed to apply the method at large scales are discussed.

Study Area
The experimental part of this paper focuses on a study site of 600,000 ha (93 km × 65 km) located in the Peruvian Amazon, as shown in Figure 1.The study site lies on the eastern side of the Andes and covers the border between the San Martin and Loreto regions, around the city of Yurimaguas.The site contains a relatively preserved mountainous area, Cordillera Escalera, in its western part, while the eastern part belongs to the Amazon lowlands.The natural vegetation comprises mostly evergreen rainforest, with the presence of some seasonal deciduous forest.However, in the lowlands, the area is heavily degraded by smallholder agriculture (rice, papaya, vegetables) as well as recent agro-industrial development in relation to the palm oil and cocoa industry, which led to significant deforestation events.
The climate is equatorial, with large amounts of rainfall throughout the year (above 150 mm per month), except from June to August (about 100 mm per month).This dry season offers the best conditions for local farmers to practice slash-and-burn agriculture and for agro-industrial companies to clear new land for oil palm plantations, and therefore, deforestation is expected to occur mostly during these months.

Sentinel-1 Data
Over the Peru study site, we have used 92 S1 images acquired between 22 October 2014 and 21 September 2017, in both ascending (43) and descending pass (49).The images were acquired in interferometric wideswath (IW) mode, defined by ESA as the pre-defined mode over land.In this mode, images are provided at a resolution of 10 m with a swath of 250 km.Over the swath, the incidence angle ranges approximately from 29 • to 41 • .In this study, we used Level-1 ground range detected (GRD) products that consist of focused SAR that has been detected, multi-looked, and projected to ground-range using an Earth ellipsoid model.The time interval between two consecutive acquisitions with the same orbit orientation was 24 days until February 2017, and 12 days afterwards, with occasional acquisition failures leading to double intervals.Images in ascending orbit are acquired 3 days after images in descending orbit.The polarization mode was single VV until the end of 2016, with a gradual switch to dual VV + VH in the first half of 2017.
In order to handle efficiently the large amount of data available, a pre-processing chain was developed in Python.The chain automatically selects the data that intersect a chosen area of interest, and downloads them from the so-called Sentinel Product Exploitation Platform (PEPS) platform, managed by the French Spatial Agency (CNES), which provides free access to data from the Sentinel satellites (https://peps.cnes.fr/).Then, the CNES Orfeo ToolBox (OTB) utilities (https://www.orfeotoolbox.org/) are used to calibrate the mono-or multi-polarized images to obtain the γ 0 backscatter coefficient, and to orthorectify the images (projection Universal Transverse Mercator) at the 10 m spatial resolution.The digital elevation model from the Shuttle Radar Topography Mission (SRTM) at 30 m resolution is used during the orthorectification process.At this stage, the equivalent number of looks (ENL) is approximately 4.4.In order to increase this ENL, which is insufficient for most applications, we applied a multi-image filter, which decreases the speckle effect while preserving the spatial resolution of the images [24,25].A 3 × 3 spatial window was chosen, and each image was filtered only with the images acquired before its own acquisition date, in order to simulate the NRT conditions.The ENL increases with the number of images in the archive used for the filtering, up to a saturation level of 39.6.With the 18 images available in the archive at the beginning of the period of interest (March to October 2016, see Sections 2.2.2 and 3), the ENL is 27.4.

Reference Data
In order to validate the deforestation mapping results, we used reference samples derived from two Sentinel-2 images acquired on 10 March 2016 and 16 October 2016, which are the only almost cloud-free Sentinel-2 images acquired in 2016.The sample selection was achieved in two steps.First, we manually selected 94 deforestation samples through visual interpretation of the two Sentinel-2 images, as well as 32 undisturbed samples.Then, these manually selected deforestation samples were used to assess the spectral signatures of deforestation events (NDVI, red, green, and blue bands at both dates), and these signatures were used to build a decision tree that automatically extracts reference samples from the Sentinel-2 images.A further visual quality check was done to ensure that the extracted samples actually corresponded to deforestation.In total, 901 reference deforestation samples have been extracted automatically, which represent deforestation events that occurred between the two Sentinel-2 acquisitions.The location of all these reference samples is shown in Figure 2, and the numbers and sizes of the reference polygons are shown in Table 1.In addition, we compared our results with the deforestation patches detected in the University of Maryland (UMD) Global Land Analysis and Discovery (GLAD) Forest Alert dataset, a Landsat-based humid tropical forest disturbance alert system [8] that is currently available for 2016 and 2017 over the Congo Basin, Peru, Brazil and insular South-East Asia (http://glad.geog.umd.edu/alerts).

The Shadowing Effect as an Indicator of Deforestation
As mentioned in the introduction, deforestation is not always characterized by a significant change of backscatter within the disturbed area.
To get around this problem, we developed an alternative methodology which consists in detecting SAR shadowing.Shadowing occurs in SAR images because of the particular side-looking viewing geometry of SAR systems.A shadow in a SAR image is an area that cannot be reached by any radar pulse, because higher objects create an obstacle between the SAR antenna and this area.Consequently, no backscatter is recorded at these locations, and shadows represent dark areas in the SAR images.This phenomenon commonly occurs in mountainous areas where high peaks create shadows, but shadows created by trees at the borders between forest and non-forest areas can also be observed in high-resolution SAR images, depending on the viewing direction, as illustrated in Figure 3.As SAR satellites have a near-polar orbit, and therefore, an almost north-south flying direction (98.18 • inclination in the case of S1), only borders oriented approximately in a north-south direction can be identified.In the particular case of S1, with a right-looking antenna, borders between forest and non-forest, from west to east, are seen as a shadow area in ascending orbit, while borders between forest and non-forest, from east to west, are seen as shadow areas in descending orbit.When the trees in a forest patch are cut, a shadow appears or disappears at some of its edges, depending on the orbit orientation (ascending or descending orbit), the patch orientation, and the presence or absence of remaining forests around the patch, as described in Figure 4. Shadows that appear are characterized by a sudden drop of backscatter in the S1 time series.Thanks to the purely geometrical nature of the shadowing effects, this decrease of backscatter is expected to be persistent over time.New shadows should consequently remain visible for a long period.On the contrary, as previously mentioned, the decrease of backscatter that is sometimes observed within a deforested area is generally not stable in time, because of environmental conditions, regrowth, or forest management.Likewise, shadows that disappear are theoretically characterized by a sudden increase of backscatter, but these shadows are replaced by bare soil, with potentially low and variable backscatter, and can therefore be challenging to detect.
When a patch is deforested inside a forest, in addition to the shadow that appears on one of its edges (for a given orbit orientation), an opposite phenomenon is sometimes observed on the opposite edge: a backscatter increase, caused by a double-bounce mechanism between the newly created bare ground and the trees at the border of the remaining forest [26].However, because this effect depends on the backscatter of the bare ground, it is dependent on the environmental conditions, and is therefore more variable and less persistent than the shadowing effect.We therefore focused only on the detection of newly created shadows, as it is believed to be the most reliable indicator of deforestation.
Figures 5 and 6 illustrate the previously described phenomena.In Figure 5, the temporal backscatter profile of an area being logged (in black) exhibits a moderate decrease (~2.5 dB) when logging occurs (in October-November 2015), but the post-disturbance backscatter then gradually increases to its original level in about 6 months, probably because of regrowth, with additional variations probably caused by soil moisture variations.In one of the edges of the logged patch, a shadow appears (in red), which is characterized by a drastic backscatter decrease (~5 dB), with no apparent evolution after the disturbance.An intact forest close to the deforestation area (in green) is also depicted, and shows a stable backscatter.It can be noted that in the case of a pure shadow, the backscatter should drop to the noise equivalent sigma zero (NESZ) value linked to the thermal noise of the instrument, which should remain below −22 dB for Sentinel-1.In practice, however, the pixels of the shadow areas associated to deforestation events are usually not pure, and are contaminated by adjacent forest or bare soil areas, leading to backscatter values higher than the NESZ but still significantly lower than the backscatter of forests, as can be seen in Figure 5.   Therefore, it appears that the shadowing effect, as described in this section, could be used as an indicator of forest removal at the stand level, either from anthropogenic causes (clearcutting for wood harvest, for conversion to farms, ranches or urban use, or for plantation renewal or extension) or from natural causes (fire, wind).In the following, we refer to these changes as deforestation, even though some of these changes would not qualify as deforestation according to the definitions of other stakeholders, depending on whether plantations are included in forest or not, and on whether deforestation is defined as a long-term or as an immediate process.In any case, other forest disturbance types, such as forest degradation through selective logging or insect outbreaks and diseases, should not be detectable through the shadowing effect at the Sentinel-1 resolution.
The deforestation detection system that we recommend to adopt when using S1 time series is composed of two steps: 1.
Detect shadows that appear or disappear in a series of images 2.
Reconstruct the deforested patches associated to the shadows These two steps, as well as the implications of the different configurations presented in Figure 4 on each of these two steps, are further discussed in Sections 2.3.2 and 2.3.3,respectively.

Detection of Shadows
First of all, in order to be detectable, shadows have to be large enough relative to the resolution of S1 images.Over flat terrain, the width of the shadows W is expressed as a function of the tree height H, and the SAR incidence angle θ: W = H tan (θ).The width of the shadow is given in Figure 7 for the range of incidence angles of S1 (between 29 • and 46 • ) and for tree heights between 0 and 40 m, with the corresponding number of pixels.For example, a shadow of 10 m (i.e., a S1 pixel) occurs when trees are greater than 10 m at 45 • and 18 m at 29 • approximately.In the case of non-flat terrain, with a local slope α oriented towards the sensor (α > 0) or away from the sensor (α < 0), the shadow width is multiplied by a factor K = cos (θ)*cos (α)/cos (θ − α).For moderate slopes below 10 • , this implies a reduction of the shadow width limited to 15% for slopes oriented towards the sensor, and an increase of the shadow width limited to 22% for slopes oriented away from the sensor.These geometrical characteristics restrict the use of this method to dense forests with high trees.According to the tree height map of Simard et al. (2011) [27], average tree height in the Peru site is above 35.5 m in the lowland and 32.5 m in the mountainous area, which ensures that shadows are detectable (more than 1.5 pixel wide) in flat terrain and over moderate slopes.Detecting the abrupt backscatter drops that characterize a newly created shadow requires adopting a change detection approach.SAR-based change detection has largely been based on two-date intensity ratios [14,17,[28][29][30][31][32], which were found to be better adapted to the statistical characteristics of SAR data than the intensity difference [33].In this study, we used the radar change ratio (RCR) [34], which consists of a ratio between the post-and pre-disturbance averaged backscatter.Compared to the more simple two-date ratio, the RCR has the benefit to be less sensitive to remaining speckle, and to ensure that the change persists over a certain period.In practice, in a time-series of N dates, the backscatter of a given pixel at date d i (i ∈ (1, N)) being noted γ 0 i , the changes that occur between date d i and date d i+1 are measured by the following radar change ratio: RCR i = M a /M b where M b is the mean backscatter in X b images before date d i (included), and M a is the mean backscatter in X a images after date d i+1 (included): γ 0 j and M a =  This change indicator can be computed for i ∈ [X b , N − X a ], and X b and X a represent the number of images on which the mean is calculated before and after the considered acquisition date.X b and X a must be chosen with respect to the following considerations:

•
Higher values of X b and X a have the advantage to reduce the speckle effect and limit the detection of areas with intrinsically variable backscatter (e.g., crops).

•
However, high values of X a will delay the effective detection of deforestation, and will therefore hamper the NRT capacity of the approach.A trade-off must therefore be found between speckle filtering on one hand, and timely provision of results on the other hand.In this study, we chose X a = 3, which proved to be sufficient in terms of speckle reduction.

•
In principle, an X b can be chosen that is as high as possible.Therefore, we chose to average all images acquired before the considered date, in order to reduce the sensibility of the change detection to seasonal and environmental effects affecting the backscatter.
Within a given S1 time series, a shadow is detected if the following rule is verified: min (RCR i ) < Y dB, where Y is a threshold to be defined, that will depend on the polarization, the local incidence angle, and characteristics of the study area.The date frame of the detected change is [d k , d k+1 ], where k verifies: RCR k = min (RCR i ). Figure 9 shows one example of the backscatter γ 0 i time series, and the radar change ratio RCR i time series over one area identified as a new shadow.The point of maximum negative change (around −6 dB) corresponds to the detected disturbance date.In this study, we identified shadow areas by applying a −4.5 dB threshold on the image formed by the temporal minimum RCR value in both ascending and descending orbits in VV polarization, over the full 2016 year.Although both VH and VV polarizations have been available, simultaneously since February 2017 (fine beam dual mode FBD), preliminary tests demonstrated that the two polarizations perform similarly.We therefore restricted ourselves to VV only, to keep a consistent time series over the full date range.Small outliers are filtered by sieving, to keep only segments of more than 4 pixels (at a 10 m pixel size).

Reconstruction of Deforested Patches
In the case when both ascending and descending orbits are available, as in our study site, and deforestation occurs within a larger forest patch (first row in Figure 4, which is the most common configuration), an efficient reconstruction of the deforested patches would consist in associating the detected shadows into pairs (one in ascending and one in descending, occurring in a given neighborhood within a reduced timeframe), and in applying a boundary operator (e.g., convex envelope) around these two edges to delineate the deforested patch, as illustrated in Figure 10.However, for a more general approach adapted to areas where only a single orbit orientation is available (see Section 4), and to other configurations than logging occurring fully within a forested area (rows 2 to 4 in Figure 4), we used a different reconstruction strategy.In this approach, the shadow areas detected in the first step are used as a starting point to detect neighboring pixels that exhibit a slighter backscatter decrease.To do so, we apply a −3 dB threshold over the same minimum RCR image than that used for the shadow detection, and small outliers are filtered by sieving, to keep only segments of more than 10 pixels (at a 10 m pixel size).Among the many potential deforestation patches generated by this step, only those that intersect with a detected shadow are kept, and constitute extended shadows.Then, the final delineation of the deforested patches is obtained by applying the MATLAB "boundary" function (based on alpha shapes [35], https://fr.mathworks.com/help/matlab/ref/boundary.html) to each extended shadow, with a shrink factor equal to 0.6.The shrink factor is a scalar in the range of [0,1] that is used to tighten or loosen the boundary around the points.A shrink factor of 0 corresponds to the convex hull of the points, and 1 corresponds to the tightest single-region boundary around the points.The value 0.6 was chosen as a trade-off value through visual assessment.This reconstruction strategy is illustrated in Figure 11.

Post-Processing: Masking Undesirable Areas
The method proposed in this study is based on the detection of persistent backscatter decreases.Other events than deforestation can cause such backscatter decreases and create false alarms.In theory, it would be advisable to use a forest/non-forest (FNF) mask in order to restrict the deforestation detection to forested areas.In practice, such FNF masks can be obtained in an up-to-date and accurate form only locally, using manual or semi-automatic methods, but obtaining them for large-scale applications remains a challenge.
In the study site, we identified that the two potential areas that can create false alarms are agricultural fields upon harvest, and rivers that can dramatically change their courses in this very flat area.Agriculture did not turn out to be a major source of error, as fields are very small and the backscatter decrease linked to harvest is of small amplitude and not persistent in time.Rivers changing their courses were masked out by applying a threshold on the minimum VV backscatter image.
Topographical heterogeneities can also be a major source of uncertainties in radar applications.However, in this study, the use of the RCR, which is a backscatter ratio, limits the effect of topography, which is essentially a multiplicative effect.Terrain-induced shadows would not provide false alarms in the deforestation detection because they are stable in time for consistent orbits, but they would result in omission errors.The use of ascending and descending orbits nonetheless reduces the amount of such areas that are not observed.Besides, deforestation does not seem to occur over very steep slopes in this part of the world.Therefore, we did not mask out areas with steep slopes in this study, although this might be necessary in other parts of the world [14].

Results
Following the approach described in Section 2, the available S1 dataset was used to produce a map of all deforested areas detected in the year 2016, which is shown in Figure 12.The deforestation map is evaluated in terms of detection rates using the 901 automatically selected deforestation reference samples described in Section 2.2.2.The detection rates correspond to the percentage of reference samples that are correctly detected by a given method, where a sample is considered detected when at least 10% of its area has been seen as deforested.It therefore corresponds to the producer's accuracy of the "disturbance" class.The detection rates of the S1-based method described in this paper are given in Figure 13 (red line) for a number of sample size ranges: 0-0.2 ha, 0.2-0.4ha, 0.4-0.6 ha, 0.6-0.8ha, 0.8-1 ha, 1-1.5 ha, 1.5-2 ha, 2-3 ha, 3-4 ha, 4-5 ha.The same detection rates were also calculated for the UMD-GLAD Forest Alert dataset (blue line).The number of reference samples in each size range is also reported (black line).For the S1-only approach, we considered only detections occurring effectively between the dates of the two Sentinel-2 images used for the sample extractions (10 March to 16 October 2016), while for the UMD-GLAD dataset, we considered all the detections occurring in 2016, because of the lower observation rate linked to the poor availability of cloud-free Landsat observations.The detection rate is always higher with the S1-based approach than in the UMD-GLAD Forest Alert dataset, except for the very large samples (over 3 ha) where both methods reach 100% detection (with very few samples though).In particular, the S1-based approach reaches a 95% detection rate already for samples in the range 0.4-0.6 ha.Values of about 90% detection rate are reached by the UMD-GLAD Forest Alert only above the 1.5-2 ha range, partly because of the coarser resolution (30 m) of the input data.We also calculated the overall sample detection rates and area detection rates of both approaches for a given minimum mapping unit (MMU), which is commonly included as a user requirement.Those detection rates correspond to the percentage of reference samples or percentage of area (of the reference samples) correctly detected as deforested when only reference samples larger than the MMU are considered.The overall sample and area detection rates are always superior for the S1-based approach with respect to the UMD-GLAD Forest Alert system, with the larger discrepancy in the area detection, as can be seen in Figure 14.In addition, we show in Figure 15 the histograms of the detection dates of the reference samples for both approaches, as well as the histogram of the difference in the detection date between the two datasets.We found that the S1-based approach provides deforestation detection on average 41.6 ± 34.7 days earlier than UMD-GLAD Forest Alert, with detections 3 to 4 months earlier not being uncommon.Our results indicate that the maximum deforestation rates occur in June, July, and August, the three driest months of the year, which is expected.In order to assess, also, the commission (false-alarm) rates, the manually selected undisturbed plots were used together with the automatically selected disturbed plots, to compose the pixel-based confusion matrix shown in Tables 2 and 3 for the S1-based approach and the UMD-GLAD Forest Alert dataset, respectively.Due to the lack of very high-resolution optical imagery available over the whole study area, and because the Sentinel-2 images cannot be interpreted unambiguously over random or systematic samples, it was not possible to adopt a probability sampling design for the selection of the validation samples, as recommended in [36].As a consequence, we could not derive the error matrix expressed in terms of estimated area proportions, but only in terms of pixel counts.The User's Accuracy (UA) and Producer's Accuracy (PA) figures presented in Tables 2 and 3 are therefore only indicative, and should be considered with caution.

Discussion
The results presented in Section 3 indicate that our method, which is based on the SAR shadow detection approach described in Section 2, provides better deforestation results than the only NRT product currently available over Peru, the Landsat-based UMD-GLAD Forest Alert dataset, both spatially (Figures 13 and 14) and temporally (Figure 15).In particular, the UMD-GLAD dataset is limited by the resolution of the Landsat data (30 m) for the detection of small deforestation samples, and by cloud cover for the timely detection of deforestation events.
For the deforested class, the confusion matrices indicate a slightly lower UA for the S1-based method compared to the Landsat-based method (99.4% vs. 99.9%),but a much higher producer's accuracy (80.3% vs. 40.2%).
This study in Peru is a proof-of-concept that the SAR shadowing effect is an indicator of deforestation that can be exploited in NRT deforestation monitoring systems.The method described in Section 2.3 should, however, not be received as a turnkey method that could be readily used globally.An operational application of this approach at large scale indeed requires further works, in particular:

•
Thresholds: The method involves several parameters and thresholds, as described in Sections 2.3.2 and 2.3.3.In particular, the thresholds that were applied on the minimum RCR image to detect shadows and potentially deforested areas (−4.5 dB and −3 dB respectively) were chosen empirically in this study, and might need to be adapted locally in order to account for the polarization, the local incidence angle, and the characteristics of the study area (for example, the type of deforestation, which will impact, more importantly, the second threshold).One way to define these thresholds automatically is to use statistics derived from reference samples, if available.The minimum RCR can be approximated as an intensity ratio.In that case, it has been shown that the optimal threshold r opt to distinguish between two classes of mean intensity ratios r A and r B (r A < r B ) has a complex expression that involves other parameters that cannot be estimated in the general case (such as the proportion of each class), but can be approximated using the particular value r 0 = √ r A r B [37].For example, the mean intensity ratio of shadows can be estimated using the mean values of the first quartile of the minimum RCR image of each reference deforestation sample, under the assumption that shadows cover approximately 25% of the deforested areas.This leads to r A = −7.2dB in the manually selected deforestation samples.The mean intensity ratio of non-shadows is approximately r B = −2.1 dB (mean value of the minimum RCR image in the scene).This leads to a threshold of r 0 = −4.65 dB, very close to the −4.5 dB value that we defined empirically.Regarding the second threshold used for the detection of potentially deforested areas, r A and r B can be calculated as the mean values of the minimum RCR image in the manually selected deforestation and undisturbed reference samples, and are found to be equal to r A = −4.6 dB and r B = −2.1 dB, leading to a threshold of r 0 = −3.35dB, again, very close to the −3 dB value defined empirically.
When reference samples are not available, a set of standard values could be defined, e.g., as a function of the local dominant drivers of deforestation.

•
Orbits: This approach, which detects shadows that appear simultaneously in descending and ascending orbits, allows obtaining a fairly good delineation of deforested patches, especially in the configuration where deforestation appears fully within a larger forest patch (first row of Figure 4).This kind of approach can be implemented in areas where both orbit orientations are available.According to the current S1 observation scenario depicted in Figure 16, these areas concern Europe, the western part of the Americas, eastern Africa, and some parts of Asia.Significant gaps for tropical deforestation include most of the Amazon and Congo river basins, which are covered by one orbit orientation only.In these cases, the detection of shadows will be less complete as only one edge can be detected in the best case, or even none in some configurations (see Figure 4).Therefore, when only one orbit orientation is available, in addition to the detection of shadows that appear, it could be worthwhile investigating also the detection of shadows that disappear, as well as the detection of backscatter increases through double-bounce, to complement the edge detection.

•
NRT: The sensitivity of the method was demonstrated for X a = 3, which represents a delay of about 1 month with a revisit frequency of 12 days.For truly NRT applications, a shorter delay is required, and this parameter involved in the calculation of the RCR needs to be reduced to 1 or 2. This would, however, cause an increase in the false-alarm rate, because of confusions caused by the speckle effect.In that case, an alert system can be proposed with X a = 1 or X a = 2, keeping in mind that its results would need to be confirmed/rejected, and complemented afterwards using higher values of X a .

•
Seasonality: In areas where natural forests exhibit a strong seasonality (e.g., deciduous forests) which are reflected in the annual C-band backscatter profiles, this seasonality should be taken into account in the calculation of the RCR.For example, in the calculation of M a and M b , γ 0 j should be replaced by γ 0 j − t j , where t j represents the mean temporal backscatter trend of natural forests, which can be retrieved from the archived S1 time series over selected undisturbed forested pixels.

•
Synergies with other sensors: The Sentinel era offers a unique opportunity to exploit the synergies between optical and SAR sensors.Whether the approach based on Sentinel-1 described in this study will be sufficient in itself to reach operational levels remains to be demonstrated, but it anyway has a strong potential to improve current and future deforestation monitoring systems, which rely mostly on optical imagery.As the continuity of the Sentinel-1 and -2 sensors are guaranteed until at least 2030, investigating, more thoroughly, the synergies between both kinds of sensors would be a well-invested effort.The S1-based approach can also be applied to future high-resolution SAR systems operating with high repetition, in particular, the future L-band NASA-ISRO Synthetic Aperture Radar (NISAR) sensor that is planned for launch in 2020-2021.

Conclusions
In this paper, we introduced a new indicator of forest disturbances that can be derived from SAR time series, in the form of shadow areas that appear at the edges of deforested patches.We described how and when these shadows appear or disappear, and how this phenomenon can be exploited to detect deforestation events.With the launch of S1, the availability of free-of-charge high-resolution SAR data, with a global coverage and high temporal repetition, has allowed testing, for the first time, the potential of this approach.We demonstrated the sensitivity of the method in a 600,000 ha test site in the Peruvian Amazon, by obtaining better detection rates than the UMD-GLAD Forest Alert dataset, a Landsat-based NRT deforestation detection system, and a better temporal characterization of deforestation events.This sensitivity shall be exploited in the future for operational applications.

Figure 1 .
Figure 1.The study site, located in Peru (red square), on an optical background image from Google Earth (© 2018 Google).Features that can be seen in the right pane include the city of Yurimaguas along Rio Huallaga in the northeastern corner, the Cordillera Escalera mountain range in the western part, and a large oil-palm plantation in the central eastern region.

Figure 2 .
Figure 2. Location of the reference samples: manually selected deforestation (yellow) and undisturbed (green) samples, and automatically extracted deforestation samples (red).Optical background image from Google Earth (© 2018 Google).

Figure 3 .
Figure 3. Illustration of the synthetic aperture radar (SAR) shadowing effect at the border between forests and deforested areas, for the descending and ascending orbit configurations.(Background image: ConservationDrones.org/OrangutanConservancy.)

Figure 4 .
Figure 4. Shadows that appear or disappear when a forest patch is deforested, depending on the orbit orientation (ascending or descending), the patch orientation, and the configuration (patch inside a forest, at the edge of a forest, or isolated).The figure is simulated for the Sentinel-1 configuration (right-looking with 98.18 • orbit inclination) which is summarized in the frame on the top right.

Figure 5 .
Figure 5. Temporal backscatter profiles of an area being logged (black), the corresponding shadow (red) that is created on one of its edges when deforestation occurs (October-November 2015), and an intact forest (green).

Figure 6
Figure6shows a subset of Sentinel-1 intensity images acquired in descending orbit on 10 May 2015 and 1 October 2016 over a rubber plantation area.The shadows created at the eastern edge of stands that have been harvested between the two dates appear as vertical dark lines, indicated by blue arrows.The backscatter within the harvested stands shows low values for one recently harvested stand, and similar values as the neighboring intact stands for the other stands harvested a few months before the second date.

Figure 6 .
Figure 6.A 2.6 km × 3.4 km subset of Sentinel-1 intensity images acquired in descending orbit on 10 May 2015 and 1 October 2016 over a rubber plantation area, centered on 11 • 16 40"N 106 • 39 40"E.Blue arrows indicate shadows that appear between the two dates because of harvest.

Figure 7 .
Figure 7. Shadow width as a function of tree height and Sentinel-1 incidence angle.The corresponding number of pixels is based on a 10 m pixel size.

Figure 8
Figure 8 illustrates how the change indicator RCR i is computed for two [d i , d i+1 ] timeframes in an area where a shadow appears, for X b = 5 and X a = 3.The first example corresponds to the RCR calculated on the exact [d i , d i+1 ] timeframe during which the shadow appears, and leads to RCR i = −5.2dB.The second example corresponds to a later [d i , d i+1 ] timeframe when no change occurs, with RCR i = 0.5 dB.

Figure 8 .
Figure 8. Temporal variation of the VV backscatter in an area detected as a new shadow (red), and illustration of the calculation of the change indicator RCR i (radar change ratio) for X b = 5 and X a = 3 for two pairs of dates: (a) the exact date interval during which change occurs; (b) a date interval when no change occurs.

Figure 9 .
Figure 9.Time series of the backscatter (red) and the corresponding RCR (blue) over a 20-pixel area where a new shadow appears.The detected logging date (black dashed line) corresponds to the date when the RCR is minimum.

Figure 10 .
Figure 10.Delineation of deforested patches in a 1.4 km × 1.9 km subset by applying automatic convex envelopes around the pairs of shadows detected in ascending (red) and descending (blue) orbits, respectively corresponding to the western and eastern edges of the deforested patch.Optical background image from Google Earth (© 2018 Google).

Figure 11 .
Figure 11.Detection and reconstruction of deforested patches in Peru in an 800 m × 800 m subset: detected shadows (red) and potentially deforested areas (blue).Isolated blue patches are rejected while blue patches connected with red patches are combined to form extended shadows, which are used to reconstruct the boundary of the deforested patch (in black).(Background optical image from Google Earth.)

Figure 12 .
Figure 12.Deforestation map for year 2016 over the study site.In the 1.5 km × 2.0 km subset (left), the automatically selected deforestation reference samples are shown as red polygons.(Background optical image from Google Earth, in greyscale).

Figure 13 .
Figure 13.Detection rate of deforestation events per reference sample size range (delimited by thin black dashed lines), for the S1-only method from CESBIO (red) and for the UMD-GLAD (University of Maryland Global Land Analysis and Discovery) Forest Alert dataset (blue).The number of reference samples available in each size range is also shown (thick black dashed line).

Figure 14 .
Figure 14.Overall sample (left) and area (right) detection rate of both approaches with respect to the minimum mapping unit (MMU) of the reference samples.

Figure 15 .
Figure 15.Number of detected reference samples per month for the S1-only method and for the UMD-GLAD Forest Alert system (left), and histogram of the difference in detection date of the reference samples (right), where negative values indicate an earlier detection in the S1-only method.

Table 1 .
Numbers and sizes of the reference samples constituting the validation database for deforestation assessment.

Table 3 .
Confusion matrix of UMD-GLAD Forest Alert dataset.