A Classification of Tidal Flat Wetland Vegetation Combining Phenological Features with Google Earth Engine

Abstract

The composition and distribution of wetland vegetation is critical for ecosystem diversity and sustainable development. However, tidal flat wetland environments are complex, and obtaining effective satellite imagery is challenging due to the high cloud coverage. Moreover, it is difficult to acquire phenological feature data and extract species-level wetland vegetation information by using only spectral data or individual images. To solve these limitations, statistical features, temporal features, and phenological features of multiple Landsat 8 time-series images obtained via the Google Earth Engine (GEE) platform were compared to extract species-level wetland vegetation information from Chongming Island, China. The results indicated that (1) a harmonic model obtained the phenological characteristics of wetland vegetation better than the raw vegetation index (VI) and the Savitzky–Golay (SG) smoothing method; (2) classification based on the combination of the three features provided the highest overall accuracy (85.54%), and the phenological features (represented by the amplitude and phase of the harmonic model) had the greatest impact on the classification; and (3) the classification result from the senescence period was more accurate than that from the green period, but the annual mapping result on all seasons was the most accurate. The method described in this study can be applied to overcome the impacts of the complex environment in tidal flat wetlands and to effectively classify wetland vegetation species using GEE. This study could be used as a reference for the analysis of the phenological features of other areas or vegetation types.

1. Introduction

Vegetation is the source of primary productivity in wetland ecosystems and has important ecological functions including water conservation, regional climate regulation, siltation promotion, biodiversity protection, and resource production [1,2,3]. China’s wetland vegetation resources are precious, and the biological invasion of Spartina alterniflora in most coastal zones threatens the growth of local native vegetation. Therefore, accurate and timely information on the distribution of wetland vegetation is very important in protecting the diversity of wetland vegetation and sustainably developing wetland ecosystems.

Satellite remote sensing images can be used to effectively monitor vegetation types such as forests and cropland [4,5]. However, the environment of wetland vegetation is complex due to biological invasion of S. alterniflora and the difference in the underlying surface water level, which are known to affect the radiation transfer of the vegetation canopy. Furthermore, the east coast of Chongming Island is more affected by sea water, while the west coast is more affected by Yangtze River water, which leads to different growth patterns of the same species. Thus, even the same species can show different spectral features in different places and can have varied spectra because of the variability in the vegetation density, water level, and soil moisture, while different types of wetland vegetation may have similar spectra. Therefore, it is difficult to identify wetland vegetation with spectral information alone. Many studies have focused on extraction of wetland vegetation information through remote sensing classification methods and data sources [6,7,8], but most of these classifications have been based on single images rather than considering a time series of images. However, vegetation classification results may vary among different growth periods. Therefore, how to effectively extract vegetation information and obtain reliable mapping results from remote sensing images in complex wetland environments is a problem that remains to be solved.

Site surveys have indicated that the vegetation in tidal flat wetlands on Chongming Island has obvious phenological characteristics such as a green period and a senescence period. Remote sensing classification based on the growth characteristics of vegetation in different specific phenological periods has been successfully applied to cropland, bamboo, and forest areas [9,10]. However, the changes in wetland vegetation are different from the regular changes observed in croplands due to harvesting as well as from the obvious seasonal changes in deciduous forests. Wetland vegetation usually grows naturally on tidal flats and exhibits typical perennial herbaceous vegetation features. Few studies have used phenological characteristics to identify wetland vegetation. Lumbierres et al. [11] used the surface phenology of wetland vegetation from the Moderate-Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) images to estimate the aboveground biomass of seasonal swamps and found that the stability of estimation combined with a maximum-value approach provided the best results. Shen et al. [12] used NDVI and climate data to study the spatiotemporal changes in the start date of vegetation growth and discussed the climate change impacts on the start date of growth in freshwater marshes in Northeast China. Li et al. [13] used Sentinel-2 and unmanned aerial vehicle (UAV) images to study the distribution of mangroves in the Yellow River estuary during the period of maximal phenological difference. At present, vegetation indices (VIs) are commonly used to analyze wetland phenological characteristics. In these approaches, a critical period with significant phenological differences is selected, and the feature collections of different phenology are constructed based on pixels or images to extract different types of wetland vegetation. However, it is unclear whether phenological features can be directly used for image classification, which requires further study of the contribution of phenological features to wetland vegetation identification.

The studies of plant phenology often require a large number of time-series satellite images to analyze plant phenological characteristics. Traditional remote sensing methods are considered to have poor timeliness and lack portability because it may require considerable time, manpower, and energy to download and preprocess the image collections. However, Google Earth Engine (GEE) is a cloud-based geospatial computing platform that includes satellite observation images and products of all levels from the past 40 years to the present. Moreover, GEE can flexibly establish machine learning classification algorithms at high speed that run in parallel and can quickly eliminate the influence of clouds and shadows in images. Some current studies have used this platform for vegetation identification. GEE can substantially simplify preprocessing and repetitive work by efficiently integrating the data and making full use of existing products [14]. Wang et al. [15] used various forest classification products provided by GEE and Landsat images to monitor and evaluate changes in forest disturbance in tropical areas over 30 years. Xie et al. [16] used GEE to analyze the changes in pasture vegetation cover over the past 30 years and provided data for the protection of degraded land. These studies indicate that GEE has significant advantages for regional classification with multitemporal features [17].

