Continuous Forest Monitoring Using Cumulative Sums of Sentinel ‐ 1 Timeseries

: Forest degradation is recognized as a major environmental threat on a global scale. The recent rise in natural and anthropogenic destruction of forested ecosystems highlights the need for developing new, rapid, and accurate remote sensing monitoring systems, which capture forested land transformations. In spite of the great technological advances made in airborne and spaceborne sensors over the past decades, current Earth observation (EO) change detection methods still need to overcome numerous limitations. Optical sensors have been commonly used for detecting land use and land cover changes (LULCC), however, the requirement of certain technical and environmental conditions (e.g., sunlight, not cloud ‐ coverage) restrict their use. More recently, synthetic aperture radar (SAR) ‐ based change detection approaches have been used to overcome these technical limitations, but they commonly rely on static detection approaches (e.g., pre and post disturbance scenario comparison) that are slow to monitor change. In this context, this paper presents a novel approach for mapping forest structural changes in a continuous and near ‐ real ‐ time manner using dense Sentinel ‐ 1 image time ‐ series. Our cumulative sum–spatial mean corrected (CUSU ‐ SMC) algorithm approach is based on cumulative sum statistical analysis, which allows the continuous monitoring of radar signal variations, derived from forest structural change. Taking advantage of the high data availability offered by the Sentinel ‐ 1 (S ‐ 1) C ‐ band constellation, we used an S ‐ 1 ground range detected (GRD) dual (VV, VH) polarization timeseries, formed by a total of 84 images, to monitor clear ‐ cutting operations carried out in a Scottish forest during 2019. The analysis showed a user’s accuracy of 82% for the (conservative) detection approach. The use of a post ‐ processing neighbor filter increased the detection performance to a user’s accuracy of 86% with an overall accuracy of 77% for areas of a minimum extent of 0.4ha. To further validate the detection performance of the method, the CUSU ‐ SMC change detector was tested against commonly ‐ used pairwise change detection approaches for the same period. These results emphasize the capabilities of dense SAR time ‐ series for environmental monitoring and provide a useful tool for optimizing national forest inventories.


Introduction
Forests are one of the major contributors to climate change mitigation, playing an integral role in the global carbon cycle by removing 2.1 Gt CO2 per year [1,2].However, accelerating changes in global climate dynamics, together with high pressures exerted by current socio-economic forces, are leading to a dramatic reduction of global forested regions, which is threatening the conservation of world forest ecosystems [3][4][5][6][7].Among the many land use and land cover changes (LULCC) that threaten global ecosystems, deforestation has been recognized as the most important human-induced land transformation during the last half century [8].In addition, prediction models based on the current global environmental trends forecast a global canopy cover reduction of 223 million hectares by 2050 [9][10][11].These facts highlight the need to develop accurate forest monitoring systems which contribute to an optimal management and conservation of the world's forests.
At present, there are numerous institutions and programs that seek to reverse the current deforestation trend through joint ecological agreements.Initiatives such as the United Nations Reducing Emissions from Deforestation and Forest Degradation (REDD+ program), seeking to develop a joint and multi-scale global forest inventory; the UK Space Agency Forest2020 project, aiming to restore 300 million hectares of tropical forest while improving forest monitoring systems in developing countries [12]; or the European Space Agency Climate Change Initiative (ESA-CCI) [13], exploiting the potential of earth observation technology for climate change response, are some examples of international projects aiming to protect forest ecosystems through the development of improved environmental monitoring systems.
Forest inventories are recognized as the standard tool for quantifying and monitoring forest characteristics [14].However, the highly dynamic nature of forests normally renders inventories inefficient tools for monitoring the continuous transformation processes that occur in these ecosystems as they are mostly generated through infrequent, sometimes one-off, ground surveys.Remote sensing (RS) techniques have emerged as valuable alternatives for overcoming this challenge, thanks to the advances that aerospace technology and geographical information systems (GIS) have accomplished over the last decades [15].
The work presented here utilised synthetic aperture radar (SAR) to address the inadequacies of forest inventories.Historically, the use of RS for forestry purposed has focused on optical systems, however this approach is hindered in areas with recurrent cloud coverage (e.g., Tropical areas) or limited solar illumination (e.g.Polar latitudes) [16][17][18].Another disadvantage of optical systems derives from the fact that they can only extract information from the very top canopy layer.This limits the use of optical methods for real time forest monitoring, since under-canopy disturbances may be unnoticed or clouds may lead to a long wait time before the disturbance is assessed.
SAR basic imaging principles rely on the measurement of the energy that is redirected back to the radar antenna after an interaction occurring between the transmitted signal and the object studied.Depending on the intensity of the backscatter, it is possible to determine diverse information about the physical properties of the target [19].Radar signals are highly influenced by three main physical parameters which are inherent of any object or surface target: the geometry, roughness, and moisture content.Any changes in these parameters in a forest may result in a noticeable variation in the backscattered intensity [20,21].
The integration of EO datasets into forest monitoring systems has generated a significant impact in forest management and policy making.Products such as the JICA-JAXA Forest Early Warning System in the Tropics (JJ-FAST) [22] or the UMD-GLAD (University of Maryland Global Land Analysis and Discovery) Global Forest change map [23]; the global forest/non-forest maps from ALOS PALSAR data [24], or the TanDEM-X forest/non-forest global map [25] are widely used by international policymakers for reviewing the annual conservation status of the worldʹs forests.However, the use of these products presents diverse limitations.Technical and environmental aspects such as the coarse spatial resolution (i.e., 100-1000m) or low temporal resolution, due to persistent cloud coverage, hinder the detection of small forest degradation events as well as the generation of sub-annual or rapid response alerts [14].
The launch of the European Space Agency's (ESA) Sentinel-1 satellite constellation in 2014, opened the door to new opportunities for environmental monitoring.With freely available data and improved sensor characteristics (e.g., dual polarization 6-or 12-day revisit period), Sentinel-1 rapidly stood out as a valuable resource for effective monitoring of highly dynamic forest ecosystems [17,26].

