Detecting Forest Changes Using Dense Landsat 8 and Sentinel-1 Time Series Data in Tropical Seasonal Forests

: The accurate and timely detection of forest disturbances can provide valuable information for e ﬀ ective forest management. Combining dense time series observations from optical and synthetic aperture radar satellites has the potential to improve large-area forest monitoring. For various disturbances, machine learning algorithms might accurately characterize forest changes. However, there is limited knowledge especially on the use of machine learning algorithms to detect forest disturbances through hybrid approaches that combine di ﬀ erent data sources. This study investigated the use of dense Landsat 8 and Sentinel-1 time series data for detecting disturbances in tropical seasonal forests based on a machine learning algorithm. The random forest algorithm was used to predict the disturbance probability of each Landsat 8 and Sentinel-1 observation using variables derived from a harmonic regression model, which characterized seasonality and disturbance-related changes. The time series disturbance probabilities of both sensors were then combined to detect forest disturbances in each pixel. The results showed that the combination of Landsat 8 and Sentinel-1 achieved an overall accuracy of 83.6% for disturbance detection, which was higher than the disturbance detection using only Landsat 8 (78.3%) or Sentinel-1 (75.5%). Additionally, more timely disturbance detection was achieved by combining Landsat 8 and Sentinel-1. Small-scale disturbances caused by logging led to large omissions of disturbances; however, other disturbances were detected with relatively high accuracy. Although disturbance detection using only Sentinel-1 data had low accuracy in this study, the combination with Landsat 8 data improved the accuracy of detection, indicating the value of dense Landsat 8 and Sentinel-1 time series data for timely and accurate disturbance detection.


