Estimating Forest Stock Volume in Hunan Province, China, by Integrating In Situ Plot Data, Sentinel-2 Images, and Linear and Machine Learning Regression Models

Abstract

The forest stock volume (FSV) is one of the key indicators in forestry resource assessments on local, regional, and national scales. To date, scaling up in situ plot-scale measurements across landscapes is still a great challenge in the estimation of FSVs. In this study, Sentinel-2 imagery, the Google Earth Engine (GEE) cloud computing platform, three base station joint differential positioning technology (TBSJDPT), and three algorithms were used to build an FSV model for forests located in Hunan Province, southern China. The GEE cloud computing platform was used to extract the imagery variables from the Sentinel-2 imagery pixels. The TBSJDPT was put forward and used to provide high-precision positions of the sample plot data. The random forests (RF), support vector regression (SVR), and multiple linear regression (MLR) algorithms were used to estimate the FSV. For each pixel, 24 variables were extracted from the Sentinel-2 images taken in 2017 and 2018. The RF model performed the best in both the training phase (i.e., R² = 0.91, RMSE = 35.13 m³ ha⁻¹, n = 321) and in the test phase (i.e., R² = 0.58, RMSE = 65.03 m³ ha⁻¹, and n = 138). This model was followed by the SVR model (R² = 0.54, RMSE = 65.60 m³ ha⁻¹, n = 321 in training; R² = 0.54, RMSE = 66.00 m³ ha⁻¹, n = 138 in testing), which was slightly better than the MLR model (R² = 0.38, RMSE = 75.74 m³ ha⁻¹, and n = 321 in training; R² = 0.49, RMSE = 70.22 m³ ha⁻¹, and n = 138 in testing) in both the training phase and test phase. The best predictive band was Red-Edge 1 (B5), which performed well both in the machine learning methods and in the MLR method. The Blue band (B2), Green band (B3), Red band (B4), SWIR2 band (B12), and vegetation indices (TCW, NDVI_B5, and TCB) were used in the machine learning models, and only one vegetation index (MSI) was used in the MLR model. We mapped the FSV distribution in Hunan Province (3.50 × 10⁸ m³) based on the RF model; it reached a total accuracy of 63.87% compared with the official forest report in 2017 (5.48 × 10⁸ m³). The results from this study will help develop and improve satellite-based methods to estimate FSVs on local, regional and national scales.

1. Introduction

The forest stock volume (FSV, m³ ha⁻¹) is the sum of the stem volumes of all the living trees per unit area, and is one of key forest variables for forest resources management and assessments on local, region and country scales [1]. FSVs have a strong relationship with the aboveground biomass (AGB) and carbon stocks [2]. To understand the spatial distribution of carbon in forests and to derive predictions for monitoring carbon stock trends, the FSV must be quantified [3]. Traditionally, the FSV is estimated by sampling several plots, which involves substantial manpower, materials, and financial resources [4]. With the development of remote sensing technology, particularly the Landsat series since 1972, satellite imagery has played an important role in forest inventory. Various studies have been performed to estimate the forest variables using low spatial resolution (LSR, LSR ≥ 30 m), moderate spatial resolution (MSR, 5 m < MSR < 30 m), and high spatial resolution (HSR, HSR ≤ 5 m) [5] data obtained by different optical sensors (e.g., Landsat [6,7,8,9,10,11], MODIS [12,13,14,15,16], SPOT [17,18,19,20], Quickbird [21,22,23,24,25,26], RapidEye [27,28,29,30]), microwave sensors [31,32,33,34], and light detection and ranging (LiDAR) sensors [35,36,37,38,39]. The successful applications of these technological tools have laid the foundation for the estimation of forest variables, such as the FSV, using remote sensing technology.

Two recent developments in the Earth Observation (EO) sector have increased the potential to improve the efficiency of retrieving global forest attributes. The Sentinel-2A and Sentinel-2B satellite series launched by the European Space Agency (ESA) through its Copernicus program in 2015 (S2A) and 2017 (S2B) provide nominal five-day revisit imagery across the globe [40]. The Sentinel-2 imagery includes 13 spectral bands with spatial resolutions ranging from 10–60 m [40]. One of the main purposes of the Sentinel-2 satellite series is vegetation analysis [4]. The Sentinel-2 satellite images provided by the operational environment monitoring system based on the European Copernicus program may be accessed freely. The different spatial resolution bands, the short revisit period, and the rich spectral information have made these images a popular source of remote sensing data for forestry research in recent years. For instance, Persson et al. [41] used Sentinel-2 data to classify common tree species in central Sweden, observing the highest overall accuracy, i.e., approximately 88.2%, using all the imagery bands in the final model. A study conducted by Ho’sciło et al. [42] in southern Poland affirmed that the Sentinel-2 series could accurately delineate tree species (e.g., beech, oak, birch, alder, and larch) with an overall accuracy above 85%. Similarly, Pandit et al. [43] explored the ability of Sentinel-2 images to estimate the forest biomass in Nepal; their biomass estimation model achieved an R² = 0.81 and RMSE = 25.57 t ha⁻¹. Zarco-Tejada et al. [44] demonstrated the potential of Sentinel-2A data to estimate the chlorophyll content in open canopy conifer forests (R² > 0.7 for June; R² > 0.4 for December). In the Brazilian Amazon, Lima et al. [45] performed comparative research on monitoring selective logging using Sentinel-2 and Landsat-8 OLI imagery. These authors found that Sentinel-2 data (43.2% detected) were more effective in detecting logging concessions than Landsat 8 data (35.5% detected). In Poland, Grabska et al. [46] used the Sentinel-2 time series to map forest stand species in the Carpathian Mountains and reported higher accuracy, i.e., a 5–10% improvement in overall accuracy compared with only using single date imagery. These studies demonstrate the potential of Sentinel-2 for forest vegetation monitoring.

