Comparison of Different Important Predictors and Models for Estimating Large-Scale Biomass of Rubber Plantations in Hainan Island, China

Abstract

Rubber (Hevea brasiliensis Muell.) plantations are among the most critical agricultural ecosystems in tropical regions, playing a vital role in regional carbon balance. Accurate large-scale biomass estimation for these plantations remains a challenging task due to the severe signal saturation problem. Recent advances in remote sensing big data, cloud platforms, and machine learning have facilitated the precise acquisition of key physiological variables, such as stand age (A) and canopy height (H), which are critical parameters for biomass estimation but have been underutilized in prior studies. Using Hainan Island—the second-largest rubber planting base in China—as a case study, we integrated extensive ground surveys, maps of stand age and canopy height, remote sensing indicators (RSIs), and geographical and climate indicators (ECIs) to ascertain the optimal method for estimating rubber plantation biomass. We compared different inputs and estimation approaches (direct and indirect) using the random forest algorithm and analyzed the spatiotemporal characteristics of rubber plantation biomass on Hainan Island. The results indicated that the traditional model (RSIs + ECIs) had low accuracy and significant estimation bias (R² = 0.24, RMSE = 38.36 mg/ha). The addition of either stand age or canopy height considerably enhance model accuracy (R² = 0.77, RMSE ≈ 21.12 mg/ha). Moreover, incorporating the DBH obtained through indirect inversion yielded even greater predictive accuracy (R² = 0.97, RMSE = 7.73 mg/ha), outperforming estimates derived from an allometric equation model input with the DBH (R² = 0.67, RMSE = 25.43 mg/ha). However, augmenting the model with stand age, canopy height, or their combination based on RSIs, ECIs, and DBH only marginally improved the accuracy. Consequently, it is not recommended in scenarios with limited data and computing resources. Employing the optimal model, we generated biomass maps of rubber plantations on Hainan Island for 2016 and 2020, revealing that the spatiotemporal distribution pattern of the biomass is closely associated with the establishment year of the rubber plantations. While average biomass in a few areas has undergone slight decreases, total biomass has exhibited significant growth, reaching 5.46 × 10⁷ mg by the end of 2020, underscoring its considerable value as a carbon sink.

1. Introduction

Biomass, defined as the accumulation of dry matter per unit time or area, represents a fundamental quantitative attribute of forest ecosystems [1,2,3]. Fluctuations in biomass mirror alterations and influences stemming from forest, climate, and human activities, rendering it an essential parameter for assessing vegetation carbon stocks and examining shifts in forest ecosystem structure and function [4,5,6]. Consequently, the accurate description and estimation of biomass are vital to a wide array of scientific endeavors [7]. In light of societal advancement, resource utilization, and global climate change, the investigation of biomass and carbon stocks has undergone significant expansion, with related research consistently receiving substantial interest [8].

Numerous methods exist for estimating biomass, with the precision and spatial scale of these approaches varying significantly [9]. Traditional field measurement remains the most accurate estimation technique, and the established allometric equation (AE) often serves as a critical foundation for alternative methods [2,5]. However, this approach is both time-consuming and spatially constrained [5]. Owing to its high spatial, temporal, and spectral resolution, as well as its capacity for continuous, long-term scientific observations, remote sensing technology has become a widely employed tool for biomass estimation [10,11,12,13]. Broadly speaking, forest biomass estimation via remote sensing has evolved through three stages: (1) estimation solely based on optical, radar, or Light Detection and Ranging (LiDAR) data [14,15]; (2) effective incorporation of multiple types of satellite imagery, mitigating the limitations of solitary satellite imagery data [16]; (3) selection of optimal variables from multi-source satellite imagery and ancillary data [5], thereby enhancing model accuracy [12,17].

Rubber (Hevea brasiliensis Muell.), a tall tree native to the Amazon basin, constitutes the primary source of natural rubber [3,6,18]. In response to the burgeoning industrial and economic growth, the tree has been extensively cultivated in Southeast Asia and the tropical regions of China since the last century [3]. Prompted by the swift escalation of natural rubber prices in the 2000s, global rubber plantation coverage has undergone rapid expansion, rendering rubber plantations significant ecosystems within numerous tropical regions. Beyond providing natural rubber, rubber plantations fulfill the ecological roles commonly associated with typical forests, such as carbon sequestration [19,20]. Studies indicate that total carbon stocks within rubber plantations surpass those of many tropical forests and agroforestry ecosystems, contributing positively to climate change mitigation efforts [21]. Accurate estimates of rubber plantation biomass are essential for predicting natural rubber yield [22], analyzing intrinsic dynamics, and understanding the carbon cycle in the tropics [21,23].