Introduction
Understanding forest dynamics provides valuable knowledge for effective forest management [1]. Remotely sensed data are widely used to characterize these dynamics. Optical satellite remote sensing is suitable for capturing current and past forest conditions over large areas because of its wide spatial coverage and frequent temporal observations [2]. To detect small-to large-scale disturbances accurately, medium spatial resolution satellite data have several desirable characteristics [3]. Recently, there has been an increase in freely available and well-calibrated satellite data (e.g.,  thanks to open data policies [4][5][6][7][8]. This situation has encouraged researchers to process data from numerous satellites to extract temporally and spatially improved forest change information [9][10][11].
Many approaches have been proposed to detect forest changes using time series medium spatial resolution satellite data. Landsat time series have been particularly used for change detection in the last combines the time series probability of forest changes predicted independently from optical and SAR models. This approach has the advantage of analyzing different sensors separately, thus preserving their sensor-specific characteristics [35]. In this context, Reiche et al. [22] developed a Bayesian approach to detect deforestation events in the tropics using Landsat and ALOS/PALSAR time series. They enhanced the conditional probabilities of non-forest using probability density functions for each sensor, and then updated the conditional probability of deforestation using both sensors based on Bayes' theorem. Reiche et al. [35] also investigated the use of Landsat, Sentinel-1 and ALOS-2/PALSAR-2 time series data in a similar probabilistic approach. Verhegghen et al. [48] simply combined the burned area mappings independently generated from Sentinel-1 and Sentinel-2.
Although several studies have investigated the combination of optical and SAR time series data, there have been relatively few studies that only used freely available satellite data (e.g., Landsat and Sentinel-1). Furthermore, even though previous studies have successfully detected forest changes, the probability of forest and non-forest mainly depended on parametric models (e.g., Gaussian models). Non-parametric models, such as machine learning, offer an alternative, as they have the ability to handle complex data [52,53] and have been widely used for forest monitoring [47]. The use of machine learning algorithms might be effective for predicting forest and non-forest probabilities in regions with various causes of forest change. However, machine learning approaches for dense time series optical and SAR data are not well developed.
The aim of this study was to evaluate the applicability of a disturbance detection approach using dense time series observations of Landsat 8 and Sentinel-1 based on a machine learning algorithm in the tropical seasonal forests of Myanmar. Disturbance probability models were built for each observation (i.e., observed pixel values at the time of data acquisition) of Landsat 8 and Sentinel-1 using a Random Forest (RF) algorithm. The predicted time series of disturbance probabilities were then used to detect forest disturbances in each pixel location. We used a simple approach to combine the disturbance probabilities of both sensors, and assessed the time needed to detect disturbances and the accuracy of disturbance detection.

Study Area
The Bago region is located in central Myanmar (Figure 1). A total of 2.3 million ha ranging from 10-600 m elevation was analyzed in this study. The study area is mainly composed of tropical mixed deciduous forests. The mean annual precipitation is approximately 2000 mm and the mean temperature is 27 • C. There are three seasons in this region: a hot season from February through May, a rainy season from May to October, and a dry season from November to January. There are several kinds of human-induced forest disturbances in this region. Logging activities (both legal and illegal) are the main cause of disturbances in the mountainous region [55,56]. Clearing for the development of timber plantations is also observed in the mountainous region. In areas near human settlements, shifting cultivation practices and conversion to permanent agriculture are major disturbance agents, especially in the eastern part of the study area. Other disturbance agents include water inundation and infrastructure development.

Data and Methods
An overall analysis of the disturbance detection process is presented in Figure 2. We also show the flowchart of overall analysis in Figure 3 to present training and test samples used for both RF disturbance probability models and disturbance detection. First, Landsat spectral indices and SAR backscatter coefficients were derived from time series observations of Landsat 8 and Sentinel-1. A harmonic regression model was fitted to represent forest structural change and vegetation seasonality in each spectral index and backscatter coefficient. RF disturbance probability models were then independently built for Landsat 8 and Sentinel-1 using the results of harmonic regression. Time series of disturbance probabilities were predicted for a two-year period, and disturbances were detected based on the combined disturbance probabilities. We assessed the accuracy of disturbance detection using only Landsat 8, using only Sentinel-1 and using the combination of both sensors. There are several kinds of human-induced forest disturbances in this region. Logging activities (both legal and illegal) are the main cause of disturbances in the mountainous region [55,56]. Clearing for the development of timber plantations is also observed in the mountainous region. In areas near human settlements, shifting cultivation practices and conversion to permanent agriculture are major disturbance agents, especially in the eastern part of the study area. Other disturbance agents include water inundation and infrastructure development.

Data and Methods
An overall analysis of the disturbance detection process is presented in Figure 2. We also show the flowchart of overall analysis in Figure 3 to present training and test samples used for both RF disturbance probability models and disturbance detection. First, Landsat spectral indices and SAR backscatter coefficients were derived from time series observations of Landsat 8 and Sentinel-1. A harmonic regression model was fitted to represent forest structural change and vegetation seasonality in each spectral index and backscatter coefficient. RF disturbance probability models were then independently built for Landsat 8 and Sentinel-1 using the results of harmonic regression. Time series of disturbance probabilities were predicted for a two-year period, and disturbances were detected based on the combined disturbance probabilities. We assessed the accuracy of disturbance detection using only Landsat 8, using only Sentinel-1 and using the combination of both sensors.   Disturbances were detected using disturbance probabilities predicted by the RF models.

Preprocessing of Landsat 8, Sentinel-1 and Topographic Data
We used the Landsat 8 Operational Land Imager (OLI) surface reflectance data available in the Google Earth Engine (GEE) [58] cloud computing platform. We only used data located in our study area (i.e., path/rows 132/47, 132/48, 133/46, 133/47 and 133/48 in World Reference System 2). All data with cloud cover of less than 70% acquired from the beginning of January 2014 to the end of December 2018 were used in the analysis. A total of 374 data were retained. The Landsat 8 surface reflectance data were atmospherically corrected using LaSRC [59]. Clouds and cloud shadows were masked using a pixel quality assessment band generated from the CFmask algorithm [60].
The Landsat 8 surface reflectance data were transformed to five spectral indices to represent the spectral response of vegetation conditions. We used the Tasseled Cap indices, which included the Tasseled Cap Brightness (TCB), Greenness (TCG), Wetness (TCW) and Angle (TCA) [61] and the Normalized Burn Ratio (NBR) [62]. A total of five time series of Landsat 8 spectral indices were developed.
Sentinel-1 data originally from Level-1 Ground Range Detected (GRD) products acquired in Interferometric Wide Swath (IW) mode from ascending and descending orbits were used in this study. We used VV and VH dual-polarized data acquired from the beginning of January 2014 to the end of December 2018, which are available in GEE. Sentinel-1A and 1B were launched in April 2014 and April 2016, respectively. As a result, the earliest scene in this study was acquired in October 2014. A total of 583 scenes were used for analysis. In GEE, Sentinel-1 GRD data have been preprocessed to backscatter coefficients (σ°) using the following steps: thermal noise removal, radiometric calibration and terrain correction with SNAP Sentinel-1 Toolbox [63]. Additionally, speckle filtering was

Preprocessing of Landsat 8, Sentinel-1 and Topographic Data
We used the Landsat 8 Operational Land Imager (OLI) surface reflectance data available in the Google Earth Engine (GEE) [58] cloud computing platform. We only used data located in our study area (i.e., path/rows 132/47, 132/48, 133/46, 133/47 and 133/48 in World Reference System 2). All data with cloud cover of less than 70% acquired from the beginning of January 2014 to the end of December 2018 were used in the analysis. A total of 374 data were retained. The Landsat 8 surface reflectance data were atmospherically corrected using LaSRC [59]. Clouds and cloud shadows were masked using a pixel quality assessment band generated from the CFmask algorithm [60].
The Landsat 8 surface reflectance data were transformed to five spectral indices to represent the spectral response of vegetation conditions. We used the Tasseled Cap indices, which included the Tasseled Cap Brightness (TCB), Greenness (TCG), Wetness (TCW) and Angle (TCA) [61] and the Normalized Burn Ratio (NBR) [62]. A total of five time series of Landsat 8 spectral indices were developed.
Sentinel-1 data originally from Level-1 Ground Range Detected (GRD) products acquired in Interferometric Wide Swath (IW) mode from ascending and descending orbits were used in this study. We used VV and VH dual-polarized data acquired from the beginning of January 2014 to the end of December 2018, which are available in GEE. Sentinel-1A and 1B were launched in April 2014 and April 2016, respectively. As a result, the earliest scene in this study was acquired in October 2014. A total of 583 scenes were used for analysis. In GEE, Sentinel-1 GRD data have been preprocessed to backscatter coefficients (σ • ) using the following steps: thermal noise removal, radiometric calibration and terrain correction with SNAP Sentinel-1 Toolbox [63]. Additionally, speckle filtering was performed to reduce speckle noise using a refined Lee filter [64] with a window size of 7 × 7. The preprocessed Sentinel-1 data were resampled to 30 m spatial resolution.
To reduce the possible effects of variations in the incident angle of Sentinel-1 backscatter coefficients [36,65], we only used observations that had almost the same incident angle at the same pixel location. Because all observations were recorded as having incident angles below 37 • or above 39 • in this study, we labeled the incident angle of each observation as "low" or "high." Only observations labeled with the majority angle class were used in the same pixel location. By doing this, the time series of backscatter coefficients were less affected by variations in the incident angle. We used observations from both ascending and descending orbits in the same time series in subsequent analysis.
The digital elevation data were obtained from the Shuttle Radar Topography Mission (SRTM) version 3 [66] in GEE. This data has 30 m spatial resolution and has been processed with a void-filling algorithm. We used SRTM data to calculate the topographic slope and elevation.

RF Disturbance Probability Models for Landsat 8 and Sentinel-1
To derive the disturbance probability of each observation in Landsat 8 and Sentinel-1, time series of the Landsat 8 spectral indices and Sentinel-1 backscatter coefficients were respectively used in the RF models. The signals of both Landsat 8 and Sentinel-1 time series data are affected by the seasonality of soil and vegetation moisture content, as well as forest structure changes. There are two main approaches to cope with seasonality: removing seasonality from signals (e.g., References [20,35]) and modeling time series seasonality (e.g., References [67,68]). In this study, we used a harmonic regression model to characterize seasonality and potential forest changes in each pixel location. We fitted an ordinary least square (OLS) model, which is almost the same as the Continuous Change Detection and Classification (CCDC) algorithm [18], to the time series of all Landsat 8 spectral indices and Sentinel-1 backscatter coefficients: where y is a predicted value at time x, x is the Julian date, a 0 represents a coefficient for overall value, a 1 and b 1 represent coefficients for intra-annual change, c 1 represents a coefficient for inter-annual change and is the residual error. The CCDC algorithm was originally developed for Landsat spectral bands; however, several studies have successfully applied it to different sensors and different spectral indices (e.g., References [28,68,69]) for change detection. The OLS model was initially fitted using observations from 1 January 2014 to 31 December 2015. After that period, the harmonic regression model was used to predict time series observations. Based on the harmonic regression model prediction and actual observation values, the root mean square error (RMSE) was calculated and used to identify break points as potential change. If three consecutive observations deviated by more than three times the RMSE, a break point was labeled at that time ( Figure 2). After identifying the break point, a new harmonic regression model was fitted from the next observation until a new break point was detected. Conversely, if three consecutive observations did not deviate by more than three times the RMSE, we assumed there was no break point at that time. The harmonic regression model was updated once new observations became available, and the RMSE was recalculated every time. Thus, different coefficients and RMSEs were obtained when new observations became available. We recorded these coefficients and the RMSE at each observation. A detailed description and illustration of the harmonic regression model can be found in Reference [18].
The coefficients and RMSE of each observation from the harmonic regression models were used as predictor variables in the RF models. Different from the original CCDC algorithm, the coefficients, RMSEs and observed values were extracted for every observation, even if there was no break point. This was because the RF model was applied to each time series observation in the subsequent procedure. Together with variables from the harmonic regression models, the topographic slope and elevation, and incident angle were also included as predictor variables. The RF models were independently developed for Landsat 8 and Sentinel-1. This means that variables from the harmonic regression models for TCB, TCG, TCW, TCA and NBR were used in the RF model of Landsat 8, and variables from VV and VH were used in the RF model of Sentinel-1. We built RF disturbance probability models for each observation of Landsat 8 and Sentinel-1. Because RF classifiers can be used to predict the probabilities of binary responses [70], several studies have used them for probability mapping with remote sensing data (e.g., References [50,71]). In this study, we predicted the probability of forest disturbances, defined as any discrete event that caused a reduction of forest canopy. Following land cover type was not considered, therefore disturbances may include conversion to other land cover types (e.g., infrastructure and agricultural land). We set disturbed and not-disturbed (stable forest) conditions as the binary outcome of the RF models. The disturbance probability was calculated as the proportion of disturbed class predicted by decision trees in the RF classifier. The RF models were independently trained using observations derived from 535 training samples collected within forested areas of the study area (see Section 3.4). In the training samples, each observation of Landsat 8 and Sentinel-1 time series data was labeled as disturbed or not-disturbed. We randomly selected five observations from each of the 535 training samples to use as training data in each RF model. Thus, a total of 2675 observations from both Landsat 8 and Sentinel-1 were used to train the RF models. Here, each observation came from any acquisition date of time series Landsat 8 and Sentinel-1 data during the monitoring period (i.e., from 1 January 2016 to 31 December 2017) and had a reference label (disturbed or not-disturbed) and predictor variables. The RF models were developed using the "randomForest" package [72] and tuned with 10-fold cross-validation of 500 trees using the "caret" package [73] in the R statistical computing environment [74].

Detection of Forest Disturbance Using Time Series Disturbance Probabilities
The developed RF models were used to predict the disturbance probabilities of time series observations, and these time series disturbance probabilities were then used to detect forest disturbance in each pixel location. We tested the disturbance detection performance using only the disturbance probabilities of Landsat 8, using only the disturbance probabilities of Sentinel-1 and using the combined disturbance probabilities of Landsat 8 and Sentinel-1 from the monitoring period 1 January 2016 to 31 December 2017. We used a simple approach to detect disturbances: if three consecutive observations exceeded a disturbance probability of 50%, then that pixel was assumed to indicate a disturbance. We required three consecutive observations to remove false positive detections. When combining the probabilities of Landsat 8 and Sentinel-1, we simply used either of the probabilities from Landsat 8 or Sentinel-1 based on the confidence of the probability. This was measured by the closeness to 0 or 1, expressed as follows: where C i represents the confidence of the disturbance probability for the i-th observation of Landsat 8 or Sentinel-1 and P(d i ) represents the disturbance probability for the i-th observation of Landsat 8 or Sentinel-1. We assumed that observations with more extreme probabilities represented more confident predictions. From the beginning of the time series, the disturbance probability was updated when new observations became available. If a new observation was obtained from a different sensor, the confidence of the probability of this new observation was compared against that of the previous disturbance probability. If the confidence was higher for the new observation, the disturbance probability was updated; if the confidence was lower for the new observation, the previous disturbance probability was retained. By doing this, we combined the disturbance probabilities from different sensors.
To understand the variability of forest disturbance detection using different thresholds for time series disturbance probabilities, we changed the threshold used to identify the occurrence of disturbance. We tested thresholds from 10-90% at intervals of 2.5%. If the disturbance probabilities of three consecutive observations exceeded the threshold, the pixel was labeled as a disturbance.
We compared disturbance detection using only Landsat 8, using only Sentinel-1 and using the combined Landsat 8 and Sentinel-1 data.

Collection of Reference Samples
We collected reference samples for the two-year monitoring period (i.e., from 1 January 2016 to 31 December 2017). Stratified random sampling was employed to collect samples. We used forest disturbance maps from 2016 and 2017 in this study area that were derived from time series analysis of annual Landsat data [75]. These disturbance maps were originally developed for the period 2000-2014, but additional Landsat data from 2000-2019 were used to generate the disturbance maps used in this study. A total of 1070 reference samples were collected in forested regions of the study area. We allocated 405 samples to areas that experienced disturbances in the monitoring period and the remaining 665 samples to areas that appeared to be stable. Reference samples were randomly distributed over each stratum in the forested regions of the study area. We used a single pixel as the sample unit.
We recorded the presence/absence of disturbance and the date of disturbance occurrence for the reference samples based on a visual interpretation of high spatial resolution satellite images. PlanetScope and RapidEye images [57] were used to identify disturbances in this visual interpretation. The PlanetScope satellites form a constellation of more than 150 CubeSats and provide images with a spatial resolution of about 3 m. The PlanetScope and RapidEye images available for reference samples varied over the monitoring period in this study ( Figure 4) largely because of the variation of satellites on orbits. The weighted average of available images for each sample was 3.3 images/month (i.e., 9.1 day interval), including cloud-contaminated images. Although this figure included cloudy observations, which might reduce the number of images available for visual interpretation, we considered this an acceptable frequency for visual interpretation. In addition, time series of Landsat 8 data and high spatial resolution images from Google Earth were used as complementary data when appropriate PlanetScope and RapidEye images were not available. We also recorded disturbance agents (i.e., agricultural development, timber plantation development, logging, shifting cultivation, infrastructure development and water inundation) in the visual interpretation to investigate potential differences in detection for different causes of disturbance ( Figure 5). Although several small-scale disturbances, such as caused by selective logging, might be omitted in visual interpretation, we considered that these omissions accounted for a small amount because most disturbances caused by selective logging operations can be identified using high spatial resolution satellite images [76].

Validation
The performance of the RF models was assessed based on the confusion matrix formed using the 2675 observations from the test samples. We set a disturbance probability of 50% as the threshold for identifying disturbances. The overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA) were calculated using the confusion matrix. In addition, the Brier score [77] was calculated to investigate the prediction performance. This is a measure of probabilistic prediction accuracy expressed as: where is the Brier score of Landsat 8 or Sentinel-1, represents the total number of observations in the RF model, represents the reference label indicating not-disturbed (0) or disturbed (1) for the i-th observation of Landsat 8 or Sentinel-1 and represents the predicted disturbance Figure 5. Distribution of samples that experienced disturbance in reference samples. We classified these disturbances for six disturbance agents.
Each Landsat 8 and Sentinel-1 observation in the reference samples was labeled as disturbed or not-disturbed based on the results of visual interpretation. We labeled observations as not-disturbed when there was no sign of disturbance in the sample. If there was disturbance in the sample, observations before the date of disturbance were labeled as not-disturbed and observations after disturbance were labeled as disturbed. The reference samples had, on average, 35.4 and 50.8 observations for Landsat 8 and Sentinel-1, respectively.
We allocated 50% of the reference samples for training and used the remainder as test samples. RF disturbance probability models were developed using the observations of training samples and then validated using the observations of test samples. In developing and validating the RF models, five observations were randomly selected from each training and test samples. Thus, we used 2675 observations for training and testing the RF models of both Landsat 8 and Sentinel-1.

