Development of spectral-phenological features for deep learning to understand Spartina alterniflora invasion

be (2) prominent spectral exist within the ﬂ ora due to phonological variations; (3) spectral between native species is often presented in the territories where S. alterni ﬂ ora intruded in. To articulate these questions, we proposed a new pixel-based phenological feature composite method (Ppf-CM) based on Google Earth Engine. The Ppf-CM method was brainstormed to battle the aforementioned three hurdles as the basic unit for extracting phonological feature is individual pixel in lieu of an entire image scene. With the Ppf-CM-derived phenological feature as inputs, we took a step further to investigate the performance of the latest deep learning method as opposed to that of the conventional support vector machine (SVM); Lastly, we strive to understand how S. alterni ﬂ ora has changed its spatial distribution in the Beibu Gulf of China from 1995 to 2017. As a result, we found (1) the developed Ppf-CM method can mitigate the phonological variation and augment the spectral separability between S. alterni ﬂ ora and the background species regardless of the signi ﬁ cant cloud coverage in the study area; (2) deep learning, compared to SVM, presented better potentials for incorporating the new phenological features generated from the Ppf-CM method; and (3) for the ﬁ rst time, we discovered a S. alterni ﬂ ora invasion outbreak occurred during 1996 – 2001.


Background
Invasive plants pose serious threats to biodiversity and environmental health by interfering with indigenous vegetation community, animal habitats, soil properties, water quality, and biogeochemical cycles (Dukes and Mooney, 2004;Hartig et al., 2002;Vitousek et al., 1997). Spartina alterniflora (S. alterniflora) is native plant species in the With the elevated human activities in coastal zones, the harms of S. alterniflora are becoming more and more obvious and aggravated . However, due to various logistical reasons, limited knowledge regarding the S. alterniflora invasion has been gained as of today (Liu et al., 2018). Consequently, no much attention from national and local conservation agencies has been paid to S. alterniflora in China. In order to develop cost-effective management strategies, monitoring the dynamics of S. alterniflora is urgently needed.
Given the fact that S. alterniflora often disperse quickly and occupy spatially-fragmented areas, it is challenging to use field methods to keep track of its invasion (Ren et al., 2019). Alternatively, remote sensing offers unique capability to monitor S. alterniflora over large areas and long time periods by acquiring the multi-temporal image acquisition from remote platforms such as satellites. In general, the following two groups of methods have been developed with regards to the intrinsic spatial resolution of the images that was adopted: high spatial resolution (finer than 10 m) and moderate spatial resolution (10-500 m).
Regarding the first category of methods (high resolution), different attempts have been made such as airborne multispectral/hyperspectral image (Beland et al., 2016;Chust et al., 2008;Hladik et al., 2013;Wang et al., 2007), SPOT-5 and Gaofen-1 , ZiYuan-3 and Gaofen-1 (Tao et al., 2017), unmanned aerial vehicle (UAV) image (Wan et al., 2014), and SPOT-6 . Regardless of various success of using such high spatial resolution images, the indispensable obstacle for applying such methods to understanding regional scales invasion are two facets: (1) limited spatial-temporal coverage; (2) high costs. Alternatively, the second category of method (moderate resolution), owing to its superior data accessibility, has greater potential to articulate the regional scale invasion problem. Although a handful studies have been carried out (Ai et al., 2016;He et al., 2010;Lu and Zhang, 2013;Ouyang et al., 2013;Ramsey et al., 2014;Wang et al., 2015), their mapping accuracy is still not satisfactory due to two major hurdles: (1) the costal zones where S. alterniflora usually intruded is inevitably annoyed by the significant coverage of cloud, which subsequently reduces the number of images that can be used; (2) poor spectral separability between S. alterniflora and its codominant native species is presented. Consequently, how to augment such spectral separability to refine the mapping accuracy is an immediate question.