In summary, due to the growth differences of wetland vegetation and the complexity of wetland environments, several limitations such as the difficulty of quickly and efficiently extracting combined wetland vegetation phenology features for mapping cannot be overcome by solely using spectral information for vegetation identification. Therefore, in this study, the GEE platform with Landsat 8 surface reflectance products was used to analyze the phenological features of wetland vegetation, and the phenological features were added as bands to the image for classification. In addition, the phenological features were compared and combined with the statistical and temporal features and obtained annual vegetation maps. Next, the optimal combination of features was used to extract wetland vegetation information from the green period and the senescence period and to estimate the area of wetland vegetation. Finally, the contribution of each feature type to the classifications and the influence of phenological features on the extraction of wetland vegetation were analyzed. The purposes of this study were to (1) determine and analyze the phenological characteristics of the VIs of wetland vegetation; (2) use phenological features to identify wetland vegetation and select the optimal feature combination for annual mapping by comparing and combining phenological features with statistical and temporal features; and (3) assess the classifications obtained by including phenological features and the classifications from the seasonal phenological periods and analyze the differences among the seasonal mapping results.

2. Materials and Methods

2.1. Study Area

Chongming Island is located on the northeastern side of Shanghai, between 31°27′00′′ N to 31°51′15′′ N and 121°09′30′′ E to 121°54′00′′ E, as shown in Figure 1. It is surrounded by a river on three sides and borders the East China Sea. It is the largest alluvial estuary island in the world [18], at 80 km in length from east to west and 18 km wide from north to south, and covers a total area of 1411 km². Chongming Island has an average annual temperature of approximately 15.3 °C and an average annual precipitation of approximately 1003.7 mm. The climate is a mild and humid northern subtropical monsoon climate. Climate disasters such as drought, heavy rains, and typhoons are common in the area [19].

Chongming Island is the most well-developed tidal flat wetland in the Yangtze River estuary. The vegetation is distributed mainly in Xisha, Beiliuyao, Niupenggang, and Dongtan, and the rest is distributed in intermittently connected patches along the dam roundabout. The development of the coastal wetlands across the island varies by area. The tidal flat wetlands on the north bank are generally lower-quality than those on the south bank and exhibit low plant diversity. There are more reclamation areas on the northern narrow tidal flats, while the distribution of vegetation on the south part of the island is continuous and interrupted only by several docks, locks, and shipyards [20]. The tidal flat wetland vegetation on Chongming Island mainly comprises Scirpus mariqueter, S. alterniflora, Phragmites australis, trees, and cropland. Most of the trees are planted trees and are present mainly on the south part of Chongming Island including Metasequoia glyptostroboides, Taxodium distichum, Salix babylonica, and Taxodium ascendens.

2.2. Data and Methods

The Landsat 8 OLI atmospherically corrected surface reflectance (SR) product provided by GEE was used in this study. The experimental process is shown in Figure 2. It is difficult to use VI to directly reflect phenological changes in vegetation because of the complex environment and the different vegetation growth conditions. Therefore, this study selected five specific VIs, namely, the NDVI, ratio vegetation index (RVI), 1640 nm shortwave infrared vegetation index (NDWI1640), 2310 nm shortwave infrared vegetation index (NDWI2310), and soil normalized vegetation index (SAVI), and used harmonic analysis to extract the phenological differences in wetland vegetation. Then, the green period (from 30 May to 15 October) and senescence period (15 November to 30 April) were determined according to the vegetation phenological characteristics. The images used in this study are shown in

Table 1. Characteristics of the Landsat 8 Operational Land Imager (OLI) images used in the experiment.

Year	Path	Row	Cloud Cover (%)	Season	Year	Path	Row	Cloud Cover (%)	Season
2018-01-15	118	38	6.63	A,S	2019-01-18	118	38	0.61	A,S
2018-03-04	118	38	41.44	A,S	2019-03-07	118	38	41.75	A,S
2018-04-21	118	38	77.97	A,S	2019-03-23	118	38	15.99	A,S
2018-05-23	118	38	18.66	A	2019-04-08	118	38	41.30	A,S
2018-06-08	118	38	87.86	A,G	2019-05-10	118	38	6.97	A
2018-06-24	118	38	61.50	A,G	2019-06-11	118	38	55.17	A,G
2018-07-10	118	38	66.41	A,G	2019-06-27	118	38	93.79	A,G
2018-07-26	118	38	5.89	A,G	2019-07-29	118	38	10.87	A,G
2018-08-11	118	38	9.25	A,G	2019-08-14	118	38	3.30	A,G
2018-08-27	118	38	27.05	A,G	2019-08-30	118	38	7.94	A,G
2018-09-28	118	38	16.82	A,G	2019-10-17	118	38	54.22	A
2018-10-30	118	38	7.32	A	2019-11-02	118	38	15.12	A,S
2018-12-01	118	38	20.58	A,S	2019-12-04	118	38	6.38	A,S
2018-12-17	118	38	0.77	A,S	2019-01-25	119	38	65.95	A,S
2018-01-06	119	38	37.59	A,S	2019-02-26	119	38	42.47	A,S
2018-01-22	119	38	42.16	A,S	2019-03-14	119	38	96.98	A,S
2018-02-07	119	38	58.59	A,S	2019-03-30	119	38	99.99	A,S
2018-02-23	119	38	0.99	A,S	2019-04-15	119	38	6.15	A,S
2018-03-11	119	38	46.00	A,S	2019-05-01	119	38	28.89	A
2018-03-27	119	38	8.95	A,S	2019-06-02	119	38	38.03	A,G
2018-04-28	119	38	3.00	A,S	2019-07-04	119	38	99.98	A,G
2018-05-14	119	38	41.22	A	2019-08-05	119	38	42.03	A,G
2018-05-30	119	38	96.37	A,G	2019-08-21	119	38	5.85	A,G
2018-06-15	119	38	21.57	A,G	2019-09-06	119	38	76.00	A,G
2018-07-17	119	38	43.71	A,G	2019-09-22	119	38	74.05	A,G
2018-08-02	119	38	47.29	A,G	2019-10-24	119	38	16.52	A
2018-08-18	119	38	69.70	A,G	2019-11-09	119	38	0.15	A
2018-09-03	119	38	70.83	A,G	2019-12-11	119	38	0.20	A,S
2018-09-19	119	38	60.92	A,G	2019-12-27	119	38	5.12	A,S
2018-10-05	119	38	64.28	A,G
2018-10-21	119	38	70.36	A
2018-11-22	119	38	4.60	A,S
2018-12-24	119	38	2.44	A,S

