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

Abstract

Mapping rubber plantations in a large area is still challenging in high-cloud-cover and complex-vegetation landscapes. Existing studies were often confined 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 classifier 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² 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 sufficient for rubber tree discrimination. These findings 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.

1. 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 introduced from the Amazon to Southeast Asia in the nineteenth century, rubber trees became deciduous in the winter or dry season with seasonal drought and cold stresses [6]. For rubber trees, most leaves fall off during the defoliation process before new leaves emerge in the refoliation process [7]. In contrast, most other tropical forests are evergreen forests, having strong resistance against influences from the physical environment, such as temperature and precipitation. Multiple phenology-based features using vegetation indices [8,9,10,11,12,13,14,15,16] have been exploited to map rubber plantations and have shown some advantages [5]. (1) It is easier to distinguish rubber trees from cropland and natural forest [17]; (2) it helps to identify and outline rubber plantations within fragmented tropical mountainous lands [12].

Hainan Province is China’s second-largest natural rubber planting area after Yunnan Province. Several studies have been published to map rubber plantations in Hainan Province [8,13,15,16,18,19,20]. Among these studies, the decision tree classifier with multiple phenology vegetation indices was the most common approach. Dong et al. [16] and Chen et al. [18] utilized the L-band PALSAR product to map forests and then extracted rubber plantations from the forest map according to their phenological features based on MODIS and Landsat data, respectively. Han et al. [13] developed a simple decision rule based on band value changes, vegetation index changes, and phenological phase changes (defoliation and foliation) to distinguish rubber plantations from other land-cover types. Cui et al. [8] explored the suitability of Planet Scope imagery for mapping rubber plantations, and the spectral features, index features, and textural features were comprehensively utilized. These studies obtained high classification accuracies of 85%, 95.2%, 92.17%, and 95.18%, respectively, but there are still three major issues that influence the rubber plantation mapping results. (1) The mapping of rubber plantations in previous studies was mainly based on MODIS [16] or Landsat data [13,15,18]. The 250 m spatial resolution of MODIS products is relatively coarse, making it difficult to delineate the border of rubber plantations in Hainan Province with fragmented vegetation landscapes. As for Landsat data, it is difficult to obtain valid cloud-free imageries due to the cloudy weather condition in Hainan Province, the short period of rubber defoliation, and the 16-day repeat cycle of Landsat satellites. (2) Previous studies [8,9,13,15,16,21] focused on the discrimination of deciduous rubber plantations from evergreen natural forests, but little consideration was given to other tropical tree species. Hainan Province is abundant with tropical tree species, such as Eucalyptus, Areca Palm, and Coconut Palm, and they are widely distributed. Some tropical economic forests were mentioned briefly [16,22], but the difference between rubber tree plantations and the economic forests has not yet been explored. (3) The defoliation and foliation periods were regarded as the optimal timing for distinguishing rubber trees from evergreen natural forests [8,11,23]. These two key periods were also called the rubber defoliation and foliation (RDF) period [18], leaf-off and leaf-on phases [24], senescence and regreening periods [22], and the wintering period [5,25]. However, the timing of the two periods was defined differently, for example, late February to April [15], January to March [10,15], January to April [8], December to late February [23], and November to March [11]. This confusion may affect the selection of imageries.

As stated above, Hainan Province is a large cloudy area with complex vegetation landscapes. Cloud-free multiphenology imageries are necessary to generate an accurate map of rubber plantations, and these imageries must be able to distinguish rubber trees from other tropical tree species. Compared to MODIS and Landsat data, Sentinel-2 can acquire optical imagery at high spatial (maximum 10 m) and temporal resolution (2–5 days revisit time at the equator). The combination of freely available satellite imagery, high resolution, a swath width of 290 km, and frequent revisit times has greatly enhanced the feasibility of mapping rubber plantations in large cloudy areas. In addition, Sentinel-2 incorporates three new spectral bands in the red-edge region, centered at 705, 740, and 783 nm, respectively. Recently, many studies have investigated the utility of Sentinel-2 imageries for precise regional-tree-species mapping [26,27,28,29,30]. The novel red-edge bands are significant for tree species classification [26,31,32,33]. Particularly, multitemporal Sentinel-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?

2. Materials and Methods

2.1. 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: ① resting phase, ② recovery growth phase, ③ vigorous growth phase, and ④ 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.

2.2. Data

2.2.1. 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² 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 near-infrared (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.

2.2.2. 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].

2.2.3. 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.

