Forest Disturbance Mapping Using Dense Synthetic Landsat / MODIS Time-Series and Permutation-Based Disturbance Index Detection

Spatio-temporal information on process-based forest loss is essential for a wide range of applications. Despite remote sensing being the only feasible means of monitoring forest change at regional or greater scales, there is no retrospectively available remote sensor that meets the demand of monitoring forests with the required spatial detail and guaranteed high temporal frequency. As an alternative, we employed the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) to produce a dense synthetic time series by fusing Landsat and Moderate Resolution Imaging Spectroradiometer (MODIS) nadir Bidirectional Reflectance Distribution Function (BRDF) adjusted reflectance. Forest loss was detected by applying a multi-temporal disturbance detection approach implementing a Disturbance Index-based detection strategy. The detection thresholds were permutated with random numbers for the normal distribution in order to generate a multi-dimensional threshold confidence area. As a result, a more robust parameterization and a spatially more coherent detection could be achieved. (i) The original Landsat time series; (ii) synthetic time series; and a (iii) combined hybrid approach were used to identify the timing and extent of disturbances. The identified clearings in the Landsat detection were verified using an annual woodland clearing dataset from Queensland’s Statewide Landcover and Trees Study. Disturbances caused by stand-replacing events were successfully identified. The increased temporal resolution of the synthetic time series indicated promising additional information on disturbance timing. The results of the hybrid detection unified the benefits of both approaches, i.e., the spatial quality and general accuracy of the Landsat detection and the increased temporal information of synthetic time series. Results indicated that a temporal improvement in the detection of the disturbance date could be achieved relative to the irregularly spaced Landsat data for sufficiently large patches.


Introduction
Forests worldwide are regarded as prime areas for highly valuable Ecosystem Services (ESS) and are as carbon sinks invaluable for human life on earth [1].Remote sensing offers ways to explicitly report on changes in forest extent and cover and has become a pivotal tool for reporting on carbon budgeting [2] supporting the development of climate policies, and project future climate change.The Normalized Difference Vegetation Index (NDVI) has been utilized in numerous satellite image based studies on forest and vegetation change in general [3].Landsat-5 Thematic Mapper (TM) and Landsat-7 Enhanced Thematic Mapper plus (ETM+) data have widely been used for land surface mapping and monitoring (e.g., [4][5][6]).However, process monitoring was always somewhat limited due to the low temporal frequency (16 days orbital repeat) and image cost.However, the release of the Landsat archive at no cost [7] has provided new opportunities for assessing historical changes in landscapes and enabling the discrimination between surface features that single date images cannot differentiate [8].Although marking a milestone in environmental monitoring is the temporal coverage with Landsat data alone still somewhat limited as cloud cover, adverse atmospheric conditions and sensor problems can cause significant data gaps [6,9].
A number of data fusion techniques have been developed to merge data from the high temporal repeat MODIS instrument and Landsat sensor data (e.g., [10][11][12]).Of note is the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM, [10]), which has been applied successfully for a range of studies (e.g., [10,13,14]).For example, Walker et al. [15] demonstrated the usefulness of the STARFM algorithm for phenological studies in the drylands of Arizona by using 12 STARFM predicted images.Schmidt et al. [16] used STARFM to predict 333 Landsat sensor images over a 7.5 years time period, with these then used to study regional ecological and phenological processes in a heterogeneous Northern Australian savanna.Bhandari et al. [17] described forest phenology in an Australian savanna by filling in fused images into a Landsat time series.Schmidt et al. [18] used 12 years STARFM time series with eight-day interval in a temporal change detection approach using the "Breaks For Additive Season and Trend" approach (BFAST, [19]).Tewes et al. [20] merged RapidEye and MODIS imagery for one year in an African savanna with 250 m MODIS data using the enhanced STARFM algorithm (ESTARFM, [12]).ESTARFM and the modified mESTARFM [21] have been successfully applied as an improvement on STARFM, although findings of [22] are inconclusive in establishing whether this improvement is evident for all environments (e.g., in areas with a dominant temporal variance, STARFM gave a superior performance).Therefore, the less complex STARFM algorithm was applied here.Hilker et al. [23] implemented a Spatial Temporal Adaptive Algorithm for mapping Reflectance Change (STAARCH) in order to detect disturbances and better parameterize STARFM.Similar to many other disturbance detection approaches [24][25][26][27], STAARCH relies on the Disturbance Index (DI, [24]), which is a linear transformation of the Tasseled Cap [28] indices.The DI was used in numerous forest clearing studies [23][24][25][26][27] and was developed to separate completely deforested areas from dark closed-canopy coniferous forests [24].
In this study, we implemented an approach that is capable of detecting disturbances in both space and time.Our method is conceptually related to the stand replacing disturbance detection algorithm, which is part of STAARCH.One criticism of STAARCH has been the fixed disturbance thresholds, which need to be given by the user.As such, we expanded on the method by implementing a more robust parameterization, where detection thresholds are permutated with random numbers for the normal distribution.Disturbances are thus detected with several threshold sets within a broader threshold confidence space.Our approach is able to process observed and irregularly spaced Landsat images, as well as synthetically derived, equidistant Landsat-like fused images.In addition, a hybrid approach was implemented that utilizes both data sources.We tested if process-based change such as clear-cut logging of forest or regrowth vegetation, stand-removing fires, multiple disturbances or subtle change can be detected by our approach.We also tested how the synthetic images perform in comparison with observed imagery.

Study Area
The study area (22,500 km²) is located within the World Reference System 2 (WRS-2) path/row 093/078 (see Figure 1).The structural formations of open forests and woodlands dominating the area are characteristic savanna environments, accounting for over 70% of Australian forests in terms of structure and biomass [29].The major tree communities in this area are formed by evergreen eucalyptus, acacia or callitris dominated woodlands to open forests [30] (see Figure 2).The analysis was restricted to the forested areas shown in Figure 2, i.e., where Foliage Projective Cover (FPC) [31] of our first year of interest exceeded 15%.The region is characterized by variable rainfall and the ecosystem itself is generally water limited.The average monthly evaporation exceeds the average monthly rainfall throughout the year [32] with a pronounced dry and wet season, where the rainfall is precipitated by short duration storms with high temporal and spatial variability [32].Wildfires are present in the study area, though they happen irregularly.They can alter the species' composition or dominance as some species (e.g., Callitris glaucophylla) are less resistant to fire than others [33,34].Anthropogenic clearing is also commonplace and mainly aims to facilitate establishment of cattle pasture [33].Commonly, bulldozers are used to fell the trees.Stem injection and selective harvesting of Callitris glaucophylla in State Forests are secondary disturbance mechanisms [29].
Remote Sens. 2016, 8, 277 3 of 22 our first year of interest exceeded 15%.The region is characterized by variable rainfall and the ecosystem itself is generally water limited.The average monthly evaporation exceeds the average monthly rainfall throughout the year [32] with a pronounced dry and wet season, where the rainfall is precipitated by short duration storms with high temporal and spatial variability [32].Wildfires are present in the study area, though they happen irregularly.They can alter the species' composition or dominance as some species (e.g., Callitris glaucophylla) are less resistant to fire than others [33,34].Anthropogenic clearing is also commonplace and mainly aims to facilitate establishment of cattle pasture [33].Commonly, bulldozers are used to fell the trees.Stem injection and selective harvesting of Callitris glaucophylla in State Forests are secondary disturbance mechanisms [29].Broad vegetation groups of the forested study area (Regional Ecosystem Mapping of 2006 [30]).The analysis was restricted to the forested areas (colors); the "other" class was compiled using  our first year of interest exceeded 15%.The region is characterized by variable rainfall and the ecosystem itself is generally water limited.The average monthly evaporation exceeds the average monthly rainfall throughout the year [32] with a pronounced dry and wet season, where the rainfall is precipitated by short duration storms with high temporal and spatial variability [32].Wildfires are present in the study area, though they happen irregularly.They can alter the species' composition or dominance as some species (e.g., Callitris glaucophylla) are less resistant to fire than others [33,34].Anthropogenic clearing is also commonplace and mainly aims to facilitate establishment of cattle pasture [33].Commonly, bulldozers are used to fell the trees.Stem injection and selective harvesting of Callitris glaucophylla in State Forests are secondary disturbance mechanisms [29].Broad vegetation groups of the forested study area (Regional Ecosystem Mapping of 2006 [30]).The analysis was restricted to the forested areas (colors); the "other" class was compiled using Figure 2. Broad vegetation groups of the forested study area (Regional Ecosystem Mapping of 2006 [30]).The analysis was restricted to the forested areas (colors); the "other" class was compiled using the 15% Foliage Projective Cover threshold, and is mainly composed of grasslands and previously cleared forests.The outlines of the National Parks are superimposed in green colors.

