The study area is in the Yellow River Delta National Nature Reserve, which is located in the northeast of Dongying City, Shandong Province, China (118°32’58’’E-119°20’27’’E, 37°34’46’’N-38°12’18’’N). It belongs to warm temperate zone and semi-humid continental monsoon climate, with four distinct seasons and rainy summers. The annual average temperature is 11.7-12.6°C, the annual average precipitation is 530-630 mm, and about 70% of the precipitation is concentrated in summer. The study area covers the intertidal zones of the Yellow River Estuary ( Figure 1 ), with an area of 923 km². P.australis, S.salsa, and S.alterniflora are the primary vegetation species in the study area (Wang et al., 2021b; Zhang et al., 2021). P.australis generally grows on both sides of the river bank; in the inner part of tidal flat, it is mixed with Tamarix Chinensis. P. australis starts to grow in April, flowering from August to September, and start senescence in October. P.australis near the river generally grows better with higher density, while P.australis in the area with higher salinity is relatively short and sparse. S.salsa is an annual herb with strong salt-tolerance. It is mostly found in mid to high tide areas and covers a wide range. It blooms red from July to October. S.alterniflora is a perennial herb native to the Atlantic coast of North America. It was introduced to the Yellow River Estuary in the 1990s. Due to its strong reproductive capacity and environmental adaptability, S.alterniflora has expanded rapidly in the tidal flat area of Yellow River Delta in recent years, resulting in degradation of the native S.salsa and seagrass bed, which has seriously affected the biodiversity in the coastal wetland. The study area was divided into four zones where Zone A and B are located in the north bank of the estuary, and Zone C and D are located in the south bank. Zone B and C are located near the river mouth.

Landsat 8 satellite is a multispectral imaging satellite launched in 2013. Its carries Operational Land Imager (OLI) sensor with 9 spectral bands from visible to shortwave infrared wavelengths. We downloaded all available Landsat 8 Level 2 Tier 1 surface reflectance images covering the study area (Row 121, Path 43) acquired during January 1, 2020 ~ December 31, 2020 from Google Earth Engine (GEE) platform. The quality assessment (QA) bands of the images were used to identify the area covered by clouds and cloud shadows. There were 18 images in total, and the average cloud coverage was 29.2%. Figure 2 illustrates the spatial distribution of the number of valid observations (no cloud/cloud shadow). The average number of valid observations is 12.7 per pixel, while the number was 11.5 per pixel over the intertidal area.

Auxiliary datasets include high-spatial-resolution images taken by DJI Phantom 4 Multispectral (P4M) Unmanned Aerial Vehicle (UAV) and Gaofen-6 satellite images, which were primarily used for reference data collection. DJI P4M UAV carries a RGB camera and a multispectral sensor with 5 spectral bands including blue, green, red, red-edge and near-infrared (NIR). In September 2020, around 11 UAV flights with an average coverage of 10.2ha were taken in the study area ( Figure 1 ). For each flight, the flight height was 50 m, resulting in 2.65 cm spatial resolution. The along path and cross path overlapping area were over 70%. Within the coverages of UAV flights, 68 field plots with size 1m × 1m were randomly selected. Location of each plot was recorded with handheld GPS RTK equipment, the vegetation species, density and growth status were also recorded at the field surveys. For all UAV images, DJI Terra software was used to generate multispectral orthophoto images. The ortho-images were then segmented into objects using multiresolution segmentation algorithm embedded in eCognition software. Base on the vegetation indexes calculated for each object, the threshold method was used to classify the objects into bare flat, S. salsa, S. alterniflora and P. australis. The classification maps were then upscaled to Landsat 30 m resolution and the fractional cover of each vegetation type within 30 m-grids were calculated using area aggregation approach. The total of 112.3 ha UAV flight coverage resulted in 826 samples.

