Mapping Annual Land Use and Land Cover Changes in the Yangtze Estuary Region Using an Object-Based Classiﬁcation Framework and Landsat Time Series Data

: A system understanding of the patterns, causes, and trends of long-term land use and land cover (LULC) change at the regional scale is essential for policy makers to address the growing challenges of local sustainability and global climate change. However, it still remains a challenge for estuarine and coastal regions due to the lack of appropriate approaches to consistently generate accurate and long-term LULC maps. In this work, an object-based classiﬁcation framework was designed to mapping annual LULC changes in the Yangtze River estuary region from 1985–2016 using Landsat time series data. Characteristics of the inter-annual changes of LULC was then analyzed. The results showed that the object-based classiﬁcation framework could accurately produce annual time series of LULC maps with overall accuracies over 86% for all single-year classiﬁcations. Results also indicated that the annual LULC maps enabled the clear depiction of the long-term variability of LULC and could be used to monitor the gradual changes that would not be observed using bi-temporal or sparse time series maps. Speciﬁcally, the impervious area rapidly increased from 6.42% to 22.55% of the total land area from 1985 to 2016, whereas the cropland area dramatically decreased from 80.61% to 55.44%. In contrast to the area of forest and grassland, which almost tripled, the area of inland water remained consistent from 1985 to 2008 and slightly increased from 2008 to 2016. However, the area of coastal marshes and barren tidal ﬂats varied with large ﬂuctuations.


Introduction
Estuaries are the interfaces between marine, freshwater, and terrestrial ecosystems and are welcomed as intermediate transitional zones with various ecosystem service functions [1]. Affected by climate change and rapidly urbanization, the estuarine and coastal regions are facing with many degradation risks, such as wetland shrinking, environmental pollution and coastal erosion [2,3]. To alleviate these problems, long-term dense monitoring of land use and land cover (LULC) change is of crucial important, because it provides essential information for depicting the history, current situation and future of LULC change, and for understanding biogeochemical processes and the mechanisms of LULC changes [4,5]. LULC products. Thus, approaches for consistent classification of image time series urgently need to develop for monitoring long-term LULC changes in estuarine and coastal regions.
In this study, a typical estuarine and coastal region, namely the Yangtze River estuary, was selected as the study area to explore these challenges. The Yangtze River estuary, one of the most famous estuaries in the world, is an ecologically and economically important region, which provides a wide range of ecosystem services, such as habitat for coastal vegetation and animals, carbon sequestration, and food supply [1,2]. With rapid urbanization and accelerated climate change, the LULC pattern in the study area has underwent unprecedented changes in past decades. Previous studies have shown that the high rate of LULC conversion in this area have created a series of environmental problems, such as coastal wetland degradation, water pollution and biodiversity loss [32]. Unfortunately, due to the lack of a system understanding of long-term LULC changes, the dynamic linkages between LULC changes and its environmental responses are poorly understood. Therefore, it is of the utmost importance for the mapping and monitoring of long-term LULC to optimize management strategies.
In this work, to explore ways to solve the challenges mentioned above, an object-based classification framework was developed to mapping long-term LULC dynamics of the Yangtze River estuary region using annual time-series of Landsat data. The specific goal is (1) to develop an object-based framework for generating consistent annual Landsat time series LULC maps by integrating four of prevailing LULC mapping techniques (object-based image analysis, hierarchical classification, and updating and backdating approaches) and (2) to quantify and assess annual LULC change in the Yangtze River estuary region from 1985 to 2016.

Study Area
The Yangtze River estuary, one of the most famous coastal and estuarine ecoregions in the world, was selected as the study area. The region is located in East China, and the study area includes the Shanghai municipality and portions of Suzhou and Nantong in Jiangsu Province. The estuary extends from 30 • 33 to 32 • 8 N in latitude and from 120 • 50 to 122 • 7 E in longitude, covering an area of 13929.75 km 2 ( Figure 1). Topographically, the region is composed of plains, except for a few hills located in the southwest. The climate is a marine subtropical monsoon climate, where clouds are frequent during the entire year. The region covers complex urban impervious materials, cropland, forest, water body and grass, as well as a variety of wetlands. However, the Yangtze River estuary has become a vulnerable region due to intense human activities and uncertain climate change.

