High-Precision Stand Age Data Facilitate the Estimation of Rubber Plantation Biomass: A Case Study of Hainan Island, China

Abstract

Rubber (Hevea brasiliensis Muell.) plantations constitute one of the most important agro-ecosystems in the tropical region of China and Southeast Asia, playing an important role in the carbon budget there. Accurately obtaining their biomass over a large area is challenging because of difficulties in acquiring the Diameter at Breast Height (DBH) through remote sensing and the problem of biomass saturation. The stand age, which is closely related to the forest biomass, was proposed for biomass estimation in this study. A stand age map at an annual scale for Hainan Island, which is the second largest natural rubber production base in China, was generated using all Landsat and Sentinel-2 (LS2) data (1987–2017). Scatter plots and the correlation coefficient method were used to explore the relationship (e.g., biomass saturation) between rubber biomass and different LS2-based variables. Subsequently, a regression model fitted with the stand age (R² = 0.96) and a Random Forest (RF) model parameterizing with LS2-based variables and/or the stand age were respectively employed to estimate rubber biomass for Hainan Island. The results show that rubber biomass was saturated around 65 Mg/ha with all LS2-based variables. The regression model estimated biomass accurately (R² = 0.79 and Root Mean Square Error (RMSE) = 14.00 Mg/ha) and eliminated the saturation problem significantly. In addition to LS2-based variables, adding a stand age parameter to the RF models was found to significantly improve the prediction accuracy (R² = 0.82–0.96 and RMSE = 4.08–10.59 Mg/ha, modeling using samples of different biomass sizes). However, all RF models overestimated the biomass of young plantations and underestimated the biomass of old plantations. A hybrid model integrating the optimal results of RF and regression models reduced estimation bias and generated the best performance (R² = 0.83 and RMSE = 12.48 Mg/ha). The total rubber biomass of Hainan Island in 2017 was about 5.40 × 10⁷ Mg. The northward and westward expansions after 2000 had great impact on the biomass distribution, leading to a higher biomass density for the inland coastal strip from south to northeast and a lower biomass density in the northern and western regions.

1. Introduction

The rubber tree (Hevea brasiliensis Muell.), which is widely planted in tropical regions such as Southeast Asia and the tropical region of China, is of great importance in producing natural rubber and in providing timber when the rubber latex harvest is finished [1,2]. The establishment of new rubber plantations—usually converted from old rubber plantations, tropical forest, and cropland [3]—and timber harvesting in old rubber plantations will greatly alter the carbon budget [4]. Therefore, accurately estimating biomass, which is an important indicator of the carbon storage capacity and the potential carbon pool size of a forest ecosystem [5,6], is a widespread concern of ecologists, plantation managers, and policy makers related to the rubber industry.

The biomass of rubber plantations at a plot level is usually estimated by allometric equation (AE) models, linked with field measurements of the Diameter at Breast Height (DBH) and Height (H), or a combination of other factors, such as the stand age and wood density [7,8,9,10,11]. In China, four rubber biomass allometric models have been developed thus far [8,10,12,13]. All these models use DBH or the combination of DBH and H as independent variables, yielding R² values greater than 0.97. Three of these models were developed in ** the global aboveground biomass, forest height, and severe forest disturbances [30]. Generally, LiDAR is the top-ranked approach for forest biomass estimation as it can remove data saturation, followed by longer wavelength radar and optical remote sensing [31]. However, the data availability is the opposite; LiDAR data are the scarcest due to limited satellites carrying LiDAR instruments, and optical remote sensing data are the most abundant.