Detecting Forest Changes With Pairwise Analysis of SAR Data
In this work we have adopted an innovative approach which uses long time series of data to detect forest change.However, to determine the effectiveness of this approach, we compared it to several methodologies which only use pairwise difference or ratio between consecutive images.As will be shown by our analysis these methods present some limitations.
Previous studies have shown the capabilities of SAR systems for forest change monitoring at lower scales [16,[27][28][29][30][31][32][33].Most of these studies were based on pairwise or single epoque image comparison for detecting the changes between two specific dates.Others investigated the use of time series for developing near real-time change detectors by using multitemporal averages (e.g., monthly average).A major limitation to this approach is that when comparing images acquired over a relative small time gap, possible changes occurring during this time gap may go undetected [32,34].For instance, direct pairwise comparison between closely acquired images could lead to an erroneous estimation of land cover changes associated with divergences in the environmental conditions.Minor variation in some land surface parameters (i.e., soil moisture, canopy moisture, plant phenology) may result in radar backscatter changes between the compared images.Another limitation is that often there is a time lag between the change happening and the SAR images showing the change.Several studies [32,[34][35][36] have reported significant delays in the forest change detection that range from weeks to months after the disturbance.This phenomenon was also observed by Watanabe et.al [36], where the reduced sensitivity to change was related to the persistence of branches left on the ground after the felling activities.Consequently, there are inherent difficulties in selecting the right time gap for performing pairwise difference analyses.Pairs with, for example, between 6-or 12-days difference may not be able to detect changes, whereas using larger time gaps can lead to the identification of changes due to seasonality.Lastly, pairwise change detection methods do not allow the advantages associated with long-time series multitemporal information to be realised, which facilitates the search for specific trends over the time domain.Therefore, building the change detector algorithm on parameters extracted from dense datasets provides a higher level of confidence when assessing the monitoring processes.
This study aims to develop a novel approach for monitoring forest structural changes in a continuous manner using Sentinel 1 image time-series data.Based on cumulative sum statistical analysis, this method was developed for exploiting the continuous variation in radar signals derived from forest changes (e.g., logging activities).This approach was validated in the Queen Elizabeth Forest Park (Scotland) in collaboration with Forest Research (FR) and Forest and Land Scotland (FLS).However, the new forest monitoring tool is designed to be applicable to any UK National Forests.
To summarize, the main novelties of this work are 1) the design of statistical tests for cumulative sums using SAR data, 2) the design of a novel SAR-based change detection algorithm for continuous forest monitoring (CUSUM-SMC) and 3) the quantitative evaluation of the change detection performance of the statistical tests for forest structural disturbances.

Study Area
The study area was located within the boundaries of the Queen Elizabeth Forest Park (19,665 ha), near Aberfoyle village, Scotland (Figure 1).This area is part of the Loch Lomond and Trossachs National Park and has been widely used as a research forest by FR.The area is characterized by large natural forests and forest plantations.Main forest types, according to the most recent (2016) National Forest Inventory (NFI), are 'conifer woodlands' (9,001.34ha) and 'broadleaved and mixed woodlands' (1,547.76ha).There is a clear predominance of plantations of Sitka spruce (Picea sitchensis) (5,351 ha) with smaller areas of Pinus contorta (468 ha) and Quercus petreae and robur (612 ha).According to historical data (1981-2010) obtained from the closest Met Office climate station at Arrochymore (56.096, -4.548), 30 m above sea level, the bioclimate of this region is considered as 'Cfb'-'oceanic or highland climate' (Köppen-Geiger classification).The mean annual temperature for the area is 12.6 °C and annual rainfall (1735 mm) is distributed evenly throughout the year.July is the hottest month with an average temperature of 19.2 °C, while January shows the lowest temperatures with an average of 6.9 °C.The soils are drifts derived from acid metamorphic rocks and Lower Old Red Sandstone sediments, classified as "brown soils" [37].To ensure the absence of any significant geometric distortion that may compromise the use of radar images, a detailed analysis of the main topographic variables (e.g.elevation, aspect, hill shade, orientation) was carried out.Results showed a predominance of lowland and foothills with generally gentle slopes, although some hills do present a significant slope, which did not condition the utilization of radar in the area of study.

SAR Data
Dual polarization (VV, VH) Sentinel-1 Level-1 ground range detected (GRD) images, were retrieved from the ESA data archive using the Sentinel Data Access Service (SEDAS) platform.S1-IW_GRD_HR products were acquired in interferometric wide swath (IW) high resolution (HR).The center frequency was 5.405 GHz (which leads to a 5.3 cm wavelength) and with around 20 x 22 meters resolution [38,39].GRD products were focused SAR data that had already been converted into intensities, multi-looked, and projected to ground range using the Earth ellipsoid model WGS84 [40].In total, the radar time-series was composed of 84 scenes, downloaded for the period between 01 January 2018 and 30 September 2019.All images were acquired in ascending mode and using a unique orbit (i.e., relative orbit 30), aiming to maintain the identical radar geometry in all acquisitions.

