Comparative study of three satellite image time-series decomposition methods for vegetation change detection

ABSTRACT Satellite image time-series (SITS) methods have contributed notably to detection of global change over the last decades, for instance by tracking vegetation changes. Compared with multi-temporal change detection methods, temporally highly resolved SITS methods provide more information in a single analysis, for instance on the type and consistency of change. In particular, SITS decomposition methods show a great potential in extracting various components from non-stationary time series, which allows for an improved interpretation of the temporal variability. Even though many case studies have applied SITS decomposition methods, a systematic comparison of common algorithms is still missing. In this study, the seasonal trend loess (STL), breaks for additive season and trend (BFAST) and multi-resolution analysis-wavelet transform (MRA-WT) were explored in order to evaluate their performance in modelling, monitoring and detecting land-cover changes with pronounced seasonal variations from simulated normal difference vegetation index time series. The selected methods have all proven their ability to characterize the non-stationary vegetation dynamics along with different physical processes driving the vegetation dynamics. Our results indicated that BFAST is the most accurate method for the examined simulated dataset in terms of RMSE, whereas MRA-WT showed a great potential for the extraction of multi-level vegetation dynamics. Considering the computational efficiency, both STL and MRA-WT outperformed BFAST.


Introduction
Vegetation change detection from remotely sensed images with high temporal resolution is important to promote better decisions for sustainable land management. The free availability of satellite image time series (SITS) datasets like Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS) and more recently Sentinel makes the development and analysis of change detection techniques an important topic in the field of remote sensing (Tewkesbury, Comber, Tate, Lamb, & Fisher, 2015). Several change detection approaches that specifically target SITS have been developed. The field of change detection started with approaches for bi-temporal datasets, including image differencing (Bruzzone & Prieto, 2000;Mas, 1999), change vector analysis (Bovolo, Marchesi, & Bruzzone, 2012), principal component analysis (Byrne, Crapper, & Mayo, 1980) or tasselled cap transformation-based methods (Minu & Shetty, 2015). Although these methods are still very popular and can also be applied to more than two images, they may be impractical for long and dense SITS and may suffer from some limitations. For example, bi-temporal images carry the risk of confounding temporary and persistent changes like interpretation of a temporary land-cover change as permanent transformation of land-cover type. Furthermore, changes may be induced by sensor change (e.g. Landsat-5 TM vs. Landsat-7 ETM+), by arbitrary thresholds and control parameters of the change-detection techniques or by disregarding the non-stationary data particularity caused by seasonal and trend variations (Coppin, Jonckheere, Nackaerts, Muys, & Lambin, 2004;Verbesselt, Hyndman, Newnham, & Culvenor, 2010a). High temporal-resolution SITS methods account for severalalthough not allof these drawbacks and, hence, are an interesting tool for tracing complex trajectories of land-cover change.
Most frequently, SITS methods are applied in combination with vegetation indices. For example, the normal difference vegetation index (NDVI) has been extensively used to trace vegetation dynamics (Fensholt & Proud, 2012;Hall-Beyer, 2003;Tucker & Sellers, 1986). The NDVI has proven to be a proxy for the status of the above-ground biomass at the landscape level for its high correlation with green-leaf density, net primary production and CO 2 fluxes (Tucker & Sellers, 1986;Running & Nemani, 1988;Wylie et al., 2003). Overall, NDVI time series are usually non-stationary due to the presence of different frequencies that reflect seasonal, gradual and abrupt vegetation changes (Martínez & Gilabert, 2009;Verbesselt, Hyndman, Zeileis, & Culvenor, 2010b). This study focuses on different approaches for decomposing NDVI time series into these components of vegetation dynamics. Such methods are widespread in the remote-sensing field for at least two reasons. Firstly, they are able to characterize the non-stationary vegetation dynamics. Secondly, they help to study the physical processes driving the vegetation dynamics. Among the different techniques used to monitor vegetation changes from multi-temporal SITS, two groups of classification methods were proposed (Huang et al., 1998;Hamilton, 1994), namely, spectral-frequency approaches and statistical approaches.
To the best of our knowledge, a systematic comparison of such methods is still missing. In this context, this study explores three popular time-series decomposition methods in order to compare their features and the particularity of their results in modelling, monitoring and detecting simulated NDVI time-series changes of vegetated land cover with pronounced seasonal variations. The three selected methods are the spectral approach, MRA-WT, and two statistical approaches: the STL and BFAST. We applied a simulated time series in this study in order to avoid unknown variation in the data. The limits and the advantages of each approach were determined using accuracy assessment statistics.