Phenological features
Incorporating phenological features derived from multi-temporal remote sensing has found to be effective for augmenting the spectral separability between invasive species and other co-existent native species (Bradley, 2014;Diao and Wang, 2018;Evangelista et al., 2009;Follstad Shah et al., 2007;Hansen et al., 2014;Hufkens et al., 2012;Peterson, 2007;Rapinel et al., 2015;Singh and Glenn, 2009). To this end, two different strategies for incorporating phenological features have been reported in the past.
(1) Scene-based methods: i.e. the basic unit for compositing phenological features is entire scene given the premise the scene is mostly cloud-free. Therefore, often time a limited number of scenes can be selected from all the avaiable multi-temporal images covering the study site. Those scenes not qualified are mainly because their cloud coverage is over a certain pre-defined threshold. As the latest development in the scene-based methods, Diao and Wang (2016) proposed an intra-annual phenological trajectory method drawn from time series of monthly Landsat images. They found that utilization of time-series of intra-annual data is more robust than a single image acquired at a particular phenological stage, e.g. senescence. Although the scene-based methods managed to succeed in some study sites, there are still two major obstacles: (1) spatial variability of phenology can be quite remarkable for the same invasive species within a scene (Fisher et al., 2006;Friedman et al., 2005;Ji and Wang, 2016). For example, Ji and Wang (2016) showed that significant phenological differences can be observed among the same invasive species pixels even within spatial vicinity. This leads to the first obstacle for the scene-based methods: spatial variability of phenology has to be accounted for before a meaningful scene-based phenological feature can be developed; (2) the second obstacle is scene-based methods, because of its basic unit being the entire scene, tend to yield limited cloud-free images. This is particularly true for study area in the coastal zone, where cloud coverage is phenomenal for the year long. Therefore, it is often impossible to derive useful phenological features in coastal zones with insufficiently selected image scenes.
(2) Pixel-based methods. Differed from scene-based method, the second strategy to incorporate phenological feature is named as pixelbased methods. In this case, the basic unit for compositing phenological features is an individual pixel given the premise that the pixel is mostly cloud-free among all the scenes available in a study site (Beckschäfer, 2017). In contrast to the scene-based methods that the entire scene will be bailed out should the image contains significant cloud coverage, cloud-free pixels within cloud-contaminated scenes can be still selected based on which meaningful phenological features for invasive species can be extracted. Consequently, pixel-based methods provide a viable solution to solve the two aforementioned obstacles associated with the scene-based methods (Griffiths et al., 2013;Roy et al., 2009;Lück and van Niekerk, 2016): (1) pixel-based methods can mitigate the spatial variability problems as the basic unit to account for phenology is individual pixel; (2) phenology of invasive species can be better characterized regardless of cloud-coverage at the scene level. This is because all cloud-free pixels in time-series images are employed. Despite the conceptual merit over the scene-based methods, pixel-based methods confront a practical challenge, i.e. one scene is literally composed of enormous number of pixels. Thus, selecting viable pixels from this huge pool of candidates over a long temporal dimension demands massive data processing and unprecedented computing resources, which cannot be afforded through traditional stand-alone remote sensing facility.
Google Earth Engine (GEE), as an emerging public high-performance cloud-based platform, opens up a new revenue to overcome this challenge (Gorelick et al., 2017). Recent studies have shown its potential for advancing our scientific understanding of various dynamic processes associated with earth systems (Hansen et al., 2013;Healey et al., 2018;Huang et al., 2017;Liu et al., 2018;Nijland et al., 2019;Schug et al., 2018). To our knowledge, pixel-based methods to incorporate phenology information with GEE are still limited (Azzari and Lobell, 2017;Massey et al., 2018). Towards this line, Azzari and Lobell (2017) proposed a season-based pixel composite method in GEE, which selected the median values of each Landsat image band and vegetation indices for each season during a four-year span. Although the method achieved a promising accuracy for cropland monitoring, it will have a hard time to be adapted for the study of invasive species, such as S. alterniflora attributed to the following two rationales: (1) the assumption of no areal changes within four years is not valid for invasive species; (2) seasonal difference is not fine enough in the temporal dimension for capturing invasive species.

Classifier
Although the pixel-based phenological feature may augment the spectral separability between invasive and native species, another contributing factor affecting S. alterniflora mapping is the development of an effective classifier. Currently, most mainstream classifiers in remote sensing applications are based on machine learning (Mather and Paul, 2009). Artificial intelligence (AI) was initially proposed in 1956 for the purpose of constructing a machine to mimic human thinking. Machine learning emerged in the 1980s to implement AI algorithms. Deep learning, as a new branch of machine learning, has achieved successful applications in various fields since 2012, e.g. computer vision, speech recognition, and remote sensing classification (LeCun et al., 2015;Zhang et al., 2016). A growing number of studies to date (e.g. Heydari and Mountrakis, 2018;Ishii et al., 2015;Li et al., 2016) have been reported on classification of remote sensing images at moderate spatial resolution with deep learning. For example, deep learning method (Stacked AutoEncoder, SAE) and traditional machine learning methods (e.g. Support Vector Machine, SVM) were compared to map African land-cover (Li et al., 2016). It was reported that SAE achieved a higher classification accuracy than SVM. However, in one of the latest developments, Heydari and Mountrakis (2018) conducted an experiment on land cover classification using the deep learning method (SAE) with 26 Landsat images. Surprisingly, the results indicated that the classification accuracy of deep learning was not as good as traditional machine learning methods, e.g. SVM. Therefore, it is still a wellworthy topic to investigate which kind of methods is better, deep learning or other alternatives. Nevertheless, one hypothesis we aim to test in this study is that whether inclusion of domain knowledge (e.g. phenological features) may play a role in the performance of deep learning or other machine learning methods for the detection of invasive S. alterniflora.

Objectives
The objectives of this study were to: (1) develop a new pixel-based phenological feature composite method (Ppf-CM) based on GEE with aims to minimize the cloud effects in the coastal zone; (2) investigate if the performance of the deep learning method can be enhanced with the proposed new pixel-based phenological feature; (3) understand how S. alterniflora has changed its spatial distribution in the Beibu Gulf of China from 1995 to 2017. As a whole, we expect that Ppf-CM can not only account for the spatial variability of phenology and thus enhance the spectral separability in S. alterniflora detection, but also maximize the use of the available Landsat data in the frequently-cloud-covered coastal zone. Moreover, we aim to investigate the full potentials in phenological knowledge-driven deep learning methods. Finally, the expected long-term and regional-scale S. alterniflora maps will serve as the intrinsic baseline data for the assessment, management, and conservation of ecosystem services in coastal zones, including biological invasion, biodiversity richness, climate change, and carbon storage.

Study area
The Beibu Gulf, located in the northwest of South China Sea, is a semi-enclosed strait surrounded by Hainan island, Leizhou Peninsula, Guangxi Zhuang Autonomous Region and Vietnam (18°47′N to 21°56′N and 107°49′E to 110°26′E). Due to the typical southern subtropical monsoon climate, the Beibu Gulf has a climate characterized by frequent cloudy and rainy days with a mean annual temperature of 22°C and precipitation of 1600 mm. The Beibu Gulf is also characterized by a mean tidal range of 2.24 m (Jia et al., 2015). The dominant vegetation species in the intertidal zone of the Beibu Gulf consist mainly of S. alterniflora and four mangrove species: Avicennia marina, Aegiceras corniculatum, Rhizophora stylosa, and Kandelia candel (Tian et al., 2017).