Because many areas in the YRD wetlands were difficult to access, the UAV flight coverages and the reference samples generated from the UAV images were limited. To supplement the reference samples, high-spatial resolution imagery acquired by Gaofen 6 satellite (GF-6) on September 4, 2020 was used to generate additional reference samples of fractional cover. GF-6 is a high-spatial-resolution satellite that was launched in 2018 as one of the series of China High-resolution Earth Observation System (CHEOS) satellites. It carries a 2-meter resolution panchromatic camera and an 8-meter multi-spectral imager with blue, green, red, and near-infrared band. We first fused the panchromatic imagery with the multispectral imagery using NNDiffuse Pan Sharpening method to obtain 2 meters-resolution multispectral image. Then, we utilized the dimidiate pixel model to estimate the fractional vegetation cover for every 2 m pixel (Song et al., 2022). As the dimidiate pixel model cannot discriminate vegetation species, we only selected those pixels that contain a single vegetation species as reference pixels. Expert knowledge and field experiences helped to determine whether a pixel contain one species. For example, S. alterniflora at the landward edge is unlikely mixed with other species (Zhang et al., 2020). As a result, 348 sample points were generated based on GF-6 images ( Figure 1 ).

To date, many cloud removal algorithms have been developed (Zhu et al., 2012a; Cao et al., 2020; Wang et al., 2022). These algorithms used one or more cloud-free imagery as reference images to predict the missing values in the cloud and cloud shadows in the target image (i.e., the cloud image). However, these algorithms had limited performance when dealing with landscapes with abrupt changes (Wang et al., 2022) such as the estuarian wetlands that are frequently inundated by tidal water.

To solve the above problems, our previous research proposed VICR, a new cloud removal algorithm based on time series reference images. VICR implements cloud removal by filling each cloud region separately. For each cloud region, it consists of three steps: (1) Virtual image construction by linear transformation using time series Landsat imagery. In this step, optimal number of reference images is determined. (2) Similar neighboring pixel selection with assist of a newly proposed temporal-weighted spectral distance. (3) Residual image estimation and cloud image reconstruction by adding residual image to the virtual image. VICR also proposed a strategy for time-series cloud image processing. Details of the model can be found in Wang et al. (2022). Following this strategy, the Landsat imagery acquired in 2020 over the study area ( Table 1 ) were sorted in the order of the cloud cover percentage; the image with the lowest cloud cover was processed first and then the cloud-removed image was used as reference image for images with larger cloud cover.

Temporal information is helpful to distinguish different salt marsh vegetation types and different fractional cover of the same vegetation type (Song et al., 2022). Our previous research found that temporal composite of spectral indices as input features help to discriminate different coastal wetland cover types (Wang et al., 2021b). We also found that the harmonic regression features improved the classification accuracies. Harmonic regression fits the time series spectral indices [such as the normalized difference vegetation index (NDVI)] using superposition of periodic curves and can well represent the phenological pattern of each vegetation species. Following our previous study, we first calculated seven spectral indexes from each Landsat 8 OLI images in 2020 (Wang et al., 2021b), including Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Soil Adjustment Vegetation Index (SAVI), Green Chlorophyll Vegetation Index (GCVI), Green Normalized Difference Vegetation index (GNDVI), Land Surface Water Index (LSWI) and Modified Normalized Difference Water Index (MNDWI). NDVI is the most common vegetation index to reflect the vegetation type and growth status (Tucker, 1979). EVI takes into account the canopy background and aerosol influences, so it is more sensitive to high biomass than NDVI (Huete et al., 2002). Compared to NDVI, SAVI is more suitable for low vegetation cover areas because it adds soil adjustment coefficient (Huete, 1988). GCVI has a larger dynamic range than NDVI and is suitable for densely vegetation areas (Grevstad et al., 2003). GNDVI has significant correlation with chlorophyll content and leaf area index (Gitelson and Merzlyak, 1998). LWSI is sensitive to canopy water content and soil moisture (Xiao et al., 2005), and MNDWI is good at identifying open water (Xu, 2006). The spectral indices were calculated using the following functions:

where ρ_blue ρ_green , ρ_red , ρ_NIR and ρ_SWIR1 are the surface reflectance in blue, green, red, near infrared and short-wave infrared 1 bands in Landsat 8 OLI images.

The annual maximum, minimum, mean, median and standard deviation of the six spectral bands (blue, green, red, NIR, SWIR1 and SWIR2) and the seven spectral indexes were calculated for each pixel based on all Landsat imagery in 2020. Therefore, a total of 39 temporal composite images were generated.

In addition, the Harmonic ANalysis of Time Series (HANTS) method was used for all spectral indexes with obvious periodicity except for MNDWI. This method is beneficial to identify plant phenology, which helps to distinguish different plants (Zhou et al., 2015). The mathematical expression of HANTS used in this study is as follows:

where A is the amplitude of the harmonic wave, which represents the fluctuation range of the spectral index time series curve; the value can reflect the difference in productivity of different vegetation types in the whole cycle. Phase φ represents the peak time of spectral index, i.e., the peak time of vegetation growth. a ₀ is the remainder value of the curve, representing the annual average value of the spectral index. In this study, the amplitude, phase and remainder of six spectral indexes constituted a total of 18 harmonic regression features.

The ELM combined Random Forest Regression (RFR), K-Nearest Neighbor Regression (KNNR) and Gradient Boosted Regression Tree (GBRT). RFR is composed of multiple regression trees based on the bagging algorithm. There is no association with each decision tree in the forest, and the final output of the model is jointly determined by each decision tree. The selection of samples and features in RFR is random, which can effectively reduce the occurrence of over fitting. In addition, RFR can evaluate the importance of different features, has strong processing ability for high-dimensional data, and has a certain anti-noise ability, which makes this method widely used in remote sensing data (Ge et al., 2020; Yang et al., 2020). KNNR is an instance-based machine learning regression model which assumes that similar samples are more proximity in the feature space (Ge et al., 2020). In the process of regression prediction, the value of k neighbors is used as the prediction result. KNNR needs to normalize all features first, and then choose a distance measurement method to calculate the similarity between pixels. In this paper, Euclidean distance was used to calculate the similarity. GBRT is also a regression-tree-based machine learning model. Different from RFR where each regression tree is independent, GBRT connects each tree (weak learner) in a linear combination to continuously reduce the residual errors by the loss function. (Di et al., 2019; Yu et al., 2021). In the training process of GBRT, weak learners are generated through multiple iterations, and each learner is trained according to the residuals of the previous learner. Through iterative improvement of each weak learner, the GBRT model is finally obtained.

The ELM developed in this study integrated three models by using GBRT model. This is because that GBRT model has the following advantages: (1) strong prediction ability for low dimensional data; (2) strong processing ability for nonlinear data; and (3) strong flexibility in handling various continuous values, discrete values, and other types of data. Specifically, the predicted values from each of the RFR, KNNR and GBRT were used as temporal-spectral features, and the same training samples for each individual model were used to train the GBRT model, which was then used to predict the fractional cover of salt marsh vegetation species. Different machine learning models all used the grid search method to determine the optimal parameters.

We aimed to investigate whether the cloud-removed imagery help to enhance the fractional cover estimation of different salt marsh species, and whether the ensemble learning regression algorithm helped to improve the accuracies. For this purpose, we designed four experimental scenarios as follows.

For each scenario, we can compare the estimation accuracies from the original imagery with those from the cloud-removed imagery; we can also compare the accuracies from each of the individual models and that from the ELM. By comparing scenario 1 with scenario 3 and by comparing scenario 2 with scenario 4, we can evaluate whether cloud removal can compensate the unavailability of observations during critical months. By comparing scenario 1 with scenario 2 and by comparing scenario 3 with scenario 4, we can evaluate the role of harmonic regression. Note that harmonic regression is essentially a gap filling algorithm which can build full time series observations, although its purpose is not recovering missing values obscured by cloud/cloud shadow.

For each scenario, ten-fold cross validation was used to evaluate the model performance. Specifically, the model was trained for ten times, at each time the model is fitted by a training data set consisting of randomly selected 90% of the total reference data, and the remaining 10% was used for validation. The accuracy assessment metrics include determination coefficient (R-square), Root Mean Square Error (RMSE) and Mean Absolute Error (MAE), and the formula are as follows:

where y_i represents the reference fractional cover measured by UAV or high-spatial-resolution imagery, y i ¯ represents the mean value of reference fractional cover, and y i ^ represents the predicted fractional cover. R-square represents the reliability of the regression model. Larger R-square indicates higher fitting accuracy. MAE can measure the average absolute difference between the fractional cover estimation and the reference values. RMSE is similar to MAE, but it can amplify larger errors.