Validation
The performance of the RF models was assessed based on the confusion matrix formed using the 2675 observations from the test samples. We set a disturbance probability of 50% as the threshold for identifying disturbances. The overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA) were calculated using the confusion matrix. In addition, the Brier score [77] was calculated to investigate the prediction performance. This is a measure of probabilistic prediction accuracy expressed as: where BS is the Brier score of Landsat 8 or Sentinel-1, N represents the total number of observations in the RF model, R i represents the reference label indicating not-disturbed (0) or disturbed (1) for the i-th observation of Landsat 8 or Sentinel-1 and P(d i ) represents the predicted disturbance probability for the i-th observation of Landsat 8 or Sentinel-1. The Brier score measures the average squared deviation of predicted probability and actual outcome [77]; lower values indicate higher prediction accuracy. The disturbance detection was validated based on 535 test samples. We calculated OA, PA and UA in each disturbance detection (i.e., only with Sentinel-1, only with Landsat 8 and combining Landsat 8 and Sentinel-1). We treated disturbances detected before the actual disturbance date as commission errors of disturbance, in accordance with other studies (e.g., References [20,21,35]). However, disturbances detected within 10 days before the actual date were regarded as correct detections, because the visual interpretation had an uncertainty of 9.1 days (the weighted average interval of available PlanetScope and RapidEye images used in the visual interpretation). We calculated the mean time lag of disturbance detection, which was defined by mean days between the reference date of disturbance and the initial date of three consecutive disturbance probabilities of greater than 50%, to investigate the temporal accuracy of detection. In this calculation, disturbances detected before the date of the actual disturbance were not included. The OA, PA, UA and mean time lag were calculated for disturbance detection with thresholds from 10-90%. Furthermore, we calculated the omission error of disturbance for each disturbance agent, which was recorded in the reference samples, using the disturbance detection results combining Landsat 8 and Sentinel-1 with a threshold of 50%. We also investigated the effects of spatial post processing to mapped disturbance pixels with a threshold of 50% by applying minimum mapping unit from 1 to 100 pixels.

