Uncertainty Analysis of Forest Aboveground Carbon Stock Estimation Combining Sentinel-1 and Sentinel-2 Images

1. Introduction

Forest ecosystem plays a vital role in understanding climate change and maintaining carbon balance [1,2,3]. Parameter such as aboveground carbon stock (AGC) measures energy cycling between the organic and inorganic environment. It is also an important indicator for investigating climate change [4,5]. Therefore, accurate estimation of forest AGC is crucial for humans to assess and protect the forest ecosystem [6,7]. However, there are certain uncertainties in estimating AGC at both large and small scales [8]. Ignoring these uncertainties will affect the accurate estimation of forest AGC. Hence, reducing the uncertainty to improve the accuracy in the estimation process has become a significant challenge [9].

The methods for estimating AGC mainly contain the traditional plot survey and remote sensing estimation. Although the plot survey can obtain accurate data, it requires a considerable investment of human and financial resources. Whereas remote sensing-based methods have the characteristics of being fast, low-cost, large-scale, and less destructive. They had been widely used in the estimation of forest aboveground biomass (AGB) and AGC [10,11,12]. Optical remote sensing data have a great advantage in obtaining structural parameters of the forest; their rich spectral information have been widely applied in studies of AGC estimation [13]. However, it is susceptible to the strong influence of canopy thickness, so it is difficult to obtain the vertical structural parameters of the forest and the information it obtains tends to be saturated in areas of high biomass and carbon stock, which has a disadvantage in estimating high AGC [14,15]. Many researchers have found that using only single optical data to estimate forest AGB or AGC has high uncertainty [16,17,18]. Synthetic Aperture Radar (SAR) data has the advantages of longer electromagnetic wavelength, strong penetration ability, and sensitivity to structural attributes [19], so it can overcome the shortcomings of optical data and can obtain valuable information concerning the vertical vegetation canopy. Therefore, combining optical with SAR data to estimate forest AGB and AGC is considered a feasible method [20,21]. For instance, Gesta et al. [22] combined optical and SAR data to estimate AGB and AGC in falcata and to assess the potential for synergy between optical and SAR data in estimating AGB and AGC. Their results showed that the AGB and AGC models generated by combining the two types of data are more robust and accurate than models derived from optical or SAR data. Given the above findings, we combine optical and SAR data to estimate forest AGC in this study. Our goal is to explore the potential of combining different types of data to reduce the uncertainty generated in the estimation process.

Uncertainty is an umbrella term for concepts that include inaccuracy, ambiguity, and inconsistency [23]. It is very common in the process of carbon stock estimation. There are three primary sources of uncertainty in the whole process of AGC estimation: measurement uncertainty, sampling uncertainty, and model uncertainty [24,25]. Shettles et al. [26] regarded model uncertainty as the main source of error, accounting for about 70% of the total uncertainty. Sources of model uncertainty mainly include the uncertainty of input variables [27], improperly setting the model function form [28], residual variability of the model [29], and parameter error of the model [30]. Uncertainty in the input variables mainly refers to the measurement error of variables such as diameter at breast height (DBH) and tree height (H), which is primarily affected by the measurement techniques, equipment, and anthropogenic measurement methods [31]. The incorrect setting of the model form is mainly due to the lack of appropriate test data or modeling techniques [32]. For the residual uncertainty of the model, the standard deviation of residual was used to measure it [33]. However, different results may be produced due to the differences in the study area and data. Compared to the residual uncertainty of the model, the uncertainty caused by the parameter error is relatively small and relates to the number of modeling samples. Current studies focused more on the uncertainty of the model. For example, Eduarda et al. [34] reduced the uncertainty of the biomass model in terms of modeling methods. Huang et al. [35] used the error transfer method to measure the uncertainty of parametric and non-parametric models. Cao et al. [36] combined the random forest and Monte Carlo simulation (RF-MC) to study the effects of models with different types of variables on the accuracy of carbon stock estimation.

In summary, optical and SAR data to measure uncertainties arising from estimation models has been widely used [37]. However, these studies neglected how different factors impact the accuracy and uncertainty of model. Furthermore, few studies have combined optical and SAR data to quantitatively assess the uncertainties generated by two different scales of plot and remote sensing model. To address these research issues, the Pinus densata in Shangri-La was used as the research tree species, Sentinel-1/2 data and field data were used to establish different estimation models. Then, error equations were used to calculate uncertainties of plot scale and remote sensing models. The main objectives of this study were to: (1) investigate the effects of different factors and plot scale on AGC estimation; (2) analyze the primary source of uncertainty at the plot scale; (3) explore the potential of SAR data in reducing uncertainty in the estimation process.

2. Materials and Methods

2.1. Study Area