Figures 4 – 7 showed the fractional cover estimation accuracies of the four scenarios using RFR, KNNR, GBRT and ELM, respectively. For all three vegetation species, the fractional cover estimation accuracies using the cloud-removed imagery were higher (greater R-square, lower RMSE and MAE) than those using the original imagery regardless of the scenarios and the machine learning models (expect for S.salsa in Scenario 2). Although the fractional cover estimation accuracies were different, all three independent models showed similar patterns as the ELM. The improvements were especially noticeable in Scenario 3 and Scenario 4 when assuming the images in March, July and October were unavailable. For example, for ELM, in Scenario 3 the average R-square increased from 0.859 to 0.922 (RMSE decreased from 8.4% to 6.2%), and in Scenario 4 the average R-square increased from 0.818 to 0.902 (RMSE decreased from 10.1% to 7.2%) when cloud removal was performed ( Figure 7 ). However, when all the original Landsat images were used, good accuracies could be achieved even without cloud removal as long as harmonic regression parameters were added as input features, and the improvement resulted from cloud removal was minimal. For example, for ELM, in Scenario 1 the average R-square was 0.881 when the original Landsat images were used, and the average R-square was 0.891 when all cloud-removed imagery were used ( Figure 7A ). When the harmonic regression parameters were not involved in the fractional cover estimation model (Scenario 2), the accuracies considerably decreased, with R-square of only 0.839 with the original imagery and 0.849 with the cloud-removed imagery. This indicates that the harmonic regression features were more important than removing clouds from the images in discriminating saltmarsh species as well as in discriminating the vegetation cover differences if the time series images were sufficient. When the images in March, July and October were not involved, the fractional cover estimation accuracies decreased significantly even when the harmonic regression features were used, especially for S.salsa and the average accuracies (Scenario 3 vs. Scenario 1 without cloud removal). For the ELM, the R-square of the estimated S.salsa fractional cover declined from 0.854 to 0.794 when images acquired during the three months were not used. In this case, cloud removal for the remaining images improved the accuracies substantially. And the R-square of the estimated S. salsa fractional cover was 0.889 (S3-CR in Figure 7 ), even higher than Scenario 1 (S1-CR in Figure 7A ). In general, cloud removal is helpful to improve the accuracy of fractional cover estimation, especially when there are few good observations.

From Figures 4 - 7 , the ELM generally achieved the best accuracies for all scenarios regardless of using the original imagery or using the cloud-removed imagery. Tables 2 – 4 respectively list the R-squares, RMSEs and MAEs of the estimated fractional cover derived from RFR, KNNR, GBRT and the ELM in Scenario 1 based on the cloud-removed images. The average R-square of the ELM estimation was 0.891, the average RMSE was 7.5% and the average MAE was 2.6%, which was higher than each individual model. Among the three individual models, RFR yielded the highest accuracies, slightly lower than those of the ELM. Compared to KNNR and GBRT, the accuracy was significantly improved when the models were integrated through GBRT, indicating that the GBRT can learn the residuals of each individual model through the integration process and effectively improve the estimation accuracy. For example, the average RMSE of the three-sub models is 8.03%, while the RMSE of the ELM is 7.48% ( Table 4 ). Especially for P.australis, the RMSE of the three sub-models is 2.87%, while the RMSE of the ELM is 8.96%, and the RMSE decreases by an average of 9.97%, which indicating that ELM significantly improved the estimation accuracy of fractional cover of P.australis. Table 2 showed the accuracy of vegetation coverage estimation of P.australis is the highest, followed by S.alterniflora, and finally S.salsa. The average R-square values of their sub-models are 0.905, 0.891, and 0.812 respectively. And for the ELM, the R-square for P.australis, S.alterniflora and S.salsa were 0.924, 0.890 and 0.859 respectively. For S. alterniflora, although the R-square of the ELM was very close to the average R-square of the three models, the ELM has obvious improvement in MAE. This also shows that the integration process can help improve the estimation accuracy.

