1. Introduction
Due to an increase in the concern regarding natural forest management and conservation, the rubber plantations in Xishuangbanna have gained worldwide attention [
1,
2]. Since the 1990s, rubber trees planted in this area have expanded rapidly due to the increasing price of rubber. As reported by a recent study, rubber plantations are being pushed to higher elevations (approximately 1400 m.a.s.l. in 2010) in Xishuangbanna, China [
3]. The rubber planting area in Yunnan Province (approximately 0.571 million ha.) exceeded that in Hainan (approximately 0.542 million ha.) and became the largest rubber plating area in China in 2015 [
4]. Xishuangbanna contains 75% of the rubber plantations in Yunnan [
5], and so therefore detailed, accurate, and timely mapping of the rubber plantations in Xishuangbanna is vitally important to satisfy the demands of rubber growers, marketers, and policymakers.
There are several methods to map rubber plantations. However, the defoliation of rubber trees during winter is a unique phenomenon for the vegetation in subtropical areas, and this phenomenon can be remotely sensed by using both optical and SAR (synthetic aperture radar) data. Consequently, phenological feature-based methods are the most commonly used methods [
6,
7]. According to the remote sensing data used, phenological feature-based rubber classification methods can be classified into two categories: methods that utilize SAR data and those that utilize optical remote sensing images. For the methods that utilize optical remote sensing images, the specific reflectance of monoculture rubber plantations and their changes could be captured by high revisit frequency remote sensors, such as MODIS (MODerate-resolution Imaging Spectroradiometer) and AVHRR (Advanced Very High Resolution Radiometer) [
8,
9]. However, the relatively coarse spatial resolutions of these time series data make it difficult to delineate and map rubber plantations in fragmented landscapes in tropical regions with high accuracy [
7]. Studies have demonstrated the promising utility of optical remote sensing data with high temporal resolution for differentiating rubber trees from other vegetation [
10]. It is even possible to map the age of a rubber plantation in Xishuangbanna using aggregated Landsat TM (Thematic Mapper) and ETM+ (Enhanced Thematic Mapper Plus) images [
11,
12,
13]. The object-level phenological features of a rubber plantation captured by Landsat NDVI data with 30 m spatial resolution were proven to be better than pixel-based rubber classifications [
14]. However, rubber trees in Xishuangbanna are planted on slopes with complex vegetation species. Considering the fragmented landscape and complex vegetation planting structure that may give rise to mixed pixels, it is challenging to use phenological features derived from time series optical satellite images with coarse spatial resolution to delineate and map the rubber plantations in this subtropical region [
6,
15]. Aerial photography and low-altitude sensing platforms (such as unmanned aerial vehicles, and maned airplane, balloon) can provide high-resolution data [
16,
17], but their endurance limits their applications in the identification of rubber plantations for larger areas. Another constraint for land cover classification in tropical areas is the frequent presence of clouds and fog cover, which makes it difficult to construct continuous mid- or high-resolution time series data with reliable quality [
18]. Although the approach presented by [
19] can monitor the expansion of rubber plantations by detecting and analyzing the shapelet structure in a Landsat-NDVI time series, and this method was proven to work well for the non-defoliation rubber identification, the 16-day revisit cycle and heavy cloud (and fog) contamination to optical remote sensing images remains a challenge for Landsat data to capture the accurate phenological change to identify defoliate rubber. Thus, a further study on using dense Landsat-like data to identify defoliation rubber is still missing.
In comparison to optical sensors, SAR can penetrate clouds and fog, and thus has advantages when mapping tropical forests. B. Chen integrated SAR and multitemporal Landsat imagery to extract the rubber plantations on Hainan Island of China, and the results demonstrated a great potential of using SAR data to map rubber plantations with high accuracy [
18]. Another strategy used SAR data to identify forests according to the different HH (Horizontal transmitting and Horizontal receiving) and HV (Horizontal transmitting and Vertical receiving) polarizations between forest and other vegetation types. The forest areas were then classified into rubber trees and other tree types following the phenological features derived from time series vegetation indices (such as NDVI, EVI, LSWI (land surface water index), and LAI (leaf area index)). These indices are commonly calculated using optical images such as MODIS, Landsat, or HJ-1A/B [
9,
20,
21]. Although SAR data are promising for the continuous acquisition of data in this frequently fog or cloud covered area, the fusion of optical data and SAR data, as well as the rugged terrain, remain obstacles to the rubber classification methods that utilize SAR data.
As the methods that utilized SAR and optical data presented promising results in terms of rubber classification, methods that integrate data from multiple sensors could be used to derive phenological and spectral reflectance features from vegetation, which would improve the accuracy of rubber identification [
22]. Other constraints to rubber classification in Xishuangbanna are terrain characteristics such as slope, aspect, and elevation [
14]. With the development of high spatio-temporal resolution optical data fusion algorithms, it is feasible to construct high spatio-temporal resolution optical data to identify rubbers [
23,
24]. Thus, this paper aims to blend high spatio-temporal resolution optical data using Landsat and MODIS data, and 5 indices (NDVI, EVI, NDMI, NBR, and TCA) will be derived based on this blended data set. Each of these indices will be used to classify rubbers to analyze how the terrain factors (TFs) affect the accuracy of rubber identification in Xishuangbanna and determine which of the blended high spatio-temporal indices would be the optimum choice for rubber classification.
3. Methods
There are mainly three parts in our proposed method: (1) the selection of optimum optical remote sensing window, (2) feature extraction based on the fused time series indices, and (3) classification of rubber trees and its accuracy assessment.
Rubbers are planted in tropical areas where the optical remote sensing images are frequently contaminated by clouds and fog. Xishuangbanna is a tropical monsoon forest climate region where rubber trees shed their leaves in February, and new leaves begin to turn green in late March. To specify the remote sensing window for the extraction of phenological features of rubber trees in this area, the analysis of the cloud coverage frequency is vitally important to the applicability of our method for rubber plantation identification. Therefore, we first used MOD035 data to analyze the average cloud coverage of the study area during four different years: 2000, 2005, 2010, and 2015, which range from June to June of the following year. This period completely covers the entire dry season of Xishuangbanna. There are many more clear days during the dry season than during the wet season. In addition, rubber trees shed their leaves during only this period when many cloud-free optical remote sensing images can be obtained. To improve the representativeness of the daily MOD035 data, we transformed the cloud coverage data to an average of every 7 days. The transformed cloud coverage data were then used to analyze the coverage frequency in Xishuangbanna to specify the proper remote sensing window.
Optical remote sensing images within the chosen remote sensing window were used to derive the changes in phenological features of rubber trees. The phenological features were derived based on the indices calculated from the fused data, including NDVI, EVI, NDMI, NBR, and TCA. Numerous vegetation indices can accurately reflect vegetation growth conditions and capture the trends in vegetation changes. Among these indices, the NDVI is a widely used vegetation index; however, there are some defects, such as its saturation to high vegetation coverage and disturbances caused by atmospheric noise and the soil background [
39]. In addition to NDVI, indices such as the EVI, NDMI, NBR, and TCA are commonly used to construct time series data for forest monitoring, and some of these indices have been proven to exhibit better performance over NDVI [
40]. This study used these 5 indices to construct the time series data, and the separability of these time series indices was compared and analyzed in rubber plantation classification, and through the fusion of time series data sets, we could obtain accurate annual phenological changes of rubber plantations and the minimum annual vegetation cover of rubber plantation could be obtained, that may be difficult to obtain from a single data (such as, Landsat data) derived phenology, the difference between our method and other methods are summarized in
Table 4.
This study finally classified the rubber plantations by applying two machine learning classifiers, the support vector machine (SVM) and random forest (RF), which were trained using the training samples. The SVM classifier is based on statistical machine learning theory to determine the location of decision boundaries that produce an optimal separation of classes. Thus, this nonparametric and sample distribution-free classifier was applied in this study. To estimate the performance of different time series indices, the RF classifier, which is another nonparametric statistical machine learning classifier, was also used. An RF is built based on the predictions from multiple decision trees, where each decision tree is constructed from a different bootstrap sample of the original data set [
41]. The features derived from the time series are high-dimensional data sets, and both SVM and RF classifiers are suitable for high-dimensional data classification. Therefore, the RF classifier was selected to be compared with the SVM classifier in this study. In the final step, the accuracies of different classification schemes were assessed based on the validation samples to determine the most suitable approach. The specific process is shown in
Figure 2, and the main parts are as follows.
3.1. Analysis of the Optical Remote Sensing Window for Xishuangbanna
This study analyzed the average cloud coverage of the study area from June to June of the following year in four different years: 2000, 2005, 2010, and 2015. Because there are as many as 365 images for each year, remarkable redundancy exists when analyzing the cloud coverage trends during different seasons, so the average cloud coverage was calculated over each 7 day period. The cloud coverage trends for each period are shown in
Figure 3. Results indicate that the lowest cloud coverage in the study area (as shown between the two red dotted lines in
Figure 3) commonly occurs from December (approximately DOY 327) to early March of the following year (approximately DOY 067). Thus, cloud-free optical remote sensing data are abundant during this period; in contrast, the highest cloud coverage in the study area appears from June to September, and more than 70% of the images are clouded (or fogged). Cloud-free MODIS images mostly occur during the dry season, and there are almost no completely clear images from June to September. Therefore, data between June and September are not suitable for classification of vegetation in Xishuangbanna using optical remote sensing time series data.
According to the field survey, the period from December to March of the following year belongs to the dry season in Xishuangbanna. With the decreases in temperature and precipitation, rubber plantations gradually enter the defoliation period; thus, vegetation indices such as EVI would reach the lowest value in February of the following year. Then, rubber plantations gradually turn green and new leaves began to grow, and the EVI value gradually increases. Therefore, from every November to March of the following year, the vegetation indices of rubber plantations are significantly different from those of other local vegetation. This period is also within the optical remote sensing window of Xishuangbanna, which indicates that the optimum time period for identifying rubber plantations using time series optical remote sensing data is from December to March of the following year.
3.2. Generation of Time Series Indices with High Spatial and Temporal Resolution
According to the optimum optical remote sensing window in Xishuangbanna, Landsat TM/ETM+/OLI, and MOD09GA data from December to March of the following year were used to fuse and extract the phenological characteristics of rubber plantations and their changes. The ESTARFM requires cloud-free and high-quality remote sensing data as inputs, so cloudless MOD09GA and TM/ETM+/OLI surface reflectance data during the remote sensing window were selected (using the MODIS QC data for MODIS, Fmask for Landsat images) as input data. Five indices were calculated based on the selected cloud-free surface reflectance data for MODIS and Landsat. Then, time series vegetation indices with high spatial and temporal resolution were generated by employing the ESTARFM algorithm. For each of the constructed time series indices, noise inevitably existed because of the directional reflectance of land cover, shadows caused by rugged terrain or atmospheric pollution. These noises are obstacles to time series feature extraction, so we employed a new noise-reduction algorithm based on Savitzky-Golay (S-G) called the Spatial-Temporal Savitzky-Golay (STSG) method to remove the noise, the new method assumes discontinuous clouds in space and employs neighboring pixels to assist in the noise reduction of the target pixel in a particular year [
42,
43]. This method is a convolution algorithm based on least squares. The main formula is as follows:
In Equation (1)
and
are vegetation indices before and after filtering; Ci is the least squares fitting coefficient of the high-order polynomial, which is the weight of the
i-th EVI value from the filter head; N is the length of the filter, its size is equal to 2m+1, and m is the window radius. The S-G filter can remove most of the noise and smooth the original time series data with good fidelity [
44]. Although a single S-G filtering fits well with the basic trend of the original data, noise still exists. Thus, a double S-G filtering (S-G filtering again based on the first S-G filtered results) was performed to remove additional noise, which can also fit the basic trend of the original data well [
45].
3.3. Analysis of Data Fusion Accuracy
3.3.1. Accuracy Analysis of the 5 Fused Time Series Indices
To assess the data fusion accuracy, it is needed to analyze the correlation between the five indices fused by ESTARFM and the corresponding indices calculated from the observed Landsat data of the same day. For example, to assess the accuracy of the ESTARFM fused indices on DOY 059 of 2015, the input MODIS and Landsat indices were DOY 003 and DOY 067 (time interval of 64 days), respectively. During the accuracy evaluation process, the correlation coefficient (R) between the fused index on DOY 059 and the index derived from the Landsat observation data was calculated. The evaluation index adopts the determination coefficient R2 and root mean square error (RMSE). The closer the determination coefficient R2 is to 1.0, the higher the accuracy of the fusion data is, and the closer it is to 0.0, the lower the accuracy of fused data is. Meanwhile, the closer the RMSE is to 0.0, the higher the accuracy of the fusion data is, and vice versa.
3.3.2. Analysis of Factors Affecting the ESTARFM Fusion Accuracy
In this study, we analyzed the following factors that may affect the ESTARFM fusion accuracy:
(1) The interval of days between input images. Time intervals of 32 days, 56 days, 64 days, and 112 days, this interval was depended on the available observed Landsat data.
(2) Terrain factors, including elevation, slope, and aspect as following:
(i) Based on DEM data analysis and field investigation, the data fusion accuracy was analyzed by dividing the elevation into seven ranges: less than 600 m.a.s.l., 601–800 m.a.s.l., 801–1000 m.a.s.l., 1001–1200 m.a.s.l., 1201–1400 m.a.s.l., 1401–1800 m.a.s.l., and greater than 1801 m.a.s.l.
(ii) The slopes were divided into 5 grades: 0–10°, 10–20°, 20–30°, 30–40°, and greater than 40°.
(iii) The slope aspect was divided into eight directions; namely, east (67.5–112.5°, semi-shady slope), southeast (112.5–157.5°, semi-sunny slope), south (157.5–202.5°, sunny slope), southwest (202.5–247.5°, sunny slope), west (247.5–292.5°, semi-sunny slope), northwest (292.5–337.5°, semi-shady slope), north (0–22.5°, 337.5–360°, shady slope) and northeast (0–22.5°, 22.5–67.5°, shady slope).
(3) Spatial heterogeneity index. Spatial heterogeneity is one of the main causes of spatial scale effects [
46] and may affect the fusion accuracy. To analyze the relationship between spatial heterogeneity and fusion accuracy, we used the spatial heterogeneity index (SHI) proposed by Zhang to quantitatively analyze the relationships between spatial heterogeneity and the three TFs. A SHI of a pixel
is defined as the sum of the absolute brightness difference between
and its eight neighboring pixels [
47]:
In Equation (2),
is the SHI of a pixel
The SHI for a region (assuming the image size is m multiplied by n) is defined as:
3.4. Comparative Analysis of Different Time Series Indices
Vegetation coverage is high in tropical areas, and the response of vegetation change in spectral space is limited and complicated. Therefore, it is necessary to compare and analyze the performance of different time series indices for rubber plantation changes in Xishuangbanna. In this study, 5 time series data sets (such as NDVI, EVI, NDMI, NBR and TCA) were selected for comparative analysis of their separability in rubber plantation classification.
3.4.1. Time Series Index-based Feature Extraction
Time series indices fused by the ESTARFM during the key phenophase of rubber plantations contain a large amount of information on the rubber growth state and changes, which are high-dimensional data. Excessive data dimensions may cause data redundancy and reduced classification accuracy. To better extract the time series features based on the fused data sets and remove redundant data, this study extracted six features based on the time series vegetation index curves using IDL (Interactive Data Language). Taking the EVI time series features as an example, as shown in
Figure 4, six classification features were extracted in this study, which were the maximum and minimum values of the time series index from December to March of the following year, the integral between the 1st day and the 45th day of the year, the integral of the overall time series index curve, the ratio of the maximum value of the time series index to the integral of the time series index curve, and the date when the minimum derivative of the time series curve could be obtained (detailed in
Table 5).
Figure 4 is a schematic diagram showing the time series EVI curve of the rubber plantations within the optical remote sensing window of Xishuangbanna. In this diagram, point a is the maximum EVI value within this remote sensing window, point c is the minimum EVI value, point b is where the minimum EVI derivative value is obtained, indicating the rate of change in deciduous rubber plantations; the time interval from a to c is the deciduous period of rubber plantations, and rubber plantations enter the regreening up period at point c. The integral of EVI corresponds to the biomass within this remote sensing window of Xishuangbanna (selected here from DOY 001 to DOY 045), which are apparently different between different vegetation types. The results of the six features extracted from the study area are shown in
Figure 5.
3.4.2. Separability Evaluation of Classification Features
To evaluate the separability of the 6 selected classification features, we need a monotonic criterion that exhibits a monotonic relationship with the classification error rate that can be measured [
48]. The most commonly used J-M (Jeffreys-Matusita) distance based on the conditional probability method was used in this study [
49]. When a feature follows the normal distribution, the J-M distance can be calculated as [
50]:
The
in formula (4) is Bhattacharyya distance [
51]:
In Equation (5),
is the mean of the
i-th category and
is the covariance matrix of the
i-th category. The range of J-M distance is [0,2]. J-M distances were calculated according to the selected sample points; a value of greater than 1.38 indicates good separability, values between 1.0 and 1.38 indicate moderately separable, and values less than 1.0 indicate poor separability [
52,
53]. Taking the EVI time series features of rubber plantations and natural forest as an example, the J-M distances are as follows:
(i) when 2 features were combined, the maximum J-M distance was 1.57, and the corresponding combination was IV and V;
(ii) when 3 features were combined, the maximum J-M distance was 1.80, and the corresponding combination was I, II and V;
(iii) when 4 features were combined, the maximum J-M distance was 1.87 and the corresponding combination was I, II, V and VI;
(ⅳ) when 5 features were combined, the maximum J-M distance was 1.93, and the corresponding combination was I, II, III, V and VI;
(ⅴ) when 6 features were combined, the maximum J-M distance was 1.98.
Then, we analyzed the J-M distances of these six features for NDVI, EVI, NDMI, NBR, and TCA time series indices, as shown in
Table 6. The results showed that when the features derived from the EVI time series index were the six features listed in
Table 5, the J-M distances of all types (natural forest, rubber plantations, other vegetation, and non-vegetation) were greater than 1.9, which indicates good separability between each of the classes, as shown in
Table 6. Therefore, these six features of the EVI time series index can not only remove data redundancy from the original time series curves but also ensure the classification accuracy.
3.5. Classification Schemes
This study analyzed which of the three data sets (monotemporal image, time series data set and selected feature data set) was the best classification data set. The classification schemes we designed are as shown in
Table 7. Each of the 3 data sets was classified by employing the SVM and RF classifiers.
5. Conclusions
This study presented a rubber plantation identification method based on vegetation phenology in cloudy and foggy tropical areas in this study. The method first analyzed the optical remote sensing window within the study area and then used spatio-temporal data fusion technology to construct a vegetation index data set with high temporal resolution within the remote sensing window. The phenological features of vegetation were captured to map rubber plantations using fused time series indices. The analysis of the presented method and several factors that may affect the high spatio-temporal resolution data fusion accuracy and rubber plantation classification show the following.
(1) Within the optimum optical remote sensing window in Xishuangbanna, fused optical remote sensing data with high spatio-temporal resolution could map the rubber distribution with high accuracy (overall accuracy of 89.51%, kappa of 0.86).
(2) Terrain factors, including the elevation, slope, and aspect, will have obvious negative effects on the accuracy of the fusion of optical data high spatio-temporal resolution: the fusion accuracy of each time series index fluctuated greatly at different elevations; the highest fusion accuracy occurred in areas with elevations from 1201–1400 m.a.s.l., and the lowest accuracy occurred in areas with elevations less than 600 m.a.s.l. The fusion accuracy decreased with the increase in slope, with the least impact on EVI and the greatest impact on NDVI; the slope aspect had less of an influence on the fusion accuracies of the five time series indices, but the fusion accuracy was lowest on the northwest slope. The SHI analysis indicated that land cover heterogeneity increased with decreasing elevation and increasing slope. The land cover types on the northwest slope tend to be heterogeneous compared to those on other slope aspects.
(3) Among the five commonly used time series indices (NDVI, EVI, NDMI, NBR and TCA), the EVI time series had a better ability to capture the different growing features of vegetation in this study. The highest classification accuracies were achieved for the best fusion results of EVI, which indicated that the EVI was the best choice among the five time series indices for time series features-used rubber classification.
The classification method presented in this study was mainly based on the phenological differences between different ground cover types; however, rubber plantations in Xishuangbanna are artificially planted vegetation, which would be affected by human activities to some extent. For example, some rubber trees may still have green leaves during the defoliation phase because of sufficient water and fertilizer conditions. These rubber plantations are more difficult to distinguish using the classification method based on phenological characteristics. The Unmanned Aerial Vehicles’ (UAV) remote sensing is a low-altitude remote sensing platform, which is not disturbed by atmospheric factors in the acquisition of images. It has the advantages that are inaccessible by traditional remote sensing technologies, such as low cost, simple operation, fast acquisition of images and high ground resolution. Therefore, the UAV-based method will become the hot spot of rubber plantations identification and change study in small area in the future. Furthermore, the method proposed by this study is a training samples-based supervised classification method in which the results and accuracies may largely depend on the training sample, which indicates the nonparametric method should be considered in further research. Another deficiency of our method would be the computation complexity, because it needs to fuse images, which is time and computing devices consuming. However, with the development of cloud computing services/devices, this may not be a barrier soon.
Through the fusion of time series data sets, the method can obtain accurate annual phenological changes of vegetation and minimum annual vegetation cover, which is difficult to obtain from a single data source (such as Landsat). This means the method is adaptable and can identify vegetation steadily in different areas. Therefore, this method can be used to identify other vegetation that have distinguishable phenological features in tropical areas. However, the extraction of non-defoliation rubber may perform better if we combine our method with the method proposed by Ye’s study [
19], because our method can obtain the rapid and unique growth process of rubber at higher temporal resolution. With the development of remote sensing satellites, higher spatial resolution data with a shorter revisiting period such as Sentinel-2 can be used to study the possible uncertainty in high spatio-temporal resolution data-based classification in future research. For the topographical factors of mountainous tropical areas, directional reflectance and the growth stage of vegetation are strong obstacles to rubber identification accuracy. Thus, it is also necessary to consider the correction of bidirectional reflectance caused by TFs. In addition, the solar zenith angle changes caused by shadows in mountainous areas are also significant influences that should be considered in future research.