Note: In the table, A represents all seasons, G represents the green period, and S represents the senescence period. The time of the green period is from 30 May to 15 October, and the senescence period is from 15 November to 30 April.

. GEE was used for pixel-based annual image analysis, and the images were composed of the pixels with the lowest cloud coverage for the period. Next, the amplitude and phase information of the VIs were used as the phenological features for identifying wetland vegetation. Furthermore, the classifications by phenological features, statistical features, temporal features, and their combinations were compared to select the best features for classification. Then, a random forest (RF) classifier was used for classification, and a confusion matrix was established to verify the wetland vegetation classification. Finally, the classifications for the two periods as determined with the harmonic model were compared with the classification for all seasons, and the impact of phenology on the classification was analyzed.

2.2.1. Data

Two field surveys were carried out from December 17 to 25 in 2018 and from November 12 to 17 2019 on Chongming Island. A total of 123 field survey sites were selected (as shown in Figure 1). The size of each site was at least 120 × 120 m and the distance between any two sites were set to far more than 30 m to ensure that there were at least four pixels in an OLI image falling in one site. The species at each site were recognized and recorded, and the locations of the sites were measured by a handheld positioning device with an accuracy of less than 0.5 m.

The identification of vegetation species of the sample data was based on either the field survey or the visual interpretation of Gaofen-2 PMS images with a spatial resolution of 4 m. All of the sample data were divided into two separate sets: one for the analyzing of the phenological differences between vegetation species, and the other for the classification. The first set consisted of a total of 465 sample pixels with good internal homogeneity that did not change within the study period including P. australis (145 pixels), S. alterniflora (110 pixels), S. mariqueter (79 pixels), trees (60 pixels), and cropland (71 pixels). All the pixels in the first set were located in the field survey sites. The second set for classification contained water bodies (178 pixels), tide flats (124 pixels), P. australis (272 pixels), S. alterniflora (135 pixels), S. mariqueter (71 pixels), trees (52 pixels), and croplands (17 pixels), which were randomly selected around Chongming Island.