Figure 8 presents the fractional cover maps of the three salt marsh vegetation types estimated by each individual model and ELM respectively. All four models show generally similar spatial patterns of dense patches of S.alterniflora and P.australis: dense coverage (fractional cover over 0.5) of S.alterniflora was mainly distributed near the river mouth, while dense coverage of P.autralis was distributed along the river bank. However, considerable differences existed in terms of the distributions of low coverage of different vegetation types. From the KNNR model, low density S. alterniflora (fractional cover between 0.1 and 0.4) was widely distributed in the supra tidal zone (Zone D), where S. alterniflora growth is impossible due to high frequency of inundation. Compared to the GBRT model, wide area of P.australis with low density was distributed in the supratidal zone in Zone D, which was also inconsistent with the reality. In addition, a small patch of high-density P.australis (fractional cover over 0.8) was found through GBRT model in the sand bar at the river mouth (Zone C), which is also unlikely to occur. In contrast, the spatial extent estimated by RFR and GBRT was similar, but there were considerable differences in the fractional cover estimations for each vegetation type. For example, the spatial extent of S.alterniflora estimated by RFR model is much smaller than that estimated by GBRT in the south coast. But it is obvious that S.alterniflora should not grow on the sea, and there are some errors in both models. On the tidal flat of the south bank, the spatial extent of P.australis estimated by RFR model is less than that estimated by GBRT, while the estimations for S.salsa by two models are obviously opposite.

Although there are some obvious errors in the fractional cover estimated by each individual models, the estimation accuracy was significantly improved by integrating the estimation results with the GBRT. For example, the over-estimation of S.alterniflora coverage along the south coast was significantly reduced, and the newly formed S.alterniflora patches can still be discovered. The estimated fractional cover of P.australis was also more reasonable. In the middle-low intertidal area with high soil salinity, the over-estimation of P.australis coverage is significantly reduced. Compared to the other two vegetation types, the final estimated fractional cover of S.salsa was low (ranging from 0.05 to 0.3) and the spatial extents of S.salsa was smaller than that estimated from the other models, which was consistent with field investigations and our previous reports (Han et al., 2022b). S.salsa is vulnerable to the tidal influence and the plants are generally sparse, therefore the estimation for S.salsa coverage is relatively difficult. By combining the three models, the fractional cover estimation for S.salsa was more robust. In general, the integration of the three individual models helps to improve the fractional cover estimation accuracy, and the spatial distribution of the estimated fractional cover of the saltmarsh vegetation species by the ELM is more reasonable.

The results from ELM showed that the three species, S. alterniflora, P. australis and S. salsa, covered 5753.97 ha, 4208.4 ha and 1984.41 ha, respectively; and the average fractional cover was 58.45%, 51.64% and 51.64%, respectively. The fractional vegetation cover maps ( Figure 8D ) showed that P.australis was mainly distributed along the river banks and along the Qingshuigou course, the old river channel before 1996 ( Figure 1 ). According to the zonal statistics ( Figure 9A ), the average fractional cover of P.australis was the highest in Zone B, which is 58.26%. The average fractional cover of P.australis in Zone D and in Zone C were 48.93% and 48.47%, respectively. Zone C demonstrated large spatial variation in P.australis coverage ( Figure 9A ). The coverage showed decreasing trend from the river banks to the tidal flats. This is consistent with the existing field-based studies (Xie et al., 2021), which reported that the biomass and coverage of P.australis decreased with increasing soil salinity and decreasing freshwater supply in the tidal flat. The coverage along the old Qingshuigou course is lower than that along the current river course, which is probably due to the insufficient water supply (Wu, 2022). With the expansion of S.alterniflora, the habitat of P.australis was invaded. As a result, the spatial extent area of P.australis was much lower than that of S.alterniflora in Zone C (1152.2 ha vs. 2604.6 ha, Figure 9E ). In Zone D, the area of P.australis was significantly larger than that of S.alterniflora (948.1 ha vs. 664.5 ha).

S.salsa had the smallest spatial extent and lowest fractional cover among the three vegetation types ( Figures 9B, E ). S.salsa was mainly distributed in the mid-high intertidal area, and the average fractional cover was around 12.6%. The dams and groins in Zone A and D blocked the tide waves (Xie et al., 2018), which affected the salinity and moisture content of the intertidal zone. Figures 9E shows the average coverages of S.salsa in zone A and D (13.07% and 10.91%, respectively) were slightly lower than those in zone B and C (14.27% and 13.44%, respectively). In the west part of Zone A, S.salsa was mixed with S.alterniflora, in the landward front of S. alterniflora invasion. In Zone B and Zone C, S. salsa was mixed with P.australis around the river banks.

