Examining Land Cover and Greenness Dynamics in Hangzhou Bay in 1985 – 2016 Using Landsat Time-Series Data

Land cover changes significantly influence vegetation greenness in different regions. Dense Landsat time series stacks provide unique opportunity to analyze land cover change and vegetation greenness trends at finer spatial scale. In the past three decades, large reclamation activities have greatly changed land cover and vegetation growth of coastal areas. However, rarely has research investigated these frequently changed coastal areas. In this study, Landsat Normalized Difference Vegetation Index time series (1984–2016) data and the Breaks For Additive Seasonal and Trend algorithm were used to detect the intensity and dates of abrupt changes in a typical coastal area—Hangzhou Bay, China. The prior and posterior land cover categories of each change were classified using phenology information through a Random Forest model. The impacts of land cover change on vegetation greenness trends of the inland and reclaimed areas were analyzed through distinguishing gradual and abrupt changes. The results showed that the intensity and date of land cover change were detected successfully with overall accuracies of 88.7% and 86.1%, respectively. The continuous land cover dynamics were retrieved accurately with an overall accuracy of 91.0% for ten land cover classifications. Coastal reclamation did not alleviate local cropland occupation, but prompted the vegetation greenness of the reclaimed area. Most of the inland area showed a browning trend. The main contributors to the greenness and browning trends were also quantified. These findings will help the natural resource management community generate better understanding of coastal reclamation and make better management decisions.


Introduction
Land cover change is one of the most important factors profoundly affecting the earth's ecological systems [1,2].Land cover changes can be related to natural processes (e.g., flooding, wildfire) and anthropogenic activities (e.g., urbanization, agriculture).The rate of change and the nature of land cover changes can also differ in time and space [3].Some regions are relatively stable (e.g., permanent forest), whereas other areas are subject to rapid and persistent transformation (e.g., urban expansion of previously vegetated areas).Previous studies showed that land cover change can lead to either increases or decreases in vegetation greenness [4][5][6].For example, the changes from agricultural land or forests to developed land usually result in a decrease in vegetation greenness [4].Conversely, returning agricultural land to forest and grassland, reforestation, and afforestation can lead to increases in greenness [5].As a key component of terrestrial ecosystems, vegetation has an irreplaceable role in regulating material and energy cycles [7,8].The associations between vegetation dynamics and land cover change have attracted great attention in recent decades [5,6,[8][9][10][11].However, most research was conducted at coarse spatial resolution by analyzing time series of satellite data (such as Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS)) at global and national levels [6,12,13].It is hard to accurately detect land cover changes occurring at finer spatial scales, and difficult to provide specific guidance for local land management.Additionally, much of the research used limited land cover data in specific years to analyze the impact of land cover change on long-time vegetation dynamics [4][5][6]12].This may not introduce large errors in more stable areas, but could be problematic for places that exhibit high rates of land cover change.Thus it is necessary to produce more dense land cover maps to minimize potential omission errors from bi-temporal data [3,4].
Coastal land cover across the globe has experienced remarkable change in recent decades due to extraordinary anthropogenic pressure, technological development, and climate change [14,15].Coastal reclamation is capable of creating abundant land to alleviate the pressure from land shortages.However, it significantly changes the land cover pattern with landscape fragmentation and modification to vegetation distribution and type [16][17][18].Compared to climatic factors, human activity-dominant reclamations affect vegetation growth more directly and efficiently [5].Positive human activities, such as establishment of wetland parks, can promote local vegetation restoration.Negative human activities, such as urban expansion, may exert pressure on the natural environment and contribute significantly to vegetative degradation [6,19].Assessing the location, extent, type, and frequency of land cover changes at finer spatial scale is critical for thoroughly understanding the consequences of reclamations.
Reclamation projects such as industrial and agricultural construction, aquaculture, maritime transport, and energy exploitation in the coastal zones have resulted in a sharp increase in the human demand for marine and coastal resources, and significant reduction of wetland areas.Reclamation activities provide large socioeconomic benefits and construction lands but jeopardize the coastal wetland ecosystems and services [16,20].Vegetation is one of the most important indicators of wetland ecosystem services, as it provides habitat, food, carbon storage, coastal erosion prevention, etc. [16,18,21], and reclamations greatly change the environment of plant growth.The rapid development of large-scale construction lands in coastal areas lead to the disappearance of plants, while barrier, aquaculture, and herbaceous marsh influence vegetation types and their regrowth ability [21].Quantifying the land cover changes and their impacts on vegetation growth are critical for evaluation of services and protection of coastal wetland ecosystems [22].
The coastal wetland is a complex ecosystem composed of water, sand, mud flat, and numerous vegetation species [23].Previous studies mostly focused on examining the coastal line or reclaimed area dynamics [24].The wetland changes and their impacts on wetland vegetation growth are rarely investigated, making it difficult to conduct a comprehensive assessment of the ecological changes in these strongly disturbed coastal areas.The limited data coverage from traditional field surveys makes it a challenge to monitor large-scale changes and their responses to these intense reclamation activities.Time-series satellite data provide the opportunity to analyze land cover change and vegetation greenness trends over large areas [25].Currently, analysis of trends in vegetation greenness is generally focused on very large areas [13], or at local and regional scales in forest [26], drylands [27], urban [4], and tundra areas [28], but the coastal reclaimed areas have not received much attention.
Remote sensing is one of the most important ways to monitor wetland dynamics [24], and many change detection techniques such as principal component analysis, image differencing, and post-classification comparison can be used [29][30][31].However, the wetland hydroperiod, the tidal level, and the acquisition time of remote sensing data may considerably affect wetland monitoring results.It is not easy to determine the optimal time to detect wetland change, especially for those areas influenced strongly by precipitation and tidal level [32], resulting in substantial uncertainty of wetland monitoring results and inconsistency among different studies.The dense Landsat time-series data can provide more characteristics of land cover [25], and will be helpful for wetland changes monitoring.By using all available observations and considering each pixel separately, the limitations like the failure of the Scan Line Corrector (SLC) in Landsat 7, clouds and shadows, and snow can be overcome.Another advantage is that relative normalization for each image is not needed because the time-series model already takes into account the effects of phenology and sun angle differences [25,33,34].However, most previous studies using dense Landsat data focused on detection of forest cover change (e.g., disturbance/regrowth) [33][34][35]; rarely has that type of research been used to investigate coastal areas exhibiting intense land cover change [24].
To meet the increasing demand of new lands for city expansion, agricultural purposes, industrial use, or port expansions, large areas of coastal zones in China have been reclaimed since the 1980s [36].For example, during the past 30 years, more than 400,000 ha of tidelands have been reclaimed in Zhejiang Province.It is equivalent to 4.1% of the total area of Zhejiang Province.Hangzhou Bay in Zhejiang Province is one of the coastal reclamation hotspots in China [37].Large-area reclamation promoted the local marine economy and guaranteed the construction land demand but caused a series of resource and ecological problems due to partial and complete land cover changes [37].How did these reclamations and land cover changes influence the vegetation growth?Did the changing trend of vegetation growth significantly differ from the non-reclamation area?Clarifying these questions is fundamental for government and managers to establish scientific and sustainable policies for coastal reclamation and wetland utilization.
This research attempts to answer the following questions: (1) Do the dense Landsat time-series data perform well for land cover change detection in coastal areas?(2) What major land cover changes occurred in Hangzhou Bay in the past three decades?(3) Do the coastal wetland reclamations provide sufficient areas for construction lands and alleviate the local cropland occupation?(4) How do the land cover changes influence the vegetation greenness trends in newly reclaimed area and what are the differences with the inland region?To answer these questions, dense Landsat time-series data from 1984 to 2016 were used to detect the land cover change and assess its influence on vegetation greenness in Hangzhou Bay, Zhejiang Province.

