Multi-Type Forest Change Detection Using BFAST and Monthly Landsat Time Series for Monitoring Spatiotemporal Dynamics of Forests in Subtropical Wetland

Land cover changes, especially excessive economic forest plantations, have significantly threatened the ecological security of West Dongting Lake wetland in China. This work aimed to investigate the spatiotemporal dynamics of forests in the West Dongting Lake region from 2000 to 2018 using a reconstructed monthly Landsat NDVI time series. The multi-type forest changes, including conversion from forest to another land cover category, conversion from another land cover category to forest, and conversion from forest to forest (such as flooding and replantation post-deforestation), and land cover categories before and after change were effectively detected by integrating Breaks For Additive Seasonal and Trend (BFAST) and random forest algorithms with the monthly NDVI time series, with an overall accuracy of 87.8%. On the basis of focusing on all the forest regions extracted through creating a forest mask for each image in time series and merging these to produce an ‘anytime’ forest mask, the spatiotemporal dynamics of forest were analyzed on the basis of the acquired information of multi-type forest changes and classification. The forests are principally distributed in the core zone of West Donting Lake surrounding the water body and the southwestern mountains. The forest changes in the core zone and low elevation region are prevalent and frequent. The variation of forest areas in West Dongting Lake experienced three steps: rapid expansion of forest plantation from 2000 to 2005, relatively steady from 2006 to 2011, and continuous decline since 2011, mainly caused by anthropogenic factors, such as government policies and economic profits. This study demonstrated the applicability of the integrated BFAST method to detect multi-type forest changes by using dense Landsat time series in the subtropical wetland ecosystem with low data availability.


Introduction
Dramatic changes have occurred in many wetland ecosystems in the previous decades and have been attributed to the influence of climatic changes and human activities [1].Research demonstrated that wetland conversion and loss on the global scale has exceeded 50% since 1990 [2].As the second largest freshwater lake in China, Dongting Lake is also one of the important wetland ecosystems.However, the land cover changes, such as land reclamation for developing agriculture and industry, and rapid urbanization have significantly threatened the ecological security of Dongting Lake wetland [3].With the drop in grain prices and profitable papermaking, a variety of land cover categories, such as paddy fields and bottomland, especially water body, have been transformed to plant the economic forests (i.e., Populus euramericana) for papermaking since the late 1990s.A large number of Populus plantations have accelerated the desiccation of wetland [4], resulting in a decrease of the wetland area [5][6][7].This situation leads to the death of wetland vegetation [8]; induces biodiversity loss [7]; and influences the community composition [4], structure [7] and function [4] of the wetland ecosystem.Therefore, accurately detecting the land cover changes related to forests and then monitoring spatiotemporal dynamics of forests in the Dongting Lake region are critical for protecting wetland ecosystems.
With the accessibility and accumulation of remotely sensed images on the spatiotemporal scale, remote sensing is a valid alternative for monitoring ecosystem changes.For monitoring spatiotemporal dynamics of forests in the Dongting Lake region with remote sensing, there are the following questions need to be answered: (1) how to accurately detect land cover change?(2) how to acquire the process information of ecosystem change, including land cover change, not just the final results?(3) how to distinguish land cover change from other change types?and (4) how to get the detailed "from what, to what" information on the premise that the land cover change has been identified?
The early common methods, which are quite simple and straightforward, include image differencing, principal component analysis, tasseled cap transform, post-classification comparison, and change vector analysis [9][10][11].The ideas of these methods are that the changes are generally detected through comparing the difference between the images in two different time phases.Nevertheless, the accuracy of change detection using these methods is dependent on the set threshold [11,12], the availability of ideal images in the same season or the classification accuracy of two images.Another disadvantage of these methods is that the changes are generally detected through comparing the difference between the images in two different time phases without considering the context information in the time dimension, lacking change law and process research [13,14].Therefore, these methods based on two-temporal images cannot answer questions 1 and 2.
Thus, some change-detection methods based on change trajectories of time series, including Vegetation Change Tracker and Landtrendr, have been developed for studying the change law and process [15,16] and have been widely used for forest disturbances and regrowth monitoring and land cover change detection [17][18][19].For example, Griffiths et al. employed the LandTrendr and an annual Landsat time series to derive annual forest disturbance maps along with recovery dynamics for answering the question of how the drastic institutional and socio-economic transformation after the collapse of socialism in Romania affected forestry [17].Kennedy et al. utilized the LandTrendr to detect change and then employed the random forest (RF) classifier (RFC) to label agents of change including urbanization, forest management, and natural change for monitoring habitat in the Puget Sound region, USA [18].The frequency of time series within these methods is typically one image per year in the same season to minimize seasonal variation and sun angle differences.However, these methods developed at yearly resolution have failed in accurately capturing the timing of the change at the sub-annual scale.This situation results in the detected timing delay of land cover change.Therefore, the optimal result that these methods could provide is generally annual or biennial change [20,21].Similarly, these methods based on multi-temporal images also cannot answer question 2.
Only the time series with high frequency (e.g., monthly) could describe the entire process of land cover change over shorter time intervals because of the abrupt characteristic of land cover change.Therefore, only change detection with dense satellite time series indeed satisfies the requirement of dynamic monitoring of land cover change [22].One of the major challenges of change detection with dense time series is to distinguish land cover change from other change types on the premise of minimizing the influences of seasonal phenological variation.The Breaks For Additive Seasonal and Trend (BFAST) meet the challenge through the iterative decomposotion of the time series into trend, seasonal and remainder components [12,23].The breaks detected in the trend component represent an ephemeral disturbance [24,25], while the breaks detected in the seasonal component represent a land cover change because of the seasonal pattern (phenology) difference between different land cover categories [25,26].BFAST has been applied in many different land cover categories, such as wetland [24], wildlife nature reserve [25], city [26], tropical dry forest [27], vegetation [28,29], agriculture [30], and abandoned energy development sites [31].Although BFAST can identify land cover changes (question 1), acquire the process information of land cover change (question 2), and distinguish land cover change from other change types (question 3), it still cannot provide the detailed "from what, to what" information (question 4).Therefore, the land cover categories should be classified before and after abrupt changes while finding the timing and location of land cover change [14,32].Zhu and Woodcock proposed a continuous change detection and classification (CCDC) algorithm for land cover change detection and land cover categories classification after change detection by using all available Landsat images [21].The advantage of CCDC is that the coefficients of the fitting time series model and the root mean square error (RMSE) of it are treated as the input parameters of the RFC for land cover classification before and after abrupt changes (question 4).
The Dongting Lake region, in which farmland is one of the major land use categories, is one of the important grain production bases in China.The inter-conversion between farmland and other land cover categories, including forest, frequently occurred because of the policy orientation and economic profit motivation in recent decades.Therefore, the phenological trajectory of all kinds of land cover categories, including farmland with complex intra-annual variations (such as double-or even triple-crop patterns within a year), must be well described for land cover change detection and classification before and after abrupt change.However, the CCDC just adopted a simple sinusoidal model for the phenological trajectory description for a variety of land cover categories, which means it cannot answer any one of four questions in Dongting Lake region with complex change processes [21,33].Although the BFAST algorithm exhibits difficulty in interpreting the land cover categories before and after the detected abrupt changes, a multi-harmonic model adopted in BFAST could be satisfactory for the complex shape of phenological trajectory and could be applied in the Dongting Lake region.Therefore, BFAST may be feasible for not only detecting the land cover changes related to forest but also classifying the land cover categories before and after abrupt changes by integrating BFAST and classification algorithm.
BFAST is initially proposed to detect change based on high temporal resolution MODIS data, which means it needs dense time series images.However, MODIS data with coarse spatial resolution have the limitation of detecting small changes because of the high fragmentation degree of land cover in the West Dongting Lake region [34].Lots of medium-resolution data could be freely captured with the Landsat data archive's opening in 2008.This event enabled change detection with dense remotely-sensed time series at finer spatial resolution [35].In comparison with the MODIS time series, the Landsat time series possess several advantages, such as higher spatial resolution (30 m) and longer archive (exceeding 47 years), which means more accurate change detection and longer terms of monitoring [36][37][38].However, the relatively low temporal resolution and frequent cloud cover and shadows result in excessive temporal gaps in the Landsat archive for certain regions or time periods [39], such as the West Dongting Lake region in South China (subtropical region with low optical data availability) in summer.Therefore, for investigating the spatiotemporal dynamics of forest in the West Dongting Lake region with dense Landsat time series by using BFAST, we should answer the new question in this study: (5) how does one acquire the dense Landsat time series with high quality for change detection using BFAST?
This research aims to investigate the spatiotemporal dynamics of forest in the West Dongting Lake region with dense Landsat time series using BFAST through answering five questions.In this research, the timing and location of forest changes, changing types (from forest to what, or from what to forest), and process are acquired by integrating BFAST and classification algorithm.Then the spatiotemporal variations of forest and the potential drivers in the West Dongting Lake region are analyzed on the basis of the information of change detection and classification.