S.alterniflora generally has the widest spatial extent and densest fractional cover among the three vegetation types. The average fractional cover of S. alterniflora is higher in Zone B and Zone C near the river mouth than those in other zones ( Figure 9C ). Although the area of S.alterniflora in zone C is larger than that in zone B (2604.6 ha vs. 1672.5 ha, Figure 9E ), the average fractional cover of zone B is higher (59.72% in Zone C vs. 70.82% in Zone B). S.alterniflora is also widely distributed along the coast of zone A, and the fractional cover in its west is higher than that in its east, which is associated with less tidal inundation due to higher elevation in Zone A. On the whole, the zonal difference of S.alterniflora coverage is associated with the invasion ages. S.alterniflora was first found in Zone B in 2008, then expanded to Zone A and Zone C, and finally expanded to Zone D in 2017 (Wang et al., 2021b). In addition, it has been reported that the live stem density of S.alterniflora is related to the invasion ages (Han et al., 2022a). Therefore, we calculated the statistics of fractional cover of S.alterniflora with different invasion ages ( Figure 9D ). Figure 9D shows that the coverage of S.alterniflora gradually increased during the first five years of the invasion (average values from 55.04% with 1 year of invasion to 75.42% with 5 years of invasion), and then kept high (around 70%). It is worth mentioning that in other studies (Wang et al., 2021b; Han et al., 2022a), they are also reported that the first five years of invasion is the key period for the expansion of S.alterniflora.

Cloud contamination is inevitable in optical remote sensing, especially for the optical imagery acquired over the cloudy coastal area. To date, many algorithms, such as mNSPI, GNSPI, WLR, ARCC, have been developed for removing clouds/cloud shadows and reconstructing missing images (Zhu et al., 2012a; Zhu et al., 2012b; Zeng et al., 2013; Cao et al., 2020). Compared to the number of algorithms that have been developed, the number of applications are limited. A few studies in recent years have applied cloud removal as preprocessing step for phenological metrics derivation (Tian et al., 2020; Zhu et al., 2021), paddy rice mapping (Zhao et al., 2021) or vegetation cover estimation (Wang et al., 2020). Zhu et al. (2021) built time-series cloud-free Landsat imagery by reconstructing cloud-contaminated imagery using NSPI algorithm and then derived dry-season phenology in tropical forest. Their study found that cloud removal could help better characterize the phenological features. Zhao et al. (2021) applied mNSPI to remove cloud from Landsat imagery, and then extracted phenological features from the time series imagery for paddy rice mapping. They mentioned that the mNSPI could not accurately restore the small and continuous boundaries on the image under the clouds. Wang et al. (2020a) is probably the only study that applied cloud removal algorithm (GNSPI algorithm) to build cloud-free Landsat imagery for green vegetation cover estimation. However, their study did not estimate vegetation cover at species level.

Most of the existing studies applied the cloud removal for constructing time-series imagery, based on which phenological features can be derived. However, cloud removal may not be the necessary step for phenological features retrieval, although few studies have discussed this issue. Time-series vegetation indices can also be reconstructed by fitting and filtering methods, such as harmonic regression (Yan and Roy, 2020) or Savitzky-Golay filtering method (Chen et al., 2021). For example, the harmonic regression utilized in our study produced a simplified continuous time-series curve that can fill the data gaps. Interestingly, our results showed that cloud removal is not necessary in all cases. When the number of Landsat imagery were sufficient and temporal features based on harmonic regression were used, cloud removal did not significantly improve the fractional cover estimation accuracy (Scenario 1) ( Figures 4A–C ). In this case, it seemed that harmonic regression played more important roles than cloud removal in fractional cover estimation, as the accuracies decreased significantly without the temporal features derived by harmonic regression. However, when the Landsat observations during the critical months were not incorporated, cloud removal for the remained imagery was very important (Scenario 3), while harmonic regression did not help to improve the accuracy. Different from green vegetation cover estimation, the fractional cover estimation at species level not only needs to build the relationship between the temporal features and the fractional cover, but also needs to discriminate among the species. Our previous research showed that the imagery in the key months was critical to represent the phenological patterns of each species and to ensure the discrimination accuracy (Wang et al., 2021b). P. autstralis starts to grow in April, reaches the maximum greenness during July and August, and enters senescence in September. S.alterniflora starts to grow in late May and early June, reaches the maximum greenness during August and September, and then enters senescence in late October to early November (Han et al., 2022a). S. salsa presents red-purple color and the abundance reached the maximum during October. When the images during the critical months are absent, harmonic regression cannot represent the correct phenological patterns. In this case, the remaining cloud-removed images provides important supplementary information. Surprisingly, we found that the utilization of all available cloud-removed imagery produced even slightly lower accuracies than that excluding the critical months. Detailed examination showed that the reconstructed imagery in October had poor visual effects because the cloud covered almost the entire the saltmarsh extent in the estuary ( Figure S1 ), which is quite challenging for all cloud-removal algorithms. This also indicates that when applying cloud removal as preprocessing step, cloud coverage and cloud-removal accuracies need to be considered.