Numerous studies have investigated biomass estimation in rubber plantations; however, the primary study areas have been concentrated on Hainan Island and in Xishuangbanna, China [24,25,26,27,28]. Methodological evolution has advanced from initial direct remote sensing inversion to the incorporation of vital physiological parameters, such as stand age [14,19,23]. Rubber plantations are typically established concurrently, resulting in uniform stand age, which correlates closely with biomass and can be estimated using dense time-series satellite imagery [20]. For example, the authors of [6] combined stand age with the random forest algorithm to estimate the biomass of rubber plantations in Hainan Island, achieving an R² of 0.82 to 0.96 and root-mean-square error (RMSE) of 4.08 to 10.59 mg/ha. However, stand age is not the most ideal parameter for biomass estimation, as rubber tree growth conditions vary across different geographical environments despite identical ages. Diameter at breast height (DBH) is the simplest parameter to measure for plantations, but its large-scale measurement and accurate remote sensing retrieval pose challenges. Fortunately, the recently launched Global Ecosystem Dynamics Investigation (GEDI) satellite can accurately obtain high-resolution forest canopy heights on a global scale, which can be utilized for precise forest biomass estimation [29]. For example, Nico et al. (2022) fused GEDI and Landsat images to obtain a global canopy height product with a resolution of 10 m for 2020 [30]; Sun et al. (2022) found high accuracy in aboveground biomass (AGB) estimation using GEDI canopy metrics (R² = 0.82, RMSE = 19.1 mg/ha) [31]. Consequently, we hypothesize that integrating stand age, canopy height, remote sensing indices (RSIs), and environment and climate indices (ECIs) may yield a more accurate rubber plantation biomass estimation. The specific aim of this study was to employ machine learning algorithms to compare the effects of different variable combinations (stand age, canopy height, RSIs, and ECIs) and estimation methodologies (direct and indirect estimates) on the accuracy of rubber plantation biomass estimation.

2. Materials and Methods

2.1. Study Area

Hainan Island, with a total area of approximately 3.4 million hectares, possesses a topography typified by peripheral lowlands and a central mountainous prominence. The region is subject to a tropical monsoon maritime climate, exhibiting annual mean temperatures between 23 °C and 25 °C. The island’s iso-precipitation line distribution indicates elevated concentrations in the central and eastern areas, whereas diminished concentrations are observed along the southwestern coast. Furthermore, Hainan Island is characterized by pronounced rainy and dry seasons, with the former extending from May to October and constituting approximately 80% of the overall annual precipitation.

Since the 1950s, extensive rubber plantations have been established on Hainan Island to promote economic, scientific, and technological advancements [27]. Hainan Island, once the largest rubber plantation area in China, has been surpassed by Yunnan province, and it now ranks second with an approximate area of 520,000 hectares. Rubber plantations on Hainan Island are predominantly located in the northwest, central, and northern regions, accounting for over 25% of the total forest area and forming the largest plantation ecosystem on the island [32,33].

2.2. Data and Processing

2.2.1. Field Inventory Data

Field data were procured from sampled plots across the island and long-term positioning monitoring plots located at the experimental farm of the Chinese Academy of Tropical Agricultural Sciences (CATAS). The field investigation took place in October 2016, during which time 70 rubber plantations with diverse stand ages on Hainan Island were selected for analysis. For each rubber tree within a 25 × 25 m zone in each plantation, DBH using a 1.3 m height was measured. Rubber stand age was determined through both field estimation and an algorithm analyzing time-series satellite imagery dating back to 1987 [34]. To correspond with the canopy height map in 2020 (elaborated in Section 2.2.3), these samples were further scrutinized using Google Earth high-resolution imagery, retaining only the sample plots with accurate geolocation and persistence in 2020. In a prior biomass estimation study, a mere 51 samples with precise geolocations and stand ages under 30 years were utilized [6]. Due to advancements in stand age monitoring technology employing satellite images, 58 samples remained accessible after excluding the rubber plantations harvested in 2020 (Figure 1). Among these samples, only a few overlapped with the GEDI footprint in 2020.

At the CATAS experimental farm, there are nine long-term positioning monitoring plots ranging from 1.39 to 2.98 hectares. From 2007 to 2014, the DBH of each tree was measured monthly during the seedling stage (≤4 years) and subsequently at quarterly or semi-annual intervals as trees grew taller. The DBH for the rubber plantations from 2007 to 2014 was calculated using Equation (1):

$$D_g=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{d}_i^2}}{N}}$$

(1)

where $D_g$ is the mean DBH, $d_i$ is the DBH of the $i\mathrm{t}\mathrm{h}$ tree in the stand, and $N$ is the total number of trees. The total biomass of each plantation ($B_{AE}$, mg/ha) was calculated using Equation (2), which was developed based on an AE model, taking into account initial planting density and tree loss due to natural disasters, such as typhoons [6].

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

(2)

where $G$ (cm) is the average girth at 1.3 m based on the field measurement.

2.2.2. Satellite Imagery

The U.S. Geological Survey (USGS) and National Aeronautics and Space Administration(NASA) provided the Landsat 5/8 optical satellite imagery, while the Japan Aerospace Exploration Agency (JAXA) contributed the Phased Array L-band Synthetic Aperture Radar (PALSAR/PALSAR-2) data. Both datasets underwent preprocessing and were implemented on the Google Earth Engine (GEE) cloud platform. The Landsat product employed was a collection 2 level-2 surface reflectance product that had been subjected to requisite preprocessing, including geometric and atmospheric correction. Clouds and shadows were masked using the quality band (QA_PIXEL) associated with each scene.