Landsat imagery
All the USGS Surface Reflectance (SR) images from the Landsat-5 TM and Landsat-8 OLI covering the study site at the spatial resolution of 30 m and temporal resolution of 16-days in GEE from January 1, 1995 to December 31, 2017 were employed as input data. Through carefully scrutinizing field observations and very high spatial resolution (VHR) images, we found that the expansion of S. alterniflora within Frames 3 and 4 ( Fig. 1) was negligible. Therefore, only the Landsat images located within Frames 1 and 2 ( Fig. 1) were used to map and monitor S. alterniflora (420 Landsat-5 TM and 105 Landsat-8 OLI, 19.96 billion pixels). For each Landsat image, six spectral features/bands were employed, including blue, green, red, near infrared (NIR), short-wave infrared 1 (SWIR 1), and short-wave infrared 2 (SWIR 2).
The number of cloud-free and low-tide Landsat images really in this study site is quite limited, as can be seen from the GEE app that we developed for exploring Landsat-8 OLI images: https://kingyan988. users.earthengine.app/view/myapp. We consider an image to be "lowtide" if all the intertidal vegetation (e.g. mangroves and S. alterniflora) were not submerged by water (Jia et al., 2019). Taking Frame 2 (see Fig. 1) in 2016-2017 as an example, only one cloud-free and low-tide Landsat image (acquired on September 16, 2016) was found. This situation became even worse in some other years, especially when taking Frame 1 into account.
To alleviate the problem caused by cloud and tide, Ppf-CM used each cloud-free and low-tide pixel selected from all Landsat scenes during a two-year period to generate a composite image in each specific mapping year (see first column in Table 1), assuming that the extent of S. alterniflora remained unchanged during this period. For example, to construct the composite image in the mapping year of 2016-2017, we will leverage all the Landsat images from January 1, 2016 to December 31, 2017 in the proposed Ppf-CM. According to our extensive field studies, we selected two years as a trade-off between building a cloudfree and low-tide Landsat composite image and detecting the coverage change of S. alterniflora.

Coastal boundary zone
As S. alterniflora is typically distributed in the intertidal areas, a coastal boundary zone covering those areas was digitized in GEE by visually interpreting Landsat images from 1994 to 1995. Moreover, the boundary was expanded by a 30-m radius buffer (the size of one Landsat pixel) to make sure that S. alterniflora were all included. It should be noted that the intertidal area has been getting smaller from 1995 to 2017 due to the expansion of coastal constructions. Therefore, the above-mentioned boundary zone was applied to the other mapping years (see the first column of Table 1) as well. All data analyses described in this study were for this area. For each year, the intertidal zone contains about 443,616 Landsat pixels. This zone is attached in the GEE code (see Section 3.1).

Reference data
To collect the reference data of the two mapping classes (S. alterniflora and Non-S. alterniflora), we conducted field surveys in summer and winter of each year from 2015 to 2018 and visually interpreted VHR images of multiple sources ranging from 2006 to 2017. Non-S. alterniflora class mainly contained mangroves, mudflats, water, and few others (e.g. impervious surface). The VHR image data included UAV images, SPOT and WorldView series images, and Google Earth VHR data, all of which have a spatial resolution finer than 10 m. As shown in Table 1, a rather comprehensive VHR data set was employed to construct the reference data set to help us identify the S. alterniflora pixels in the Landsat images from 2006 to 2017. Particularly, the coverage proportion of VHR data reached 76% in the mapping year of 2016-2017. The number of available VHR images, however, was very limited before 2003. For those mapping years (1994-1995, 1996-1997, 1998-1999, 2000-2001, and 2002-2003), besides two scenes of SPOT 4 images, we also visually interpreted the Landsat images to collect the reference data.
The reference data were randomly sampled following the principles that they were as evenly distributed as possible throughout the study area and were constructed in such a way that the number of S. alterniflora and Non-S. alterniflora (mangroves, mudflats, water, and others) pixels were comparable. Moreover, we make sure that the reference data covered at least 2.2% (about 10,000 Landsat pixels) of the whole study area for any given year. In each year, the reference data were split for tranining (about 60%, or 6000 Landsat pixels) and accuracy assessment (about 40%, or 4000 Landsat pixels).

Methods
We analyzed the S. alterniflora invasion through three steps. First, we constructed composite images of different phenological stages using the proposed Ppf-CM. For comparison purposes, we also constructed five other types of Landsat images. More details can be found in Section 3.1. In addition, to assess the effectivity of the proposed Ppf-CM, we compared Ppf-CM image with the other five types of images by comparing their intra-class and inter-class variance. More details can be found in Section 3.2. Second, we input each types of image to the classifier to map S. alterniflora in 10 different mapping years. Two classifiers SAE and SVM are employed and compared. More details can be found in Section 3.3. Third, the area of S. alterniflora derived from the classifiers are compared across the 10 mapping years to analyze its change over time.

Pixel-based phenological features composite method
As shown in Fig. 2, our proposed Ppf-CM can be divided into three stages: First, the cloud, cloud shadow, and water pixels in the Landsat Surface Reflectance (SR) images were masked by the C Function of Mask (CFMask) method (see Section 3.1.1). Second, two key phenological periods (green and senescence) of S. alterniflora were identified by analyzing the pixel-level Normalized Difference Vegetation Index (NDVI) temporal profiles (see Section 3.1.2). Third, a green composite image in the green period and a senescence composite image in the senescence period were constructed, each containing 6 spectral bands. The two composite images were then combined to yield a new image (Ppf-CM based image, 12 bands) (see Section 3.1.3). The method can be reproduced with the following GEE code URL in the mapping year of  Table 1 Landsat images, high spatial resolution data and field survey description during 1995-2017 over Frame 1 and 2 in the Beibu Gulf. Coverage proportion is defined as the total area of high spatial resolution data and field survey data divided by the whole study site area. The preprocessing of Landsat SR images includes two major steps: first, all the cloud, cloud shadow, and water pixels were masked by the "pixel_qa" band provided in the Landsat SR images, by excluding pixels with flags in "pixel_qa" bits 2, 3, and 5. The "pixel_qa" band contains pixel quality attributes generated from the CFMask algorithm (Foga et al., 2017), with each bit representing one type of pixel quality (e.g. bit 2: water, bit 3: cloud shadow, bit 5: cloud). Second, considering that CFMask cannot filter all cloud and cloud shadow pixels, especially when the cloud coverage proportion is high, we also excluded Landsat images where cloud and cloud shadow (pixel_qa bit 3 and 5 flagged) covered > 70% of the image within the specific coastal boundary zone. The threshold was set to 70% in this work following the previous studies (Griffiths et al., 2013;Hermosilla et al., 2016;White et al., 2014).