Accuracy of RF Disturbance Probability Models
The RF models for each observation of Landsat 8 and Sentinel-1 had an OA of 88.9% and 83.3%, respectively, when observations that had disturbance probabilities of greater than 50% were labeled as disturbances. In the RF model for Landsat 8, the PA of disturbances was 78.6% (Table 1). In the RF model for Sentinel-1, however, the PA of disturbances was 58.6%, somewhat lower than that of the RF model for Landsat 8 ( Table 2). The PA of no disturbance was slightly higher for Sentinel-1 (94.5%) than for Landsat 8 (93.5%). The UAs of disturbance and no disturbance for Landsat 8 were higher than those for Sentinel-1. The Brier scores were 0.086 and 0.128 for the Landsat 8 and Sentinel-1 models, respectively, indicating that the RF model for Landsat 8 provided better predictions.

Accuracy of Disturbance Detection Using Time Series Disturbance Probabilities
The OAs of forest disturbance detection using only Landsat 8, using only Sentinel-1 and using a combination of Landsat 8 and Sentinel-1 were 78.3%, 75.5% and 83.6%, respectively (Figure 6a), when a threshold of 50% was used to identify disturbances. The disturbance detection with Landsat 8 and Sentinel-1 had the highest OA, followed by that using only Landsat 8 and that using only Sentinel-1. In terms of the UA of disturbance, the highest accuracy was again achieved by combining Landsat 8 and Sentinel-1 (71.5%). Conversely, for the PA of disturbance, using only Landsat 8 achieved slightly higher accuracy (84.7%) than the combination of Landsat 8 and Sentinel-1 (82.2%). The mean time lags were 63.1 days and 86.9 days for disturbance detection using only Landsat 8 and using only Sentinel-1, respectively (Figure 6b). When combining Landsat 8 and Sentinel-1, the mean time lag was 52.8 days, some 10.3 days less than detection using only Landsat 8. A total of 200 samples were detected as disturbance and 335 samples were predicted as stable in test samples during the monitoring period. Figure 7 shows an example of disturbance detection using Landsat 8, Sentinel-1 and the combination of Landsat 8 and Sentinel-1. Visually, the detected disturbances were largely reduced by combining Landsat 8 and Sentinel-1 observations compared with using only Landsat 8 or Sentinel-1 observations. The OAs of forest disturbance detection using only Landsat 8, using only Sentinel-1 and using a combination of Landsat 8 and Sentinel-1 were 78.3%, 75.5% and 83.6%, respectively (Figure 6a), when a threshold of 50% was used to identify disturbances. The disturbance detection with Landsat 8 and Sentinel-1 had the highest OA, followed by that using only Landsat 8 and that using only Sentinel-1. In terms of the UA of disturbance, the highest accuracy was again achieved by combining Landsat 8 and Sentinel-1 (71.5%). Conversely, for the PA of disturbance, using only Landsat 8 achieved slightly higher accuracy (84.7%) than the combination of Landsat 8 and Sentinel-1 (82.2%). The mean time lags were 63.1 days and 86.9 days for disturbance detection using only Landsat 8 and using only Sentinel-1, respectively (Figure 6b). When combining Landsat 8 and Sentinel-1, the mean time lag was 52.8 days, some 10.3 days less than detection using only Landsat 8. A total of 200 samples were detected as disturbance and 335 samples were predicted as stable in test samples during the monitoring period. Figure 7 shows an example of disturbance detection using Landsat 8, Sentinel-1 and the combination of Landsat 8 and Sentinel-1. Visually, the detected disturbances were largely reduced by combining Landsat 8 and Sentinel-1 observations compared with using only Landsat 8 or Sentinel-1 observations.
(a) (b) The omission error of disturbance for each disturbance agent is shown in Figure 8. Water inundation achieved the lowest omission error (0.0%), although the sample size for this agent was small (n = 28) during the monitoring period. Infrastructure development also achieved a low omission error (4.8%). The highest omission error occurred for the detection of disturbances caused by logging (41.7%). The development of timber plantations, shifting cultivation and agricultural development had similar omission errors (15.0%, 15.9% and 19.4%, respectively).
The OAs of forest disturbance detection using only Landsat 8, using only Sentinel-1 and using a combination of Landsat 8 and Sentinel-1 for each minimum mapping unit are shown in Figure 9. The highest OA occurred at minimum mapping unit of 32, 3 and 7 pixels for Landsat 8 (80.0%), Sentinel-1 (76.3%) and the combination of Landsat 8 and Sentinel-1 (84.9%), respectively. The omission error of disturbance for each disturbance agent is shown in Figure 8. Water inundation achieved the lowest omission error (0.0%), although the sample size for this agent was small (n = 28) during the monitoring period. Infrastructure development also achieved a low omission error (4.8%). The highest omission error occurred for the detection of disturbances caused by logging (41.7%). The development of timber plantations, shifting cultivation and agricultural development had similar omission errors (15.0%, 15.9% and 19.4%, respectively).
The OAs of forest disturbance detection using only Landsat 8, using only Sentinel-1 and using a combination of Landsat 8 and Sentinel-1 for each minimum mapping unit are shown in Figure 9.       Figure 10 shows the accuracy of disturbance detection using only Landsat 8, using only Sentinel-1 and combining Landsat 8 and Sentinel-1 for different thresholds. In terms of the OA, PA and UA of disturbance, Sentinel-1 generally achieved a lower accuracy for all thresholds, except for the UA of disturbance for thresholds greater than 40% (Figure 10a). Combining Landsat 8 and Sentinel-1 produced higher OAs than disturbance detection using only Landsat 8 for all threshold values. We found that the PA of disturbance was higher for disturbance detection combining Landsat 8 and Sentinel-1 with thresholds greater than 52.5%, whereas disturbance detection using only Landsat 8 generally achieved higher accuracy with thresholds of 50% and below. The UA of disturbance detection combining Landsat 8 and Sentinel-1 was higher than the others for thresholds below 60%, but disturbance detection using only Landsat 8 or only Sentinel-1 was generally higher when the threshold was over 62.5%. The highest OA occurred at thresholds of 72.5%, 62.5% and 77.5% for Landsat 8 (82.4%), Sentinel-1 (78.3%) and the combination of Landsat 8 and Sentinel-1 (86.0%), respectively. The PA and UA of disturbance achieved almost the same balanced accuracy at thresholds of 67.5%, 55.0% and 70.0% for Landsat 8 (PA: 70.7%, UA: 72.7%), Sentinel-1 (PA: 62.9%, UA: 62.2%) and the combination of Landsat 8 and Sentinel-1 (PA: 79.5%, UA: 79.5%), respectively. In terms of the mean time lag, Sentinel-1 took longer to detect disturbances for all threshold values (Figure 10b). Combining Landsat 8 and Sentinel-1, the mean time lag of disturbance detection was generally shortest for thresholds above 47.5%. For disturbance detection below that threshold, the mean time lag was shorter when using Landsat 8 only.  Figure 10 shows the accuracy of disturbance detection using only Landsat 8, using only Sentinel-1 and combining Landsat 8 and Sentinel-1 for different thresholds. In terms of the OA, PA and UA of disturbance, Sentinel-1 generally achieved a lower accuracy for all thresholds, except for the UA of disturbance for thresholds greater than 40% (Figure 10a). Combining Landsat 8 and Sentinel-1 produced higher Oas than disturbance detection using only Landsat 8 for all threshold values. We found that the PA of disturbance was higher for disturbance detection combining Landsat 8 and Sentinel-1 with thresholds greater than 52.5%, whereas disturbance detection using only Landsat 8 generally achieved higher accuracy with thresholds of 50% and below. The UA of disturbance detection combining Landsat 8 and Sentinel-1 was higher than the others for thresholds below 60%, but disturbance detection using only Landsat 8 or only Sentinel-1 was generally higher when the threshold was over 62.5%. The highest OA occurred at thresholds of 72.5%, 62.5% and 77.5% for Landsat 8 (82.4%), Sentinel-1 (78.3%) and the combination of Landsat 8 and Sentinel-1 (86.0%), respectively. The PA and UA of disturbance achieved almost the same balanced accuracy at thresholds of 67.5%, 55.