Landsat Data
All available Landsat-5 TM data of WRS-2 path/row 093/078 for the three years period of 2007-2009 were considered in this study.Landsat-7 ETM+ data were ignored because of substantial data loss [35] due to the Scan Line Corrector (SLC) failure.All Landsat data were standard terrain corrected (Level 1T), processed by the Level 1 Product Generation System (LPGS) maintained by the USGS.Images that were not corrected to L1T precision were discarded in order to guarantee a coherent co-registration among the images.The Landsat data were corrected for atmospheric, topographic and bi-directional effects [36] by deriving nadir BRDF adjusted reflectance (NBAR; based on a sun zenith angle of 45 ˝).Automatic cloud screening was performed by a two-step method that utilizes the Fmask algorithm [37] and a subsequent outlier detection based on the results of the first step [38].

MODIS Data
In order to match the Landsat data and to account for MODIS NBAR reflectance being the best product for STARFM predictions with Landsat-5 [15], we analogously produced a MODIS NBAR dataset based on a fixed sun zenith angle of 45 ˝.We acquired all data of the MODIS Bidirectional Reflectance Distribution Function (BRDF) model parameters (MCD43A1) and the BRDF/Albedo quality product (MCD43A2).The products are provided at 8-day temporal resolution and utilize multi-angular data from 16-day periods to enable the inversion of the BRDF model [39].The quality product was parsed, such that only good or very good NBAR retrievals were used in predicting the synthetic time series [18].The quality masks and NBAR images were reprojected and super-sampled with nearest neighbor resampling to match the spatial characteristics of the Landsat data.

Reference and Auxiliary Datasets
We used the initial estimate of overstory Foliage Projective Cover (i.e., the derived product for 2007) as reference population for the Disturbance Index (see Section 4.2).Overstory FPC is the vertically projected percentage of photosynthetic foliage from tree and shrub life forms and is a biophysically meaningful descriptor of Australian vegetation as trees and shrubs with sparse foliage and irregular crown shapes are predominant [5].Annual estimates are operationally generated by the Department of Science, Information Technology and Innovation, Queensland Governance, on the basis of Landsat-5 TM and Landsat-7 ETM+ imagery.
The performance of the implemented disturbance detection algorithm was assessed by comparison with a well-established clearing dataset.The Queensland's Statewide Landcover and Trees Study (SLATS) [40] provides information on Queensland's forest and woodland extent, and annually reports clearing activities in support of state legislation; i.e., the Vegetation Management Act 1999 [41].The clearing history dates back to the late 1980s.The approach follows a semi-automatic procedure, where an automatic bi-temporal change detection on the basis of annual FPC estimates is performed in the first step.Afterwards, extensive manual editing and quality checking is performed to isolate the woody vegetation change due to anthropogenic reasons, including field inspection where required.Natural biomass losses, e.g., due to fire, are not reported because the resulting land cover change is regarded as non-permanent.

Methods
A schematic workflow summarizing the major parts of this study is given in Figure 3. MODIS and Landsat NBAR data are used to predict and validate a dense synthetic time series using STARFM.A disturbance detection approach was implemented and applied to: (i) the Landsat time series; (ii) the synthetic time series; and (iii) to a hybrid setup.Detection results are compared to reference data and to each other.

Prediction of a Dense Time Series
We produced a dense synthetic Landsat-like time series using STARFM [10].STARFM predicts surface reflectance at Landsat spatial resolution and MODIS temporal frequency, thus making it useful for applications that require high resolution in both time and space.STARFM imagery can either be predicted by using one or two bracketing pairs of corresponding MODIS and Landsat images plus one MODIS image for the desired prediction day [10].A pair of MODIS and Landsat observations is hereafter referred to as base pair.We operated STARFM in the two-base-pair-mode because STARFM is only able to track land cover changes if the disturbance is at least recorded in one of the base images [27] and the overall prediction quality is expected to increase by integrating additional temporal information [10].
Each Landsat image with more than two-thirds of high quality (HQ) observations (i.e., clear-sky pixels as derived by [38]) was used for the base images.The remaining images were discarded.For every Landsat date, the closest MODIS image (given by the middle of the 16-day generation period of MCD43A1) was chosen to complete the base pair.Synthetic Landsat-like images were predicted with a temporal resolution of 8 days.
Figure 4 depicts all available Landsat (a) and MODIS (b) images during the selected three-year period.The Landsat and MODIS base pairs are indicated by the point signature, the target MODIS images are drawn with the line signature.A synthetic image was predicted for each target date on basis of the bracketing base pairs (31 image pairs).The MODIS base images were of high quality (worst scene: 91.95% HQ observations; mean HQ percentage: 98.76%).The band uncertainties for STARFM processing were set to 0.002 and 0.005 for the visible and infrared bands, respectively, and the maximum search distance for finding spectrally similar pixels was set to 750 m [10,13,15].
The predictive quality of STARFM was assessed in various other studies (e.g., [10,13,15]), in part with a similar setup in the same environmental setting (e.g., [16][17][18]).The predicted imagery was generally described to be useful, and as such, we only present a limited quality assessment.We used a cross-validation-like strategy, where a synthetic equivalent to each existing Landsat image was predicted by applying STARFM to the bracketing base pairs.Band-wise linear regressions were performed.Poisson Disc sampling [42] (~10,000 samples) was used to retrieve a random sample in order to avoid under-or overrepresentation of certain areas and to minimize spatial autocorrelation.The coefficient of determination R 2 was used to describe the prediction quality.