S. alterniflora phenology analysis
To capture the optimal intra-annual phenological windows for S. alterniflora detection, an average annual NDVI temporal profile was generated to determine phenological features using all Landsat-8 OLI images from January 1, 2014 to December 31, 2017 (Fig. 3), where NDVI was computed by the formula: NDVI = (NIR -red)/(NIR + red) (Tucker, 1979). As shown in Fig. 3, on each date the NDVI distribution were calculated at the pixel level with 500 pure S. alterniflora Landsat pixels evenly distributed throughout the whole study area, in which the pixels were selected by using VHR remote sensing data visual interpretation method and field measurement boundaries. The annual NDVI temporal profile of S. alterniflora was similar to the results by (Ai et al., 2017), thus confirming the reliability of our Landsat-based phenological characteristics. Fig. 3 shows that the mean NDVI values before day of year (DoY) 132 were all smaller than 0.3, whereas the mean NDVI values between DoY 164 and 260 were much higher than 0.4. Therefore, we can identify two clearly different phenological stages of S. alterniflora, i.e. the senescence period from DoY 1 to 142 and the green period from DoY 154 to 270. It should be noted that we expanded the time range of the green period and the senescence period because the dates may vary across different years during 1995-2011.
Through the observation of field measurements, it was found that the spatial variability of phenology did exist in S. alterniflora, e.g. prominent phenological differences can be observed for pixels belonging to the same invasive species across a specific-date image. For example, as shown in Fig. 4, two stands of S. alterniflora were at rather different phenological stages on the same date, associated with clearly different spectral features. S. alterniflora distributed in Fig. 4(a) was much greener than that of Fig. 4(b).

Green and senescence images composite method
From the two key phenological periods identified in Section 3.1.2, i.e. green and senescence periods, we first generated two composite Landsat images for each mapping year, namely one green composite image and one senescence composite image. Through those composite processes, the spatial variation in plant phenology can be smoothed. The green composite image with six spectral bands was constructed when Enhanced Vegetation Index (EVI) (Van Leeuwen and Huete, 1999) reached the maximum on a per-pixel basis during the green period. Similarly, the senescence composite image was constructed when Plant Senescence Reflectance Index (PSRI) reached maximum during the senescence period. PSRI, ranging from −1 to 1, is a vegetation index which is sensitive to senescence-induced reflectance changes (Merzlyak et al., 1999). An increase in PSRI represents the increasing senescence of vegetation, where PSRI was calculated using the formula: PSRI = (ρ 680nm − ρ 500nm ) / ρ 750nm , in which ρ stands for the spectral reflectance. Replacing the ρ 680nm , ρ 500nm , and ρ 750nm respectively with red, blue, and NIR bands in Landsat, the formula of PSRI becomes: PSRI = (red -blue)/NIR in this work. Although the red, blue, and NIR bands of Landsat imagery have wider spectral band widths, the reflectance of S. alterniflora in the three Landsat bands resembles those at the three specific wavelengths (Fig. 5), which suggested that the Besides the Ppf-CM based image, we also constructed other relevant Landsat images (i.e. three scene-based and two pixel-based images) through leveraging the phenological information for comparison purposes (Table 2). These various scenarios represent the common means to capture phenological characteristics in remote identification of plant species. The definition of each constructed image can be found in Table 2.

Intra-class variability and inter-class separability evaluation
Intra-class variability and inter-class spectral separability are the key factors in determining the detection accuracy of S. alterniflora. A comparative experiment was conducted between Ppf-CM and five other types of Landsat images, including SI, GI, GSI, SCI, and GCI (see Table 2). The six types of Landsat images were constructed in 2005-2006, 2007-2008, 2010-2011, 2013-2014, and 2016-2017 (see the top five rows of Table 1). Those mapping years were chosen for conducting the comparative experiment because the coverage proportion of the reference data were high (> 31%). Parts of the testing data (about 50%) were chosen by using a random selection method ("ran-domColumn" function provided in GEE) for evaluation, including S. alterniflora and Non-S. alterniflora (mangroves and mudflats). Mangroves and mudflats in the Non-S. alterniflora class were chosen because they are difficult to distinguish from the S. alterniflora and they occupy most of the area besides S. alterniflora within the coastal boundary zone.
The Standard Deviation (SD) of S. alterniflora, mangrove, and mudflat were the first statistic to assess its intra-class variability. Smaller SD means smaller intra-class variance, which is better for classification (Dalponte et al., 2009;Shao et al., 2016). Second, the Jeffries-Matusita (J-M) distance, as a widely recognized spectral separability metric (Tolpekin and Stein, 2009), was employed to quantify the respective spectral separability between S. alterniflora and mangrove, mudflat. J-M distance is asymptotic to the value of 2.0 with increasing spectral separability, with larger J-M distance representing larger inter-class distance between S. alterniflora and others, which is better for classification (Schmidt and Skidmore, 2003). Both SD and J-M values were calculated by the ENVI 5.1.