A typical traditional remote sensing image processing approach includes data downloading and local computer processing, which render image processing computationally demanding, and also hinder the processing capability for large datasets. Regarding imagery processing, the development of cloud computing technology has been fundamentally changing traditional remote sensing image processing. The Google Earth Engine (GEE) portal is a powerful cloud-based computing platform for image processing [47,48]. The GEE archives massive, publicly-available remote sensing data, provides a programming environment, programming tools, and virtual machines for users with relatively simple code, and can process imagery data online. The GEE greatly improves the processing efficiency when using substantial amounts of remote sensing data. In recent years, the GEE was used in land cover mapping [49,50,51,52,53,54,55,56,57,58], agricultural applications [59,60,61,62,63], disaster management, and earth sciences studies [64,65,66]. This remote sensing data processing cloud platform makes the rapid processing of Sentinel-2 images covering large areas possible.

In the FSV estimation field, some studies have examined assessments of the FSVs using remote sensing. For instance, Condés et al. [67] found the model prediction of the plot-level growing stock volume using satellite images and field data to be useful; the result showed that the adj-R² increased from 0.19 to 0.42. Using the random forests (RF) regression algorithm, Chrysafis et al. [68] estimated the FSV based on Sentinel-2 image, which provided relatively better results (R² = 0.63, RMSE = 63.11 m³ ha⁻¹) than Landsat-8 OLI images (R² = 0.62, RMSE = 64.40 m³ ha⁻¹). However, some studies were conducted that combined optical images and microwave data to estimate forest variables [69,70,71]. A noteworthy study was conducted by Mauya et al. [72]; these authors assessed the multiple linear regression (MLR) models built by Sentinel-1, Sentinel-2, and ALOS PALSAR-2 images to predict the FSV, and found that Sentinel-2 images performed best with an RMSEr = 42.03% and a pseudo-R² = 0.63. For predicting forest variables, Pham et al. showed that machine learning algorithms were likely to become more attractive in remote sensing [73]. These authors suggest that future studies using more methods, large areas, and Sentinel-2 data to predict the FSV should be conducted.

However, ground plot survey data are still indispensable for remote sensing modeling [74]. The costs of ground plot surveys have always been high, which has presented some obstacles to the estimation of the provincial FSV by remote sensing. In addition, the traditional sample location survey technology often produces serious positional deviations, which may impact the modeling accuracy of the plot data-based remote sensing estimations and predictions [75]. Thus, traditional sample location technology is also an important reason for the inaccurate matching of sample plots and pixels, which results in estimation bias. Furthermore, to the best of our knowledge, little to no research has been conducted to compare the results of the RF, support vector regression (SVR), and MLR using Sentinel-2 images to predict FSVs. In addition, no FSV mapping has been conducted in Hunan Province, and this research gap directly affects forest policy making and management. Moreover, Hunan Province is located in southern China, where cloudy conditions frequently occur; these conditions are a great challenge to mapping the FSV in Hunan Province using remote sensing images.

Hence, Sentinel-2 data on the GEE platform, 459 sample plots, and three algorithms were used in this study to achieve the following objectives: (1) to identify and select the most important variables of the Sentinel-2 images for FSV estimation using plot-level tree measurements from a large number of in situ sites in Hunan Province, southern China, where forest species and stand structures are complex; (2) to assess and understand the performance of machine learning algorithms and the typical MLR models for FSV estimation; and (3) to map the FSV in Hunan Province. This study will help move towards the overall goal of developing and improving GEE-based remote sensing approaches to estimate the FSV on local, regional and national scales.

2. Materials and Methods

2.1. Study Area

The study was conducted in Hunan Province, covering an estimated area of 211,854.69 km². This area is a transition zone between the middle and lower ridges of the Yangtze and Xiangjiang river basins and the Yungui Plateau. The area has a subtropical monsoon climate with four distinct seasons. The annual precipitation ranges from 1200 to 1700 mm, and the mean annual temperatures vary between 16 °C and 18 °C [76]. These favorable climatic conditions make Hunan Province one of the major forest areas in China. In the study area, the Global PALSAR-2/PALSAR Forest/Non-Forest Map product was used as ancillary data to determine the forest areas and show the distribution of the training plots and test plots (Figure 1).