Prediction of a Dense Time Series
We produced a dense synthetic Landsat-like time series using STARFM [10].STARFM predicts surface reflectance at Landsat spatial resolution and MODIS temporal frequency, thus making it useful for applications that require high resolution in both time and space.STARFM imagery can either be predicted by using one or two bracketing pairs of corresponding MODIS and Landsat images plus one MODIS image for the desired prediction day [10].A pair of MODIS and Landsat observations is hereafter referred to as base pair.We operated STARFM in the two-base-pair-mode because STARFM is only able to track land cover changes if the disturbance is at least recorded in one of the base images [27] and the overall prediction quality is expected to increase by integrating additional temporal information [10].
Each Landsat image with more than two-thirds of high quality (HQ) observations (i.e., clear-sky pixels as derived by [38]) was used for the base images.The remaining images were discarded.For every Landsat date, the closest MODIS image (given by the middle of the 16-day generation period of MCD43A1) was chosen to complete the base pair.Synthetic Landsat-like images were predicted with a temporal resolution of 8 days.
Figure 4 depicts all available Landsat (a) and MODIS (b) images during the selected three-year period.The Landsat and MODIS base pairs are indicated by the point signature, the target MODIS images are drawn with the line signature.A synthetic image was predicted for each target date on basis of the bracketing base pairs (31 image pairs).The MODIS base images were of high quality (worst scene: 91.95% HQ observations; mean HQ percentage: 98.76%).The band uncertainties for STARFM processing were set to 0.002 and 0.005 for the visible and infrared bands, respectively, and the maximum search distance for finding spectrally similar pixels was set to 750 m [10,13,15].

Disturbance Detection
Our multi-temporal disturbance detection implementation is designed to detect disturbance caused by a variety of stand-replacing processes in a binary sense, i.e., stable vs. disturbed.Disturbance masks are generated for each time step of the input time series.The algorithm consists of four steps: (i) data transformation; (ii) temporal filtering; (iii) multi-temporal disturbance detection; and (iv) spatial filtering.In addition, we also (v) developed a hybrid approach to combine the assets of observed and synthesized data.The predictive quality of STARFM was assessed in various other studies (e.g., [10,13,15]), in part with a similar setup in the same environmental setting (e.g., [16][17][18]).The predicted imagery was generally described to be useful, and as such, we only present a limited quality assessment.We used a cross-validation-like strategy, where a synthetic equivalent to each existing Landsat image was predicted by applying STARFM to the bracketing base pairs.Band-wise linear regressions were performed.Poisson Disc sampling [42] (~10,000 samples) was used to retrieve a random sample in order to avoid under-or overrepresentation of certain areas and to minimize spatial autocorrelation.The coefficient of determination R 2 was used to describe the prediction quality.

Disturbance Detection
Our multi-temporal disturbance detection implementation is designed to detect disturbance caused by a variety of stand-replacing processes in a binary sense, i.e., stable vs. disturbed.Disturbance masks are generated for each time step of the input time series.The algorithm consists of four steps: (i) data transformation; (ii) temporal filtering; (iii) multi-temporal disturbance detection; and (iv) spatial filtering.In addition, we also (v) developed a hybrid approach to combine the assets of observed and synthesized data.

(i)
Data transformation.Our approach is based on the analysis of temporal trajectories of the Disturbance Index (DI) transformation, which was developed to separate completely deforested areas from forest regions [24].It is a linear transformation of the rescaled Tasseled Cap indices [28,43,44]: where the rescaled indices X r are computed by normalizing the respective index X (brightness B, greenness G and wetness W) with the mean µ and standard deviation σ of a reference population: The initial state (2007) of overstory FPC was used to generate a wooded vegetation mask (FPC > 15%) from which the rescaling statistics were retrieved.The disturbance detection was also confined to this area.

(ii)
Temporal filtering.The DI time series might be noise affected, e.g., due to undetected cloud remnants.We eliminate positive spikes above a DI of 0.5 and thereafter negative spikes below or equal to zero.A spike is considered an observation that exceeds the threshold, whereas the bracketing data points are inconspicuous.Disturbances normally trigger a persistent DI response for a longer period.Positive spikes result in false positives if not removed (e.g., remnant clouds).
Similarly, observations following a negative spike can also result in false positives.In the case of dense synthetic imagery, we apply an additional Savitzky-Golay filter [45] to smooth the DI time series and to lessen the noise that results from prediction errors.(iii) Multi-temporal disturbance detection.A threshold-dependent rule set is used to identify potential disturbance events.A potential disturbance at time t is flagged if This query is similar to the STAARCH approach, though [23] tested whether the DI gets beyond DI T , while starting below.Instead, we test if the DI(t) is above DI T and has increased by at least ∆DI T .This change was introduced because some disturbances are characterized by a very high increase in DI, but have a larger base level value than the given threshold from the beginning.The remaining constraints originate from STAARCH and are used to enforce some disturbance related physical properties: recently cleared stands are normally characterized by an increase in brightness, a decrease in NDVI and a decrease in wetness [23].
As it is likely that there are a large number of close misclassifications (i.e., pixels with values slightly smaller than one of the chosen thresholds), we expanded this approach by permutating the five thresholds (DI T , ∆DI T , ∆B r,T , ∆W r,T and ∆NDVI r,T ) with random numbers for the normal distribution.For each threshold, 100 random numbers are generated (the mean µ of the normal distribution is defined by the user and σ = 0.5). Figure 5 illustrates the generation of random numbers for the example of DI T .Randomly generated thresholds are shown in Figure 5a for each permutation and Figure 5b shows the normally-distributed histogram.A random number vector is generated for each threshold and they are subsequently combined to form a multi-dimensional threshold confidence space.A three-dimensional example is shown in Figure 5c, where darker colors indicate values that are closer to the 3-dimensional centroid.We generated a larger amount of combinations for illustration purposes only, to better visualize this three-dimensional threshold confidence space.In our practical implementation, we generate such a point cloud in the five-dimensional threshold space with 100 permutations.
DIT, while starting below.Instead, we test if the DI(t) is above DIT and has increased by at least ΔDIT.This change was introduced because some disturbances are characterized by a very high increase in DI, but have a larger base level value than the given threshold from the beginning.The remaining constraints originate from STAARCH and are used to enforce some disturbance related physical properties: recently cleared stands are normally characterized by an increase in brightness, a decrease in NDVI and a decrease in wetness [23].
As it is likely that there are a large number of close misclassifications (i.e., pixels with values slightly smaller than one of the chosen thresholds), we expanded this approach by permutating the five thresholds (DIT, ΔDIT, ΔBr,T, ΔWr,T and ΔNDVIr,T) with random numbers for the normal distribution.For each threshold, 100 random numbers are generated (the mean μ of the normal distribution is defined by the user and σ = 0.5). Figure 5 illustrates the generation of random numbers for the example of DIT.Randomly generated thresholds are shown in Figure 5a for each permutation and Figure 5b shows the normally-distributed histogram.A random number vector is generated for each threshold and they are subsequently combined to form a multi-dimensional threshold confidence space.A three-dimensional example is shown in Figure 5c, where darker colors indicate values that are closer to the 3-dimensional centroid.We generated a larger amount of combinations for illustration purposes only, to better visualize this three-dimensional threshold confidence space.In our practical implementation, we generate such a point cloud in the five-dimensional threshold space with 100 permutations.
Consequently, the multi-temporal detection (above rule set) is performed with 100 different threshold sets.The count of a pixel being potentially disturbed at time t is considered as the disturbance probability P(t) = [0,100].A moving maximum filter of width w is applied to the P(t) time series, in order to only retain the most certain disturbance in a given period and a disturbance is eventually found if P(t) exceeds a given threshold PT. (iv) Spatial filtering.After identifying disturbance events in the temporal sequence, each disturbance mask (binary image; 1 = disturbed, 0 = stable forest) at time t is spatially filtered.Consequently, the multi-temporal detection (above rule set) is performed with 100 different threshold sets.The count of a pixel being potentially disturbed at time t is considered as the disturbance probability P(t) = [0,100].A moving maximum filter of width w is applied to the P(t) time series, in order to only retain the most certain disturbance in a given period and a disturbance is eventually found if P(t) exceeds a given threshold P T .
(iv) Spatial filtering.After identifying disturbance events in the temporal sequence, each disturbance mask (binary image; 1 = disturbed, 0 = stable forest) at time t is spatially filtered.Disturbed pixels are rejected if none of their eight neighbors are disturbed in order to reduce outlier-induced noise [23].Undisturbed pixels are set to disturbed if at least five neighbors are disturbed because they are likely also disturbed.Finally, the disturbance masks are segmented and each disturbed object with an area smaller than a given threshold n min is removed.This filter is suitable for further noise reduction, but has the contradicting effect on the identification of small and narrow linear disturbances.The latter two filters originate from the Fmask cloud and cloud shadow detection algorithm [37] which in terms of image processing is related to disturbance detection because it is essentially also a binary classification method.
(v) Hybrid detection.We additionally implemented a hybrid detection that makes use of observed and synthesized data.The Landsat detection result is used as the starting point and provides irregularly spaced disturbance masks.For each identified disturbance, the part of the synthetic time series between the potential event and the foregoing date (denoted as time slice) is investigated in detail.Analogous to the standalone detections, the (i) DI transformation; (ii) temporal filtering; and (iii) permutation-based detection of disturbances is performed for the synthetic images in the time slice.If there is a synthetic image with a higher disturbance probability P(t), this date is recorded instead.For simplicity, we assume that only one disturbance may occur per time slice.
The disturbance detection approach was applied to: (1) the Landsat time series; (2) the synthetic time series; and (3) a hybrid setup.The parameterization of the Landsat and STARM detection is given in Table 1.The synthetic time-series was processed with smaller DI and P detection thresholds in response to prediction-related uncertainty.The window size for identifying the observation with the highest disturbance probability (w) was increased as a consequence of the increased temporal density.The minimum area threshold was decreased because the synthetic time series was characterized by more noise, both in temporal and spatial terms.Thus, patches were more likely to be split apart into several smaller objects.The hybrid detection used the Landsat detection result as processing mask and the temporal improvement was solely based on the highest disturbance probability of the synthetic time series in a given time slice.Informative areas for a variety of disturbance processes are shown, and the products are qualitatively compared.A quantitative assessment was performed using stratified Poisson Disc Sampling [42], based on the occurrence of disturbances in chips of 5 ˆ5 pixels; as perforated detections in any product prevent the usage of a purely pixel-based assessment (see the examples in the results section).The temporal resolution of the Landsat detection was reduced to the annual SLATS reporting periods prior to the quantification.The overall accuracy, producer's accuracies and user's accuracies for the non-disturbed, 2008 and 2009 disturbance classes were computed.