Discussion
In this study, we investigated the detection of disturbances in the tropical seasonal forests of Myanmar using dense time series of Landsat 8 and Sentinel-1 data based on RF models. Previous studies mainly used parametric models to combine time series data from optical and SAR satellites for change detection, and other approaches, such as using machine learning algorithms, have not been well developed. This study presented a new approach to detect disturbances by combining disturbance probabilities of Landsat 8 and Sentinel-1 predicted by RF models. The RF disturbance probability models for each observation of Landsat 8 and Sentinel-1 were separately developed using predictor variables derived from a harmonic regression model. Combining the time series disturbance probabilities of Landsat 8 and Sentinel-1 observations improved the accuracy of disturbance detection compared with using only Landsat 8 or Sentinel-1.

Discussion
In this study, we investigated the detection of disturbances in the tropical seasonal forests of Myanmar using dense time series of Landsat 8 and Sentinel-1 data based on RF models. Previous studies mainly used parametric models to combine time series data from optical and SAR satellites for change detection, and other approaches, such as using machine learning algorithms, have not been well developed. This study presented a new approach to detect disturbances by combining disturbance probabilities of Landsat 8 and Sentinel-1 predicted by RF models. The RF disturbance probability models for each observation of Landsat 8 and Sentinel-1 were separately developed using predictor variables derived from a harmonic regression model. Combining the time series disturbance probabilities of Landsat 8 and Sentinel-1 observations improved the accuracy of disturbance detection compared with using only Landsat 8 or Sentinel-1.
The RF disturbance probability model of Landsat 8 produced better results than that of Sentinel-1. The RF model of Landsat 8 achieved a relatively high OA and low Brier score. In this study, we used a harmonic regression model based on the CCDC algorithm to account for both seasonality and disturbance-related changes. This approach provided model coefficients and RMSEs representing inter-annual and intra-annual changes and random variability. These variables effectively allowed the RF model to detect disturbances. Furthermore, the harmonic regression model in this study could reduce the effects of outliers caused by undetected clouds in cloud masking because three consecutive observations were used to identify potential change. We consider this was one reason for relatively high accuracy. The RF disturbance probability model of Sentinel-1 had a low PA of disturbance (58.6%), indicating large omission errors of disturbance. This result was not surprising because C-band SAR is generally less sensitive to forest structures [42,43]. Although we applied careful preprocessing to remove noise from the SAR backscatter coefficients, such as terrain correction, the mountainous regions in the study area might have affected the accuracy of detection, because rugged terrain affects the geometric and radiometric conditions of SAR observations. In GEE, Sentinel-1 GRD data is processed to backscatter coefficients σ • ; however, using backscatter coefficients γ • might be another option because of its performance for topographic variation [78]. Although we cannot directly use Sentinel-1 data in GEE when backscatter coefficients γ • are generated, it could improve the accuracy of detection. Another solution for improving detection accuracy is to use other preprocessing methods for Sentinel-1 data, such as a multi-temporal SAR filter to reduce speckle noise. This has a possibility to contribute better results for change detection.
We found that using the disturbance probabilities of both Landsat 8 and Sentinel-1 achieved highly accurate disturbance detection. The combination of Landsat 8 and Sentinel-1 improved the OA and UA of disturbances using a threshold value of 50% (Figure 6a). Although the PA of disturbances was slightly lower than disturbance detection using only Landsat 8, the quantitative difference was small. In terms of temporal accuracy, the combination of Landsat 8 and Sentinel-1 reduced the time required to detect disturbances compared with detection using only Landsat 8. Considering the temporal distribution of the reference samples, most of the forest disturbances in this region occurred in the dry season ( Figure 5). Forest disturbances in the dry season can be observed by optical sensors (i.e., Landsat 8) immediately after the event; however, the combination of Landsat 8 and Sentinel-1 provided better results for timely and accurate detection through the increased number of observations. These results revealed the value of the frequent observations achieved by using both Landsat 8 and Sentinel-1. A similar contribution from high temporal observations using multiple sensors was demonstrated in previous studies [35]. Although Sentinel-1 alone did not perform well in terms of disturbance detection, its combination with Landsat 8 improved the accuracy of disturbance detection. In this context, combining Landsat 8 and Sentinel-1 observations can contribute to successful disturbance detection while reducing commission errors of disturbance.
Increasing the thresholds required to identify disturbances led to a decrease in the PA of disturbances and an increase in the UA of disturbances. The mean time lag of detection also generally increased with the threshold value. This tradeoff confirms similar time series change detection results from previous studies [22,35,79]. No single threshold achieved the highest PA and UA of disturbance and the shortest mean time lag of detection. Although we used a threshold of 50% to identify disturbances for analysis and mapping, it is possible to use other thresholds according to the objective of disturbance detection. To detect most forest disturbances, lower thresholds could achieve low omissions of disturbances, although false positive detections will increase. In contrast, higher thresholds (e.g., 70%) could be used to achieve balanced PA and UA of detection and mean time lag of detection. Applying a minimum mapping unit to mapped disturbances also might be used for balancing PA and UA of detection because this procedure was effective to eliminate false positive detection in this study.
Although the combination of Landsat 8 and Sentinel-1 achieved relatively high accuracy, the omission of disturbances was considerable for several disturbance agents. The omission errors for each disturbance agent showed that conversion to other land cover type, such as water inundation and infrastructure development, can easily be detected, whereas subtle disturbances such as logging were difficult to detect. In this study, we detected forest disturbances that included subtle or small events as well as deforestation. Logging activity in this region is not usually conducted by clearing large-scale forest stands, but by small-scale tree-based harvesting [80] (i.e., selective logging). This might increase the omission error of disturbances because it is difficult to detect small-scale disturbances such as those caused by selective logging using medium spatial resolution images (e.g., Landsat) [81]. One possible solution would be to use Sentinel-2 data combined with Sentinel-1 data. Small-scale disturbances might be detected using Sentinel-2 data because we could conduct disturbance detection in 10 m spatial resolution and with higher temporal resolution. Although Sentinel-2 data was not used in this study due to the availability of bottom of atmosphere reflectance products during the monitoring period in GEE, further studies are needed to investigate the effectiveness of Sentinel-2 data for detecting small-scale disturbances. We found that the forest disturbances that caused land cover changes were easily distinguished from stable forests by using Landsat 8 and Sentinel-1. This can be attributed to the different surface properties of water and city compared with forest areas.
The approach presented in this paper has demonstrated the utility of disturbance probabilities predicted by a machine learning algorithm (i.e., RF) for accurate and timely disturbance detection.
Although several studies have improved the accuracy of forest monitoring using a Bayesian approach for optical and SAR data, this approach requires several premises, such as probability density functions for Bayesian updating, to combine the conditional probabilities from both sensors (e.g., [22,35]). Conversely, we can easily implement the proposed approach, because it simply uses either of the disturbance probabilities of Landsat 8 or Sentinel-1 without prior knowledge. Compared with the Bayesian approach in previous studies, such as Reiche et al. [35], the accuracy of disturbance detection was low in this study. Reiche et al. [35] detected deforestation with the PA and UA of 88.9% and 88.0%, respectively, for a high confidence in change area case (82.2% and 71.5% in this study). The mean time lag of detection in this study was also longer compared with those of Reiche et al. [35] (i.e., 11 days in the same definition of mean time lag) and another study by Reiche et al. [44] using Sentinel-1 data (i.e., 7.5 days). The difference of forest change patterns might partially explain this lower accuracy in this study because our study included small-scale disturbances that are usually difficult to detect. Although some small-scale disturbances were omitted in our approach, most forest disturbances were successfully captured. The simplicity and accuracy of the proposed approach are suited to operational forest monitoring. Thus, the approach presented here, which enhances the disturbance probability predicted by a machine learning algorithm with variables from a harmonic regression model, represents a promising alternative for forest monitoring.
Although the presented approach produced encouraging results, it has several limitations. The RF classifier was applied to each observation of Landsat 8 and Sentinel-1, therefore it was highly dependent on computational resources. Even though RF is a highly efficient processing algorithm [52], it takes a considerable time for large-area mapping. The availability of reference images was another constraint for accurate detection. To examine the exact date of disturbance occurrence, adequate high spatial resolution reference images are necessary, especially for detecting subtle or small-scale changes. In this study, we used PlanetScope and RapidEye images that had an average interval of 9.1 days as reference images; however, more frequent reference images might be needed for near-real-time scenarios, because temporal accuracy is important for near-real-time monitoring. Recently, however, PlanetScope satellites have provided more frequent, near-daily observations as the number of satellites has increased. This will aid the reliability of data sources for reference images in time series analysis.