Classifier selection and parameterization
Besides the procedure applied to generate the composite image, the classifier selection is another key factor in determining the detection accuracy of S. alterniflora. As a deep learning model, SAE was chosen as the classifier in this work because it has been successfully applied on Landsat image classification in several studies (Heydari and Mountrakis, 2018;Li et al., 2016). In addition, SVM, a traditional machine learning method, was also selected considering the comparison results of different classifiers in He et al. (2015) and Heydari and Mountrakis (2018). The detailed description of the two algorithms can be found in literature, for example (Chen et al., 2014;Vincent et al.,   The parameter setting method of SVM in this study was done according to the sensitivity analysis method in (Heydari and Mountrakis, 2018). The "kernelType" was set to "RBF" (Gaussian function). For the "gamma" parameter in the RBF kernel, we tested values of 0.01, 0.1, 0.5, 1, 2, 5, 10, 25, 50, 100, and 300. For the "cost" parameter in SVM, we tested values of 0.1, 0.5, 1, 2, 5, 10, 25, and 50. We randomly selected 50% reference data for parameter optimization. By selecting the highest overall accuracy achieved from different parameter combinations, this study set "gamma" and "cost" to 2 and 10 respectively.
The SAE parameter setting for classification of multi-spectral images was described by several authors, e.g. (Heydari and Mountrakis, 2018), but the room for improvements remains significant, especially for the selection of the optimal number of iterations. For example, Heydari and Mountrakis (2018) set the iteration number 100 without further analysis. It should be noted that 100 iterations may not be sufficient to obtain a stable classification, at least in this work. Therefore, we set a large maximum number of iterations (10,000) with the optimal number of iterations determined by the criteria of model performance. We applied two criteria to estimate the optimal number of iterations, i.e. the mean error and the number of misclassified pixels. In the training data set (see Section 2.2.3), we labeled S. alterniflora and Non-S. alterniflora pixels with 1 and 0, respectively. The mean error is the mean difference between predicted values derived from SAE and labeled values. In contrast, the numbers of misclassified points are the sum of S. alterniflora misclassified as Non-S. alterniflora and Non-S. alterniflora misclassified as S. alterniflora. The optimal number of iterations should yield both a low mean error and few misclassified pixels. In the example (the mapping year in 2016-2017) shown here (Fig. 6), the optimal number of iteration should be located in the red circle, which has almost the smallest mean error and smallest number of misclassified pixels simultaneously. Other SAE key parameter values or range applied in this study are given in Table 3. Fig. 7 shows the SD of S. alterniflora reflectance. From the singledate images (GI and SI) to the composite images (GCI and SCI), the SD in each spectral band decreased. This was true in each mapping year. This indicates that the bands in composite images have successfully mitigated the intra-class variability in spectral reflectance. This comparison was also applied to the mangroves and mudflats (see Table S1), which achieved consistent results with S. alterniflora. Fig. 8 shows spectral separability statistics between S. alterniflora and mangrove, mudflat for the six types of Landsat images. For S. alterniflora and mangroves, both GI and GCI acquired in the green periods yielded consistently low J-M values (ranging from 0.93 to 1.12), whereas both SI and SCI acquired in the senescence periods got higher J-M values (from 1.18 to 1.37) than GI and GCI. Compared to statistics in the green periods (GI and GCI), the increased spectral separability using senescence images (SI and SCI) indicated that the senescence periods are more suitable for distinguishing S. alterniflora from mangroves. Conversely, for S. alterniflora and mudflats, the GI and GCI achieving a higher J-M values than SI and SCI indicated that the green periods are more suitable for distinguishing S. alterniflora from mudflats. Upon further comparison, both GSI and Ppf-CM achieved much higher J-M values (ranging from 1.46 to 1.81) than the other four types of Landsat images. This confirms that the combination of two images acquired at the two distinctive phenological stages increases the interclass distance, hence enhancing the spectral separability between S. alterniflora and mangroves, mudflats. As a whole, Ppf-CM achieved the Fig. 5. The mean reflectance of S. alterniflora in green period (May 31, 2017) and senescence period (January 15, 2018). Each spectral line was the mean value of ten spectral curves. The blue, red, and NIR band widths of Landsat 8 are 450-515 nm, 630-680 nm, and 845-855 nm, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) highest J-M value (ranging from 1.69 to 1.81) among the six types of Landsat images.

Classification performance assessment
A comparison experiment was carried out to evaluate the performance of SAE and SVM, and to choose an optimal strategy to map S. alterniflora. The six types of Landsat images (Table 2) in the ten different mapping years (Table 1) served as input of these two classifiers. The results were evaluated with four different metrics: Producer's Accuracy (PA), User's Accuracy (UA), Overall Accuracy (OA), and Kappa Statistic (KS). The classification performance was rather similar for each type of Landsat images in different mapping years from 1995 to 2017. The spatial coverage of the reference data was much higher in 2016-2017 than that in the other years. Accordingly, we focused our analysis on the results obtained for 2016-2017 (Table 4). The OA and KS statistic results in the remaining nine mapping years were shown in Table 5 and S2. The comparison results show that the classification accuracy was improved in three aspects.
(3) Integrating phenological features of the two distinctive stages improved the performance of SAE more than that of SVM. As shown in Table 4, although SAE and SVM performed similarly when using data from only one phenological stage (i.e. GI, GCI, SI, and SCI), SAE achieved higher accuracy than SVM when combining data from two phenological stages (i.e. GSI and Ppf-CM). Taking the mapping year 2016-2017 as an example, for SVM, the OA and KS increased by 4.28% 11.57% and 0.12-0.31, respectively; for SAE, the OA and KS increased by 5.33%~12.92% and 0.15-0.35, respectively.
Besides the mapping year in 2016-2017, the same conclusion can be drawn in the other mapping years as well (Table S2). Moreover, Table 4 also shows that SAE with Ppf-CM based composite image achieved the highest accuracy among all classification procedures, where the values of OA and KS were 96.22% and 0.90, respectively. Moreover, the OA and KA of the nine other mapping years achieved consistently higher values for SAE with Ppf-CM (Table 5 and S2). Therefore, SAE with Ppf-CM will be applied as the optimal procedure to derive S. alterniflora maps in the Beibu Gulf in the ten mapping years (see the first row in Table 1).