Many studies have used remote sensing to estimate the biomass of rubber plantations. The relationship between rubber biomass, spectral bands, and vegetation indices derived from Landsat Thematic Mapper (TM) was explored in **shuangbanna Prefecture, Yunnan Province, China, and different biomass models using stepwise regression [32,33] and Random Forest (RF) [34] were developed (R² = 0.43–0.67 and Root Mean Square Error (RMSE) = 30.27–46.05 Mg/ha). The total rubber biomass in **shuangbanna was about 1.97 × 10⁷ Mg, and the average biomass was 75.90 Mg/ha [32]. Yasen and Koedsin found that an artificial neural network based on WorldView-2 multispectral imagery could estimate biomass more accurately than stepwise multiple linear regression in Phuket Province, Thailand [2]. Charoenjit et al. [35] first estimated tree girth—using spectral bands, textural parameters, and vegetation indices from Thaichote multispectral imagery by a linear multiple regression model—and then used it with an AE model to estimate the rubber biomass in eastern Thailand. Avtar et al. [36] found that the Horizontal transmitting with Vertical receiving (HV) band of the Phase-Array L-band Synthetic Aperture Radar (PALSAR) positively correlated with the biomass of rubber plantations. The correlation was higher than for the Horizontal transmitting with Horizontal receiving (HH) band, but it began to saturate when the biomass was about 50 Mg/ha. Trisaongko and Paull estimated the biomass of rubber plantations with a very high accuracy (R² = 0.98 and RMSE = 1.88 Mg/ha), indicating that the acquisition time of PALSAR-2 and the choice of AE model significantly affect the accuracy [37]. By summarizing this literature, two main approaches for estimating the biomass of rubber plantations on a large scale can be obtained: (1) establishing direct relationships between biomass and spectral/SAR bands, vegetation indices, and textural parameters and (2) obtaining proxy factors such as DBH from remote sensing data and then using them as the parameters of AE models. In fact, all these approaches use canopy features from optical or SAR satellite imagery to estimate biomass. As the relationship between rubber biomass and canopy features gradually weakens with an increase in the stand age [38], the estimated biomass of middle-aged and old rubber plantations often exhibits large deviations.

There is also a close relationship between rubber biomass and stand age [39]. Good models (R² = 0.97–0.99) of the biomass/carbon stock of rubber plantations and stand age have been developed in **shuangbanna, Yunnan Province [8,40,41] and in Thailand [9]. Compared to DBH, the rubber stand age is easier to acquire accurately from time-series satellite imagery. Thanks to the free access of Landsat data policy [42], the establishment year (same as the stand age) of rubber plantations can be accurately traced back to the end of the 1980s at an annual resolution and RMSE of about 2 years [3,43], demonstrating great potential for biomass estimation. In addition to the stand age, the rubber biomass is also closely related to the geographical environment in which it grows, such as light conditions, soil moisture, and fertility. Differences in the growth condition of rubber plantations may affect the intensity of the spectral reflectance that can be used to estimate the biomass [35]. Therefore, we assume that combining the stand age and spectral characteristics can more accurately estimate the biomass of rubber plantations. The objectives of this study are (1) to explore the relationship between rubber biomass and stand age and LS2-based variables (e.g., spectral bands and vegetation indices); (2) to develop biomass models using stand age and LS2-based variables; and (3) to generate a biomass map of rubber plantations in Hainan Island, China and to investigate its spatial distribution pattern.

2. Materials and Methods

2.1. Study Area

Hainan Island (19°20′N–20°10′N, 108°21′E–111°03′E, Figure 1) is the second largest island in China and has an area of about 3.4 million ha. The topography of the island is mountainous in the interior, with lowlands along the coast. Wuzhi Mountain (1867 m AMSL) is the highest mountain on the island [44]. The climate varies from tropical to subtropical. The annual mean temperature is approximately 23–25 °C, and the monthly temperature varies between approximately 16 °C in January and approximately 30 °C from May to July. The average annual precipitation is 1500 to 2000 mm, about 80% of which occurs during the rainy season from May to October. However, the spatial distribution of precipitation can be as high as 2400 mm in the central and eastern area and as low as 900 mm in the coastal areas of the southwest [45].

Hainan Island is the second largest rubber plantation base in China. In the 1950s, a large area of rubber trees was planted to meet China’s national defense and economic needs. Statistics indicate that Hainan Island had about 5.43 × 10⁵ ha of rubber plantations by the end of 2017, accounting for more than a quarter of the total forest area, and is the largest plantation ecosystem there [44].

2.2. Data and Processing

2.2.1. Field Inventory Data

Field data consist of two parts: the main field data sampled in 2016 and supplementary data derived from [38] or collected from the experimental farm of the Chinese Academy of Tropical Agricultural Sciences (CATAS) from 2007 to 2014. The main field data are new and have a uniform acquisition time, so are used to explore the relationship between rubber biomass and remote sensing variables, to develop Random Forest (RF) models, and to evaluate their accuracy. Both the main field data and supplementary data were used to fit a biomass model using the stand age as the independent variable. The supplementary data were mainly collected for increasing samples in the 1990s—a period established with limited rubber plantations, which made it difficult to acquire wealthy field data by random sampling.