2.2.4. 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². In this study, this dataset was used to assess the accuracy of the resultant rubber plantation map at the county level.

2.2.5. 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. Number of training and validating samples.

Coded Name	Tree Species	Training Samples	Validating Samples	Total
RT	Rubber tree	150	200	350
LC	Lychee	150	225	375
CA	Casuarina	150	191	341
CP	Coconut Palm	150	218	368
EU	Eucalyptus	150	221	371
AP	Areca Palm	150	212	362
NF	Natural forest	150	287	437
	Total	1050	1554	2604

.

2.3. 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, phenology-based 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.

2.4. Methodology

2.4.1. 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 ① Resting, ② Recovery, ③ Vigorous, and ④ 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 double-stage 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 ①②③④ for short.

2.4.2. 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]:

$$J{M}_{ij}=2\left(1-e^{-B_{ij}}\right)$$

(1)

where the $B_{ij}$ is Bhattacharyya distance as

$$B_{ij}=\frac{1}{8}{\left({\mu }_i-{\mu }_j\right)}^T{\left(\frac{c_i+c_j}{2}\right)}^{-1}\left({\mu }_i-{\mu }_j\right)+\frac{1}{2}ln\left(\frac{1}{2}\frac{\det \left(c_i+c_j\right)}{\sqrt{\det \left(c_i\right)\det \left(c_j\right)}}\right)$$

(2)

where ${\mu }_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.

2.4.3. 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.

2.4.4. 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.

3. Results

3.1. 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.

3.2. Separability between Rubber Trees and Other Tree Species

Table 2. Class separability between rubber trees and other tree species based on J m distance.

Tree species	LC	CA	CP	EU	AP	NF
①	1.50	1.98	1.96	1.21	1.62	1.70
②	1.58	1.97	1.98	1.57	1.67	1.77
③	1.47	1.98	1.97	1.41	1.63	1.74
④	1.59	1.98	1.98	1.41	1.58	1.68
①②	1.88	1.99	1.99	1.84	1.93	1.96
①③	1.85	1.99	1.99	1.80	1.91	1.93
①④	1.86	1.99	1.99	1.84	1.89	1.95
②③	1.86	1.99	1.99	1.87	1.90	1.94
②④	1.88	1.99	1.99	1.89	1.88	1.96
③④	1.87	1.99	1.99	1.86	1.86	1.93
①②③	1.96	1.99	1.99	1.97	1.97	1.99
①②④	1.97	1.99	1.99	1.97	1.98	1.99
①③④	1.96	1.99	1.99	1.96	1.97	1.99
②③④	1.96	1.99	1.99	1.96	1.98	1.99
①②③④	1.99	1.99	1.99	1.99	1.99	1.99

