Capability of Phenology-Based Sentinel-2 Composites for Rubber Plantation Mapping in a Large Area with Complex Vegetation Landscapes

: Mapping rubber plantations in a large area is still challenging in high-cloud-cover and complex-vegetation landscapes. Existing studies were often conﬁned to the discrimination of rubber trees from natural forests and rarely concerned other tropical tree species. The Sentinel-2 constellation, with improved spatial, spectral, and temporal resolution, offers new opportunities to improve previous efforts. In this paper, four Hainan Sentinel-2 composites were generated based on the detailed phenological stages delineation of rubber trees. The random forest classiﬁer with different phenological stage combinations was utilized to discuss the capability of Sentinel-2 composites to map rubber plantations. The optimal resultant rubber plantation map had a producer’s accuracy, user’s accuracy, and F1 score of 81%, 84.4%, and 0.83, respectively. According to the rubber plantation map in 2020, there was a total of 5473 km 2 rubber plantations in Hainan, which was 2.93% higher than the statistical data from the Hainan Statistical Yearbook . According to the Hainan Statistical Yearbook , the area-weighted accuracy at the county level was 82.47%. The mean decrease in accuracy (MDA) was used to assess the feature importance of the four phenological stages. Results showed that the recovery growth stage played the most important role, and the resting stage was the least important. Moreover, in terms of the combinations of phenological stages, any dataset group with two phenological stages was sufﬁcient for rubber tree discrimination. These ﬁndings were instrumental in facilitating the rubber plantation mapping annually. This study has demonstrated the potential of Sentinel-2 data, with the phenology-based image-compositing technique, for mapping rubber plantations in large areas with complex vegetation landscapes.


Introduction
The rubber tree (Hevea brasiliensis) is the primary source of natural rubber, an essential industrial raw material. With the increasing demand for natural rubber, rubber plantations have expanded to almost every tropical forest region in the world. The expansion of rubber plantations has contributed to local economic development. However, the transition from natural forests to rubber plantations has significant ecological impacts on the water balance, carbon cycle, biodiversity, and ecosystem function [1][2][3][4]. Knowledge of the spatial extent and dynamic of rubber plantations is significant to ensure the sustainability of the natural rubber industry, including plantation management, rubber futures, national economic policy, and ecological conservation.
Remote sensing technology is an important tool in mapping rubber plantations at local and regional scales [5]. Since 2012, rubber's unique phenological characteristics have been widely exploited in the delineation of rubber plantations [5]. When rubber trees were 2 [26], dense Sentinel-2 time series [28], and phenology/seasonal Sentinel-2 composites [29] can improve the classification of tree species. For rubber plantation mapping, Xiao et al. showed the significant potential of Sentinel-2 imageries with red-edge spectral indices [9]. However, the nonrubber tree species distinguished from rubber plantations were still only natural forests.
This study's main objective was to investigate the utility of phenology-based Sentinel-2 composites for rubber plantation mapping in a large cloudy area with complex vegetation landscapes. Specifically, the following questions were addressed: Which phenological stage of rubber trees is important for rubber plantation mapping? What combinations of rubber tree phenological stages contribute to increasing mapping results? Which bands of Sentinel-2 are important for rubber plantation mapping? How does the band importance of Sentinel-2 vary over the phenological stages of rubber trees?