Simulated data
The analysis based on simulated time series aimed to robustly evaluate three decomposition approaches in a controlled environment. NDVI time series were simulated at 16-day temporal resolution based on the methodology described in  and improved with components that characterize real SITS. Figure 1 shows examples of three simulated time series with different timing and magnitude of change, mean NDVI values and abrupt changes' values for both seasonal and trend components. Additionally, different noise levels (0.01-0.08 σ) were added to both dormant and growing seasons.

Methodology
We simulated multiple time series by varying the parameters (see Figure 1 for examples). Then each of the datasets was analysed with three time-series decomposition methods and the different temporal components were determined. The values of the obtained components were then compared to the known values of the simulated datasets and the root mean square error (RMSE) was calculated. For a quantitative evaluation, the average RMSE was computed using approximately 300 iterations for the simulation of time series by varying parameter combinations. Hereafter, brief background information and the implementation parameters of the three methods compared in this study are described.

Season trend loess (STL)
STL is a seasonal trend decomposition algorithm. It decomposes the time-series signal into three components: trend, seasonal and remainder (Brandt et al., 2014). It is based on a locally weighted regression smoother filter developed by Cleveland, Cleveland, McRae, and Terpenning (1990). Two recursive loop procedures are implemented in order to estimate and update the seasonal and the trend component. The inner loop consists of six steps: (1) detrending, (2) cycle sub-series smoothing, (3) low-pass filtering of smoothed cycle subseries, (4) detrending of smoothed cycle sub-series, (5) deseasonalizing and (6) trend smoothing. The outer loop passes through the inner loop and computes robustness weights (Cleveland et al., 1990). As a consequence, STL is flexible enough to reveal nonlinear patterns and is robust to outliers (Jacquin et al., 2010). Additionally, it handles any type of seasonality (Jiang et al., 2010). STL assumes that the decomposition model is additive, such that (Equation (1) where S t is the decomposed NDVI time series, Ct t is the trend component, Cs t is the seasonal component and Ca t is the remaining variation after the estimation of trend and seasonal components, named random component. Six inputs must be chosen in order to decompose the data (Cleveland et al., 1990;Jiang et al., 2010;Sanchez-Vazquez et al., 2012): the frequency of observation (periodicity) Δt, the number of robustness iterations n o (depending on the presence of aberrant behaviour being equal to 0 if the data are stable), the number of inner loop's passes n i (a recommended value equal to 2 is selected), the low-pass filter's smoothing n l parameter (greater than or equal to n p ), the seasonal component's smoothing parameter n s (the recommended value is greater than or equal to 7 and near to the periodic length) and the trend component's smoothing parameter n t (ranging between 1.5n p and 2n p ).
Generally, some of these parameters already have recommended values while others are directly related Multi-resolution analysis-wavelet transform MRA-WT Uses a wavelet that looks at different components of the time series with windows adjusted according to scale. (Martínez & Gilabert, 2009) Empirical mode decomposition EMD Gets an ensemble of adaptive orthogonal components called intrinsic mode functions. (Kong et al., 2015) Season trend loess STL Decomposition of time series into three components: trend, seasonal and remainder (Jacquin et al., 2010) Break for additive trend and season BFAST Decomposition of time series into three components using an additive model (Verbesselt et al., 2010a) Detecting breakpoints and estimating segments in trend DBEST Uses a novel segmentation algorithm that simplifies the trend into linear segments. to the input time series. In our case, Δt = 23 images per year. The STL parameters were determined according to the characteristics of the NDVI time series since the original NDVI time series are easily affected by noise. Finally, the STL parameters for the NDVI time series were defined as follows: n p = 23 samples/year, n i = 2, n o = 4, n 1 = 24 and n s = 23.

Break for additive season and trend (BFAST)
BFAST is a generic algorithm for detecting the structural change in SITS (Jamali et al., 2015;Verbesselt et al., 2010aVerbesselt et al., , 2010b. It integrates the decomposition of time series into three components, namely, a trend component, a seasonal component and a remainder component (Watts & Laffan, 2014). The first estimation of the seasonal component uses the STL algorithm. In addition, with the embedded methods of detecting structural changes within the trend and the seasonal components, BFAST is able to detect the number of changes and their dates (Hutchinson et al., 2015). The first of these methods is the ordinary least squares residuals-based moving sum (OLS MOSUM) test (Zeileis, 2005). It determines the occurrence or absence of statistically significant abrupt changes. The second one is the Bayesian information criterion (Bai & Perron, 2003), which determines the optimal number and dates of abrupt changes (Forkel et al., 2013;Hutchinson et al., 2015;Verbesselt et al., 2010a). The user-defined parameters are the number of images of the time series (N), the model seasonality (i.e. "harmonic" or "dummy"), a parameter noted h that determines the minimum time interval between breakpoints (BPs), the number of iterations I and the maximum number of abrupt changes (max BP ) (Hutchinson et al., 2015). BFAST handles only the additive decomposition model (Equation (1)). The trend component is represented by a piecewise linear function which is in turn used to represent both linear and nonlinear variations (Verbesselt et al., 2010a). The seasonal component is modelled by a harmonic function using a finite sum of trigonometric terms (Verbesselt et al., 2010b). BFAST can be used to analyse various types of SITS and can be applied to other disciplines dealing with seasonal or non-seasonal series, such as hydrology, climatology and econometrics (Verbesselt et al., 2010a).
In this study, the harmonic model for the seasonal component was chosen following the recommendation by Verbesselt et al. (2010b). With respect to max BP , several values were tested until the detected number of BP remained unchanged. For the parameter h, we maintained the recommended value used in De Jong, Verbesselt, Schaepman and De Bruin (2012) and Verbesselt et al. (2010a).

Multi-resolution analysis-wavelet decomposition (MRA-WT)
The multi-resolution analysis (MRA) is devoted to a discrete signal as it decomposes the signal into m decomposition levels by successively translating and convolving the elements of a high-pass filter and a low-pass scaling filter associated with the mother wavelet implemented in a hierarchical algorithm. This process is known as the pyramid algorithm (Mallat, 1989). The filters retain the small-and large-scale components of the signals, also known as detail (D) and approximation (A) series. The scaling parameter is expressed as 2 to the power of a = 2 j , j = 1,..,m, where m is the highest decomposition level.
The detail component (D m (t)) reflects the periodic fluctuations at different scales resulting through highpass filtering where the period is related to the m level. However, the sum of the detail component from 2 to the level m describes the total inter-annual variability (V(t)). The first step is the noise reduction level. Consequently, it is eliminated from the total seasonal variability. It pursues the rapid changes caused by high frequencies. Besides, the inter-annual variation is modelled by the approximation component (A m (t)), which is a low-pass filtered component. The original signal can be recovered by applying In fact, the connection between the scale and the frequency of WT is known by period. The period p is correlated with the scale a, the sampling period Δt and the centre frequency of the mother wavelet as follows (Equation (3)): The period refers to the global temporal scale over the wavelet. In this study, ΔT is equal to 16 days and v c = 0.6634. This latter is given by the frequency of Meyer wavelet as recommended by Martínez and Gilabert (2009). According to Equation (3), p/ 2 = 385 days for a = 2 5 days. Hence, A 5 describes a long-term variability of the semi-period of 385 days. V = Σ i D i (i = 2,.., 5) describes the total inter-annual variability of the time series. The four different detail components D 2 (24-48 days), D 3 (49-96 days), D 4 (97-192 days) and D 5 (193-385 days) roughly correspond to monthly, bimonthly, seasonal and semiannual variations, respectively. In order to interpret and highlight the major properties of the studied canopy using the MRA-WT components, we calculated various NDVI metrics. These metrics basically reflect the vegetation phenology and rely on the statistical functions applied to one or more of these NDVI time series: the inter-annual variability, the intra-annual variability, and the detail component that reflects the semi-annual variations (Martínez & Gilabert, 2009). NDVI mean indicates the average amount of vegetation, NDVI min depicts the minimum NDVI's level, ΔNDVI reflects the amplitude of the annual phenological cycle and t max is a temporal metric. This latter is an indicator of the timing of the maximum NDVI. For this purpose, the maximum of the semi-annual variation is identified. In our study, it is presented by the fourth detail component (D 4 ). Figure 2 shows the trend components obtained by the decomposition of simulated time-series examples presented in Figure 1(a-c). Figure 3 illustrates how the decomposition methods fit the seasonal variability compared to the model seasonality for the same simulated time series (Figure 1(a-c)). The results depicted in Figures 2 and 3 are summarized in three characteristics as follows:

Results
• Capturing the details: According to the simulated dataset, the most detailed trend component was provided by the STL method followed by the WT inter-annual variation (Figure 2). Both contain small portions of seasonal variation that may lead to the false interpretation of the trend component variability. In contrast, BFAST estimated one linear segment trend with a small negative slope (in the case of (a) the simulated time series) or a succession of linear segments with good BP detection (in the case of (b) and (c)). However, the details of the seasonal component were better estimated by BFAST and MRA-WT approaches (black and red lines in Figure 3). • Boundary effect: As shown in Figure 2 case (a), both STL and MRA-WT trends diverge. The first goes down and the second goes up close to the end, which demonstrates the boundary effect. As a solution, the SITS methods were replicated at the beginning and the end of the series, obtaining a better fitting in the other two simulated time series (b) and (c). • Influence of sharp changes: The change observed in the simulated seasonal component (a) was only detected by BFAST and MRA-WT decompositions. In contrast, STL decomposition tracks out this dip in the trend component and not in the seasonal one. Additionally, V t shows shortduration deflections after the BP. Its impact spreads over all levels of wavelet coefficients affecting the seasonal coefficients (Figure 3(a,c)). The lowest average RMSE for both trend and seasonal estimated components was obtained using the BFAST approach (Table 2) followed by MRA-WT and STL methods. The STL method assumes an invariant seasonal component, which may explain this behaviour. Despite that, the original trend had  The analyses of the three methods have been performed on a personal computer with one 2.50-GHz CPU and 6.0-GB memory. In this case, STL and MRA-WT are relatively fast methods for the simulated time series with an average execution time of 354 s for STL, 500 s for MRA-WT and 677 s for BFAST. As such, STL was the fastest of the three methods.

Discussion
In this paper, we have explored three methods (BFAST, STL and MRA-WT) to decompose a nonstationary NDVI time series. These methods produced different accuracies according to the simulated NDVI time series. BFAST gave the most accurate results in terms of RMSE, whereas MRA-WT proved to be more informative due to the multi-level decomposition (details of components from D 2 to D 5 ). We assume that this is mainly due to the hierarchical extraction strategy of multiple frequency dynamics within the time series. STL was by design not able to detect abrupt changes in the seasonal component of the simulated NDVI time series, because of the static seasonal component embedded within it.
Based on the obtained results, we made the following observations. Firstly, the decomposition nature of the three methods varies under different conditions, but they are all suitable for the decomposition of NDVI time series. This resulted in comparable decomposition results. Generally, BFAST and MRA-WT described the long-term behaviour almost identically, although the BFAST method gave better RMSE values in the decomposition of the simulated data compared to STL and MRA-WT. It appeared able to detect changes with good temporal precision. Moreover, it appeared appropriate for the estimation of the trend component since it provides detailed information about the trend, including the number of BPs, their magnitude and the corresponding date (Browning, Maynard, Karl, & Peters, 2017). However, the BFAST model is computationally slow compared to other algorithms. Furthermore, it would benefit from additional tools for noise component investigation. In fact, BFAST implements STL for the initial estimation of the seasonal component, for which it is slower by design. The runtime of satellite image time-series decomposition methods is an important issue because the remote-sensing community deals with large datasets. A faster implementation of BFAST would therefore be of advantage. Alternatively, the MRA-WT algorithm gives a very useful multi-level decomposition that contains more information about the plant growing cycle.  According to a recent study investigating the sensitivity of multi-scale temporal vegetation dynamics to climatic variability, different components react differently to climatic changes with various correlation levels (Qiu et al., 2016). The MRA-WT decomposition method can yield smoother seasonal component profiles than the other two models. Also, the seasonal component generated by the MRA-WT model seems to be most reasonable and natural, in terms of reflecting the dynamics in the main vegetation stages, such as the start and the end of the season. Nevertheless, a deep research has to be done to enhance the semi-period (p/2) of each component (Percival & Walden, 2006).
The STL model is very versatile and robust. For this reason, it is widely used in the economic field. However, our results showed a limited efficiency in decomposing NDVI time series compared to other tested models.
Overall, the three models generated similar multiscale temporal dynamics. Precisely, BFAST described the general BPs in both seasonal and trend components, while MRA-WT provided more informative multi-temporal components. STL maintained a greater part of the variability in the remainder component; however, it reflected the overall changes well and had computational advantages.
The following ideas can be considered for future studies: (1) A more in-depth study of the seasonal variability could be done using temporal and NDVI-derived metrics. Therefore, in a follow-up study, we will aim to explain the interesting total variability of MRA-WT dynamics.
Our field survey could contribute to such a study by improving the interpretation of species-specific temporal dynamics, especially in forested area. (2) The residual component of the BFAST and STL decomposition represents the inherent noise (Jamali et al., 2015) and is often ignored by researchers. Nevertheless, this noise may encompass valuable data, which has already been proven in others fields (Grieser, Trömel, & Schönwiese, 2002). Therefore, we recommend the use of statistical tools such as the extreme events theory to make use of and interpret this component.
(3) Decomposition of NDVI time series with high variability into sub-series allows more detailed analysis of potential drivers of changes compared to the original vegetation index time series. A variety of decomposition models has been developed and used in the literature. However, it remains an active research topic and new models continue to be proposed such as DBEST (Jamali et al., 2015). Then, the comparison could be extended to include other algorithms.

Conclusion
NDVI time series are usually non-stationary since they present different frequency components, such as seasonal variations, long-and short-term fluctuations, which may not be limited to the mean of the series but also affect its overall variance structure. In this context, the aim of this work was to examine and compare the effectiveness of three time-series decomposition methods (MRA-WT, STL and BFAST) that are appropriate for non-stationary time-series analysis. The study was conducted using a simulated 16-day NDVI time series. Based on the RMSE values, the best results of decomposition were given by BFAST and MRA-WT with small difference between them. However, we acknowledge that the results may vary for different sites or with different parameter sets required by these methods. The STL method, on the other hand, has computational advantages and may be the suitable choice for analysis of large datasets and in circumstances where trend breaks are not anticipated.