2.2. In situ Sample Plot Data Collection

Based on the historic forest survey data, 10 forest types were chosen for this study. The forest types were first analyzed using the 8th National Forest Inventory (NFI) data, and then the top seven tree species were selected according to the stock volume. These species included Cunninghamia lanceolata, Pinus massoniana, Quercus sp., Pinus elliottii, Populus sp., Cinnamomum camphora, and Cupressus funebris. In addition, three mixed forest types were selected, including broad-leaved mixed forests, coniferous and broad-leaved mixed forests, and coniferous mixed forests. These ten types account for over 98% of the total FSV in Hunan (

Table 1. The proportion of the FSV from the ten forest types in Hunan Province, China.

Sample Plot Type	Stock Volume (m³)	Proportion (%)
Total forest stock volume	330,992,700	100.00
Cunninghamia lanceolata	109,035,700	32.94
Pinus massoniana	46,395,900	14.02
Quercus sp.	5,098,800	1.54
Pinus elliottii	3,020,400	0.91
Populus sp.	2,022,900	0.61
Cinnamomum camphora	2,244,100	0.68
Cupressus funebris	1,329,900	0.40
Broad-leaved mixed forests	84,515,000	25.53
Coniferous and broad-leaved mixed forests	36,681,900	11.08
Coniferous mixed forests	34,136,400	10.31
Total	324,481,000	98.00

Note: The Latin name of the tree species is in italic.

).

A total of 15 survey levels were generated in each forest type. Each forest survey type was further divided into five tree height levels and three canopy density degrees (0.2–0.39, 0.4–0.69, and 0.7–1.0). At each surveyed level, four sample plots were designed and investigated. Finally, 459 useful circle sample plots with a diameter of 30 m were obtained from 15 October 2017, to 30 April 2018. The tree heights (measured by a VERTEX LASER VL5 manufactured in Haglof, Sweden) and diameters at breast height (DBHs) were measured for all individual trees (measured by a caliper with a DBH greater than 5 cm) in the circle sample plot. TBSJDPT (Figure 2) was put forward and used to locate the sample plot center and the individual trees [39].

To match the pixels of the Sentinel-2 imagery, the sample plot data were resampled using the ArcGIS 10.2.2 software (Esri, Redlands, CA, USA). By using the center position of each sample plot, the 25 m × 25 m area sample plots were generated by using the ArcTools “Buffer and Feature Envelope To Polygon”. The generated shapefile was then used to extract the trees in the sample plot. Using the formula built by Hunan Forestry Industry Survey and Design Research Institute, we calculated the FSV of each tree, summed it to obtain the FSV of the sample plots, and then scaled it up to obtain the FSV per hectare (

Table 2. The FSV calculation formula of the seven major tree species in Hunan Province [78].

Tree Specie	Formula
Cunninghamia lanceolata	V = 0.000058777042D1.9699831H0.89646157
Pinus massoniana	V = 0.000062341803D1.8551497H0.95682492
Quercus sp.	V = 0.000050479055D1.9085054H0.99076507
Pinus elliottii	V = 0.000086791543D(1.6638000575+0.0094299757(D+10H))H(0.9693404868-0.0292030826(D+2.5H))
Populus sp.	V = 0.000041028005D1.8006303H1.13059897
Cinnamomum camphora	V = 0.000050479055D1.9085054H0.99076507
Cupressus funebris	V = 0.000058777042D1.9699831H0.89646157

Note: The FSV calculation formula of the other tree species are supplied in supplementary 1.

). In addition, we used the “boxcox” function in the R software to transform the FSV data to a normal distribution (f(y) = y^λ) [77].

2.3. Sentinel-2 Images Preprocessing and Variable Calculation