Multi-year spatial dynamics of S. alterniflora
The distribution of S. alterniflora in the Beibu Gulf was mapped with Ppf-CM based SAE in the ten mapping years during 1995-2017. Keyhole Markup Language (KML) format files of these mapping results Fig. 6. Mean error and number of misclassified points to estimate the optimal number of iterations in SAE, in which the number in colour bar stands for the density of points.

Table 3
The key parameter values/range setting of SAE classifier.

Parameter
Parameter value/ range Parameter in this work # of nodes in 1st hidden layer 6 to 60 (step of 6) 36 # of nodes in 2nd hidden layer 0 to 60 (step of 6) 36 # of nodes in 3rd hidden layer 0 to 60 (step of 6) 0 Batch size 10 to 100 (step of 10) 10 # of iterations 10,000 7132-9756 J. Tian, et al. Remote Sensing of Environment 242 (2020) 111745 were loaded into Google Earth. Fig. 9 shows the S. alterniflora distribution in the Beibu Gulf in the mapping year of 2016-2017. In addition, Fig. S1, a partial zoom in figure in R1 (Fig. 9), was employed to show how difficult it is to detect S. alterniflora. Although the distributions of S. alterniflora are very dense and relatively easy to classify than some other geographic locations, we cannot easily identify them by visual interpretation in the two Landsat images. In the mapping year of 2016-2017, S. alterniflora covered a total area of 1365 ha. The highest spatial coverage of S. alterniflora in the Beibu Gulf was in Guangxi Province, particularly in the Dandou Sea (R1, 539 ha) and Lianzhou Bay (R2, 291 ha) (Fig. 9). The number of patches (i.e. groups of connected S. alterniflora pixels) was 478. The average and maximum area of all the patches were 2.9 ha and 83 ha, respectively. Ever since it was intentionally introduced in China from 1979, the area of S. alterniflora has increased from about 1 ha to 1365 ha. From 1979 to 2017, the invasion speed of S. alterniflora was very high, with severe impacts on ecosystems and biodiversity. Fig. 10 shows the area and Compound Annual Growth Rate (CAGR) (Taubenböck et al., 2014) of S. alterniflora from 1979 to 2017 in the Beibu Gulf. A very small area growth (i.e. from 1 ha to 2 ha) occurred within such a long period during 1979-1995, while the area of S. alterniflora began to increase rapidly from 1996 to 2001 (i.e. from 2 ha to 84 ha) with high CAGR values during this period. The rapid increase indicated the outbreak time of S. alterniflora, particularly in 1996-1997, with CAGR as high as 112%. From 2002 to 2017, the area increased at a relatively low and stable CAGR (< 30%), with an area coverage from 137 ha to 1365 ha. However, it should be noted that the area of S. alterniflora kept growing at a very fast speed with about 82 ha per year.
The Dandou Sea (R1 in Fig. 9) is one of the main distribution regions of S. alterniflora in the Beibu Gulf. Fig. 11 shows the R1 extents during 1995 to 2017. The area and CAGR of S. alterniflora had a similar trend to that of the whole Beibu Gulf. As shown in Fig. 11, the diffusion modes of S. alterniflora included both point-source and boundarysource, in agreement with the reproduction through both seeds and roots. In 2016-2017, S. alterniflora has almost covered the whole mudflats where mangrove forests were not present. The results show that the potential habitats of mangrove forests were seriously invaded by S. alterniflora, which agrees with the field observations (Zhang et al., 2012).

Discussion
In this paper, we developed a new composite method to extract effective phenological feature at the pixel level (i.e. Ppf-CM) to monitor invasive S. alterniflora dynamics from 1995 to 2017 in a cloudy coastal zone, using 525 Landsat images in GEE. This study presents that: (1) Ppf-CM not only accommodates the spatial phenology variability of S. alterniflora as well as enhances the spectral separability, but also eases the problems caused by the scarcity of cloud-free Landsat scenes; (2) incorporation of the new phenological feature yielded from Ppf-CM as input data can improve the performance of the deep learning classifier; (3) This is the first time that invasion outbreak of the S. alterniflora was discovered at the regional scale.
Besides the aforementioned successes, we expect that Ppf-CM in GEE will advance the phenology-based composite image analysis to detect various vegetation species, despite the severe cloud contamination in satellite images. In addition, as revealed by our attempt to incorporate phenological information in deep learning, we have witnessed great potential for developing phenological knowledge based deep learning methods. As a whole, the derived long-term and regionalscale S. alterniflora maps will serve as the intrinsic baseline data for the assessment, management, and conservation of ecosystem services in coastal zones.

Phenological features (Ppf-CM)
The results revealed that Ppf-CM achieved the best performance in the spectral separability statistics (Figs. 7-8), which can be attributed to two factors: First, Ppf-CM minimizes the intra-class variability in spectral reflectance by mitigating the spatial variability of phenology in S. alterniflora. Second, Ppf-CM maximizes the inter-class distances between S. alterniflora and background by identifying and combining two distinctive phenological features (i.e. GCI and SCI images). Particularly, this study site extremely lacks cloudless Landsat images, which has hindered S. alterniflora detection research in the past. The proposed Ppf-CM performed excellently in cloud-minimizing by capturing the phenological features at the pixel level.