SAR Data Pre-Processing
All images were pre-processed using ESA-SNAP 7.0 and ESA-NEST 4.1 software and following a four step procedure: 1. subset, 2. radiometric calibration, 3. coregistration, and 4. geometric terrain correction.First, the subset process was carried out to eliminate those areas located out of the boundaries of the region of interest.Second, radiometric calibration was applied to convert backscatter values to sigma-naught (0).It should be noted that in cases where the analysis considers multiple orbits, the use of sigma0 as calibration will not be effective.A calibration process able to remove geometrical distortions associated to steep slopes will be required (e.g., gamma0).The use of a multi-orbit approach was beyond the scope of the work presented here, however; the benefits of this approach will be explored in the future.. Third, data coregistration served to precisely stack together all images into a single multi-layer (composite) product.Fourth, terrain correction allowed the geolocation of all pixels, using the 30 m Shuttle Radar Topography Mission (STRM) digital elevation model (DEM) and subsequently projected to the (WGS84, UTM 30N) coordinates system.The final output resulted in a data stack formed by 168 images (84 co-polarization, 84 crosspolarization).

Validation Data: Forest Mask and Logging Activities Reference Data
The forest mask was obtained from the 2018 United Kingdom National Forest Inventory [41].The forest definition used for the forest mask includes both conifer and broadleaves areas, regardless of their pure or mixed typology.We did not discriminate between forest species because of the mixed nature of the forest and the impossibility of differentiation, at the pixel level, between typologies of forested areas.
All of the developed algorithms were tested for detecting clear-cut activities carried out in the Queen Elizabeth forest park between 1 January and 30 September 2019 (Figure 2b).We worked with the harvesting forest office of the FLS with the aim of gathering and verifying information about the start and end dates of each logging activity.
To obtain the initial draft of the change map, we looked for any clear optical images (cloud coverage <10%) acquired over Aberfoyle over this period.A pair of Landsat 8 (OLI) shortwave infrared images acquired on 25 June 2018 and 28 June 2019 (a year later) provided the best characteristics for the required period (Figure 2a).Shortwave infrared is a band combination highly recommended for studying vegetation health and stress, change detection, and disturbed soils [42].It combines the Short-wave Infrared (SWIR), Near infrared (NIR), and red bands to generate images where vegetation can be clearly differentiated from bare soils.Using a visual interpretation process, we manually identified forest areas which had experienced significant transformation.Finally, once all the three land classes (change-forest, stableforest, non-forest) were mapped as vector polygons, we extracted the pixel's boundaries of each particular area by overlapping these polygons over a terrain corrected Sentinel 1 image, extracted from the cube.This process generated a total of 30 areas, 10 for each land class, of various sizes (from 2.5 to 50 ha).These areas were used for the signal behavior analysis, as shown in the following section.

Cumulative Sum Change Detection Analysis
In this study we adapt cumulative sum (CUSUM) statistical methods for monitoring forest disturbances over dense time series.CUSUM has been widely used by the financial sector, being recognized as an excellent big data analytic tool for detecting changes in the market [43].This method however has not been extensively used in scientific and environmental applications.Recent research [43][44][45][46] showed the potential capabilities of this methodology for monitoring environmental variations.Manogaran et al. [43], proposed a novel climate change detection algorithm to monitor changes in the seasonal climate.The preliminary results showed a high robustness of the method when detecting variations in climate.Kellndorfer [44] explores the use of CUSUM for monitoring deforestation and forest degradation processes, however no quantitative validation of the method is performed.In this work, we developed a novel methodology which is based on the concepts of the CUSUM.
CUSUM is used to monitor the deviation of a variable from its mean based on measurements acquired over a given time frame.It stands as a powerful statistical method for analysing multitemporal processes, since it allows the detection of both slow and abrupt variations in the mean value of a quantity of interest [43,46].In most cases, these changes in the mean are associated with temporal or permanent changes in the variable analysed [47], a fact which makes this method especially suitable for environmental monitoring purposes, where we have seasonal trends that we must remove.Contrary to what the name may suggest, the cumulative sum is not the cumulative sum of the values, but the cumulative sum of the differences between the values and the mean.This is: where N is the total number of temporal samples, Ii represents a sample take at time i and E[I] is the expected value of the samples over a training time frame.The expected value E[] is approximated by using the sampling average and assuming a large enough sample size.

CUSU-SMC
In this section we present the modification to the CUSUM algorithm to make it more suitable for work with SAR data, which can show seasonal trends.We call this proposed change detector CUSU-SMC (cumulative sum-spatial mean corrected).
As the main goal of this research focuses on detecting the forest changes that occurred from the beginning of 2019 onwards, we used the information compiled in all 2018 images as control samples.As is shown later, the selection of the mean plays an important role.We used spatio-temporal datasets, where data are stored in cubes with spatial dimensions (range and azimuth) and a temporal dimension.Accordingly, we can calculate the mean value E[I] in the spatial domain or the temporal domain.
1. Temporal mean: previous approaches, followed by [43,44], calculate the average backscatter value for each pixel, using all the images acquired over the reference time, in this case, the year 2018.
where Nt is the number of images in the reference time.
This will generate a single mean value  ̅ for each specific pixel.This approach presents certain theoretical disadvantages for dense time series over vegetation.A unique mean value for each pixel does not take into account any temporal variation that might be associated with the seasonal change in vegetation (e.g.phenology variations) or meteorological events (e.g., precipitation or temperature).We understand that seasonal variations (due for instance to phenology) represent genuine changes to the forest canopy.Additionally, a drought or a rain event which produces changes all over the forested areas may be misinterpreted as structural changes.
2. Spatial mean: To minimize the possible influence that seasonal variation may have on the detection of forest structural changes, we used a spatial mean (instead of a temporal one) which follows the seasonal variation of forest ecosystems.This approach uses the mean estimated over the forest areas at each specific date in the timeseries.
where Ns is the number of pixels included in forest areas and Ii are the pixels in forested areas.In other words, each time step has a characteristic  ̅ .
To do this, we needed to use a forest inventory or produce a rough forest classification mask using SAR or optical images.