To obtain a temporal match between the investigated sample plot data and the satellite data for the whole of Hunan Province, the GEE filtering options were used to select the Sentinel-2 Level-1C (TOA, top-of-atmosphere) satellite images between 1 May 2017, and 31 October 2017, and between 1 May 2018, and 31 October 2018. The Level-1C product is an ortho-image in UTM/WGS84 projection. Level-1C processing includes radiometric and geometric corrections, including orthorectification-Digital Elevation Model (DEM, processed from Shuttle Radar Topography Mission 90 m Digital Elevation Data Version 4) to project the images in cartographic geometry. Per-pixel radiometric measurements are provided in TOA reflectances along with the parameters to transform them into radiances, and spatial registration was processed on a global reference system with subpixel accuracy. In addition, the Level-1C products provide a bitmask band (QA60) to show cloud information. More details about the Level-1C products could be found in the Sentinel-2 User Handbook [79]. The Level-1C product stored in the GEE has been shown to be suitable for direct use for land mapping by Hird et al. [57]. The cloud mask steps were processed in this study. The QA60 band was first used to identify and mask out the flagged cloud and cirrus pixels. Then, the two thresholds of Band 1 ≥ 0.15 and Band 2 ≤ 0.25 were used to mask the remaining clouds and noise [57,80]. A total of 3150 images covering the entire study area over the study period were processed. After that, the Level-1C TOA reflectance data were resampled with a resolution of 25 m to match the spatial resolution of survey plots. Then, we extracted the median values of the characteristic variables using the GEE cloud computing platform [58]. The characteristic variables included the multispectral bands (B2, B3, B4, B5, B6, B7, B8, B8A, B10, B11, and B12) and the vegetation indices selected according to the research of Astola et al. [4], Wittke et al. [81], and Xia et al. [82]. The vegetation indices were calculated as follows (

Table 3. Vegetation indices calculated from Sentinel-2 data.

Characteristic Variable	Index Short Name	Calculation Method
Vegetation indices	NDVI_B5	(B5 − B4)/(B5 + B4)
	NDVI_B6	(B6 − B4)/(B6 + B4)
	NDVI_B7	(B7 − B4)/(B7 + B4)
	NDVI_B8	(B8 − B4)/(B8 + B4)
	NDVI_B8A	(B8A − B4)/(B8A + B4)
	SAVI	1.5*(B8 − B4)/(B8 + B4 + 0.5)
	RVI	B8/B4
	MSI	B8/B11
	EVI	2.5*(B8 − B4)/(B8 + 6*B4 − 7.5*B2 + 1)
	EVI2	2.5*(B8 − B4)/(B8 + 2.4*B4 + 1)
	TCW	0.1509*B2 + 0.1973*B3 + 0.3279*B4 + 0.3406*B8 + 0.7112*B11 + 0.4572*B12
	TCB	0.3037*B2 + 0.2793*B3 + 0.4734*B4 + 0.5585*B8 + 0.5082*B11 + 0.1863*B12
	TCG	− 0.2848*B2 − 0.2435*B3 − 0.5436*B4 + 0.7243*B8 + 0.0840*B11 − 0.1800*B12

).

2.4. Selection of Relevant Variables for FSV Estimation

For the MLR method, the “redun” function was used to avoid collinearity. The parameter R² was set to 0.9. Then, the variables were further filtered based on their significance (i.e., P < 0.05) and variance inflation factor (i.e., VIF > 10). For the machine learning methods, in order to obtain the best fit model, the selection of the relevant variables for the FSV estimation was carried out using the “variable selection using random forests” (VSURF) package [83]. All the extracted data were analyzed in order to obtain the variable importance measures (VIM) based on the VSURF package in the R 3.5.3 software (R Core Team, Vienna, Austria). Importance type 1 represents the percentage increase in the mean square error (PercentIncMSE), and importance type 2 represents the increase in the NodePurity (IncNodePurity); each had the ability to evaluate the importance of variables [83]. The selected variables were then used in RF modeling and SVR modeling.

2.5. Statistical Models for Estimating the FSV

The MLR is a statistical technique that uses several explanatory variables (independent variables) to predict the outcome of a response variable (dependent variable). The MLR is a widely used method to construct regression models for various applications [72,84,85,86,87]. The goal of the MLR is to model the linear relationship between the explanatory (independent) variables and response (dependent) variable (FSV). The formula for the MLR is

$$\mathrm{y}=a+b_1{x}_1+b_2{x}_2+b_3{x}_3+\dots +b_i{x}_i+\epsilon$$

(1)

where i is the number of variables, $\mathrm{y}$ is the dependent variable, x_i is the independent variable, a is the intercept, b_i is the slope coefficient for the explanatory variable, and ε is the error term of the model.

The SVR maps the output data to the high-dimensional feature space by defining the kernel function and constructs an optimal classification hyperplane in this space. The SVR can be regarded as a linear algorithm in high-dimensional space [88]. The SVR has been proven to be a popular method for regression modeling [89,90,91,92]. In this study, it was constructed using the following SVR function from the training process:

$$\mathrm{FSV}={{\displaystyle \sum }}_{i=1}^n{\alpha }_{i}k\left(x_i;x\right)+b$$

(2)

where $k\left(x_i;x\right)$ is the kernel function, ${\alpha }_i$ is the Lagrange multiplier, $x_i$ is the training vector, and $b$ is the bias term in the regression.

Using the e1071 package in the R 3.5.3 software, the parameters of the SVR were optimized and cross-validated. The radial basis function (RBF) kernel was selected following previous studies [71,93]. Then, the best regression model was optimized by the tune function (the SVM-type is eps-regression, the SVM-Kernel is the radial, cost is set as 2^c (c is 2-9 sequence with an interval of 1), epsilon is set as a 0-1 sequence with an interval 0.01, and gamma is default).

The RF is a decision tree algorithm and an effective machine learning model for predicting a forest of variables. Based on its powerful modeling capabilities, the RF regression has been widely used in scientific research [94,95,96,97,98,99]. The principle of the RF algorithm is to use the bootstrap method to randomly extract multiple samples to generate a group of regression trees (ntree) from the original sample population. To build each tree, a randomly-chosen subset of predictors is used at each splitting node (mtry) by the RF, and there is no need to prune each tree grown. For each tree grown, a procedure named the “out-of-bag” (OOB) error was used to independently build each tree based on the training data. The variable importance (VI) and OOB error were calculated by the RF algorithm [100].

