Estimating Forest Aboveground Biomass by Combining Optical and SAR Data: A Case Study in Genhe, Inner Mongolia, China

Estimation of forest aboveground biomass is critical for regional carbon policies and sustainable forest management. Passive optical remote sensing and active microwave remote sensing both play an important role in the monitoring of forest biomass. However, optical spectral reflectance is saturated in relatively dense vegetation areas, and microwave backscattering is significantly influenced by the underlying soil when the vegetation coverage is low. Both of these conditions decrease the estimation accuracy of forest biomass. A new optical and microwave integrated vegetation index (VI) was proposed based on observations from both field experiments and satellite (Landsat 8 Operational Land Imager (OLI) and RADARSAT-2) data. According to the difference in interaction between the multispectral reflectance and microwave backscattering signatures with biomass, the combined VI (COVI) was designed using the weighted optical optimized soil-adjusted vegetation index (OSAVI) and microwave horizontally transmitted and vertically received signal (HV) to overcome the disadvantages of both data types. The performance of the COVI was evaluated by comparison with those of the sole optical data, Synthetic Aperture Radar (SAR) data, and the simple combination of independent optical and SAR variables. The most accurate performance was obtained by the models based on the COVI and optical and microwave optimal variables excluding OSAVI and HV, in combination with a random forest algorithm and the largest number of reference samples. The results also revealed that the predictive accuracy depended highly on the statistical method and the number of sample units. The validation indicated that this integrated method of determining the new VI is a good synergistic way to combine both optical and microwave information for the accurate estimation of forest biomass.


Introduction
Forest aboveground biomass (AGB) accounts for the dominant share of terrestrial biomass stocks [1]. There is a strong need for estimating forest AGB across large spatial scales. For example, estimates of forest AGB support sustainable forest management, bioenergy production and the detection of land-use change. Furthermore, the assessment of carbon stocks for global climate change modeling and initiatives such as Reducing Emissions from Deforestation and Forest Degradation (REDD) and REDD+ rely on forest AGB estimates [2,3]. During the last decades valuable remote sensing tools have been developed to measure AGB in time and space.
Various types of remotely sensed data have been used to estimate AGB [4][5][6][7][8]. Multispectral images are among the most spatiotemporally available, since multispectral spaceborne sensors are nonparametric methods. This shortage of knowledge is particularly evident when combining various remote sensing predictors from different sensors (i.e., ranging from spectral indices to predictors related to SAR backscatter intensity) [34].
In this study, we tested the synergic use of optical (Landsat 8 OLI) and radar (RADARSAT-2) remote sensing data to improve the estimation accuracy of forest biomass. First, five widely-used statistical prediction methods were applied to optical data and microwave data separately, as well as to the simple combination of independent optical and microwave data, across a range of sample sizes. In order to examine the effect of predictor variable selection on biomass estimation accuracy, the optimal variables for optical and microwave data were identified by a stepwise selection procedure and then applied to the same prediction algorithms and sample sizes used in first step. Finally, a combined VI (COVI) incorporating both optical and microwave information was proposed and the model performance based on the new COVI was evaluated by comparison with those based on all other tested variables. Our aim is to present an enhanced synergic inversion method that accurately determines forest aboveground biomass from remote sensing data.
Sensors 2016, 16,834 3 of 16 when combining various remote sensing predictors from different sensors (i.e., ranging from spectral indices to predictors related to SAR backscatter intensity) [34].
In this study, we tested the synergic use of optical (Landsat8 OLI) and radar (RADARSAT-2) remote sensing data to improve the estimation accuracy of forest biomass. First, five widely-used statistical prediction methods were applied to optical data and microwave data separately, as well as to the simple combination of independent optical and microwave data, across a range of sample sizes. In order to examine the effect of predictor variable selection on biomass estimation accuracy, the optimal variables for optical and microwave data were identified by a stepwise selection procedure and then applied to the same prediction algorithms and sample sizes used in first step. Finally, a combined VI (COVI) incorporating both optical and microwave information was proposed and the model performance based on the new COVI was evaluated by comparison with those based on all other tested variables. Our aim is to present an enhanced synergic inversion method that accurately determines forest aboveground biomass from remote sensing data.