. LC: Lychee, CA: Casuarina, CP: Coconut Palm, EU: Eucalyptus, AP: Areca Palm, NF: Natural Forest. ①: Resting stage, ②: Recovery growth stage, ③: Vigorous growth stage, ④: Slowdown growth stage.

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 ①③ (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.

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.

3.3. 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. Accuracy assessment of rubber plantation mapping results.

Tree	LC	CA	CP	EU	AP	NF	RT
Species	EC/EO	EC/EO	EC/EO	EC/EO	EC/EO	EC/EO	TP	PA%	UA	F1
①	10/5	1/0	1/1	18/26	14/8	14/12	142	71.00%	73.10%	0.72
②	6/5	0/1	0/0	16/13	14/14	9/8	155	77.50%	79.10%	0.78
③	13/9	1/1	1/2	18/15	11/11	6/2	150	75.00%	78.90%	0.77
④	14/5	2/0	1/1	22/19	6/9	10/16	145	72.50%	73.20%	0.73
①②	8/5	0/0	1/1	12/13	14/14	8/11	157	77.50%	78.10%	0.78
①③	8/3	0/0	2/1	14/12	14/16	7/9	156	78.00%	78.80%	0.78
①④	12/2	1/0	1/0	22/14	11/13	5/11	148	74.00%	78.70%	0.76
②③	9/3	0/0	0/0	15/14	15/13	7/7	154	77.00%	80.60%	0.79
②④	9/3	0/0	0/0	14/17	14/12	6/6	157	78.50%	80.50%	0.79
③④	15/6	0/0	1/0	18/17	6/4	10/5	150	75.00%	82.40%	0.79
①②③	7/3	0/0	0/1	13/14	16/12	7/5	157	78.50%	81.80%	0.80
①②④	8/3	0/0	0/0	12/15	16/10	6/3	158	79.00%	83.60%	0.81
①③④	10/2	0/0	0/1	20/12	10/10	7/8	153	76.50%	82.30%	0.79
②③④	8/2	0/1	0/0	13/15	16/8	8/2	155	77.50%	84.70%	0.81
①②③④	7/3	0/0	0/0	13/12	12/12	6/3	162	81.00%	84.40%	0.83

EC: errors of commission, EO: errors of omission, TP: true positive, PA: producer’s accuracy, UA: user’s accuracy, F1: F1 score; LC: Lychee, CA: Casuarina, CP: Coconut Palm, EU: Eucalyptus, AP: Areca Palm, NF: Natural Forest, RT: Rubber tree; ①: Resting stage, ②: Recovery growth stage, ③: Vigorous growth stage, ④: Slowdown growth stage.

. 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 (①②③④) 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 (①②④), 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 (②④), 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.

As to the four mono-stage dataset groups, the recovery growth stage (②) performed best with a UA, PA, F1, and true positive of 77.5%, 79.1%, 0.78, and 155, while the resting stage (①) performed worst with a UA, PA, F1, and true positive of 71%, 73.1%, 0.72 and 142, respectively.

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.

3.4. Rubber Plantation Map and Statistical Data Validation

The result of the rubber plantation map based on the all-stage combination dataset group (①②③④) 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. Accuracy assessment at the county level using statistical data.

ID	County	Our Work	Yearbook	Difference	Accuracy	Area Weight	Area-Weighted
		(km2)	(km2)	(km2)	(%)		Accuracy
3	Danzhou	1055.07	870.84	184.23	78.84	0.10	7.88
18	Qiongzhong	595.91	559.12	36.79	93.42	0.08	7.47
13	Baisha	579.16	632.81	−53.65	91.52	0.06	5.49
11	Chengmai	483.52	508.47	−24.95	95.09	0.06	5.71
15	Ledong	356.14	315.88	40.26	87.25	0.08	6.98
12	Lingao	322.54	219.91	102.63	53.33	0.04	2.13
5	Qionghai	314.61	350.46	−35.85	89.77	0.05	4.49
10	Tunchang	279.10	364.67	−85.57	76.53	0.04	3.06
7	Wanning	228.17	265.62	−37.45	85.90	0.06	5.15
14	Chanjiang	205.56	153.60	51.96	66.17	0.05	3.31
17	Baoting	187.25	224.79	−37.54	83.30	0.03	2.50
4	Wuzhishan	186.59	161.38	25.21	84.38	0.03	2.53
9	Dingan	174.33	205.62	−31.29	84.78	0.04	3.39
2	Sanya	156.48	122.54	33.94	72.30	0.05	3.62
1	Haikou	136.16	161.38	−25.22	84.37	0.06	5.06
8	Dongfang	94.62	99.89	-5.27	94.72	0.07	6.63
16	Lingshui	81.22	53.52	27.70	48.24	0.03	1.45
6	Wenchang	37.04	46.22	−9.18	80.14	0.07	5.61
	Total	5473.47	5316.72	156.75	97.05		82.47

. 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², which is about 2.93% higher than the 2020 statistical data (5316 km²). 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. Feature importance for S2 bands and rubber tree phenological phases.

		Resting	Recovery	Vigorous	Slowdown	Sum
Visible	B2	17.64	16.65	7.67	14.01	55.97
	B3	13.74	16.03	8.23	11.13	49.14
	B4	13.62	13.79	4.97	12.74	45.11
Red edge	B5	8.12	8.39	4.53	7.26	28.29
	B6	3.91	9.64	13.72	8.52	35.79
	B7	3.08	11.09	12.98	5.53	32.69
NIR	B8	3.4	8.92	11.06	5.18	28.55
	B8A	3.77	9.15	10.44	4.65	28.01
SWIR	B11	4.96	8.55	14.21	13.28	41.02
	B12	4.46	7.67	10.95	9.63	32.71
Sum		76.71	109.88	98.76	91.92

, the accuracy at the county level was calculated as follows:

$$\mathrm{Accuracy}=\left(1-\frac{abs\left(Difference\right)}{Yearbook}\right)\times 100$$

(3)

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.

3.5. Feature Importance

Table 5 shows the feature importance for the rubber tree discrimination that employed all four phenological stages (dataset group ①②③④). 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.

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.

4. Discussion

4.1. 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 three-monthly 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.

4.2. 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.

4.3. 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].

4.4. 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.

5. 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. Rule-based 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.