The main field data were collected in October of 2016 by randomly selecting rubber plantations with the aid of Google high-resolution satellite images. In each plantation, the DBH of each rubber tree at a 1.3 m height in an area of 25 × 25 m was measured, and the mean and standard deviation of the DBH were calculated. The stand age was obtained by interviewing farmers or was evaluated according to the height of tapped bark—a relatively fixed height of bark is consumed each year for latex harvest—and double checked with the plantation establishment year map of 2017 (described in Section 2.3). Field photos were taken by a Casio EX-H20G camera embedded with a Global Positioning System (GPS) receiver. The GPS information of the photos was extracted using these steps: (1) converting it into a Keyhole Markup Language (KML) file format by PhotoKML (v1.3) (https://www.softpedia.com/get/Multimedia/Graphic/Graphic-thers/PhotoKML.shtml), (2) importing the KML file in ArcGIS (v10.2, www.esri.com) and converting it into a shapefile format, and (3) extracting the latitude and longitude of each plot. A total of 51 rubber plantations were collected during the field campaign. The spatial distribution of field plots is shown in Figure 1 and descriptive statistics on the stand age and DBH are presented in

Table 1. Descriptive statistics of the stand age and Diameter at Breast Height (DBH) for the 51 samples acquired in 2016 (main field data set).

	Mean ± Std	Min	25th Percentile	50th Percentile	75th Percentile	Max
Stand Age	14.49 ± 6.69	5.0	10.0	13.0	16.5	30.0
DBH (cm)	17.83 ± 3.31	12.05	15.33	17.49	20.06	24.60

. Among the 18 counties where rubber trees have been cultivated, 16 of them have sampling plots.

The DBH and stand age of rubber plantations included in the supplementary data were collected in a similar way, but accurate GPS information of some plantations was not clear. Actually, the stand age is the primary concern when building a stand age-based biomass model. The supplementary data consist of 19 rubber plantations with different stand ages. The total biomass (below + aboveground) for a single rubber tree was calculated by the AE model (Equation (1)) [13]:

$$\ln B_{G\_tree}=2.4021\ln G-4.5803$$

(1)

where $B_{G\_tree}$ is the biomass (kg/tree) and $G$ is the girth at 1.3 m (cm). This model has a good accuracy (R² = 0.99, systematic error ≤ 1.469% [13]). On Hainan Island, rubber trees were generally planted at a density of 3 × 7 m or 4 × 5 m in columns and rows, resulting in about 495 trees/ha at the initial stage of plantation establishment. However, the tree density decreased with the year because of typhoons, which frequently hit the island and caused huge damage in the rubber industry [46]. Tree density data of more than 80 state farms (not distinguishing the stand age), from 1952 to 2001, were collected by Hainan State Farm Bureau [47], showing an average density of 381.45 trees/ha. Equation (1) was converted into Equation (2) at a hectare scale by multiplying the average tree density and by dividing it by 1000 (converting kg to Mg):

$$B_G=3.9105G^{2.4021}·{10}^{-3}$$

(2)

where $B_G$ is the biomass (Mg/ha). The girth and biomass of each plantation were stored in the attribute table of the shapefile. According to the size of biomass, the field data were split into eight groups (30, 40, 50, 60, 70, 80, 110, and 140 Mg/ha), which were further used as parameters for stratified k-fold cross-validation. Stratified k-fold cross-validation is similar to k-fold cross-validation, but it conducts stratified sampling instead of random sampling to ensure that relative class frequencies are approximately preserved in each train and validation fold. Then, the shapefile was uploaded to the Google Earth Engine (GEE, https://earthengine.google.com/) platform for subsequent analysis and modeling.

2.2.2. Satellite Imagery

Landsat 7/8 and Sentinel-2 multispectral imagery was used in this study. Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI) imagery were obtained from the U.S. Geological Survey (USGS, www.usgs.gov), and the Sentinel-2 MultiSpectral Instrument (MSI) imagery was obtained from the European Space Agency (ESA, https://sentinel.esa.int). All the top-of-atmosphere (TOA) reflectance data of Landsat 7/8 collection 1 Tier 1 and Sentinel-2 Level-1C were acquired during the day of year (DOY) between 120 and 330 (April 1 to November 30) and from 2016 to 2017. Surface Reflectance (SR) data were not used because Sentinel-2 SR (Level-2A) for the study area was not available at the GEE platform when we first carried out this study, and the mix use of Sentinel-2 TOA and Landsat SR may have the problem of cross-calibration between different sensors. Landsat 7/8 data have a 30-m spatial resolution, and Sentinel-2 MSI has a 10-m spatial resolution. These TOA reflectance data had already been orthorectified with a high geometric accuracy and reduced radiometric uncertainties between different acquisition dates and different sensors, such as ETM+ and OLI. All these data were available on the GEE cloud computing platform and could be readily accessed.

The cloud and cloud shadow of Landsat ETM+ and OLI were masked using bitmask information from the associated quality assessment band, which was generated by the widely used C Function of Mask (CFMask) algorithm [48]. Google’s cloudScore algorithm, together with the Temporal Dark Outlier Mask (TDOM), were used for masking the cloud and shadow of Sentinel imagery [49,50]. Regression models [50] were used to minimize TOA reflectance differences among the sensors of Sentinel-2 MSI, Landsat 7 ETM+, and Landsat 8 OLI, using Landsat 8 OLI TOA reflectance as a benchmark.

A median value composite (similar to a maximum value composite [51], but the median value in each band over time was selected) was applied to Landsat 7/8 and Sentinel-2 (LS2) to generate 20-m resolution cloud-free imagery. Then, two commonly used vegetation indices—the Normalized Difference Vegetation Index (NDVI) [52] and Enhanced Vegetation Index (EVI) [53]—were calculated using Equations (3) and (4), respectively:

$$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{NIR}}+{\rho }_{\mathrm{red}}}$$