Study Area
The study site focused on the Genhe area (50°48′N, 121°34′E, Figure 1) in the northeast of the Inner Mongolia Autonomous Region, China. It covers an area of approximately 625 km 2 (25 km × 25 km). Cold temperate climatic conditions prevail, with a rainy season occurring from June to August

Study Area
The study site focused on the Genhe area (50˝48 1 N, 121˝34 1 E, Figure 1) in the northeast of the Inner Mongolia Autonomous Region, China. It covers an area of approximately 625 km 2 (25 kmˆ25 km).
Cold temperate climatic conditions prevail, with a rainy season occurring from June to August with a mean annual rainfall of 464 mm. The mean annual temperature is approximately´5.3˝C. The forests of this study site mainly consist of native coniferous and broad-leaved species spreading along a topographically gentle terrain dominated by Betula and Larix spp.

Field Measurements and Aboveground Biomass Computation
The field data collection was conducted from 25 June to 29 August 2013. We measured tree DBH and tree height (H) and identified all tree species on 181 plots. DBH and H were used to generate in-situ AGB estimates to train the remote AGB-retrieval models and assess estimation accuracy, and they were measured using the Haglof Digitech Calliper and Vertex IV laser instrument, respectively. In this study, trees with DBH ě3 cm were measured. The measurements were done using a grid-based systematic sampling technique, with approx. 30 mˆ30 m sample plot, distributed systematically within the stand. Based on the DBH and H, the aboveground biomass was calculated for each composition (stock, branch and leaf) of the individual tree by applying species-specific allometric equations developed for the main forest in the northeast of China [40]. The biomass value for a single tree was estimated by summing up the biomass values of each tree component.

Remote Sensing Data Acquisition and Pre-Processing
The study utilized both Landsat 8 OLI and RADARSAT-2 satellite images. Landsat images are among the most frequently used medium spatial-resolution data for remote forest AGB estimation at local and regional scales [9,41]. RADARSAT-2, carrying a fully polarized SAR operating at C-Band (5.3 GHz) with a wavelength of approximately 5.6 cm, has opened up new opportunities for remote forest AGB estimation by the use of backscattering coefficients at different polarizations [21,22]. Landsat 8 OLI and RADARSAT-2 satellite images that covered the region of interest were acquired during the time that most closely coincided with the field measurement dates (i.e., the 25 August 2013 and the 30 June 2013, respectively). The Landsat 8 OLI image was obtained with 8% cloud cover, sun azimuth angle of 155.84 and sun elevation angle of 47.07. It has a resolution of 30 mˆ30 m for the non-panchromatic OLI bands. RADARSAT-2 data, which was acquired under conditions of a sunny and clear sky with less cloud cover than Landsat 8 OLI, was in a single look complex (SLC) format with a mean incident angle of 36.6. The fine quad-polarization wide (FQW) mode with 8 m nominal spatial resolution was selected for this study.
The Landsat 8 OLI image was obtained in digital number (DN) and needed to be converted to reflectance values. Preprocessing was performed in ENVI 5.1 software. Firstly, conversion from DN to Top-Of-Atmosphere (TOA) spectral radiances was implemented, following the approach described on the USGS website. Then the atmospheric correction of Landsat 8 OLI image to surface reflectance was performed using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercube (FLAASH). Prior to the analysis, the image was finally scrutinized for accurate geo-referencing using field-based GPS reference points (e.g., road intersections, bridges, rock formations). Further correction of Landsat 8 OLI image was not necessary.
Preprocessing of RADARSAT-2 data was achieved using the Next European Space Agency (ESA) SAR Toolbox (NEST) software. RADARSAT-2 image was first radiometrically calibrated to obtain the linear radar backscattering coefficients transformed from the DN. A 3ˆ3 boxcar filter was then applied to the backscattering coefficients data to reduce speckle noise. Boxcar averaging adds N adjacent samples, divides the sum by N, and then writes that value into the Nth sample location and is suitable in homogeneous areas such as large forest areas. Subsequently, the image was geometrically corrected using the pre-processed Landsat 8 OLI image with an overall error of approximately 0.65 pixels. After that, RADARSAT-2 data was resampled to the same spatial resolution of Landsat 8 OLI image (30 m). Landsat 8 OLI and RADARSAT-2 needed to cover the entire and accurate extent of the study area. The two images were thus resized, resulting in a complete and precise coverage of the test site (25 kmˆ25 km).