Study Area
The Donting Lake region, the second largest freshwater lake and one of the major grain production areas in China, is located in northern Hunan Province (include three prefectural-level divisions of Changde, Yueyang and Yiyang) along the southern bank of the Yangtze River (Figure 1a,b).The major land cover categories in this study area include cropland, water body, forest, shoal (including mudflats, marshland and reed land), and urban (including settlements and roads) in this study area.Land cover changes frequently occurred in this region in recent decades due to natural disasters and human activities.As one of the major land cover categories in this area, the land cover changes related to forest, especially economic forests (i.e., Populus euramericana), are significant.It takes only 3-4 years to make the Populus euramericana with fast growth and strong adaptability grow into useful timber for papermaking.The mass planting of Populus euramericana began in the late 1990s in West Donting Lake (Figure 1c).Many primary and secondary forests, such as Metasequoia and willow, have been replaced by the Populus euramericana.A variety of land cover categories, including farmland and wetland, have been transformed to plant the Populus euramericana due to decline in grain prices and the increasing profitability of paper-making.Parts of the planting areas have disappeared because of flood, urbanization, decline in timber prices, government policy, and construction of the Three Gorges Dam (TGD).As the study area in this work, West Dongting Lake is artificially divided into two parts, namely, the core zone (yellow border in Figure 1d,e) and the transition zone, according to the distance from the main water body and demand for wetland protection (Figure 1d).The core zone includes West Dongting Lake Nature Reserve and buffer zone of the trunk streams flowing into West Dongting Lake.The transition zone, which has some distance from the main water body, is in the adjacent area where some small lakes and tributaries are dispersedly distributed.The value range of Global Digital Elevation Model (GDEM) is mostly from 20 to 40 m in West Dongting Lake, except for the southwestern mountains with GDEM greater than 50 m (Figure 1e).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 33 spatiotemporal variations of forest and the potential drivers in the West Dongting Lake region are analyzed on the basis of the information of change detection and classification.

Study Area
The Donting Lake region, the second largest freshwater lake and one of the major grain production areas in China, is located in northern Hunan Province (include three prefectural-level divisions of Changde, Yueyang and Yiyang) along the southern bank of the Yangtze River (Figure 1a,b).The major land cover categories in this study area include cropland, water body, forest, shoal (including mudflats, marshland and reed land), and urban (including settlements and roads) in this study area.Land cover changes frequently occurred in this region in recent decades due to natural disasters and human activities.As one of the major land cover categories in this area, the land cover changes related to forest, especially economic forests (i.e., Populus euramericana), are significant.It takes only 3-4 years to make the Populus euramericana with fast growth and strong adaptability grow into useful timber for papermaking.The mass planting of Populus euramericana began in the late 1990s in West Donting Lake (Figure 1c).Many primary and secondary forests, such as Metasequoia and willow, have been replaced by the Populus euramericana.A variety of land cover categories, including farmland and wetland, have been transformed to plant the Populus euramericana due to decline in grain prices and the increasing profitability of paper-making.Parts of the planting areas have disappeared because of flood, urbanization, decline in timber prices, government policy, and construction of the Three Gorges Dam (TGD).As the study area in this work, West Dongting Lake is artificially divided into two parts, namely, the core zone (yellow border in Figure 1d,e) and the transition zone, according to the distance from the main water body and demand for wetland protection (Figure 1d).The core zone includes West Dongting Lake Nature Reserve and buffer zone of the trunk streams flowing into West Dongting Lake.The transition zone, which has some distance from the main water body, is in the adjacent area where some small lakes and tributaries are dispersedly distributed.The value range of Global Digital Elevation Model (GDEM) is mostly from 20 to 40 m in West Dongting Lake, except for the southwestern mountains with GDEM greater than 50 m (Figure 1e).

Methods
Figure 2 shows an overview of the general workflow followed in this study.The monthly Landsat NDVI time series were generated by using two spatiotemporal interpolation methods named Neighborhood Similar Pixel Interpolator (NSPI) and its modified version MNSPI, and a spatiotemporal data fusion algorithm named enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM).The BFAST was then adopted to detect the abrupt changes with the monthly NDVI time series.The land cover categories before and after abrupt changes were identified through setting the coefficients of the time series model as the input parameters of the RF algorithm.To determine forest extent, we created a forest mask for each image in the time series and merged these masks to produce an 'anytime' forest mask for extracting all the forest pixels.After accuracy assessment of change detection and classification, the spatiotemporal distributions of the number and timing of abrupt changes related to forest were analyzed.Lastly, we focused on describing the variation of forest area in West Dongting Lake and discussing the potential drivers.The critical procedures will be introduced in the following sections in detail.

Methods
Figure 2 shows an overview of the general workflow followed in this study.The monthly Landsat NDVI time series were generated by using two spatiotemporal interpolation methods named Neighborhood Similar Pixel Interpolator (NSPI) and its modified version MNSPI, and a spatiotemporal data fusion algorithm named enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM).The BFAST was then adopted to detect the abrupt changes with the monthly NDVI time series.The land cover categories before and after abrupt changes were identified through setting the coefficients of the time series model as the input parameters of the RF algorithm.To determine forest extent, we created a forest mask for each image in the time series and merged these masks to produce an 'anytime' forest mask for extracting all the forest pixels.After accuracy assessment of change detection and classification, the spatiotemporal distributions of the number and timing of abrupt changes related to forest were analyzed.Lastly, we focused on describing the variation of forest area in West Dongting Lake and discussing the potential drivers.The critical procedures will be introduced in the following sections in detail.

Datasets and Preprocessing
The remotely sensed time series with coarse spatial resolution, such as MODIS and NOAA/AVHRR, are often used for land cover change detection because of their high temporal resolution.However, they have the limitation of detecting small changes because of the high fragmentation degree of land cover in the West Dongting Lake region.Considering the temporal resolution (16 days) of Landsat images and time span of abrupt change and its influence, the time interval between different images in time series is set as 1 month for accurate change detection.Therefore, a monthly Landsat time series with 19 years, from 2000 to 2018, have been adopted for change detection and analyzing spatiotemporal variations of forest in West Dongting Lake region.
Although the Landsat data records were continuously captured from 1972 to present, the relatively low temporal resolution and frequent cloud cover and shadows result in excessive temporal gaps in the Landsat archive for certain regions or time periods [39], such as the West Dongting Lake region in South China (subtropical region with low optical data availability) in summer.Moreover, the permanent failure of the scan-line corrector (SLC) of the ETM+ on Landsat 7 caused approximately 22% of the pixels to be unscanned in ETM+ image since 2003 [40].In spite of three sensors with 16-day revisit cycle operation, only 197 Landsat images, including 56 TM, 103 ETM+ and 38 OLI, with less than 60% of pixels contaminated by clouds, cloud shadows, and SLC-off gaps were collected in the West Dongting Lake region between 2000 and 2018 (Table 1, Figure 3).The region was not covered by the Landsat images with contaminated pixels less than 60% within 59 months (Table 1).Two spatiotemporal interpolation methods named NSPI for SLC-off gaps filling and its modified version MNSPI for cloud removal were integrated into an iterative process for interpolating gaps pixels and cloudy pixels, respectively [40][41][42].This method is effective for the time series without extremely persistent clouds, such as the 197 Landsat images with less than 60% of contaminated pixels.However, it is difficult to fill the contaminated pixels with acceptable accuracy when the clouds cover most (over 60%) of the image.Therefore, the missing images in 59 time phases need to be reconstructed in another way to generate a monthly Landsat time series between 2000 and 2018, a total of 228 images (12 scenes per year multiply by 19 years).A spatiotemporal data fusion algorithm named ESTARFM was used to reconstruct the missing Landsat images by fusing available Landsat and MODIS images [43].Seventy-three MODIS surface reflectance 8-day composite images for data fusion were obtained (Table 2).The ESTARFM uses two pairs of cloud-free fine-resolution Landsat and coarse-resolution MODIS images at base date tm and tn, and one cloud-free coarse-resolution MODIS image at prediction date tp to construct the fine-resolution image at tp.The ESTARFM (1) uses a