Three widely used vegetation indices, the normalized difference vegetation index (NDVI) [35], enhanced vegetation index (EVI) [36], and land surface water index (LSWI) [37], were calculated using Equations (3)–(5), respectively.

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

(3)

$$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)

$$LSWI=\frac{{\rho }_{NIR}-{\rho }_{SWIR1}}{{\rho }_{NIR}+{\rho }_{SWIR1}}$$

(5)

where ${\rho }_{\mathrm{blue}}$, ${\rho }_{red}$, ${\rho }_{NIR}$, and ${\rho }_{SWIR1}$ are the blue, red, near-infrared (NIR), and shortwave infrared (SWIR1) bands of the Landsat 5/8 imagery, respectively.

The 25 m PALSAR/PALSAR-2 imagery utilized consisted of an annual global mosaic composed of images acquired during the growing season. The mosaic data encompassed horizontal transmit horizontal receive (HH) and horizontal transmit vertical receive (HV) bands. Referring to prior analogous studies [38], the gray level co-occurrence matrix (GLCM) was calculated for the HV band. Subsequent to GLCM computation, 18 bands were derived, but a limited number of key variables were screened using correlation analysis [39]. The selection criterion entailed retaining only one band with a high correlation coefficient; ultimately, HV_savg and HV_shade were the final bands retained. The original optical spectral bands (blue, green, red, NIR, SWIR1, SWIR2) and the NDVI, EVI, and LSWI, as well as SAR variables (HH, HV, HV_savg, HV_shade), are collectively referred to as remote sensing indicators (RSIs).

2.2.3. Terrain and Climate Data

Elevation data possessing a spatial resolution of 1 arc-second (approximately 30 m resolution) were obtained from the National Aeronautics and Space Administration’s Jet Propulsion Laboratory (NASA JPL) and originated from the Shuttle Radar Topography Mission (SRTM) Version 3 (V3). Slope was calculated using the application programming interface (API) of the GEE. Furthermore, the annual average precipitation, with an approximate resolution of 11 km, was a composite derived from the fifth generation of the European ReAnalysis (ERA5) dataset, which represents an enhanced global dataset for the land component generated by the European Centre for Medium-Range Weather Forecasts (ECMWF) [40]. These variables, including elevation, slope, annual average temperature, and precipitation, are collectively referred to as environmental and climate indicators (ECIs).

2.3. Biomass Estimation

2.3.1. Work Flow

The workflow diagram is depicted in Figure 2, outlining a three-step process: (1) preparation of biomass inversion variables encompassing direct variables, such as RSIs and ECIs, and indirect variables, like canopy height (H) and DBH; (2) evaluation and comparison of various biomass models based on these variables; and (3) identification of the optimal model for estimating rubber plantation biomass in Hainan Island, accompanied by an analysis of its temporal and spatial characteristics. The canopy height for 2016 (the year of the island-wide field investigation) was estimated through the development of a height and stand age model. At present, there is no readily available remote sensing-based model for rubber DBH estimation. We modeled and compared field observations, stand age, RSIs, and ECIs, selecting the optimal model to estimate DBH.

2.3.2. Key Independent Parameter Estimation

Stand Age (A)

The 2020 and 2016 maps of rubber plantations were generated through the utilization of dense optical and radar satellite imagery coupled with a phenology-based decision tree algorithm. The implementation of these techniques yielded an overall classification accuracy of 90%, as reported by [33]. The stand age maps for the years 2016 and 2020 relied on these rubber plantation maps and a comprehensive analysis of long-term time-series images from Landsat and Sentinel-2 satellites beginning in 1987. This analysis followed an algorithm previously outlined [34]. The algorithm’s core principle involves utilizing distinctive features of land-use and land-cover changes during the establishment phase of plantations and the subsequent progression of canopy closure from rubber seedlings to fully mature dense plantations. By employing these parameters, the algorithm can effectively identify the precise year of establishment for rubber plantations. The target year minus the establishment year is the stand age of a rubber plantation in a specific year.

Canopy Height (H)

The canopy height data for the entirety of Hainan Island in 2020 were produced in [30]. This study is considered a milestone in the field as it produced the first global 2020 wall-to-wall canopy height map at 10 m spatial resolution through the integration of GEDI and Sentinel-2 imagery with a probabilistic deep learning model. Given the strong relationship between stand age and canopy height [41], it was judged to be feasible to predict the canopy height of rubber plantations in 2016 based on this relationship and the wall-to-wall canopy height map from 2020. The primary steps included: (1) reclassifying the 2020 stand age map of rubber plantations in 1-year increments, followed by masking of the 2020 canopy height map; (2) performing sampling for each age class map with 50 sample points and extracting the corresponding canopy height for these sample points; (3) computing the percentile mean and standard deviation for the canopy height within a range of 20 to 80 percent (excluding misclassified rubber pixels and canopy height outliers); and (4) constructing a regression model with stand age as the independent variable and average canopy height as the dependent variable.

DBH

The DBH is also closely related to the stand age and the growth environment of rubber trees. By combining the field DBH, stand age (A), canopy height (H), RSIs, and ECIs, we developed and compared three random forest (RF)-based DBH models using Equation (6).