Regression Algorithms
Different statistical and machine-learning techniques have been applied to model complex relationships between remote sensing signals and forest biomass [32][33][34][35][36]. In our study, we compared five parametric or nonparametric prediction models in terms of their accuracy for biomass estimation: KNN, SVM, BPNN, RF and Stepwise Linear Regression (LMSTEP).
RF is an ensemble learning algorithm [42], designed to improve the classification and regression trees (CART) method by integrating a large set of decision trees based on a deterministic technique, by selecting a random set of variables and a random sample from the training dataset. RF technique has furthermore become popular in remote sensing as a nonlinear and non-parametric alternative, with promising predictive capabilities for high-dimensional datasets [43]. Fassnacht et al. [34] showed that the RF was a suitable modeling technique for the estimation of forest aboveground biomass using hyperspectral and LiDAR data.
Generally speaking, in a KNN method the weighted mean of the response value from the most similar neighbor(s) is assigned to a target unit of interest, where the similarity is defined in a feature space consisting of candidate predictor variables [44]. The KNN technique has received considerable attention due to its ease of availability, and some previous studies revealed that the KNN technique has the potential to increase the precision of vegetation parameters estimates [45][46][47]. The KNN was used and compared with other models by Tian et al. [37] with respect to their predictive power of forest aboveground biomass and a relatively good performance of KNN was found in comparison to other models. Support vector regression (SVR) is a training approach based on the framework of statistical learning theory. The main idea behind SVR is to estimate the linear dependency between pairs of n-dimensional input vectors and 1-dimensional target variables by fitting an optimal approximating hyperplane to a set of training samples. The method has proven its robustness to dimensionality and generalization ability [48]. As applied to remote sensing data, SVR has shown great capability to predict biophysical/chemical plant parameters by relating spectral information to in situ samples [49][50][51].
BPNN is one of the most popular and proven neural networks [52,53]. This machine-learning algorithm has great potential for estimating vegetation parameters due to its ability to find and learn complex linear and nonlinear relationships between radiometric data and vegetation parameters [52,54]. Zheng et al. [53] found that BPNN was better at predicting forest growing stock volume compared to LMSTEP based on Landsat Thematic Mapper (TM) and field data.
LMSTEP is a classical parametric method. The method aims at automatically choosing a set of predictive explanatory variables among all possible variables, and the partial F-statistic is a common criterion, with α as a threshold for adding or removing a predictor variable. We therefore include this technique in order to assess whether more flexible and complex nonparametric methods could outweigh the more frequently-used LMSTEP technique, which accounts for linear relationships between response and predictor variables.

Experiments
Using all the algorithms, we conducted three experiments (Experiments 1-3) with different data type combinations (Table 1). Predictors used in these experiments varied, depending upon available sensor spectral (or backscatter coefficients) data. Experiments were as follows: Experiment 1 was composed of three subtests: (i) Landsat 8 OLI image bands and image-derived variables; (ii) RADARSAT-2 backscatter coefficients and image-derived variables; and (iii) the joined variables from subtests (i) and (ii); Experiment 2 that also encompassed three subtests: (i) optical optimal variables chosen by the predictor selection method; (ii) SAR optimal variables; and (iii) the joined variables from subtests (i) and (ii); and Experiment 3: the COVI and optical and SAR optimal variables. Methods for the predictor selection and the COVI construction are detailed in Sections 2.6 and 2.7, respectively. The following is an overview about the modeling set-up in Experiments 1-3. Four main processing steps were included: After preparing an original dataset of response and predictor variables, the dataset of size n was ordered according to ascending biomass values and was then subdivided in five equal-size subsets of size ns = n/5. (II) Using the stratified dataset derived from step I, we developed subsamples for each stratum using bootstrapping [55]. In addition to the model diagnostics, predictive models were fitted to obtain spatial prediction maps of biomass for visual comparison. All statistical calculations were implemented using the R statistical package [56].