STARFM Quality Assessment
Figure 6 displays the coefficient of determination R 2 for each observed and predicted image pair as well as the Landsat poor quality percentage.In most cases, R 2 is larger than 0.87 (25% quantile; mean R 2 for all bands is 0.90; mean R 2 for band 4 is 0.88; mean R 2 for band 3 is 0.93).The prediction quality is clearly wavelength-dependent.The SWIR and red bands were predicted with superior quality.The blue band was predicted with least precision and the NIR band was also not perfectly predicted.The green band is in between the blue and the red band.There is little evidence that the prediction was strongly affected by the temporal distance of the base pairs to the target image.Some of the R 2 drops seem to coincide with a high portion of poor quality pixels in the observed Landsat images.
quality.The blue band was predicted with least precision and the NIR band was also not perfectly predicted.The green band is in between the blue and the red band.There is little evidence that the prediction was strongly affected by the temporal distance of the base pairs to the target image.Some of the R 2 drops seem to coincide with a high portion of poor quality pixels in the observed Landsat images.

Landsat Disturbance Detection
In order to test the functionality of the implemented disturbance detection with respect to a number of known disturbance processes, we first applied the approach to the observed Landsat images.The results were compared with the SLATS land cover change dataset [40].The overall accuracy is 83.14%.The producer's accuracy is 97.11% for the non-disturbed class, and 67.43% and 71.51% for 2008 and 2009 detections, respectively.The user's accuracy is 75.85% for the non-disturbed class, and 93.66% and 97.71% for 2008 and 2009 detections, respectively.
Clear-cut events recorded in SLATS are also detected with our approach.The area depicted in Figure 7 was mapped as two-stage clearing (2008: blue, 2009: yellow) (see Figure 7a).Our approach (Figure 7b) reveals that there were actually three distinct clearings ( Figure 8 depicts an area where we identified more patches (Figure 8b) compared to SLATS (Figure 8a).The two road clearings (16 September 2008 and 27 March 2009) were detected by both approaches, though the southern feature was identified imperfectly in our approach.There is one particular date (7 February 2009: green) at which excessive disturbances were mapped in Figure 8b.DI time series are shown in Figure 8g for fire and regrowth clearing; marked as ( 1)-( 2) in Figure 8b.A tremendous DI increase is evident for both pixels.Multispectral imagery for the pre-and post-disturbances are shown in Figure 8c-f.

Landsat Disturbance Detection
In order to test the functionality of the implemented disturbance detection with respect to a number of known disturbance processes, we first applied the approach to the observed Landsat images.The results were compared with the SLATS land cover change dataset [40].The overall accuracy is 83.14%.The producer's accuracy is 97.11% for the non-disturbed class, and 67.43% and 71.51% for 2008 and 2009 detections, respectively.The user's accuracy is 75.85% for the non-disturbed class, and 93.66% and 97.71% for 2008 and 2009 detections, respectively.
Clear-cut events recorded in SLATS are also detected with our approach.The area depicted in Figure 7 was mapped as two-stage clearing (2008: blue, 2009: yellow) (see Figure 7a).Our approach (Figure 7b) reveals that there were actually three distinct clearings ( Figure 8 depicts an area where we identified more patches (Figure 8b) compared to SLATS (Figure 8a).The two road clearings (16 September 2008 and 27 March 2009) were detected by both approaches, though the southern feature was identified imperfectly in our approach.There is one particular date (7 February 2009: green) at which excessive disturbances were mapped in Figure 8b.DI time series are shown in Figure 8g for fire and regrowth clearing; marked as ( 1)-( 2) in Figure 8b.A tremendous DI increase is evident for both pixels.Multispectral imagery for the pre-and post-disturbances are shown in Figure 8c-f.Figure 9 displays an area that was affected by fire and/or clearing.Analogous to Figure 8, the fire (7 February 2009: green) was detected in our approach (Figure 9b).but is not evident in SLATS (Figure 9a).The fragmented clearing (1 July 2009: blue) was detected by both approaches, though some smaller features are missing in the Landsat detection.Multispectral imagery for pre-and post-disturbance are shown in Figure 9c-f.DI time series are shown in Figure 9g for representative pixels in the burnt area, in the cleared area and in an area affected by both disturbance agents.The cleared and burnt pixels were successfully detected and are characterized by a single increase in DI, whereas the twofold disturbance was undetected; two more subtle increases in DI are evident for this pixel.Figure 9 displays an area that was affected by fire and/or clearing.Analogous to Figure 8, the fire (7 February 2009: green) was detected in our approach (Figure 9b).but is not evident in SLATS (Figure 9a).The fragmented clearing (1 July 2009: blue) was detected by both approaches, though some smaller features are missing in the Landsat detection.Multispectral imagery for pre-and post-disturbance are shown in Figure 9c-f.DI time series are shown in Figure 9g for representative pixels in the burnt area, in the cleared area and in an area affected by both disturbance agents.The cleared and burnt pixels were successfully detected and are characterized by a single increase in DI, whereas the twofold disturbance was undetected; two more subtle increases in DI are evident for this pixel.Figure 9 displays an area that was affected by fire and/or clearing.Analogous to Figure 8, the fire (7 February 2009: green) was detected in our approach (Figure 9b).but is not evident in SLATS (Figure 9a).The fragmented clearing (1 July 2009: blue) was detected by both approaches, though some smaller features are missing in the Landsat detection.Multispectral imagery for pre-and post-disturbance are shown in Figure 9c-f.DI time series are shown in Figure 9g for representative pixels in the burnt area, in the cleared area and in an area affected by both disturbance agents.The cleared and burnt pixels were successfully detected and are characterized by a single increase in DI, whereas the twofold disturbance was undetected; two more subtle increases in DI are evident for this pixel.

STARFM Disturbance Detection
We applied the disturbance detection algorithm to the synthetic time series.Figure 10 depicts the previously introduced subsets demonstrating the strengths and weaknesses of this approach.
There is more detection noise when using synthetic data (compare Figure 10d,e, blue box).There is more temporal variability within the patches opposed to the Landsat results.Some patches are characterized by noise in the disturbance timing (e.g., red box in Figure 10b), whereas other patches are more clearly separated into regions of distinct clearing stages (yellow box in Figure 10b).
Rather small patches are more strongly affected by the temporal intermixing: for example, the burn scar (green box in Figure 10e

STARFM Disturbance Detection
We applied the disturbance detection algorithm to the synthetic time series.Figure 10 depicts the previously introduced subsets demonstrating the strengths and weaknesses of this approach.
Remote Sens. 2016, 8, 277 13 of 22 Figure 10.Disturbance detection for the previously introduced subsets; see Figure 7 for (a-c), Figure 8 for (d-f), Figure 9 for (g-i).The Landsat detection, STARFM detection, and hybrid detection is shown in columns 1-3, respectively.See Figure A1 for the Landsat random color map and Figure A2 for the STARFM and hybrid random color map.(j) Temporal profiles of the Disturbance Index for the pixel marked in (g,h).
As a consequence of the spatial deficiencies in the STARFM detection (i.e., more noise, problems at edges and omission of small structures) but promising additional temporal information, we developed and applied the hybrid detection approach to the Landsat and STARFM time series.Results of the hybrid detection are shown in the third column in Figure 10 and are similar to the Landsat and STARFM results.The detection is characterized by coherent patches as in the Landsat detection, and also features the temporal richness of the STARFM detection.
Figure 11 displays histograms of the temporal improvement for the three largest time slices in the Landsat time series (see Figure 4a).All cleared pixels in the respective time slices were used to generate the histograms.The temporal improvement using the synthetic images is measured relative to the Landsat detection.In the following, this shift in disturbance timing is simply referred to as an Figure 10.Disturbance detection for the previously introduced subsets; see Figure 7 for (a-c), Figure 8 for (d-f), Figure 9 for (g-i).The Landsat detection, STARFM detection, and hybrid detection is shown in columns 1-3, respectively.See Figure A1 for the Landsat random color map and Figure A2 for the STARFM and hybrid random color map.(j) Temporal profiles of the Disturbance Index for the pixel marked in (g,h).
There is more detection noise when using synthetic data (compare Figure 10d,e

, blue box).
There is more temporal variability within the patches opposed to the Landsat results.Some patches are characterized by noise in the disturbance timing (e.g., red box in Figure 10b), whereas other patches are more clearly separated into regions of distinct clearing stages (yellow box in Figure 10b).
Rather small patches are more strongly affected by the temporal intermixing: for example, the burn scar (green box in Figure 10e) is mainly decomposed to three different dates with the largest time difference being 78 days (blue: 24 November 2008; violet: 10 October 2008; brown: 10 February 2009).The edges and the general appeal of the shapes are less defined in comparison to the Landsat product.This mainly results in the omission of disturbed pixels at the edges of disturbed areas.Similarly, small and thin structures are less pronounced in the STARFM result, which results in a less exact estimation of timing and extent.Nevertheless the shape of the structures is still evident and the objects are not completely undetected, e.g., the fragmented clearing pattern (white box in Figure 10h).DI time series are shown in Figure 10j for the marked point in Figure 10g,h, which reveals that the sharp DI increase in the discrete Landsat imagery was interpolated as gradual change by STARFM.
As a consequence of the spatial deficiencies in the STARFM detection (i.e., more noise, problems at edges and omission of small structures) but promising additional temporal information, we developed and applied the hybrid detection approach to the Landsat and STARFM time series.Results of the hybrid detection are shown in the third column in Figure 10 and are similar to the Landsat and STARFM results.The detection is characterized by coherent patches as in the Landsat detection, and also features the temporal richness of the STARFM detection.
Figure 11 displays histograms of the temporal improvement for the three largest time slices in the Landsat time series (see Figure 4a).All cleared pixels in the respective time slices were used to generate the histograms.The temporal improvement using the synthetic images is measured relative to the Landsat detection.In the following, this shift in disturbance timing is simply referred to as an "improvement" (regardless if this is true or not; no validation data are available to conform if this was indeed an improvement).No or little improvement could be achieved for a number of pixels (red bars), whereas the disturbance timing of many pixels could be improved, either being concentrated on one date (Figure 11a) or indicating several clearing stages (Figure 11b,c).
Remote Sens. 2016, 8, 277 14 of 22 was indeed an improvement).No or little improvement could be achieved for a number of pixels (red bars), whereas the disturbance timing of many pixels could be improved, either being concentrated on one date (Figure 11a) or indicating several clearing stages (Figure 11b,c).This area is show in Figure 12; the same color-coding was applied.
Figure 11.Improvement in disturbance timing using the hybrid detection.The improvement (in 8-day steps) is measured relative to the Landsat detection (dark red bar, i.e., no improvement).The three slices that are characterized by the largest data gaps in the Landsat time series (see Figure 4) are shown in (a-c), respectively.The majority of the disturbances in (c) are attributed to only one disturbed area.This area is show in Figure 12; the same color-coding was applied.

STARFM Quality Assessment
The blue band prediction was of least quality, which was amended by non-clear-sky conditions.Cloud pixels were not used in the quality assessment, but there is a chance that some clouds were missed by the automatic cloud detection methods [38] and/or there was generally a higher aerosol loading.The blue band is strongly affected by atmospheric scattering [46].Accurate measures of aerosol optical depth (AOD) are difficult to obtain in Australia [47] and the Landsat data were corrected by using a fixed and low AOD [36].This might have caused some inaccuracies, especially in images where AOD was actually higher.The influence of AOD decreases exponentially with increasing wavelength and the gradual increase of the prediction quality from blue to green, red and the SWIR bands (where the influence of AOD is negligible) supports the assumption of this being an aerosol-induced effect.As remaining atmospheric contamination in the reference Landsat image is known to decrease quality measures that compare synthetic and real images [15], it is likely that atmospheric scattering was not completely removed by the atmospheric correction of either input source.
The NIR band is an exception to this as it is also affected by comparably low prediction quality.The relative performance of the NIR and the red band was already discussed in the literature with respect to the Landsat radiometric processing level and the environmental setting.A summary of several comparable studies is given in Table 2. Walker et al. [15] found similar values (also using NBAR MODIS imagery) in a dryland forest/woodland area in Arizona.They found similar R 2 values of 0.94 (red) and 0.84 (NIR).Hilker et al. [13] and Watts et al. [14] observed the reverse behavior.Hilker et al. [13] proposed that the relative poorer red estimates are due to a greater effect of atmospheric contamination on shorter wavelengths.Walker et al. [15] additionally noted that the location of [13]'s study area in a temperate Canadian forest would support this assumption as consequence of a generally stronger contaminated atmosphere.This could also explain why R 2 values of [13] are consistently lower (especially in the shorter wavelength bands; red: 0.27 < R 2 < 0.67; NIR: 0.73 < R 2 < 0.82) than in this study and in [15,17].On the contrary, the study area of [14] is located in north-central Montana with semi-arid climate (similar to the climate in [15] and in this study).However, Landsat data in their study were only TOA corrected, thus atmospheric contamination was still present in their data.Bhandari et al. [17] stated that the weaker performance of NIR predictions in a dryland environment, compared to the red band, could also be due to differences in vegetation canopy, branch and trunk structure between broadleaved and coniferous forests, resulting in different BRDF properties at MODIS and Landsat scales.Regardless of the band-specific performance, a strong relationship between synthetic and observed images is generally existent (mean R 2 = 0.90).Moreover, our retrieved R 2 values are approximately in the same dimension or even higher compared to similar studies (e.g., [13,17]).As such, we consider the generated synthetic time series of sufficient quality to be used for further analyses.

Disturbance Index
Reference statistics were taken from a single class (i.e., all woody vegetation pixels) as the area under investigation was assumed to be quite homogeneous with only evergreen forest communities.The DI transformation might fail if both deciduous and evergreen forests are present within the normalization population [27].Nevertheless, the results might be improved by rescaling each forest/woodland class separately.This would be especially useful if this approach would be applied to more heterogeneous landscapes [23].
We used the initial state of our time series to generate the reference population, assuming that the land cover was merely stable.If this approach is to be applied to longer time series, the reference population should be updated occasionally, especially if a great proportion of the reference pixels experienced land cover change.
The Disturbance Index was initially designed to highlight stand replacing disturbances, developed and tested in dark closed-canopy coniferous forests [24,25].Masek et al. [26] state, that the DI transformation performs weaker in sparse forests.Australian forests and woodlands are characterized by trees and shrubs with sparse foliage and irregular crown shapes [5].Despite not being intended for this environment, the disturbance detection algorithm was indeed able to detect stand-replacing disturbances reliably.The more subtle disturbance mechanisms (like selective logging, gradual change or rapid multiple disturbances) were detected imperfectly, but the DI was never intended to do so.This also applies to low FPC stands where the contrast of low biomass forest vs. cleared area can be smaller than the assumed detection thresholds, as shown for some missed fire and logging disturbances in Figure 8.Given that the DI was designed to highlight the transition from dark conifers to bright soils, we recognize that some disturbances may be omitted, although we consider the omission rates acceptable under the given environmental settings.

Landsat Disturbance Detection
In general, the same clear-cut objects were found in the SLATS and the Landsat detection, which points out that both approaches identify stand-removing disturbances reliably.Nevertheless, the increased temporal resolution of our result indicates the feasibility for high temporal resolution time series applications.The spatially more coherent appeal of our approach in comparison to SLATS could be due to the implemented spatial filtering or potentially be caused by the permutation of the detection thresholds.Thresholding techniques with only a simplex of thresholds can cause perforated results if the spectral response of the disturbance is similar to the used threshold.Using multiple threshold sets and creating some sort of confidence area, as implemented in this study, can partially mitigate this effect.
The DI profile of the first clearing in Figure 7g depicts a signal sequence that would potentially not have been detected with an approach scanning for the DI to get beyond a certain hard threshold, as used in STAARCH.The DI is high from the beginning, but nevertheless shows a significant increase at the actual clearing date.This observation necessitated an adaption in the thresholding ruleset compared to STAARCH.
Several disturbance mechanisms were present in our study area, among them extensive clear-cut logging (Figure 7), fragmented clearing (Figure 9), regrowth clearing (Figure 8) and fire (Figure 8).Differences between our detection and the SLATS reference dataset are mainly attributed to the disturbance mechanism.The SLATS assessment features extensive manual editing to isolate the woody vegetation change due to anthropogenic reasons.As such, non-permanent disturbances, e.g., caused by fire, were not available in the reference dataset.The regrowth clearing patches (see Figure 8) were potentially falsely classified as natural disturbance by the SLATS image interpreter.In general, there were many excess patches in one specific time-step in our product.Fires are known to be a common feature in the study area [33] and as they are frequently induced by lightning in combination with the high flammability of eucalyptus species [48], it is very likely that many of these excess patches are burn scars.We conclude that both datasets are generally similar and of high quality, if the different thematic scopes are considered.
More subtle disturbances and multiple disturbances in rapid succession were subject to omission; e.g., pixel (3) in Figure 9b.The effect of the first disturbance (fire) was not strong enough to be identified by our approach or FPC was potentially too low to trigger detection (see also previous section).The second disturbance (chaining-based logging) might have been detected if the forest would not have experienced the first disturbance.The rapid sequence of disturbances prevented the automatic identification as the potentially strong increase in DI was masked by the fire disturbance.The clearing was recorded in the SLATS dataset, which is likely due to the bi-temporal approach, where annual FPC estimates are used.In this special case, a reduction in the temporal resolution was beneficial.
In a binary sense, our algorithm was able to reliably detect disturbances caused by a variety of different processes.Although not attempted in this study, additional spectral information, e.g., fractional ground cover estimates [49], may be utilized in a future version of this algorithm to approach a disturbance type classification (clear cut logging, regrowth clearing, selective harvesting, fire, etc.).
The quantitative detection accuracy was of medium quality, and although the statistics are not of highest quality, visual investigation points to a successful disturbance detection (as discussed in the previous paragraphs).Based on the user's and producer's accuracies, false classifications only occur between the non-disturbed and disturbed classes, but not within the disturbed classes, i.e., years.This points to a successful identification of the disturbance timing, though spatial and thematic differences between the products may likely have caused omissions and commissions of disturbances.Thus, several factors have decreased the quantitative assessment: the SLATS data was produced with a completely different methodology, being based on an automatic bi-temporal change detection of annual FPC estimates.In a second step, extensive manual editing was performed to ensure quality and to remove non-anthropogenic change (like fire).Thus, there are disturbances that are not reported in SLATS.In addition, manual quality checking ensures that detection noise-a common feature of automatic approaches-is strongly decreased.On the other hand, manual editing is a subjective process, and may result in false removals.Subtle disturbances were subject to omission in our approach, and potentially, the manual quality checking of SLATS ensured that these patches were correctly identified by their approach.The different temporal resolutions might also have influenced the results, where, e.g., gradual change like selective logging might be easier to detect with a bi-temporal annual approach as changes from time to time are small, but are high in larger time steps.
Due to the aforementioned discussion, we conclude that the implemented algorithm produced reasonable results with irregularly spaced Landsat data for most disturbance mechanisms, but mainly has an edge in enabling the disturbance history assessment with higher temporal resolution.Although it would be preferable to also edit our results manually, we assume that this fully automatic algorithm can be regarded as functional, and processing synthetic data with high temporal resolution can be approached.In addition, the disturbance detection is in good agreement with the SLATS dataset (or was found to be reasonable in cases of discrepancies) and is therefore suited to assess the quality of the disturbance products that were generated with synthetic imagery.

(i)
Synthetic time series stand-alone detection.The higher temporal diversity in the STARFM product could indicate that there were indeed several clearing stages with high temporal frequency.This would be a substantial asset for state legislation as a detailed temporal documentation of the clearing succession is needed in cases of illegal logging [50].On the contrary, the partially noisy appeal of disturbance timing could rather point to uncertainties in the STARFM prediction.It is hard to assess which process is true or if both factors are evident.In a local temporal window around each clearing event, the latter process seems to be certainly true, though.The boundaries of adjacent dates are somewhat diffuse and tend to become blurred and mixed up within one patch (e.g., Figure 10b red box; 25 June 2008: dark green, 3 July 2008: light-green, 11 July 2008: pink).This implies that the synthetic images might be useful to detect disturbance/clearing patches, but not at the full temporal resolution.One option might be to declare a detection uncertainty in temporal terms, e.g., ˘8 days or ˘16 days.Secondly, one could aggregate the results to a lower temporal frequency, e.g., to fortnightly or monthly time steps, which in many cases would still be an improvement to sparse Landsat data.On the contrary, there seem to be disturbed patches that could indeed be of different dates (e.g., Figure 10b yellow box; green: 9 June 2008; pink: 11 July 2008), because there is a sharp temporal and spatial boundary between the two polygon parts.As no reference data were available at an appropriate temporal scale, this proposition cannot be falsified or verified at the moment.
Nevertheless, this effect appears to be related to patch size (and maybe even patch shape or complexity, as a square shape is more likely to appear less mixed with its surrounding in a coarse resolution MODIS observation) whether the distribution to different dates could be true or just a data-blending artifact.For example, the burn scar (green box) in Figure 10e is strongly affected by noise in the disturbance timing with large temporal differences.This effect is being caused by diffuse spatio-temporal blending artifacts in the synthetic images, which is likely due to the inability of STARFM to identify pure homogeneous coarse resolution pixels for that rather small area.Therefore the observed change in the bracketing Landsat base images and the intermixed reflectance change in the MODIS image of the prediction days (the stable forest is changing as well due to phenology) cannot be exactly assigned to the correct location of the change.This results in some kind of unpredictable spatial assignment of the change, visible as a blurring with the surrounding.
The edges of changing shapes were identified to be less defined in comparison to the Landsat product.According to Gao et al. [10], differences are indeed greatest at the border of shapes, because the spatial resolution of the STARFM prediction cannot exactly match that of the high resolution image.This results in some blending or blurring of the object of interest with its surroundings.Therefore, small structures that are comprised of a large proportion of border in relation to core area or object edges in general are subject to larger prediction errors or blending artifacts.This is evident in the underlying synthetic images as well as in the final results (note especially the more rugged shape of the larger disturbances in Figure 10b compared to Figure 10a).
Although the disturbance algorithm is somewhat inaccurate in the detection of small-scale structures (especially the exact date and at the outline of the shape), the shape of structures is nevertheless evident in the results.This supports the visual observation that the spatial detail of the Landsat scenes was conserved in the synthetic time series appropriately (with respect to the known limitations of STARFM).The striping pattern in Figure 10h (white box) was at least partly detected with synthetic images.It was not detected to the full extent, as there occurred some spatial blending of the surrounding areas (subsidiary) and especially blending in the time domain, transitioning from an undisturbed to a disturbed state, as can be seen in Figure 10j.STARFM is capable of estimating the transition from one state to another [10] which is a double-sided feature.On the one hand, it enables the detection of a disturbance at all-if the disturbance is at least recorded in one of the base pairs [10].On the other hand, it would be desirable if the change would be captured as a pronounced distinct change event.Instead, the disturbance signal was predicted as gradual change, which somewhat limits the applicability of a STARFM-based disturbance detection.

(ii)
Hybrid detection.The hybrid detection was characterized by coherent patches as in the Landsat detection, but additionally favored from increased temporal detail of the synthetic time series.Some pixels were characterized by no or little shift in disturbance timing.Thus, the timing either could not be improved with dense synthetic images or the discrete Landsat images were captured at the optimal point in time by chance.Figure 12 displays the major disturbances for the time slice shown in Figure 11c.The non-improved pixels are merely situated at the edge of patches or are small objects, which indicate that STARFM was not able to extract useful information (regarding land cover change) from the subtle MODIS reflectance change at a coarse spatial scale-as discussed in Figure 12i.Some Landsat images were not used for processing because of strong overall cloud contamination (see Figure 4).Fortunately, parts of the disturbance were visually recognizable as can be seen in Figure 12b,c.A fire front is visible for patch B (red box), but not for A and C, which were obviously burnt some time before.The left-out images indicate that patch B was burnt on 5 December 2008 and patches A and C burnt anytime between this date and 3 November 2008.The hybrid detection dates the patches A and C primarily to 24 November 2008 (blue) and subsidiary to 10 December 2008 (green), whereas patch B is dated to 10 December 2008 (green), 18 December 2008 (light-yellow) and 26 December 2008 (yellow), which is in quite good agreement with the left-out Landsat observations.As such, the general decomposition of the disturbance dates seems to be reasonable and underpins the potential applicability of improving the temporal resolution of fire history assessments using dense synthetic imagery.The estimated temporal accuracy of the hybrid approach in comparison to the manually revealed approximate dating in the Landsat poor quality images also backs the aforementioned suggestion of a detection uncertainty of 8-16 days or alternatively an aggregation of the results to coarser time steps.
We consider the hybrid approach superior to the standalone Landsat and STARFM disturbance detections, as the hybrid approach benefits from the general quality of the Landsat detection and its well-defined spatial characteristics, as well as the temporal improvement in disturbance timing for sufficiently large shapes as a consequence of Landsat/MODIS data-fusion.

Conclusions
We successfully implemented a binary disturbance detection approach capable of processing irregularly spaced time series and/or equidistant dense synthetic imagery.It was tested whether disturbances caused by a variety of known processes could be detected.Process-based stand-replacing change, e.g., extensive clearing and fire, was reliably identified in observed Landsat imagery.Subtle changes were detected imperfectly as a consequence of the employed Disturbance Index not being intended for gradual change.Multiple disturbances in rapid succession-resulting in relatively subtle single changes-were also problematic due to the multi-temporal approach.Permutating the detection thresholds allowed for a more robust parameterization and a spatially more coherent detection.The results were in good agreement with reference data, or were found to be reasonable in cases of discrepancies.Differences were mainly identified to be dependent on the disturbance mechanism as a result of different product scopes.The standalone usage of the synthetic time series revealed deficiencies with regards to noise and spatio-temporal inaccuracies at the edge of disturbances, as well as higher uncertainty regarding small objects.Nevertheless, the additional temporal information proved to be promising and resulted in the development of a hybrid detection approach that utilized the assets of both data sources, i.e., the spatial coherence and general accuracy of the Landsat detection and the temporal improvement capabilities inherent in the synthetic time series-limited by known object-size related data-blending limitations.

Figure 1 .
Figure 1.Study area: The Queensland study area is displayed as the white rectangle.The corresponding full Landsat frame (WRS-2 path/row 093/078) is indicated by the dark grey box.The largest settlements in and near to the study area are displayed.

Figure 2 .
Figure 2.Broad vegetation groups of the forested study area (Regional Ecosystem Mapping of 2006[30]).The analysis was restricted to the forested areas (colors); the "other" class was compiled using

Figure 1 .
Figure 1.Study area: The Queensland study area is displayed as the white rectangle.The corresponding full Landsat frame (WRS-2 path/row 093/078) is indicated by the dark grey box.The largest settlements in and near to the study area are displayed.

Figure 1 .
Figure 1.Study area: The Queensland study area is displayed as the white rectangle.The corresponding full Landsat frame (WRS-2 path/row 093/078) is indicated by the dark grey box.The largest settlements in and near to the study area are displayed.

Figure 2 .
Figure2.Broad vegetation groups of the forested study area (Regional Ecosystem Mapping of 2006[30]).The analysis was restricted to the forested areas (colors); the "other" class was compiled using series; (ii) the synthetic time series; and (iii) to a hybrid setup.Detection results are compared to reference data and to each other.

Figure 3 .
Figure 3. Schematic workflow of this study with references to the sections and sub-sections.

Figure 3 .
Figure 3. Schematic workflow of this study with references to the sections and sub-sections.

Figure 4 .
Figure 4. Landsat (a) and MODIS (b) input images for generating the dense synthetic Landsat-like time series.A high quality (HQ) observation is an observation, which passed the quality tests.Landsat images with more than two-third HQ pixels observations were used for the prediction.The three largest data gaps are indicated by arrows.

Figure 4 .
Figure 4. Landsat (a) and MODIS (b) input images for generating the dense synthetic Landsat-like time series.A high quality (HQ) observation is an observation, which passed the quality tests.Landsat images with more than two-third HQ pixels observations were used for the prediction.The three largest data gaps are indicated by arrows.

Figure 5 .
Figure 5. Illustration of the threshold permutation with random numbers for the normal distribution: (a) sequentially generated random numbers; (b) histogram of the random numbers; and (c) three-dimensional example of the combined threshold confidence area for a large number of permutations (1000).Darker colors represent data points that are located in proximity to the confidence centroid.

Figure 5 .
Figure 5. Illustration of the threshold permutation with random numbers for the normal distribution: (a) sequentially generated random numbers; (b) histogram of the random numbers; and (c) three-dimensional example of the combined threshold confidence area for a large number of permutations (1000).Darker colors represent data points that are located in proximity to the confidence centroid.

Figure 6 .
Figure 6.Coefficient of determination R 2 , inferred from linear regressions between observed Landsat images and synthetic equivalent imagery.The Landsat poor quality percentage is shown as vertical bars.
16 September 2008: blue; 27 March 2009: violet; 28 April 2009: beige).The shapes are very similar in Figure 7a,b, with better defined edges in Figure 7b.The three clearing stages are clearly visible in the underlying multispectral imagery (see Figure 7c-f).DI time series are shown in Figure 7g for three representative pixels marked as (1)-(3) in Figure 7b.Large DI increases-up to 4 × ΔDI(t)-are evident at the clearing events.

Figure 6 .
Figure 6.Coefficient of determination R 2 , inferred from linear regressions between observed Landsat images and synthetic equivalent imagery.The Landsat poor quality percentage is shown as vertical bars.
16 September 2008: blue; 27 March 2009: violet; 28 April 2009: beige).The shapes are very similar in Figure 7a,b, with better defined edges in Figure 7b.The three clearing stages are clearly visible in the underlying multispectral imagery (see Figure 7c-f).DI time series are shown in Figure 7g for three representative pixels marked as (1)-(3) in Figure 7b.Large DI increases-up to 4 ˆ∆DI(t)-are evident at the clearing events.