Previous studies have confirmed that machine learning algorithms have great potential in green vegetation cover estimation. For example, Wang et al. (2018) reported high accuracy (RMSE=8.5%) of RFR for green vegetation cover estimation based on Sentinel-2 imagery. Yang et al. (2020) used RF soft classification method to estimate fractional abundance of halophytic species based on high-spatial resolution WorldView-2 imagery over Venice lagoon, Italy, and also achieved high accuracy (RMSE ranging from 0.06 to 0.19). However, the integration of multiple machine learning models into one ensemble model has not been introduced into vegetation cover estimation, especially at species level. Our results showed that the performance of the model can be ranked as the following order: ELM > RFR > GBRT > KNNR.

However, it was also found that the accuracy improvement of the ELM, which was measured by R-square and RMSE seemed not considerable compared to RFR. For example, the increase in the average R-square was only 0.002. Note that the accuracy assessment metrics were calculated based on reference samples, whose spatial locations might influence the evaluation. And when we looked at the spatial distribution of high to low coverage of each salt marsh species, the ELM apparently yielded more reasonable results. Some examples are shown in Figure 10 , illustrating the results in zoomed-in area in Figure 8 . Compared to ELM, RFR overestimated the spatial extent of S. alterniflora in the bare tidal flat close to the sea, and some S. alterniflora even appeared in the seawater (first column in Figure 10 ). In addition, it was unlikely that S. salsa grew on the sand bar near the river mouth (third column in Figure 10 ). Although both RFR and ELM over-estimated S. salsa cover in this area, ELM generally produced lower error than RFR. Previous research in other fields also reported that similar RMSEs or R-squares from different methods does not necessarily mean similar performance in every location (Di et al., 2019; Requia et al., 2020). Di et al. (2019) applied the ELM to estimation PM2.5 concentration across the contiguous United States. They found that each individual model did not perform equally well in every location or at all PM2.5 concentration levels although the overall R-squares are similar; however, the ensemble model complemented each other and produced more spatially balanced results. By integrating individual models in a non-linear manner, the model that performs better at some locations contribute more to the ensemble model, which improves the overall performance of the ELM.

The performance of the ensemble-learning-based fractional cover estimation depends on at least three factors: (1) whether the reference samples can represent the reality of fractional cover (2) whether the predictor indicators (temporal composite features used in our study) can be associated with the variability of response variable (fractional cover in our study), and (3) whether the model can capture the relationship between predictor indicators and the fractional cover of each vegetation species. Our study attempts to improve (2) and (3) by using time-series cloud-removed imagery and by developing ELM, respectively. As field surveys in coastal wetlands are difficult, in this study we relied on UAV images and GF-6 high-spatial-resolution imagery to create reference samples (Di et al., 2019; Yang et al., 2020; Song et al., 2022). Although high-spatial-resolution imagery has been widely used for reference sample collection for fractional cover estimation model, uncertainties still existed. First, the spatial coverage of each UAV flight was very small compared to the whole study area, and there might be little difference in the vegetation fractional cover over each UAV image. We need to conduct many UAV flights to generate enough reference samples, which is time and labor consuming. Second, using GF-6 and other high spatial resolution satellites would still have a certain mixed pixel effect. In the future, more efforts need to be taken to overcome the problem of sample selection. Future research will explore deep learning models for sample augmentation for the fractional cover estimation. In addition, in YRD several wetland restoration projects are being implemented in recent years. Continued monitoring of the fractional cover change in the coastal wetlands are necessary for better evaluating the effectiveness of the restoration.