Monitoring the Vegetation Dynamics in the Dongting Lake Wetland from 2000 to 2019 Using the BEAST Algorithm Based on Dense Landsat Time Series

: The dynamic monitoring and analysis of wetland vegetation play important roles in revealing the change, restoration and reconstruction of the ecosystem environment. The increasing availability of high spatial-temporal resolution remote sensing data provides an unprecedented opportunity for wetland dynamic monitoring and change detection. Using the reconstructed dense monthly Landsat time series, this study focuses on the continuous monitoring of vegetation dynamics in Dongting Lake wetland, south China, in the last two decades (2000–2019) by using the Bayesian estimator of abrupt change, seasonal change, and trend (BEAST) method. Firstly, the ﬂexible spatiotemporal data fusion (FSDAF) model is applied to blend Landsat and moderate-resolution imaging spectroradiometer (MODIS) images on the basis of the input image pair selection strategy named “cross-fusion” to generate the monthly time-series normalized di ﬀ erence vegetation index (NDVI) with the spatial resolution of 30 m. Then, the abrupt changes, trend, and seasonality of the vegetation in the study area as well as the uncertainties of change detection are estimated by the BEAST method. Results show that there is a close relationship between the ground true data and the estimated changepoints. A high overall accuracy (OA) of 87.37% and Kappa coe ﬃ cient of 0.85 were achieved by the proposed framework. Additionally, the temporal validation got the interval intersection of 86.57% and the absolute di ﬀ erence of mean interval length of 6.8 days. All of the results demonstrate that the vegetation changes in the Dongting Lake wetland varied spatially and temporally in the last two decades, because of extreme weathers and anthropogenic factors. The presented approach can accurately identify the vegetation changes and time of disturbance in both the spatial and temporal domains, and also can retrieve the evolution process of wetland vegetation under the inﬂuence of climate changes and human activities. Therefore, it can be used to reveal potential causes of the degradation and recovery of wetland vegetation in subtropical areas.


Introduction
Wetland is an open, complex and interactive component of the ecosystem [1], as well as the co-occurrence of environmental services and biological diversity [2,3]. However, the wetland vegetation Appl. Sci. 2020, 10 is extremely sensitive to climate and environmental changes. In the last several decades, due to the global climate change and human activities, the wetland ecosystem has faced many degradation risks, such as area shrinkage [4], ecological degradation [5], invasion of alien species [6] and pollution [7]. Therefore, a long-term dynamic monitoring of wetland vegetation is of practical significance and is important and essential to sustainable development.
With the increasing availability of high spatial-temporal resolution data, remote sensing technology provides spatially comprehensive information for monitoring and detecting vegetation disturbance and land cover changes [8][9][10]. Coarse spatial resolution data, such as MODIS and advanced very high-resolution radiometer (AVHRR) images, have been used to monitor the surface change in regional or large area, because of the high temporal resolution and wide spatial coverage. For example, the multi-temporal MODIS NDVI 16-day composite data have been used to detect the land cover change in the Albemarle-Pamlico Estuary System (APES) region in the US [11]; the MODIS imagery was employed to monitor the land cover change in Borneo from 2002 to 2005 [12]; the AVHRR data were used in monitoring the land cover dynamics in Bangladesh [13]. However, these moderate-resolution images caused mixed pixels. Furthermore, these images cannot be used to extract the detailed land cover changes or vegetation disturbance, due to their low spatial resolution. Remote sensing images of a 30-m resolution can identify most vegetation disturbances and land cover changes [14]. Thus, the free Landsat archive with a spatial resolution of 30 m provides an opportunity for the continuous dynamic monitoring of land cover changes. Landsat time series at high temporal resolution have been used to monitor the areas of land change at a regional scale, which has accurately reflected the forest [15], wetland [16] and coastline [17] changes in the long term. In addition, the Landsat time series could provide information on the ecosystem cover and condition, which helps the understanding of the ecosystem function, resilience and dynamics [18]. Nevertheless, the research of the long-term evolution of the wetland ecosystem in subtropical regions has been seriously hampered by the insufficient high-quality and long-term continuous wetland change monitoring data, because of frequent cloudy and rainy weather.
To get high spatial and temporal resolution data, a number of Spatial and Temporal satellite Image Fusion (STIF) models have been developed to blend the MODIS (coarse resolution and short interval) with Landsat images (finer resolution and long interval), such as the Flexible Spatiotemporal Data Fusion (FSDAF) model [19], the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) model [20], the enhance STARFM (ESTARFM) model [21] and the Robust Adaptive Spatial and Temporal Image Fusion model [22]. The STARFM has been one of the most widely used spatial-temporal fusion models, which generates high spatial and temporal resolution images at the prediction time based on the input image pairs on the base date. It considers the influence of the spatial distance and spectral similarity of neighboring pixels, as well as the temporal variation [23]. But the STARFM may cause prediction errors in heterogeneous landscapes with rapid land surface dynamics (e.g., Dongting Lake wetland), as the base images of the STARFM cannot capture transient or sudden changes of landcover. Besides, the STARFM requires the temporal variation information of pure pixels for fusion. If there is no pure pixel in the moving window, the fusion quality will be poor. The ESTARFM and its modified version (i.e., mESTARFM) have been proposed to reduce the prediction errors induced by searching spectrally similar neighbor pixels in the STARFM [21,23]. However, the second problem results in collinearity in the equation system, which makes the fused images cannot be used in continuous dynamics monitoring of heterogeneous wetland. The FSDAF model integrates unmixing based methods, spatial interpolation, and the STARFM into one framework, which has better capability for monitoring dynamic disturbance in heterogeneous area and provides a new possibility for long-term dynamic monitoring [19].
How to accurately and effectively identify the vegetation disturbance with different durations and spatial distribution is challenging and critical. Conventionally, post-classification comparison [24][25][26][27], image differencing- [28] and transformation-based methods [29], were used to monitor the vegetation disturbance and land cover change based on one pair of remote sensing images across a large timescale.
Most of these methods entail high quality images and complex data preprocessing, which cannot be always available in the vegetation monitoring, because of the seasonal phenology of vegetation. Furthermore, these methods ignore ecosystem changes and temporal context and regularity, as they extract the changes by comparing two time-phase images. Therefore, many temporal decomposition techniques have been proposed to track the temporal behaviors of vegetation and account for the abrupt, seasonal and trend changes in the ecosystem. For instance, the Landsat-based detection of trends in disturbance and recovery (LandTrendr) method uses Landsat time series data to segment the trend signals for change detection based on a piecewise linear model [30]. Breaks for additive seasonal and trend (BFAST) [31] and detecting breakpoints and estimating segments in trend (DBEST) [32] algorithms were used successively to extract the periodic variations and trend of time series and identify the abrupt changes by the Seasonal-Trend decomposition procedure based on Loess. Moreover, other monitoring approaches such as the vegetation change tracker (VCT) [33], continuous change detection and classification (CCDC) method [34], and exponentially weighted moving average change detection (EWMACD) [35] were proposed. The VCT algorithm reconstructs recent forest disturbance history based on the spectral-temporal properties of land cover and forest change processes. The EWMACD algorithm uses the exponentially weighted moving average chart and Landsat images to detect the multitemporal forest change and can achieve an accuracy of 85% [35]. By the CCDC approach, land cover change can be detected when the difference between the observed and predicted (based on a time series model) pixels exceeds a threshold three times consecutively and then the kinds of land cover change could be identified by the CCDC method (based on random forest classifier). Although these methods have been widely used, they have some defects. First, these methods should preset a range of parameters. Second, some of these algorithms cannot estimate the seasonal changes. Third, these methods ignore the real problems in ecosystem evolution, because they only adopt the single best-fitting model (i.e., linear model or piecewise linear model) in a signal estimation. Foremost, the single best-fitting model may lead to opposite results that vegetation greening is suggested by linear model but this trend only reflects in a previous segmentation of time series when use piecewise linear modeling [36][37][38]. To solve these problems, the BEAST method was proposed to utilize a Bayesian-ensemble-based approximator to fit arbitrarily complex trends of time series, which can quickly and accurately characterize individual components and credible uncertainty within the time series [39]. To some extent, this ensemble strategy makes the method more ecologically meaningful and is beneficial to the mechanistic interpretations of the terrestrial ecosystem and climate-biosphere interactions as well as the monitoring performance [39].
The STIF and time-series decomposition algorithm were often used separately, and have been seldomly used for long-term monitoring in wetland areas. On one hand, the cloudy and rainy weather in the subtropical regions make the acquisition of continuous remote sensing data with high spatial and temporal resolution difficult, which hampers the studies related to wetland monitoring [17,40]. On the other hand, previous time-series algorithms focus on seeking a single best model for change detection, but ignore the interpretation of the ecosystem change mechanism, so they may result in controversy because of the opposite results [16,41]. In this research, we integrate the FSDAF model with the BEAST algorithm to continuously monitor the 20-year vegetation changes of the Dongting Lake wetland in south China, for a better understanding of the long-term evolution process of wetland ecosystem under global warming. The rest of this paper is structured as follows: Section 2 briefly introduces the Dongting Lake wetland and the dataset used in this study. Section 3 describes how to fuse time series using the FSDAF model based on the cross-fusion strategy and how the BEAST algorithm is applied to detect the abrupt changes, trend, and seasonality of wetland vegetation. Section 4 demonstrates the fusion and change detection results. Section 5 discusses the implication of the results and some ideas for further work. In the final Section, some conclusions are presented.