The OOB error can be estimated as follows:

$$OO{B}_{error}= \frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2$$

(3)

where $y_i$ is the measured FSV, ${\widehat{y}}_i$ is the predicted FSV, and n is the total number of OOB samples.

After selecting the variables, the number of grown trees (ntree) and the number of variables (mtry) at the time of node splitting were determined using the minimum mean square error as the criterion. First, the RF model’s ntree parameter was set to 500, and the mtry parameter was set to the total number of variables. Then, the mean square error was used to determine the best number of variables (mtry) and the best number of trees (ntree). After calculating the best parameters, the RF regression model was established and tested. In our three models, 70% of the total samples were used for training samples, and 30% of the overall samples were used for testing.

Two statistical indicators were selected to evaluate the performance of the MLR, SVR, and RF models [88,101]. The coefficient of determination (R², Equation (4)) is the proportion of the variance in the dependent variable that can be explained by the independent variables. The root mean square error (RMSE, Equation (5)) is the standard deviation of the residuals, which is the difference between the surveyed data and the fitted model.

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(4)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{n}}$$

(5)

where $y_i$ is the measured FSV, $\overline{y}$ is the mean measured FSV, ${\widehat{y}}_i$ is the predicted FSV, i is the same index, and n is the sample plot number.

3. Results

3.1. Characteristics of the in Situ FSV Data

For this study, 321 out of the 459 total sample plots (70%) were used as the training data, and 138 out of the 459 total sample plots (30%) were used as the test data.

Table 4. Summary statistics for the FSVs of the training data and test data.

Descriptive Statistics	Training Data	Transformed Training Data	Test Data	Transformed Test Data
Mean	121.11 (m3 ha−1)	3.98	120.53 (m3 ha−1)	3.95
Median	103.33 (m3 ha−1)	4.08	98.37 (m3 ha−1)	4.02
Minimum value	1.42 (m3 ha−1)	1.11	4.25 (m3 ha−1)	1.55
Maximum value	577.49 (m3 ha−1)	6.87	450.11 (m3 ha−1)	6.37
Variance	9019.13	1.15	9053.02	1.24
Kurtosis	2.73	−0.23	0.22	−0.89
Skewness	1.40	−0.18	0.89	−0.11
Number of sample plots	321	321	138	138

lists the summarized statistics of the plot-level training data and test data. For the training data, the minimum FSV per hectare was only 1.42 m³, and the maximum FSV per hectare was 577.49 m³. For the test data, the FSV per hectare ranged from 4.25 m³ to 450.11 m³. These large data ranges were the reason for the large variances in the training and test data. For the training data and test data before transformation, it is clear that the training data and test data did not follow a normal distribution (Figure 3, the green histograms with red curves). A λ = 0.3030303 was calculated and used to transform the training data and test data, and it is clear that the data were transformed into an approximately normal distribution (Figure 3, the red histograms with green curves).

3.2. Major Variables Related to the FSV Data

Figure 4 shows the results of the numbers of selected variables using the RF algorithms provided by the VSURF package. First, no variables were eliminated based on the VI mean and VI standard deviations. Based on the OOB error, eight variables were selected: B2, B3, B4, B5, B12, TCW, TCB, and NDVI_B5. The PercentIncMSE (percent increase in the mean square error) and IncNodePurity (increase in the NodePurity) estimated from the RF OOB data were used to rank all predictor variables according to their abilities to estimate the FSV; the higher the value, the more important the variable (Figure 5).

For the MLR algorithm, 14 out of the 24 variables were eliminated by the “redun” function to avoid multicollinearity; these eliminated variables were NDVI_B8, EVI2, SAVI, NDVI_B7, B7, B8A, TCB, NDVI_B6, B6, B3, B11, TCG, NDVI_B8A, and B4. The remaining variables were reselectedusing the P significance in linear modeling. Finally, B5 and MSI were selected for the MLR model.

3.3. Optimal Regression Model for the RF, SVR, and MLR

For the RF model, before selecting the best RF regression model, the mtry and ntree parameters must first be determined. Based on the number of selected variables, an error rate loop algorithm was used to calculate the smallest error rate. In Figure 6a, it is clear that mtry = 5 resulted in the smallest error rate. In Figure 6b, ntree = 257 was determined to be the best parameter.

In the SVR algorithm, 909 models were trained. The best model that was selected had cost = 8, gamma = 0.125, and epsilon = 0.68. Then, the best SVR regression was used to build the regression model for the FSV prediction. In the MLR algorithm, the results of the MLR model are shown in

Table 5. Coefficient of the MLR model.

	Estimate	Std. Error	t Value	P
(Intercept)	5.6320995	0.9240733	6.095	3.17e − 09 ***
B5	−0.005478	0.0004089	−13.397	< 2e − 16 ***
MSI	1.9034603	0.4423448	4.303	2.24e − 05 ***

*** P values at significance level of 0.001.

.