Datasets and Preprocessing
The remotely sensed time series with coarse spatial resolution, such as MODIS and NOAA/AVHRR, are often used for land cover change detection because of their high temporal resolution.However, they have the limitation of detecting small changes because of the high fragmentation degree of land cover in the West Dongting Lake region.Considering the temporal resolution (16 days) of Landsat images and time span of abrupt change and its influence, the time interval between different images in time series is set as 1 month for accurate change detection.Therefore, a monthly Landsat time series with 19 years, from 2000 to 2018, have been adopted for change detection and analyzing spatiotemporal variations of forest in West Dongting Lake region.
Although the Landsat data records were continuously captured from 1972 to present, the relatively low temporal resolution and frequent cloud cover and shadows result in excessive temporal gaps in the Landsat archive for certain regions or time periods [39], such as the West Dongting Lake region in South China (subtropical region with low optical data availability) in summer.Moreover, the permanent failure of the scan-line corrector (SLC) of the ETM+ on Landsat 7 caused approximately 22% of the pixels to be unscanned in ETM+ image since 2003 [40].In spite of three sensors with 16-day revisit cycle operation, only 197 Landsat images, including 56 TM, 103 ETM+ and 38 OLI, with less than 60% of pixels contaminated by clouds, cloud shadows, and SLC-off gaps were collected in the West Dongting Lake region between 2000 and 2018 (Table 1, Figure 3).The region was not covered by the Landsat images with contaminated pixels less than 60% within 59 months (Table 1).Two spatiotemporal interpolation methods named NSPI for SLC-off gaps filling and its modified version MNSPI for cloud removal were integrated into an iterative process for interpolating gaps pixels and cloudy pixels, respectively [40][41][42].This method is effective for the time series without extremely persistent clouds, such as the 197 Landsat images with less than 60% of contaminated pixels.However, it is difficult to fill the contaminated pixels with acceptable accuracy when the clouds cover most (over 60%) of the image.Therefore, the missing images in 59 time phases need to be reconstructed in another way to generate a monthly Landsat time series between 2000 and 2018, a total of 228 images (12 scenes per year multiply by 19 years).A spatiotemporal data fusion algorithm named ESTARFM was used to reconstruct the missing Landsat images by fusing available Landsat and MODIS images [43].Seventy-three MODIS surface reflectance 8-day composite images for data fusion were obtained (Table 2).The ESTARFM uses two pairs of cloud-free fine-resolution Landsat and coarse-resolution MODIS images at base date t m and t n , and one cloud-free coarse-resolution MODIS image at prediction date t p to construct the fine-resolution image at t p .The ESTARFM (1) uses a linear mixture model to establish the relationship between the reflectance of a mixed coarse-resolution MODIS pixel and the fine-resolution Landsat pixels within it at base date t m or t n , and then (2) employs a conversion coefficient to convert the reflectance changes of a mixed coarse resolution pixel to the fine-resolution pixels within it based on spatial and spectral similarity.The missing Landsat pixels at t p are reconstructed based on the two relationships.Either the pair of images at t m or t n can be used as the reflectance at base date to compute the fine resolution reflectance at t p , and a more accurate reflectance at t p can be obtained by a weighted combination of the two prediction results.The validity of the ESTARFM was verified by our previous experiment [44].After the data reconstruction, NDVI was calculated from Landsat surface reflectance data, considering the NDVI robustness.Finally, an available monthly NDVI time series per pixel with 30-m resolution during this time period for the responding variable in the following change detection process was established.

Breaks for Additive Seasonal and Trend
BFAST is an iterative algorithm, decomposing time series into three components: trend, seasonal, and remainder components, for change detection [12,23].The simple form of decomposition model (Equation ( 1)) is as follows:   [46].

High-resolution images
Images for validation of change detection and land cover classification, and determination of land cover category of each reference pixel with Landsat images together.
Google earth GDEM ASTER GDEM version 2.0 data with 1 arc-second grid used for analyzing spatiotemporal dynamics of forest.
Remote Sens. 2020, 12, 341 9 of 33 In addition, the high-resolution images were acquired from Google earth for validation of change detection and land cover classification, and determination of land cover category of each reference pixel with Landsat images together (Table 2).The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) GDEM data were collected for analyzing spatiotemporal dynamics of forest (Table 2).

Breaks for Additive Seasonal and Trend
BFAST is an iterative algorithm, decomposing time series into three components: trend, seasonal, and remainder components, for change detection [12,23].The simple form of decomposition model (Equation ( 1)) is as follows: where Y t is the observed data at time t, T t , S t , and e t , which separately represent the trend, seasonal, and remainder component.p represents the number of observations.T t can be fitted as the piecewise linear models (Equation ( 2)) with specific slopes β i and intercepts α i on different m + 1 segments: where i = 1, . . ., m, and m represent the number of breakpoints in the trend component.The seasonal component S t can be fitted as the piecewise harmonic models (Equation ( 3)) on different n + 1 segments: where j = 1, . . ., n, and n is the number of breakpoints in the seasonal component.Variable k is the order of harmonic function and its set as 3 to more accurately detect complex phenological changes, such as double-or even triple-crop patterns within a year in this study area.As a result, not only the land cover categories with simple phenology patterns, such as forest with yearly phenological cycle, but also those with the complex phenology patterns could be accurately described and then classified.The frequency f is set as 12 for annual observations for a monthly Landsat time series in this research.
The iteration process is initialized by estimating the seasonal component Ŝt through utilizing the seasonal-trend decomposition (STL) algorithm [48].First, the residual-based moving sum (MOSUM) test [49] is used to test whether breakpoints τ* 1, . . ., τ* m occur in the adjusted trend component Y t − Ŝt .The position and number of breakpoints in the trend component are then determined through minimizing the Bayesian Information Criterion (BIC) [50] and the residual sum of squares (RSS), respectively.Then robust regression on the basis of M-estimation [51] is used to estimate coefficients α and β in the trend component for each segment, and the trend estimate Tt is computed on the basis of Equation ( 2).Similarly, the MOSUM test is utilized to test whether breakpoints τ# 1, . . ., τ# n occur in the de-trended seasonal component Y t − Tt , and the minimization of the BIC and RSS is conducted to determine the position and number of the breakpoints in the seasonal component.Coefficients γ j,k and θ j,k of the seasonal component in Equation ( 3) for each segment are also estimated on the basis of M-estimation, and then the seasonal estimate Ŝt is computed on the basis of Equation (3).The above steps are iteratively performed for breakpoint detection and parameter estimation until the number and position of detected breakpoints remains stable.
Considering the frequent land cover changes that have resulted from the active human activities in this study area, we set 1 year as the minimal segment size between successive abrupt changes for fine change detection.Both changes will be detected as long as their interval exceeds 1 year.
The phenological characteristics among land cover categories are different, and they can be effectively captured by the seasonal component.The breakpoint in the seasonal component, which means the phenological change before and after breakpoint, suggests a potential land cover change.In this study, land cover changes related to forest, including conversion from forest to other land cover categories (hereinafter referred to as FTO) and from other land cover categories to forest (hereinafter referred to as OTF), are detected.Results calculated by using BFAST in sample areas of forest, an area converted from forest to shoal (Figure 4a) or from cropland to forest (Figure 4b), indicate the phenological characteristic differences among different land cover categories.Although the forest and shoal have a significant cyclic intra-annual variation with one peak in each year, the seasonal component of the shoal possesses small fluctuations of intra-annual NDVI values.Cropland has an intra-annual variation with two peaks in each year compared with the phenological characteristic of one peak within 1 year for forest.The breakpoints in the seasonal component indicate the structural change of the phenological characteristic and therefore suggest a permanent land cover change.Specifically, the monthly NDVI time series decomposed by BFAST can detect the land cover changes related to forest in the subtropical wetland.
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 33 The phenological characteristics among land cover categories are different, and they can be effectively captured by the seasonal component.The breakpoint in the seasonal component, which means the phenological change before and after breakpoint, suggests a potential land cover change.In this study, land cover changes related to forest, including conversion from forest to other land cover categories (hereinafter referred to as FTO) and from other land cover categories to forest (hereinafter referred to as OTF), are detected.Results calculated by using BFAST in sample areas of forest, an area converted from forest to shoal (Figure 4a) or from cropland to forest (Figure 4b), indicate the phenological characteristic differences among different land cover categories.Although the forest and shoal have a significant cyclic intra-annual variation with one peak in each year, the seasonal component of the shoal possesses small fluctuations of intra-annual NDVI values.Cropland has an intra-annual variation with two peaks in each year compared with the phenological characteristic of one peak within 1 year for forest.The breakpoints in the seasonal component indicate the structural change of the phenological characteristic and therefore suggest a permanent land cover change.Specifically, the monthly NDVI time series decomposed by BFAST can detect the land cover changes related to forest in the subtropical wetland.The short-term disturbances, such as flooding and replantation post-deforestation, result in conversion from forest to forest (hereinafter referred to as FTF).Although these ephemeral changes, including disturbance to forest and then recovery performed within a short term, could not result in the permanent change related to forest, such as the land cover change (FTO or OTF), they still cause short-term forest change and spatiotemporal variation of forest.Therefore, the description of flooding and replantation post-deforestation is as important as the detection of land cover changes related to forest for monitoring the spatiotemporal variation of forest in West Dongting Lake.Results calculated by using BFAST over sample areas of forestland, an area flooding (Figure 4b), or replantation post-deforestation (Figure 4c), showed that these interruptions only cause abrupt changes in the trend component.Two breakpoints (dash lines in Figure 4) closely adjacent to each other in trend component represent the disturbance (the former) and recovery (the latter), respectively.The major difference between flooding and replantation post-deforestation is whether the value of NDVI time series is greater than zero when the disturbance occurred.These results showed that the different ephemeral disturbances, even the recovery process, could also be detected and distinguished by BFAST.