Shangri-La is located in northwestern Yunnan Province in southwest China, with a total area of 11,613 km² (Figure 1b). It has a large undulating topography, with an elevation difference of 4,042 m (the lowest elevation is 1,503 m) and an average elevation of 3,459 m (Figure 1d). The rainfall time concentrated from June to September, mean annual precipitation is 607 mm and the average annual evaporation is 1643.6 mm [38]. It has four towns and seven villages with a high forest cover of 75%. It is one of the world’s top ten species gene pools [39]. Its main vegetation types are cold-temperate coniferous forest and rubber forest. Pinus densata is one of its dominant tree species (Figure 1c), which covers an area of 1848.18 km² and accounts for 16.18% of the area of Shangri-La [40].

The methods of this paper are divided into three stages (Figure 2). In the first stage, remote sensing images were processed and remote sensing features were extracted, AGC at each field plot were calculated from the field survey. AGC remote sensing-based estimation models were established in the second stage using the RF method, measurement uncertainty and uncertainty of the mono-carbon stock model at the plot scale were calculated. In the third stage, we used error equations to calculate the remote sensing-based estimation model uncertainty and the total uncertainty at the plot scale. Then, we combined two uncertainties to obtain the total uncertainty of the AGC estimation.

2.2. Data Acquisition and Processing

2.2.1. Field Survey

A total of 60 plots were surveyed in 2019 and 2021. Twenty plots were surveyed in December 2019 and the other forty were surveyed in May 2021. Plot sizes of 10 m x 10 m were considered in this study to ensure correspondence between the field measurement and pixel size in the imagery [41]. Sample points were randomly distributed in pure Pinus densata forest areas; each plot was spaced more than 3 km apart. The field survey mainly recorded DBH and H of pure stand of Pinus densata with DBH of 5 cm or more. The binary lumber volume table was used to calculate the capacity of individual trees in each plot [42]. The AGB of Pinus densata was calculated using the allometric growth equation [43]:

$$W=0.073\times D^{1.739}\times H^{0.880}$$

(1)

where W is AGB, D is the diameter at breast height, and H is the tree height. Then the carbon stock is calculated as:

$$C=W\times f_c$$

(2)

where $C$is the AGC of Pinus densata, and $f_c$ is the carbon content coefficient, which is 0.5131 [44]. The field survey data are shown in Table 1.

2.2.2. DEM

The Digital Elevation Model (DEM) used in this study was downloaded from the official website of the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov). It is the Shuttle Radar Topography Mission (SRTM) Version 3 product, with a resolution of 30 m. We resampled it to 10 m using the nearest neighbor method.

2.2.3. Remote Sensing Images and Processing

In this study, we selected two scenes of Sentinel-1 Interferometric Wide Mode (IW) Ground Range Detected (GRD) images acquired in November 2021 with a resolution of 10 m (Table 2). Pre-processing operations such as orbit correction, thermal noise removal, radiometric calibration, speckle filtering, and terrain correction were performed in the SNAP software provided by the European Space Agency (ESA) [45]. Owing to the significant altitude disparities in the study area, we utilized the altitude, slope, and aspect extracted from DEM to correct the image for alleviating the effect of elevation differences on the images [46]. Concerning the influence of the wave propagation path on the images, the alteration in the wave propagation path can be analyzed by extracting terrain information (altitude, slope, and aspect) from DEM; we then used SARscape software to correct this impact through the alteration [47]. Finally, the processed Sentinel-1 images were exported to ENVI format, stitching and cropping operations were performed.

Sentinel-2 is a multispectral imaging satellite that carries a Multispectral Imager (MSI) with 13 spectral bands. The image is unique in that it has three special red-edge bands, making it particularly effective in monitoring vegetation information. We selected six scenes of Sentinel-2 Level-2A (L2A) images in November 2021 (Table 2). It has been atmospherically and orthographically corrected [48]. Then, all 20 m spectral bands were resampled to 10 m using the nearest neighbor method [49]. Finally, the band synthesis, splicing, and cropping were completed in the ENVI software.

2.2.4. Extraction and Selection of Remote Sensing Factors

Remote sensing factors were extracted from the Sentinel-1/2 images after pre-processing through SNAP, SARscape, and ENVI software. Studies have revealed that backscatter coefficients and texture features of SAR were crucial for AGB estimation [50,51]. Hence, we extracted SAR backscatter coefficients (VV, VH) and SAR texture features from Sentinel-1 images. Key parameters such as spectral textures (excluding atmospheric bands B1, B9, and B10), vegetation indices, and texture features were extracted from Sentinel-2 images according to the relevant research [52]. The parameters extracted in this study contained 2 SAR backscatter coefficients, 10 spectral textures, 9 vegetation indices, 16 SAR texture features, and 152 optical texture features, with a total of 189 factors (Table 3). Then, we used the SPSS software to analyze the correlation of remote sensing factors. Finally, high correlation factors were used to construct AGC estimation models. These factors mainly included 3 spectral textures, 4 vegetation indices, 5 optical texture features, and 2 SAR texture features.