Predictor Selection
To examine the effect of predictor selection on model performances and identify more suitable predictors to obtain higher model accuracy, a stepwise procedure was used to exclude relatively poor predictors and highly intercorrelated predictors. We emulated the predictor selection method used in the study of Kattenborn et al. [35]. In each step, n models (n = number of predictors) were run using all but one predictor variable (n´1).The model validation described above in Section 2.5 was applied to each model. Therefore, the model which had the best performance (R 2 ) indirectly indicates the (excluded) predictor with the lowest explanatory power. In the next step, this predictor was removed from the set of candidate predictors. Model performance was evaluated using average values of the different model methods (KNN, SVM, BPNN, RF and LMSTEP). This procedure was continued until convergence of the R 2 value was reached ( Figure 2). The stepwise procedure was used to generate optical and SAR optimal predictor variables, respectively (i.e., OOVs in Experiment 2 (i) and SOVs in Experiment 2 (ii)). used in the study of Kattenborn et al. [35]. In each step, n models (n = number of predictors) were run using all but one predictor variable (n − 1).The model validation described above in Section 2.5 was applied to each model. Therefore, the model which had the best performance (R 2 ) indirectly indicates the (excluded) predictor with the lowest explanatory power. In the next step, this predictor was removed from the set of candidate predictors. Model performance was evaluated using average values of the different model methods (KNN, SVM, BPNN, RF and LMSTEP). This procedure was continued until convergence of the R 2 value was reached ( Figure 2). The stepwise procedure was used to generate optical and SAR optimal predictor variables, respectively (i.e., OOVs in Experiment 2 (i) and SOVs in Experiment 2 (ii)).

Combined Vegetation Index Construction
As stated before, both optical and SAR data have advantages and shortcomings in monitoring vegetation parameters. In order to overcome the limitations associated with the use of a single data type, a combined VI incorporating weighted optical and microwave information was designed and its equation can be written as: where COVI is the combined VI, OVI and MBI are the most important optical and microwave predictor variables chosen, respectively. A is the weighting factor and its calculation is elaborated below.
(1) After preparing a matrix of response and the most important optical (OVI) and microwave (MBI) predictor variables, the dataset of size n was first ordered according to ascending biomass values. OVI and MBI then needed to be normalized to allow for direct comparison between optical and microwave variables, e.g., with different scales and dynamic ranges. (2) With limited and discrete samples, the sensitivity value of the predictor variable to the response variation could be expressed via first difference quotient method. This method consists of forward, backward and central difference quotients, as an approximation of the first derivative method [57,58]. In the study, the central difference quotient approach was applied to calculate the sensitivity value of OVI (or MBI) to biomass changes. Take OVI as an example, the processing can be carried out as: where OVI_ni is the normalized value of OVI and OSi is the corresponding sensitivity value of case i.
(3) Based on the above analysis, the weights of optical data (i.e., OVI) can be formulated as follows:

Combined Vegetation Index Construction
As stated before, both optical and SAR data have advantages and shortcomings in monitoring vegetation parameters. In order to overcome the limitations associated with the use of a single data type, a combined VI incorporating weighted optical and microwave information was designed and its equation can be written as: where COVI is the combined VI, OVI and MBI are the most important optical and microwave predictor variables chosen, respectively. A is the weighting factor and its calculation is elaborated below.
(1) After preparing a matrix of response and the most important optical (OVI) and microwave (MBI) predictor variables, the dataset of size n was first ordered according to ascending biomass values. OVI and MBI then needed to be normalized to allow for direct comparison between optical and microwave variables, e.g., with different scales and dynamic ranges. (2) With limited and discrete samples, the sensitivity value of the predictor variable to the response variation could be expressed via first difference quotient method. This method consists of forward, backward and central difference quotients, as an approximation of the first derivative method [57,58]. In the study, the central difference quotient approach was applied to calculate the sensitivity value of OVI (or MBI) to biomass changes. Take OVI as an example, the processing can be carried out as: where OVI_n i is the normalized value of OVI and OS i is the corresponding sensitivity value of case i. (3) Based on the above analysis, the weights of optical data (i.e., OVI) can be formulated as follows: where A i is the weight value of OVI in case i, MS i is the sensitivity value of MBI under the i'th condition Equation (2).
In order to obtain the quantitative expression of the weight of optical (or microwave) data (i.e., A), the relationship between the weight value of optical data and corresponding OVI value was then developed by statistical regression analysis.