Land Cover Classification after Breakpoint Detection
BFAST could identify abrupt change types in the remotely sensed time series: Abrupt changes in the seasonal component indicate a permanent land cover change; the ephemeral changes, such as flooding and replantation post-deforestation, could be detected in the trend component.However, land cover categories before and after the breakpoints could not be identified by using the BFAST algorithm.Similar to the CCDC algorithm, the coefficients of the time series models (Equations ( 2) and ( 3)) fitted by BFAST are used as input parameters to RFC [52] for land cover classification, which is the main reason why the monthly Landsat time series and BFAST are adopted in this study.The subsequent classification after change detection could also reduce the pseudo change (false positive in change detection).Therefore, the land cover categories at any given time and location could be effectively acquired by integrating BFAST and RF algorithms.
Each pixel has its own time series models (the trend model in Equation ( 2) and seasonal model in Equation ( 3)) before and after any breakpoint.For example, n + 1 time series models exist when n breakpoints are detected in the seasonal and trend components by BFAST.Given that each land cover category has its own specific shape of the time series models [53], the coefficients of each time series model used to describe the shape could be utilized as input parameters to the trained RFC for classifying land cover category.Accordingly, the land cover category of each pixel for each time series model in the entire period is determined.Two coefficients (α and β in Equation ( 2)) in the trend component and six coefficients (γ k , θ k , and K = 3 in Equation ( 3)) in the seasonal one were selected as the classification features in this study.Variable α represents the overall NDVI value of each time series model, and β indicates inter-annual differences, the combination of which describe the mean value and inter-annual trend of each land cover category.Meanwhile, γ k and θ k provide information about intra-annual seasonal change pattern primarily caused by phenological change.
The land cover reference data for classification, a total of 26,137 reference pixels for five land cover categories, were collected on four dates within the period covered by the Landsat time series (Table 3).An approximately equal number (around 6500) of reference pixels were collected on each date.The number of reference pixels for different land cover categories on each date is different, considering the difference among areas, intra-class heterogeneity and data availability of each land cover category.The land cover category of each reference pixel at the given time was determined through careful visual interpretation by using the Landsat time series and the images with high spatial resolutions from Google Earth; and the reference pixels of each land cover category are equably collected at representative regions.Random sampling with replacement design was used to divide the land cover reference data into two parts of training data (for training the classifier) and test data (for validating the classifier) with sampling probabilities of 0.7 and 0.3, respectively.Due to the lack of remotely sensed images with higher temporal and finer spatial resolution than Landsat images over the last few decades covered by Landsat, the primary and optimal source for accuracy assessment of long-term change detection using historical dense time series of Landsat is the Landsat images themselves [21,[54][55][56].Meanwhile, the images with high spatial resolution from Google Earth could help manually interpret the land cover categories and change.A random stratified sample design proposed by Zhu and Woodcock [21] was used to assess the detection accuracy in this study.A total of 600 reference pixels were selected.Three hundred example pixels represent the areas where forest was persistent between 2000 and 2018 and another 300 pixels represent those where forest change occurred, which were identified by the BFAST and RF algorithm.In each pixel, we visually examined the NDVI time series data and images with high spatial resolution acquired from Google Earth to determine whether and when forest change occurred for that specific location.
A randomly selected subset, a total of 30% (7837 reference pixels) of the land cover reference pixels, were used to assess the accuracy of classification with the RF algorithm (Section 2.2.3).A confusion matrix was produced to estimate user's and producer's accuracy of each land cover category, and overall accuracy.

Spatiotemporal Dynamics Analysis of Forest
Eleven indicators were utilized to analyze the spatiotemporal dynamics law of forest in West Dongting Lake region from the aspect of duration, number and timing of abrupt change, and areas (Table 4).D f is used to measure the duration of forest plantation within each pixel between 2000 and 2018.Two indicators (P s , P c ) are utilized to measure the proportion of stable and changed forests, respectively.N c (i), P c (i) are separately used to measure the number and proportion of changed pixels with different numbers of abrupt changes.M L is calculated to illustrate the proportion of abrupt changes in each region, while two indicators (i K , R L ) are adopted to describe the change frequency of each forest pixel and each region, respectively.T l and AP K,T are utilized to describe the timing of the last abrupt change and the annual proportion of all abrupt changes, respectively.In the end, forest areas A L,t in every time phase within time series are calculated.Five indicators (P s , P c , M L , R L , A L,t ) are calculated in each region, while six indicators (i K , N c (i), P c (i), M L , T l , AP K,T ) are calculated based on different change types.The ratio of the number of changes in a certain year to the total number of changes.N c (i, T, K) represents the number of changed pixels with i times change type n f (L, t) represents the number of forest pixels in region L in t. s is the area of one Landsat pixel (900 m 2 )

Accuracy Assessment for Change Detection
For changed pixels, the results showed the producer's accuracy of 90.7% and user's accuracy of 84.3% (Table 5).For stable pixels, the results showed the producer's accuracy of 85.4% and user's accuracy of 91.3%.The change detection with an overall accuracy of 87.8% indicated that the changes related to forest were well captured by BFAST algorithm.The accuracy of the timing of change detection was further assessed for 253 reference pixels in which the change was identified with the Landsat time series data or high spatial resolution images acquired from Google Earth and correctly detected by BFAST (Table 6).The proportion of pixels in which the change timing detected by BFAST is the same as that identified with reference data is 77.47%.The proportions of the pixels in which the detected timing is later than that identified with reference data within 2 months or over 2 months are 13.04% and 9.49%, respectively.If we consider that the actual date of change is within a quarter of a year of the detected change date, the accuracy of timing is as high as 90.51% (77.47% + 13.04%).The reference Landsat images before and after abrupt change and results of change detection with BFAST for each conversion type were displayed to further intuitively illustrate the algorithm ability (Figure 5).A typical area converted from forest to urban in 2008 was found by BFAST algorithm, in which wavelet characteristics of the seasonal component were weakened to a small fluctuation of intra-annual NDVI values after the breakpoint (Figure 5a).The slow urban expansion resulted in partial forest cut within the pixel, in which this spectral change was insufficiently large to be detected by the BFAST algorithm until it was totally changed.Accordingly, the timing of change detected by BFAST was a little later (2 months) than the actual timing observed from the reference data.Another typical area converted from shoal to forest in 2005 was also found by BFAST algorithm, in which the land cover before and after breakpoint and breaking time detected by BFAST and RF algorithms could correspond to those acquired from the reference Landsat images (Figure 5b).Intra-annual changes, such as the flooding and replantation post-deforestation, were difficult to be detected by algorithms, such as Landtrendr, with the frequency of time series of one image per year.The accuracies of flooding and replantation post-deforestation detected by BFAST algorithm were high compared with the reference Landsat images before and after the breakpoint (Figure 5c,d).Moreover, the confidence interval of each detected breaking time (two dashed lines in the trend component) was narrow.

Accuracy Assessment for Classification
The overall accuracy and Kappa coefficient of the land cover classification with five categories are 92.96% and 0.9026, respectively (Table 7).Not only producer accuracies but also user accuracies for the different land cover categories are mostly greater than 90%.These results indicate that the RF classifier is sufficient for classifying the land cover categories in Dongting Lake region.The largest confusion is between forest and shoal, with the producer accuracy for shoal being 73.44%.The possible reason is that other vegetation growing in shoal, such as the annual shrubs or herbs, are similar to forest in the spectral characteristics at the spatial resolution of 30 m.Previous studies mostly focused on the regions that are known to have been forested at the beginning of the first monitoring period and neglected the forest regions in which the forest has not been planted in the initial period [27,57], such as the pixels, the land cover of which changed from another land cover category in the initial period to forest in the later period.Considering the characteristic of frequent change between forest and the other land cover categories in West Donting Lake, we created a forest mask of each image in the time series and merged these masks to produce an 'anytime' forest mask for extracting all forest pixels.The distribution map of the duration proportion of forest plantation between 2000 and 2018 (D f , Table 4) is shown in Figure 6.The forests are mainly distributed in the core zone surrounding the water body and southwestern mountains.The D f of most forest regions is 1 (stable forests), whereas the rest of the forests, the D f of which are less than 1 (changed forests), are dispersedly distributed.Especially, the forest areas with D f less than 0.5 are mainly distributed in the south of the West Donting Lake region, and most changed forests with D f greater than 0.5 are concentrated in the core zone.In contrast with the forests in South China, which are mainly distributed in mountainous regions because of the limited flatlands with low elevation utilized for croplands or settlements, a large part of forests in the West Dongting Lake region are distributed in the flatlands with low elevation, except for the forests in the southwestern mountains (Figures 1e and 6).
Further comparison of the spatial distribution between stable and changed forests for each DEM range and region showed considerable differences (Table 8).The proportions of stable and changed forests (P s and P c , Table 4) are 72.08% and 27.92% for all forest regions, a total of 1,180,893 forest pixels.Over 60% (44.45%/72.08% = 61.67%) of the stable forests are distributed in the transition zone, including those with DEM greater than 50 m (27.55%); the stable forests in the core area are mainly distributed in the region with DEM ranging from 20 to 40 m.This finding indicates that the stable forests are mainly distributed in relatively high elevation regions of the transition zone and low elevation regions of the core zone.The proportion of changed forests with DEM ranging from 30 to 40 m is the highest (11.12%), and those with DEM ranging from 20 to 40 m is over 60% of the changed forests, especially including those in the core zone.In addition, the proportion of changed forests in the transition zone with DEM greater than 50 m is also high (6.88%).These results showed that the changes related to forest were mainly detected in the low elevation region surrounding the water body and relatively high elevation region of the transition zone.Although the proportion of forests in the core zone (11.69% + 27.63% = 39.32%) is lower than that in the transition zone (16.23% + 44.45% = 60.68%), the proportion of changed forests in all of the forests in the former (11.69%/39.32%= 29.73%)was slightly higher than that in the latter (16.23%/60.68%=26.75%).Meanwhile, the proportion of changed forests in the region with each DEM range decreased with the DEM increase.These results further showed that the forest changes in the core zone and low elevation region are more prevalent.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 33 forests, especially including those in the core zone.In addition, the proportion of changed forests in the transition zone with DEM greater than 50 m is also high (6.88%).These results showed that the changes related to forest were mainly detected in the low elevation region surrounding the water body and relatively high elevation region of the transition zone.Although the proportion of forests in the core zone (11.69% + 27.63% = 39.32%) is lower than that in the transition zone (16.23% + 44.45% = 60.68%), the proportion of changed forests in all of the forests in the former (11.69%/39.32%= 29.73%)was slightly higher than that in the latter (16.23%/60.68%=26.75%).Meanwhile, the proportion of changed forests in the region with each DEM range decreased with the DEM increase.These results further showed that the forest changes in the core zone and low elevation region are more prevalent.Df is calculated as the ratio of the number of time phases when the land cover category within the pixel is forest to the total number of time phases (or images) within the monthly Landsat time series.All the forest pixels were extracted by using the 'anytime' forest mask, and then the proportion and spatial distribution of stable and changed forests were also analyzed in this section.However, the spatiotemporal characteristic of forest changes must be analyzed to further explain the variation characteristic of forests.Therefore, the number and timing of abrupt changes related to forest, and the variation of forest areas will be discussed in the following sections.

Number of Abrupt Changes Related to Forest
The number (Nc(i), Table 4) and proportion (Pc(i), Table 4) of pixels with single (or multiple) occurrence(s) (iK, Table 4) of abrupt change for each change type were compared to show considerable differences (Table 9).In all the changed pixels related to forest (27.92%), the number of abrupt changes of 25.35% out of 27.92%, which exceeds 90% (25.35%/27.92%= 90.80%) of all the changed pixels related to forest, is from 1 to 3. The number of FTO and OTF, which separately occurred in  All the forest pixels were extracted by using the 'anytime' forest mask, and then the proportion and spatial distribution of stable and changed forests were also analyzed in this section.However, the spatiotemporal characteristic of forest changes must be analyzed to further explain the variation characteristic of forests.Therefore, the number and timing of abrupt changes related to forest, and the variation of forest areas will be discussed in the following sections.

Number of Abrupt Changes Related to Forest
The number (N c (i), Table 4) and proportion (P c (i), Table 4) of pixels with single (or multiple) occurrence(s) (i K , Table 4) of abrupt change for each change type were compared to show considerable differences (Table 9).In all the changed pixels related to forest (27.92%), the number of abrupt changes of 25.35% out of 27.92%, which exceeds 90% (25.35%/27.92%= 90.80%) of all the changed pixels related to forest, is from 1 to 3. The number of FTO and OTF, which separately occurred in 13.36% and 11.89% of the forest pixels, is basically 1.The proportion of the forest pixels with FTF is 15.70%.In contrast with FTO and OTF, the proportion of forest pixels with one time FTF (flooded or replantation post-deforestation) is generally equal to that with multiple times FTF.Results indicated that not only the proportion of forest pixels with FTF but also the number of FTF is higher than those of FTO and OTF.The distribution map of the number of breakpoints is shown in Figure 7.Most areas with two or more times changes are in the core zone surrounding the water body, and the areas with one change are dispersedly distributed in the whole study area, especially the southwest mountains (Figure 7a,b).FTO or OTF, the number of which is generally 1, are principally distributed in the areas surrounding the water body (Figure 7c-f).The distribution of areas with different numbers of FTF is almost the same as that with the estimated number of all abrupt changes (Figure 7g,h).
The number of abrupt changes is summarized for each DEM range and region to facilitate further comparison (Figure 8).Most changes occurred in the areas with DEM ranging from 20 to 40 m, the proportion (M L , Table 4) of which in all of detected changes is over 65% (Figure 8a).In these areas with DEM ranging from 20 to 40 m, the proportion of detected changes in the core zone is higher than that in the transition zone regardless of the change types (Figure 8b-d).The detected changes in the region with DEM over 40 m are generally distributed in the transition zone, the proportion of which is much less than that with DEM ranging from 20 to 40 m.This is because the elevation in West Dongting Lake is generally less than 40 m and the main land cover categories in the transition zone with DEM less than 40 m are cropland and settlement (Figures 1 and 6).The proportion of OTF (18.59%) and FTO (18.59%) in the transition zone with DEM over 50 m (southwest mountains) is relatively low compared with the proportion of OTF (42.89%) and FTO (42.78%) in the core zone with DEM less than 40 m (Figure 8b,c).The proportion of FTF in the southwest mountains (24.64%) is also lower than that in the core zone with DEM less than 40 m (43.33%; Figure 8d).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 33 13.36% and 11.89% of the forest pixels, is basically 1.The proportion of the forest pixels with FTF is 15.70%.In contrast with FTO and OTF, the proportion of forest pixels with one time FTF (flooded or replantation post-deforestation) is generally equal to that with multiple times FTF.Results indicated that not only the proportion of forest pixels with FTF but also the number of FTF is higher than those of FTO and OTF.The distribution map of the number of breakpoints is shown in Figure 7.Most areas with two or more times changes are in the core zone surrounding the water body, and the areas with one change are dispersedly distributed in the whole study area, especially the southwest mountains (Figure 7a,b).FTO or OTF, the number of which is generally 1, are principally distributed in the areas surrounding the water body (Figure 7c-f).The distribution of areas with different numbers of FTF is almost the same as that with the estimated number of all abrupt changes (Figure 7g,h).
The number of abrupt changes is summarized for each DEM range and region to facilitate further comparison (Figure 8).Most changes occurred in the areas with DEM ranging from 20 to 40 m, the proportion (ML, Table 4) of which in all of detected changes is over 65% (Figure 8a).In these areas with DEM ranging from 20 to 40 m, the proportion of detected changes in the core zone is higher than that in the transition zone regardless of the change types (Figure 8b-d).The detected changes in the region with DEM over 40 m are generally distributed in the transition zone, the proportion of which is much less than that with DEM ranging from 20 to 40 m.This is because the elevation in West Dongting Lake is generally less than 40 m and the main land cover categories in the transition zone with DEM less than 40 m are cropland and settlement (Figures 1 and 6).The proportion of OTF (18.59%) and FTO (18.59%) in the transition zone with DEM over 50 m (southwest mountains) is relatively low compared with the proportion of OTF (42.89%) and FTO (42.78%) in the core zone with DEM less than 40 m (Figure 8b,c).The proportion of FTF in the southwest mountains (24.64%) is also lower than that in the core zone with DEM less than 40 m (43.33%; Figure 8d).The ratio (RL, Table 4) of the total number of changes to the total number of changed forest pixels was proposed to further reflect the frequency of forest change in different regions.Although the proportion of change number in the core zone is slightly less than that in the transition zone (47.49% versus 52.51% of 665,744), the ratio in the core zone (2.29) is higher than that in the transition zone (1.82) considering the proportion of the changed forests in both zones (11.69% versus 16.23% of 1,180,893, Table 8).Moreover, the ratio in the core zone is slightly higher than that in the transition zone regardless of the change types, especially FTF (2.15 versus 1.66).Specifically, the proportion of  The ratio (RL, Table 4) of the total number of changes to the total number of changed forest pixels was proposed to further reflect the frequency of forest change in different regions.Although the proportion of change number in the core zone is slightly less than that in the transition zone (47.49% versus 52.51% of 665,744), the ratio in the core zone (2.29) is higher than that in the transition zone (1.82) considering the proportion of the changed forests in both zones (11.69% versus 16.23% of 1,180,893, Table 8).Moreover, the ratio in the core zone is slightly higher than that in the transition zone regardless of the change types, especially FTF (2.15 versus 1.66).Specifically, the proportion of The ratio (R L , Table 4) of the total number of changes to the total number of changed forest pixels was proposed to further reflect the frequency of forest change in different regions.Although the proportion of change number in the core zone is slightly less than that in the transition zone (47.49% versus 52.51% of 665,744), the ratio in the core zone (2.29) is higher than that in the transition zone (1.82) considering the proportion of the changed forests in both zones (11.69% versus 16.23% of 1,180,893, Table 8).Moreover, the ratio in the core zone is slightly higher than that in the transition zone regardless of the change types, especially FTF (2.15 versus 1.66).Specifically, the proportion of changed forests (Section 3.2.1)and the ratio in the core zone are higher than those in the transition zone.This finding further indicates that the forest change in the core zone is more frequent than that in the transition zone.

Timing of Abrupt Changes Related to Forest
The spatial distribution of the timing (T l , Table 4) of the last abrupt change detected from the monthly NDVI time-series is shown in Figure 9.Most FTO in the core zone occurred between 2011 and 2015, whereas that in the transition zone is between 2012 and 2017 (Figure 9a,b).The major years of OTF in the core zone region are 2000, 2012, and 2016, while most transition zones had OTF in 2012, 2016, and 2017 (Figure 9c,d).FTF mainly occurred in 2017 for the southwest mountains, while large parts of the core zone region had FTF between 2006 and 2017 (Figure 9e,f).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 33 changed forests (Section 3.2.1)and the ratio in the core zone are higher than those in the transition zone.This finding further indicates that the forest change in the core zone is more frequent than that in the transition zone.

Timing of Abrupt Changes Related to Forest
The spatial distribution of the timing (Tl, Table 4) of the last abrupt change detected from the monthly NDVI time-series is shown in Figure 9.Most FTO in the core zone occurred between 2011 and 2015, whereas that in the transition zone is between 2012 and 2017 (Figure 9a    The changed forest pixels with nonFTO (or nonOTF, or nonFTF) represent that no FTO (or OTF, FTF) occurred in these changed pixels.
In contrast with the analysis of the timing of last change in Figure 9, the timing of all changes was summarized each year (APK,T, Table 4) to facilitate further comparison (Figure 10).Large parts of FTO occurred since 2011, while a relatively low proportion of FTO was detected before 2011 (Figure 10a).Several clear peaks of OTF were detected during 2000, 2012, 2016, and 2017 (Figure 10b).The changed forest pixels with nonFTO (or nonOTF, or nonFTF) represent that no FTO (or OTF, FTF) occurred in these changed pixels.
In contrast with the analysis of the timing of last change in Figure 9, the timing of all changes was summarized each year (AP K,T , Table 4) to facilitate further comparison (Figure 10).Large parts of FTO occurred since 2011, while a relatively low proportion of FTO was detected before 2011 (Figure 10a).Several clear peaks of OTF were detected during 2000, 2012, 2016, and 2017 (Figure 10b).The frequency distribution of the timing of FTF was relatively even for each year from 2000 to 2018, except for one peak in 2017 (Figure 10c).
The frequency distribution of the timing of FTF was relatively even for each year from 2000 to 2018, except for one peak in 2017 (Figure 10c).

Variation of Forest Areas
The spatial distribution of forest pixels and the forest areas (AL,t, Table 4) at any given time could be effectively acquired with the application of integrating BFAST and RF algorithms to the monthly NDVI time series (Figure 11).The variation of forest areas in West Dongting Lake experienced three steps (Figure 11a The increased forest areas in the core zone are significantly higher than that in the transition zone during the period of 2000 through 2005 (Figure 11b).This finding indicates that forest plantation in the core zone is the major contribution for increased forest areas between 2000 and 2005.However, the increasing rate of forest area in the transition zone is higher than that in the core zone since 2003.The forest areas in the core zone generally remained steady between 2003 and 2011.The forest areas in the transition zone first also had the same trend as that in the core zone between 2006 and 2008 and then increased between 2008 and 2011.The overall trend of the forest areas in the core zone

Variation of Forest Areas
The spatial distribution of forest pixels and the forest areas (A L,t , Table 4) at any given time could be effectively acquired with the application of integrating BFAST and RF algorithms to the monthly NDVI time series (Figure 11).The variation of forest areas in West Dongting Lake experienced three steps (Figure 11a The increased forest areas in the core zone are significantly higher than that in the transition zone during the period of 2000 through 2005 (Figure 11b).This finding indicates that forest plantation in the core zone is the major contribution for increased forest areas between 2000 and 2005.However, the increasing rate of forest area in the transition zone is higher than that in the core zone since

Discussion
The accuracy assessment results of change detection show higher user accuracy for stable pixels and higher producer accuracy for change pixels (Table 5), indicating more commission errors than omission errors in the detected change pixels.The errors are possibly due to the following reasons: (1) noise; (2) partial change; (3) and influence of classification accuracy.Although the validity of the NSPI, MNSPI and ESTARFM are verified in this study and our previous experiment [44], the noise cannot be removed completely, which is mistakenly detected as change (commission error).The partially changed pixels are difficult to be detected by BFAST due to the insufficiently large spectral change (omission error).The similarity between shoal covered by other vegetation and forest in the spectral characteristics within the Landsat image with 30 m resolution may result in misclassification, and then the inter-conversion between them cannot be identified even if the conversions have been detected by BFAST before classification (omission error).
The timing errors of change detection are possibly due to the following reasons: (1) slow change; (2) reconstruction effect of time series; and (3) parameter setting of BFAST.As shown in Figure 5a, the slow urban expansion resulted in a partial forest cut within the pixel, in which this spectral change was insufficiently large to be detected by the BFAST algorithm until it was totally changed.Accordingly, the timing of change detected by BFAST was a little later (2 months) than the actual timing observed from the reference data.NSPI, MNSPI and ESTARFM algorithms interpolated the contaminated pixels and reconstructed the missing images without considering possible land cover change at that time phase, which means that the land cover category of the interpolated or reconstructed pixel with actual land cover change may have mistakenly remain unchanged at that time phase.This situation may have resulted in the land cover change detected timing delay of one or several months.We set 1 year as the minimal segment size between successive change detections in BFAST to reduce pseudo changes.However, this strategy caused the interval of two breakpoints closely adjacent to each other in the trend component, representing the short-term disturbances (such as flooding and replantation post-deforestation), to be just 1 year.In fact, the cycle of the short-term disturbances is generally less than 1 year, and it may delay the timing of the recovery process.
On the basis of the error analysis, there are some feasible strategies as follows: (1) change detection through integrating temporal and spatial information and (2) classification using deep learning.Similar to most change detection methods utilizing time-series data, the BFAST adopted in this study is mainly focused on the information on the temporal domain and independently processes each pixel, without considering the change characteristics generated on the spatial domain [14,58].This results

Discussion
The accuracy assessment results of change detection show higher user accuracy for stable pixels and higher producer accuracy for change pixels (Table 5), indicating more commission errors than omission errors in the detected change pixels.The errors are possibly due to the following reasons: (1) noise; (2) partial change; (3) and influence of classification accuracy.Although the validity of the NSPI, MNSPI and ESTARFM are verified in this study and our previous experiment [44], the noise cannot be removed completely, which is mistakenly detected as change (commission error).The partially changed pixels are difficult to be detected by BFAST due to the insufficiently large spectral change (omission error).The similarity between shoal covered by other vegetation and forest in the spectral characteristics within the Landsat image with 30 m resolution may result in misclassification, and then the inter-conversion between them cannot be identified even if the conversions have been detected by BFAST before classification (omission error).
The timing errors of change detection are possibly due to the following reasons: (1) slow change; (2) reconstruction effect of time series; and (3) parameter setting of BFAST.As shown in Figure 5a, the slow urban expansion resulted in a partial forest cut within the pixel, in which this spectral change was insufficiently large to be detected by the BFAST algorithm until it was totally changed.Accordingly, the timing of change detected by BFAST was a little later (2 months) than the actual timing observed from the reference data.NSPI, MNSPI and ESTARFM algorithms interpolated the contaminated pixels and reconstructed the missing images without considering possible land cover change at that time phase, which means that the land cover category of the interpolated or reconstructed pixel with actual land cover change may have mistakenly remain unchanged at that time phase.This situation may have resulted in the land cover change detected timing delay of one or several months.We set 1 year as the minimal segment size between successive change detections in BFAST to reduce pseudo changes.However, this strategy caused the interval of two breakpoints closely adjacent to each other in the trend component, representing the short-term disturbances (such as flooding and replantation post-deforestation), to be just 1 year.In fact, the cycle of the short-term disturbances is generally less than 1 year, and it may delay the timing of the recovery process.
On the basis of the error analysis, there are some feasible strategies as follows: (1) change detection through integrating temporal and spatial information and (2) classification using deep learning.Similar to most change detection methods utilizing time-series data, the BFAST adopted in this study is mainly focused on the information on the temporal domain and independently processes each pixel, without considering the change characteristics generated on the spatial domain [14,58].This results in the accuracy of change detection relying heavily on the quantity and quality of time-series data in the temporal domain, such as the reconstruction accuracy of the monthly Landsat time series in this study.In fact, the change essentially occurs in the spatial domain and is reflected in the temporal domain, in other words, change determines the necessity of defining 'time' by human.Therefore, time series analysis based on an independent pixel is insufficient to accurately and roundly describe a series of change characteristics.Several algorithms recently eliminated the noise and other change (such as seasonal change) interference through integrating spatial information, which improved the change detection accuracy [18,[59][60][61][62].These methods have less reliance on the quantity and quality of time-series data.For example, Hamunyela et al. used spatially normalized NDVI values on the basis of setting the various window sizes to reduce phenological differences and lower latency in change detection [59], and Guttler et al. developed the object-based change detection method [62].They may be used to resolve the problem mentioned in the previous paragraph.Therefore, additional algorithms that use the spatiotemporal data analytical techniques (such as spatiotemporal geo-statistic methods) [63,64], which consider the correlativity and heterogeneity of spatiotemporal data, are anticipated.In addition, the improvement of the classification accuracy can also contribute to change detection.Deep learning methods have recently progressed considerably and have played a key role in addressing remote sensing image classification [65,66].The validity of classification using Convolutional Neural Network (CNN) in South China has been verified in our previous studies [67][68][69].Compared with traditional classification methods, including support vector machine and decision tree, the CNN improves the accuracy of classification through extracting deep features of spatial and spectral trajectories.Specifically, the abandoned land and even paddy rice field with different rice-cropping systems are classified well.
Although the West Dongting Lake region is small in spatial scale, the generation, storage, and processing of the monthly Landsat Time series between 2000 and 2018 and the change detection with BFAST still requires expensive storage and computation.Therefore, time series data analysis in the large regional scale (such as global scale) and long temporal scale (such as the 47-year archive of Landsat) needs more expensive storage and computation because of the large amount of data and its processing [21].Computing capability, such as the cloud-computing technology with the Google Earth Engine (GEE) platform [70], has been greatly improved in the last few years.GEE synchronously archives the Landsat data from USGS and effectively performs data processing by using the cloud computing technology through utilizing millions of servers worldwide.This technology shows potentials and prospects of the emerging GEE platforms in large regional and long temporal scales for land cover change research, and several studies have proved its feasibility [31,71,72].Therefore, the data storage, processing, and analysis will be conducted by using the GEE platform in the next exploration.
Although the proportion of forests in the core zone is lower than that in the transition zone, the proportion of changed forests in all the forests in the former is higher than that in the latter, especially in the areas with low elevation.The forests in the core zone are mainly distributed in the region with DEM ranging from 20 to 40 m, while over two-thirds of the forests in the transition zone were planted in the region with DEM exceeding 40 m (Table 8), and these forests are concentrated in the southwestern mountain with DEM over 50 m.Therefore, the forests with low elevation surrounding the water body in the core zone are easily flooded, which increases the number of FTF.In addition, the variation of timber price caused other land cover categories, such as cropland or shoal, to be switched to forest and replantation post-deforestation when prices rose, or the forest to be switched to other land cover when the price fell because of economic profit motivation, increasing the number of OTF or FTO.Meanwhile, the southwestern mountain with high elevation is suitable for tree plantation and unsuitable for growing crop and settlement, which leads to a relatively low proportion of OTF and FTO compared with the proportion of OTF and FTO in the core zone.These situations result in a high change frequency in the core zone.
The low grain price and high timber price, policy incentives and market demand resulted in a rapid expansion of forest plantation in the West Dongting Lake region from 2000 to 2005.The early expansion mainly occurred in the core zone, and then in the transition zone since 2003.One potential factor is the transferal of parts of forest plantation from the core zone surrounding the water body to the transition zone due to TGD raising the water level of Dongting Lake.A possible reason for forest areas keeping relatively steady from 2006 to 2011 is a recovery in the grain price and stable timber price.Since 2011, a series of policies for protecting wetlands and restoring ecology made by the local government, especially the project of returning forest to wetland after 2013, forbade poplar plantation in the Dongting Lake region.It is the key driver of the difference of changes between before and after 2011.These policies have led to a continuing decline of the forest areas since 2011, except for the forest areas in the core zone sharply fluctuating with the fluctuating price of timber in 2014 and 2016.The policy of cleaning the poplar plantation in 2017 further accelerated the transferal of forest plantation from the core zone to the transition zone, resulting in a decline of forest areas in the core zone and an increasing trend in the transition zone during 2017.It also explained the high proportions of OTF and FTF in 2017 (Figure 10b,c).Therefore, the variation trend of forest areas indicated that anthropogenic factors, such as government policies and economic profits, are the dominating drivers of forest plantation in the West Dongting Lake region.
Results in this study effectively showed not only the final results of forest variations, but also the spatiotemporal dynamic process of forests and the potential drivers in West Dongting Lake, providing abundant knowledge for wetland ecosystem protection.

Conclusions
This study confirmed the applicability of the BFAST method to detect changes related to forest in the subtropical wetland ecosystem using the monthly Landsat NDVI time series and investigated the spatiotemporal variations of forest in the West Dongting Lake region between 2000 and 2018.By using the monthly Landsat time series generated by the spatiotemporal interpolator and data fusion algorithms, the forest changes, including FTO, OTF, and FTF, and land cover categories before and after breakpoints, could be effectively captured by integrating BFAST and RF algorithms, with an overall accuracy of 87.8%.The spatiotemporal variations of forest were described from the aspects of duration and distribution of forest, the number and timing of forest changes, and the variation of forest areas, respectively: (1) The forests are mainly distributed in the core zone surrounding the water body and the southwestern mountains.Over 60% of the stable forests are distributed in the transition zone, whereas the changed forests in the core zone and low elevation region are more prevalent.(2) Not only the proportion of pixels with FTF but also the number of FTF is higher than those of FTO and OTF.The forest change in the core zone is more frequent than that in the transition zone.Most changes occurred in the areas with DEM ranging from 20 to 40 m, in which the proportion of detected change number in the core zone is higher than that in the transition zone regardless of the change type.(3) Large parts of FTO occurred since 2011 due to the policy for protecting wetlands and restoring ecology.Some clear peaks of OTF were detected during 2000, 2012, 2016, and 2017 because of the jumping price of timber.FTF was relatively even for each year, except for one peak in 2017.(4) The variation of forest areas in the West Dongting Lake region experienced three steps: rapid expansion of forest plantation from 2000 to 2005, to which forest plantation in the core zone is the major contribution; relatively steady from 2006 to 2011, except for a small increase in the transition zone between 2008 and 2011; and continuous decline since 2011 with small sharp fluctuations in the core zone and a laggard downtrend in the transition zone.These variations were mainly caused by anthropogenic factors, such as the government policies and economic profits.

Figure 1 .
Figure 1.Geographic location, Landsat and DEM image of the study area: (a) general location within China, (b) within Hunan Province, and (c) within the Dongting Lake region; (d) Landsat 8 OLI image (RGB = near-infrared, red, and green band); and (e) Global Digital Elevation Model (GDEM) image of the study area.

Figure 1 .
Figure 1.Geographic location, Landsat and DEM image of the study area: (a) general location within China, (b) within Hunan Province, and (c) within the Dongting Lake region; (d) Landsat 8 OLI image (RGB = near-infrared, red, and green band); and (e) Global Digital Elevation Model (GDEM) image of the study area.

Figure 2 .
Figure 2. Flowchart of investigating the spatiotemporal variations of forest in the West Dongting Lake region.NSPI, MNSPI and ESTARFM represent the spatiotemporal interpolation method named Neighborhood Similar Pixel Interpolator and its modified version, and a spatiotemporal data fusion algorithm named enhanced spatial and temporal adaptive reflectance fusion model, respectively.BFAST is the change detection algorithm named Breaks For Additive Seasonal and Trend, and RF represents random forest algorithm.

Figure 3 .
Figure 3. Proportions of 197 available Landsat images with different percentages of contaminated pixels between 2000 and 2018 in the West Donting Lake region.

Figure 3 .
Figure 3. Proportions of 197 available Landsat images with different percentages of contaminated pixels between 2000 and 2018 in the West Donting Lake region.

Figure 4 .
Figure 4. Detected changes in the seasonal and trend component of monthly NDVI time series between 2000 and 2018 by using BFAST: (a) The breakpoint in the seasonal component St represents the conversion from forest to other land cover categories (FTO); (b) the breakpoint in the seasonal component St represents the conversion from other land cover categories to forest (OTF), and two breakpoints in the trend component Tt represent flooding and recovery (FTF); (c) two breakpoints closely adjacent to each other in the trend component Tt represent replantation post-deforestation (FTF).Yt represents the original observed data, and e t represents the remainder component.The dash line represents the breakpoint and indicates when abrupt change occurred.The confidence intervals of the estimated date of abrupt changes are shown (red).
represents the number of changed pixels with i times abrupt changes in region L n i=1 i × N c (i) represents the total number of abrupt changes The ratio R L of the number of abrupt changes to the number of changed pixels in region L R L = n i=1 i×N c (i,L) N c (L) N c (L) represents the number of change pixels in region L Timing Timing (T l ) of last change [25] T l = 2000, 2001, . . ., 2018 Annual proportion AP K,T of abrupt changes [24] AP K,T = n i=1 i×N c (i,T,K) n i=1 i×N c (i) T = 2000, 2001, . . ., 2018

Figure 5 .
Figure 5. Examples of accuracy assessment for the detected changes and corresponding Landsat images (RGB = near-infrared, red, and green band) before and after abrupt changes: (a) Forest to other land cover (urban), (b) Other land cover (shoal) to Forest, (c) Flooding, and (d) Replantation post-deforestation.Yt is the original observed data; Tt, St, and et separately represent the trend, seasonal, and remainder component decomposed by BFAST.The dashed line represents the breakpoint and indicates when abrupt change occurred.The confidence intervals of the estimated date of abrupt changes are shown (red).

Figure 5 .
Figure 5. Examples of accuracy assessment for the detected changes and corresponding Landsat images (RGB = near-infrared, red, and green band) before and after abrupt changes: (a) Forest to other land cover (urban), (b) Other land cover (shoal) to Forest, (c) Flooding, and (d) Replantation post-deforestation.Y t is the original observed data; T t , S t , and e t separately represent the trend, seasonal, and remainder component decomposed by BFAST.The dashed line represents the breakpoint and indicates when abrupt change occurred.The confidence intervals of the estimated date of abrupt changes are shown (red).

Figure 6 .
Figure 6.Distribution map of the duration proportion (Df) of forest plantation between 2000 and 2018.Df is calculated as the ratio of the number of time phases when the land cover category within the pixel is forest to the total number of time phases (or images) within the monthly Landsat time series.

Figure 6 .
Figure 6.Distribution map of the duration proportion (D f ) of forest plantation between 2000 and 2018.D f is calculated as the ratio of the number of time phases when the land cover category within the pixel is forest to the total number of time phases (or images) within the monthly Landsat time series.

Figure 7 .
Figure 7. Spatial distribution of the number (iK) of abrupt changes between 2000 and 2018 for different forest change types and their close-up subsection: (a,b) all changes; (c,d) conversion from forest to other land cover categories (FTO); (e,f) conversion from other land cover categories to forest (OTF); and (g,h) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).

Figure 8 .
Figure 8. Proportion (ML) of the number of different forest change types for each DEM range and region: (a) all changes; (b) conversion from forest to other land cover categories (FTO); (c) conversion from other land cover categories to forest (OTF); and (d) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).ML is calculated as the ratio of the number of abrupt changes in a certain region or DEM range to the total number of abrupt changes in the whole study area between 2000 and 2018 for each change type.

Figure 7 . 33 Figure 7 .
Figure 7. Spatial distribution of the number (i K ) of abrupt changes between 2000 and 2018 for different forest change types and their close-up subsection: (a,b) all changes; (c,d) conversion from forest to other land cover categories (FTO); (e,f) conversion from other land cover categories to forest (OTF); and (g,h) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).

Figure 8 .
Figure 8. Proportion (ML) of the number of different forest change types for each DEM range and region: (a) all changes; (b) conversion from forest to other land cover categories (FTO); (c) conversion from other land cover categories to forest (OTF); and (d) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).ML is calculated as the ratio of the number of abrupt changes in a certain region or DEM range to the total number of abrupt changes in the whole study area between 2000 and 2018 for each change type.

Figure 8 .
Figure 8. Proportion (M L ) of the number of different forest change types for each DEM range and region: (a) all changes; (b) conversion from forest to other land cover categories (FTO); (c) conversion from other land cover categories to forest (OTF); and (d) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).M L is calculated as the ratio of the number of abrupt changes in a certain region or DEM range to the total number of abrupt changes in the whole study area between 2000 and 2018 for each change type.
,b).The major years of OTF in the core zone region are 2000, 2012, and 2016, while most transition zones had OTF in 2012, 2016, and 2017 (Figure 9c,d).FTF mainly occurred in 2017 for the southwest mountains, while large parts of the core zone region had FTF between 2006 and 2017 (Figure 9e,f).