Conclusions
This study assessed the effectiveness of using dense time series Landsat 8 and Sentinel-1 data for the detection of forest disturbances in tropical seasonal forests. The disturbance probability of each observation was predicted from a machine learning algorithm (Random Forest) using predictor variables from harmonic regression models. Combining the time series of disturbance probabilities from Landsat 8 and Sentinel-1 data based on a simple approach, we detected disturbances with an OA of 83.6%, which was higher than disturbance detection using only Landsat 8 (78.3%) or only Sentinel-1 (75.5%). In addition, the combination of Landsat 8 and Sentinel-1 data effectively reduced the mean time lag for detecting disturbances. Although disturbance detection using only Sentinel-1 data achieved a relatively low accuracy, Sentinel-1 data were effective for accurate and timely detection when combined with Landsat 8 data. We showed the tradeoff between reducing the omission error and commission error of disturbances by changing the thresholds required to identify disturbances from the time series disturbance probabilities. This suggests that the optimal thresholds vary depending on the objective of analysis. Subtle or small-scale disturbances caused by logging activity were difficult to detect; however, the presented approach can detect most large-scale disturbances, such as land cover change from forest. Overall, using a machine learning algorithm was effective for detecting disturbances in an accurate and timely manner. A simple approach to combining the disturbance probabilities of multiple sensors is a promising approach for tropical seasonal forests.