2.3. Modeling Method

Random Forest (RF) [53] is an integrated learning model that is widely used in a variety of fields such as classification and nonparametric regression prediction [54]. The implementation of the RF model in this study is based on the “Random Forest Regressor” algorithm provided by the “Scikit-learn” package for the Python language. In this study, different remote sensing-based estimation models were constructed by the RF. The components of each model are BIM: spectral features, VIM: vegetation indices, TIM: optical texture features, S-2M: selected optical remote sensing factors, and S-1/2M: SAR texture features and selected optical factors (Table 4).

2.4. Accuracy Evaluation

70% (42 groups) of the plots were randomly selected for modeling, the remaining 30% (18 groups) were used for validation. Cross-validation was performed in the process of modeling. The evaluation indexes for the accuracy of the estimation model are coefficient of determination (R²), root-mean-square error (RMSE), relative root-mean-square error (rRMSE), and prediction accuracy (P). To ensure that the model results were as objective as possible, each model was fitted ten times in this study to allow the mean values of evaluation indicators to be used for comparison [55]. The formulas are as follows:

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

(3)

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

(4)

$$rRMSE={\displaystyle \frac{RMSE}{\overline{y}}}\times 100\%$$

(5)

$$P={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left(1-\left|{\displaystyle \frac{y_i-{\widehat{y}}_i}{{\widehat{y}}_i}}\right|\right)\times 100\%$$

(6)

where $y_i$ is the observed value, ${\widehat{y}}_i$ is the predicted value, $\overline{y}$ is the mean of the observed values and $n$ represents the plot number.

2.5. The Analysis of Uncertainty at the Plot Scale

2.5.1. The Uncertainty of Measurement

The observed values of diameter at breast height and tree height were denoted as D and H, their errors from the normal distribution were denoted as ${\epsilon }_D,{\epsilon }_H$, and their standard deviation of error were ${\sigma }_D,{\sigma }_H$. The observed carbon stock model is $C=f(E,a)$, $E$ are the independent variable (DBH, H) and $\alpha$ denotes the parameter of the model. According to Taylor’s first-order expansion principle, the error of the mono-carbon stock model is:

$$C_{bia}\approx f\left(e\right)-f\left(E\right)\approx f^{\prime }\left(E\right)(e-E)$$

(7)

where $f^{\prime }(E)={\displaystyle \frac{\partial f(E,\alpha )}{\partial E_1}}$ is the partial derivative of the function $f\left(x\right)$ concerning the variable $E_1$. The variance matrix of prediction error for the model can be expressed as:

$${\sigma }_{C_{bia}}^2\approx {G\left(f^{\prime }\right(E)·(e-E\left)\right)}^2=f^{\prime }\left(E\right){Var}_{e-E}{f}^{\prime }{(E)}^T$$

(8)

where ${Var}_{en-En}$ is the covariance matrix of $e_n$-$E_n$. If the error between the independent variables is uncorrelated, there is:

$${Var}_n=\left[\begin{array}{ccc}Var(e_1-E_1)& 0& \begin{array}{cc}\dots & 0\end{array}\\ 0& Var(e_2-E_2)& \begin{array}{cc}\dots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ \dots \end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ 0\end{array}& \begin{array}{c}⋮\\ Var(e_n-E_n)\end{array}\end{array}\end{array}\right]$$

(9)

For the $f(E,a)$, there are two types of models: unary model and binary model. The form of $f(E,a)$ used in this article is binary model ($f(E,a)=a{D}^{\alpha }{H}^{\beta }\times f_c$), and according to the above equation:

$${\sigma }_{C_{bia}}^2=f^{\prime }\left(E\right)Varf\left(E\right)=\left(\begin{array}{cc}{\displaystyle \frac{\partial C}{\partial D}}& {\displaystyle \frac{\partial C}{\partial H}}\end{array}\right)\left(\begin{array}{cc}{\sigma }_D^2& {\sigma }_{DH}\\ {\sigma }_{DH}& {\sigma }_H^2\end{array}\right)\left({\displaystyle \frac{\frac{\partial C}{\partial D}}{\frac{\partial C}{\partial H}}}\right)$$

(10)

$$=C^2\left({\alpha }^2\right.{\displaystyle \frac{{\sigma }_D^2}{D^2}}+{\beta }^2{\displaystyle \frac{{\sigma }_H^2}{H^2}}+2\alpha \beta {\displaystyle \frac{{\sigma }_{DH}^2}{DH}})$$

(11)

where ${\sigma }_D^2$ is the variance of error in DBH, ${\sigma }_H^2$ is the variance of error in H, and ${\sigma }_{DH}^2$ is the covariance of ${\sigma }_D$ and ${\sigma }_H$. $\alpha$, $\beta$ are the parameters of the mono-carbon stock model. The grouping method [32] was used to calculate ${\sigma }_D^2$ : 1) The means of the observed DBH were ranked in ascending order; 2) The observed DBH were grouped on the basis of 10 trees by groups; 3) Calculate the mean value of DBH in each data group as well as the standard deviation between measurements; 4) Measurement error in DBH was calculated by fitting a linear model between standard deviation of measurement error and DBH. The formulas are as follows:

$$\overline{D_I}={\displaystyle \frac{1}{n}}{\sum }_{j=1}^n\overline{D_j}$$

(12)

$${Dif}_D=D_1-D_2$$

(13)

$${\sigma }_D={\displaystyle \frac{1}{n}}\sqrt{{\sum }_{j=1}^n{({Dif}_{D,J}-\overline{{Dif}_{D,j}})}^2}$$

(14)

where $\overline{D_I}$ is the mean value of DBH in each group, $\overline{D_j}$ is the mean value of the jth sample tree in each group, ${Dif}_D$ is the difference between $D_1$ and $D_2$, ${\sigma }_D$ is the standard deviation of the difference in each group of data, and the square is ${\sigma }_D^2$ in equation (11). $\overline{{Dif}_{D,j}}$ is the mean value of ${Dif}_{D,J}$ in each group. The linear model can be fitted by the relationship between the $\overline{D_I}$ and the corresponding ${\sigma }_D$ in each group. The expression is given below:

$${\sigma }_D=a+b\overline{D_I}$$

(15)

Since tree height can only be measured by the altimeter, the difference between the measured values and the true values is large in several cases. The average of multiple measurements is still far from the true value [56]. Therefore, the absolute error between the measured value and the true value is used as ${\sigma }_H$ and the curve is fitted to estimate the error of single tree height. The expression is as follows:

$${\sigma }_H=a{H}^2+bH+c$$

(16)

where ${\sigma }_H$ is the error of tree height, $H$ is measured value of tree height, and $a$, $b$, $c$ are parameter of the model. ${\sigma }_{DH}^2$ can be calculated by ${\sigma }_D$ and ${\sigma }_H$. The ${\sigma }_{DH}^2$ can be expressed by the following equation:

$${\sigma }_{DH}^2={\displaystyle \frac{\sum (D_i-\overline{D})(H_i-\overline{H})}{N-1}}$$

(17)

where ${\sigma }_{DH}^2$ is the covariance of ${\sigma }_D$ and ${\sigma }_H$, $D_i$ is the $i$th measured value of DBH, $\overline{D}$ is mean value of DBH, $H_i$ is the $i$th measured value of $H$, $\overline{H}$ is mean value of $H$, and $N$ is dimensionality.

2.5.2. Uncertainty in the Mono-Carbon Stock Model

(1) Residual Uncertainty of the Mono-carbon Stock Model

To calculate the residual uncertainty of the mono-carbon stock model, we used the observed data of Pinus densata carbon stock to fit the mono-carbon stock model. Then the uncertainty caused by residual variation is calculated by the difference between the predicted and observed values. Let the observed value of carbon stock be $C$ and the prediction model of carbon stock be $C^{\prime }=f(E,\widehat{\alpha })$, the residual $\epsilon$ can be expressed by the following equation:

$$\epsilon =f(E,\widehat{\alpha })-C$$

(18)

where $\epsilon$ denotes the difference between the observed and predicted values, $\widehat{\alpha }$ is the predicted value of the parameter in the carbon stock model and $E$ is the variables of model (DBH, H).

The standard deviation of residual can be used to measure the uncertainty caused by the variation of model residual. It was found that the standard deviation of residual was linearly related to the AGC [33]. Therefore, the standard deviation can be calculated by fitting a linear relationship between the standard deviation and carbon stock. The equation is as follows:

$${\sigma }_{\epsilon }=\beta f(E,\widehat{a})$$

(19)

where ${\sigma }_{\epsilon }$ is the standard deviation of residual, $\beta$ is the fitting parameter. The six-step approach [57] was used to calculate the ${\sigma }_{\epsilon }$. The steps include: 1) rank the observed carbon stock values ($y$) of the plots in ascending order; 2) calculate the residuals of model ($C_{\epsilon }$), which are the difference between the observed and predicted values; 3) divide the modeling data into groups of N, with a total of 10 groups; if the last group is less than N, the remaining data is counted in the previous group; 4) calculate the mean values of the predicted carbon stock ($\overline{\widehat{y}}$) and the standard deviation of residual (${\sigma }_{\epsilon }$) in each group. The mean value, residual, and standard deviation of residual were calculated as follows:

$$\overline{\widehat{y}}={\displaystyle \frac{1}{n}}\sum _{j=1}^n{\widehat{y}}_j$$

(20)

$$C_{\epsilon }=y-\widehat{y}$$

(21)

$${\sigma }_{\epsilon }={\displaystyle \frac{1}{n-1}}\sqrt{\sum _{j=1}^n{\left(C_{{\epsilon }_j}-{\overline{C}}_{{\epsilon }_j}\right)}^2}$$

(22)

where ${\widehat{y}}_j$ denotes the predicted value of carbon stock in the $j$th plot, $C_{{\epsilon }_j}$, ${\overline{C}}_{{\epsilon }_j}$ denote the residual and the residual mean value of the $j$th plot and $n$ denotes the number of samples.

5) fitting the predicted mean value $\overline{\widehat{y}}$ and standard deviation of residual ${\sigma }_{\epsilon }$. The relationship can be expressed as:

$${\sigma }_{\epsilon }=\theta \left(\overline{\widehat{y}}\right)$$

(23)

6) substitute the predicted values of carbon stock into the fitted formula. Calculate the standard deviation of the residuals for each plot. Then, divide them by the sum of the observed carbon stock values to obtain the model residual’s uncertainty.

(2) Parameter Uncertainty of the Mono-carbon Stock Model

The Taylor series first-order expansion was used to measure the uncertainty of model parameter. The carbon stock model is first-order expanded by a Taylor series as follows:

$$C^{\prime }=f(E,\widehat{a})\approx f(E,a)+{\displaystyle \frac{\partial f\left(E,\widehat{\alpha }\right)}{{\partial }_{a_j}}}(\widehat{\alpha }-\alpha )$$

(24)

where $f(E,\alpha )$ is the observed value of carbon stock, $E$ is the independent variable, $\widehat{\alpha }$ is the simulated value of parameter, $\alpha$ denotes the true value of parameter, ${\displaystyle \frac{\partial f(E,\alpha )}{\partial \alpha j}}$ is the partial derivative of model parameter ${\alpha }_j$, ${\displaystyle \frac{\partial f(E,\alpha )}{\partial \alpha j}}(\widehat{\alpha }-\alpha )$ is the error (${\sigma }_{\gamma }$) of model parameters. ${\sigma }_{\gamma }$ can be approximated as:

$${\sigma }_{\gamma }^2\approx Z_{jk}var(a)Z_{jk}^T$$

(25)

where $Z_{jk}$ denotes the matrix with $Z$ as $j\times k$, $Z_{jk}^T$ is the transpose matrix of $Z_{jk}$ and $var\left(\alpha \right)$ is the covariance matrix of estimated parameter $\alpha$ in the carbon stock equation.

2.5.3. Uncertainty Synthesis

The error in indirect measurement can be expressed by synthesizing each direct measurement error [58]. Therefore, the total error ($M_Z$) can be described as:

$$M_Z=\sqrt{M_{{\mathrm{j}}_1}^2+M_{{\mathrm{j}}_2}^2+M_{{\mathrm{j}}_3}^2+\dots +M_{{\mathrm{j}}_{\mathrm{i}}}^2}$$

(26)

where $M_Z$ denotes the total error, $M_{{\mathrm{j}}_{\mathrm{i}}}$ denotes the error of the ith variable. Based on this equation, the total uncertainty (${\sigma }_t$) is calculated as:

$${\sigma }_t=\sqrt{{\sigma }_{\epsilon }^2+{\sigma }_{\gamma }^2+{\sigma }_c^2}$$

(27)

where ${\sigma }_{\epsilon }$ denotes the residual uncertainty of the mono-carbon stock model, ${\sigma }_{\gamma }$ denotes the uncertainty of model parameters and ${\sigma }_c$ indicates the measurement uncertainty. The error transfer from the mono-carbon stock model to the plot scale uncertainty (${\sigma }_p$) is calculated as follows:

$${\sigma }_p=\sqrt{{\sum }_{i=1}^N{\sigma }_{ti}^2/S}$$

(28)

where $N$ is the number of sample wood, ${\sigma }_{ti}$ is the total uncertainty of the mono-carbon stock model for the $i$th plot wood plant in the plot, and $S$ is the area of the plot.

2.6. Uncertainty Analysis of Remote Sensing-Based Estimation Models

2.6.1. The Uncertainty of Model Residual

We also used the six-step approach [57] to calculate the residual uncertainty in remote sensing-based estimation models; the detailed steps and formulas are shown in the six-step approach section of Part 2.5.2, Point 1. Finally, the uncertainty of the remote sensing model (${\sigma }_R$) can be expressed as:

$${\sigma }_R={\sum }_{i=1}^n{\sigma }_{{\epsilon }^{\text{'}}}/y_t$$

(29)

where ${\sum }_{i=1}^n{\sigma }_{{\epsilon }^{\prime }}$ is the sum of the standard deviation of the residuals in all plots, and $y_t$ is the sum of the measured AGC.

2.6.2. Total Uncertainty in Carbon Stock Estimation

The total uncertainty ($\sigma$) of AGC estimation in Pinus densata was calculated as follows:

$$\sigma =\sqrt{{\sigma }_p^2+{\sigma }_R^2}$$

(30)

where ${\sigma }_p$ is the uncertainty at the plot scale; ${\sigma }_R$ is the uncertainty of the remote sensing model.

3. Results and Analysis

3.1. Comparison of Accuracy for Estimation Models

Figure 3-4 show the scatter plots of the prediction accuracy for each model and the distribution of R² and RMSE values for different models. It showed that the single-band factor directly involved in modeling has the lowest accuracy with R² of 0.551, RMSE of 23.049 t/ha, and rRMSE of 45.158% (Figure 3a). The prediction accuracies of the models were all improved after we gradually introduced different remote sensing factors during the modeling. Since optical image is prone to light saturation in areas with high biomass and carbon stock, it has a disadvantage in estimating high AGC. Therefore, after introducing Sentinel-1 SAR data based on the S-2M, the model has the highest prediction accuracy with R² of 0.872, RMSE of 10.821 t/ha, and rRMSE of 21.201%. Sentinel-1 SAR data can overcome shortcomings of the optical data to a certain extent according to the results, which have an important role in improving the estimation accuracy [59]. The above results can be more intuitively shown in Figure 4. R² of the models gradually increases, and the RMSE gradually decreases with the introduction of different remote sensing factors; the model with joint Sentinel-1 SAR and Sentinel-2 optical data has the largest R² and the smallest RMSE. This provides some references for improving the accuracy of carbon stock estimation and the rational selection of variables.

3.2. Uncertainty at the Plot Scale

3.2.1. The Uncertainty of Measurement

Scatter plots of measurement error of DBH can be obtained by equations (12) and (14) (Figure 5a). It showed that the measurement error of DBH increases with the increase of DBH, which is basically a linear relationship. The error model of DBH is:${\sigma }_D=0.1913D-0.4964$, $R^2=0.665$. Then, the true values of DBH were substituted into the fitted model to obtain the error of DBH. For the measurement error of H, we can calculate by fitting a model through the relationship between the H and the measurement error of H. As shown in Figure 5b, similar to the measurement error of DBH, the error of H is also positively correlated with the H and its error form can be expressed as:${\sigma }_H=0.0095H^2-0.075H+1.3359$, $R^2=0.265$. After obtaining the error of H, ${\sigma }_{DH}^2$ is calculated according to equation (17), which is 0.70. Finally, the measurement uncertainty is calculated by equation (11). The calculation shows that for the binary carbon stock model, the uncertainty in carbon stock estimation due to measurement error is 3.02%.

3.2.2. Uncertainty in the Mono-Carbon Stock Model

(1) Residual Uncertainty

Figure 6a shows the fit of the binary carbon stock model for Pinus densata, the equation is $C=0.101D^{1.859}{H}^{0.635}\times f_c$ and the R² of the model is 0.984. Figure 6b shows the scatter plot of the error equation for the mono-carbon stock model, from which it can be seen that the standard deviation of residual increased with the increase in the predicted value of carbon stock. The error equation of the binary carbon stock model is y=0.0613x-0.369 and R²=0.944, which has an excellent fitting effect. The predicted values of carbon stock were then substituted into the error equation to calculate the total error value. Finally, the residual uncertainty of model was obtained by dividing the total error value into the total observed value of carbon stock, which was 5.96%.

(2) Parameter Uncertainty

For the parameter uncertainty of the binary carbon stock model, we first used Taylor’s first-order expansion to obtain the variance-covariance matrix of parameters. Then, the parameter uncertainty of the model was calculated using the formula. By calculating, we can know that the parameter covariance matrix of the binary carbon stock model is $\left(\begin{array}{ccc}0.00000075& -0.00003400& 0.00009900\\ -0.00003400& 0.01100000& -0.01100000\\ 0.00009900& -0.01100000& 0.04100000\end{array}\right)$, which substituted into equation (25) to conclude that the parameter uncertainty of binary carbon stock model accounts for about 6.86%. Synthesizing the two uncertainties, the total uncertainty of the binary carbon stock model is 9.09%.

3.2.3. Total Uncertainty at the Plot Scale

The total uncertainty is 9.58% by synthesizing the measurement uncertainty and the uncertainty of the mono-carbon stock model through the equation (27). When the uncertainty of the mono-carbon stock model is superimposed on the plot scale, we used the equation (28) to calculate the uncertainty of the plot scale, which is 3.71%. Table 5 shows the various sources of uncertainty in carbon stock estimation at the plot scale.

3.3. The Uncertainty of Remote Sensing-Based Estimation Models

Figure 7-9 illustrate each model’s error equations, error value of plots and uncertainties. The uncertainties of BIM, VIM, TIM, S-2M, and S-1/2M are 21.98%, 13.63%, 12.01%, 10.14%, and 7.57% by substituting the average values of carbon stock into error equations of each model (Figure 9). Through the Figure 7, the width of the confidence ellipse for each model can reflect the magnitude of the model’s R² to a certain extent, R² gradually increases as the ellipse becomes narrower. To more intuitively reflect the difference between the predicted values and the measured values for each model, we compared the observed value and the corresponding predicted value of 60 field sites (Figure 8). BIM, VIM, and TIM have phenomenon of the overestimation of low values and the underestimation of high values and the maximum absolute errors of three models are 61.14 t/ha, 41.847 t/ha, and 36.047 t/ha, which suggests that models modeled using only a single variable have larger error. Then, we filtered the three types of remote sensing factors to establish S-2M. As we can see the error curve of this model is significantly slower compared to the previous three models and its error value is also reduced (Figure 8b), as well as the absolute value of its maximum error is 28.477 t/ha. This indicates that the joint use of different types of remote sensing factors is necessary and makes up for the disadvantage of single-factor modeling, which plays an important role in and reducing uncertainty [60]. However, optical data suffers from light saturation, which impacts the accurate estimation of AGC [61]. Therefore, we added SAR data to S-2M in the following experiments to further explore the potential of SAR data in reducing uncertainty. The trend of the error curve of the model after adding SAR data is further slowed down and its error is also reduced compared with S-2M (Figure 8b); the maximum absolute value of error is 19.808 t/ha and most of the absolute values of error are below 10 t/ha, which indicates that the predicted values are close to the measured values. The comparative analysis of the errors for different models shows that the single-factor model has the largest error, and the joint use of different factors is necessary. SAR data plays an important role in improving estimation accuracy and reducing errors compared to optical data.

3.4. Total Uncertainty in Carbon Stock Estimation

The total uncertainty of AGC estimation for Pinus densata is shown in Table 6. After synthesizing the uncertainties of plot scale and estimation models by equation (30), the total Uncertainty of BIM, VIM, TIM, S-2M, and S-1/2M are 22.29%, 14.12%, 12.56%, 10.79%, and 8.43%. The BIM has the highest uncertainty, while S-1/2M with the introduction of sentinel-1 SAR data has the lowest uncertainty compared with the S-2M. Thus, SAR data have some advantages in reducing the uncertainty of AGC estimation.

4. Discussion

4.1. Analysis of the Uncertainty at the Plot Scale

Although many studies have been conducted on forest carbon stock estimation, there is still a great deal of uncertainty in the current estimation process. In this study, the uncertainty of measurement and residual variation and parameter error of the mono-carbon stock model were calculated by error equations. The uncertainties in the measurement and the mono-carbon stock model were then combined to obtain the uncertainty of plot scale. It was found that the error of the mono-carbon stock model is the most significant source of uncertainty among the uncertainties at the plot scale. Even though the measurement uncertainty is so small that its contribution to the total uncertainty is practically negligible, we identified a small systematic effect that should be considered. This error may have a greater impact if we use different personnel or equipment to collect the field data [62,63,64]. The error of the mono-carbon stock model mainly originated from parameter variation and the residual uncertainty was relatively low, consistent with the results of Huang et al. [35]. However, the parameter error of the model was higher than those of Chen et al. [65] and Wang et al. [66]. The reason is that Chen et al. [65] used 4004 plot trees for modeling, which was much higher than the number of plot trees in this study. When the number of plot trees is large enough, the uncertainty due to the model parameter can converge to zero. Wang et al. [66] used grouping and combined with the Monte Carlo method to measure the errors of unitary and binary biomass models. This method can effectively reduce the variability of parameter covariance matrices. In this paper, we only calculated the parameter uncertainty and residual uncertainty of the binary model; the uncertainty of the unitary carbon stock model needs to be further investigated in the future.

4.2. Influence of Plot Size and Forest Stand Factors on AGC Estimation and Uncertainty

The plot size used in this study is 10 m x 10 m to keep up with the spatial resolution of remote sensing images. Although several studies have shown that uncertainty decreases with the plot size increases [67,68], we need to consider the available space and representative sampling, etc. during the selection of plots. Enough plots were needed to arrange in the limited space to meet the need of the modeling because of the limitation of the study area, which can ensure the representativeness and accuracy of the Pinus densata samples. Even though estimation accuracy and uncertainty can be affected by small-scale plots, we can effectively improve this condition by increasing the number of plots and modeling samples [69]. Secondly, representative sampling is one of the important reasons why we chose small-scale plots. The difference in vegetation density and age structure also affects the estimation of AGC due to the heterogeneity of the distribution of Pinus densata in the study area. Therefore, setting more small-scale plots on this basis gives a better analysis of the factors affecting the accuracy of AGC estimation, thus improving the accuracy and reducing the uncertainty. Canopy density and average age of forest stand factors are also important factors that affect the accuracy and uncertainty of AGC estimation [70]. Previous study has shown that with the gradual increase of stand age, the forest carbon stock also presents a trend of increase [71]. Still, the growth rate decreases slowly with the maturity of the forest stand and tends to balance in the later stage of development [72]. Expect the forest stand age, canopy density is also an important factor affecting the change of AGC and can determine forest structure and vegetation growth to a large extent [73]. An increased canopy opening increased the potential of height growth, thereby increasing the AGC. Therefore, we should consider these two factors to analyze their specific effects on the change and uncertainty of AGC in the subsequent research.

4.3. Uncertainty of Different Remote Sensing Models

In this paper, we constructed different remote sensing models based on Sentinel-1/2 images. Data saturation in optical data is a key issue limiting AGC estimation accuracy, especially in dense forests with high AGB or AGC [74]. The reason is that optical data provide only spectral and horizontal spatial characteristics [75], limiting the effectiveness of spectral reflectance and vegetation indices for estimating high biomass or carbon stock. This is also confirmed by the high uncertainty of BIM and VIM in this study. To alleviate the effect of optical data saturation on the uncertainty of carbon stock estimation, we introduced Sentinel-1 SAR data on the basis of optical data. SAR can better obtain vertical structure parameters of forest vegetation because of its strong penetration ability to forest canopy and sensitivity to water content [76]. Furthermore, the combination of optical and SAR data has been proven to reduce the saturation in AGC estimation [77,78]; the basic principle of fusion is the synergy of spectral and backscattering features. Therefore, the uncertainty of S-1/2M effectively reduced after introducing Sentinel-1 SAR data, which shows that the combined use of optical and SAR data is an effective way to reduce uncertainty in the estimation process.

4.4. Limitations and Future Research

Firstly, since the study area of this paper covers only one county and Pinus densata is a single layer of pure forest [79], this conclusion will be altered if compared with the different types of forest structure, tree species, etc. in other regions. Consequently, the applicability of this study’s method in multi-layered stands or forests of different ages with complex structures needs to be further explored. Secondly, the uncertainty in this paper is only related to the plot scale and remote sensing-based estimation model. Still, during the total process of data collection, processing, model building, and carbon stock estimation, it inevitably generates various sources of error due to the differences in the measurement, data processing and calculation, model building methods and methods of carbon estimation. Therefore, we need to consider the effects of climate, light saturation, topography, error transmission, and the forest carbon cycle on the estimation of carbon stock. Finally, allometric growth equation may lead to uncertainty in the estimation of carbon stock. We attempted to improve the estimation accuracy by using an allometric growth equation for each tree species in the study area. However, some uncertain factors such as tree density, soil texture, and climatic conditions can affect the growth of H and DBH, thereby affecting the calculation of carbon stock. In the future, researchers should consider how precipitation, soil, and climate affect the growth of DBH and H and the calculation of carbon stock. Taking these factors into account will improve the accuracy and reduce the uncertainty in the process of carbon stock estimation.

5. Conclusion

In this study, the random forest was used to establish different remote sensing-based estimation models based on Sentinel-1/2 images and field data. The uncertainties resulted from the plot scale, and different estimation models were calculated using error equations. We then analyzed the impact of uncertainty at different scales on carbon stock estimation. (1) Different remote sensing factors have certain effects on the carbon stock estimation and the uncertainty of the model. In particular, the BIM has the lowest accuracy; the introduction of vegetation and texture factors improves the accuracy of the model, while the S-1/2M combining optical and SAR data has the lowest uncertainty. (2) The error of the mono-carbon stock model is the main uncertainty at the plot scale, which is 9.09 %. The uncertainty at the plot scale is 3.71%, which is lower than the uncertainty of all remote sensing-based estimation models. Therefore, the uncertainty of the estimation model is the primary source of uncertainty in the process of carbon stock estimation compared with the uncertainty at the plot scale. (3) The R² of S-1/2M combining optical and SAR data was improved by 0.046, and the uncertainty was reduced by 2.36% compared with the S-2M. Adding SAR data to optical data can effectively reduce the uncertainty in the estimation process. This study analyses the effect of different factors on AGC estimation and evaluates the uncertainty of plot scale quantitatively, which provides references for reducing the uncertainty in the process of carbon stock estimation.