A basic hypothesis test was made for the linear model (Figure 7). In Figure 7, the top left graph, the residuals versus the fitted values, was used to test the linearity assumptions. The scatter points were concentrated near a straight line, which indicated that the linear relationship was good. The top right graph was a normal Q-Q plot, which was used to test the normality. The scatter points were mainly concentrated along the straight line, which indicated that the residual normality was good. The bottom left graph was the scale-location graph, which was used to test the homoscedasticity. The points were randomly distributed around the curve, which indicated that the homoscedasticity assumption was accepted. The bottom right graph used the residual and leverage to determine the outliers, and no high leverage points or strong influence points were found.

3.4. Comparison of the Predicted FSV Estimates among the Three Models (MLR, SVR, and RF)

In the training phase, the MLR model had the worst performance with the smallest R² = 0.38 and highest RMSE = 75.74 m³ ha⁻¹ (Figure 8a). This model was followed by the SVR model with an R² = 0.54 and RMSE = 65.60 m³ ha⁻¹ (Figure 8c). The best performance among the three algorithms was that of the RF model, which had the highest R² = 0.91 and smallest RMSE = 35.13 m³ ha⁻¹ (Figure 8e). In the test phase, it was clear that the RF model also performed the best among the three kinds of models (R² = 0.58, RMSE = 65.03 m³ ha⁻¹) (Figure 8f); the next best was the SVR model, with an R² = 0.54 and RMSE = 66.00 m³ ha⁻¹ (Figure 8d); the MLR model performed the worst, with the smallest R² = 0.49 and highest RMSE = 70.22 m³ ha⁻¹ (Figure 8b).

3.5. Modeling Results Comparison between Selected Variables and all Variables

Table 6. The performance comparison of machine learning methods under different variables.

	Methods	Best Model Parameters	R2.training	RMSE.training (m3 ha−1)	R2.test	RMSE.test (m3 ha−1)
Selected variables	RF	mtry = 5	0.91	35.13	0.58	65.03
		ntree = 257
	SVR	cost = 8	0.54	65.60	0.54	66.00
		gamma = 0.125
		epsilon = 0.68
All variables	RF	mtry = 7	0.92	34.83	0.58	66.04
		ntree = 495
	SVR	cost = 4	0.61	60.58	0.51	67.86
		gamma = 0.04166667
		epsilon = 0.57

shows the modeling performance of all variables and selected variables under the machine learning methods. In the training phase, the result of the RF model with all variables (e.g., R².training = 0.92, RMSE = 34.83 m³ ha⁻¹) is basically the same as that of the RF model with the selected variables (e.g., R².training = 0.91, RMSE = 35.13 m³ ha⁻¹). And the results of the RF model are better than those of the SVR model with all variables (e.g., R².training = 0.61, RMSE = 60.58 m³ ha⁻¹) and selected variables (e.g., R².training = 0.54, RMSE = 65.60 m³ ha⁻¹). In the test phase, the result of the RF model with selected variables (e.g., R².training = 0.58, RMSE = 65.03 m³ ha⁻¹) is basically consistent with the result of the RF model with all variables (e.g., R².training = 0.58, RMSE = 66.04 m³ ha⁻¹). Similarly, the results of the RF model in the test phase are better than those of the SVR models with all variables (e.g., R².training = 0.51, RMSE = 67.86 m³ ha⁻¹) and selected variables (e.g., R².training = 0.54, RMSE = 66.00 m³ ha⁻¹). In the SVR model, the performance of all variables in the training phase is slightly better than that of the selected variables; and in the test phase, the performance of the selected variables is slightly better than that of the all variables.

3.6. Map of the FSV Estimation in Hunan Province in 2017

First, the Global PALSAR-2/PALSAR Forest/Non-Forest Map product was used to extract the forest area in Hunan Province. A total of 143,258,710 pixels with a 25 m spatial resolution were defined as forest areas in this product (Figure 9. left), totaling 8.95 × 10⁶ ha of forest area. Then, the selected variables and the best RF model were used to estimate the FSV in the study area. The right of Figure 9 shows that the FSV in Hunan Province ranged from 12.7421 m³ ha⁻¹ to 269.649 m³ ha⁻¹. The mean FSV per hectare was 39.09 m³, and the total FSV in Hunan Province was 3.50 × 10⁸ m³.

4. Discussion

The main purpose of this study was to evaluate the potential variables of the Sentinel-2 data using different algorithms to predict the FSV based on reliable field survey data, and to map the FSV for the first time in the southern province of China. The FSV is an important variable of forest management reports at the provincial and national levels. The use of free remote sensing imagery (e.g., Sentinel-2) and cloud processing platforms (e.g., GEE) to process and build prediction models for estimating and mapping the FSV is especially important in southern China. One of the reasons for this is that the area is covered by clouds many days each year, which seriously affects provincial forest mapping research that uses remote sensing images. Another reason is that the provinces of southern China are an important part of China’s forestry. Through this research, we can effectively use the Sentinel-2 data with the GEE platform and apply the RF algorithm to spatially map the FSV in Hunan Province. In addition, this approach is also an effective way to conduct forest carbon monitoring, which is sensitive and important to climate change.