Figure 7 .
Figure 7. SLATS reference dataset, Landsat disturbance detection result for a clearcut area: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for three selected pixels marked in (b).

Figure 8 .
Figure 8. SLATS reference dataset, Landsat disturbance detection result for an area affected by regrowth, clearing and burning: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for two selected pixels marked in (b).

Figure 8 .
Figure 8. SLATS reference dataset, Landsat disturbance detection result for an area affected by regrowth, clearing and burning: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for two selected pixels marked in (b).

Figure 8 .
Figure 8. SLATS reference dataset, Landsat disturbance detection result for an area affected by regrowth, clearing and burning: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for two selected pixels marked in (b).

Figure 9 .
Figure 9. SLATS reference dataset, Landsat disturbance detection result for an area affected by clearing and burning: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for three selected pixels marked in (b).
) is mainly decomposed to three different dates with the largest time difference being 78 days (blue: 24 November 2008; violet: 10 October 2008; brown: 10 February 2009).The edges and the general appeal of the shapes are less defined in comparison to the Landsat

Figure 9 .
Figure 9. SLATS reference dataset, Landsat disturbance detection result for an area affected by clearing and burning: (a) SLATS forest clearing (blue: 2008; yellow: 2009); (b) Landsat disturbance detection (see Figure A1 for the random color coding); (c-f) Landsat images depicting the clearing stages (RGB = NIR; SWIR1, red); and (g) temporal profiles of the Disturbance Index for three selected pixels marked in (b).