Ppf-CM minimizes the intra-class variability
Ppf-CM minimizes the intra-class variability in spectral reflectance by mitigating the spatial variability of phenology. Previous studies (e.g.  Fisher et al., 2006;Friedman et al., 2005;Ji and Wang, 2016) have shown that the spatial variability of phenology is nontrivial in invasive species. For example, prominent phenological differences (i.e. green and senescence periods) can be observed for pixels belonging to the same invasive species across a specific-date Landsat image (Fig. 4). This phenomenon results in large variability in the spectral reflectance of the studied invasive species, i.e. intra-class variability. To address this issue, as the latest development, Ji and Wang (2016) developed a pixelbased method to estimate the timing of coloration across different geographic locations in Landsat imagery with the help of the MODIS End of Season-Time (EoST) product. However, this method was limited by the accuracy of the EoST product and the discrepancies caused by the different spatial resolution of MODIS and Landsat. In this study, in lieu of inclusion of additional data (e.g. MODIS EoST product), our proposed method (Ppf-CM) aims to borrow help from time series of Landsat images. The essence of Ppf-CM is to generate a respective green (maximum EVI) and senescence (maximum PSRI) composite image, in which the phenological variability is largely ameliorated. As a result, the intra-class variability for the invasive species is minimized. Despite the merits of our developed Ppf-CM method, it has to be noted that two limitations still exist. First, PSRI has difficulty in capturing the senescence timing when cirrus cloud is present. Therefore, an enhanced PSRI or a more robust cloud mask method should be developed to address this issue. Second, Ppf-CM cannot capture the green and the senescence timing finer than 16 days due to the moderate temporal resolution of Landsat image. To solve these problems, further work based on temporally refined Landsat image is worthwhile to investigate (Hilker et al., 2009;Zhu et al., 2010). Alternatively, applying Ppf-CM on Sentinel 2 imagery at the 5 days' revisit frequency can be tested as well.

Ppf-CM maximizes the inter-class separability
Ppf-CM maximizes the inter-class distance by identifying and combining the two distinctive phenological features in S. alterniflora detection. Prior studies (e.g. Evangelista et al., 2009) have documented the effectiveness of incorporating time-series images acquired at multiple phenological stages in improving spectral separability of invasive species. However, these studies were scene-based methods which shared two common problems: (a) the spatial variability of phenology cannot be mitigated within an image scene; (b) selecting entire scenes according to cloud coverage reduces the number of available images sharply, which is a particularly prominent phenomenon in the cloudy coastal zone. Problem (a) has been discussed in Section 5.1.1. Problem (b) is a very challenging issue in this study site. As mentioned in Section 2.2.1, taking Frame 2 (see Fig. 1) in 2016-2017 as an example, there was only one cloud-free and low-tide Landsat image (acquired on September 16, 2016) available within the coastal boundary zone. This situation became even worse in some other years, especially when taking both Frame 1 and Frame 2 into account.
In this study, Ppf-CM identified and combined two distinctive phenological features (i.e. green composite image and senescence   Tian, et al. Remote Sensing of Environment 242 (2020) 111745 composite image) to better monitor S. alterniflora, which guided the selection of the spectral data to be composited at the pixel level. In the senescence period, mangrove forest and the other evergreen species are much greener than S. alterniflora, making them easily distinguishable (Ai et al., 2017). However, some mixed pixels of S. alterniflora are extremely difficult to distinguish from the mudflats (Chust et al., 2008). In the green period, although S. alterniflora is hard to be distinguished from the evergreen vegetation, it becomes easy to distinguish S. alterniflora from the mudflats (Ouyang et al., 2013). Therefore, the combination of green and senescence periods can improve the spectral separability in S. alterniflora detection. Compared to scene-based methods, Ppf-CM will optimize the use of each available pixel in Landsat imagery, hence enlarging the amount of remote sensing data for time-series analyses, which successfully addressed the issues of lacking entire scene cloud-free remote sensing images. In addition, only two phenological stages are employed in Ppf-CM rather than stacking multiple phenological stages, reducing demands for the amount of data. Moreover, the two specific stages theoretically accommodate most of the effective phenological information and abandon the redundant phenological information for S. alterniflora detection. As a result, the inter-class separability between invasive species and others are maximized, notwithstanding the scarcity of entire cloud-free image scenes. Although our proposed Ppf-CM in GEE is designed for S. alterniflora detection, it should be beneficial to all phenology-based composite image analysis, especially for long-term and large-scale mapping, with an excellent performance in the aspect of cloud-minimizing. For the next step, Ppf-CM can be expanded to detect various vegetation species, by compositing spectral data from their corresponding phenological features. For example, when applying Ppf-CM to map paddy rice, we expect to find at least three dominant phenological stages: water (sowing), green (emergence and full growth) and bare soil (harvest). Capturing these three stages with spectral data should further enhance the spectral separability between paddy rice and others.

Deep learning with Ppf-CM
This is the first study to our knowledge to investigate the phenological features in improving the performance of deep learning. As revealed in our results, we have seen great potential in developing phenological knowledge based deep learning methods. Recent literature has documented the applications and benefits of deep learning, a new branch of machine learning, in various fields (LeCun et al., 2015). However, it is still a controversial topic whether deep learning is more suitable than traditional machine learning methods for Landsat image classification (Heydari and Mountrakis, 2018;Li et al., 2016). Nevertheless, this study aims to evaluate the hypotheses that applying phenological knowledge as input data of deep learning can improve the performance in the detection of invasive S. alterniflora. A comparison experiment between SAE and SVM with six different types of Landsat images was conducted (see Table 2). The best performance was achieved with Ppf-CM based SAE method, suggesting SAE benefited more than SVM from combining data in the two phenological stages using Ppf-CM (Table 4). This was mainly because, compared to the traditional machine learning method (SVM), deep learning method (SAE) can easily capture and exploit more complex information with its deep network structure and strong self-learning ability (Penatti et al., 2015;Xie et al., 2015).
There exist some future works that are worth noting. First, only the SAE model was employed in this study. It treats each pixel in Landsat imagery as an individual input unit and neglects the spatial information Table 5 Overall accuracy (%) of the six types of Landsat images in nine different mapping years with SAE and SVM.