$$DB{H}_{RF}=f(RSIs,ECIs,A|H)+\epsilon$$

(6)

where $DB{H}_{RF}$ is the DBH estimated by the random forest algorithm; $A|H$ represents the stand age and canopy height, which were additional independent variables during the test; and $\epsilon$ is the error in the estimated DBH. The Boruta algorithm was used to select optimal independent variables from the RSI and ECI sets for the DBH and subsequent biomass estimation, as the feature selection process is fundamental to any machine learning project [42]. The RF algorithm was selected for DBH estimation due to its general immunity to overfitting and data noise and its widespread use in the remote sensing community [43,44]. Ultimately, the optimal model was chosen to estimate the DBH for rubber plantations for different nominal years.

2.3.3. Biomass Estimation Models

The RF algorithm utilized several independent variables, including RSIs, ECIs, A, H, and DBH. The primary reason for selecting RSIs and ECIs as variable sets was their ease of acquisition and their extensive use in traditional biomass estimation via remote sensing. Moreover, the models were assessed by incorporating A, H, DBH, and their combinations (eight combinations in total; Figure 2) using Equation (7).

$$B_{RF}=f(RSIs,ECIs,A\&|H\&|DB{H}_{RF})+\epsilon$$

(7)

where $B_{RF}$ is the biomass estimated, and $A\&|H\&|DB{H}_{RF}$ means adding one or a combination of $A$, $H$, and $DB{H}_{RF}$ to the $RSIs$ and $ECIs$ during model estimation. The RF model’s hyperparameters were set to 100 for the number of trees and 1 for the random state. In addition to the RF models, biomass estimates based on the AE model ($B_{AE}$) using Equations (2) and (8) were compared, utilizing the best $DB{H}_{RF}$ as an input parameter.

$$B_{AE}=f(DB{H}_{RF})+\epsilon$$

(8)

2.4. Accuracy Assessment

Considering the limited availability of field data, this study adopted a k-fold (k = 4 was used herein) cross-validation approach for modeling and evaluation. Plots were stratified into diverse categories based on their biomass (10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, and greater than 110 mg/ha), which served as the grouping variable. The sample was subsequently divided into four smaller subsets referred to as folds, utilizing a hierarchical k-fold cross-validation method. This technique refines the conventional k-fold cross-validation method by incorporating a grouping variable, ensuring a roughly even distribution of biomass across folds. During each iteration, k−1 folds were employed for training, while the remaining fold was designated for validation. To evaluate model performance, the R² and RMSE were implemented as metrics. The performance of various band combinations was compared by calculating their mean values within the k-fold loop.

2.5. Spatial Analysis

The biomass of rubber plantations in different administrative regions (provinces and counties/cities) and topographic environments (with differing elevation and slope) on Hainan Island was examined using the GEE cloud platform. Elevation was segmented into seven groups with bin values of 0, 50, 100, 200, 300, 400, and 600 m. Slope was categorized into seven levels with bin values of 0, 5, 10, 15, 20, and 25°, slightly adjusted based on the technical standard Comprehensive control of soil and water conservation—General rule of planning (GB/T 15772-2008) [45]. Additionally, the spatiotemporal characteristics of biomass changes between 2016 and 2020 were analyzed.

3. Results

3.1. Estimation of Canopy Height and DBH

The canopy height (H)–stand age model is presented in Figure 3a. With increasing stand age over the first 10 years, the height of rubber trees rose dramatically. The growth rate slowed between 10 and 15 years, and after 15 years, the trees’ height increased very slowly. According to logistic curve fitting, the R² reached 0.95. Figure 3b–d show the DBH models established using three different methods. DBH estimated via RSIs and ECIs displayed a high degree of accuracy for rubber trees around 17 years but had clear tendencies for overestimation for trees younger than 15 years and underestimation for trees older than 20 years. The linear fit demonstrated an R² of 0.36 and an RMSE of 4.38 cm. When the canopy height or stand age was introduced, the prediction accuracy significantly improved: the scatter points for the observed and predicted DBH were mainly distributed along the 1-to-1 line, while the R² reached 0.76 and 0.77 and the RMSE decreased to 2.65 cm and 2.60 cm, respectively. Considering that stand age is easier to obtain using remote sensing than canopy height, we finally chose the model with RSIs, ECIs, and stand age as independent variables to estimate DBH.

3.2. Comparison of Different Biomass Estimation Models