Experiments 1 and 2 Performances
The stepwise selection method retained seven predictors of optical sensor type, i.e., SWIR bands; SR, NDVI, OSAVI, MSI and NDWI. The microwave optimal variables chosen by the selection procedure included: VV, HV, HH/HV and RVI. In Experiments 1 and 2 performed with sole optical or SAR metrics (i.e., Experiments1 (i), (ii) and 2 (i), (ii)), the results were always not satisfactory with RMSE > 30 Mg/ha and R 2 < 0.6, regardless of what prediction algorithm or sample size was used. Additionally, the employment of optical variables generally led to slightly smaller RMSE and higher R 2 values compared to the microwave-only results. The differences in RMSE and R 2 between the best optical result and the best SAR result were about 5 Mg/ha and 0.08, respectively. Figures 3 and 4 summarize the model performances based on optical and microwave integrated dataset from Experiments 1 (iii) and 2 (iii), respectively. The use of the combined optical and microwave dataset resulted in better AGB estimates compared to the use of either the sole optical data or the sole microwave data, regardless of the prediction algorithm or sample size that was used. Furthermore, a relatively stable increase in R 2 and decrease in RMSE were observed along with increasing sample size (colored horizontal lines from left to right). RF performed best compared to all other tested methods, especially when large numbers of sample units were available, and LMSTEP, in most cases, turned out slightly higher RMSE and lower R 2 values when compared to all other methods.
From Figures 3 and 4, we can also see that using the selected optimal variables resulted in relatively higher predictive performances of models (Experiment 2) compared to those using original tested predictors (Experiment 1). For instance, for RF, the differences in RMSE and R 2 between the models based on combined optical and microwave optimal variables (i.e., OOVs + SOVs in Experiment 2 (iii) in Figure 4) and those using a combination of original optical and microwave variables (i.e., OVs + SVs in Experiment 1 (iii) in Figure 3) reached approximately 5 Mg/ha and 0.085.
From Figures 3 and 4, we can also see that using the selected optimal variables resulted in relatively higher predictive performances of models (Experiment 2) compared to those using original tested predictors (Experiment 1). For instance, for RF, the differences in RMSE and R 2 between the models based on combined optical and microwave optimal variables (i.e., OOVs + SOVs in Experiment 2 (iii) in Figure 4) and those using a combination of original optical and microwave variables (i.e., OVs + SVs in Experiment 1 (iii) in Figure 3) reached approximately 5 Mg/ha and 0.085.

Experiment 3 Performance
According to the results from previous tests, the simple combination of optical and microwave dataset (i.e., OVs + SVs in Experiment 1 (iii) or OOVs + SOVs in Experiment 2 (iii)) improved results compared to either the optical-only or microwave-only data, regardless of what sample size or prediction algorithm was used, both before and after predictor variables selection. In order to further improve biomass estimation accuracy, a new VI incorporating optical and microwave information was developed. In subsequent analyses, we focused on the preferable RF-based methodology with optical and microwave optimal variables and the largest sample size.
From the previous experiments, the most important optical and microwave variables detected by RF were OSAVI and HV, which were therefore employed in the COVI computation. More details on the method of constructing the COVI were provided in Section 2.7. When setting up the new VI, we developed a subsample using stratified bootstrapping as well: the dataset of size n (181) ranked on the basis of rising biomass values was divided into ten equal-size subsets of size ns = n/10; one-sixth of each of the ten subsets was then extracted and was subsequently gathered up to obtain the subsample. Lastly, the 100 input datasets with the largest sample size (Class 4), whose predictor variables contained the COVI and residual optical and microwave optimal variables (i.e., SWIR bands, SR, NDVI, MSI, NDWI and VV, HH/HV, RVI), were fit by RF prediction method.