Based on the spectral bands and vegetation indices extracted from Sentinel-2 data, this study has shown that B5 (Red-Edge 1) was the most important variable when estimating the FSV using both the machine learning methods and MLR method, which had been confirmed in recent studies concerning forest prediction [102] and tree species classification [103]. In the gross primary productivity field (GPP), Lin et al. found that the red-edge band was useful for estimating the GPP, and noted that the red-edge reflectance was sensitive to the leaf chlorophyll content [104]. In addition, the leaf chlorophyll content was an important forest variable. Except for B5, the modeling variables selected by the machine learning models and the MLR were not the same. However, the accuracies of the different model verification results were not very different, which shows that B5 has a substantial advantage in estimating the FSV. In the study conducted by Chrysafis et al. [68], when using the RF algorithm and Sentinel-2 images to estimate the FSV in a Mediterranean forest ecosystem, they found that the most important variable was B11 (SWIR 1), which was different from our findings in this study. Our research was consistent with a study conducted by Astola et al. [4], which showed that B5 was the most important variable with which to predict the FSV. Regarding the FSV prediction performance, we compared our R² that measured our predictive capability with that of Astola et al. [4], who conducted research using Sentinel-2 and multilayer perceptron and regression trees to estimate the FSV; we found that our results (R² = 0.58) were slightly better than their best results with a multilayer perceptron model (R² = 0.56). This observation may occur because the RF algorithm usually performs better than the multilayer perceptron [71], or the sample survey data based on our advanced positioning technology improved the estimation accuracy. Regarding the variable selection in machine methods, our results also showed that the traditional vegetation indices did not perform well when estimating the FSV in this study, and Lu et al. [105] found the same trend, i.e., that the original band performed better than vegetation indices when estimating the FSV. Table 6 shows that the training and test results before and after selecting the variables of the RF modeling were basically the same, which indicated that the VSURF package was a good tool with which to select variables for RF modeling. However, Table 6 also shows that the selected variables had a certain impact on SVR modeling. This is because the VSURF package is a variable selecting tool based on the RF model, and is not fully applicable to the SVR model.

Figure 8 shows that the three algorithms have some saturation problems in the training and test phases. In the training phase, the maximum value of the RF estimation (363.63 m³ ha⁻¹) was larger than the other two algorithms (SVR, 336.79 m³ ha⁻¹; MLR, 308.39 m³ ha⁻¹), and the minimum value of the RF estimation was smaller (4.97 m³ ha⁻¹) than the other two algorithms (SVR, 6.52 m³ ha⁻¹; MLR, 5.52 m³ ha⁻¹). In the test phase, although the maximum value estimated by MLR (290.26 m³ ha⁻¹) was the largest, the minimum value (7.48 m³ ha⁻¹) was the smallest, and the data range of the RF estimated value was the largest (7.48–272.23 m³ ha⁻¹). This indicated that the RF model was the best model with which to estimate the FSV. In addition, the RF exhibited the best performance among the three algorithms according to the R² and RMSE. When using the RF model to predict the FSV in Hunan Province (Figure 9), we found that the smallest FSV per hectare was 12.7421 m³, which was larger than the smallest FSV measured in the field (1.42 m³ ha⁻¹), and the largest FSV per hectare was 269.649 m³, which was smaller than the highest FSV measured (577.50 m³ ha⁻¹). This result indicated that the RF model overestimated the low FSV values and underestimated the high FSV values, which may be due to the common saturation problem in optical remote sensing vegetation analyses [83]. The overestimation problem could be caused by the understory vegetation (e.g., shrub and grass), which typically impacts the reflectance values. The high FSV areas often have complex canopy structures, which may affect the reflectance values. A similar study by Ou et al. [106] also found this was a common problem in the estimation of the FSV or biomass using multispectral remote sensing data. However, if the forest area with a low FSV and the forest area with a high FSV reach a certain ratio, the underestimation and overestimation problems of the RF will reach a certain balance. For this reason, when using the RF to estimate the FSV over a large area, it may offset some errors caused by the defects of the model or the image data. Regarding the SVR model, Figure 8c shows two distinct trend lines consisting of points. This trend was caused by the large data volume (n = 321 in the training phase) and parameter optimization. Due to the parameter optimization, as many samples as possible were within this hyperplane; therefore, many data points were concentrated on the edge.