The nine biomass models developed with different variables and field data from 2016 are presented in Figure 4. Notably, the model exclusively employing RSIs and ECIs exhibited the lowest accuracy in estimation (R² = 0.24 and RMSE = 38.36 mg/ha) along with the most substantial estimation bias, particularly for rubber plantations with biomass greater than 100 mg/ha (Figure 4a). However, with the introduction of stand age, a significant improvement in prediction accuracy was observed: R² increased to 0.77 and RMSE decreased to 21.16 mg/ha. Additionally, scatter points demonstrated a closer alignment along the 1-to-1 line (Figure 4b). Likewise, when solely incorporating canopy height, the model’s prediction accuracy remained comparable to that achieved with stand age (R² = 0.77 and RMSE = 20.91 mg/ha), while the slope of the linear fit improved from 0.79 to 0.80 (Figure 4c). However, the inclusion of DBH led to a substantial enhancement in model accuracy: R² increased from 0.77 to 0.97, RMSE decreased to 7.73 mg/ha, and the slope improved from 0.80 to 0.93 (Figure 4d). Furthermore, when simultaneously incorporating two out of three variables (A, DBH, and H), the prediction accuracy of both the models incorporating DBH surpassed that of the model without DBH. The R² for the former models reached 0.97, while the latter model lagged at a value of 0.77 (Figure 4e–g). Incorporating DBH, A, and H together yielded an identical R² of 0.97 but with a slightly reduced RMSE of 7.43 mg/ha. However, it is worth noting that the biomass calculated using the AE and DBH (from the best RF model) exhibited lower precision compared to the RF models, with an R² of 0.67 and an RMSE of 25.43 mg/ha (Figure 4i). Consequently, the RF model utilizing RSIs, ECIs, DBH, A, and H as independent variables was selected as the final biomass model. Utilizing this model, the biomass of rubber plantations in both 2016 and 2020 was determined.

3.3. Spatiotemporal Pattern for Rubber Biomass at City/County Scale in Hainan Island

Figure 5 provides an overview of the spatial distribution of rubber plantation biomass in 2016 and 2020. In 2016, the rubber plantations with a biomass of about 70 mg/ha were mainly located in the northwestern and northern regions of Hainan Island, including Lingao and Changjiang county and northwest of Danzhou city. On the other hand, the central region and southeastern coastal regions had higher biomass for the rubber plantations, typically ranging between 70 mg/ha and 85 mg/ha (Figure 5a). In 2020, there was a significant increase in both the area and biomass of rubber plantations. The rubber acreage increased throughout the island, especially in the northwest and northern regions. Rubber plantations with a biomass of approximately 65 mg/ha or less were found to be more widely distributed in the western region of Danzhou city, the northern region of Changjiang and Ledong county, Haikou city, Ding’an county, and Wenchang city. In most of the remaining areas, the biomass of rubber plantations increased to around 70 mg/ha. The areas with the highest biomass were often distributed in the border areas of different cities and counties, such as the border areas between Chengmai county and Tunchang county, Qiongzhong county, Wuzhishan city, Baisha county, and Wanning city, where the biomass of many of the rubber plantations exceeded 120 mg/ha (Figure 5b).

The average biomass of rubber plantations across different cities and counties changed greatly between 2016 and 2020 (Figure 6a). In 2016, the average biomass was approximately 75 mg/ha. Wuzhishan city, Wanning city, and Sanya city had the highest average biomass of about 84 mg/ha. Conversely, the average biomass in Changjiang county, Dongfang county, Ledong city, and Dingan city was relatively low, approximately 69 mg/ha. With the exception of Baisha, Chengmai city, Danzhou county, and Qiongzhong, the average biomass of rubber plantations in the other 14 cities and counties showed a decreasing trend. Notably, Chengmai city had the highest increase of 5.70 mg/ha, far exceeding the other cities and counties. The cities and counties with remarkable decreases in biomass were Wanning county, Tunchang city, Lingao city, and Dongfang county. In 2020, Qiongzhong county had the highest average biomass of approximately 85 mg/ha, followed by Wuzhishan city and Chengmai county, Baisha city, and Wanning county, which had average biomass exceeding 79 mg/ha. The average biomass in Sanya city, Baoting county, Danzhou city, Lingao county, and Tunchang county was concentrated around 73 mg/ha. Dongfang city and Changjiang county had the lowest average biomass for rubber trees among the 18 cities and counties, with values of 63.81 mg/ha and 63.57 mg/ha, respectively.

Figure 6b demonstrates that the total biomass of the rubber plantations in different cities and counties on Hainan Island varied significantly. Danzhou city had the highest total biomass in both 2016 and 2020, with a value of 1.06 × 10⁷ mg in 2020, which was significantly higher than that of other cities and counties. The growth of rubber plantation biomass in Danzhou city was also very rapid, increasing by about 1.93 × 10⁶ mg. In 2016, the total biomass of rubber plantations in Baisha county, Chengmai county, and Qiongzhong county exceeded 4.00 × 10⁶ mg, with the first two cities and counties also demonstrating significant increases. Additionally, Ledong county, Qionghai city, Tunchang county, and Lingao county had total rubber plantation biomass exceeding 2.00 × 10⁶ mg in 2016. Conversely, Dongfang county and Linshui city had smaller total biomass values of 0.63 × 10⁶ mg and 0.67 × 10⁶ mg, respectively, in 2020. The total biomass of rubber plantations on Hainan Island in 2016 and 2020 was 4.52 × 10⁷ and 5.46 × 10⁷ mg, respectively.

3.4. Topographic Pattern of the Latest Rubber Biomass on Hainan Island

Rubber plantations exhibited notable variability in average biomass across differing elevations and slopes in 2020, as demonstrated in

Table 1. The average and total biomass of rubber plantations at different elevations and slopes on Hainan Island in 2020.