Images types
Two classifiers 2013-2014 2010-2011 2007-2008 2005-2006 2002-2003 2000-2001 1998-1999 1996-1997 1994-1995 (Finn et al., 2016). Therefore, some other deep learning models (e.g. Fully Convolutional Networks) that accommodates the relationship between one pixel and its neighbors can be tested (Hu et al., 2015). Second, we did not implement SAE on GEE in this study due to the excessive programming efforts to migrate the codes. Fortunately, deep learning methods will become easier to implement on GEE in the near future because Google is trying to import AI models to GEE. Third, in addition to phenological information, some other types of prior   Fig. 9) from 1995 to 2017. The red pixel represents S. alterniflora, and the base map is a Pleiades image (three-channel composite: red, green, blue) at a spatial resolution of 0.5 m acquired in November 2017.
knowledge (e.g. spatial structure, spatial contexture, and texture) should be considered as well to improve the performance of deep learning.

Invasion outbreak of S. alterniflora
This is the first time that the invasion outbreak of the S. alterniflora was discovered at the regional scale, which benefitted from the following three factors: (1) The proposed Ppf-CM eases the problems caused by the scarcity of entire scene cloud-free Landsat scenes; (2) The deep learning classifier optimizes the use of the phenological information yielded from Ppf-CM; (3) GEE provides an ideal platform for handling massive amount of data required by the pixel-based method (i.e. Ppf-CM) with its powerful data storage and high computational capacity.
Prior studies (Beland et al., 2016;Ramsey et al., 2014;Rosso et al., 2006;Wang et al., 2015;Liu et al., 2017) have successfully employed remote sensing techniques to monitor the S. alterniflora dynamics. However, these studies were limited by the spatial or temporal coverage due to the available remote sensing data imagery. In this study, we derived S. alterniflora maps during 1995-2017 in the Beibu Gulf, China. This long-term regional-scale dataset can serve as the intrinsic baseline data for the assessment, management, and conservation of ecosystem services in coastal zones.
The invasion outbreak time of S. alterniflora in the Beibu Gulf, China was found to be during 1996-2001 ( Fig. 10) from 1979 to 2017. This finding is close to but earlier than the conclusion drawn by local forest agencies through long-term observations at a local site in the Dandou Sea. The difference can be attributed to two key reasons: First, the scope of human eyes is limited by the terrain and other environmental factors (Zeleke et al., 2013). Second, although S. alterniflora expanded rapidly during 1996-2001, it can be easily omitted by human observers because it covered only a small area at that time. Compared to the field surveys, remote sensing demonstrates its capability to remotely monitor a very large area in a rapid and responsive manner (O'Connell et al., 2017). From 2002 to 2017, although the CAGR of S. alterniflora was relatively low, the total area kept growing at a very fast speed (about 82 ha per year). As a whole, the areas of S. alterniflora increased from 1 ha to 1365 ha during 1979-2017, which posed significant threats to the biodiversity and environmental health in the coastal zone of China. Therefore, it becomes more critical for the conservation agencies to develop cost-effective control and management strategies, and to restore the coastal areas invaded by S. alterniflora.
Adaptation of our proposed methods from regional scales to national ones still faces a significant uncertainty. The main associated species of S. alterniflora is the mangrove forest in this study site. However, some other associated species will appear in different geographical zones in China, such as Phragmites australis (Ai et al., 2017). To distinguish S. alterniflora from those species may require different phenological knowledge.

GEE provides great potentials for pixel-based analysis
This study developed a new pixel-based method (Ppf-CM) in GEE using 525 full Landsat scenes (19.96 billion pixels) to monitor S. alterniflora dynamics. We found that, as mentioned in Section 5.1, Ppf-CM not only enhances the spectral separability between S. alterniflora and others, but also accommodates the problems caused by the scarcity of entire cloud-free Landsat scenes. These findings are in general agreement with prior GEE-supported pixel-based studies (e.g. Azzari and Lobell, 2017) and confirm that pixel-based methods outperform scenebased methods to monitor S. alterniflora. Although GEE-supported methods have been increasingly leveraged in various fields in the last three years, such as agriculture (Phalke and Özdoğan, 2018) There still exist some limitations of GEE. For example, we cannot save a geometry file larger than 512 KB, and it would be desired to make the deep learning methods much easier to implement on GEE platform. In the near future, it is foreseeable that a growing number of large-scale mapping studies with GEE will emerge, which benefits a lot for us to understand, manage, and protect the earth. However, we also should pay attention to some issues: (1) We call for worldwide researchers to release their derived maps in GEE so that others can test, use, or improve them; (2) GEE should establish sound accuracy assessment, feedback, and modification mechanism to improve the accuracy of the released products; (3) GEE should encourage all the users around the world to upload their high quality ground observation data to advance large-scale mapping studies.

Conclusions
This paper proposed a new pixel-based cloud-minimizing phenological feature (Ppf-CM) for deep learning to understand Spartina alterniflora invasion in GEE. Ppf-CM mitigated the spatial variability in plant phenology, improved the spectral separability between S. alterniflora and background, and accommodated the impact of scarce entire scene cloud-free Landsat imagery. It will advance all the phenology-based composite image analysis at the pixel level to detect various vegetation species, notwithstanding the large impact of clouds on the quality of full Landsat scenes. Moreover, we found that using phenological features yielded from Ppf-CM improved the performance of a deep learning method (i.e. SAE) in detecting S. alterniflora. Hence, we have seen great potential in developing prior knowledge based deep learning methods. Finally, long-term and regional-scale maps of S. alterniflora were derived using the most powerful detection strategy (i.e. Ppf-CM based SAE) during 1995-2017 in the Beibu Gulf, which discovered the outbreak of S. alterniflora during 1996-2001 for the first time at the regional scale. The long-term and regional-scale S. alterniflora maps can serve as the intrinsic baseline data for the assessment, management, and conservation of ecosystem services in this study area, including biological invasion, biodiversity richness, climate change, and carbon storage.