Experiment 3 Performance
According to the results from previous tests, the simple combination of optical and microwave dataset (i.e., OVs + SVs in Experiment 1 (iii) or OOVs + SOVs in Experiment 2 (iii)) improved results compared to either the optical-only or microwave-only data, regardless of what sample size or prediction algorithm was used, both before and after predictor variables selection. In order to further improve biomass estimation accuracy, a new VI incorporating optical and microwave information was developed. In subsequent analyses, we focused on the preferable RF-based methodology with optical and microwave optimal variables and the largest sample size.
From the previous experiments, the most important optical and microwave variables detected by RF were OSAVI and HV, which were therefore employed in the COVI computation. More details on the method of constructing the COVI were provided in Section 2.7. When setting up the new VI, we developed a subsample using stratified bootstrapping as well: the dataset of size n (181) ranked on the basis of rising biomass values was divided into ten equal-size subsets of size ns = n/10; one-sixth of each of the ten subsets was then extracted and was subsequently gathered up to obtain the subsample. Lastly, the 100 input datasets with the largest sample size (Class 4), whose predictor variables contained the COVI and residual optical and microwave optimal variables (i.e., SWIR bands, SR, NDVI, MSI, NDWI and VV, HH/HV, RVI), were fit by RF prediction method. Figure 5a shows the sensitivities of OSAVI and HV to biomass change. From this figure we can see that the HV sensitivity values rose in the first stage, reached a maximum of about 0.004, and then slowly declined with the increasing biomass. Meanwhile, as the biomass increased, there was a dramatic decrease in the sensitivity of OSAVI. Furthermore, the OSAVI was found to have a higher sensitivity to biomass variation than the microwave backscatter coefficient (HV), but the saturation problem of the optical data was also obvious. When biomass values were less than approximately 100 Mg/ha, OSAVI was more sensitive to the change of biomass than HV and played a more important role in the COVI. In contrast, the sensitivity values of HV outperformed those of OSAVI when biomass values exceeded about 100 Mg/ha, and thus the COVI mainly depended on HV information. Therefore, the scatterplot of OSAVI (or HV) weight values and biomass values (Figure 5b) indicates that the biomass value of about 100 Mg/ha could be perceived as a threshold: optical information largely controlled the COVI when biomass values were less than 100 Mg/ha and microwave information became increasingly important accompanied by increasing biomass values. The quantitative regression relationship was finally developed using the weight of optical data (i.e., A) as the dependent variable and corresponding OSAVI as the independent variable. Overall, the mathematical expression of the COVI can be shown as:

Combined Vegetation Index
where OVI and MBI are OSAVI and HV, respectively.

Model Performances
According to the results in Experiments 1 and 2, the models (i.e., NBM in Table 2) based on a simple combination of optical and microwave optimal variables (i.e., OOVs + SOVs in Experiment 2(iii)) with RF method and the largest sample size (Class 4) performed best compared to all other tested combinations. Table 2 shows that the models in Experiment 3 (i.e., BM) using the predictor set that included the COVI and residual optical and microwave optimal variables with the same prediction method (RF) and sample size (Class 4) generated better RMSE and R 2 compared to NBM. The best-performing model reached a mean R 2 of 0.82 with a mean RMSE of 15.95 Mg/ha and RMSEr of 14.17%, and the latter achieved a mean R 2 of 0.75 with a mean RMSE of 21.03 Mg/ha and RMSEr of 18.68%. Table 2. Performances of the models (BM) based on COVI and residual optical and microwave optimal variables in Experiment 3 and the models (NBM) which performed best in Experiments 1 and 2. The quantitative regression relationship was finally developed using the weight of optical data (i.e., A) as the dependent variable and corresponding OSAVI as the independent variable. Overall, the mathematical expression of the COVI can be shown as: COVI " pp´0.846qˆOVI 2´1 .008ˆOVI`1.293qˆOVI`p1´pp´0.846qˆOVI 2´1 .008ˆOVI`1.293qqˆMBI (4) where OVI and MBI are OSAVI and HV, respectively.