Elevation (m)	0–50	50–100	100–200	200–300	300–400	400–600	>600
Biomass (mg/ha)	69.08	75.17	79.02	84.48	85.74	86.52	82.34
Biomass (×106 mg)	4.83	14.05	23.3	9.7	3.95	1.85	0.06
Slope (°)	0–5	5–10	10–15	15–20	20–25	>25
Biomass (mg/ha)	74.48	78.41	83.33	86.09	87.54	89.38
Biomass (×106 mg)	24.09	16.35	8.13	4.84	2.72	1.61

. The most substantial average biomass was observed in rubber plantations situated at elevations between 400 and 600 m above sea level, which exhibited an average biomass of 86.52 mg/ha. Comparatively, plantations distributed at elevations between 300 and 400 m and 200 and 300 m above sea level displayed average biomasses of 85.74 mg/ha and 84.48 mg/ha, respectively. In contrast, the rubber plantations situated at elevations between 0 and 50 m above sea level displayed the lowest average biomass, with a value of 69.08 mg/ha. The average biomass fluctuated slightly among the various slope regions. Specifically, plantations situated in the flattest areas (<5°) exhibited an average biomass of 74.48 mg/ha, while the highest average biomass, 89.38 mg/ha, was observed in plantations situated in regions with slopes greater than 25°. The average biomass of rubber plantations increased marginally with increases in slope, from 74.48 mg/ha in plantations with slopes less than 5° to 89.38 mg/ha in plantations with slopes above 25°.

This study investigated the total biomass of rubber plantations situated at varying elevations and slopes. The findings reveal that the plantations at elevations of 100–200 m above sea level exhibited the highest total biomass, exceeding 2.3 × 10⁷ mg. Rubber plantations at elevations of 50–100 m and 200–300 m possessed total biomass of 1.41 × 10⁷ mg and 0.97 × 10⁷ mg, respectively. Conversely, plantations at elevations above 600 m had considerably lower total biomass amounting to 0.6 × 10⁵ mg. Furthermore, the total biomass of the rubber plantations progressively diminished with increasing slope. For instance, in the plains with slopes less than 5°, the total biomass approached 2.41 × 10⁷ mg, but it decreased sharply to 1.64 × 10⁷ mg in areas with slopes of 5–10°. Notably, regions with slopes exceeding 25° exhibited significantly lower total biomass for rubber plantations amounting to 1.61 × 10⁶ mg.

4. Discussion

4.1. Biomass Estimation Using Different Variables

Signal saturation is a common problem in estimating forest biomass using optical and SAR remote sensing data [12,46,47]. It has been reported that data saturation may occur when biomass density reaches 100–150 mg/ha in moist tropical forests, depending on the complexity of the forest stand structures caused by biophysical environments [12]. For uniformly planted rubber plantations, many studies have found that the saturation point for the biomass is lower than that in moist tropical forests. For example, saturation occurred at 50 mg/ha with the PALSAR HV band, 40 mg/ha with the WorldView-2 images, and 65 mg/ha with the LS2-based spectral indices [6,48,49]. Therefore, it is difficult to accurately estimate the biomass of rubber plantations solely using traditional predictors, such as spectral bands, vegetation indices, geo-environmental indicators, or their combinations. One of our models demonstrated the a estimation bias using RSIs and ECIs: the biomass of rubber plantations was seriously underestimated when it was greater than about 70 mg/ha, with an R² of 0.24 and RMSE of 38.36 mg/ha (Figure 4a).

The rate of forest biomass accumulation is closely related to stand age, especially in the transitional stage of forest succession [50,51,52]. For specific species, such as rubber trees, the economic life cycle is about 30 years. During the growth and development stage (<10 years), the biomass increases rapidly with the increase in the stand age and gradually slows down with the increase in latex harvesting time [53]. The stand age of rubber trees in the same plantation is consistent, which brings great convenience to the accurate estimation of biomass and alleviates the problem of signal saturation to a great extent. When the stand age was added to the independent variable sets, the estimation accuracy for biomass was significantly improved (R² = 0.77, RMSE = 21.16 mg/ha; Figure 4b). However, obtaining a high-precision map of rubber stand age in advance is essential for estimating biomass in large areas. Thanks to the development of remote sensing big data and cloud computing platforms, such as GEE [54], it is now possible to quickly access and analyze decades of remote sensing images and obtain more accurate stand age information for rubber plantations by monitoring changes in land use and rubber tree growth [34,55].

Each rubber plantation also has a similar canopy height, which is another key variable, in addition to DBH, in estimating aboveground biomass on the ground [56,57,58]. Estimating tree canopy heights from traditional optical and radar imagery was also challenging until the GEDI instrument was successfully launched on the International Space Station (ISS) [31]. By producing high-resolution laser-range observations of the Earth’s three-dimensional structure, the GEDI instrument significantly increases the accuracy of the worldwide estimation of tree canopy heights [59]. When canopy height was added in the RF model, the estimation accuracy for rubber plantation biomass also significantly improved (R² = 0.77 and RMSE = 20.91 mg/ha). Compared to the stand age, the accuracy of biomass estimation was slightly higher after adding the canopy height variable (Figure 4b,c), which could explain the differences in canopy height even though the stand age was the same in the same rubber plantation. The height difference between trees in the same plantation can be obtained through the GEDI, while the stand age monitored through time-series images tends to be the same.