(3)

$$\mathrm{EVI}=2.5\times \frac{{\rho }_{\mathrm{NIR}}-{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{NIR}}+6\times {\rho }_{\mathrm{red}}-7.5\times {\rho }_{\mathrm{blue}}+1}$$

(4)

where ${\rho }_{\mathrm{blue}}$, ${\rho }_{\mathrm{green}}$, ${\rho }_{\mathrm{red}}$, and ${\rho }_{\mathrm{NIR}}$ are the blue, green, red, and near-infrared (NIR) bands of the composite image (

Table 2. Spectral bands and their wavelengths of Landsat 7/8 and Sentinel-2 used in this study.

	Blue	Green	Red	NIR
Landsat 7 ETM+	B1 (0.45–0.52 µm)	B2 (0.53–0.61 µm)	B3 (0.63–0.69 µm)	B4 (0.78–0.90 µm)
Landsat 8 OLI	B2 (0.45–0.52 µm)	B3 (0.53–0.60 µm)	B4 (0.63–0.68 µm)	B5 (0.85–0.89 µm)
Sentinel-2 MSI	B2 (0.44–0.54 µm)	B3 (0.54–0.58 µm)	B4 (0.65–0.69 µm)	B8 (0.77–0.91 µm)

). Previous studies indicated that principal component analysis (PCA) of Landsat images had stronger relationship with forest biomass than individual bands [31]. Therefore, PCA was computed based on the LS2 composite image.

2.2.3. Terrain Data

The Shuttle Radar Topography Mission (SRTM) digital elevation data (V3) at a resolution of 1 arc-second (approximately a 30-m spatial resolution) were provided by the National Aeronautics and Space Administration (NASA) Jet Propulsion Laboratory (JPL) [54]. The elevation was directly obtained from SRTM V3 data, and the corresponding slope was calculated using GEE’s Application Programming Interface (API).

2.3. Rubber Plantation Biomass Estimation

Most AE biomass models are built on DBH, H, or their combination [12,13]. However, the biomass of rubber plantations also displays a close relationship with stand age—the R² of the established biomass model exceeds 0.93 [8,18]. A stand age-based biomass model is more popular in the remote sensing community because the stand age can be estimated more accurately by dense time-series satellite images than DBH or H [3,43], which are more easily and precisely obtained on the ground. In addition to the stand age, factors such as geography and climatic will affect the growth, meaning that rubber plantations of the same age have different canopy characteristics and biomasses. We made a hypothesis that integrating the stand age and spectral signatures from satellite imagery acquired in the growing season (April to December) can improve the estimation accuracy of the biomass model. Two methods were used to estimate the biomass of rubber plantations, and their differences were compared. The first was a regression model solely fitted by stand age, and the second was a Random Forest (RF) model built on LS2-based variables and/or the stand age. To generate a wall-to-wall biomass map of Hainan Island, both models require a reliable stand age map of rubber plantations in advance.

2.3.1. Map** the Stand Age of Rubber Plantations