Study Area
The Dongting Lake wetland, located in central-south China, is a part of the Yangtze River basin. There are two national nature reserves (i.e., Xi Dongting Lake (Mupinghu) nature reserve and east Dongting Lake nature reserve) in the Dongting Lake wetland. As the second largest freshwater lake wetland in China, the Dongting Lake wetland plays a significant role in local ecology protection and economic development. In this research, the east Dongting Lake wetland (Latitude/Longitude: 28 • 59 N to 29 • 38 N, 112 • 43 E to 113 • 15 E, 12,000 km 2 ) was selected as the study area ( Figure 1). This region has the subtropical monsoon climate with an annual precipitation of 1300 mm and an average temperature of 17 • C. The wetland area is a sedimentary plain with an elevation of 20-45 m. It is composed of fluvial and lacustrine deposits, estuarine deltas as well as outer lakes.

Study Area
The Dongting Lake wetland, located in central-south China, is a part of the Yangtze River basin. There are two national nature reserves (i.e., Xi Dongting Lake (Mupinghu) nature reserve and east Dongting Lake nature reserve) in the Dongting Lake wetland. As the second largest freshwater lake wetland in China, the Dongting Lake wetland plays a significant role in local ecology protection and economic development. In this research, the east Dongting Lake wetland (Latitude/Longitude: 28°59′ N to 29°38′ N, 112°43′ E to 113°15′ E, 12,000 km 2 ) was selected as the study area ( Figure 1). This region has the subtropical monsoon climate with an annual precipitation of 1300 mm and an average temperature of 17 °C. The wetland area is a sedimentary plain with an elevation of 20-45 m. It is composed of fluvial and lacustrine deposits, estuarine deltas as well as outer lakes. With sufficient moisture and heat, the Dongting Lake wetland is a habitat of biogeographic and ecological interest, and has many endemic, rare, and relict species [42]. The wetland vegetation in the Dongting Lake wetland includes herbs (dominated by reeds and sedges), economic crops (dominated by single season rice and double cropping rice), deciduous broad-leaved forests (dominated by planted Populus nigra forests) and evergreen broad-leaved mixed forests [43,44]. During the past 20 years, the vegetation structure has changed frequently, resulting in the increasing instability of the ecological environment [43]. Besides, the combination influence of global warming and human factors increases the sedimentation rate, which results in the reduction of the lake area and ultimately decreases biological diversity [43,45]. In the past, double cropping rice accounts for the majority of the paddy rice in this region. With the drop in grain prices, a lot of double-cropping rice paddies have been converted to single-season rice paddies, so single-season rice dominates the paddy rice nowadays [46]. In addition, some paddy fields and bottomland have been planted poplar, because of the high market demand since 1990 [45]. However, since 2011, the poplar forests have been intensely deforested to plant native tree species for comprehensive ecological improvement [47,48]. Under the With sufficient moisture and heat, the Dongting Lake wetland is a habitat of biogeographic and ecological interest, and has many endemic, rare, and relict species [42]. The wetland vegetation in the Dongting Lake wetland includes herbs (dominated by reeds and sedges), economic crops (dominated by single season rice and double cropping rice), deciduous broad-leaved forests (dominated by planted Populus nigra forests) and evergreen broad-leaved mixed forests [43,44]. During the past 20 years, the vegetation structure has changed frequently, resulting in the increasing instability of the ecological environment [43]. Besides, the combination influence of global warming and human factors increases the sedimentation rate, which results in the reduction of the lake area and ultimately decreases biological diversity [43,45]. In the past, double cropping rice accounts for the majority of the paddy rice in this region. With the drop in grain prices, a lot of double-cropping rice paddies have been converted to single-season rice paddies, so single-season rice dominates the paddy rice nowadays [46]. In addition, some paddy fields and bottomland have been planted poplar, because of the high market demand since 1990 [45]. However, since 2011, the poplar forests have been intensely deforested to plant native tree species for comprehensive ecological improvement [47,48]. Under the coupling effects of factors mentioned above, the dynamics of wetland vegetation in Dongting lake changes dramatically.