Landsat 8 OLI atmospherically corrected SR product images with a spatial resolution of 30 m obtained through GEE were selected for this study. These images needed to be quality-filtered before analysis. The pixel quality and radiant saturation properties were used to remove cloud cover areas and cloud shadows [21,22]. The original reflectance value was divided by 10,000 and converted to a value between 0 and 1. Images of Chongming Island from January 1, 2018 to December 31, 2019, were used in this study. A total of 62 images met the requirements of this study, as shown in Table 1. The experiment used five main VIs (

Table 2. Vegetation index (VI) calculation formulas and associated Google Earth Engine (GEE) code.

VIs	Calculation Formula	GEE Code
NDVI [30]	(NIR-RED)/(NIR + RED)	normalizedDifference([‘B5’,’B4’])
NDWI2130 [29]	(NIR-SWIR2130)/(NIR + SWIR2130)	normalizedDifference([‘B5’,’B6’])
NDWI1640 [29]	(NIR-SWIR1640)/(NIR + SWIR1640)	normalizedDifference([‘B5’,’B7’]
SAVI [26]	1.5 × (NIR-RED)/(NIR + RED + 0.5)	Code equation ‘1.5 × (nir-red)/(nir + red + 0.5)’
RVI [28]	NIR/RED	Code equation ‘nir/red’

) including NDVI [23,24,25], which is sensitive to vegetation information; SAVI [25,26,27], which is sensitive to soil information and can be used to distinguish vegetation from other ground objects; RVI [28], which is sensitive to green vegetation; and NDWI1640 and NDWI2310 [27,29], which are sensitive to leaf water content and can reflect the differences between different vegetation types.

2.2.2. Phenological Analysis Method

The VI data had many data gaps in the time series due to the influence of clouds. It is difficult to use VIs to directly reflect the phenological differences among various vegetation types due to their fluctuation and interference caused by growth differences and the environment in the time series. This study first used the mean NDVI as an example to compare the capacity of two filtering methods, namely, a harmonic analysis of the time series (Hants) and the Savitzky–Golay (SG) smoothing method, in terms of the phenological characters extracted and to compare the fitting value with the real value using correlation analysis.

The fitted values of the Hants and SG filtering methods for the NDVI of the three types of wetland vegetation are shown in Figure 3. Both Hants and SG filtering performed well on P. australis, and the Pearson correlation coefficients reached 0.9659 and 0.9442, respectively, as shown in Figure 3a. However, the SG filter had a poor fitting effect for the S. alterniflora time-series images, with many data gaps. In Figure 3b, the SG method overestimated the NDVI in the lower-NDVI regions where the NDVI value was less than 0.5 and underestimated NDVI in the higher-NDVI (0.5–0.9) regions. Figure 3c shows that the two methods estimated the NDVI of S. mariqueter well, but that the fitting effect of the SG filtering method was slightly lower than that of Hants. Remote sensing images of Chongming Island are affected by its local microclimate; most areas are covered by high clouds throughout the year, resulting in few clear images. Therefore, there were many gaps in the time series. In general, Hants does not require high data integrity, but SG has strict data requirements. On Chongming Island, where time-series gaps are prone to occur, the fitting effect of Hants for the NDVI was better than that of the SG method. The four other vegetation indices were compared in this study, and similar results were obtained.

The wetland vegetation phenological characteristics extracted by the Hants and SG filtering methods were compared. The SG filter was found to simulate the NDVI value of wetland vegetation well, but was susceptible to the raw value distribution and had multiple inflection points (Figure 4), which affected the distinction of vegetation phenology periods. The Hants method combines Fourier transform with least-squares fitting. It can decompose the time spectrum of the vegetation index into several sine and cosine curves of different frequencies and superimpose the curves that can reflect the characteristics of the time series. The fitting of sequences of different time intervals can be effectively applied to the extraction of the phenological features of vegetation [31].

In summary, the Hants filtering method was used in this study to extract the phenological characteristics of wetland vegetation, and follow-up research was conducted based on the results. Harmonic models were used to determine the changes in the mean value of five VIs (in Table 2) in the Landsat time-series images and to analyze the phenological features of vegetation. Then, the periodic fluctuations of the VIs were simulated, and the amplitude and phase of the harmonic models were obtained to effectively reflect phenological differences among vegetation. The calculation formula is as follows:

$$Y={\mathit{a}}_{\mathbf{0}}+{\mathit{b}}_{\mathbf{1}}cos\left(\frac{2\pi t}{\mathit{T}}\right)+{\mathit{b}}_{\mathbf{2}}sin\left(\frac{2\pi t}{\mathit{T}}\right)+c_{1}t$$

(1)

$$Amplitude = {\left({\mathit{b}}_{\mathbf{1}}{}^{\mathbf{2}}+ {\mathit{b}}_{\mathbf{2}}{}^{\mathbf{2}}\right)}^{½}$$

(2)

$$Phase = atan\left({\mathit{b}}_{\mathbf{1}}/{\mathit{b}}_{\mathbf{2}}\right)$$

(3)

Note: ${\mathit{a}}_{\mathbf{0}}$ is the coefficient of the overall value for the VI; ${\mathit{b}}_{\mathbf{1}}$,${\mathit{b}}_{\mathbf{2}}$ are the intra-annual changes of the VI;$c_1$ is the inter-annual changes of the VI;$t$ is time; $Y$ is the fitted value corresponding to $t$; $\mathit{T}$ is the frequency (365 days). The experiment used the amplitude and phase of the harmonic model to represent the phenological features of vegetation.

2.2.3. Feature Selection

Multitemporal Landsat images are often used to extract the phenological characteristics of different vegetation growth periods. Previous studies have shown that the use of multiseason images based on vegetation phenological information can effectively increase the differentiation between vegetation types with similar spectra, thereby improving the classification accuracy for various wetland vegetation types [32]. Three groups of features were used in this study including three statistical variables (median, maximum, minimum), five temporal variables (25% and 75% of the time-series quantile [33] and average intervals of the time series at 25–50%, 50–75% and 25–75%), and two phenological variables (amplitude and phase). These three groups of features and their combinations were used to analyze the impact of each feature group on classifying wetland vegetation. At most, 10 feature variables were selected for each VI, and a total of 50 feature variables were tested. The feature variables and their codes in GEE are shown in

Table 3. Selection of vegetation index (VI) features.

Group	Name	Description	GEE Code
Statistical features	Median	Median of time-series	$\mathrm{ee}.\mathrm{Reducer}.\mathrm{median}()$
	Max	Maximum of time-series	$\mathrm{ee}.\mathrm{Reducer}.\max ()$
	Min	Minimum of time-series	$\mathrm{ee}.\mathrm{Reducer}.\min ()$
Temporal features	Percentile 25	The value at the 25% quartile of the time series	$\mathrm{ee}.\mathrm{Reducer}.\mathrm{percentile}\left(25,75\right)$
	Percentile 75	The value at the 75% quartile of the time series
	Mean Interval (25,75)	The mean value of the time series from 25% to 75%	$\mathrm{ee}.\mathrm{Reducer}.\mathrm{intervalMean}\left(25,75\right)$
	Mean Interval (25,50)	The mean value of the time series from 25% to 50%	$\mathrm{ee}.\mathrm{Reducer}.\mathrm{intervalMean}\left(25,50\right)$
	Mean Interval (50,75)	The mean value of the time series from 50% to 75%	$\mathrm{ee}.\mathrm{Reducer}.\mathrm{intervalMean}\left(25,75\right)$
Phenological features	Amplitude	Amplitude of VI time series filtered by Hants over a period of time	${\mathrm{Code} \mathrm{equation}}^{ }\mathrm{based} \mathrm{on} \left(2\right)$
	Phase	Phase of VI time series filtered by Hants over a period of time	${\mathrm{Code} \mathrm{equation}}^{ }\mathrm{based} \mathrm{on} \left(3\right)$

.

2.2.4. Classification and Validation

The RF classifier is an ensemble of multiple decision trees. For the input images, n trees will have n classification results. The classification of the RF classifier is determined by the mode value of classifications from decision trees. Since the results of multiple decision trees are integrated, an RF classifier can accept high-dimensional sample input and can evaluate the importance of each feature in the classification with better accuracy than other methods [34]. The number of decision trees needs to be set according to the actual situation when using the RF method for classification. The use of more decision trees provides higher classification accuracy, but also requires more calculations and a longer calculation time [35,36]. This study used a RF classifier constructed with 500 trees.

In total, 849 sample points were used for classification, of which 70% of the data by stratified sampling were used for calibration and 30% for verification [37]. Then, the accuracy of the classifications was evaluated. Classification accuracy refers to the proportion of pixels in the classified image that are correctly classified. This experiment used the confusion matrix method to verify the classification accuracy including the overall accuracy (OA), producer’s accuracy (PA), and user’s accuracy (UA).

3. Results

3.1. Analysis of Wetland Vegetation Phenological Characteristics

The predicted values of the wetland vegetation harmonic model from 2018 to 2020 were compared with the actual average values of the vegetation index, as shown in Figure 5, and indicated that the changes of phenology in croplands were weak, and its amplitude was the smallest. Compared with those of wetland vegetation and trees, the periodic changes in cropland phenology were slight. Thus, with suitable phenological features, cropland can be rapidly distinguished from wetland vegetation. The VI fluctuation trends of P. australis and trees were similar, but the VI values of trees were larger than those of P. australis. Comparing S. mariqueter and S. alterniflora revealed that the NDVI and SAVI values of S. alterniflora were always higher than those of S. mariqueter, while the fluctuation ranges of the NDWI1640, NDWI2310, and RVI values of S. alterniflora were always greater than those of S. mariqueter. Therefore, phenological information can be used to distinguish these two types of vegetation. Considering the similarity in amplitude between P. australis and trees and between S. alterniflora and S. mariqueter, the differences in the VIs of each pair were compared, and the phases were found to be obviously different, as shown in Figure 5.

3.2. Potential Separability of the Classes Using the Derived Variables

The above analysis indicated that there are obvious phenological differences between different types of vegetation. Thus, the amplitude and phase extracted from the harmonic model were used to represent the phenology of the vegetation, and they were added as a group of features for classification. The NDVI was used as an example for comparing the distribution characteristics of the vegetation samples in terms of the three feature groups (10 variables). Heterogeneous statistical and temporal features were observed within the same vegetation type in different areas of Chongming Island due to the complex environment, as shown in Figure 6. However, the phenological features reflected the differences between the different vegetation types well, reduced the influence of environmental interference, and overcame the issues caused by the complex wetland environment. In terms of statistical features and temporal features, the values of P. australis and S. mariqueter in the 25% and 75% quartiles overlapped greatly, but these two types of vegetation could be clearly distinguished by their phenological features. Although there was some overlap in the phenological features of P. australis and those of trees, the statistical and temporal features of P. australis and trees were clearly different. These results indicate that combining phenological features can help to distinguish vegetation types outside the dam.

3.3. Classifications under Different Feature Combinations Considering Phenological Features

This study found that when the three groups of features were separately added to the image, the best classification with the least number of features was obtained with phenological features (

Table 4. Classifications under three kinds of features and their combinations.

Accuracy	SF	TF	PF	SF + TF	SF + PF	TF + PF	SF + TF +PF
Number of features	18	25	10	43	28	35	53
OA (%)	79.78	73.62	80.32	79.22	82.84	83.94	85.54
P. australis PA (%)	85.18	85.71	90.91	90.18	83.81	95.04	88.46
P. australis UA (%)	78.63	68.29	81.48	73.19	82.24	80.99	82.14
S. alterniflora PA (%)	79.31	82.35	89.29	76.79	85.71	86.79	83.64
S. alterniflora UA (%)	80.70	66.67	83.33	91.49	90.57	90.20	92.00
S. mariqueter PA (%)	57.58	32.26	77.14	25.00	66.67	60.71	60.00
S. mariqueter UA (%)	61.29	50.00	77.14	37.50	66.67	65.38	66.67
Tree PA (%)	54.17	23.08	38.09	27.78	56.52	47.62	86.67
Tree UA (%)	76.47	54.54	100.00	55.56	61.90	76.92	81.25
Cropland PA (%)	62.50	11.11	100.00	75.00	100.00	100.00	83.33
Cropland UA (%)	100.00	50.00	66.67	100.00	100.00	100.00	100.00
Mudflat PA (%)	86.87	82.69	61.54	89.13	80.00	83.64	92.68
Mudflat UA (%)	73.58	78.18	76.19	78.85	83.02	80.70	79.17
Water PA (%)	87.50	87.18	81.71	88.89	94.03	81.69	89.19
Water UA (%)	92.10	95.77	78.82	95.52	87.50	96.67	98.51

Note: In the table, SF represents statistical features, TF represents temporal features, and PF represents phenological features.

). When phenological features were used for classification, the PA and UA of P. australis, S. alterniflora, and S. mariqueter were high, and the wetland vegetation could be identified well. However, using only phenological information for classification was insufficient to distinguish between P. australis and trees, and the distinction between water and tidal flats was poor. Therefore, the study further combined features to extract the wetland vegetation. In this study, combining only statistical features and temporal features helped distinguish between croplands and tidal flats, but did not improve the accuracy of wetland vegetation identification. When phenological features were combined with statistical features, S. alterniflora, S. mariqueter, and cropland could be better identified than when only statistical features were used, and the ability to distinguish between vegetation types was effectively improved. Compared with using only phenological features, using statistical features with phenological features could more effectively distinguish water bodies and tidal flats. When temporal features and phenological features were both used for classification, the OA, PA, and UA for the five types of vegetation were improved compared to those using only temporal features. Moreover, including temporal features also overcame the limitations of identifying tidal flats and water bodies only by phenological features. When the three feature types were combined for mapping in this study, the highest OA was achieved. With this combination, the phenological information obtained from the harmonics could be used to effectively distinguish wetland vegetation types, and statistical features and the temporal features could be used to distinguish between tidal flats and water. When combined, the three kinds of features could effectively distinguish wetland vegetation from trees and cropland and increase the distinction between trees and cropland.

In summary, classification with all three kinds of features exhibited the best performance in the identification of wetland vegetation outside the dam on Chongming Island. Thus, the vegetation distribution of Chongming Island was analyzed according to this result (Figure 7). Considering the administrative regions of Chongming Island, P. australis was found to be distributed throughout the tidal flats and wetlands outside the Chongming Island dam; it is mainly concentrated in Niupenggang, as shown in Figure 7a. The wetland vegetation in this area is mainly P. australis and has not been disturbed by the invasive species S. alterniflora. However, in the middle area of Niupenggang, croplands have been formed by artificial reclamation. There are rows of trees along the dam, distributed below the area of Figure 7a in Xisha, that are mainly mangroves planted to block winds and stabilize the bank. S. alterniflora is distributed mainly at Beiliuyao on Chongming Island, as shown in Figure 7b. S. alterniflora is distributed mainly near the coast at the middle and low tide levels. P. australis is distributed mainly at the edge of the dam at the high tide level and exhibits a banded distribution. In this area, there is ecological competition between the native species P. australis and the invasive species S. alterniflora. Another major distribution area of P. australis on Chongming Island is located at southwestern Dongtan on Chongming and extends to central Dongtan along the dam. S. mariqueter is concentrated mainly in the southeast corner of Dongtan. Most of the northern part of Dongtan is an ecological restoration area dominated by S. alterniflora that has been flooded for a long time. The distribution of vegetation communities at Dongtan is shown in Figure 7c.

3.4. Features That Contribute to Wetland Vegetation Identification

Table 4 shows that the classification accuracy of wetland vegetation identification was the highest when a combination of statistical variables, temporal variables, and phenological variables was used. Then, the influence of each VI and each kind of feature on the classification as well as the importance of each component were analyzed. This study was divided into two parts. In the first part, the five VIs were unchanged, and the feature types were removed from the variables one at a time to analyze the influence of each feature type on the classification, which is shown in Table 4. In the second part, the three kinds of features were unchanged, and the VIs were removed from the variables one at a time to analyze the influence of each VI on the result (Figure 8). The results of the first part of the analyses indicated that the OA of the classifications decreased when any one of the three kinds of features was removed, and the phenological features had the greatest impact on OA among the three kinds of features with the greatest decreases in accuracy. The removal of phenological features had the greatest impact on S. mariqueter classification, as its PA and UA were significantly reduced. When phenological features were removed, the PA and UA of S. alterniflora decreased; although the PA of P. australis slightly increased, the UA of P. australis was significantly reduced. When statistical features or temporal features were excluded, the PAs of the three types of wetland vegetation were slightly increased, but the corresponding UAs and OAs decreased. Therefore, the PA and UA of ground objects will be negatively impacted if temporal or statistical features are not included for classification.

The OA of the classification decreased when any kind of VI was excluded, and the removal of SAVI had the greatest impact on the classifications, followed by the removal of NDVI and NDWI2310. In this study, removing any VI reduced the PA of P. australis; although the UA of P. australis was improved by removing VIs, the reliability of the results decreased due to the reduced PA. Therefore, reducing the number of VIs used will lead to a reduction in the accuracy of vegetation classification. Analyzing the changes in S. alterniflora and S. mariqueter classification when VIs were removed revealed that their UAs decreased. Although their PAs increased when certain VIs were removed, their OAs decreased. Of the vegetation types, P. australis was the most sensitive to SAVI and NDVI. When SAVI and NDVI were removed, the PA and UA of P. australis decreased. S. alterniflora was sensitive to NDVI, which means that NDVI is critical for identifying S. alterniflora. When NDVI was removed, the PA and UA of S. alterniflora were significantly reduced. S. mariqueter was sensitive to SAVI, NDVI, and NDWI2310. When SAVI, NDVI, and NDWI2310 were removed, the PA of S. mariqueter increased slightly, while its UA decreased significantly.

3.5. Classification of the Green Period and Senescence Period with the Optimal Feature Combination

According to the phenological information obtained with the harmonic model and the vegetation growth characteristics (Figure 5), the images in this study were divided into two characteristic periods: the green period and the senescence period. The optimal feature combination was used to identify the vegetation types. The differences between the green period, senescence period, and all-season classifications were compared in

Table 5. Comparison of classifications among all seasons, the senescence period and the green period.

Accuracy	All Seasons	Senescence Period	Green Period
Number of images	63	24	26
OA (%)	85.54	81.57	78.59
P. australis PA (%)	88.46	88.99	93.41
P. australis UA (%)	82.14	78.86	72.03
S. alterniflora PA (%)	83.64	72.92	76.00
S. alterniflora UA (%)	92.00	87.50	90.48
S. mariqueter PA (%)	60.00	50.00	51.43
S. mariqueter UA (%)	66.67	75.00	75.00

, and the reasons for these differences were analyzed. Compared with those for the senescence period and the green period, the all-season classification was better, and the OA of the senescence period was higher than that of the green period. The PA for P. australis increased slightly, but the PAs for S. alterniflora and S. mariqueter both decreased. Therefore, the reliability of mapping S. alterniflora and S. mariqueter in the senescence period and the green period would be lower than that in all seasons. Although the PA for P. australis increased, its UA declined in the senescence and green periods. Therefore, the classifications of P. australis in the senescence and green periods were not as good as all seasons on the whole.

Comparing the classifications in the senescence period (Figure 9) and the green period (Figure 10) with that in all seasons (Figure 7) revealed that the classifications were different when they were performed in different periods. Considering the accuracy of the classifications shown in Table 5 revealed that the classification in the senescence period was similar to that in all seasons, but classification in all season could better distinguish wetland vegetation from trees and cropland. Moreover, the error of the senescence period lay mainly in the identification of areas of S. alterniflora close to the beach as P. australis. Analyzing the phenological difference between P. australis and S. alterniflora showed that P. australis mostly turned yellow during the senescence period, while S. alterniflora mostly remained green. This characteristic can be used to distinguish P. australis from S. alterniflora; however, S. alterniflora near the coast (at the low tide level) was affected by tidal water stress, and its senescence period was earlier than that of S. alterniflora at the mid-tide level. Therefore, the characteristics of images of S. alterniflora near the coast were close to those of P. australis during the senescence period, which led to some misclassifications. The classification error for the green period was due mainly to misclassification between S. alterniflora and trees. The main reason for these misclassifications was that S. alterniflora and trees in the greening season both appeared as saturated bright green in the image and their VI characteristics were similar. The harmonic model could effectively distinguish between the phenological features of S. alterniflora and trees. However, in some areas, S. alterniflora remained green during the green period, similar to the phenological characteristics of the tree canopy, so there were some partial misclassifications.

The areas of vegetation were further analyzed based on the classification of vegetation outside the Chongming Island dam in different periods in this study.

Table 6. Comparison of vegetation-type areas outside the Chongming Island dam extracted for all seasons, the senescence period, and the green period.

Type	All Seasons (km2)	Senescence Period (km2)	Green Period (km2)
P. australis	55.30	57.47	50.10
S. alterniflora	61.95	52.68	57.81
S. mariqueter	20.51	10.85	18.16
Tree	5.09	8.84	2.55
Cropland	1.56	2.07	1.53
Total	144.41	131.91	130.15

indicates that the all-season classification effectively avoided underestimating the area of vegetation better than the classifications in the green and senescence periods. The areas of P. australis, trees, and cropland during the senescence period were estimated to be higher than those in the all-season classification, while the areas of S. alterniflora and S. mariqueter were estimated to be lower. The total area of vegetation in the green period was estimated to be lower than that in the all-season classification.

3.6. Comparisons with Other Methods

The results from mapping the annual wetland vegetation by using phenological features in this study were compared with the results from other related studies that used phenology to structure feature collections in the study area. Tian et al. [38] recently published a method for mapping S. alterniflora using spectral-phenological characteristics for comparison. They proposed a new method for annual S. alterniflora mapping based on the phenological characteristics of the pixels and compared it with classifications based on a single image and a composite image of pixels in the green period or the senescence period. The results showed that the method based on the phenological composition of the pixels was the best method to reflect the characteristics of S. alterniflora over a year. Therefore, the feature combination method considering phenology proposed in this paper was compared with the optimal pixel-based phenological composition method described by Tian et al. [38] to extract wetland vegetation for a certain year. The two phenological methods were applied to the samples from the study area in this article, and a RF classifier was used to classify the wetland vegetation. The benefits and drawbacks of the compared wetland vegetation mapping methods were then analyzed.

Considering the classification accuracy data in

Table 7. Comparison of wetland vegetation extraction by the feature combination method considering phenology and Tian et al.’s method based on pixel phenological composition.

Accuracy	The Result of Tian et al.	The Result of the Feature Combination
OA	78.26	85.54
P. australis PA (%)	84.91	88.46
P. australis UA (%)	77.59	82.14
S. alterniflora PA (%)	73.85	83.64
S. alterniflora UA (%)	88.89	92.00
S. mariqueter PA (%)	80.00	60.00
S. mariqueter UA (%)	80.00	66.67

, the method of combining phenological features studied in this paper was better able to extract the wetland vegetation types outside the Chongming Island dam than the method of Tian et al. [38]. The classification differences can be observed by comparing Figure 7 and Figure 11. Although the method of Tian et al. [38] was sensitive to S. mariqueter and could distinguish S. mariqueter from P. australis well, the method in this paper resulted in higher OA, UA, and PA for S. alterniflora and P. australis. With regard to the biological invasion of S. alterniflora, the method in this paper could better distinguish S. alterniflora from the native vegetation P. australis and S. mariqueter. Tian’s method was mainly designed on the basis of the phenological period of S. alterniflora and did not consider other vegetation types. The phenological characteristics of S. alterniflora may partially overlap those of other vegetation types, resulting in low discriminations between P. australis and S. alterniflora and flourishing S. alterniflora vegetation may be classified as trees in some areas. Moreover, the differentiation between cropland and P. australis also needs to be improved in Tian’s method.

4. Discussion

4.1. Advantages of the Harmonic Model in Extracting Vegetation Phenological Features

S. alterniflora exhibited a significant difference in VI between the green period and the senescence period according to the study by Zhang et al. [39]. In the appropriate phenological period, S. alterniflora in the tidal flat wetland could be identified well with the VI. However, the phenological periods in this study could only be approximated by phenological analysis because the raw time series data lacked some of the VI information [40,41]. The SG filtering method is commonly used to fill in missing values and improve images for analysis, but it encounters difficulty in reflecting the phenological information about vegetation and does not reveal the periodic characteristics of vegetation growth. To solve these problems and analyze the phenological differences in wetland vegetation on Chongming Island, this paper used the harmonic model to effectively complete and simulate the phenological information in the time series.

The harmonic model not only effectively obtained periodic signals from time series data despite noise interference, but also identified some critical time points in the phenological period such as the beginning of the green period, the beginning of the senescence period, and the duration of the green period. Harmonic models can comprehensively and effectively reflect the phenological characteristics of wetland vegetation during the growth period of a year, and they can also be used in croplands and bamboo forests [42,43]. Li et al. [36] used VIs and harmonic models to extract mangroves and many kinds of herbs and found that harmonic models could extract phenological information for vegetation well and that mangroves and herbs had significantly different phenology. The length of the time series for building a harmonic model in this study was selected as two years. This choice was made because when different time lengths were tried, the harmonic model that was based on a long time length (the longest was from 2014 to 2019) caused greatly increased amplitude and phase variations in the training pixels of the same species. The reasons include both natural and non-natural aspects such as the complex tidal environment, interspecies competition, land cover changes due to human activities and ecological governance of invasive species. After several comparisons, a 2-year time length was chosen in this study. This rationale addresses a trade-off between building a harmonic model and decreasing the complex impacts on the same pixel.

4.2. The Features Combined with Phenological Features for Wetland Vegetation Identification

Although many studies have revealed that during a suitable phenological period, VI data could effectively differentiate wetland vegetation such as mangroves and S. alterniflora [13,17], there are few studies that have conveyed how to apply phenological features to images for classification and investigate their impacts. In this study, the amplitude and phase obtained by the harmonic model were used to represent phenological features and analyze their influence on the classification. Compared with statistical features and temporal features, this approach indicated that phenological features provided the highest accuracy when only one kind of feature was used to classify wetland vegetation, and they contributed the most to improving classification accuracy when combined classification was used.

Furthermore, according to the results of the harmonic analysis and the phenological characteristics of vegetation, wetland vegetation phenology was divided into two periods: the green period and the senescence period. The vegetation canopy characteristics in remote sensing images from the green period and the senescence period were different. Zhang et al. [39] showed that identifying suitable seasonal phases according to the phenological differences in vegetation was conducive to the extraction of different types of vegetation. In that study, analyzing the classifications from the two periods revealed that the different types of wetland vegetation were distinguished more effectively in the senescence period than in the green period. Therefore, performing classification with images from the senescence period could be useful for distinguishing and identifying wetland vegetation. The results of this study are consistent with the above conclusions. Moreover, the classification based on all seasons was superior to that based only on the senescence and green periods. The classification method in this paper was also compared with the method of Tian et al. [38] and found that the OAs for the classification and identification of P. australis and S. alterniflora were better than those obtained with the method of Tian et al. [38]. The reliable annual mapping results obtained from this study are of great significance for the development of long-term annual wetland vegetation mapping.

4.3. A Typical Useful Method for Identifying Wetland Vegetation Based on Google Earth Engine (GEE)

GEE can perform pixel-based calculations, and using high-density GEE remote sensing images for mapping can allow for the full utilization of the high-quality pixels in the image through cloud detection codes, which greatly increase the utilization of data [44]. This study used GEE to flexibly and effectively calculate the features of Landsat images [45] and mapped the wetland vegetation outside the Chongming Island dam. GEE can extract the time series features of wetland vegetation through simple codes and obtain Landsat images for any time and region with the removal of clouds. GEE’s code and methods are easy to share and can effectively reduce work duplication. Moreover, the open-source Python interface can be connected to GEE, which is conducive to applying open-source methods [46]. The GEE platform includes a variety of classification methods. In future research, GEE image data will be used on the cloud platform with deep learning methods to perform classification. However, GEE use also has some limitations such as the limit on the number of user calculations. A calculation number that is too large will exceed the permissible calculation scope in GEE.

5. Conclusions

Wetland vegetation is easily affected by the complexity of the wetland environment, and it is difficult to directly use raw spectral data and VIs to classify wetland vegetation. However, the harmonic model can be used to extract changes in the VI over time and can effectively analyze the phenological changes in wetland vegetation. In this study, the GEE platform was used to analyze three feature types: statistical features, temporal features, and phenological features. The combination of the three features provided the best classification of wetland vegetation, reaching a classification accuracy of 85.54%. Including phenological features in the classification combinations resulted in the most obvious improvement in classification accuracy, with the overall accuracy increasing by 6.32%. The phenological characteristics after harmonic transformation can be used to effectively distinguish wetland vegetation from croplands and trees and to improve discrimination between P. australis and S. alterniflora. In this study, the harmonic model was used to analyze the characteristics of wetland vegetation in two periods. The growth period of wetland vegetation was divided into a green period and a senescence period, and the optimal feature combination method was used to extract the wetland vegetation in the two periods separately. The senescence period classification was more accurate than the green period classification, and the all-season classification was better than that from either period alone. Analyzing changes in vegetation phenology can help to better distinguish wetland vegetation. This study provides an accurate reference for the development of remote sensing inversions in different phenological periods and an accurate method for the annual mapping of wetland vegetation. This phenological method can also be used in other areas to extract other vegetation types that exhibit phenological characteristics.