After mapping the FSV for the whole of Hunan Province, a statistic using the ENVI 5.3 “Quick Stats ” tool (Exelis Visual Information Solutions, Boulder, Colorado) was calculated, and the mean FSV was 39.09 m³ ha⁻¹. For comparison, we calculated the mean FSV of Hunan Province based on the Analysis Report of Hunan Forestry Statistics Annual Report in 2017, and found that the value was 42.15 m³ ha⁻¹ [107]. The estimation accuracy of the mean FSV reached 92.74%. However, for the sample plot data, the mean values were 121.11 and 120.53 m³ ha⁻¹ for the training and test, respectively. This result shows that among the sample plots we selected, there were too few samples with a small FSV. Using these data to directly model may cause some errors in the model prediction results. However, before modeling, a normal transformation was performed, and the mean values were 3.97 (Table 4) and 3.04 (39.09^λ) before and after transformation, respectively (Figure 3). The normal transformation narrowed this gap, which may compensate for the impact of uneven data sampling to some extent [108]. Finally, the total FSV that we predicted in Hunan Province was 3.50 × 10⁸ m³ in forest areas. At the end of 2017, the Hunan provincial government published a report on the FSV stating that the total value was 5.48 × 10⁸ m³ [107]. Based on this observation, the RF model reached an accuracy of 63.87% in predicting the FSV. We also focused on the forest area reported on the government report in 2017; this showed that the forest area was 1.3 × 10⁷ ha [107], which was different from that extracted from the PALSAR-2/PALSAR Forest/Non-Forest product in 2017 (8.95 × 10⁶ ha). This difference is also an important factor that affects the estimation result of the FSV. In addition, we paid attention to the study conducted by Shen et al. [83], who predicted the forest biomass in Guangdong Province using the RF model and remote sensing data, which achieved slightly lower accuracy (58.88%) than that observed in this study. However, one important difference is that Shen et al. estimated the biomass, and we estimated the FSV.

Since Hunan Province is located in southern China, the greatest challenge was to obtain images with little or no cloud data that covered the whole province. To reduce the effects of clouds and noise on the images, we used three masking methods. To reduce or even eliminate the holes caused by the masks on the images, we processed a total of 3150 images and extracted the median values of the overlapping image parts. As all the sample plots were in a relatively homogeneous forest environment (at least within 30 m of the sample plot boundary is homogeneous), the sample plots vector file was used to extract the mean values of the spanned pixels. Through this method, high-quality remote sensing data for southern China can be effectively obtained, and the potential of using Sentinel-2 to map the FSV in southern China is increased. Such a potential is truly important in southern China, where forest plantations grow rapidly in order to meet the demand for wood products, and for ecological protection [109]. The dates of the selected remote sensing image are mainly based on two considerations. First, the trees in the forest grow slowly in one year. The difference between the acquisition time of the sample plot data and the remote sensing data is less than one year, which has little effect on the results. In addition, we also use the median value of the two-year growing season to reduce this effect. Furthermore, we consider that the spectral information between May and October is better, but the image in one year could not cover the entirety of Hunan Province, so we chose to use the data for two year periods.

Although our prediction results were good, some limitations of the study should be noted. The first are the limitations of the Sentinel-2 Level-1C TOA data. While three mask steps were used to process the Sentinel-2 Level-1C TOA data, this approach could not eliminate the bidirectional reflectance effects from changing the sun–sensor–surface geometries between the acquisitions. Further processing steps that transform the TOA data to bottom-of-atmosphere (BOA) surface data would improve the results; however, the GEE platform did not offer an online tool for this purpose. The local tools that can convert large-scale Sentinel-2 TOA data to BOA data are limited. According to Hird et al. [57], although the TOA data were used, the annual composites of Sentinel-2 input variables would suffice to minimize the bidirectional effects to an acceptable degree. In this study, since approximately two years of Sentinel-2 data were used to extract the band characteristics and vegetation indices, it would be acceptable for this purpose. The second limitation is the modeling algorithm. Although the RF algorithm has been proven to exhibit outstanding performance, the phenomena of the underestimation of high values and the overestimation of low values always exists. The innovation of the algorithms to solve this problem will be a future research direction. The third limitation is the value distribution of the sample data. Although we conducted data transformation, this approach still cannot avoid the modeling bias caused by the uneven distribution of the data values of the sample. The fourth limitation is that not all the data analysis work can be implemented on the GEE platform. Due to the lack of some regression functions, we had to use the R software for auxiliary analyses, which seriously reduced the efficiency of the data analysis. As the Sentinel-2 BOA data are being continuously processed and uploaded to the GEE platform, we believe that the global data products will soon be available on the GEE platform. At the same time, we also noticed the good performance of deep learning algorithms in forest variable estimations [110]. Based on these limitations, we suggest using Sentinel-2 BOA data, deep learning algorithms, and more reasonable sample plot data to conduct FSV estimation research in the future.

5. Conclusions

In this study, based on the Sentinel-2 derived variables, three modeling algorithms (RF, SVR, and MLR) were compared to explore the possibilities of using the Sentinel-2 imagery variables for FSV estimation; the RF model was found to be the best method. The red edge B5 was found to be the most important variable in both the machine learning methods and the linear regression method for FSV estimations. Finally, the best RF model was used to predict the FSV and assess the FSV prediction accuracy. The FSV mapping filled the gaps in FSV study in Hunan Province, China, and provided a reference for the distribution of the FSV in the region. Furthermore, the study also indirectly proved the potential to use Sentinel-2 data combined with the RF model to estimate FSVs, laying the foundation for plotting the FSV in southern China. Future research will estimate FSVs using Sentinel-2 BOA data and deep learning algorithms.