Landsat Image Preprocessing
A consistent Landsat data preprocessing process was developed to process all the Landsat images. The process included atmospheric correction, image co-registration, poor-quality observation detection, gap filling, mosaicking and data clipping. Landsat TM (Landsat Thematic Mapper)/ETM+ (Landsat Enhanced Thematic Mapper Plus) bands 1-5 and 7 and bands 2-7 for Landsat OLI (Landsat 8 Operational Land Imager) were used in this work.

Image Selection
The Landsat images used in this study were acquired from the United States Geological Survey (USGS, http://glovis.usgs.gov/). All Landsat data covering the study area were assessed from 1984 to 2016. Two issues were considered in selecting high-quality images (Table 1): (1) the cloud contamination of a selected image was < 10% (mainland cover parts) as much as possible in the study area; (2) the selected images should cover the three key growing seasons of natural vegetation in each year to improve the classification accuracy [33]. These seasons include the early growing season ("spring," March to May), peak growing season ("summer," June to September), and late growing season ("autumn," October to December). However, due to cloud contamination, cloud-free images were not available for some years during all three key phenological phases. Thus, images of adjacent years were used as alternatives. In addition, one scene at the lowest tide for each year was selected for tidal flat extraction (Table 1).

Landsat Image Preprocessing
A consistent Landsat data preprocessing process was developed to process all the Landsat images. The process included atmospheric correction, image co-registration, poor-quality observation detection, gap filling, mosaicking and data clipping. Landsat TM (Landsat Thematic Mapper)/ETM+ (Landsat Enhanced Thematic Mapper Plus) bands 1-5 and 7 and bands 2-7 for Landsat OLI (Landsat 8 Operational Land Imager) were used in this work.

Image Selection
The Landsat images used in this study were acquired from the United States Geological Survey (USGS, http://glovis.usgs.gov/). All Landsat data covering the study area were assessed from 1984 to 2016. Two issues were considered in selecting high-quality images (Table 1): (1) the cloud contamination of a selected image was < 10% (mainland cover parts) as much as possible in the study area; (2) the selected images should cover the three key growing seasons of natural vegetation in each year to improve the classification accuracy [33]. These seasons include the early growing season ("spring," March to May), peak growing season ("summer," June to September), and late growing season ("autumn," October to December). However, due to cloud contamination, cloud-free images were not available for some years during all three key phenological phases. Thus, images of adjacent years were used as alternatives. In addition, one scene at the lowest tide for each year was selected for tidal flat extraction (Table 1).

Atmospheric Correction and Geometric Correction
For all selected Landsat images, the Fast Line of Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm was chosen to perform atmospheric corrections, which can accurately correct the radiation errors caused by various reasons [34]. All images were geometrically registered into a UTM projection (zone 51, North) with geometric errors of less than one pixel. Thus, all images have the same coordinate system.

Gap Filling and Cloud Removal
The failure of the Scan Line Corrector (SLC) on Landsat 7's ETM+ sensor in 2003, caused to approximately 22% of the pixels per scene not being scanned. The SLC failure has reduced the capacity of those images for mapping LULC dynamics. Thus, methods for gap filling were needed to address the missing values before image consequent processing. The gap-filling approach developed by Garcia [35] was used to solve the missing values in this study. This method is a three-dimensional spatiotemporal interpolation algorithm. The algorithm is a penalized least-squares regression based on a three-dimensional discrete cosine transform. It can maintain the original signal where data exist and minimize smoothing effects where new data are interpolated [36]. This algorithm has been successfully applied in previous studies [37]. Compared with the existing gap-filling methods [38], Garcia's method did not require ancillary data to simulate the missing pixel values and can accurately restore the values of unscanned pixels.
In this study, clouds and cloud shadows were visually detected and removed through the assembly of cloud-free observations from multi-temporal imagery of the same year.

Calculation of Spectral Indices
Several spectral indices were selected for LULC mapping in this study: the Enhanced Vegetation Index (EVI) [39], Modified Normalized Difference Water Index (MNDWI) [40], and Normalized Difference Build-Up Index (NDBI) [41]. The EVI was used to separate vegetation types from not-vegetation types, the MNDWI was used to extract water; and impervious and pervious surfaces were mapped with the NDBI. To minimize classification errors due to phenological differences, the effect of individual precipitation events and image temporal differences [42], the average EVI (EVIave), the average MNDWI (MNDWIave), and the average NDBI (NDBIave) were calculated. These spectral indices were calculated using the Landsat surface reflectance based on the following equations and spectral bands: where ρ BLUE , ρ GREEN , ρ RED , ρ NIR , and ρ SWIR1 are the surface reflectance values of the blue band, green band, red band, near infrared band, and shortwave infrared band (band 6 for Landsat OLI, band 5 for Landsat TM/ETM+), respectively.

Classification Scheme
Based on the research aims, features of the acquired images and field surveys, a nine-class scheme was designed for the study area (Table 2) based on a modified version of the scheme proposed by Gong et al. [43].

Object-Based Classification Framework
An object-based framework was designed to integrate object-based image analysis, hierarchical classification, updating and backdating approaches for LULC mapping in this study, as shown in Figure 2. In the framework, multi-temporal imagery was first segmented to generate image objects. Then, a rule-based hierarchical classification method was implemented to generate the initial LULC maps. In this stage, lots of prior knowledge obtained from field survey data and historical maps, such as spectral characteristics, phonological metrics, and spatial features for each class, was incorporated into the rule-based classifier to improve the classification accuracy. Trial and error methods are used to select the classification parameters and determine the optimal threshold values to classify each types, but this approach may increase the chance for classification inconsistencies. To overcome these issues, updating and backdating approaches were used to reduce the misclassification errors and generate more consistent LULC maps. After these steps, an accuracy assessment was performed for the LULC products. Finally, time series LULC products were used for change analysis.

Impervious
Built-up land, roads, transportation, residential, commercial services, industrial areas and settlements in villages Cropland Areas cultivated with crops such as rice, beans, and maize

Multi-Temporal Segmentation
In this study, image segmentation was implemented in eCognition®Developer 8.7 software. The multiresolution segmentation algorithm was selected as the image segmentation algorithm, which was proposed by Baatz and Schäpe [44]. The algorithm merges pixels based on a heterogeneity criterion, which can be adjusted by selecting the scale parameter, weights of spectral bands, spectral factor, and compactness factor [45,46]. Following Drǎguţ et al. [47], the ESP tool was used to aid in determining appropriate parameter settings for image segmentation. To guarantee the homogeneity of objects across the 32 years of Landsat time series, all time-series images were processed using the same parameters that were sufficiently detailed to map all the LULC classes in this study. Namely, multi-temporal images ("spring," "summer," and "autumn") were segmented in each year using 18 bands with equal weighting, a shape parameter of 0.1, compactness of 0.6, and scale factor of 20.

Object-Based Hierarchical Classification Stage
After image segmentation, an object-based hierarchical classification approach was implemented to generate the initial LULC products. An overview of the LULC mapping rule set is shown in Figure 2. In Figure 2, four steps are applied in the rule set. In the first step, water was extracted using the MNDWI. Considering the seasonal variations of water bodies, the summer images were used as reference images because summer is the rainy season in the study area. In the second step, a hybrid method was used within each land class to separate "coastal wetlands" from "not coastal wetlands." The algorithm named "assign class by thematic layer" within eCognition was first used to extract the 'coastal wetlands', of which the input vector was the coastal wetland layer from the previous year. The final "coastal wetland" boundaries were modified by visual interpretation. In the water category, the algorithm named "find enclosed by class" within eCognition was used to separate inland water and nearshore water. In the third step, within the "not coastal wetlands" class, "vegetation" was separated from "not vegetation" using EVI summer . Within the "coastal wetlands," coastal marshes were distinguished from barren tidal flats using EVI summer . In the final step, within the "not vegetation" type, harvest cropland was separated from the other two categories through analysis of the multi-temporal EVI and NDBI. Notably, high EVI and low NDBI values were obtained when harvest cropland was in the replantation phase. Built-up land was separated from other barren land using the NDBI. Within the vegetation category, a method similar to that adopted to classify "coastal wetlands" was applied to separate forest and grass from cropland. After hierarchical classification, error analysis showed that a certain amount of error could still be observed in the initial classified maps. To reduce these errors, updating and backdating approaches were used to generate more accurate LULC products (Figure 2). In this stage, classification error detection analysis was implemented, and 15 major misclassification scenarios occurred (Figure 3) (for coastal marshes, barren tidal flats and nearshore water, good classification results were achieved from the previous stage, so these three classes did not require updating/backdating). To perform updating/backdating, the overlay analysis method was first used to determine "correct classification objects" and "uncertain objects." Objects with identical classification types (attributes) between the reference map (T0) ( Table 1) and the initial classified map (Ti) were assigned as "correct classification objects." Objects with different classification types between Ti and T0 were labeled as "uncertain objects." updating/backdating). To perform updating/backdating, the overlay analysis method was first used to determine "correct classification objects" and "uncertain objects." Objects with identical classification types (attributes) between the reference map (T0) ( Table 1) and the initial classified map (Ti) were assigned as "correct classification objects." Objects with different classification types between Ti and T0 were labeled as "uncertain objects." As shown in Figure 3, the "correct classification objects" maintained the original classified results when performing updating/backdating. Among the "uncertain objects," only the misclassified objects required updating/backdating for classification correction, and the remaining "uncertain objects" maintained the original classified results (Figure 3). A combination of rule-based methods and a visual interpretation approach was used to identify and reclassify the misclassified objects from the "uncertain objects" ( Table 3). The visual interpretation approach of updating/backdating error objects was used through manual editing of the initial classified layer. The rule-based approach is a semi-automatic reclassified process. For example, if one misclassification object requires an update to impervious land, two criteria must be met: (a) the attribute of the object in T0 must be impervious, and (b) the values of EVIave and NDBIave in Ti must be less than 0.2 and greater than −0.1, respectively. In this stage, the works of Yu et al., Linke and McDermid. and Linke et al. are recommended as references for understanding more details of updating and backdating methods [48][49][50]. The reference maps used in this stage were generated by visual interpretation of the initial classified maps based on multi-temporal Landsat imagery and Google Earth maps.   As shown in Figure 3, the "correct classification objects" maintained the original classified results when performing updating/backdating. Among the "uncertain objects," only the misclassified objects required updating/backdating for classification correction, and the remaining "uncertain objects" maintained the original classified results (Figure 3). A combination of rule-based methods and a visual interpretation approach was used to identify and reclassify the misclassified objects from the "uncertain objects" ( Table 3). The visual interpretation approach of updating/backdating error objects was used through manual editing of the initial classified layer. The rule-based approach is a semi-automatic reclassified process. For example, if one misclassification object requires an update to impervious land, two criteria must be met: (a) the attribute of the object in T0 must be impervious, and (b) the values of EVIave and NDBIave in Ti must be less than 0.2 and greater than −0.1, respectively. In this stage, the works of Yu et al., Linke and McDermid. and Linke et al. are recommended as references for understanding more details of updating and backdating methods [48][49][50]. The reference maps used in this stage were generated by visual interpretation of the initial classified maps based on multi-temporal Landsat imagery and Google Earth maps.

Accuracy Assessment
Orthorectified aerial photos, high-resolution imagery from Google Earth, and field survey data were combined to produce reference data. Since available historical high-resolution imagery cannot cover the whole study area for the period from 1985 to 1998, reference data were also created from multi-temporal Landsat imagery and historical maps through visual interpretation in those period. The method proposed by Cochran was used to determine the total sample size [51]. Once the total sample size was determined, the approach proposed by Olofsson et al. [52] was used to determine the allocation of samples for specific categories. Namely, we allocated 50 samples per rare class (forest, grass, coastal marshes and unused land), and a sample size of n-200 was then proportionally allocated to the remaining classes. Accordingly, we used a random sampling scheme to generate the checkpoints for each class. Because the accuracy assessment required very intensive visual analysis and long-term time series and high spatial resolution data were not all available in this study, accuracy assessments were only conducted for the LULC maps from 1985,1992,1998,2005,2011, and 2016. The sampling design scheme of this study is shown in Table 4. Accuracy was assessed using a confusion matrix to calculate the overall accuracies, user's and producer's accuracies, and the Kappa coefficient [53]. An accuracy assessment of the LULC products indicated that the map accuracies from the object-based classification framework approach were high (Supplementary Table S1). The overall accuracies (OA) were 87.82%, 87.43%, 85.77%, 86.30%, 86.74%, and 87. 85% in 1985, 1992, 1998, 2005, 2011, and 2016, respectively, and the Kappa coefficients were 0.77, 0.78, 0.76, 0.78, 0.79, and 0.82, respectively. Cropland, inland water, coastal mashes and barren tidal flats were accurately classified, with user's accuracies (UA) and producer's accuracies (PA) of greater than 80% for all tested maps. For impervious cover, PA ranged from 77.63% to 87.21%, and UA ranged from 61.98% to 87.85%. The relatively low classification accuracy of impervious cover in 1985 and 1992 was mostly caused by low-density impervious regions. Such regions were often confused with cropland due to the limited spatial resolution and image quality. The rare LULC types, including forest, grass and other barren land, were relatively poorly classified, and the associated UAs were generally lower than 75%.

Long-Term LULC Dynamics
The spatial distributions of LULC types are clearly demonstrated in Figure 4, and the long-term changes are illustrated across each panel. The LULC classification results showed that great change was observed across the study area, especially change near the coastlines and peri-urban region. Specifically, cropland was the most extensive LULC type, but it continuously and rapidly decreased from 1985 to 2016. In contrast, impervious cover continuously and significantly expanded due to rapidly urbanization. Figure 4 also indicates that forest and grass substantially increased during the study period, although the proportions are still very small. The spatiotemporal distribution of coastal wetlands (including coastal marshes and barren tidal flats) varied and exhibited large fluctuations. Other barren land displayed an upward trend due to large-scale reclamation and the formation of increasing amounts of reclaimed lands from 2000 to 2016.
Quantitatively speaking, the land area of the study region increased by 730.13 km 2 from 1985 to 2016, with an average annual increase of 23.55 km 2 . This increase was due to land reclaimed from coastal wetlands (coastal marshes, bare tidal flat and nearshore water). For specific categories, the cropland area continuously decreased from 8076.63 km 2 (80.61%) in 1985 to 5959.33 km 2 (55.44%) in 2016 ( Figure 5), with an average annual decrease of 68.3 km 2 , and had a trend of continuous decrease (Slope = −65.96, R 2 = 0.99). The impervious coverage increased from 642.84 km 2 (6.42%) in 1985 to 3140.88 km 2 (22.55%) in 2016, with an average annual increase of 80.58 km 2 , thus resulting in a trend to continuous increase (Slope = 80.95, R 2 = 0.99). This expansion was mainly at the expense of cropland. The area of forest and grass increased by more than triple the original area during this period due to the long-term implementation of urban greening. The proportion of inland water was stable from 1985 to 2008 and then experienced a slight increase from 2008 to 2016. The area of coastal marshes fluctuated. Barren tidal flat areas underwent continued declining trends due to long-term reclamation activities (Slope = −6.41, R 2 = 0.55). Sustainability 2020, 12, 659 12 of 18

Monitoring Gradual Changes with Long-Term Time Series of LULC Products
Compared with most previous studies that used bi-temporal or sparse time series images [28], this work contributed to the analysis of LULC changes by using annual, long-term Landsat time series data. With annual time series LULC maps, the long-term (including gradual changes) and short-term trends (e.g., sudden change) in LULC dynamics can be monitored. Figure 4 clearly reflects the continuous spatiotemporal process of LULC change in the Yangtze River estuary from a time series perspective. More importantly, gradual changes that would not be observed using bi-temporal or sparse time series maps can be detected with LULC products produced in this paper. Here, we consider LULC gradual changes to be those that are occurring over a period of several years or more, usually relating gradual LULC transitions over that time period. At a conceptual level, some of the most notable example of gradual change involve to LULC undergoing a gradual shift through time, such as the process of re-growing saltmarsh vegetation after damage by typhoon. To present the details of these gradual changes, the generated maps in a selected region (Chongming Dongtan wetland natural reserve) were enlarged and illustrated, as shown in Figure 6. Figure 6 not only shows the sudden change of coastal marshes by disease or natural disasters or large-scale reclamation (denoted by the red ellipse region in Figure 6) but also shows regrowth and long-term recovery from these events (denoted by the red rectangle and circle regions in Figure 6). It is difficult to map such gradual changes using bi-temporal or sparse time series maps.

Monitoring Gradual Changes with Long-Term Time Series of LULC Products
Compared with most previous studies that used bi-temporal or sparse time series images [28], this work contributed to the analysis of LULC changes by using annual, long-term Landsat time series data. With annual time series LULC maps, the long-term (including gradual changes) and short-term trends (e.g., sudden change) in LULC dynamics can be monitored. Figure 4 clearly reflects the continuous spatiotemporal process of LULC change in the Yangtze River estuary from a time series perspective. More importantly, gradual changes that would not be observed using bi-temporal or sparse time series maps can be detected with LULC products produced in this paper. Here, we consider LULC gradual changes to be those that are occurring over a period of several years or more, usually relating gradual LULC transitions over that time period. At a conceptual level, some of the most notable example of gradual change involve to LULC undergoing a gradual shift through time, such as the process of re-growing saltmarsh vegetation after damage by typhoon. To present the details of these gradual changes, the generated maps in a selected region (Chongming Dongtan wetland natural reserve) were enlarged and illustrated, as shown in Figure 6. Figure 6 not only shows the sudden change of coastal marshes by disease or natural disasters or large-scale reclamation (denoted by the red ellipse region in Figure 6) but also shows regrowth and long-term recovery from these events (denoted by the red rectangle and circle regions in Figure 6). It is difficult to map such gradual changes using bi-temporal or sparse time series maps.

Advantages and Disadvantages of the Proposed Framework
In this study, an object-based classification framework was developed to generate a consistent Landsat time series for the long-term monitoring of LULC changes in the Yangtze River estuary region. This objective was completed through a combination of object-based image analysis, hierarchical classification, and updating and backdating approaches. Several advantages have been made to generate consistent long-term LULC maps in the proposed framework. First, compared to the traditional supervised classification algorithms [17], the proposed framework can make better use of spatial information such as shape and context, which can improve classification accuracy for complex LULC mapping. Second, compared to popular machine learning and deep learning algorithms, the proposed method does not need training samples, which is difficult to select training samples for these popular algorithms [54]. Third, object-based classification framework can largely reduce the errors caused by spatial misregistration, which is a common phenomenon and difficult to completely eliminate when using pix-based methods [55]. Finally, the proposed framework (updating and backdating stages) allow for a more straightforward utilization of prior knowledge. For example, the study area has not obviously appeared the phenomenon of counter-urbanization, thus we can use updating approach to improve the mapping accuracy of built up class (The classification result of reference image can be directly updated as the result of corresponding misclassified objects from the "uncertain objects").
However, the approach has some disadvantages that must be addressed. Data gaps (due to cloud contamination and gap filling) and temporal difference problems in images are the main sources of classification error in this study, and novel methods, such as data fusion, should be introduced or proposed to address these problems. Another issue in the proposed method is related to the mixed pixel problem. Due to complex LULC types in estuarine and coastal regions, mixed pixels are very common and may result in the inaccurate estimation of LULC. Hence, sub-pixel analysis may be useful to overcome these mixed pixel problems in future studies. In addition, compared to LULC mapping, using a Google Earth Engine-like platform [56,57], the proposed method is still a time-consuming and laborious work.

Implications for Ecological Restoration
Under the combined effects of human activities and climate change, most of estuarine and coastal ecosystems are in a degraded state. Ecological restoration is widely used to reverse the degraded state. Understanding annual LULC changes is the first step to have some insight on ecological degeneration and conducting restoration projects. Analyzing a long-term time series of LULC products is critical for predicting further LULC changes in management of ecosystem services [16,23], which is important for decision making during the implementation of ecological restoration projects. Information provided by the annual LULC products can be used to assess restoration effects associated with ecosystem functions and services. On top of that, annual LULC products provide basic data for understanding the mechanisms and processes of ecosystem restoration.

Conclusions and Future Work
In summary, an object-based classification framework that integrates object-based image analysis, hierarchical classification, and updating and backdating approaches was developed for mapping long-term LULC changes in the Yangtze River estuary region from 1985-2016 using annual Landsat time series images. The results showed that the annual time series of LULC maps could clearly reflect the continuous spatiotemporal change process of LULC in the study area that would not have been observed using bi-temporal or sparse time series maps. The impervious area rapidly increased from 6.42% to 22.55% of the total land area from 1985 to 2016, and the cropland area dramatically decreased from 80.61% to 55.44%. Additionally, the area of forest and grass increased by more than three times the original area. The proportion of inland water remained stable from 1985 to 2008 and exhibited a slight increase from 2008 to 2016. The area of coastal marshes and barren tidal flats displayed fluctuating change, and barren tidal flats exhibited a rapidly decreasing trend. Furthermore, the time series of LULC products revealed gradual changes that are valuable for both ecological and environmental management. An accuracy assessment of single-year classifications indicated that the designed framework provides an efficient method with OAs of over 86%.This study also shows the limitations of the proposed framework such as not full automatic and mixed pixel problem. Further work will focus on the change mechanisms among different land cover classes and more robust and automatic classification algorithms require further development to map long-term LULC change.