CUSU-SMC Implementation
The implementation was done in Python, downloading data using the SEDAS Sentinel Data Service Access Hub.
Data cubes were first created using the Sentinel-1 images.Three data cubes were created, using the polarization channels VV, VH, and their ratio VH/VV.The cubes were then overlaid with forest inventories masks.The forest mask was previously extracted from the 2018 Scottish National Forest Inventory [41].Forest can be classified differently based on different precision ranges.In this case, it was decided that forest must contain indistinct conifer and broadleaf masses, no matter if they are pure or mixed stands.Once the forest cube is calculated, it is possible to get the spatial forest mean (SFM) for each specific date and implement it into the CUSU-SMC analysis.
In terms of numerical implementation, the CUSU is calculated with the following iteration: 2) Spatial mean | CUSU-SMC where  = 0 and i is an iteration over the temporal axis of the cube.Nt is the total number of acquisitions used in the training dataset.The image I can be the VV, VH polarization channels or their ratio.The comparison of results using the three observables is provided in the following.

Z Score Calculation
The result of the CUSUMs (refers to both CUSUM and CUSU-SMC approaches) is either a cumulate intensity or ratio.In order to produce a detector a threshold value was needed for both CUSUM algorithms..In this work we proposed a rigorous statistical framework to set thresholds for both CUSUM algorithms.
The first step was to understand which type of distribution the CUSUMs have.In the following we will show that after all the processing the CUSUM in the training dataset was distributed as a Gaussian zero mean with standard deviation .This allowed us to use a Z test for setting thresholds in a statistically rigorous way, with some significance value attached to the results.The first step was to calculate the standard score, which will provide information about the probability of a point having a set Z score.The Z score needs knowledge of the population mean  ̅ and standard deviation  : where, σ is calculated using all the images acquired during the reference year (2018, with 54 images in our case). is the cumulative sum at time i and  is the standard deviation of  at the time i.
In the future, we will explore the use of alternative methodologies such as autoregression models associated with the output of the CUSUMs.
The null hypothesis is that ʺno changeʺ occurred.This is the hypothesis where the  , under test belongs to the same population as the training dataset.The alternative hypothesis is that a change happened, or more formally the  under test does not belong to the same population as the training dataset:

H1 = Change
In the test the p-value is the probability that the CUSUM sample belongs to the null hypothesis (H0).This means that a small p-value will suggest that the H0 is rather unlikely and therefore we can reject H0 and argue that we have a large confidence that H1 is true.

Logging Area Detection Accuracy Assessment
In this work, the accuracy of the algorithm was assessed through an error matrix (or confusion matrix), an approach which evaluated the quality of the classifier, based on the correspondence between the classification results and a reference image (Table 1).From the confusion matrix, overall accuracy (OA, Equation ( 7)), producer accuracy (PA, Equation ( 8)), user accuracy (UA, Equation ( 9)) and F-measure (F-score, Equation (10)) were individually calculated for each data cube.F-measure is one of the best measures for evaluating the detector performance since it summarizes the balance between (PA) and (OA), prioritizing the identification of real detections over the detection of true negatives [48,49].

Post-Processing and Sieve Filtering
The high sensitivity of electromagnetic signals to many variations occurring in target parameters (e.g., geometry, moisture contents, or surface roughness) makes SAR technology a valuable tool for detecting changes.However, this sensitivity also leads to significant misclassifications and estimation errors when performing change detection analyses.Tackling high false alarm rates is a common issue when working with satellite imagery.A methodology that can be used at the expense of spatial resolution is post-processing with morphological filters, such as the sieve filter [50].The sieve filter removes all polygons smaller than a given size and merges them with their largest neighbors [50,51].Beside the threshold on the number of connected pixels, a sieve filter differentiates among two types of pixel connectivity.A type 4 connectivity is established for adjacent pixels that share a side, while a pixel with a connectivity of 8 shares any side or a vertex with its adjacent [51].
The filtering impact on classification accuracy was evaluated using three different thresholds of 10, 15, and 25 pixels; with a connectivity of type 8, as this is less restrictive in terms of shape.The analysis was performed using the geospatial tools of ArcMap 10.5 (Int, Region group), and QGIS software (Gdal-Sieve filter).The final accuracy assessment was carried out in the Python environment, making use of the open source Geospatial Data Abstraction Library (GDAL).