Study Area
Hainan Province is located in the southern part of China ( Figure 1). It is characterized by a marine tropical monsoon climate. The annual average precipitation is more than 1600 mm, and the mean annual temperature is 23~25 • C. The annual minimum temperature is generally above 5 • C. Before the 20th century, Hainan Province was dominated by tropical natural forests, and rubber trees were introduced from British Malaysia. Since the 1950s, rubber trees have been planted in large quantities through the deforestation of natural forests. After half a century of development, Hainan Province has become the second-largest rubber-planting area in China after Yunnan Province, and natural rubber production accounts for about 40.5% of China's total production.
The deciduous habit of rubber trees is one of its most critical phenological characteristics, allowing annual leaf renewal [34]. Generally, the defoliation periods of rubber trees last for about one month between February and March [21]. In this study, according to rubber trees' phenology and growth rhythms [35], the growth cycle of rubber trees was divided into four phenological stages: 1 resting phase, 2 recovery growth phase, 3 vigorous growth phase, and 4 slowdown growth phase. The photos in Figure 2 show the status of rubber trees at the four phenological stages, respectively. Detailed descriptions of the four phenological stages and their timing will be presented in Section 3.1. In addition to rubber trees, there are a variety of tropical economic tree species in Hainan Province, such as Eucalyptus, Coconut, Areca nut, and Litchi. All the economic forests, including rubber trees, are widely distributed in the study area, with complex and fragmented landscapes.

Sentinel-2 Data
The Sentinel-2 mission is based on a constellation of two satellites, Sentinel-2A (S2A) and Sentinel-2B (S2B). Combining S2A and S2B can provide a five-day revisit time at the equator and better temporal resolution at higher latitudes. In this study, we collected all the Sentinel-2 Level 2A products of the study area acquired from 1 January 2019 to 31 December 2021. The Sentinel-2 Level 2A products were downloaded from the Copernicus Open Access Hub of the European Space Agency (ESA). The Level-2A product is composed of 100 × 100 km 2 tiles, and Hainan Province is covered by eight tiles, as shown in Figure 3.
Hainan Province has cloudy weather all year round. The average percentages of Sentinel-2 images with cloud cover (CC) less than 10% for the four phenological stages were 29.57%, 14.36%, 14.21%, and 16.65%, respectively. Therefore, for each phenological stage, it is tough to obtain less-cloud or cloud-free Sentinel-2 images for all eight tiles. In this study, to obtain four cloud-free composites for the entire study area in Section 2.4.1, Sentinel-2 images with an estimated CC ≤ 75% were selected for further preprocessing.
The preprocessing of Sentinel-2 imageries included cloud and cloud shadow masking. The cloud masks were converted from the layers of Sentinel-2 cloud probability, and the cloud shadow masks were generated based on the intersection of cloud projection and nearinfrared (NIR) band darkness. Finally, the 20 m bands were resampled to 10 m resolution using nearest neighbor resampling, and the 10 Sentinel-2 10/20 m bands (B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12) were merged into a single raster dataset. Hainan land cover from "ChinaCover2020" and Sentinel-2 tiles for the study area.