Study Area
Hangzhou Bay, an important ecological monitoring zone and world-famous trumpet-shaped tidal estuary, is located in the northeast of Zhejiang Province, China (see Figure 1).The trumpet shape of the bay leads to strong deformation of tidal waves and strengthens the tidal current.These two characteristics considerably influenced the geomorphological processes and resulted in erosion in the north side of Hangzhou Bay and sedimentation in the south side [17].Since the 1960s, artificial seawalls have built up the north side and stopped the natural changes, while the coastline on the south side has changed rapidly due to the strong sedimentation induced by the bay's geographic characteristics and coastal reclamation.The southern coastal zone is experiencing unprecedented development due to the surge of economic activities in China and is the main area of coastal reclamation in Zhejiang Province.
Although great attention has been paid to the protection of the ecological environment in Hangzhou Bay in recent years, the coastal wetland has shrunk rapidly and a large coastal area has been reclaimed due to rapid urbanization and industrialization [37].The inner baseline of the study area was determined by the coastal line in 1984 (hereafter 1984-coastline) with a buffer distance of 15 km.The outer boundary was profiled based on the 2016 Landsat Operational Land Imager (OLI) image to include the current coastal wetland such as shallow sea, mud flats, Scirpus mariqueter (hereafter Scirpus), and Spartina alterniflora (hereafter Spartina) and a 15-km strip of land along the 1984-coastline.The regions on the land side and on the bay side of the 1984-coastline were defined as inland and reclaimed area, respectively.The two parts together were called coastal area of Hangzhou Bay.