CUSUMs Framework
All Sentinel-1 GRD images acquired for our period of study considered the VV and VH polarization channels.Channels were calibrated as normalised radar cross sections.We used these to compute the ratio of VH over VV.The timeseries was stored in 3 data cubes, one for VV, one for VH, and one for ratio.In the following step, we reduced noise by applying a boxcar filter (3 x 3 window size).The filtering helped reduce speckle variation.Subsequently, unitless backscatter values were converted to dB using a logarithmic transformation.This pre-processing was identical for CUSUM and CUSU-SMC.
where DN is digital number.
The following step marks the separation between the CUSUM and CUSU-SMC approaches.CUSUM uses a unique mean value for each pixel calculated over the reference timeseries (year 2018).CUSU-SMC uses a unique mean for each acquisition, and this mean is updated for each acquisition.Z-score calculation allowed us to evaluate the probability of change, rejecting or accepting our hypothesis under different levels of significance.Change detection maps were subject to a first accuracy assessment, making use of reference logging activity masks, provided by the FLS as the validation dataset.Finally, a sieve morphological filter, with multiple windows sizes, was applied to reduce the number of random flagged pixels, and the accuracy was assessed again to evaluate the effects of the morphological filter on the final result (Figure 3).

Trends and Biases of CUSUMs
In this section we want to validate the CUSUM [46] and compare this with the new CUSU-SMC.As explained, the main difference is in the use of a spatial mean instead of a temporal one.To demonstrate the results, in the following we will plot the CUSUM and CUSU-SMC using the polarization ratio only.Other polarization channels show similar behavior.Figure 4a presents the CUSUM over three plots for three classes, "change", "conifer", and "broadleaves".In this test the training was done using the 2018 data (54 dates).As can be seen in Figure 4a, all the lines seem to converge to zero after exactly one year and then diverge afterwards.This is due to a mathematical constraint that forces a value of zero when the CUSUM is calculated for exactly 54 dates.The equation used to calculate the CUSUM is linear and we get zero when the added part equals the mean: As a comparison Figure 4b presents the proposed CUSU-SMC over the same regions for the three classes, "change", "conifer", and "broadleaves".This figure shows that CUSU-SMC seems to slowly, but constantly, diverge for each of the classes.This observation can be made for almost all the pixels in the image.This is because some pixels will naturally have a backscattering which is higher (or lower) than the mean  ̅ due to different forest structures.The constant change will act as a bias and accumulate over time forming a ramp.Since most of the pixels present different forest structures they will all show a ramp.Since the differences act as a constant bias the trend is approximately linear (i.e., this is the function with a constant gradient).We can easily estimate the line of best fit for the temporal trend during the training dataset and remove it for the rest of the time series.After removal, the CUSU-SM pixels will be fluctuating around 0, unless there is a change.Figure 4c shows the result of CUSU-SMC (after correction), where the ramp did not remove the change because it is estimated over the training dataset where no changes occurred.Moreover,It should be noted that the corrected trends are not constrained to be zero at the end of the training interval.The change that they show is what was actually accumulated compared to the rest of the forest.
One rational for using the spatial mean is that in theory this will help remove anomalies due to meteorological conditions or forest phenology.Since the weather and phenology will affect all the surrounding forest in a roughly equal way, this will modify the spatial mean that we calculate and therefore the weather signal would be removed.This procedure may be seen as an attempt to calibrate the data based on the weather condition.Of course, there will still remain some random fluctuation, due to the fact that some forest species may not change backscattering in the same way compared to other forest stands.However, this should be a second order error, that could be eventually corrected if we have, either an accurate forest inventory or a good model to predict more accurately the change of the forest behavior due to weather conditions.

Testing the Assumptions
In this section we want to test the assumption that the CUSUMs obey a Gaussian zero mean distribution.We decided to determine normality by investigating graphically the output of data histograms and quantile-quantile plots.The results obtained for each of the three classes "forest", "conifer", and "broadleaves" provided enough evidences to ensure the data is normally distributed.As can be seen in Figure 5, histograms highlighted a close approximation of the data to the theoretical bell curve distribution; while in the Q-Q plots data points appeared closely distributed along the diagonal line.The data used for generating the histograms was obtained by masking out each region ("forest", "conifer", "broadleaves") over a CUSU-SMC ratio image (29.07.2019).Data used for Q-Q plots correspond to temporal values of unique pixels, randomly selected, for each region.

CUSUM Logging Activities Detection
The results obtained from the accuracy assessment showed an optimal change detection performance, for both the CUSUM and the CUSU-SMC.The accuracy assessment was carried out by using multiple significance levels in order to test the performance of the methods.Despite different results being obtained for each level of significance, we decided to set alpha =0.05 as the most appropriate value for comparing the tested products.Results are shown in Table 2 from which we can conclude:

1) Best polarization channel:
Assessing the best polarimetric channel to monitor logging for temperate forest is crucial since previous studies have shown some discrepancy regarding this.Previous work [16,52,53] suggested that co-polarized backscatter (HH, VV) was more sensitive to forest gaps; with other work [44,54] demonstrating higher sensitivity of cross-polarized (VH, HV) for forest structural changes, as a result of the disruption of the high backscatter values associated with the complex geometries of tree canopies.
The temporal trend for VV, VH, and the ratio are showed in Figure 6.Data were filtered with a 3 x 3 boxcar filter.The last row of Figure 6 shows the backscatter differences between the stable forest areas and the logging areas for the VV, VH and its ratio.The CUSUM method [46], showed the best results for the VH polarization channel (OA = 73.6%;F = 69.5%),providing a medium to high performance in the detection of true positives (TP = 60.2%) and a significantly high performance in the rejection of false positives (FP = 13.0%).Similarly, the proposed CUSU-SMC also showed its higher performance with the VH polarization channel (OA = 65.7% and F = 68.0%),offering a high performance in detecting true positives (TP = 72.4%),but a relatively high false alarm rate (FP = 40.6%).The ratio (VH/VV) combination provided a good performance for both CUSUM (F = 68.7%)and CUSU-SMC (F = 67.4%)methods.Finally, the impact of ground on the co-polarized channel (VV) provided the lowest separation between forest and deforested areas.
2) Detection strategies: Despite the similar overall performance observed for both CUSUM and CUSU-SMC, we found a notable difference in terms of true positives and false alarms.Specifically, the CUSU-SMC showed a much higher power of the statistical test, which means it is much more capable of detecting true positives, but it also had an increase of true negatives.Given that the total OA is similar, we can identify two main detection strategies as: a) conservative strategy (*C); this leads to low false positive rate, but also lower true positive; b) tolerant strategy (*T); this leads to higher true detection rates at the expense of having more false alarms.
In Table 2 we identified the algorithm and polarization channel that would lead to best results for one of the two strategies.
Table 2. CUSUM and CUSU-SMC accuracy assessment.Accuracy rates were calculated using 2019 logging areas as reference, the same acquisition date (29 September 2019) and a significance level of (α = 0.05).TP = true positives, TN = true negatives, FP = false positives, FN = false negatives, OA = overall accuracy, PA = producer accuracy, UA = user accuracy.(*C = conservative approach, *T = tolerant approach, and *NC = results obtained without the application of the fitted ramp correction).The correction for the constant bias applied in the CUSU-SMC approach, provided a significant improvement in the detection performance (OA = +5%, F= +4.3%), as a result of a considerable reduction of the false alarms rates (FP = -6.7%)(see Table 2, *NC).It is therefore important to apply this when using the CUSU-SMC.

CUSUM
The similar F-score obtained for CUSUM and CUSU-SMC (F-score differences of 2.1%) may inform us that correcting for seasonal trends does not improve the detection.This could be explained by the stable behavior of the radar backscatter signal observed for forest areas.This behavior was also observed by Quegan et al. [55], who found that old forest exhibits a lower seasonal variability than younger forests.To check this point in Figure 7, we plotted the single temporal mean (CUSUM method) and the spatial forest mean (CUSU-SMC method) as time trends, where the single temporal mean was a constant in the plot.The different curves represent conifer and broadleaf forest stands.It can be seen how conifer and broadleaf plots present opposite behaviors, in terms of time trends.The spatial forest mean averaged over all the trends, and therefore these temporal trends are averaged out.If the forest inventory is large enough, calculating the spatial forest mean over different forest types may improve the final detection.

Comparison with Pairwise Analysis
In order to compare the CUSUMs performance against pairwise combination we considered two different approaches: 1. Continuous pairwise.This approach aims to generate a continuous change analysis based on the comparison of each image present in the timeseries to the image acquired three times before.
where, i is an iteration over the temporal axis of the cube.Nt is the total number of acquisitions used in the training dataset.The image I used for the analysis corresponds to a channel, or the ratio after applying a boxcar filter (window size = 3).
2. Static pairwise.This follows the common direct pairwise change detection, by which two acquisitions separated in time by a specific event (clear-cut activities in our study case) are directly compared.Considering that all logging activity studies were carried out between 1 December 2018 and 20 June 2019, we have used the image acquired on 3 November 2018 as the reference image (pre-event), and the last image of the timeseries acquired on 29 September 2019 as the change image (post-event).
Table 3. Pairwise vs CUSUMs detection performance comparison.The change detection accuracy assessment was performed using a significance level of (α = 0.05) and the same image (29 September 2019) for both the CUSUMs and pairwise approaches.The CUSUMs approaches were the ones for conservative (*C) and tolerant (*T) strategies.TP = true positives, TN = true negatives, FP = false positives, FN = false negatives, OA = overall accuracy, PA = producer accuracy, UA = user accuracy.The results obtained for both pairwise methods are shown in Table 3.It is clear thatthe pairwise method had a significantly lower performance in the detection of logging activities when compared to the proposed CUSUM and CUSU-SMC methods.The low detection performance may have been due to the high signal variation that SAR images may have on single dates.Minor changes to weather variables (e.g., rainfall, temperature, wind) or to vegetation phenology may lead to significant misclassification errors.

ROC Curves
Although the Gaussian model fits very well the CUSUM values for no-change, it is still important to check that it is not the statistical model adopted that is producing the different performance, but the framework to process the SAR data.We therefore analysed the receiver operating characteristic (ROC) curves [56], which are independent of the selection of the threshold (and therefore the significance level).Figure 8 presents the ROC curves of both continuous and static pairwise approaches, CUSUM, and CUSU-SMC using VH polarization and ratio (VV has not been shown due to its significantly lower performance).ROC curves plot the true positive rate against the false-positive rate while the threshold is varied and have been widely used for change detection performance comparisons in multiple scientific disciplines [57].The closer the ROC curve is to the upper left corner, the higher the overall accuracy of the test [58].In our test, detection performance was tested over multiple significance levels (0.4, 0.3, 0.2, 0.15, 0.1, 0.05, 0.01, 0.001).A higher detection performance was obtained by both CUSUM and CUSU-SMC, which show a similar pattern.Using the ratio, the pairwise approaches always show lower detection accuracy, while using the VH channel the pairwise gives similar performance when the probabilities of false alarms is very high (higher than 0.5).