Model Performances
According to the results in Experiments 1 and 2, the models (i.e., NBM in Table 2) based on a simple combination of optical and microwave optimal variables (i.e., OOVs + SOVs in Experiment 2(iii)) with RF method and the largest sample size (Class 4) performed best compared to all other tested combinations. Table 2 shows that the models in Experiment 3 (i.e., BM) using the predictor set that included the COVI and residual optical and microwave optimal variables with the same prediction method (RF) and sample size (Class 4) generated better RMSE and R 2 compared to NBM. The best-performing model reached a mean R 2 of 0.82 with a mean RMSE of 15.95 Mg/ha and RMSE r of 14.17%, and the latter achieved a mean R 2 of 0.75 with a mean RMSE of 21.03 Mg/ha and RMSE r of 18.68%. 3.3. Wall-to-Wall Predictions Figure 6 illustrates the wall-to-wall mean biomass predictions as obtained from the best-performing models (i.e., BM) based on the COVI and residual optical and microwave optimal variables using the largest number of sample units (Class 4) and RF method. The mean biomass estimates were reasonably well spread over the original range of reference biomass values.
The predicted values varied from 14.68 Mg/ha to 249.67 Mg/ha, with a mean value of 106.45 Mg/ha. In most of the study area, the biomass values were in the range from 70 to 120 Mg/ha. Furthermore, many areas with low predicted AGB were situated near villages or major roads. The rivers and cloud cover areas have been masked out and other non-forested areas could be identified as homogeneous areas of low predicted biomass estimates.

Performance of Sensor, Statistical Model and Sample Size
The results of Experiments 1 and 2 show that the sole use of either microwave data or optical data did not provide satisfactory results, and the simple combinations of independent optical and microwave data improved results compared to either the optical-only or microwave-only data, regardless of which prediction algorithm or sample size was used, before and after predictor variables selection. This is in line with earlier findings [15,59,60]. Attarchi et al. [29] evaluated 11 different multiple linear regression models using optical and L-band SAR variables, and demonstrated that only the joint use of SAR and multispectral data allowed a good estimation of AGB in those regions. The use of both optical and microwave information confirms our initial Figure 6. Wall-to-wall map of mean biomass estimates as obtained from the 100 bootstrapped model runs, using the COVI and residual optical and microwave optimal variables, random forest and largest sample size (models in Experiment 3).

Performance of Sensor, Statistical Model and Sample Size
The results of Experiments 1 and 2 show that the sole use of either microwave data or optical data did not provide satisfactory results, and the simple combinations of independent optical and microwave data improved results compared to either the optical-only or microwave-only data, regardless of which prediction algorithm or sample size was used, before and after predictor variables selection. This is in line with earlier findings [15,59,60]. Attarchi et al. [29] evaluated 11 different multiple linear regression models using optical and L-band SAR variables, and demonstrated that only the joint use of SAR and multispectral data allowed a good estimation of AGB in those regions. The use of both optical and microwave information confirms our initial hypothesis that a combination of the two types of data enhances biomass estimations. This is because in the area with relatively dense vegetation coverage, the efficiency of multispectral data is affected by the saturation phenomenon.
In contrast, C-band SAR data has penetration depth of few centimetres and is sensitive to top canopy architecture [21]. In addition, the sensitivity of C-band backscatter to forest AGB could also be partly attributed to heterogeneity in the forest canopy structures [21]. Cougo et al. [22] evaluated the relationships between backscattering of a RADARSAT-2 image and the structural attributes of regenerating mangrove vegetation, and found significant relationships between backscattering coefficients in VH, HH and VV and the average height, DBH and AGB which was mainly in the range from 40 to 200 Mg/ha. Chand et al. [21] also showed reasonable correlations between backscatter values derived from C-band SAR data and forest biometric parameters. In our study, according to the sensitivity analysis of optical VI (OSAVI) and SAR backscattering (HV) to biomass variation (Figure 5a), the sensitivity values of HV outperformed those of OSAVI when biomass values exceeded about 100 Mg/ha. Overall, we therefore conclude that C-band SAR signal, which is sensitive to canopy patterns, could provide important information that is not present in optical data, and their combination is valuable for accurate biomass estimation.
Although the objective of this study was not to compare algorithms, the study has shown that among the considered methods, RF outperformed the other tested approaches, particularly when many sample units were used. The good performance has also been found in other studies. Fassnacht et al. [34] showed that RF was better at predicting forest aboveground biomass compared to other tested methods based on hyperspectral and LiDAR data. The RF was also used and compared with other models by Latifi et al. [36] for the estimation of aboveground forest biomass and it proved to be superior to the other examined methods. We attribute the outstanding predictive accuracy to the flexibility and robustness of the RF approach, which markedly differs from all other tested methods due to its conceptual design. A problem with RF may be that the applied subsampling may lead to considerable variance of the estimates when applied to a small number of sample units (compare results regarding sample size). In most cases, LMSTEP performed worse than other tested algorithms, especially when applied to optical data. This phenomenon can be explained by the fact that relationships between optical predictors and observed biomass are likely nonlinear and thus not well fitted by LMSTEP.
Our results also indicate that the number of sample units is important for explaining the variance in R 2 and RMSE on the test site. This is in accordance with other studies [34]. Figures 3 and 4 show that a practically relevant effect on R 2 and RMSE can be attributed to sample size (compare colored stripes), especially for RF.