In a previous study, we developed a new algorithm to identify the establishment year (same as the stand age, where the stand age equals the target year minus the establishment year) of rubber plantations based on all Landsat images of Hainan Island [3]. Specifically, signatures that topsoil exposed during plantation establishment and canopy closure linearly increased during nonproductive periods (rubber seedling to mature plantations) were employed to identify the establishment year. Annual images of the minimum and average Land Surface Water Index (LSWI) from 1987 to 2015 were respectively composited from all Landsat images (about 2000 scenes in total) acquired in the growing season of the corresponding year and were then used to generate an establishment year map of rubber plantations for Hainan Island in 2015. The establishment year map at an annual resolution since 1987 shows a very high accuracy (R² = 0.85/0.99 and RMSE = 2.34/0.54 year at a pixel/plantation scale). The establishment year established before 1987 was uniformly set as 1986 because the earliest collection of Landsat 5 imagery of Hainan Island started in 1987. We used the same algorithm and thresholds to generate the stand age map in 2017, using the newly obtained map of rubber plantations in 2017 as a benchmark and adding annual composited images of 2016 and 2017 to the previous annual composited data set.

2.3.2. Regression Biomass Model

Both the main field data (51 plots) and supplementary field data (19 plots) were used to build a regression biomass model. The field data used were obtained from rubber plantations of less than 30 years old. It is difficult to accurately obtain a stand age of over 30 years from both remote sensing and ground surveys—the former case is because of a lack of time-series images for map** (the earliest collection of Landsat 5 images of Hainan Island began in 1987), while the latter case is due to difficulties in determining the stand age from blurred tap** bark of old trees. The biomass was grouped according to stand age, and their average and standard deviation values were calculated. The stand age was used as the independent variable, and the average biomass was the dependent variable. An “S”-shaped sigmoidal biomass model was built using an online curve-fitting utility named MyCurveFit (mycurvefit.com). The biomass estimated by the regression model was referred to as $B_{RG}$.

2.3.3. Random Forest (RF) Biomass Model

The RF machine learning algorithm is generally immune to data noise and overfitting [55] and has been widely used within the remote sensing community such as classification [56,57,58,59], biomass estimation [25,60], disease monitoring [61,62], and hazards detection [63]. The good performance of the RF model over other methods, such as maximum likelihood and single decision trees, has been demonstrated by many studies and was summarized by Belgiu et al. [64].

Three tests were conducted during RF model building: (1) biomass saturation (sensitivity to biomass is lost up to a point) with different variables from LS2 composite images, (2) model building using LS2-based variables, and (3) model building using LS2-based variables and the stand age. Based on previous studies on forest biomass [31], we selected variables such as B3, B4, B5, B7, PC1, PC2, PC3, NDVI, and EVI to construct an independent variable data set. A scatter plot and the Pearson correlation coefficient between different variables and biomass were used to assess the biomass saturation. We calculated multiple correlation coefficients using samples of different biomass sizes, starting from 40 Mg/ha, with 5 Mg/ha as the step until 80 Mg/ha.