Although canopy height was comparable as a biomass estimator variable to DBH, its contribution was found to be significantly lower than that of DBH when tested in this study. Compared with stand age and canopy height, the RF model with only the DBH variable added had the highest estimation accuracy (Figure 4d; R² = 0.97 and RMSE = 7.73 mg/ha). In the other four RF models using multivariate combinations of stand age, canopy height, and DBH, the accuracy of the model with the DBH variable was also significantly higher than that of the model without DBH (Figure 4e–g). This could be primarily explained by the fact that DBH is the most important variable in forest biomass estimation [3,60,61]. The biomass used for verification in this study came from Equation (2), a DBH-based AE model. Moreover, there are still many uncertainties regarding the current canopy height products. The GEDI data used for canopy height retrieval were discrete footprint samples with a diameter of about 25 m, and the other areas were obtained from indirect estimation with data from multiple sources, such as Landsat imagery [62]. In addition, high geolocation uncertainty in GEDI data is likely to introduce high uncertainty in canopy height estimation. The horizontal geolocation error for each 25 m footprint center for the first-release GEDI data is within 10 m, but a recent evaluation of GEDI early-mission data for about 8 months showed a horizontal geolocation error of 23.8 m [63,64]. Furthermore, canopy height data are currently only available for the nominal years of 2019 and 2020 [62,65]. To match the 2016 DBH data, the canopy height used to estimate rubber plantation biomass in 2016 was indirectly estimated by the regression model developed with stand age (Figure 3a). Due to this indirect estimation, there were more dependencies among the DBH, stand age, and canopy height, which is also an important reason why the model accuracy could not be significantly improved when multiple variable combinations were added at the same time.

4.2. Spatial Distribution Pattern for Rubber Plantations and Biomass on Hainan Island

Rubber plantations constitute over a quarter of the forested area on Hainan Island, representing the largest forest plantation ecosystem in the region [33]. Therefore, quantifying the biomass of these plantations holds considerable significance for the island. The spatial distribution of rubber biomass on Hainan Island is intricately linked to the timeline of the establishment and the location of rubber plantations (Figure 1). Following a rapid increase in natural rubber prices after 2003, rubber planting areas on the island expanded rapidly, particularly in the western and northern regions [33,66]. Although the southeastern and southern parts of Hainan offer favorable temperature and rainfall conditions for rubber cultivation, the former is prone to frequent severe typhoons, and the latter, as a genuine tropical area, generates more substantial economic returns from the production of tropical fruits rather than rubber trees [67,68,69]. This phenomenon explains the high average biomass in Sanya city in the south during 2016, which subsequently declined significantly in 2020 (Figure 6). There are few new rubber plantations in these areas, resulting in most of the old rubber plantations gradually being felled to plant other tropical fruits. Consequently, the establishment of new rubber plantations along the southeastern coast and in the south is substantially lower in comparison to the western and northern regions.

With regard to alterations in spatial distribution, the extent of rubber plantations in 2020 was markedly greater than in 2016 (Figure 5). This observation may be ascribed to the limitations of existing algorithms, which solely detect rubber plantations exhibiting forest morphological characteristics via remote sensing approximately five years post-establishment [34]. Consequently, the rubber plantations that principally expanded between 2016 and 2020 were established between 2010 and 2015. Despite a precipitous decline in natural rubber prices following their peak in 2011, the enlargement of rubber planting areas persisted due to a temporal disconnect between plantation growth and price oscillations [66,70,71]. For example, statistics indicate that, between 2010 and 2015, newly established rubber plantations on Hainan Island encompassed 128,000 ha (averaging 21,000 ha annually), while this figure diminished to approximately 3200 ha in 2016 [72]. As a result, by the end of 2020, a substantial portion of Hainan Island’s rubber plantations consisted of trees less than 15 years old, displaying relatively low biomass density. The rise in the proportion of young rubber plantations explains the reduced average biomass for rubber plantations in numerous cities and counties in 2020 compared to 2016 (Figure 6). However, the augmented total area of rubber plantations in 2020 contributed to the elevation of total biomass in the majority of cities and counties. The total biomass of rubber plantations island-wide also increased from 4.52 × 10⁷ mg in 2016 to 5.46 × 10⁷ mg in 2020.

Rubber plantations characterized by higher biomass (exceeding 120 mg/ha) are predominantly situated in central mountainous regions and the border areas of different cities and counties (Figure 5). This distribution may be attributed to the elevated expenses associated with short-term, intensive agricultural development in mountainous locales, rendering the cultivation of perennial rubber trees a more favorable option [73]. The considerable costs involved in updating old rubber plantations in these regions facilitate the preservation of older rubber plantations [34]. For example, the top four rubber plantations with the highest average biomass in 2020 were Qiongzhong county, Wuzhishan city, Chenmai county, and Baisha county, three of which have typical mountainous terrain (Figure 6). In addition, border areas often comprise underdeveloped regions with unfavorable geographic conditions, which, in turn, increases the likelihood of maintaining older rubber plantations [74,75].