Figure 9 .
Figure 9. Spatial distribution of the timing (Tl) of last change of different forest change types and their close up of the subsection detected from monthly NDVI time-series: (a,b) conversion from forest to other land cover categories (FTO),(c,d) conversion from other land cover categories to forest (OTF), and (e,f) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).The changed forest pixels with nonFTO (or nonOTF, or nonFTF) represent that no FTO (or OTF, FTF) occurred in these changed pixels.

Figure 9 .
Figure 9. Spatial distribution of the timing (T l ) of last change of different forest change types and their close up of the subsection detected from monthly NDVI time-series: (a,b) conversion from forest to other land cover categories (FTO),(c,d) conversion from other land cover categories to forest (OTF), and (e,f) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).The changed forest pixels with nonFTO (or nonOTF, or nonFTF) represent that no FTO (or OTF, FTF) occurred in these changed pixels.

Figure 10 .
Figure 10.Annual proportion (APK,T) of different forest change types detected from Landsat NDVI time series between 2000 and 2018: (a) conversion from forest to other land cover categories (FTO), (b) conversion from other land cover categories to forest (OTF), and (c) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).APK,T is calculated as the ratio of the number of changes in a certain year to the total number of changes between 2000 and 2018 for each change type.
): a rapid expansion of forest plantation from 2000 to 2005, except for a small drop between 2001 and 2002; the forest areas were relatively steady from 2006 to 2011, except for one valley in 2008; the forest areas continued to decline since 2011.