Remote Sensing Data and Preprossing
The NDVI can accurately capture and reflect the growth and nutrition information of the vegetation, and is used to characterize the growth conditions of wetland vegetation for detecting the vegetation dynamic changes in this research. For acquiring the NDVI time series of the last 20 years, all available standard Level 1 terrain corrected and orthorectified Landsat TM/ETM +/OLI Top of Atmosphere (TOA) reflectance images (Path/Row: 124/40, available pixels ≥ 30%), and the MODIS 13Q1 NDVI products (MODIS tile number: h27, v06) covering the whole study area from 2000 to 2019 were downloaded from the USGS website (https://glovis.usgs.gov/).  Figure 2. The preprocessing steps of the Landsat TOA reflectance images include: (1) convert TOA reflectance data to surface reflectance via radiometric calibration and FLAASH atmospheric correction, (2) generate monthly Landsat images using simple synthesis algorithm and quality bands of the Landsat image acquisition [49], (3) detect and mask the remaining cloudy area in monthly Landsat images using the Fmask algorithm [50], (4) calculate NDVI using the B3 and B4 bands of Landsat TM/ETM + sensor and the B4 and B5 of Landsat OLI sensor, (5) clip images using the study area vector.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 24 coupling effects of factors mentioned above, the dynamics of wetland vegetation in Dongting lake changes dramatically.

Remote Sensing Data and Preprossing
The NDVI can accurately capture and reflect the growth and nutrition information of the vegetation, and is used to characterize the growth conditions of wetland vegetation for detecting the vegetation dynamic changes in this research. For acquiring the NDVI time series of the last 20 years, all available standard Level 1 terrain corrected and orthorectified Landsat TM/ETM +/OLI Top of Atmosphere (TOA) reflectance images (Path/Row: 124/40, available pixels ≥ 30%), and the MODIS 13Q1 NDVI products (MODIS tile number: h27, v06) covering the whole study area from 2000 to 2019 were downloaded from the USGS website (https://glovis.usgs.gov/).  Figure 2. The preprocessing steps of the Landsat TOA reflectance images include: (1) convert TOA reflectance data to surface reflectance via radiometric calibration and FLAASH atmospheric correction, (2) generate monthly Landsat images using simple synthesis algorithm and quality bands of the Landsat image acquisition [49], (3) detect and mask the remaining cloudy area in monthly Landsat images using the Fmask algorithm [50], (4) calculate NDVI using the B3 and B4 bands of Landsat TM/ETM + sensor and the B4 and B5 of Landsat OLI sensor, (5) clip images using the study area vector. A total of 460 geo-rectified MODIS NDVI images with the spatial resolution of 250 m and temporal resolution of 16 days were downloaded. Firstly, the maximum value composite method was used to generate the monthly MODIS NDVIs and reduce the cloud contaminated pixels [51]. The monthly MODIS NDVI data were re-projected into the WGS 84/UTM zone 49 N projection and resampled to 30 m by the nearest neighbor method to match the Landsat images. Then, a Landsat image was selected as the reference to register and clip the MODIS NDVIs. Moreover, the Savitzky-Golay filter with the window size of 4 was used to optimize the low-quality pixels in MODIS NDVI time series, which can keep the shape and width of the original signal while filter the noise [8,9]. Therefore, the bad quality data in MODIS NDVI can be eliminated and the real dynamic disturbance can be preserved as much as possible. Finally, the high-quality, monthly MODIS NDVI time series was acquired for subsequent spatial-temporal fusion. A total of 460 geo-rectified MODIS NDVI images with the spatial resolution of 250 m and temporal resolution of 16 days were downloaded. Firstly, the maximum value composite method was used to generate the monthly MODIS NDVIs and reduce the cloud contaminated pixels [51]. The monthly MODIS NDVI data were re-projected into the WGS 84/UTM zone 49 N projection and resampled to 30 m by the nearest neighbor method to match the Landsat images. Then, a Landsat image was selected as the reference to register and clip the MODIS NDVIs. Moreover, the Savitzky-Golay filter with the window size of 4 was used to optimize the low-quality pixels in MODIS NDVI time series, which can keep the shape and width of the original signal while filter the noise [8,9]. Therefore, the bad quality data in MODIS NDVI can be eliminated and the real dynamic disturbance can be preserved as much as possible. Finally, the high-quality, monthly MODIS NDVI time series was acquired for subsequent spatial-temporal fusion.

Reference Data
The reference dataset was collected by National Forestland Change Survey (NFCS) of China in 2018 and TimeSync-plus v4.6 tool [52]. The NFCS produces a forestland database annually for recording forest management activities and the spatial-temporal changes of forestland utilization and land cover types by two steps: (1) visual interpretation based on the remote sensing images with the resolution higher than 5 m (e.g., GF-1, GF-2 and ZY-3) and (2) verifying the results by field survey of each change plot (accuracy over 95%). The changes of the forestland larger than 667 m 2 (smaller than one Landsat pixel) caused by the following reasons are included in the database: afforestation and regeneration, forest harvesting, planning adjustment, requisition and occupation, deforestation, fire, geologic hazard, other disasters and other natural factors. A total of 8960 regions with changes in the study area were extracted from the NFCS database of 2018. In addition to NFCS data, 420 90 × 90 m 2 sampling plots (3 × 3 Landsat pixels) with a spacing of 5 km were extracted by TimeSync based on time-series Landsat and Google Earth images for the validation of the detected changepoints in the last 20 years. The spatial and temporal changes information of every sample plot was recorded monthly. These changes were classified as vegetation degradation (browning), recovery (greening) and stability. Figure 3 depicts the distribution of sampling plots and the forestland changes in 2018 over the study area.

Reference Data
The reference dataset was collected by National Forestland Change Survey (NFCS) of China in 2018 and TimeSync-plus v4.6 tool [52]. The NFCS produces a forestland database annually for recording forest management activities and the spatial-temporal changes of forestland utilization and land cover types by two steps: (1) visual interpretation based on the remote sensing images with the resolution higher than 5 m (e.g., GF-1, GF-2 and ZY-3) and (2) verifying the results by field survey of each change plot (accuracy over 95%). The changes of the forestland larger than 667 m 2 (smaller than one Landsat pixel) caused by the following reasons are included in the database: afforestation and regeneration, forest harvesting, planning adjustment, requisition and occupation, deforestation, fire, geologic hazard, other disasters and other natural factors. A total of 8960 regions with changes in the study area were extracted from the NFCS database of 2018. In addition to NFCS data, 420 90 × 90 m 2 sampling plots (3 × 3 Landsat pixels) with a spacing of 5 km were extracted by TimeSync based on time-series Landsat and Google Earth images for the validation of the detected changepoints in the last 20 years. The spatial and temporal changes information of every sample plot was recorded monthly. These changes were classified as vegetation degradation (browning), recovery (greening) and stability. Figure 3 depicts the distribution of sampling plots and the forestland changes in 2018 over the study area.

Auxiliary Data
The 30-m Shuttle Radar Topography Mission Digital Elevation Model downloaded from USGS website (https://glovis.usgs.gov/) was used to analyze the spatial-temporal vegetation dynamics from

Auxiliary Data
The 30-m Shuttle Radar Topography Mission Digital Elevation Model downloaded from USGS website (https://glovis.usgs.gov/) was used to analyze the spatial-temporal vegetation dynamics from the elevation aspect. The annual average temperature and precipitation data over the Dongting Lake wetland were applied to investigate the relation between climate changes and vegetation variation.

Methods
This research attempts to propose a strategy for the long-term monitoring of the wetland ecosystem and detecting the abrupt wetland vegetation changes under global warming using the fused observations. The method consists of three main steps: (1) collecting and preprocessing raw Landsat and MODIS observations of the period 2000-2019 to predict high spatial and temporal resolution time series; (2) blending Landsat and MODIS NDVI images to obtain the monthly cloud-free fused time series (30-m resolution) over the study area using the FSDAF algorithm; (3) using the BEAST method to decompose the NDVI time series, trend and seasonality of wetland vegetation over the east Dongting Lake, before analyzing the vegetation dynamics from the aspects of duration, distribution, elevation and climate. The flowchart of the proposed approach is shown in Figure 4.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 24 the elevation aspect. The annual average temperature and precipitation data over the Dongting Lake wetland were applied to investigate the relation between climate changes and vegetation variation.

Methods
This research attempts to propose a strategy for the long-term monitoring of the wetland ecosystem and detecting the abrupt wetland vegetation changes under global warming using the fused observations. The method consists of three main steps: (1) collecting and preprocessing raw Landsat and MODIS observations of the period 2000-2019 to predict high spatial and temporal resolution time series; (2) blending Landsat and MODIS NDVI images to obtain the monthly cloudfree fused time series (30-m resolution) over the study area using the FSDAF algorithm; (3) using the BEAST method to decompose the NDVI time series, trend and seasonality of wetland vegetation over the east Dongting Lake, before analyzing the vegetation dynamics from the aspects of duration, distribution, elevation and climate. The flowchart of the proposed approach is shown in Figure 4.

Generation of Monthly Cloud-Free NDVI Time Series Based on the Cross-Fusion Method and FSDAF Algorithm
In order to suppress the influence of cloudy and rainy weather on wetland vegetation dynamics detection, the Flexible Spatiotemporal Data Fusion (FSDAF) method was adopted to generate the consecutive monthly NDVI time series. The FSDAF method blends the coarse-and fine-resolution images at time t1 with the coarse-resolution image at the prediction time t2 to predict a fine-resolution image at t2. Five main steps are included in the FSDAF method: (1) classify the fine-resolution image at t1 using the iterative self-organizing data analysis technique algorithm (an unsupervised classification algorithm: cluster the similar pixels based on distance function to accomplish image classification) [53]; (2) detect the temporal change of each class in the coarse-resolution images between t1 and t2; (3) predict the fine-resolution image at t2 and calculate residuals at each coarse pixel; (4) distribute the temporal prediction residuals to fine pixels by the Thin Plate Spline interpolator; (5) obtain the final prediction of fine-resolution image at t2. Detailed descriptions of FSDAF can be found in [19].
The performance of the fused data is greatly affected by the selected input images [54]. Some studies selected the base image pair acquired at the closest date to the prediction image date to ensure

Generation of Monthly Cloud-Free NDVI Time Series Based on the Cross-Fusion Method and FSDAF Algorithm
In order to suppress the influence of cloudy and rainy weather on wetland vegetation dynamics detection, the Flexible Spatiotemporal Data Fusion (FSDAF) method was adopted to generate the consecutive monthly NDVI time series. The FSDAF method blends the coarse-and fine-resolution images at time t 1 with the coarse-resolution image at the prediction time t 2 to predict a fine-resolution image at t 2 . Five main steps are included in the FSDAF method: (1) classify the fine-resolution image at t 1 using the iterative self-organizing data analysis technique algorithm (an unsupervised classification algorithm: cluster the similar pixels based on distance function to accomplish image classification) [53]; (2) detect the temporal change of each class in the coarse-resolution images between t 1 and t 2 ; (3) predict the fine-resolution image at t 2 and calculate residuals at each coarse pixel; (4) distribute the temporal prediction residuals to fine pixels by the Thin Plate Spline interpolator; (5) obtain the final prediction of fine-resolution image at t 2 . Detailed descriptions of FSDAF can be found in [19]. The performance of the fused data is greatly affected by the selected input images [54]. Some studies selected the base image pair acquired at the closest date to the prediction image date to ensure the highest similarity between them. Some research determined the input base images by the correlation coefficient between the coarse images acquired on the base date and that on the prediction date [55,56]. The "cross-fusion" algorithm takes both similarity and correlation criteria into consideration, which is superior to other input image pair selection methods [54]. In order to improve the fusion result, the cross-fusion method was exploited to quantitatively consider the consistency between the input images acquired at the base time and that acquired at the prediction time, which can automatically determine the optimal input image pairs for FSDAF fusion. In the cross-fusion, a number of base image pairs were selected as candidates from all available images in the light of similarity criterion. All candidate base image pairs were used to predict the fused image by weighted average.
In this study, we generated a monthly cloud-free NDVI time series (30-m resolution) by cross fusing Landsat and MODIS images using the FSDAF algorithm. Then, the fusion results were applied to enhance the data credibility and reconstruct NDVI time series for dynamic monitoring.

Decomposition Model of BEAST
In this research, the BEAST approach was utilized to decompose the reconstructed NDVI time series to obtain abrupt changes as well as seasonal and trend signals for monitoring the dynamic disturbances of wetland vegetation [36]. By the BEAST algorithm, the data space (y) is transformed to the ecological or biophysical space with interspersed abrupt changes over time (t) and remainder (ε) by the decomposition modelŷ (t i ) = f (t i ; Θ) (Equation (1), f (·) donates the function of the decomposition model).
where i is the index of time dimension, Θ s and Θ T represent the parameters related to seasonality (S) and trend (T), respectively, and ε is the residual.
In the BEAST algorithm, the definition of abrupt changes, trend and seasonality is broader than that in other decomposition algorithms. Seasonality (S (t)) is considered as a seasonal signal embedded with disturbance that indicates the phenological change, and it is fitted by a piecewise harmonic model (Equation (2)).
where p is the time-series period and the p separates time scale into k segments; [ξ k , ξ k+1 ] is the interval of the (p + 1)-th segment; t represents the time; L k is the harmonic order at segment k, a k, l and b k, l specify the parameters for sins and cosines, respectively. The changepoints, including all turning points, are deviated from previous regular signals. Trend (T (t)) is defined as a complex nonlinear curve consisting of numerous piecewise lines (the time span was divided by knots (τ j , j = 0, · · · m) into (m + 1) intervals [τ j , τ j+1 ]; a and b are the coefficients of the piecewise linear function) (Equation (3)).
When parameters Θ T and Θ s are assumed as: {Θ T , Θ s } = M, β M , M is the model structure to be estimated (Equation (6)), including the total number and occurrence time of changepoints, harmonic orders of seasonal and trend signals; β M represents the coefficients in M (Equation (7)). Thus, Equation (1) can be rewritten as Equation (8):

Bayesian Ensemble Model of BEAST
According to Equation (8), the vegetation time series can be written asŷ (t i ) = f (t i ; Θ). Conventional models seek the best single linear, piecewise linear, or polynomial model to estimate parameter Θ. BEAST applies general linear models to fit signals and detect abrupt changes based on Bayesian model averaging. Each model in BEAST has a probability to be the true estimation. The Bayesian model averaging (BMA) scenario (Equation (9)) [57] is exploited to parameterize the informative weight for quantifying uncertainties (likelihood P β M , σ 2 , M | y ). The model reduces overfitting and alleviate model misspecification in conventional methods (i.e., single best methods).
where P y | β M , σ 2 , M is the probability of observing data (y)-this probability is simply Gaussian P y | β M , σ 2 , M -β M is the given model's parameters; π β M , σ 2 , M is the prior distribution of the model parameters and it is considered the normal-inverse Gamma distribution; π(M) is the prior on model structure M; σ denotes the data noise. The posterior distribution P β M , σ 2 , M | y includes the information for inferring ecosystem dynamics. The Markov Chain Monte Carlo (MCMC) sampling was utilized to obtain random samples for posterior inference as the intractable analysis of the posterior distribution. MCMC generates various sample-based statistics to test the BMA-based hypotheses. The details of the MCMC method are available in [58]. Additionally, the information on the posterior inference and uncertainties estimation of change detection can be found in [36].

Parameters for BEAST
In general, the key prior parameters for BEAST are easy to determine and most of them depend on the attributes of time-series data (e.g., temporal interval). In our experiment, the maximum number of knots, which is the maximum number of changepoints that can be detected, was set as 8, according to the actual condition of vegetation change in the Dongting Lake wetland; the separation distance in both seasonality and trend was set as 12, as we used monthly data. The minimum harmonic order of seasonality and trend were set as 1 corresponding to the linear model and 0 (which is a constant), respectively. The maximum harmonic orders were set as 6 and 1, respectively. Thus, the trend harmonic order (0 or 1) represents a constant or a linear model, while that of seasonality (range from 1 to 6) can approximate the seasonal signal as linear or nonlinear model arbitrarily well. Based on the BMA strategy, each NDVI pixel (y i ) in a time series was fitted for 10,000 times by MCMC sampling to make the number of samples large enough to ensure the quality of random sampling, and quantify the probability distributions of the number and timing of changepoints.

Validation and Accuracy Assessment
The change points detected by the proposed method were verified by the NFSC data and the manually collected samples from TimeSync-plus tool. Spatially, accuracies of both seasonal and trend change points were evaluated in terms of Kappa coefficients, user's and producer's accuracy (UA and PA) and overall accuracy (OA). The Kappa coefficient evaluates the consistency between the detected change points and validation data. UA and PA determine the commission and omission errors for seasonal and trend change points. OA reflects the ratio of correctly detected change points to the total validation samples. Temporally, we used interval intersection to measure the proximity between the timing of changes and the validation data to manifest the BEAST-derived dates as being close to the actual events.

Spatial-Temporal Vegetation Dynamics Analysis
The vegetation dynamics were analyzed from the aspects of distribution, duration, elevation and climate. The year with the maximum NDVI, timing and the number of abrupt changes estimated by BEAST were used to analyze the vegetation dynamics of the last 20 years. In order to understand the evolution of the wetland system, DEM and climate factor (annual average temperature and precipitation) were considered in the investigation. The proportion of the pixels with single/multiple seasonal and trend changepoint(s) and the proportion of the stable pixels and of the changed pixels were calculated to descript the vegetation changes. Additionally, annual proportion of pixels with the maximum NDVI value and the most significant seasonal and trend changes were extracted from the annual average temperature and precipitation to investigate the relation between vegetation variation and climate and human factors.

Fused Monthly NDVI Time Series
The base image pairs for each month of the last 20 years were selected as candidates to fuse NDVI images by the FSDAF algorithm. The fused NDVI images were used to fill the missing pixels or replace the contaminated pixels of the original Landsat data, but they are not directly utilized for time-series analysis. On one hand, this filling process ensured the best consistency between the derived and the actual data. On the other hand, fusion was only performed in the assigned region (cloudy and bad pixels), which saves the time of FSDAF fusion. Consequently, the final synthetic NDVI image ought to have higher accuracy and quality than the fused image. To avoid confusion, the fusion results were termed as fused images and the final synthetic images were referred to as reconstructed images. The 3D volumetric view of the fused NDVI time series stack in a sub-region of the study area in the past 20 years is shown in Figure 5. The fused image has abundant visual characteristics and details in the heterogeneous region.
In order to assess the consistency of NDVI trajectories between MODIS and the fused NDVI, the NDVI values of three vegetation types (forest, wetland herbs and crop) were calculated and compared. Figure 6 depicts the NDVI time-series relationship between MODIS and the fusion results of the forest, wetland herb and cropland. The trends between MODIS and FSDAF NDVI time series in three vegetation types are similar. The trend of the higher-resolution FSDAF data changes more dramatically than that of the lower-resolution MODIS NDVI in some periods, because higher resolution images capture more detailed landscape changes. The pixel mixing problem of the surface target in coarse-resolution MODIS images may also cause the difference [9]. Obviously, the finer spatial resolution data derived from the FSDAF algorithm provides an opportunity to monitor detailed landscape changes, including less discernable deforestation, degradation or other disturbances.
NDVI image ought to have higher accuracy and quality than the fused image. To avoid confusion, the fusion results were termed as fused images and the final synthetic images were referred to as reconstructed images. The 3D volumetric view of the fused NDVI time series stack in a sub-region of the study area in the past 20 years is shown in Figure 5. The fused image has abundant visual characteristics and details in the heterogeneous region.  In order to assess the consistency of NDVI trajectories between MODIS and the fused NDVI, the NDVI values of three vegetation types (forest, wetland herbs and crop) were calculated and compared. Figure 6 depicts the NDVI time-series relationship between MODIS and the fusion results of the forest, wetland herb and cropland. The trends between MODIS and FSDAF NDVI time series in three vegetation types are similar. The trend of the higher-resolution FSDAF data changes more dramatically than that of the lower-resolution MODIS NDVI in some periods, because higher resolution images capture more detailed landscape changes. The pixel mixing problem of the surface target in coarse-resolution MODIS images may also cause the difference [9]. Obviously, the finer spatial resolution data derived from the FSDAF algorithm provides an opportunity to monitor detailed landscape changes, including less discernable deforestation, degradation or other disturbances. Besides, MODIS and Landsat NDVI images in December 2017 were used to assess the fusion accuracy on the corresponding dates (Figure 7). The fused dense NDVI stack preserves a same spatial resolution with Landsat images and a temporal resolution of one month. Generally, the fusion data show a close correspondence with the actual NDVI image and possess the resembling spatial characteristics of Landsat NDVI. The correlation analysis was also performed between the fused and actual Landsat NDVI in December 2017. As shown in the scatter density plot of Figure 8, high correlation (p < 0.001, RMSE (Root Mean Squared Error) = 0.053, R 2 = 0.93) between the fused and actual NDVI was achieved, despite the land cover change occurred between the base and predicted dates. The results demonstrate that the cross-fusion and the FSDAF method can capture the differences between the base and predicted dates and estimate the interpolated NDVI time series. Therefore, they can be used in the following change detection. Besides, MODIS and Landsat NDVI images in December 2017 were used to assess the fusion accuracy on the corresponding dates (Figure 7). The fused dense NDVI stack preserves a same spatial resolution with Landsat images and a temporal resolution of one month. Generally, the fusion data show a close correspondence with the actual NDVI image and possess the resembling spatial characteristics of Landsat NDVI. The correlation analysis was also performed between the fused and actual Landsat NDVI in December 2017. As shown in the scatter density plot of Figure 8, high correlation (p < 0.001, RMSE (Root Mean Squared Error) = 0.053, R 2 = 0.93) between the fused and actual NDVI was achieved, despite the land cover change occurred between the base and predicted dates. The results demonstrate that the cross-fusion and the FSDAF method can capture the differences between the base and predicted dates and estimate the interpolated NDVI time series. Therefore, they can be used in the following change detection.

Abrupt Changes, Trend, and Seasonality of Wetland Vegetation Derived by the BEAST Method
The BEAST algorithm was used to decompose the NDVI time series into abrupt changes, trend and seasonality components to monitor vegetation dynamics. Seasonality is the periodic variation with abrupt changes caused by intra-annual climatic change or phenological factors. Trend signal is considered as the gradual change resulted from long-term environmental trends, successional dynamics and chronic disturbances, reflecting the stable, browning and greening trend of vegetation. Abrupt changes are related to severe disturbances in both seasonality and trend. Specifically, most disturbance results in the turning of time series, such as flooding, insects, fire, logging, weeding, replantation post-deforestation, urbanization, extreme weather, crop rotation and climate change,

Abrupt Changes, Trend, and Seasonality of Wetland Vegetation Derived by the BEAST Method
The BEAST algorithm was used to decompose the NDVI time series into abrupt changes, trend and seasonality components to monitor vegetation dynamics. Seasonality is the periodic variation with abrupt changes caused by intra-annual climatic change or phenological factors. Trend signal is considered as the gradual change resulted from long-term environmental trends, successional dynamics and chronic disturbances, reflecting the stable, browning and greening trend of vegetation. Abrupt changes are related to severe disturbances in both seasonality and trend. Specifically, most disturbance results in the turning of time series, such as flooding, insects, fire, logging, weeding, replantation post-deforestation, urbanization, extreme weather, crop rotation and climate change,

Abrupt Changes, Trend, and Seasonality of Wetland Vegetation Derived by the BEAST Method
The BEAST algorithm was used to decompose the NDVI time series into abrupt changes, trend and seasonality components to monitor vegetation dynamics. Seasonality is the periodic variation with abrupt changes caused by intra-annual climatic change or phenological factors. Trend signal is considered as the gradual change resulted from long-term environmental trends, successional dynamics and chronic disturbances, reflecting the stable, browning and greening trend of vegetation. Abrupt changes are related to severe disturbances in both seasonality and trend. Specifically, most disturbance results in the turning of time series, such as flooding, insects, fire, logging, weeding, replantation post-deforestation, urbanization, extreme weather, crop rotation and climate change, can be detected as abrupt changes by BEAST. Additionally, the possibility and uncertainty of the distribution of these three components were quantified in the study area. Figure 9 shows an example of an NDVI time-series decomposition in a forest stand. In this example, the BEAST algorithm has detected two seasonal changepoints (Figure 9a) and four trend changepoints (Figure 9b) with the 95% credible interval of occurrence time, which possibly demonstrates the unstable characteristics of forest growth. The probabilities of both the seasonal and trend change points were estimated using the BMA strategy. The threshold of the probability of abrupt changes can be controlled manually to determine the location of change points according to the prior knowledge. BEAST is beneficial from repeating sampling (10,000 times in this experiment) and the BMA method. The standard deviations of fitting in BEAST can be controlled by the BMA strategy (Figure 9e,f). In other words, BEAST improves the accuracy and stability of fitting in seasonality and trend. Figure 9g,h show the estimated harmonic order for fitting seasonal and trend signals, respectively. In Figure 9i,j, the probabilities distribution of seasonal and trend change points are present. In order to understand the evolution of wetland vegetation during the last two decades, the slope of trend with the 95% confident interval and its standard deviations (Std) were measured (Figure 9k,l). Figure 9m depicts the probability of having a positive rate of change in trend that assists in identifying the disturbance types (e.g., deforestation and recovery). can be detected as abrupt changes by BEAST. Additionally, the possibility and uncertainty of the distribution of these three components were quantified in the study area. Figure 9 shows an example of an NDVI time-series decomposition in a forest stand. In this example, the BEAST algorithm has detected two seasonal changepoints (Figure 9a) and four trend changepoints (Figure 9b) with the 95% credible interval of occurrence time, which possibly demonstrates the unstable characteristics of forest growth. The probabilities of both the seasonal and trend change points were estimated using the BMA strategy. The threshold of the probability of abrupt changes can be controlled manually to determine the location of change points according to the prior knowledge. BEAST is beneficial from repeating sampling (10,000 times in this experiment) and the BMA method. The standard deviations of fitting in BEAST can be controlled by the BMA strategy (Figure 9e,f). In other words, BEAST improves the accuracy and stability of fitting in seasonality and trend. Figure 9g,h show the estimated harmonic order for fitting seasonal and trend signals, respectively. In Figure 9i,j, the probabilities distribution of seasonal and trend change points are present. In order to understand the evolution of wetland vegetation during the last two decades, the slope of trend with the 95% confident interval and its standard deviations (Std) were measured (Figure 9k,l). Figure 9m depicts the probability of having a positive rate of change in trend that assists in identifying the disturbance types (e.g., deforestation and recovery). The Landsat images corresponding to the detected changes are shown in Figure 10 to help understand the changing process. A total of six abrupt changes were detected (Figure 9a,b) and one abrupt change (approximately occur in 2012) was observed by both the seasonal and trend signals.  : (a,b) represent the best fitted seasonal and trend signals (red line), respectively. The 95% credible interval (gray envelope) of corresponding signals, the most possible changepoints (dotted line) as well as the possible locations of the break time (lemon envelope) were also calculated by BEAST. In addition, several parameters, such as the occurrence probability of being changepoint, the standard deviations (Std) of the seasonal/trend components, the estimated harmonic order used to approximate signals, the probability distribution of the total number of changepoints, the slope of trend and its standard deviations were shown in (c-l). (m) represents the estimated probability of having a positive rate of change in trend.
The Landsat images corresponding to the detected changes are shown in Figure 10 to help understand the changing process. A total of six abrupt changes were detected (Figure 9a,b) and one abrupt change (approximately occur in 2012) was observed by both the seasonal and trend signals. Firstly, a sudden change occurred in the seasonality of September 2006, which might be caused by phenological shifts (Figure 9a). In Figure 10, the area has no vegetation on 5 May 2006, but has forest on 27 May 2007. In this duration, the stand recovered and this change has been accurately detected by BEAST. The vegetation was a stable from May 2007 to January 2009, but a sudden drop of NDVI value occurred in February 2010 and caused another abrupt change (clear cutting). The second trend changepoint was found in September 2011, reflecting the NDVI value jump caused by replantation post-deforestation. Thereafter, the forest stand has experienced canopy closure, selective cutting and the second canopy closure. All these abrupt events caused by managements and recovery were correctly detected with a high temporal accuracy.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 24 Firstly, a sudden change occurred in the seasonality of September 2006, which might be caused by phenological shifts (Figure 9a). In Figure 10, the area has no vegetation on 5 May 2006, but has forest on 27 May 2007. In this duration, the stand recovered and this change has been accurately detected by BEAST. The vegetation was a stable from May 2007 to January 2009, but a sudden drop of NDVI value occurred in February 2010 and caused another abrupt change (clear cutting). The second trend changepoint was found in September 2011, reflecting the NDVI value jump caused by replantation post-deforestation. Thereafter, the forest stand has experienced canopy closure, selective cutting and the second canopy closure. All these abrupt events caused by managements and recovery were correctly detected with a high temporal accuracy.

Spatial Validation
In the spatial validation, a total of 2673 pixels of browning disturbance (degradation), 1552 pixels of greening vegetation (recovery) and 4264 pixels of stable vegetation extracted with the NFCS 2018 database and TimeSync-plus tool were used to verify the seasonal and trend changepoints in the research duration. The changepoint detection accuracies of the presented method are shown in Table 1. The promising results of seasonal and trend changepoints were achieved with an OA of 87.37% and Kappa coefficient of 0.85. The highest accuracy was achieved for stable vegetation, with PA = 88.50% and UA = 86.89%. The PA and UA of the three kinds of vegetation are all higher than 85%. The detection accuracies of greening pixels are higher than that of browning pixels, which means that the proposed

Spatial Validation
In the spatial validation, a total of 2673 pixels of browning disturbance (degradation), 1552 pixels of greening vegetation (recovery) and 4264 pixels of stable vegetation extracted with the NFCS 2018 database and TimeSync-plus tool were used to verify the seasonal and trend changepoints in the research duration. The changepoint detection accuracies of the presented method are shown in Table 1. The promising results of seasonal and trend changepoints were achieved with an OA of 87.37% and Kappa coefficient of 0.85. The highest accuracy was achieved for stable vegetation, with PA = 88.50% and UA = 86.89%. The PA and UA of the three kinds of vegetation are all higher than 85%. The detection accuracies of greening pixels are higher than that of browning pixels, which means that the proposed method is more sensitive to vegetation greening than browning events. Additionally, all the PA values are higher than the UA values, indicating more commission errors are in change detection.

Temporal Validation
To evaluate the temporal accuracy in dynamic detection, 76 sampling plots experiencing vegetation disappearance in the research duration were used for validation. We derived the spatial and temporal information of those plots with the TimeSync-plus tool, which can extract the deforestation issues within a maximum of 30 days from all available Landsat and Google Earth images. The temporal interval intersection and mean interval length are shown in Table 2. The interval length is the duration between the sudden rise and drop of NDVI in time series. High temporal accuracy for vegetation degradation events (84.71% interval intersection) and recovery events (88.42% interval intersection) were achieved (Table 2). Additionally, the mean interval length of validation (113.1 days in average) and the detected results (119.9 days in average) were calculated. The detection results highlight the high consistency between validation and detection at time scale, which achieved the interval intersection of 86.57% in average and a small absolute difference (6.8 days) of the mean interval length.

Spatial-Temporal Dynamics of Wetland Vegetation in the Dongting Lake Area
In this study, the total number of seasonal and trend changepoints and the slope of trend were estimated by the BEAST algorithm in the period 2000-2019 over the study area. Figure 11a,b visualize all the seasonal and trend changepoints in Dongting Lake wetland, suggesting that the change in trend is more sharply than that in seasonality. The cropland (e.g., the paddy rice field in the central part of the study area) has less changes in seasonality than in trend, indicating that it has stable intra-annual phenology but frequent crop rotations. The wetland herb (the northeastern part of the study area, surrounding the Dongting Lake) has similar change features with cropland, but the trend has a negative slope (~−0.1/20 a) (Figure 11c). This may be caused by the climate change. Most changepoints in seasonality was detected in the Pinus massoniana-dominated area in the southwestern part of the study area, because that area suffered from pine moth. The changepoints of the permanent water body were extracted precisely, which are the least in both seasonality and trend. This further demonstrates the feasibility of the BEAST method. In Figure 11a,c, the stable area has positive slope of trend, and the area with changes has negative slope of trend. The area with stable vegetation has a greening trend, which is possible due to the global warming. The natural wetland has most changing areas and faces the risk of degradation, because of the frequent human activities and environment changes. We calculated the proportion of pixels with single/multiple changepoint (s) (Figure 12a) and the proportion of stable and changed pixels in different DEM ranges (Figure 12b). In general, the seasonal changepoints are less than the trend changepoints, which means that permanent changes (e.g., deforestation) are more than ephemeral changes (e.g., flooding). Most areas have two or less seasonal changepoints but at least three trend changepoints. In addition, 80% pixels are changed (Figure 12b), indicating the dramatic vegetation changes over the whole Dongting Lake wetland during 2000-2019. Besides, the changed pixels mainly distribute in the area lower than 20 m (22% of the total pixels) and between 30~40 m (18% of the total pixels), which are the core wetland and built-up area. Forested wetland (southwestern of the study area) higher than 40 m, especially the region higher than 50 m, is the relatively stable vegetation area. The regions with poplar and wetland herb surrounding the water body have negative slopes in trend during the research period (Figure 11c), and they have the most frequent changes. We calculated the proportion of pixels with single/multiple changepoint (s) (Figure 12a) and the proportion of stable and changed pixels in different DEM ranges (Figure 12b). In general, the seasonal changepoints are less than the trend changepoints, which means that permanent changes (e.g., deforestation) are more than ephemeral changes (e.g., flooding). Most areas have two or less seasonal changepoints but at least three trend changepoints. In addition, 80% pixels are changed (Figure 12b), indicating the dramatic vegetation changes over the whole Dongting Lake wetland during 2000-2019. Besides, the changed pixels mainly distribute in the area lower than 20 m (22% of the total pixels) and between 30~40 m (18% of the total pixels), which are the core wetland and built-up area. Forested wetland (southwestern of the study area) higher than 40 m, especially the region higher than 50 m, is the relatively stable vegetation area. The regions with poplar and wetland herb surrounding the water body have negative slopes in trend during the research period (Figure 11c), and they have the most frequent changes.
To better understand the evolution mechanism of wetland vegetation, we extracted the time with the most significant disturbances in seasonality and trend (Figure 13a,b), and measured the year with the maximum NDVI ( Figure 13c). As Figure 13a shows, many significant changes intensively occurred in 2008, because the most serious frozen disaster in more than a hundred years struck southern China and led to the phenological shift and vegetation death. The land cover that was transformed from vegetation to build-up area in dispersed study regions in the early years of the research duration has been detected by BEAST (Figure 13b). In addition, the most significant disturbance of trend was observed in 2011, when the clearing of poplar policy was carried out. The extensive deforestation of poplar also brought vegetation death and degeneration of ecosystem. Figure 13c shows that the greenhouse effect may benefit the growth of vegetation and crops, because the NDVI peak of the stable area occurred in the years after large areas of poplars have been removed since 2011. The maximum NDVI value was captured in 2011 in some areas, because the extensive deforestation of poplar led to the gradual decrease in NDVI since 2011. All of the results demonstrate that the vegetation changes vary spatially and temporally.
Besides, the changed pixels mainly distribute in the area lower than 20 m (22% of the total pixels) and between 30~40 m (18% of the total pixels), which are the core wetland and built-up area. Forested wetland (southwestern of the study area) higher than 40 m, especially the region higher than 50 m, is the relatively stable vegetation area. The regions with poplar and wetland herb surrounding the water body have negative slopes in trend during the research period (Figure 11c), and they have the most frequent changes. To better understand the evolution mechanism of wetland vegetation, we extracted the time with the most significant disturbances in seasonality and trend (Figure 13a,b), and measured the year with the maximum NDVI ( Figure 13c). As Figure 13a shows, many significant changes intensively occurred in 2008, because the most serious frozen disaster in more than a hundred years struck southern China and led to the phenological shift and vegetation death. The land cover that was transformed from vegetation to build-up area in dispersed study regions in the early years of the research duration has been detected by BEAST (Figure 13b). In addition, the most significant disturbance of trend was observed in 2011, when the clearing of poplar policy was carried out. The extensive deforestation of poplar also brought vegetation death and degeneration of ecosystem. Figure 13c shows that the greenhouse effect may benefit the growth of vegetation and crops, because the NDVI peak of the stable area occurred in the years after large areas of poplars have been removed since 2011. The maximum NDVI value was captured in 2011 in some areas, because the extensive deforestation of poplar led to the gradual decrease in NDVI since 2011. All of the results demonstrate that the vegetation changes vary spatially and temporally. The proportion of pixels with abrupt changes and that with the maximum NDVI value are shown in Figure 14a. The annual average temperature and precipitation over the study area in the research duration were calculated (Figure 14b) to investigate the relation between vegetation variation and climate change. As Figure 14 shows, the climate changes are not the main reason for vegetation variation, though the annual average temperature and precipitation in recent years are slightly different from that of previous years. The seasonal and trend change both peak in 2008 (Figure 14a), when the frozen disaster happened and resulted in a lots of vegetation death [59]. The change of Dongting Lake wetland in the later period of the study duration is larger than that in the early period, and the change accelerated gradually, because of the environmental pollutions and rapid urbanization in recent years. The maximum NDVI value appears in 2011, because the poplar cleaning policy in the Dongting Lake area came into force. Most poplars were cut in 2011, and the effect of the policy lasted for several years (five years since 2011). In general, extreme weathers and policy incentives are the main reasons for the vegetation change in the last 20 years. The proportion of pixels with abrupt changes and that with the maximum NDVI value are shown in Figure 14a. The annual average temperature and precipitation over the study area in the research duration were calculated (Figure 14b) to investigate the relation between vegetation variation and climate change. As Figure 14 shows, the climate changes are not the main reason for vegetation variation, though the annual average temperature and precipitation in recent years are slightly different from that of previous years. The seasonal and trend change both peak in 2008 (Figure 14a), when the frozen disaster happened and resulted in a lots of vegetation death [59]. The change of Dongting Lake wetland in the later period of the study duration is larger than that in the early period, and the change accelerated gradually, because of the environmental pollutions and rapid urbanization in recent years. The maximum NDVI value appears in 2011, because the poplar cleaning policy in the Dongting Lake area came into force. Most poplars were cut in 2011, and the effect of the policy lasted for several years

Discussion
The accurate change detection of the ecosystem and vegetation remains a challenge, particularly in sub-tropical and tropical regions. A few of researches attempted to address this problem ( Table 3). The BFAST algorithm was utilized to detect the forest loss in the Aberdare National Park, Kenya, using the 500-m resolution MODIS EVI data [60]. However, it only detected 72% of forest loss accurately, due to a large amount of mixing and missing pixels in the MODIS pixels. Furthermore, it cannot perform the temporal detection. Another study focused on the accurate detection of forest change in Sierra Madre Occidental mountain and achieved a high detection accuracy (97.6%) using the VCT method [61]. The discontinuous Landsat data used in [61] result in the inability of obtaining high accuracy in time scale. Thus, the change region was only analyzed annually in the research. The BFAST algorithm was also used to detect the forest and grassland change over the homogeneous area in Queensland, Australia, taking the STARFM fused time series as input [9]. Despite the 96.2% changepoint estimation accuracy of the clearing event and 83.6% of the re-clearing event were achieved, the timing accuracy was seriously delayed (40 days). For improving image fusion and increasing the detection precision of change types, the ESTARFM model and random forest classification algorithm were adopted [42]. But the uncertainty in the detection procedure was not estimated and the slow disturbance event (such as diseases and pests) could not be detected. In general, the methods in [9,42] are not suitable for monitoring the vegetation dynamics in the Dongting Lake area as the high heterogeneity in this region where the STARFM and ESTARFM algorithms could not provide accurate fused images.

Discussion
The accurate change detection of the ecosystem and vegetation remains a challenge, particularly in sub-tropical and tropical regions. A few of researches attempted to address this problem ( Table 3). The BFAST algorithm was utilized to detect the forest loss in the Aberdare National Park, Kenya, using the 500-m resolution MODIS EVI data [60]. However, it only detected 72% of forest loss accurately, due to a large amount of mixing and missing pixels in the MODIS pixels. Furthermore, it cannot perform the temporal detection. Another study focused on the accurate detection of forest change in Sierra Madre Occidental mountain and achieved a high detection accuracy (97.6%) using the VCT method [61]. The discontinuous Landsat data used in [61] result in the inability of obtaining high accuracy in time scale. Thus, the change region was only analyzed annually in the research. The BFAST algorithm was also used to detect the forest and grassland change over the homogeneous area in Queensland, Australia, taking the STARFM fused time series as input [9]. Despite the 96.2% changepoint estimation accuracy of the clearing event and 83.6% of the re-clearing event were achieved, the timing accuracy was seriously delayed (40 days). For improving image fusion and increasing the detection precision of change types, the ESTARFM model and random forest classification algorithm were adopted [42]. But the uncertainty in the detection procedure was not estimated and the slow disturbance event (such as diseases and pests) could not be detected. In general, the methods in [9,42] are not suitable for monitoring the vegetation dynamics in the Dongting Lake area as the high heterogeneity in this region where the STARFM and ESTARFM algorithms could not provide accurate fused images. We used the cross-fusion strategy and FSDAF model to generate the dense Landsat time series stack, which made the detection more sensitive and robust. For example, the Pinus massoniana forest suffered pine moth was accurately identified. Compared with the BFAST algorithm, the BEAST can estimate more extra parameters, including probability distribution of timing, number of changepoints and the rate of change in trend signal, so it can achieve satisfying performance in both spatial and temporal domains. The spatial and temporal validation shows high OA of 87.37% and Kappa coefficient of 0.85%, and a satisfactory timing accuracy with an interval intersection of 86.57% and absolute difference of mean interval length of 6.8 days. However, there are commission and omission errors in the detection process over the time series. The commission errors are possibly due to the noise from the Savitzky-Golay filter in MODIS processing and the prediction errors in the FSDAF fusion. The Savitzky-Golay filter was used to smooth the NDVI time series by reducing noise and cloud-contaminated pixels, but it influenced the abrupt changes in time series and induced commission errors [9]. In this study, the cross-fusion and FSDAF algorithm were used to blend MODIS and Landsat observations to generate a continuous monthly NDVI time series. However, the fusion results have uncertainties, despite the high correlation between the fused and actual images. On one hand, there are some prediction errors. Indistinguishable changes cannot be reflected in the base coarse-resolution MODIS images [8]. On the other hand, the residual cloud pixels in a monthly MODIS NDVI are ubiquitous, particularly in subtropical and tropical areas, which may also result in fusion errors that are mistaken as changes. Additionally, the main reason for omission errors are the slight and continuous transformation of vegetation to other land cover types or the inter-conversion between vegetation, which induce indistinguishable spectral variation that cannot be detected by BEAST.
We focused on the detection of abrupt changes in the vegetation area caused by disturbance, vegetation recovery and regrowth, and also the phenology and species succession as well as the diseases and insect pests. The phenological shift caused by species succession, plant diseases and insect pests would only show as slow changes that cannot be immediately detected by the proposed framework. The slow changes caused by other issues (e.g., urban expansion and forest cutting) would also delay the change detection. Besides, the quality of the reconstructed NDVI time series is a key factor to influence the detection accuracy at temporal scale. A 250-m pixel with small changes that might not be observed in the base MODIS image may be taken as an unchanged pixel to extend the detection interval and lead to an inconsistent intersection. The maximum separation distance in trend was set as 12 months (referred to one year) to reduce pseudo changes and this setting is suitable for most pixels. That means the setting only achieved local optimum but not global optimum, which would lead to changepoints that were too close. Thus, the delay of detection timing is also attributed to this situation, because the real cycle of some disturbances is shorter than 12 months.
In summary, the proposed method preforms well in extracting the vegetation changes from monthly time series, as it benefits from the synthetic use of FSDAF and the BEAST algorithm. However, there are several limitations in practical applications, which might be solved by the following strategies.
Data: This study only used NDVI to characterize the vegetation growth for detecting the disturbance. Other vegetation indices (e.g., red-edge band vegetation index) may also be useful and flexible for vegetation monitoring. Besides, SAR data are immune to the cloudy and rainy weather and sensitive to the structure, moisture and dielectric properties of ground targets, so can be employed for continuous monitoring of vegetation. The red-edge and shortwave infrared reflectance is critical for the early detection of vegetation diseases and insect pests [62]; the temporal rate of change of the red-edge chlorophyll index and the structure-related NDVI can be used to detect forest loss [63]; multivariate SAR images have been employed to detect vegetation changes and have got promising results [64]. Thus, exploring other time-series resources that are more sensitive and robust to vegetation variation is one of the key points in future work.
Algorithm: The BEAST algorithm has solved the following three issues of the conventional decomposition method. Ecologically, the definition of the change components (abrupt changes, seasonality and trend) in BEAST is more inclusive, which helps to detect the indiscernible abnormality in time series and investigate the potential driver of phenological shift. Mathematically, the BEAST algorithm uses BMA and MCMC strategies to fit the seasonal and trend signals and quantify the uncertainty of fitting and detection, but not to seek a single beast model that may result in poor fitting.
Computationally, more parameters about vegetation variation can be derived by the BEAST algorithm, such as the probability of being a changepoint, probability distribution of total changepoint number, the estimated harmonic order and the probability of having a positive rate of change in trend. However, the BEAST algorithm has larger computation loads than other algorithms. When it is used to monitor vegetation dynamics over a large scale using high-resolution data, the huge computation load would limit the application. The Google Earth Engine (GEE) platform combines a range of satellite imagery and geospatial datasets with cloud-based computation power for environmental data analysis, which is conducive to implement the proposed method. Besides, BEAST cannot consider both spatial and temporal information simultaneously, and it only uses the information of temporal variation to seek abrupt changes in a pixel. Therefore, the detection accuracy depends on the quality and interval of the time-series images in the temporal domain. Actually, there are spatial correlations between the vegetation change and natural or human factors. Furthermore, the time series decomposition algorithm based on individual pixel cannot characterize a series of vegetation changes in the spatial domain. Developing an object-based time series decomposition algorithm may link the spatial and temporal information of vegetation changes and solve the problem.

Conclusions
The study presents a strategy integrating FSDAF and BEAST algorithms to detect and characterize the long-term vegetation changes in wetland areas by using the fused time series. The strategy considers both similarity and consistency in selecting base image pairs for FSDAF fusion. A high correlation was found between the FSDAF-fused and the actual NDVIs, with R 2 = 0.93 and RMSE = 0.051. The fused NDVI time series is suitable for monitoring dynamic disturbance in heterogeneous areas. Besides, the BEAST algorithm robustly detected abrupt changes, seasonality and trend as well as uncertainties of the vegetation in the study area in the past two decades. In addition to the indiscernible deforestation and degradation, the plant diseases and insect pests and climate-driven gradual vegetation change events were reflected by the seasonality and trend. Both the spatial and temporal accuracy of the detection were evaluated, with 87.37% OA and 86.57% interval intersection. In the research period, the wetland vegetation in the Dongting Lake area changed dramatically in distribution and durations, because of the combined impacts of climate changes and human activities. Results indicate that the frozen disaster in early 2008 over southern China is the main factor leading to phenological shift and vegetation death in Dongting Lake wetland, and the ecological conservation policies instituted by the governments (i.e., poplar cleaning and replacement) are the key points to explain the abrupt vegetation changes in trend. Generally, the extreme weather, (frozen disaster) and anthropogenic factors (government policies, variation of economic profits, land reclamation for developing agriculture and industry as well as rapid urbanization) are the primary reasons for vegetation changes. Additionally, the stable area (e.g., cropland and relatively high-altitude forest) had a greening trend during the last 20 years, which is an evident of the global warming bringing vegetation greening. The proposed method fills the gap of continuous vegetation monitoring in heterogeneous wetland areas using long-term data with high spatial and temporal resolutions. This method can reveal the potential drivers for the degradation and recovery of wetland vegetation in subtropical areas, which can be extended to other large-scale vegetation monitoring.