Final Detection Mask
The differences between CUSUM and pairwise are also clearly noticeable in the final detection maps.Since the detector we developed is statistical, we can also attach probabilities to the change.This therefore helps evaluate how confident we are that a change is genuine.
Figure 9 compares the change detection map of static pairwise against the CUSU-SMC method for two study areas affected by clear-cuts.The change map for the CUSUM is very similar to the CUSU-SMC and therefore has not been presented here.As can be seen, CUSU-SMC offered the best performance, detecting almost the entire affected area in a homogeneous way, with only a lower detection capability observed at the edges.On the contrary, the pairwise analysis identified a significantly smaller area of change, being those changed areas detected in a scattered way.Additionally, the pairwise approach showed overall lower probability values, a fact which manifests the lower power of the test.

Post-Processing Sieve Filter
Table 4 shows that the sieve filter produces a significant increase in the detection accuracy as a result of a reduction of false alarms.This spatial filter added between 3.6% to 5% to the overall accuracy and 3.5% to 7% to the F-score; this improvement being greater with a smaller filter size of 10 connected pixels.Contrary to Rüetschi et al. [16], our results suggest that the use of wider windows would result in a more inefficient performance of the detector, since, considering the accuracy/resolution trade-off, the slight accuracy improvements would not compensate for the loss in spatial resolution.We believe that the best filter size is a function of plot size and logging extent, and therefore needs to be set using knowledge of the area.The sieve filter also provided a considerable improvement in visual interpretation, since cleared areas could be distinguished more readily due to the high, and continuous, aggregation of ʹchangeʹ pixels and the significant enhancement of the polygon edges.Finally, Figure 10 shows the influence of the sieve filter on the final change detection accuracy of both the CUSUM and CUSU-SMC methods using ROCs.

Forest Logging Detection Maps with Sieve Filter
In Figure 11, the final change detection maps show the capabilities of the cumulative sum to detect logging activities.Despite the high presence of isolated false positives in forest areas, all clearcut areas are distinguishable.The implementation of the pixel connectivity sieve filter significantly helped to remove those flagged pixels, which were likely caused by environmental factors (e.g., meteorological condition variation).This step contributed to the refinement of the final product by generating a better representation of the disturbed areas, while improving the overall detection accuracy rate.

Discussion
The results of this study indicate that the proposed method, which is based on the analysis of dense SAR timeseries by the cumulative sum test, detects forest logging activities in a continuous way, with similar spatial accuracy to previous forest change detectors [16,31,54,59].The analysis of SAR temporal backscatter signal showed significant differences between felled areas and intact forest areas.Despite these differences, observed in both polarizations, the cross-polarization (VH) channel showed the higher sensitivity to forest structural variations, providing the best change detection performances.From the various parameters and cumulative sum approaches tested, the combination of the CUSUM method with a 3 x 3 mean filter generated the best results, offering an OA = 73.6% and F = 69.5%.The subsequently added sieve-pixel connectivity filter (10 pixels) increased the detection accuracy values to UA = 85.9%,OA = 77.2, and F = 74.0.

Near Real Time
The proposed method emerges as a valuable forest monitoring tool when it is compared to the existing forest change detector products.The exploitation of dense SAR timeseries for the CUSUMs allows the monitoring of forest structural disturbances in a rapid, continuous manner, which contrasts with the current monitoring methods based on optical systems.Both CUSUM and CUSU-SMC methods allow the gathering of information about changes that extends over long periods of time.This may be significantly helpful when analysing diverse patterns and trajectories associated with those disturbances which can be crucial to forecast future alterations, as in the case of illegal deforestation, agriculture expansion, or risk management.When developing a near real-time (NRT) monitoring tool, low latency is considered as the driving factor for obtaining high temporal accuracies [16,60].This stands out as the main disadvantage associated with the use of spaceborne optical sensors, since data acquisition is highly limited by cloud coverage.The CUSUM methods have been tested using a repeat-track orbit and single acquisition mode (ascending).For the 22-month study period (1 January 2019 -30 September 2020), a total of 84 observations were used (~3.8 img/month).Despite the Sentinel-1 image acquisition strategy allowing a lower latency for the study area (~1 img/3 days), the irregular topography of the region hindered the use of a multi-aspect approach.The utilization of multiple acquisition modes would generate a significant increment in the false alarm rate, caused by the geometrical and signal divergences between different viewing angles.The implementation of CUSU-SMC in flatter areas, would allow the combined use of tracks and orbits (ascending/descending).This approach needs to be further explored in order to assess the contribution of the multi-aspect approach to time detection n performance.

Adaptability to Other Sensors and Applications (Synergies)
The CUSUM-SMC method enhances the capability of the SAR system for monitoring environmental disturbances over natural ecosystems, thus standing out as a reliable alternative to the current optical-based deforestation monitoring methods [31].The single requirement of having access to a continuous data acquisition program gives CUSUMs approaches the ability to be easily adaptable to other SAR missions and applications.Apart from the ensured continuity of the Sentinel-1 mission until 2030, there are other SAR missions whose image acquisition programs would allow the implementation of this method.The C-band RADARSAT Constellation Mission (RCM), operational from November 2019 and with a lifespan of 7 years, will provide additional C-band data with a revisit time of 4-days [61,62], a fact which will contribute to shorten the times between acquisitions [16].Moreover, the recent announcement of the future free-data policy of ALOS-2/(ScanSAR), the upcoming ALOS-4 L-band (launch date: 2021), and NASA-ISRO (NISAR) L,S-band (launch date: 2022) SAR missions might serve to enhance the performance of the CUSUMs, due to the higher sensitivity that longer wavelengths have to structural forest disturbances [17].
With respect to adaptability to other environmental monitoring applications, this method shows high versatility.The ability of CUSUMs to detect disturbances in a continuous manner makes it potentially suitable for other applications (e.g., flood mapping, coastal erosion, land degradation) characterized by high dynamicity in the environment.