Every possible variable combination by gradually adding a distinct variable from the independent variable set—B3 and B4, and B3 and B5 as combination examples—was used to build RF models. We also tested the model performance under conditions of with/without stand age as a fixed independent variable. We kept the best 10 RF models each time and finally compared their performance under different variable combinations. Considering the problem of biomass saturation, model building was also conducted for samples from different biomass groups (smaller than mean values of 50, 60, 70, 80, 110, and 140 Mg/ha, respectively). Due to the limited sample size, we did not perform model testing on samples with a biomass of less than 50 Mg/ha, even using k-fold cross-validation (K = 4 was used here). In order to improve the efficiency of the test, we downloaded the sample data containing LS2-based variables from GEE to a local computer and then used the scikit-learn (https://scikit-learn.org/stable/index.html) python package for model testing. The number of trees in the forest was set as 100 and random_state was set as 1. The biomass estimated by the RF model was referred to as $B_{RF}$.

2.4. Accuracy Assessment

Stratified k-fold cross-validation (StratifiedKFold) was used to evaluate the model accuracy. The 51 samples were split into four smaller sets (each set is called a fold) using grou** field (values close to 30, 40, 50, 60, 70, 80, 110, and 140 Mg/ha) for stratification. The model was trained using k-1 of the folds as training data and was validated with the remaining part. The performance measure reported by k-fold cross-validation was then the average of the values computed in the loop. Model performances were evaluated by the determination of coefficients (R²) and Root Mean Square Error (RMSE). The scatter of $B_G$ against $B_{RG}$ and $B_{RF}$ was plotted, respectively. Matplotlib (https://matplotlib.org) and Seaborn, which is a Python data visualization library based on matplotlib (https://seaborn.pydata.org), were used for plotting.

3. Results

3.1. Variable Sensibility to the Biomass of Rubber Plantations

The scatter plots and correlation coefficients of biomass against different spectral bands, vegetation indices, and principal components are presented in Figure 2. The reflectance of B3, B5, and B7 decreases with biomass when the biomass is lower than 65 Mg/ha and then increases slowly. Among them, B3 drops more significantly, followed by B7 and B5. PC1, PC2, PC3, NDVI, and EVI first increase with the increase in biomass, but there are almost no changes after the biomass is greater than 65 Mg/ha, except for PC2, showing an obvious decrease in a high biomass range. The relationship between B4 and biomass is the weakest, even if the biomass is divided into two parts at 65 Mg/ha. When the biomass is lower than 65 Mg/ha, its correlation coefficient with B3, B7, PC2, and NDVI is quite high among these variables, ranging from −0.27 to −0.59, −0.20 to −0.72, 0.37 to 0.64, and 0.45 to 0.62, respectively. When the biomass continues to increase, the correlation coefficient between it and all variables weakens.

3.2. Biomass Estimated by the Regression Model

A nonlinear symmetrical sigmoidal model, named the four parameter logistic (4PL) regression model, was used to establish a regression model between the biomass of rubber plantations and stand age (Equation (5), 95% confidence) (Figure 3):

$$B_{RG}=145.40-\frac{115.27}{1+{\left(\frac{x}{17.46}\right)}^{3.23}}$$

(5)

where $B_{RG}$ is the biomass of rubber plantations (Mg/ha) and $x$ is the stand age. The R² is 0.963, adjusted R² is 0.967, RMSE is 7.12 Mg/ha, and Sum of Square Error (SSE) is 911.12. A slight overestimation is observed in immature plantations for which the stand age is less than 5 years. The maximum biomass for this model is about 145.40 Mg/ha. When combined with the biomass model and latest stand age map, the biomass of rubber plantations on Hainan Island could be obtained. Validation by main ground reference data collected in 2016 shows that R² reaches 0.79 and that RMSE is 14.00 Mg/ha (Figure 4a). The scatter plot of $B_{RG}$ against $B_G$ demonstrates that most points distributed along the 1-to-1 line, but the model exhibits a slight overestimation when rubber biomass is less than about 80 Mg/ha (Figure 4a).

3.3. Biomass Estimated by RF Models

R² of the best RF models developed by variable combinations of B3, B4, B5, B7, PC1, PC2, PC3, NDVI, and EVI ranges from 0.38 to 0.81, and the RMSE is between 4.31 and 24.34 Mg/ha (Figure 5). The RF model trained with samples of less than or equal to group 50 (specifically, all samples with a biomass < 55 Mg/ha) holds the highest R² (0.81) and lowest RMSE (4.31 Mg/ha). R² decreases to about 0.65 when sample groups less than 60, 70, and 80 are used for modeling and then gradually decreases to about 0.38 when the samples of a higher biomass (groups 110 and 140) are added to the models. However, the RMSE of the RF models generally increases from 4.31 to 24.34 Mg/ha with the addition samples of a higher biomass.

When the stand age is a fixed independent variable, the prediction accuracies of RF models are significantly improved. R² ranges from 0.81 to 0.96, and the RMSE is between 4.08 and 10.59 Mg/ha (Figure 5). The accuracy of the RF model built with samples of less than and equal to the 50 group is the highest, where R² reaches 0.96 (0.15 higher than that without using the stand age) and RMSE drops to 4.08 Mg/ha. When the sample biomass is high, the increase in R² is more obvious. For example, R² of a model built with samples of less than or equal to the 140 group increases from 0.38 to 0.85 and the RMSE decreases from 24.34 to 10.59 Mg/ha. The RMSE of the RF models generally increases with the addition of higher biomass samples.

The scatter plot of the best RF model developed by stand age and LS2-based variables with k-fold cross-validation is shown in Figure 4b. The points are more densely distributed near the 1-to-1 line when the biomass of rubber plantations is less than about 65 Mg/ha, and most points are then more scattered afterwards. The R², RMSE, and slope of the fit line are 0.70, 16.41 Mg/ha, and 0.78, respectively. A slight underestimation is observed for plantations which exhibit a higher biomass.

Since the RF models have a higher prediction accuracy for plantations with a lower biomass (Figure 5) and the regression model presents a high accuracy for plantations with a higher biomass (Figure 4a), a hybrid model integrating the optimal results of these two models can better predict the accuracy. The scatter plot of $B_G$ against results obtained from the hybrid model (B_{RG + RF}) is shown in Figure 4c. The prediction results are more closely distributed around the 1-to-1 line. The slope of the fit line, R², and RMSE are 0.92, 0.83, and 12.48 Mg/ha, respectively.

3.4. The Spatial Distribution of Rubber Plantations and Biomass

The establishment year map of rubber plantations on Hainan Island by the end of 2017 is shown in Figure 6a. Immature plantations with a stand age of less than four years were excluded from this map because of larger map** uncertainties. Rubber plantations are densely distributed in northern and center-to-east regions. Danzhou city possesses the largest area of rubber plantations, whilst other counties, such as Baisha, Chengmai, Qiongzhong, and Qionghai, also have dense distributions. The establishment year map shows an interesting distribution: A considerable number of rubber plantations have been established in the western and northern regions since 2000, whereas the old rubber plantations are mainly distributed along the coastal strip from south to northeast.

The spatial distribution of the biomass of rubber plantations generated by the hybrid model is shown in Figure 6b. The spatial distribution pattern is quite consistent with the spatial distribution of the establishment year. The biomass of newly established rubber plantations in the western and northern regions is generally below 60 Mg/ha, while the biomass of aging plantations distributed along the coastal strip from south to northeast and some border regions such as Danzhou city and Baisha county usually exceeds 110 Mg/ha. The biomass of some rubber plantations discretely distributed in Danzhou city and counties of Baisha and Qiongzhong is greater than 90 Mg/ha.

The average biomass of rubber plantations varies in different counties or cities (Figure 7a). Wanning city possesses an average biomass of about 105 Mg/ha and tops the rankings, followed by the cities of Qionghai, Haikou, Wenchang, Sanya, and Wuzhishan, all of which have an average biomass of about 100 Mg/ha. Baoting county has an average biomass of about 95 Mg/ha, followed by the counties of Dingan, Qiongzhong, and Ledong, possessing an average biomass that is slightly greater than 90 Mg/ha. The average biomass of Danzhou city, representing the largest production base on Hainan Island, is about 80 Mg/ha. Changjiang and Lingao county exhibit the lowest values among the 18 counties or cities, displaying an average biomass of about 72 Mg/ha.

The total biomass of rubber plantations varies dramatically among different counties or cities (Figure 7b). Danzhou city possesses a total biomass of 9.77 × 10⁶ Mg, which is far higher than other counties or cities, followed by the counties of Baisha, Qiongzhong, Chengmai, and Qionghai city, all of which have a total biomass of rubber plantations greater than 4.50 × 10⁶ Mg. The total biomasses of rubber plantations in Tunchang county, Wanning city, Ledong county, Dinan city, and Lingao county are about 3.21 × 10⁶, 3.05 × 10⁶, 2.78 × 10⁶, 2.59 × 10⁶, and 2.45 × 10⁶, respectively, while the remaining counties or cities have values of less than 2.00 × 10⁶ Mg. Dongfang city and Lingshui county have the smallest biomass of rubber plantations, possessing values of 6.88 × 10⁵ and 6.37 × 10⁵ Mg, respectively. The total biomass of rubber plantations on Hainan Island in 2017 was 5.40 × 10⁷ Mg.

The average biomass of rubber plantations varies greatly at different elevations and slopes (Figure 8a,b). The highest average biomass is located at an elevation between 200 and 300 m, which is about 93 Mg/ha, followed by rubber plantations distributed in an elevation range from 300 to 400 m and from 400 to 600 m. The lowest average biomass is located at an elevation between 50 to 100 m, which is about 85 Mg/ha. In areas with a slope below 12 degrees, the average biomass of rubber plantations increases with an increase in the slope, gradually increasing from about 83 Mg/ha within five degrees to about 94 Mg/ha at regions of 8–12 degrees. After the slope exceeds 12 degrees, the average biomass of rubber plantations is maintained at about 96 Mg/ha.

The total biomass of rubber plantations on Hainan Island is mainly distributed in low elevation and flat areas (Figure 8c,d). At an elevation of 100 to 200 m, the total rubber biomass is the highest, exceeding 2.16 × 10⁷ Mg, followed by areas of 50 to 100 and 200 to 300 m, with values of about 1.31 × 10⁷ Mg and 8.95 × 10⁶ Mg, respectively. The total biomass in areas above a 600-m elevation is extremely low. The total biomass of rubber plantations gradually decreases with an increasing slope. In flat areas of less than five degrees, the total biomass is close to 2.28 × 10⁷ Mg, followed by areas of 5–8 degrees, for which the total biomass is 1.12 × 10⁷ Mg.

4. Discussion

4.1. Biomass Saturation with Different LS2-Based Variables

Many studies have reported that the saturation point between the biomass of rubber plantations and different variables from remote sensing was very low. In this study, the saturation point of rubber plantations with B3, B4, B5, B7, PC1, PC2, PC3, NDVI, and EVI from LS2-based composite images was about 65 Mg/ha, which was higher than in studies conducted in Cambodia [36] and Thailand [2]. Avtar et al. [36] found that the PALSAR HV band positively correlated with the aboveground biomass of rubber plantations but was saturated at 50 Mg/ha. Yasen and Koedsin [2] found that the aboveground biomass of rubber plantations was saturated with WorldView-2 images at 40 Mg/ha in Phuket Province, Thailand. The main reason for the higher saturation point we monitored was that it contains belowground biomass. If the underground biomass was deducted at a rate of 20% [13], the saturation point of the aboveground biomass was also around 50 Mg/ha, which was consistent with these studies. Since the most LS2-based variables were saturated when the rubber biomass was very low, it was difficult to accurately estimate the biomass of rubber plantations from these variables. For example, only one model we developed using samples less than 55 Mg/ha had an R² of 0.81. The R² of other models built on samples with a biomass of less than 85 Mg/ha was about 0.65 and decreased to 0.38 when all samples were used for modeling (Figure 5a).

4.2. Biomass Estimation Using Stand Age and Remote Sensing Parameters

The regression biomass model fitted by stand age had a high accuracy (R² = 0.96 and RMSE = 7.12 Mg/ha), which was consistent with previous studies [8,10,18]. The high accuracy of the model could be explained by the fact that rubber plantations have a shorter economic life cycle (about 25–30 years). The growth of biomass only became slow after more than 20 years [38], so the biomass had a close relationship with stand age. The maximum biomass from this regression model (145.40 Mg/ha) was similar to studies such as [9,65], reporting 63.70 and 65.10 Mg C/ha of carbon stock for rubber plantations with rotation periods of 38 and 25 years, respectively. However, this value was much higher than the value from Trisasongko’s model—the maximum biomass of a 30-year-old rubber plantation was about 70 Mg/ha [37]. In addition to stand age, the planting density also affects biomass accumulation greatly [10]. We have incorporated planting density in the regression model, according to 50 years of extensive surveys from Hainan State Farm Bureau.

Obtaining a high-precision stand age map was essential for estimating the biomass of rubber plantations over a large area. Currently, most studies use one or several satellite images to estimate stand age, obtaining a stand age map of different age groups (e.g., <7, 5–15, and >15) [66,67,68]. In this study, the stand age map of rubber plantations in 2017 was generated at annual steps using all LS2 images acquired from 1987 to 2017 and a novelty algorithm developed in a previous study [3]. With this reliable stand age map, the biomass map through the regression model of Hainan Island exhibited high accuracy (R² = 0.79 and RMSE = 14.00 Mg/ha), better than models without the stand age—for example, higher than the regression model (R² = 0.43) developed in Cambodia [36]; artificial network model (R² = 0.66 and RMSE = 11.97 Mg/ha) established in Phuket Province, Thailand [2]; and RF model (R² = 0.43 and RMSE = 34.05 Mg/ha) built in ** rubber plantations and their stand age [3]. This is why we used a single pixel instead of the average value of the 3 × 3 pixel window of the field data for modeling.

4.4.3. Tree Density

Due to being frequently affected by severe typhoons, the tree density of rubber plantations on Hainan Island gradually decreases as the stand age increases [46,70]. However, we could not establish a model of the tree density and stand age because of the lack of abundant field data. Instead, a fixed density of 381.45 trees/ha, according to field data collected by Hainan State Farm Bureau, was used to estimate the biomass of rubber plantations. This tree density was smaller than the 418.99 tree/ha value used in Thailand, which is rarely hit by typhoons [9]. Using a fixed tree density underestimates the biomass of young and middle-aged plantations but overestimates the biomass of aged plantations. However, the biomass of rubber plantations established before 1987 was underestimated due to the underestimation of stand age (uniformly set to 1986), which, to some extent, offsets the overestimated biomass of old rubber plantations.