Figure 11 .
Figure11.Improvement in disturbance timing using the hybrid detection.The improvement (in 8-day steps) is measured relative to the Landsat detection (dark red bar, i.e., no improvement).The three slices that are characterized by the largest data gaps in the Landsat time series (see Figure4) are shown in (a-c), respectively.The majority of the disturbances in (c) are attributed to only one disturbed area.This area is show in Figure12; the same color-coding was applied.

Figure 11 .
Figure 11.Improvement in disturbance timing using the hybrid detection.The improvement (in 8-day steps) is measured relative to the Landsat detection (dark red bar, i.e., no improvement).The three slices that are characterized by the largest data gaps in the Landsat time series (see Figure 4) are shown in (a-c), respectively.The majority of the disturbances in (c) are attributed to only one disturbed area.This area is show in Figure12; the same color-coding was applied.

Figure 12 .Figure 12 .
Figure 12.Hybrid detection and observed Landsat images: (a) hybrid detection, the improvement in disturbance timing relative to the Landsat detection was color-coded as in Figure 11c; and (b-d) Figure 12.Hybrid detection and observed Landsat images: (a) hybrid detection, the improvement in disturbance timing relative to the Landsat detection was color-coded as in Figure 11c; and (b-d) observed Landsat images (RGB = NIR; SWIR1, red).The 3 November 2008 and 5 December 2008 images were not used in the Landsat detection, nor were they used in generating the dense synthetic time series due to poor overall quality.(b) Pre-disturbance image; (c) Disturbance in progress; and (d) Post-disturbance image.

Figure A1 .
Figure A1.Color table for the Landsat disturbance detection results.

Figure A2 .
Figure A2.Color table for the STARFM and hybrid disturbance detection results.

Table 1 .
Disturbance detection parameterization for the standalone Landsat and STARFM detections, as well as the hybrid detection, which follows the parameterization of the standalone detections.

Table 2 .
Literature-based comparison of STARFM red and NIR band prediction performance as a function of environmental setting and Landsat radiometric processing level.