Limitations and Future Work
The area of study used for this work may have led to underestimations in terms of the potential capabilities of the algorithm.Aspects such as the irregular topography of the region or the low backscatter variability, showed by old forests throughout the whole timeseries, may have led to lower accuracies in the detection of forest structural disturbances.
In the future we are planning to test the algorithm more globally, over very different forest types and conditions.To achieve this in an efficient way we are working on the transcription of the code to the Google Earth Engine (GEE) platform.GEE will allow us to test the CUSU-SMC over other forest regions around the globe, assessing the influence of the diverse environmental conditions and minimizing the pre-processing times.The resultant GEE product will provide a worldwide, free, and time efficient forest monitoring tool accessible to foresters.
Additionally, we will try some more complex statistical methodologies to apply to the CUSUMs outputs, such as autoregression models [63].
Another very interesting future work is to try to relate the forest change to diverse physical parameters of the forests, such as leaf area index [64] or biomass, thus exploring the influence of forest variability on the CUSUMs detection performance [65].

Conclusions
This study introduced a new forest monitoring method for detecting forest disturbances from the analysis of dense SAR time series.Prolonged differences in the radar backscatter signal (before and after logging activities) between stable and disturbed forest areas highlighted the sensitivity of C-band SAR systems to forest structural changes.The high spatial and temporal resolution, global coverage, and free data policy of the ESA-Sentinel 1 mission has allowed the exploitation of dense SAR time-series for environmental monitoring.The cumulative sum method makes use of dense time-series, allowing the detection of forest changes in a continuous and rapid manner.The main outcome of this study was the demonstration of the efficient detection capabilities of SAR systems over other more expensive (e.g., LiDAR , manual surveys) or weather dependent monitoring systems (e.g., optical).
The method presented was tested for logging activities in Scottish forests, providing faster detection rates (early warning/rapid response) than current optical methodologies used for the national forest inventory.With the single requirement of having access to continuous data acquisition, the proposed method can be easily implemented on other SAR sensors such as the recently launched RCM constellation or the ALOS-Palsar mission (open policy recently agreed), and on other global forests.CUSU-SMC stands out as a valuable monitoring tool for foresters/forest managers because it reduces the current response times, while also cuts the costs and efforts, when evaluating forest damage.

Figure 1 .
Figure 1.Reference map showing the location and areas covered by the Queen Elizabeth Forest Park.

Figure 2 .
Figure 2. Forest logging activity identification.(a) Landsat 8 shortwave infrared image used for the visual identification of forest changes.Acquired: 28 June 2019; pink tone areas represent land cover transformation from forest to bare ground.(b) Colored polygons represent Forest and Land Scotland (FLS)-verified forest logging activities carried out between 1 January 2019 and 30 September 2019.

Figure 3 .
Figure 3. Proposed methods for the cumulative sum (CUSUM) and cumulative sum-spatial mean corrected (CUSU-SMC) forest change detectors.

Figure 5 .
Figure 5. Assessing if CUSUM and CUSU-SMC are normally distributed over the reference year.The data used for generating the histograms was obtained by masking out each region ("forest", "conifer", "broadleaves") over a CUSU-SMC ratio image (29.07.2019).Data used for Q-Q plots correspond to temporal values of unique pixels, randomly selected, for each region.

Figure 6 .
Figure 6.This figure illustrates the sensitivity of different polarization channels to forest structural changes by comparing the annual backscatter timeseries of logging areas with undisturbed forest.Temporal differences extracted from the comparison of logging areas and undisturbed forest (bottom) showed the higher sensitivity of the VH polarization channel.

Figure 7 .
Figure 7. Synthetic aperture radar (SAR) backscatter signal comparison for mean reference values used for the CUSUM and CUSU-SMC change detection methods.The lack of significant differences in the values obtained for the forest mask for both cases tested could explain the similarity in the final results.

Figure 8 .
Figure 8.This figure compares the performance of the pairwise approaches with the CUSUM and CUSU-SMC methods using a ROC curve representation, for both the Ratio (VH/VV) and VH channels.PW = pairwise.ROC curves were generated using the results obtained for the accuracy assessment, performed for the image acquired on 29 September 2019, using significance levels of (α=0.4,0.3, 0.2, 0.15, 0.1, 0.05, 0.01, 0.001).

Figure 9 .
Figure 9. Logging activities change detection maps obtained for pairwise (static approach) (left) and CUSU-SMC (right) methods.The images show the probability of change (being P = 1 considered as change, and P = 0 as no change), obtained for the 29 September 2019 and ratio data cube, for the areas affected by clear-cuts ʺUnknown 4ʺ (top), and ʺ33020ʺ (bottom).

Figure 10 .
Figure 10.Contribution of the tested sieve filter to detection performance of the CUMSUM and CUSU-SMC approaches.Accuracy test carried out using a significance level of ( 0.05 .

Figure 11 .
Figure 11.Forest changes detected in the Queen Elizabeth Forest Park using the conservative strategy (CUSUM VH).(a) Shows the logging activities reference dataset obtained for the period of study.(b) Illustrates the changes detected (red) for the CUMSUM VH (avg 3 x 3) detector.Subsets (c) and (d) show a closer look of the effect of the sieve filter on the final change detection maps [(c) 10 pixels; (d) 25 pixels].

Table 4 .
Results showing the spatial filtering effects on detection performance/accuracy.