Combined Vegetation Index
The best-performing model achieved in this study was based on the COVI and residual optical and microwave optimal variables. In fact, some promising results have been obtained by combining optical with SAR data to estimate biomass and other vegetation parameters [29][30][31]. Usually, a simple combination of the two types of data is used. For instance, the common method is that independent optical and SAR variables are simultaneously considered as predictor variables. Parametric or nonparametric relationships are developed between these variables and the predictive parameter, which is similar to the way of OVs + SVs in Experiment 1 (iii) (or OOVs + SOVs in Experiment 2 (iii)). However, the theoretically combined mechanism has not been considered in this case. In this study, a combined VI was proposed based on the optical reflectance and SAR backscatter mechanism, which determined the contributions of optical and microwave information under different biomass conditions. In other words, vegetation indices with fixed expressions are improved by the new VI with variable weights of optical and microwave data as biomass variation. From the comparison analysis of the 5-fold cross-validation (Table 2), it can be assumed that the models (i.e., BM) based on the COVI and residual optical and microwave optimal variables outperformed the corresponding models (i.e., NBM) using a simple combination of optical and microwave optimal variables (i.e., OOVs + SOVs) and obtained the best predictive performance. The only difference between BM and NBM is the combined method of OSAVI and HV. OSAVI and HV was integrated through a simple combination in NBM, while in BM, OSAVI and HV was combined using the proposed synergic method resulting in the new VI (i.e., COVI). We therefore conclude that the COVI can yield a significant improvement in the biomass estimation compared to a simple combination of independent SAR (i.e., HV) and optical variables (i.e., OSAVI). The high accuracy in estimating AGB can be attributed to the more specific mathematical and physical significance of the COVI due to the reasonable way of incorporating weighted optical and microwave information. Therefore, the COVI has valuable potential to enhance aboveground biomass estimation accuracy. However, further tests and validation are needed before the wide application of COVI.

Overall Performance of the Biomass Model
We concede that C-band SAR data has weaker penetrability than longer wavelength radar, such as L-and P-band data, which could lead to a negative effect on capturing vertical vegetation structure information, especially in dense vegetation coverage areas. We therefore suspect that one of the reasons why the models in the study have good performances was that the reference data in this study was concentrated on low-medium biomass values, with few field samples in high value ranges (more than 200 Mg/ha). The accuracy of the biomass estimation (R 2 = 0.82, RMSE = 15.95 Mg/ha, RMSE r = 14.17%) using the models based on the COVI and residual optical and microwave optimal variables (Experiment 3) in our study is comparable to or better than those of other similar studies. For example, Kattenborn et al. [35] evaluated the potential of data fusion based on hyperspectral and InSAR data using RF for temperate forest biomass estimation with R 2 of 0.73 and RMSE r of 15%. In the study of Fassnacht et al. [34], the data combination based on hyperspectral and airborne LiDAR data was performed to estimate forest biomass with a mean R 2 of 0.72 and RMSE r of 28.46%. Laurin [61] estimated tropical forest biomass using the integration of airborne LiDAR metrics with hyperspectral bands with R 2 of 0.7 and RMSE r of 35.83%.
In conclusion, the combined VI (COVI) provides a good synergistic way to integrate optical with microwave information. Therefore, the COVI presented could be useful in biomass and other vegetation parameters estimation. This study may be valuable in guiding further research on quantifying the different responses of optical reflectance and microwave backscatter to biomass and other vegetation variables. For the future, our synergistic method should be expanded to other forest types and sensor systems, especially the introduction of longer wavelength SAR data, for further validation and wide application.