Collection of Landsat Imagery and Calculation of Vegetation Indices
Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and OLI images (Worldwide Reference System (WRS) Path 118 and Row 39) with cloud coverage of less than 80% for the years 1984 to 2016 were downloaded from the United States Geological Survey (USGS) through bulk ordering (https://espa.cr.usgs.gov).A total of 431 Landsat L1T (terrain corrected data) scenes (TM: 223 images; ETM+: 164 images; OLI: 44 images) covering the study area (Figure 2) were collected.Geo-referencing, automatic atmospheric correction [38], and detection of clouds and shadows [39,40] of the Landsat L1T products had been done by USGS.Therefore, the Landsat timeseries data could be directly used in this research for calculating Normalized Difference Vegetation Index (NDVI).For more details about the products, readers can refer to the website https://landsat.usgs.gov/landsat-data-access.
Many vegetation indices can be derived from Landsat surface reflectance data [41].The NDVI is the common index due to its robust and easily interpretable [42], and insensitivity to observation frequency in time-series analysis [43].Because of the spectral confusion between water/wetland and some dark impervious surfaces (e.g., roads, building roofs) [44], Normalized Difference Water Index (NDWI) [45] and Modified Normalized Difference Water Index (MNDWI) [46] are widely used for improving the separation between water and impervious surfaces [47].Previous research has

Collection of Landsat Imagery and Calculation of Vegetation Indices
Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and OLI images (Worldwide Reference System (WRS) Path 118 and Row 39) with cloud coverage of less than 80% for the years 1984 to 2016 were downloaded from the United States Geological Survey (USGS) through bulk ordering (https://espa.cr.usgs.gov).A total of 431 Landsat L1T (terrain corrected data) scenes (TM: 223 images; ETM+: 164 images; OLI: 44 images) covering the study area (Figure 2) were collected.Geo-referencing, automatic atmospheric correction [38], and detection of clouds and shadows [39,40] of the Landsat L1T products had been done by USGS.Therefore, the Landsat time-series data could be directly used in this research for calculating Normalized Difference Vegetation Index (NDVI).For more details about the products, readers can refer to the website https://landsat.usgs.gov/landsatdata-access.
Many vegetation indices can be derived from Landsat surface reflectance data [41].The NDVI is the common index due to its robust and easily interpretable [42], and insensitivity to observation frequency in time-series analysis [43].Because of the spectral confusion between water/wetland and some dark impervious surfaces (e.g., roads, building roofs) [44], Normalized Difference Water Index (NDWI) [45] and Modified Normalized Difference Water Index (MNDWI) [46] are widely used for improving the separation between water and impervious surfaces [47].Previous research has demonstrated that MNDWI can enhance water features and extract water bodies with higher accuracy than NDWI [17].Therefore, NDVI and MNDWI were used in this research.demonstrated that MNDWI can enhance water features and extract water bodies with higher accuracy than NDWI [17].Therefore, NDVI and MNDWI were used in this research.

Collection of Land Cover Reference Data
Field surveys for collection of different land cover types was conducted in 2015 and 2016 with GPS recording the locations (see ground photos for typical land cover types in Figure S1).The Gaofen-1 satellite imageries with spatial resolution of 2 m × 2 m (launched in 26 April 2013 by China), which were acquired on 9 January 2015 and 2 August 2015, were used to collect more land cover sample plots through visual interpretation.A total of 1950 samples were selected, of which some were randomly selected for use as training samples in the Random Forest classifier and the rest as test samples for evaluation of the classification results.All of the samples were cross-checked using the time-series NDVI (a higher NDVI indicates the higher probability of vegetation) and MNDWI (a higher MNDWI suggests the higher probability of water).Based on field investigation, China's wetland classification standard, and the research objective, a classification scheme with 10 categories was designed (see Table 1).

Collection of Land Cover Reference Data
Field surveys for collection of different land cover types was conducted in 2015 and 2016 with GPS recording the locations (see ground photos for typical land cover types in Figure S1).The Gaofen-1 satellite imageries with spatial resolution of 2 m × 2 m (launched in 26 April 2013 by China), which were acquired on 9 January 2015 and 2 August 2015, were used to collect more land cover sample plots through visual interpretation.A total of 1950 samples were selected, of which some were randomly selected for use as training samples in the Random Forest classifier and the rest as test samples for evaluation of the classification results.All of the samples were cross-checked using the time-series NDVI (a higher NDVI indicates the higher probability of vegetation) and MNDWI (a higher MNDWI suggests the higher probability of water).Based on field investigation, China's wetland classification standard, and the research objective, a classification scheme with 10 categories was designed (see Table 1).

Time Series Model and Change Detection Analysis
The land surface changes can be generally divided into three categories: (1) intra-annual change caused by vegetation phenology driven by seasonal patterns of environmental factors like temperature and precipitation; (2) gradual inter-annual change caused by climate variability, vegetation growth or gradual change in land management or degradation; and (3) abrupt change caused by land cover change due to different factors such as floods, fire, insects, and urbanization [25,48].The time-series model has components of seasonality, trend, and breaks and can capture all three categories of surface changes (Equation ( 1)) [25,33].
The initial time-series regression model can be written as follows: where the dependent variable ∧ ρ at a given time t (Julian date incorporating year) is expressed as the sum of an intercept β, a sum of different frequency harmonic components representing seasonality a and b, the slope of potentially temporal trend in the NDVI time-series data c; T is number of days per year (T = 365).After segmentation, Equation ( 1) can be re-written as: where i is the ith segment, τ * k is the kth breakpoint.The Breaks For Additive Seasonal and Trend (BFAST) algorithm [33,34,48] was used to detect abrupt changes based on NDVI time-series data.This algorithm can detect multiple breakpoints in time-series regression models using structure change detection [49].This method used a triangular residual sum of squares (RSS) matrix to test the hypothesis that the regression coefficients remain constant over time against the alternative hypothesis that at least one of the coefficients changes in the given regression model (Equation ( 1)).The optimal number of partitions is then obtained by minimizing the Bayesian information criterion, while the position of the breakpoints is determined by minimizing the RSS among all possible partitioning schemes given by the RSS matrix [50,51].
The unpartitioned time-series regression (Equation ( 1)) can therefore be re-written for each potential partitioning following Equation (2).The model coefficients were estimated by the ordinary least squares method based on the NDVI of each segment.In the time-series model, the overall value of the NDVI was captured by β.Coefficients a i and b i were used to estimate the intra-annual changes caused by phenology and sun angle differences of the NDVI in the ith segment.The inter-annual change of the NDVI in the ith segment was captured by c i .The coefficients fitted in the time-series models, mean, median and standard deviation calculated for each segment were used as the inputs for land cover classification.The Random Forest classifier was applied to each segment to obtain land cover information at any given time covered by the time-series model.
An important parameter in the multiple breakpoints detection algorithm is the minimum segment size relative to the sample size, referred to as h [33,49].The sample size refers to the number of observations in the given time series.A large h value may increase the stability of the model, but may result in short segments being missed for the frequently changed area.However, a too small value of h may result in the detection of spurious breakpoints due to the noise that might be contained in the data.Based on iteration experiments from 0.05 to 0.20 with the step 0.01, an h value of 0.10 was chosen, which had the highest accuracy of breakpoint detection.This value corresponds to approximately 25 samples in this case, as most pixels have the 60% available data of 431 images.The maximum times of land cover change in the study area was set to six based on field surveys.It should be noted that the detected breakpoint mainly started in 1987 and ended in 2015 because the h value determined a minimum number of clear observations through the moving window.

Continuous Land Cover Classification
After the multiple breakpoints detection process, each pixel had its own time-series models before and after any detected breakpoints.By classifying each segment based on the time-series model coefficients, the land cover type for the entire period for each time-series model could be acquired.
The main basis of the land cover classification was that different land cover classes have distinctive shapes for the estimated time-series models.Figure 3 illustrates the estimated time-series models for the 10 land cover types in the study area.The sea and mud tidal flat did not have obvious seasonal change, and water and construction land had slight seasonal change due to the mixed pixels containing vegetation.All the vegetation-related land covers showed obvious seasonal change trends.Scirpus showed the least amplitude of variation, and paddy and herb had the highest amplitude of variation.The peak time of NDVI for herb was slightly earlier than Spartina.Forest had the highest NDVI, followed by dryland and paddy field.The different characteristics among these land covers imply that the information contained in the NDVI time-series model is helpful for land cover classification.

Continuous Land Cover Classification
After the multiple breakpoints detection process, each pixel had its own time-series models before and after any detected breakpoints.By classifying each segment based on the time-series model coefficients, the land cover type for the entire period for each time-series model could be acquired.
The main basis of the land cover classification was that different land cover classes have distinctive shapes for the estimated time-series models.Figure 3 illustrates the estimated time-series models for the 10 land cover types in the study area.The sea and mud tidal flat did not have obvious seasonal change, and water and construction land had slight seasonal change due to the mixed pixels containing vegetation.All the vegetation-related land covers showed obvious seasonal change trends.Scirpus showed the least amplitude of variation, and paddy and herb had the highest amplitude of variation.The peak time of NDVI for herb was slightly earlier than Spartina.Forest had the highest NDVI, followed by dryland and paddy field.The different characteristics among these land covers imply that the information contained in the NDVI time-series model is helpful for land cover classification.The four variables (β, a, b, and c) generated from time-series model estimation were used as inputs for land cover classification.Figure 4 shows the mean values of these four variables for each land cover category based on land cover reference data.It is clear that different land covers show distinct shapes in the plots of the four variables.The mean, median, and standard deviation values of NDVI for each segment were also calculated.Thus, seven variables were used as the input for classification.The Random Forest classifier was used to perform land cover classification because of its relatively high accuracy and computational efficiency [25,52].The four variables (β, a, b, and c) generated from time-series model estimation were used as inputs for land cover classification.Figure 4 shows the mean values of these four variables for each land cover category based on land cover reference data.It is clear that different land covers show distinct shapes in the plots of the four variables.The mean, median, and standard deviation values of NDVI for each segment were also calculated.Thus, seven variables were used as the input for classification.The Random Forest classifier was used to perform land cover classification because of its relatively high accuracy and computational efficiency [25,52].

Continuous Land Cover Classification
After the multiple breakpoints detection process, each pixel had its own time-series models before and after any detected breakpoints.By classifying each segment based on the time-series model coefficients, the land cover type for the entire period for each time-series model could be acquired.
The main basis of the land cover classification was that different land cover classes have distinctive shapes for the estimated time-series models.Figure 3 illustrates the estimated time-series models for the 10 land cover types in the study area.The sea and mud tidal flat did not have obvious seasonal change, and water and construction land had slight seasonal change due to the mixed pixels containing vegetation.All the vegetation-related land covers showed obvious seasonal change trends.Scirpus showed the least amplitude of variation, and paddy and herb had the highest amplitude of variation.The peak time of NDVI for herb was slightly earlier than Spartina.Forest had the highest NDVI, followed by dryland and paddy field.The different characteristics among these land covers imply that the information contained in the NDVI time-series model is helpful for land cover classification.The four variables (β, a, b, and c) generated from time-series model estimation were used as inputs for land cover classification.Figure 4 shows the mean values of these four variables for each land cover category based on land cover reference data.It is clear that different land covers show distinct shapes in the plots of the four variables.The mean, median, and standard deviation values of NDVI for each segment were also calculated.Thus, seven variables were used as the input for classification.The Random Forest classifier was used to perform land cover classification because of its relatively high accuracy and computational efficiency [25,52].We established a set of rules to improve the classification results.First, the change points were deleted if the classes of the two adjacent segments were the same.This step can exclude the noises of NDVI time series, even the breakpoint induced by the NDVI outlier.Second, the mean MNDWI value of each segment was calculated and the threshold method was used to check water and construction land classes (the threshold was determined by the mean of the referenced data).If a pixel was classified as construction land and the MNDWI was larger than 0 in the segment, the pixel was reclassified to water.If a pixel was classified as water and the MNDWI was less than −0.1, the pixel was reclassified to construction land.This step reduced the classification error between water and construction land.

Trends Analysis of NDVI Time-Series Data
NDVI is an indicator of chlorophyll abundance and correlates to vegetation growth condition and photosynthetic capacity [53].The NDVI changes can be divided into three categories (seasonal, gradual, and abrupt changes) through the time-series models (Section 3.2).Seasonal change occurs when the land surface phenology changes, e.g., earlier onset of greening in spring driven by temperature or rainfall.This type of change will not affect the NDVI trend component.The NDVI trend analysis focused on gradual changes and abrupt changes.Gradual change reflects the slowly acting environmental processes such as climate change, land management practices, and plant life cycle on vegetation growth or degradation.Abrupt change shows the impacts of fast-acting processes such as land use/cover change, wildfire, typhoon, and flood events on vegetation growth and coverage.By summing the gradual and abrupt changes for each pixel, the total change indicates the net greening or browning during the long period.
For each pixel, the time-series model of each segment was estimated and the start and end times of each model was also recorded.Equation (3) defined the overall NDVI value for each pixel at the start and end of each segment.Equations ( 4)-( 6) defined the accumulated gradual change (G), the accumulated abrupt change (A), and the total change in greenness (Tot).Annual mean NDVI change magnitude of gradual and abrupt changes for each year for all pixels in the study area was also calculated. NDV where, a i and c i : coefficients for overall value and for inter-annual change (slope) for the ith segment; t start,i and t end,i : start (Julian date) and end (Julian date) of the ith segment; NDVI start,i and NDVI end,i : overall NDVI values at the start and end of the ith segment; K: total number of segments estimated for a pixel.

Impacts of Land Cover Change on Vegetation Greenness Dynamics
Positive and negative changes of NDVI through time can be referred to as greening and browning, respectively [12].If there is little or no land cover change in the study area, the greenness trends can be easily quantified through the slope coefficient based on a simple linear regression of the NDVI.However, for areas characterized by frequent land cover changes or trend changes, this approach may provide misleading or incomplete results [4,12].Therefore, this study separated the abrupt change from gradual change.Using the land cover information for each pixel at any given time, the detected gradual and abrupt NDVI changes per land cover class can be summarized.Note that the summarized abrupt NDVI change was based on the posterior date class to account for the effects of land cover change on the greenness trends.All described preprocesses and analyses were performed using R statistical software (R Development Core Team, 2011) with the packages "raster", "bfast", "bfastSpatial", "strucchange", "randomForest", etc.

Assessment of Breakpoints Detection Results
A stratified random sample design was used for assessing the change detection accuracy.A total of 1000 sample plots were selected, of which 400 were from unchanged areas throughout the time of the analysis (1984-2016) and 600 were selected from the changed area where they were determined by the multiple breakpoints detection algorithm.High and medium spatial resolution satellite data were often used for collection of reference data [25,54,55].In this research, the reference data for evaluation of change times and change years were manually interpreted by the expert using time-series Landsat data, assisted by high spatial resolution imagery in Google Earth and Tiandi Maps (http://www.tianditu.cn/).This method to collect validation data is consistent with approaches recommended or used in previous Landsat time-series studies [25,54].By carefully examining the time-series data, the expert can identify whether the pixels were changed and when the change occurred.The acquisition date of the image in which the breakpoint was detected was used as a surrogate date when the land cover change occurred.If there were multiple changes within one pixel, the number of change times and the corresponding years were recorded.Only the pixels with high interpretation confidence were selected.
The normalized root mean square error (NRMSE), normalized mean absolute error (NMAE), and bias were used to assess the accuracy of multiple breakpoints detection results.NRMSE and NMAE reflect the relative amount of residual variance.Bias indicates whether the model predictions over-or under-estimated.The NRMSE, NMAE, and bias can be computed using the following equations: where y max − y min is the range of the observation of the observed variables, ŷi is the corresponding predicted value for that observation, n is the total number of observations used for validation, and y i is the ith observation of the observed variables.

Assessment of Land Cover Classification
The land cover reference data previously introduced for continuous land cover classification were used as the basis for assessing the accuracy of classification.As the land cover reference data were collected based on the Gaofen-1 data acquired in 2015, this study only assessed the accuracy of land cover classification at this specific time interval (9 January 2015-2 August 2015).The k-fold cross-validation method was used for the model training and test with the k value set as 10.This study selected the best model for the final classification.In this approach, each reference pixel was used only once for training and accuracy assessment.The user's accuracy, producer's accuracy, and overall classification accuracy were used to assess the land cover classification result.

Accuracy Assessment of Change Detection and Classification
The accuracy assessment for change detection results (Table 2) showed that the land cover change was successfully detected with an overall accuracy of 88.7%.We further evaluated the successfully detected 467 changed samples and obtained the NRMSE, NMAE, and bias values of 0.27, 0.13, and −0.08, respectively, indicating that the multiple breakpoint detection algorithm provided reliable information in the times of land cover change detection in the coastal area without strong systematic errors.For these 467 changed pixels, 713 breakpoints were interpreted by the expert analysts, but 677 breakpoints were detected by the multiple breakpoints detection algorithm.According to the 677 breakpoints detected by the algorithm, the proportion of the breakpoints in the same year of change for the algorithm results and the interpreted results was 86.1% (583 breakpoints).The accuracy assessment result for wetland classification (Table 3) indicated that the overall accuracy for ten land cover classes reached 91.0%.The user's accuracies were over 85.9% and producer's accuracies were over 80.4% for all land cover classes.The major confusions were between water and Scirpus, and between construction land and dryland.

Intensity and Types of Land Cover Change
The area statistics of different change times between 1987 and 2015 (see Table 4) indicated that a large proportion of the study area has changed, accounting for 39%.Of these changed areas, the once and twice changed areas accounted for 73.5% and 22.6% of the total changed area.The land cover change had obvious spatial patterns (see Figure 5), for example, the unchanged areas were mainly located inland and the changed areas were distributed in the reclaimed region.Such patterns are expected because most of the reclaimed areas were experienced the change from sea/mud flat to vegetation coverage or aquaculture.In addition, due to the rapid economic development, large amounts of area close to villages and roads were experienced obvious land cover change from dryland to construction land.The land covers of Hangzhou Bay have experienced intense and complex change in the past 32 years (Table 5).The prominent changes were the changes from dryland to construction land, from water to herb, and from forest to dryland.The largest lost and gained areas occurred in water and construction land, respectively.The areas of Scirpus, Spartina, herb, and construction land showed net increases, and the areas of forest, mud flat, dryland, water, and sea showed net decreases.To simplify the analysis of wetland changes, this study collapsed the initial categories of wetland-related land covers into three classes: CW (coastal wetland, including sea, mud, Scirpus, Spartina), SW (swamp wetland, including herb), and AW (artificial wetland, including water and paddy field).The largest wetland-related change was from CW to AW, followed by the change from AW to SW (Figure 6).The largest loss of wetland was due to the increase of construction land and dryland.    1 The numbers in parentheses are the percentages of changed areas at different change times to total changed area.
The land covers of Hangzhou Bay have experienced intense and complex change in the past 32 years (Table 5).The prominent changes were the changes from dryland to construction land, from water to herb, and from forest to dryland.The largest lost and gained areas occurred in water and construction land, respectively.The areas of Scirpus, Spartina, herb, and construction land showed net increases, and the areas of forest, mud flat, dryland, water, and sea showed net decreases.To simplify the analysis of wetland changes, this study collapsed the initial categories of wetland-related land covers into three classes: CW (coastal wetland, including sea, mud, Scirpus, Spartina), SW (swamp wetland, including herb), and AW (artificial wetland, including water and paddy field).The largest wetland-related change was from CW to AW, followed by the change from AW to SW (Figure 6).The largest loss of wetland was due to the increase of construction land and dryland.    1 The numbers in parentheses are the percentages of changed areas at different change times to total changed area.
The land covers of Hangzhou Bay have experienced intense and complex change in the past 32 years (Table 5).The prominent changes were the changes from dryland to construction land, from water to herb, and from forest to dryland.The largest lost and gained areas occurred in water and construction land, respectively.The areas of Scirpus, Spartina, herb, and construction land showed net increases, and the areas of forest, mud flat, dryland, water, and sea showed net decreases.To simplify the analysis of wetland changes, this study collapsed the initial categories of wetland-related land covers into three classes: CW (coastal wetland, including sea, mud, Scirpus, Spartina), SW (swamp wetland, including herb), and AW (artificial wetland, including water and paddy field).The largest wetland-related change was from CW to AW, followed by the change from AW to SW (Figure 6).The largest loss of wetland was due to the increase of construction land and dryland.The land cover changes of the Hangzhou Bay showed obvious spatial variations.As most pixels experienced land cover changes once or twice, here we just showed the change directions for the earliest and the latest land cover changes (Figure 7a,b).For the inland region, land cover was mainly changed to construction land and dryland, and was located near villages, newly built roads, and forests.In the reclaimed region, the land cover was mainly changed to herb, water, and construction land.In the earliest land cover change, the middle of the reclaimed area was mainly changed to water for aquaculture, and the western and eastern parts mostly changed to herb.In the latest land cover change, large areas in the western part were changed to herb due to the establishment of Hangzhou Bay National Wetland Park, and many places in the middle and eastern parts were changed to construction land due to the establishment of Hangzhou Bay New District.The land cover changes of the Hangzhou Bay showed obvious spatial variations.As most pixels experienced land cover changes once or twice, here we just showed the change directions for the earliest and the latest land cover changes (Figure 7a,b).For the inland region, land cover was mainly changed to construction land and dryland, and was located near villages, newly built roads, and forests.In the reclaimed region, the land cover was mainly changed to herb, water, and construction land.In the earliest land cover change, the middle of the reclaimed area was mainly changed to water for aquaculture, and the western and eastern parts mostly changed to herb.In the latest land cover change, large areas in the western part were changed to herb due to the establishment of Hangzhou Bay National Wetland Park, and many places in the middle and eastern parts were changed to construction land due to the establishment of Hangzhou Bay New District.The years of the land cover changes also showed spatial variations.The spatial distributions of the earliest and latest years of land cover changes (Figure 7c,d) indicated that the changes in the inland area occurred mainly before 2005, and in the reclaimed area after 2005.For the earliest land cover changes, the western part of the reclaimed area occurred mainly before 2000, and changes in the middle and eastern areas occurred mostly after 2000 (Figure 7c).For the latest land cover changes, major changes occurred after 2005, especially in 2013 in the west and middle parts of the study area.The difference between the earliest and the latest change years ranged from 1 year (excluding areas that changed only once during the period) to 28 years.
The annual change area of each land cover category showed distinct trends during 1987-2015.The temporal dynamics of each land cover category (Figure 8) based on the lost (changed to other categories), gained (changed from other categories), and net changed area indicated that the major lost land cover categories were water, dryland, and forest, which peaked at 2013, 2001, and 2001, respectively (Figure 8a).The lost areas of dryland and forest showed increase trends from 1997 and peaked in 2001.Mud flat showed obvious decrease starting in 2005.The most obviously gained category was construction land, which occurred mainly after 1999 (Figure 8b).Dryland area also showed considerable increase at the same time due to the compensation mechanism for occupied dryland.Therefore, the net changed area of dryland was little during the study period (Table 5 and Figure 9c).Construction land was the only category that kept a consistent net increase during the period of 1987-2015.
Remote Sens. 2018, 10, 32 13 of 22 The years of the land cover changes also showed spatial variations.The spatial distributions of the earliest and latest years of land cover changes (Figure 7c,d) indicated that the changes in the inland area occurred mainly before 2005, and in the reclaimed area after 2005.For the earliest land cover changes, the western part of the reclaimed area occurred mainly before 2000, and changes in the middle and eastern areas occurred mostly after 2000 (Figure 7c).For the latest land cover changes, major changes occurred after 2005, especially in 2013 in the west and middle parts of the study area.The difference between the earliest and the latest change years ranged from 1 year (excluding areas that changed only once during the period) to 28 years.
The annual change area of each land cover category showed distinct trends during 1987-2015.The temporal dynamics of each land cover category (Figure 8) based on the lost (changed to other categories), gained (changed from other categories), and net changed area indicated that the major lost land cover categories were water, dryland, and forest, which peaked at 2013, 2001, and 2001, respectively (Figure 8a).The lost areas of dryland and forest showed increase trends from 1997 and peaked in 2001.Mud flat showed obvious decrease starting in 2005.The most obviously gained category was construction land, which occurred mainly after 1999 (Figure 8b).Dryland area also showed considerable increase at the same time due to the compensation mechanism for occupied dryland.Therefore, the net changed area of dryland was little during the study period (Table 5 and Figure 9c).Construction land was the only category that kept a consistent net increase during the period of 1987-2015.

Analysis of Long-Term NDVI Trends and Their Relationship to Land Cover Change
The NDVI change trends showed distinct spatial variations between inland and reclaimed areas (Table 6 and Figure 9).The inland region showed mainly a browning trend, and reclaimed area showed a considerable greening trend according to NDVI values.The means of the total NDVI changes were −0.052 in inland and 0.209 in reclaimed area.Similar to total NDVI change, the mean of gradual NDVI change was negative, but positive in the reclaimed region.Overall, more than half of Hangzhou Bay showed a greening trend according to the gradual NDVI.The pixels with abrupt NDVI change showed that only 11.4% appeared as positive trends in the inland region, indicating most land cover changes resulted in decreased greenness, especially for nearby villages (Figure 9c).For the reclaimed region, most of the land cover changes showed an increased NDVI following abrupt changes.Generally, abrupt changes played a negative effect on the total greenness in Hangzhou Bay due to the larger contribution from the inland region.The annual mean gradual NDVI changes were all positive and the magnitude marginally increased after 2001 (Figure 10).The area that showed positive gradual NDVI changes varied between 53% and 63% of the total area during the study period.The annual mean abrupt NDVI changes were generally negative and varied in time.The magnitude was higher in 1991 and 2003 than in other years.The mean annual abrupt NDVI changes quickly increased after 2003.The percentage of positive area showed increase after 1992 and reached maximum around 2008.

Analysis of Long-Term NDVI Trends and Their Relationship to Land Cover Change
The NDVI change trends showed distinct spatial variations between inland and reclaimed areas (Table 6 and Figure 9).The inland region showed mainly a browning trend, and reclaimed area showed a considerable greening trend according to NDVI values.The means of the total NDVI changes were −0.052 in inland and 0.209 in reclaimed area.Similar to total NDVI change, the mean of gradual NDVI change was negative, but positive in the reclaimed region.Overall, more than half of Hangzhou Bay showed a greening trend according to the gradual NDVI.The pixels with abrupt NDVI change showed that only 11.4% appeared as positive trends in the inland region, indicating most land cover changes resulted in decreased greenness, especially for nearby villages (Figure 9c).For the reclaimed region, most of the land cover changes showed an increased NDVI following abrupt changes.Generally, abrupt changes played a negative effect on the total greenness in Hangzhou Bay due to the larger contribution from the inland region.The annual mean gradual NDVI changes were all positive and the magnitude marginally increased after 2001 (Figure 10).The area that showed positive gradual NDVI changes varied between 53% and 63% of the total area during the study period.The annual mean abrupt NDVI changes were generally negative and varied in time.The magnitude was higher in 1991 and 2003 than in other years.The mean annual abrupt NDVI changes quickly increased after 2003.The percentage of positive area showed increase after 1992 and reached maximum around 2008.Gradual NDVI changes indicated the vegetation growth or degradation process.The impacts of different land covers on the gradual change trends were different (Figure 11).All land cover categories except sea and dryland showed greenness trends according to the mean gradual NDVI changes, with herb showing the largest magnitude, following Scirpus and Spartina (Figure 11a).The mean gradual NDVI change in mud was the lowest in magnitude.The water and construction land also showed increased gradual NDVI changes.The reason might be the increased growth of phytoplankton from eutrophication of water, mixed pixel with vegetation, and the impact of NDVIs from different remote sensors (discussed in Section 4.1).Due to the large magnitude and high frequency of herb and water (Figure 11b), they were the main contributors to the overall greenness trend in terms of gradual change.Gradual NDVI changes indicated the vegetation growth or degradation process.The impacts of different land covers on the gradual change trends were different (Figure 11).All land cover categories except sea and dryland showed greenness trends according to the mean gradual NDVI changes, with herb showing the largest magnitude, following Scirpus and Spartina (Figure 11a).The mean gradual NDVI change in mud was the lowest in magnitude.The water and construction land also showed increased gradual NDVI changes.The reason might be the increased growth of phytoplankton from eutrophication of water, mixed pixel with vegetation, and the impact of NDVIs from different remote sensors (discussed in Section 4.1).Due to the large magnitude and high frequency of herb and water (Figure 11b), they were the main contributors to the overall greenness trend in terms of gradual change.Gradual NDVI changes indicated the vegetation growth or degradation process.The impacts of different land covers on the gradual change trends were different (Figure 11).All land cover categories except sea and dryland showed greenness trends according to the mean gradual NDVI changes, with herb showing the largest magnitude, following Scirpus and Spartina (Figure 11a).The mean gradual NDVI change in mud was the lowest in magnitude.The water and construction land also showed increased gradual NDVI changes.The reason might be the increased growth of phytoplankton from eutrophication of water, mixed pixel with vegetation, and the impact of NDVIs from different remote sensors (discussed in Section 4.1).Due to the large magnitude and high frequency of herb and water (Figure 11b), they were the main contributors to the overall greenness trend in terms of gradual change.Abrupt NDVI changes reflected the impacts of land cover changes on vegetation regrowth.When one land cover was changed to Spartina, herb, forest, and paddy field, positive mean abrupt NDVI changes were observed.On the other hand, if one land cover was changed to such classes as sea, water, mud flat, Scirpus, and construction land, negative mean abrupt NDVI changes were observed.The largest negative change magnitude was associated with areas where land covers were changed to construction land, mainly at the cost of dryland.The construction land and water were the most important contributors to the overall negative abrupt NDVI change trend, as many abrupt changes ended up in these classes (Figure 11b) and with a high magnitude of mean NDVI (Figure 11a).

Uncertainties and Limitations of This Study
This research used the BFAST algorithm to successfully extract land cover change and greenness dynamics in the coastal wetland ecosystem based on dense Landsat time series data.The accuracies are compatible with other studies (e.g., [25,33,34,56]).However, there are some uncertainties and limitations that need to be noted.
The temporal density of the constructed NDVI time-series data is a critical factor determining the probability of breakpoint detection.Although NDVI is one of the most robust indexes when observation frequency varies [43], the limitation in data availability (such as in the 1980s and 1990s) is likely to result in missing breakpoints for detecting fast dynamics of some land cover types.Therefore, land cover change detection based on Landsat images before 2000 might be underestimated.The minimum segment size (h) is another important temporal density-relevant factor.This study set h as 0.10 through a preliminary test based on the validation data.As h is defined as a minimum observation through the moving window, there should be at least 25 (most pixels in the 431 images had about 60% available data) clear observations for the beginning (first segment) and ending (last segment) breakpoints, and two consecutive breakpoints.Generally, continuous changes in a short duration seldom occurred in the inland area, because the main change was the conversion from dryland to construction land.On the other hand, some breakpoints in the reclaimed area where frequent land cover change existed might be missed.Both data density and h value mostly induced omission errors (with the bias −0.08 in Hangzhou Bay), while commission errors were partially overcome by intelligent post-processing (Section 3.3).The omission errors can be improved if the frequency of images increases by incorporation of other similar sensors such as HuanJing satellite (HJ-1A/B CCD), Advanced Land Observation Satellite (ALOS), Sentinel-2, or data fusion with high temporal resolution data (e.g., MODIS).
Different sensors (TM sensor on Landsat 4 and 5, ETM+ sensor on Landsat 7, and OLI sensor on Landsat 8) might influence the greenness trend analysis due to reflectance inconsistencies introduced by variable sun-surface-sensor geometry, sensor degradation, atmosphere, and sensor calibration [41,57].Previous studies showed that TM and ETM+ sensors were calibrated well with each other [58], but the consistency of data between ETM+ and OLI is still under further research.A previous study found that there were subtle differences between two sensors in NDVI, and had average difference values less than ±0.05 and a smaller standard deviation [59].Additionally, correlation analysis of the vegetation indices indicated that both sensors had a very high linear correlation coefficient, with R 2 greater than 0.96 [59].The subtle differences and high correlation of vegetation indices demonstrated that ETM+ and OLI imagery could be used as complementary data.However, it has been reported that the NDVI differences between the two sensors can be as large as 5% [57], and the NDVI values from OLI were larger than those from TM/ETM+ sensors [4], this might lead to a positive gradual NDVI change trend.However, this will not change the absolute different NDVI change trends between inland area and reclaimed area.Nevertheless, further studies are needed to quantify the impacts of the different sensors on NDVI trend analysis.
We assumed that the simple sinusoidal model and NDVI would be able to estimate all kinds of land cover types and breakpoints.A more complex time series model is needed for these highly seasonal land covers, such as agriculture with multiple rotations.Choosing more vegetation indices such as NBR (Normalized Burn Ratio) and NDMI (Normalized difference Moisture Index) may benefit the land cover change detection accuracy [43].Developing change detection approaches that take advantage of all the spectral information available should therefore constitute a priority for future work.In addition, vegetation indexes or bands that are less biased between different sensors also should be noticed.

Drivers of Land Cover Change Dynamics in Hangzhou Bay
Following the abrupt change detection and posterior and prior classification of breakpoints, this study produced the history of land cover and land cover change at Hangzhou Bay.It captured more details of the land cover change than using limited land cover maps in given years.The area experienced land cover change amounts to 39% across the study area, which is a little higher than Guangzhou City, China [4].The land cover change had distinct characteristics between inland and the reclaimed area.For the inland area, urbanization was the main driver of the land cover change, with most occurring before 2005.Inspired by the fast-growing economy and urban population growth, large amounts of area close to towns and roads were changed from dryland to construction land.Such a land cover change pattern was common around the world [1,[60][61][62][63], especially in these rapid economic development regions of China [4,64].Another dominant inland land cover change was from forest to dryland, which was driven by the strict cultivated land protection policies of China, because the acceleration of urbanization at the expense of arable land threated food supply and demands [64].
In the reclaimed region, the land cover was mainly changed to herb, water, and construction land.Aquaculture (shrimp farms, crab farms, and fish ponds) was the main purpose of reclamation in most of the Hangzhou Bay area [65], because the cost of this kind of reclamation was relatively small, and humans can derive substantial value from aquaculture outputs [14].Another important reason was the higher levels of soil salinity, which limited the land as a resource for agriculture during the early stage of the reclaimed land [66].During 1984-2016, the sea and mud flat areas were mainly changed to water for aquaculture in Hangzhou Bay.The water area showed an obvious increase trend in 2005, mainly from the conversion of mud flats.This conformed to the reality that more reclamation projects started in 2005 [64].In the latest land cover change, large areas in the middle and eastern Hangzhou Bay changed to construction land accompanied by the establishment of Hangzhou Bay New District (planning an area of 145 km 2 ) in 2001.This land cover change pattern was common in coastal areas of China [14].Before the 1990s, coastal reclamation was aimed at increasing salt, agricultural, and fishery production in China.In recent decades, rapid expansion of the coastal economy and accelerated coastal population growth have caused a sharp increase in industry, urban expansion, and infrastructure.Coastal reclamation has resulted in the cumulative loss of coastal wetlands in Hangzhou Bay.The largest wetland-related change was from CW to AW, followed by AW to SW.The largest loss of wetland was the conversion to construction land and dryland.This undoubtedly severely affected the wetland ecosystem function and structure.Aquaculture creates artificial wetland and retains some functions of a wetland ecosystem.However, industrial building and port construction are no longer wetlands, and cause permanent damage to them.Fortunately, the government has recognized the importance of ecological function in coastal wetland.With the support of the Global Environment Facility and World Bank, the Zhejiang provincial government has established Hangzhou Bay National Wetland Park (43.5 km 2 ) in the study area.This was the main driver of the increased herb area and increased vegetation coverage in the western part of Hangzhou Bay.

Impacts of Land Cover Change on Vegetation Greenness
Land cover change has significant impacts on vegetation variation [5].It can change the vegetation coverage directly and bring subsequent effects on revegetation and regrowth.The negative effects of land cover change such as urbanization has proved to degrade vegetation cover [4], whereas positive land cover change such as returning farmland to forest and grassland lead to sustained vegetation growth [5].However, previous studies mostly focused on the impact of land cover change in long NDVI change trends at national and global levels or on specific, large ecological restoration programs [5,12,13].This study used abrupt and gradual NDVI changes to reflect the direct impacts of land cover change and subsequent effects on vegetation greenness, respectively.The abrupt NDVI change can be directly related to human activities or natural processes.The gradual NDVI change reflected coupled effects of land cover, regrowth process, climate change, rising atmospheric CO 2 concentration, and nitrogen deposition.The results showed that NDVI changes revealed distinct spatial variations between inland and reclaimed areas.The abrupt and gradual changes decreased mean NDVI in the inland area, but increased mean NDVI in the reclaimed area.The main reason was the different land cover change patterns between the two areas.For the inland area, the dominant land cover changes were dryland to construction land and forest to dryland, both of them inducing decreased greenness.This negative effect of land cover changes on vegetation greenness has been widely observed [4,13].
Most the reclaimed area in Hangzhou Bay showed increased greenness based on the analysis of gradual and abrupt NDVI.Herb and water were the main contributors to the overall increase of NDVI in terms of gradual change.The water was mainly changed from sea, mud flat, and herb in the reclaimed region.When mud flat and sea were changed to water, the greenness increased as water is a suitable (steady) environment for vegetation growth around water ponds.Herb had higher NDVI and showed the largest gradual NDVI increase trend once they had a stable living environment.Herb area has been significantly restored with the establishment of Hangzhou Bay National Wetland Park and increased strict protection policies of coastal wetland in recent years, and will play an important role in improving the vegetation coverage of the coastal area [14].Spartina was also an important driver of both increased gradual NDVI and abrupt NDVI.With its great capacity for reducing tidal wave energy, mitigating erosion, and trapping sediments, Spartina has been widely introduced to many coastal and estuarine regions of the world as a species for ecological engineering [67].This exotic plant was introduced to Hangzhou Bay in the 1990s and has been expanding rapidly.However, evidence has recently been presented that Spartina may out-compete native plants, threaten the native ecosystems and coastal aquaculture, and cause declines in native species richness [67].
One previous study indicated that wetland reclamation could cause nearby wetlands to become vulnerable because of the changes in wetland habitat environment and structure [2].Reclamation can cause the quantity, density, and biomass of vegetation to decrease in the tidal flat [37].According to our study, both herb and water had strong regrowth ability and were the main contributors to the overall increase of NDVI.Climate change, eutrophication, and the protection of coastal wetlands might also promote vegetation greenness, as most land cover showed increased gradual NDVI in the study area.A highly intense land use change contributed to a pronounced vegetation greenness decline, while low intense land use change might result in vegetation greenness increase.This was also suggested in a previous study which focused on the species diversity and biomass of macrozoobenthos [68].However, the vegetation sprawl over the reclaimed area should be considered as a potential threat to loss of mudflat habitats.Meanwhile, the influences of the reestablishment of marsh vegetation on other plant lives are less known.Increased vegetation greenness may potentially restore the ecosystem's food web, but at very high densities it may also obstruct waterbird foraging and movement [69].Hence the benefits of the increased greenness from gradual NDVI change and land cover changes remain open to further study [70].

Conclusions
This study integrated dense Landsat time-series data, BFAST algorithm, and the Random Forest classification method to study the land cover change and vegetation greenness variations of Hangzhou Bay during 1987-2015.The overall accuracies of 88.7% for land cover change and 91.0% for ten land cover classifications were obtained.The area experienced land cover change amounts to 39% of the study area, and showed spatial-temporal dynamics.The intense changes mainly occurred in 1995, 2001, and 2013, coincident with the local policy and the establishment of Hangzhou Bay New Zone and Hangzhou Bay National Wetland Park.For the inland area, urbanization was the main driver of the land cover change and mostly occurred before 2005.The land cover change reduced mean vegetation greenness because the conversions from construction land to dryland and from forest to dryland were the dominant land cover changes.In the reclaimed area, the land cover was mainly changed to herb, water, and construction land.Land cover change has led to the increase of the mean vegetation greenness of the reclaimed area, benefitting from the increased area of herb, water, and Spartina.In summary, construction land and water were the most important contributors to the overall negative abrupt NDVI change trend.Herb and water were the main contributors to the overall increase of NDVI in terms of gradual change.
Using all available Landsat data can provide detailed information about land cover change and vegetation growth in the frequently changed coastal area.This method can take full advantage of vegetation phenology information, avoid the hydroperiod effects, and quantify the abrupt and gradual NDVI changes.Such information of land cover change and impacts on coastal vegetation greenness will help the research and natural resource management communities to generate better understanding of coastal ecosystems and make better management decisions.However, further studies are still ongoing to thoroughly understand the impacts of reclamation and coastal wetland reduction on environmental functions, including biodiversity, rising sea level, and carbon cycle.

Figure 1 .
Figure 1.Location of the study area: Hangzhou Bay, Zhejiang Province, China (The red dash line is the boundary of the study area.The white dash line is the 1984 coastline defined from visual interpretation of the 1984 Landsat imagery.The regions on the land side and on the bay side of the 1984 coastline are defined as inland and reclaimed area, respectively.The color composite using the shortwave infrared, near-infrared, and red spectral bands is from a 2016 Landsat 8 Operational Land Imager image).

Figure 1 .
Figure 1.Location of the study area: Hangzhou Bay, Zhejiang Province, China (The red dash line is the boundary of the study area.The white dash line is the 1984 coastline defined from visual interpretation of the 1984 Landsat imagery.The regions on the land side and on the bay side of the 1984 coastline are defined as inland and reclaimed area, respectively.The color composite using the shortwave infrared, near-infrared, and red spectral bands is from a 2016 Landsat 8 Operational Land Imager image).

Figure 2 .
Figure 2. The number of Landsat scenes between 1984 and 2016 used in this research.

Figure 2 .
Figure 2. The number of Landsat scenes between 1984 and 2016 used in this research.

Figure 3 .
Figure 3. Fitted lines of each land cover type based on the parameters from time-series models (note: Dryl., dryland; Cons., construction land).

Figure 4 .
Figure 4.The values of the four parameters (β intercept, c trend, a harmonic cos, b harmonic sin) in the time-series model for different land cover classes based on the reference samples for each land cover category (note: the trend value c was enlarged by 10 4 ; Dryl., dryland; Cons., construction land).

Figure 3 .
Figure 3. Fitted lines of each land cover type based on the parameters from time-series models (note: Dryl., dryland; Cons., construction land).

Figure 3 .
Figure 3. Fitted lines of each land cover type based on the parameters from time-series models (note: Dryl., dryland; Cons., construction land).

Figure 4 .
Figure 4.The values of the four parameters (β intercept, c trend, a harmonic cos, b harmonic sin) in the time-series model for different land cover classes based on the reference samples for each land cover category (note: the trend value c was enlarged by 10 4 ; Dryl., dryland; Cons., construction land).

Figure 4 .
Figure 4.The values of the four parameters (β intercept, c trend, a harmonic cos, b harmonic sin) in the time-series model for different land cover classes based on the reference samples for each land cover category (note: the trend value c was enlarged by 10 4 ; Dryl., dryland; Cons., construction land).
amounts of area close to villages and roads were experienced obvious land cover change from dryland to construction land.

Figure 5 .
Figure 5. Times of land cover change on the south side of Hangzhou Bay in Zhejiang Province, China (the black dash line represents the coastline in 1984).

Figure 5 .
Figure 5. Times of land cover change on the south side of Hangzhou Bay in Zhejiang Province, China (the black dash line represents the coastline in 1984).
Remote Sens. 2018, 10, 32 11 of 22 amounts of area close to villages and roads were experienced obvious land cover change from dryland to construction land.

Figure 5 .
Figure 5. Times of land cover change on the south side of Hangzhou Bay in Zhejiang Province, China (the black dash line represents the coastline in 1984).

Figure 7 .
Figure 7. Land cover change directions and corresponding change years at Hangzhou Bay in Zhejiang Province, China, during 1987-2015.(a,b) are the earliest and latest land cover change directions, respectively; (c,d) are the years of the earliest and latest land cover changes, respectively.If a pixel changed only once, (a) is the same as (b), and (c) is the same as (d) (note: Dryl., dryland; Cons., construction land).

Figure 7 .
Figure 7. Land cover change directions and corresponding change years at Hangzhou Bay in Zhejiang Province, China, during 1987-2015.(a,b) are the earliest and latest land cover change directions, respectively; (c,d) are the years of the earliest and latest land cover changes, respectively.If a pixel changed only once, (a) is the same as (b), and (c) is the same as (d) (note: Dryl., dryland; Cons., construction land).

Figure 8 .
Figure 8. Annual change area of each land cover at Hangzhou Bay in Zhejiang Province, China, during 1987-2015.(a) lost area (area changed to other categories); (b) gained area (area changed from other categories); and (c) net changed area (note: Dryl., dryland; Cons., construction land).

Figure 8 .
Figure 8. Annual change area of each land cover at Hangzhou Bay in Zhejiang Province, China, during 1987-2015.(a) lost area (area changed to other categories); (b) gained area (area changed from other categories); and (c) net changed area (note: Dryl., dryland; Cons., construction land).

Figure 10 .
Figure 10.Annual mean and percentage of positive gradual and abrupt NDVI changes at Hangzhou Bay in Zhejiang Province, China.

Figure 11 .
Figure 11.Magnitude of mean gradual changes (green bar in (a)) and mean abrupt changes (red bar in (a)) for different land cover categories, and the frequency of time series models classified into the same land cover category (green bar in (b)), and the land cover categories they changed to (red bar in (b)) (Dryl., dryland; Cons., construction land).

Figure 10 .
Figure 10.Annual mean and percentage of positive gradual and abrupt NDVI changes at Hangzhou Bay in Zhejiang Province, China.

Figure 10 .
Figure 10.Annual mean and percentage of positive gradual and abrupt NDVI changes at Hangzhou Bay in Zhejiang Province, China.

Figure 11 .
Figure 11.Magnitude of mean gradual changes (green bar in (a)) and mean abrupt changes (red bar in (a)) for different land cover categories, and the frequency of time series models classified into the same land cover category (green bar in (b)), and the land cover categories they changed to (red bar in (b)) (Dryl., dryland; Cons., construction land).

Figure 11 .
Figure 11.Magnitude of mean gradual changes (green bar in (a)) and mean abrupt changes (red bar in (a)) for different land cover categories, and the frequency of time series models classified into the same land cover category (green bar in (b)), and the land cover categories they changed to (red bar in (b)) (note: Dryl., dryland; Cons., construction land).

Table 1 .
Description of the designed land cover classification system in this research.

Table 2 .
Accuracy assessment of change detection in the spatial domain.

Table 3 .
Accuracy assessment of wetland classification result.

Table 4 .
Areas of different change times at Hangzhou Bay in Zhejiang Province, China, between 1987 and 2015.
1The numbers in parentheses are the percentages of changed areas at different change times to total changed area.

Table 4 .
Areas of different change times at Hangzhou Bay in Zhejiang Province, China, between 1987 and 2015.

Table 4 .
Areas of different change times at Hangzhou Bay in Zhejiang Province, China, between 1987 and 2015.TypesChange Times Area (ha) Percentage (%)

Table 5 .
Land cover change matrix at Hangzhou Bay in Zhejiang Province, China, during the period of 1987-2015 (ha).Note: Dryl., dryland; Cons., construction land; Sum_in., increased area from other classes; Net chg., net changed area of each class; Sum_de., decreased area due to conversion to other classes.

Table 5 .
Land cover change matrix at Hangzhou Bay in Zhejiang Province, China, during the period of 1987-2015 (ha).Note: Dryl., dryland; Cons., construction land; Sum_in., increased area from other classes; Net chg., net changed area of each class; Sum_de., decreased area due to conversion to other classes.

Table 6 .
The mean and positive percentage of total, gradual and abrupt NDVI at Hangzhou Bay in Zhejiang Province, China.

Table 6 .
The mean and positive percentage of total, gradual and abrupt NDVI at Hangzhou Bay in Zhejiang Province, China.