Figure 10 .
Figure 10.Annual proportion (AP K,T ) of different forest change types detected from Landsat NDVI time series between 2000 and 2018: (a) conversion from forest to other land cover categories (FTO), (b) conversion from other land cover categories to forest (OTF), and (c) conversion from forest to forest (FTF, such as flooding and replantation post-deforestation).AP K,T is calculated as the ratio of the number of changes in a certain year to the total number of changes between 2000 and 2018 for each change type.
): a rapid expansion of forest plantation from 2000 to 2005, except for a small drop between 2001 and 2002; the forest areas were relatively steady from 2006 to 2011, except for one valley in 2008; the forest areas continued to decline since 2011.
2003.The forest areas in the core zone generally remained steady between 2003 and 2011.The forest areas in the transition zone first also had the same trend as that in the core zone between 2006 and 2008 and then increased between 2008 and 2011.The overall trend of the forest areas in the core zone began to fall since 2011, except for the forest areas in the core zone sharply fluctuating in 2014 and 2016.The forest areas in the transition zone began to fall since 2013, which is laggard compared with the downtrend in the core zone; the decreasing trend and rate of forest areas in the transition zone almost did not change.The forest areas in the core zone sharply decreased during 2017.By contrast, the forest areas in the transition zone showed an increasing trend in 2017.Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 33 began to fall since 2011, except for the forest areas in the core zone sharply fluctuating in 2014 and 2016.The forest areas in the transition zone began to fall since 2013, which is laggard compared with the downtrend in the core zone; the decreasing trend and rate of forest areas in the transition zone almost did not change.The forest areas in the core zone sharply decreased during 2017.By contrast, the forest areas in the transition zone showed an increasing trend in 2017.