4.3. Uncertainty Analysis and Potential Application Prospects

The precision of rubber plantation biomass estimation primarily depends on the accurate determination of three variables: DBH, stand age, and canopy height. Remote sensing offers a comparatively more feasible approach for acquiring stand age and canopy height data than obtaining DBH measurements. The estimation of stand age and canopy height involves a relatively comprehensive algorithm that demonstrates notably high accuracy, with an RMSE of approximately 2 years, when sufficient image data are available [34,55]. Nevertheless, the availability of Landsat 5 images prior to 2000 was limited, leading to a scarcity of data sources and consequently affecting the estimation accuracy. With regard to canopy height estimation, significant improvements in accuracy have been achieved through the use of GEDI data [30,62], but it still has high uncertainty due to intrinsic geometric positioning errors associated with the GEDI and the availability of footprint points only. While the indirect estimation of DBH through RSIs, ECIs, stand age, and canopy height displays high accuracy, the employed algorithm lacks clear theoretical significance. Furthermore, tropical regions are characterized by severe landscape fragmentation [25,76], and many rubber plantations are managed by smallholders with limited land area that is not easy to detect, which exacerbates the problems of collecting stand age and canopy height at the spatial resolution level of satellite imagery. Hence, a more advanced nonlinear regression fitting analyzer [77] that can synthetically consider all disturbance factors related to the rubber tree biomass assessment could be developed to raise the accuracy of the final results.

Despite the existence of limitations, the algorithm employed in this study demonstrates promising potential and can be readily applied to other rubber cultivation regions, such as Xishuangbanna in Yunnan province, China, and various Southeast Asian countries. Firstly, the utilization of platforms such as GEE and the availability of extensive time-series big data, such as from Landsat, Sentinel, and SDGSAT-1, enable accurate estimation of rubber stand age [54,78]. Secondly, the canopy height products from 2020 are on a global scale, and canopy height in recent years can be further mapped using the latest LiDAR data acquired from GEDI and GF-7, as well as L-band SAR data, like PALSAR-2 [13]. By incorporating ground survey data to refine the DBH estimation model and biomass allometric equations specific to different regions, it is anticipated that satisfactory accuracy in biomass estimation can be achieved across diverse geographical areas.

Against the backdrop of persistently low natural rubber prices for over a decade and China’s explicit commitment to carbon peaking and carbon neutrality [79,80], the practical significance of large-scale biomass/carbon storage estimation in rubber plantations is becoming increasingly profound. Accurately determining the carbon sink capacity of rubber plantations holds the potential to facilitate carbon trading initiatives and enhance the economic benefits associated with rubber cultivation. For example, implementation of the Verified Carbon Standard (VCS) in forestry carbon sequestration projects within rubber plantations can yield tangible economic advantages for various stakeholders, including rubber farmers, businesses, and government regulators [81]. Hence, with the continuous upgrading of satellite-borne remote sensing techniques and the latest image processing algorithms, the ideological elaboration of carbon stock surveying from the regional scale to the sample-plot scale and then the individual-tree scale could be achieved in the near future [82].

5. Conclusions

Estimating biomass in rubber plantations holds significant scientific and practical implications, but achieving accurate estimations on a large scale remains challenging. This study employed remote sensing images, ground surveys, and key physiological factors, such as stand age and canopy height, to investigate the impact of different predictors and estimation methods on biomass estimation accuracy in rubber plantations. The results demonstrate that incorporating stand age or canopy height factors into the random forest model based on RSIs and ECIs significantly improves biomass prediction accuracy. This approach effectively addresses the signal saturation problem encountered in traditional remote sensing inversion, increasing the R² value from 0.24 to 0.77. Furthermore, when the DBH—the most widely used variable for biomass estimation using the AE—was added to the model, the prediction accuracy surpassed that obtained when adding stand age or canopy height, resulting in an R² value of 0.97. However, incorporating two or all of the variables (stand age, canopy height, and DBH) into the model did not yield significant improvements in prediction accuracy as compared to the more straightforward approach using RSIs, ECIs, and DBH.

This study applied the optimal model to estimate biomass in rubber plantations on Hainan Island, China’s second-largest rubber planting region, for the years 2016 and 2020. The results reveal that spatiotemporal biomass distribution patterns closely correlate with rubber plantation distribution and establishment years. Since the onset of the 21st century, the increasing price of natural rubber has prompted the establishment of numerous rubber plantations across the island, particularly in the northwest and northern regions. Consequently, these areas exhibit relatively low biomass densities (e.g., <50 mg/ha). In contrast, older rubber plantations in the central mountainous region and the border areas of different cities and counties feature higher biomass densities (e.g., >120 mg/ha). Although the average biomass in over half of the cities and counties in 2020 was marginally lower than in 2016, the total biomass of rubber plantations in most cities and counties exhibited a clear growth trend due to the expansion of plantation areas. By the end of 2020, the total biomass of rubber plantations on Hainan Island reached 5.46 × 10⁷ mg, emphasizing its substantial value for the carbon cycle on Hainan Island.