MODIS Data
Due to the frequent cloudy weather in Hainan Province, it is tough to construct consistent year-long Sentinel-2 time series with reliable data quality. Therefore, in this study, we used coarse spatial resolution MODIS data to delineate the phenological stages of rubber trees. The timing of the lowest NDVI according to the MODIS-NDVI time series curve has been used to determine the periods of rubber tree defoliation and foliation [8,13,36]. This study used two MODIS-NDVI products with 250 m resolution, MOD13Q1-NDVI and MYD13Q1-NDVI. The two products showed high compatibility [37,38], and it has been shown that they can be used in combination [25,39]. These products were obtained through the online MODIS data repository (https://modis.gsfc.nasa.gov/data/dataprod/mod13. php, accessed on 3 March 2022).
The two NDVI products for the study area from January 2019 to December 2021 were rearranged by the property of start time. Subsequently, the combined NDVI time series were smoothed using the locally linear regression smoothing technique [40].

Land Cover Data
The main objective of this study was to distinguish rubber plantations from other tropical tree species. Therefore, we used the Hainan land cover dataset from "ChinaCover2020" to mask non-forest land ( Figure 3). The "ChinaCover" product was generated based on an object-based approach [41,42], and the latest version is "ChinaCover2020" with a 10 m resolution. The classification accuracy of the Hainan land cover dataset was evaluated using 1200 independent ground survey samples (712 samples for forest and 488 for non-forest) that were collected using a random sampling approach [43]. The producer's accuracy and user's accuracy for the forest category were higher than 95%. Therefore, the forest map can serve as a reliable base dataset for rubber plantations' delineation.

Rubber Plantations' Area Data from Statistical Yearbook
The Statistical Bureau of Hainan Province (SBHP) published annual reports on the rubber plantations' statistical data [44]. The data on rubber plantations were from the Statistical Survey System on Agriculture, Forestry, Animal Husbandry and Fishery of Hainan Province. Statistical departments at all levels collected data by means of sample surveys, key surveys, or comprehensive surveys and reported them layer by layer according to the actual local conditions. The area of rubber plantations was recorded for each of the 18 counties in Hainan Province. According to the statistical yearbook, in 2020, the rubber plantation area was 5316.72 km 2 . In this study, this dataset was used to assess the accuracy of the resultant rubber plantation map at the county level.

Field Survey Data
Two field works were conducted in August 2017 and October 2020. The primary purpose of the field surveys was to collect samples for updating the "ChinaCover" product. As shown in Figure 1, the field surveys cover all counties in Hainan Province. The reference samples were selected following two rules: the samples were located at the center of homogeneous forest patches, and all the trees in a plot were mature forests.
The typical tropical tree species in Hainan Province are Rubber tree, Eucalyptus, Litchi, Coconut Palm, Areca Palm, and Casuarina. In addition, the natural forest is mainly distributed in the central mountainous areas of Hainan, which the field surveys did not cover, so we collected the natural forest samples from the high-resolution images in Google Earth, as shown in Figure 1.
In total, 2604 reference samples were collected. All samples were divided into two sets, one designed for training and the other for assessing classification accuracy. For convenience, coded names were given for each tree species, as shown in Table 1.

Experimental Design
A comprehensive overview of this study is shown in Figure 4. First, the temporal behaviors based on MODIS-NDVI were studied to delineate the phenological features of rubber trees and determine the timing of the four phenological stages. Next, phenologybased image compositing was used to generate four Sentinel-2 composites, and 15 dataset groups were generated with different phenological stage combinations. Then, using the random forest classifier and Jeffries Matusita distance, the classification accuracy and class separability were analyzed for the 15 dataset groups. After that, the spatial distribution of the resultant rubber plantation map was further validated based on the Statistical Yearbook of Hainan Province in 2020. At last, the mean decrease in accuracy was used to assess the feature importance of Sentinel-2 bands and rubber trees' phenological stages for rubber plantation mapping.

Phenology-Based Image Compositing
Image compositing is an approach to reduce a series of images into a single image [45]. Recently, it has been introduced to generate cloud-free Sentinel-2 composite datasets at large scales [46,47]. Additionally, the phenological and seasonal Sentinel-2 composites have been used for tree species mapping [29].
In this study, we adopted the technique of median compositing to generate compositing images for each Sentinel-2 tile and phenological stage. The temporal intervals in the compositing process were set following the timing of the corresponding phenological stages, which will be presented in Section 3.1. When there was no good-quality observation during the phenological stage, we filled the data gap with the nearest good-quality observation in the Sentinel-2 time series. At last, four Hainan Sentinel-2 composites were generated based on the seamless mosaic tool in ENVI 5.3. The four composites were named as 1 Resting, 2 Recovery, 3 Vigorous, and 4 Slowdown, respectively, corresponding to the four phenological stages.
To further evaluate the importance of the different phenological stages and the potential of multistage combinations for rubber plantation mapping, the four Hainan Sentinel-2 composites were rearranged as 15 dataset groups: four mono-stage images, six doublestage combinations, four three-stage combinations, and one all-stage combination. For convenience, the 15 dataset groups were named by the serial numbers of the corresponding phenological stages. Taking the all-stage combination as an example, it was called 1 2 3 4 for short.

Class Separability Based on Jeffries Matusita Distance
The performance of the 15 dataset groups in separating rubber trees from other tree species was assessed using the Jeffries Matusita (JM) distance in this study. The JM distance is the average distance between two class density functions, and it can be calculated as [48]: where the B ij is Bhattacharyya distance as where µ i and c i are the mean and covariance matrix of the i-th category. The JM distance is a parametric test, ranging between 0 and 2. It provides an easy comparison of class separability, with 0 indicating no separability and 2 complete separation.

Random Forest Classifier and Feature Importance
In this study, the random forest (RF) classifier was used to map rubber plantations. RF is a supervised machine learning algorithm for classification and regression. It is an ensemble of decision trees, which are constructed based on the bootstrap aggregated sampling (bagging). RF takes the majority vote of the decision trees for classification. In this study, the randomForest package in the software R was employed to build the RF model [49]. Unified parameters were set in the RF model for all 15 dataset groups. The number of features for each split (parameter mtry) was set as the square root of the total feature number. The number of decision trees (ntree) in the RF model was set to 500, as this value was commonly used for remote sensing classification [29].
In addition, the mean decrease in accuracy (MDA) was used to assess the importance of the 10 Sentinel-2 bands and four phenological stages for rubber plantation discrimination. MDA, also known as permutation importance, is one of the most efficient feature importance measures for the random forest.

Accuracy Assessment
To assess the performance of Sentinel-2 images on rubber plantations discrimination, we carried out two levels of classification accuracy. (1) Confusion matrix based on the validating samples, with the estimated user's accuracy (UA), producer's (PA) accuracy, and F1 score; (2) statistical data validation at the county scale.

Phenological Stages Delineation of Rubber Trees
The temporal profile of rubber trees from 2019 to 2021 was plotted based on the MODIS-NDVI, as shown in Figure 5. The blue and orange curves represent the original and smoothed NDVI time series. The troughs of the smoothed curve are indicated with red arrows and the corresponding dates.
In Figure 5, the four phenological stages are labeled by the corresponding numbers, and the vertical blue dotted lines separate the timings of the stages. The resting stage lasts for about 50 days [50]. During this stage, the old leaves of rubber trees have turned red and yellow, and concentrate on falling off the rubber trees. According to the long-term collection of fallen leaves on a field, the amount of fallen leaves accounts for more than 70% of the whole year [50]. In this study, the resting stage's start date and end date were set 25 days before and after the dates of troughs. In the recovery growth stage, the rubber trees sprout new leaves, and the first awning leaves complete the process of leaf development. During this stage, the number of leaves accounts for about 70% of the whole year [20]. This stage lasts for about 50 days after the resting stage. The slowdown growth stage lasts for about 50 days before the resting stage. In this stage, affected by cold air masses, the rubber trees gradually enter the defoliation period. The rubber growth tends to stop, and the old leaves begin to turn red. The vigorous growth stage is the longest stage in the rubber trees' growth cycle, lasting for about 210 days, from middle May to late November. In this stage, the rubber trees' leaf area index maintains a high level. The main causes of NDVI fluctuation are atmospheric conditions and precipitation.  Table 2 shows the J m distance between rubber trees and other tree species. In this study, we separated the J m distance into four levels and marked it with different colors: strong (1.9-2.0, blue), good (1.8-1.9, green), weak (1.7-1.8, yellow), and poor (<1.7, pink). Taking the dataset group 1 3 (combination of the resting stage and vigorous growth stage) as an example, the J m distance between rubber trees and Lychee was 1.85, corresponding to the level of good separability.

Separability between Rubber Trees and Other Tree Species
The J m distance in Table 2 varied greatly with tree species and dataset groups (phenological stage combination). Casuarina and Coconut Palm were most easily distinguishable from rubber trees among the six tree species. For all the 15 dataset groups, the J m distances were larger than 1.9, which indicates strong separability. As for natural forest, when the dataset groups contained two or more phenological stages, it was easily distinguishable from rubber trees with strong separability. However, when only one phenological stage was contained, the separability between natural forest and rubber trees was less than 1.8, corresponding to weak or poor separability. As for Lychee, Areca Palm, and Eucalyptus, when three or more phenological stages were contained, the J m distance was at the level of strong separability. When a single phenological stage was contained, the J m distance was less than 1.7, and these three tree species were difficult to distinguish from rubber trees.
A J m distance larger than 1.8 is regarded as satisfactory discrimination between different classes. In general, based on Sentinel-2 data, any combination of two phenological stages was good enough for distinguishing rubber plantations from the other six tree species.

Accuracy Assessment Using Survey Data
The rubber plantation mapping results based on the 15 dataset groups were validated using the 1554 validating samples. The measures UA, PA, F1 score, and true positive (TP) for rubber plantations are shown in Table 3. In addition, the sources of commission and omission errors for rubber plantations are also included in Table 3. Judging by the UA, PA, and F1 score, it was apparent that the dataset group of all-stage combination ( 1 2 3 4 ) was superior to other 14 dataset groups. The UA, PA, and F1 score for rubber plantations were 81.0%, 81.4%, and 0.83, respectively. The best dataset group with three phenological stages was the combination of the resting, recovery growth, and slowdown growth stage ( 1 2 4 ), with a UA, PA, F1, and true positive of 79%, 83.6%, 0.81, and 158, respectively. The best dataset group with two phenological stages was the combination of recovery growth and slowdown growth stage ( 2 4 ), with a UA, PA, F1, and true positive of 78.5%, 80.5%, 0.79, and 157, respectively. In general, the more phenological stage in the dataset group, the higher the F1 score. However, when the dataset group contained more than three phenological stages, the difference in mapping accuracy was slight.
The commission and omission errors were mainly from Eucalyptus and Areca Palm, accounting for about two-thirds of the total errors. The rest of the classification errors were from Lychee and the natural forest. The classification errors from Casuarina and Coconut Palm were negligible.

Rubber Plantation Map and Statistical Data Validation
The result of the rubber plantation map based on the all-stage combination dataset group ( 1 2 3 4 ) is shown in Figure 6. The distribution of rubber plantations was in good agreement with the resultant maps in previous studies [13,15,16,18]. The codes labeled in Figure 6 correspond to the 18 counties in Hainan Province, and the one-to-one correspondence is shown in Table 4. The densest rubber plantation area was located in the northwest of Hainan Province, including the counties of Danzhou (3), Chengmai (11), Lingao (12), and Baisha (13), where there is the largest natural rubber production base in Hainan Province. The numbers behind the counties correspond to the codes in Figure 6.
The spatial distribution of the resultant map was further validated based on the Statistical Yearbook of Hainan Province in 2020. This study estimated the rubber plantation area of Hainan in 2020 was approximately 5473 km 2 , which is about 2.93% higher than the 2020 statistical data (5316 km 2 ). Additionally, a county-level comparison of rubber plantation area between this study and the statistical yearbook is shown in Table 4. The 18 counties were sorted based on the rubber plantation area estimated in this work. In Table 5, the accuracy at the county level was calculated as follows: where the operator abs() denotes the absolute value, and the area weight is the proportion of each county's total Hainan Province area. The area-weighted accuracy in total is the sum of the area-weighted accuracy for each county. The accuracy of the rubber plantations area at the county level ranged from 48.24% in Lingshui (16) to 95.09% in Chengmai (11). In total, the accuracy was as high as 97.05%, and the area-weighted accuracy was 82.47. Therefore, it could be concluded that the resultant map of rubber plantations had high accuracy in county-level spatial distribution.  Table 5 shows the feature importance for the rubber tree discrimination that employed all four phenological stages (dataset group 1 2 3 4 ). The colored cells denote the MDA values for each feature. The bottom row in Table 5 shows the sum of the feature importance for each phenological stage, and the last column offers the sum of the feature importance for each Sentinel-2 band.

Feature Importance
The recovery growth stage was the most important phase, with a cumulative MDA of 109.88, followed by the vigorous and slowdown growth stages, with an MDA of 98.76 and 91.72, respectively. The resting stage was the least important phase with an MDA of 76.71.
The visible bands were much more important in the resting, recovery growth, and slowdown growth stages than in the vigorous growth stage. The SWIR bands showed much higher feature importance in the vigorous growth and slowdown growth stages than in the resting and recovery growth stages. Two red-edge bands (B6, B7) and the two NIR bands showed higher feature importance in the recovery growth and vigorous growth stages than in the resting and slowdown growth stages. Only the red-edge band B5 always had low importance values during all four phenological stages.
The importance of each Sentinel-2 band was uneven when summarized over the four phenological stages. The three visible bands (B2, B3, and B4) were the most important, followed by the two SWIR bands (B11 and B12) and red-edge band B6. The remaining bands in the red-edge region (B5 and B7) and the NIR region (B8 and B8A) showed low importance values.

Capability of Sentinel-2 for Rubber Plantation Mapping
Xiao et al. [9] have presented the potential for mapping rubber plantations at a regional scale from three advantages of Sentinel-2 imagery, i.e., the three red-edge bands, 10/20 m spatial resolution, and 5 days' revisiting time. The study area of this study [9] is a relatively small area, which facilitates the selection of cloud-free imagery. However, for a large area such as Hainan Province, with complex vegetation landscapes and frequent cloudy weather, the capability of Sentinel-2 for rubber plantation mapping needs to be further explored.
New in this study was the phenology-based image compositing for rubber tree discrimination. Recent studies have reported the advantages of time series images [28,29] or phenological features [51] for tree species classification. Compared to these studies, our approach is more suitable for rubber plantation mapping in a large area.
First, image compositing is necessary due to the high probability of CC in Hainan Province. The convolution filter to generate gap-free time series [28] is not recommended if the study area has high CC values as the Hainan Province. Due to climatological gradients and orbit swath overlaps, large areas exhibit higher variability in data availability [52,53]. In this case, image compositing with the same time intervals for different Sentinel-2 tiles can ensure the effectiveness of the rubber tree identification model for the whole area.
Second, the time intervals of image composition are set according to the timings of the rubber trees' phenological stages. The timings of the phenological stages are determined based on the smoothed annual MODIS-NDVI curve. For tree species classification, it has been highlighted that input features should be chosen according to the considered species and their phenological characteristics in order to include key phenological stages that enhance their separability [29]. However, the phenology of rubber trees is highly sensitive to climate change, particularly rainfall and temperature [25,54,55]. The changes in climatic factors may delay or advance the defoliation period. As shown in Section 3.1, the timings of the phenological stages change over the three years. Additionally, the durations of the slowdown growth, resting, and recovery growth stages are short. Therefore, the threemonthly or seasonal compositing techniques [29] are unsuitable for rubber tree mapping.
Third, phenological features (e.g., the start time of the phenological stage and length of the growing season) have proven helpful in tree species classification [51]. However, these metrics have rarely been utilized in rubber tree mapping. The main reason is still the frequent cloudy weather in the distribution areas of the rubber trees. With poor temporal accuracy, the derived phenological features could have little impact on rubber tree mapping accuracy.

Phenological Stage Importance
Regarding the importance of phenological stage for rubber tree mapping, based on the J m distance in Table 2, the MDA in Table 5, and the F1-score in Table 3, it can be concluded that the recovery growth stage plays the most important role, while the resting phase plays the least important.
The importance of the phenological stage is directly associated with phenological events (or status) of rubber trees and other tree species. In the recovery growth stage, for rubber plantations, the canopy recovers with new leaves' emergence and formation, while for other tree species, many leaves fall off trees due to spring drought. Shi et al. [56] analyzed the seasonal dynamics of litterfall production of tropical natural forests in Hainan Province. Results showed that spring (March and April) was one of the two peak seasons for litterfall production, and the main influencing factor was drought. In general, new leaves have higher chlorophyll contents than old leaves. The differences in phenological traits between leaf on and leaf off, coupled with the spectral differences between fresh leaves and old leaves, are the main reasons for the most important role of the recovery growth stage.
Unexpectedly, the resting stage plays the least important role among the four phenological stages. As almost all the leaves of rubber trees drop off in the resting stage, it was regarded as the optimal time window or the key phenological stage for mapping rubber trees [8,9,23]. The reason for the least important role of the resting stage could be twofold: understory vegetation in rubber plantations or unsynchronized phenology of rubber trees. Chen et al. [18] also noted they were the two important factors that affect the rubber tree mapping accuracy. In the resting stage, the canopy's opening accelerates the growth of understory vegetation. The mixture of understory vegetation and withered leaves makes the spectrum signature of rubber plantations uncertain. As for the unsynchronized phenology, Hu et al. reported that the spatial characteristics of rubber tree phenology were consistent with the spatial distribution of terrain and elevation in Hainan Province [57]. Chen et al. [18] noted that the phenology of rubber trees was closely related to the factors of stand age, site-specific environment, etc.
The least important role of the resting stage may change if the texture features based on high spatial resolution images are considered for rubber tree mapping. The forest texture is closely related to the canopy structure, which is defined by foliage properties and branch arrangement [58]. In the resting phase, the defoliation of rubber trees exposes crown branches and increases canopy roughness. While for other tree species, leaves still dominate the crown. These differences enhance the separability between rubber trees and other tree species.

Band Importance
The feature importance of Sentinel-2 bands for rubber tree mapping partly agrees with previous studies using Sentinel-2 for tree species classification [26,29,31,33,59]. It should be noted that in previous studies, band importance was calculated for all tree species, but in this study, feature importance was quantified specifically for rubber trees. In this study, the band importance in different phenological stages can be associated with phenological events of rubber trees.
In the vigorous growth stage, all the tree species have lush foliage. The red-edge bands B6 and B7, the two NIR bands, and the two SWIR bands play important roles with high MDA values, which agree well with the results in [26,31,33,59]. Immitzer et al. [33] and Nelson et al. [59] reported that the red-edge bands B6 and B7, the narrow NIR band, and the SWIR 2 band were among the most important bands. Persson et al. [26] also reported that the two SWIR bands were among the highest-ranked bands. The red band (B4) and red-edge band (B5) showed the lowest importance, which partly agrees with the previous results that the red band of a summer acquisition ranked considerably lower than most other summer bands [29,31,32].
In the slowdown growth stage, the leaves of rubber trees gradually turn red and yellow, but for other tree species, the leaves are still green. The red-edge bands B6, B7, and the two NIR bands turn less important, while the SWIR band B11 still shows the second-highest importance. This result is also confirmed by Koller et al., that the SWIR region is especially sensitive to leaf water content, implying that species-specific dynamics of leaf water content during phases of phenological transition help to discriminate tree species [29]. The three visible bands (B2, B3, and B4) play important roles in the slowdown growth stage, since the difference between the yellow (or red) leaves and the green leaves is visible to the naked eye.
In the resting stage, as discussed in Section 4.2, the spectral signatures of rubber plantations are actually from the mixture of understory vegetation and withered leaves. Due to the uncertainty of the mixing ratio, the bands from red-edge to SWIR further decline in importance; however, the three visible bands show higher importance.
In the recovery growth stage, new green leaves are gradually restored in the canopy. Due to the clear difference between fresh leaves and old leaves in the visible region, the three visible bands remain the top three important features. This result is also confirmed by Persson et al., that the blue and green bands from May were in the top seven highest-ranked bands [26]. In addition, Immitzer et al. showed the importance of the blue band [33], and Nelson reported that the red band was more important, probably since an image from early May was included [59].

Research Limitations and Prospects
(1) Previous studies [28,29,60] have confirmed the usefulness of vegetation indices (such as NDVI and EVI) for tree species classification. In the scope of rubber tree mapping, the spectral-indices-based decision trees have been used to discriminate rubber trees from natural forests. It can be expected that the spectral indices are helpful for rubber tree mapping. In addition, several red-edge-related spectral indices can be obtained based on Sentinel-2 imagery [61], and the importance of red-edge bands can be explored deeply. (2) In this study, based on the RF algorithm, the MDA was used to evaluate the feature importance of the phenological stages and the Sentinel-2 bands for rubber plantation discrimination. However, the selection of features with high importance does not warrant that this is the best set of features for a given problem [62]. Features with high correlation may reduce the reliability of the RF-based MDA importance, and have a negative effect on the feature selection. Different solutions have been proposed to overcome some of the known flaws of MDA [63], and several methods have been proposed to select the optimal subset of features [64]. These techniques have been used in forest parameter monitoring, such as tree species diversity [65], growing stock volume [66], and forest stand parameters [66]. However, the impacts of dependent input features on the RF-based MDA importance for tree species classification have not been discussed. In this study, there is no doubt that there is a high correlation between the Sentinel-2 imageries with adjacent phenological stages. In addition, the adjacent bands of Sentinel-2 data are possibly correlated to each other because of the continuity of the bands. Therefore, the feature importance of the phenological stages and the Sentinel-2 bands for rubber tree discrimination needs to be explored deeply in the future. (3) In this study, the four Hainan Sentinel-2 composites were generated based on imageries of three years due to the high cloud cover. This technique is time-consuming and easily affected by the frequent deforestation and reforestation in Hainan Province. Our results showed that any dataset group with two phenological stages was sufficient for rubber tree mapping. Actually, for each Sentinel-2 tile, it is achievable to composite two cloud-free images annually with different phenological stages. Therefore, based on the four Hainan Sentinel-2 composites, future studies should tap the potential of dataset groups with double phenological stages. Standardizing the rule sets for identifying rubber plantations, especially in decision trees, could facilitate the rubber plantation mapping annually. (4) As we know, the first Sentinel-2 satellite was launched in 2015. If we want to map the rubber plantations in Hainan Province before the Sentinel-2 data was available, we still have to resort to Landsat series data. Given the high cloud cover in Hainan Province and the 16-day repeat cycle of Landsat satellites, it is tough to generate four Hainan compositing Landsat images corresponding to the four phenological stages. Recent studies have shown the harmonization of Landsat and Sentinel-2 data [67], and the combination of Sentinel-2 and Landsat for land surface phenology characterizing [68,69]. Therefore, simulated Hainan Landsat images can be generated by adjusting Sentinel-2 radiometry to replicate the spectral bandpasses of Landsat 5/TM or 8/OLI for the bands common to both sensors. In the future, the rubber tree mapping model based on the simulated Landsat images should be explored and will be used to monitor the dynamics of rubber plantations during 1990-2020.

Conclusions
Knowledge of the spatial extent of rubber plantations at a regional scale is significant to ensure the sustainability of the natural rubber industry and ecological conservation. This work used three years of Sentinel-2 data to map rubber plantations in Hainan Province, China. Unlike previous studies, six tropical tree species were involved in being distinguished from rubber trees, and the growth cycle of rubber trees was divided into four stages: resting, recovery growth, vigorous growth, and slowdown growth stage. A detailed MODIS-NDVI curve was used to determine the timings of the four phenological stages, and accordingly, four Sentinel-2 composites were generated. The random forest classifier and the four phenological stage composites were used to identify rubber trees. The capability of phenology-based Sentinel-2 composites for mapping rubber plantations was shown in two levels of classification accuracy: the F1 score based on survey samples was 0.83, and the area-weighted accuracy at the county level based on the Statistical Yearbook of Hainan Province was 82.47%.
Although our study demonstrated the capability of the phenology-based Sentinel-2 image compositing approach for mapping rubber plantations, it is time-consuming and easily affected by frequent deforestation and reforestation. In the future, the utilization of vegetation indices and red-edge-related spectral indices will be specifically addressed. Rulebased decision trees should be built for any two phenological stages to facilitate the rubber plantation mapping annually. Moreover, based on the four compositing images in this work, the potential of Landsat data in mapping rubber plantations deserves further exploration.