Figure 11 .
Figure 11.Variation of forest areas (AL,t) between 2000 and 2018 for (a) entire study area and (b) core and transition zones.

Figure 11 .
Figure 11.Variation of forest areas (A L,t ) between 2000 and 2018 for (a) entire study area and (b) core and transition zones.

Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 33 Figure 2. Flowchart
of investigating the spatiotemporal variations of forest in the West Dongting Lake region.NSPI, MNSPI and ESTARFM represent the spatiotemporal interpolation method named Neighborhood Similar Pixel Interpolator and its modified version, and a spatiotemporal data fusion algorithm named enhanced spatial and temporal adaptive reflectance fusion model, respectively.BFAST is the change detection algorithm named Breaks For Additive Seasonal and Trend, and RF represents random forest algorithm.

Table 1 .
Number of available and missing images for generating a monthly Landsat time series between 2000 and 2018.

Table 2 .
Datasets used in this study.All 197 available Landsat images (Table1) used as inputs to NSPI and MNSPI for SLC-off gaps filling and cloud removal.One-hundred-and-sixty-nine out of 197 images and 59 prediction images reconstructed with ESTARFM, a total of 228 images, utilized to generate a monthly Landsat time series for change detection with BFAST.

Table 3 .
Land cover reference data collected at four dates for classification.

Table 4 .
Indicators for describing spatiotemporal dynamics of forest.
c (i) and proportion P c (i) of changed pixels with i times abrupt changes

Table 5 .
Accuracy assessment of change detection.

Table 6 .
Accuracy assessment of the timing of change detection.

Table 7 .
Accuracy assessment of classification with the RF algorithm.

Table 8 .
Proportion (%) of stable (PS) and changed (Pc) forests for each DEM range and region.

Table 8 .
Proportion (%) of stable (P S ) and changed (P c ) forests for each DEM range and region.

Table 9 .
Number (N c (i)) and proportion (P c (i)) of pixels with single (or multiple) occurrence(s) of abrupt change between 2000 and 2018.
1 TCP represents Total of changed pixels, and 2 TNB represents Total number of breakpoints.

Table 9 .
Number (Nc(i)) and proportion (Pc(i)) of pixels with single (or multiple) occurrence(s) of abrupt change between 2000 and 2018.
1 TCP represents Total of changed pixels, and 2 TNB represents Total number of breakpoints.