*Article* **Estimation of Forest Biomass in Beijing (China) Using Multisource Remote Sensing and Forest Inventory Data**

#### **Yan Zhu 1, Zhongke Feng 1,\*, Jing Lu <sup>1</sup> and Jincheng Liu <sup>2</sup>**


Received: 2 December 2019; Accepted: 28 January 2020; Published: 31 January 2020

**Abstract:** Forest biomass reflects the material cycle of forest ecosystems and is an important index to measure changes in forest structure and function. The accurate estimation of forest biomass is the research basis for measuring carbon storage in forest systems, and it is important to better understand the carbon cycle and improve the efficiency of forest policy and management activities. In this study, to achieve an accurate estimation of meso-scale (regional) forest biomass, we used Ninth Beijing Forest Inventory data (FID), Landsat 8 OLI Image data and ALOS-2 PALSAR-2 data to establish different forest types (coniferous forest, mixed forest, and broadleaf forest) of biomass models in Beijing. We assessed the potential of forest inventory, optical (Landsat 8 OLI) and radar (ALOS-2 PALSAR-2) data in estimating and mapping forest biomass. From these data, a wide range of parameters related to forest structure were obtained. Random forest (RF) models were established using these parameters and compared with traditional multiple linear regression (MLR) models. Forest biomass in Beijing was then estimated. The results showed the following: (1) forest inventory data combined with multisource remote sensing data can better fit forest biomass than forest inventory data alone. Among the three forest types, mixed forest has the best fitting model. Forest inventory variables and multisource remote sensing variables can match each other in time and space, capturing almost all spatial variability. (2) The 2016 forest biomass density in Beijing was estimated to be 52.26 Mg ha−<sup>1</sup> and ranged from 19.1381–195.66 Mg ha<sup>−</sup>1. The areas with high biomass were mainly distributed in the north and southwest of Beijing, while the areas with low biomass were mainly distributed in the southeast and central areas of Beijing. (3) The estimates from the RF model are better than those from the MLR model, showing a high *R*<sup>2</sup> and a low root mean square error (RMSE). The *R*<sup>2</sup> values of the MLR models of three forest types were greater than 0.5, and RMSEs were less than 15.5 Mg ha−1, The *R*<sup>2</sup> values of the RF models were higher than 0.6, and the RMSEs were lower than 13.5 Mg ha<sup>−</sup>1. We conclude that the methods in this paper can help improve the accurate estimation of regional biomass and provide a basis for the planning of relevant forestry decision-making departments.

**Keywords:** forest biomass estimation; forest inventory data; multisource remote sensing; random forest; biomass density

#### **1. Introduction**

Forest ecosystems are an important component of the terrestrial ecosystem. Forests store 76%~98% of the organic carbon in terrestrial ecosystems [1] and play an irreplaceable role in mitigating global warming caused by the increase in atmospheric carbon dioxide [2]. Forest biomass reflects the material cycle of forest ecosystems and is an important indicator for measuring changes in forest structure

and function. Additionally, forest biomass is closely related to the carbon sources and sinks in forest ecosystems [3]. Because the monitoring of forest biomass resources is expensive and time consuming, most countries do not have effective monitoring systems. Therefore, accurate estimations of forest biomass can effectively replace forest monitoring systems and are an important basis for assessing ecosystem processes, the carbon balance of ecosystems and climate change [4].

Meso-scale forest biomass estimations are usually obtained from forest inventory data [5]. In many countries, the use of large-scale forest inventories is considered an effective method for estimating biomass accurately [6]. China conducts a large-scale forest resource survey every five years to provide good data for statistical forest resources. Using these inventory data, forest biomass can be estimated at provincial or national scales [7]. However, with the continuous change in the forest resource structure of China, there have been some problems with these inventory data in regional biomass estimation [8]. To obtain the total volume or biomass of forest, the volume or biomass of one tree is calculated, and then the volumes or biomasses of all the trees in the sample plot are added together. Obviously, this method requires a high amount of manpower and material resources [9]. Moreover, inventory data cannot fully reflect forest information [6,10]. Therefore, we need data that cover a wide area and contain a high amount of vegetation information to supplement forest inventory data.

With developments in technology, remote sensing has increased the possibilities for forest biomass research [11,12]. The use of remote sensing data in the research of meso-scale biomass is an important technical method. Various remote sensing indicators based on optical sensors, such as the normalized difference vegetation index (NDVI) and other factors obtained by image transformations, have been shown to be well correlated with the ground vegetation, providing reliable information for forest biomass estimation [13–15]. However, applications with optical data are often limited due to the complexity of biomass in time and space and limitations in the spatial and spectral characteristics of satellite data [16]. More abundant remote sensing data are needed to depict detailed forest information. Lidar can penetrate dense forests, provide accurate three-dimensional information of trees, and then be used to obtain forest biomass [17]. However, because of its limited coverage, high cost and inconvenience to transport, Lidar is not suitable for forest biomass estimation at the meso-scale [18]. Synthetic Aperture Radar (SAR) data, such as L-band Advanced Land Observing Satellite/Phased Array L-band Synthetic Aperture Radar (ALOS/PALSAR) [19] and X-band TerraSAR-X data, are widely used in the estimation of forest biomass [20,21]. SAR is not affected by illumination and climate conditions and it can penetrate vegetation to obtain information, covering relatively large areas in a short period of time [20].

At the meso-scale, many studies have demonstrated the potential of optical and radar remote sensing-derived indicators to estimate forest biomass [22,23]. However, there is a large range and many uncertainties of remote sensing. For example, the resolution of remote sensing images might be insufficient, and the vertical structure information of forest canopies cannot be obtained, which has certain limitations in high biomass areas. Therefore, at the regional scale, the accuracy of forest biomass estimation using remote sensing data is low [24].

Optical and radar remote sensing data can match forest inventory data in time and space [25]. In addition, these data can provide forest attributes and structural information that are missing from inventory data. Therefore, combining multisource remote sensing data with forest inventory data for regional forest biomass research provides a more consistent spatial and temporal analyses than forest inventory data alone. Generally, the uncertainty in the estimation can be reduced by this combination [25]. Furthermore, this combination can promote the application of forest biomass estimation and other ecological research at the meso-scale. However, there is little information on the potential of the combination of sample plot survey data with multisource remote sensing data to estimate and map biomass. Most forest biomass estimation studies focus on the impact of environmental variables on forest biomass [26]. Therefore, it is necessary to better assess and understand the modeling potential of sample survey factors and remote sensing factors to provide decision makers with information on forest resources.

In this study, by synthesizing the existing technical means, we combined forest inventory data with multisource remote sensing data to estimate forest biomass and improve the accuracy of biomass estimation at the regional scale. There are three objectives of the present study: the primary objective is to assess the potential of forest inventory data combined with multisource remote sensing data in modeling and mapping forest biomass. The second objective is to estimate the biomass of different forest types in Beijing in 2016 and provide data support for regional biomass estimation. The third objective is to estimate biomass using multiple linear regression (MLR) and the random forest (RF) model and compare the performances of the two models.

#### **2. Materials and Methods**

#### *2.1. Data Collection*

#### 2.1.1. Forest Inventory Data

The National Forest Resources Continuous Inventory system is a method of forest resource investigation that aims to understand the status and dynamics of macroforest resources and periodically reviews both with fixed sample plots. It is an important part of the comprehensive monitoring system for forest resources and ecological conditions in China. China's Ministry of Forestry has carried out eight consecutive national surveys and inventories of forest resources [27]. According to the technical regulations of the national forest resources continuous inventory, systematic sampling is used to lay out fixed sample plots, the size of which is 4 km × 4 km, and the sample plots are laid out at the intersection point of the kilometer network of the newly compiled 50,000 or 100,000 topographic map of the country. To ensure that the sample points are not repeated and missed, computer technology such as GIS, is used as far as possible [28]. In recent years, the collection of forest inventory data depends primarily on manual work, and it is supported by high-tech survey instruments that can automatically collect data to improve the accuracy of inventory results [29].

In this study, we used the Ninth Beijing Forest Inventory data of 2016, which involve 1431 sample plots located in all districts and counties of Beijing, as shown in Figure 1, covering coniferous, broadleaf and mixed coniferous-broadleaf forest types. The dataset describes in detail plot locations, measurement dates and forest compositions. For each plot, multiple attributes were collected, including the mean diameter at breast height (DBH), mean tree height, mean age, crown density, volume, land use and cover, and ecological conditions. The biomass, mean DBH, mean tree height, mean age and crown density were used as the inventory variables to establish the model, as shown in Table 1. The real biomass value was calculated using the equation for the biomass-volume relationship of the stand type and age group [30]. The stem volume of each tree was provided by FID, and the stand volume of each fixed sample plot was the sum of all tree volumes. The area of each plot was 0.0667 hectares.


**Table 1.** Statistics of the main forest inventory dataset (plot number and biomass of the three forest types in Beijing, China).

There are significant differences in topography and administrative functions among different districts in Beijing. From the topographic point of view, the northwest is a mountainous area with a higher terrain, and the southeast is a plain; from the administrative function point of view, the central part of Beijing is the capital and core functional area, the Northwest Mountainous Area and the southwest are ecological conservation functional areas, and the plain is a densely populated scientific and technological innovation and economic development area, which also leads to differences in forest biomass distribution. More than 80% of Beijing's forest resources are distributed in mountainous counties in the west and north of the city. The forest coverage in mountainous areas of Beijing has reached more than 50%, but the forest area in the plain of southeast Beijing accounts for less than 20% of the whole city.

**Figure 1.** Spatial distribution map of the forest sample plots in Beijing.

The forest ecosystem contains arbors, shrubs and herbs, but the amount of biomass from shrubs and herbs is less than the amount of arboreal biomass [31]. Therefore, this study considers only arboreal biomass and does not consider shrub and herb biomass.

#### 2.1.2. Remote Sensing Data and Preprocessing

We used optical (Landsat 8 OLI) and radar (ALOS-2 PALSAR-2) remote sensing sources.

#### Landsat 8

Landsat 8 Operational Land Imager (OLI, which developed by Bauer Aerospace and Technology Corp, Colorado, USA) images included a 15 m panchromatic band, with a spectral range from 0.500 to 0.680 μm and eight 30 m multispectral bands, with a spectral range from 0.433 to 2.300 μm. They were selected for biomass estimation due to their suitability in terms of their resolution ratio; a spatial resolution of approximately 30 m by 30 m is adequate to assess information at the forest stand level. We selected eight Landsat 8 OLI scenes of Beijing with low cloud cover as the research images. The image range was 122 to 124 paths and 31 to 33 rows. The image acquisition time used in this study was June-August 2016, and the time phase was basically the same as the time phase of the Beijing forest inventory data.

The image preprocessing steps included geometric correction, radiation correction, atmospheric correction, and image clipping. Because the downloaded images were Level-1 data products, the geometric accuracy was high, so only radiation and atmospheric corrections were needed.

Based on the preprocessed Landsat 8 OLI data, we acquired the surface reflectance for 6 bands of the Landsat 8 OLI (Band 2-Band 7) and then acquired the vegetation indices, namely, normalized difference vegetation index (NDVI), difference vegetation index (DVI) and ratio vegetation index (RVI) (Table 2), through band processing and the Landsat 8 OLI image calculation. The NDVI, DVI, and RVI are commonly used vegetation indices that are sensitive to vegetation, as shown in Equations (1)–(3). Another dataset was derived by image transformation from the original satellite band, which involved tasseled cap transformation (TCT) and texture features as shown in Table 2. Tasseled cap transformation, also known as a K-T transform, is an image enhancement method for vegetation information extraction. It can enhance vegetation information of images. After the K-T transform, the same number of components as the number of bands can be obtained, and the second

component is the green index, which has a strong relationship with the vegetation coverage and biomass on the ground [32]. Therefore, based on the TCT coefficients of the OLI sensor onboard Landsat 8, we chose the second band generated from the TCT, which was marked as the greenness [32].

**Table 2.** Remote sensing factors calculated from the Landsat 8 and ALOS-2/PALSAR-2 images.


In addition, we extracted the texture factor of the image. Texture is an important feature of remote sensing images and can be extracted by using the gray-level cooccurrence matrix (GLCM). Previous research has shown that Band 2 of a Landsat 8 OLI image contains much information about the image; thus, we extracted the texture feature of Band 2. The larger the selected window is, the greater the information content will be [33]. According to the sample area, five texture eigenvalues were extracted from the 15 × 15 window, namely, the mean, variance, contrast, correlation and second moment [34], as shown in Equations (4)–(8):

Normalized Difference Vegetation Index (NDVI):

$$\text{NIDVI} = \frac{\text{NIR1} - R}{\text{NIR1} + R} \tag{1}$$

Difference Vegetation Index (DVI):

$$\text{RVI} = \frac{\text{NIR1}}{\text{R}} \tag{2}$$

Ratio Vegetation Index (RVI):

$$\text{DVI} = \text{NIR1} - \text{R} \tag{3}$$

Mean(ME):

$$\text{ME} = \sum\_{i,j=0}^{N-1} iP\_{ij} \tag{4}$$

Variance(VA):

$$\text{VA} = \sum\_{i,j=0}^{N-1} P\_{ij}(i - \text{ME})^2 \tag{5}$$

Contrast(CO):

$$\text{CO} = \sum\_{i,j=0}^{N-1} iP\_{ij}(i-j)^2 \tag{6}$$

Correlation (CC):

$$\text{CC} = \sum\_{i,j=0}^{N-1} iP\_{ij} \left| \frac{(i - \text{ME})(j - \text{ME})}{\sqrt{VA\_iVA\_j}} \right| \tag{7}$$

Second Moment(SM):

$$\text{SM} = \sum\_{i,j=0}^{N-1} iP\_{ij}^{-2} \tag{8}$$

Using bilinear interpolation, the average values of the remote sensing factors at and near sampling points can be extracted. This method effectively solves the problem that occurs when sampling sites do not match the image completely and can cover areas that the inventory data cannot fully cover, thus improving the estimation accuracy.

#### ALOS-2/PALSAR-2

We downloaded six images of ALOS-2/PALSAR-2 (L-band) taken in 2016 from Japan Aerospace Exploration Agency (JAXA (http://www.eorc.jaxa.jp/ALOS/en/index.htm)). The PALSAR data had a 25-m spatial resolution and contain two polarized bands, HH and HV. The preprocessing of PALSAR data was completed by the JAXA. The digital numbers (DN) of the PALSAR signal amplitude were extracted and converted to gamma naught backscattering coefficients (dB) in decimal units using the following equation [35,36]:

$$I^0 = 10 \times \log\_{10}{DN^2} - \text{CF} \tag{9}$$

where Γ<sup>0</sup> is the backscattering coefficient, DN is the digital number value of pixels, and CF is the calibration factor, which equals −83 [36]. Then, we calculated the sum, difference and ratio values using the backscattering coefficients of HH and HV, as shown in Equations (10)–(12):

$$\mathbf{I}\cdot\mathbf{sum} = I\_{HH}^{0} + I\_{Hv}^{0} \tag{10}$$

$$\text{difference} = I\_{HH}^0 - I\_{Hv}^0 \tag{11}$$

$$\mathbf{r}\mathbf{r}\mathbf{i}\mathbf{i}=\boldsymbol{\Gamma}\_{HH}^{0}/\boldsymbol{\Gamma}\_{Hv'}^{0}\tag{12}$$

where Γ<sup>0</sup> *HH* and <sup>Γ</sup><sup>0</sup> *Hv* are the backscattering coefficients of HH and HV in decibels.

#### *2.2. Multiple Regression Model*

The allometric growth equation is the most widely used model for estimating forest biomass. Many studies have confirmed the advantages of the allometric growth equation for estimating forest biomass [37–39]. This model regresses a correlated variable (biomass) based on one or more independent variables. The DBH and tree height, as the two most relevant factors of biomass, are often used in biomass prediction in the form of single or compound variables. Based on the allometric model and previous research results, we introduced new variables (Landsat 8 data and backscattering coefficients) into the model to explore its ability to estimate forest biomass, as shown in Equation (13).

$$\ln(\mathbf{B}) = \beta\_0 + a \ln(d^2 H) + \beta\_1 \mathbf{x}\_1 + \beta\_2 \mathbf{x}\_2 + \dots + \beta\_j \mathbf{x}\_j \tag{13}$$

where B is the biomass of the sample plot, each *xj* is an independent variable (*j* = 1, 2, 3 ...), β*<sup>j</sup>* is the regression coefficient of *xj*, β<sup>0</sup> is a constant, and a is the regression coefficient of the model.

However, the remote sensing variables were highly collinear. To overcome this problem, we used a stepwise regression analysis method, which gradually screens variables and leaves highly correlated variables that are not collinear in the model, to retain a model that was not very complex and to reduce the number of calculations. The basic idea of stepwise regression is to introduce variables into the model one by one. After introducing an explanatory variable, we need to conduct F-test and t-test for the selected explanatory variables one by one. When the original explanatory variables are no longer significant due to the introduction of later explanatory variables, they will be deleted. To ensure that only significant variables are included in the regression equation before each new variable is introduced. This is a repeated process until neither significant explanatory variables are selected into the regression equation nor insignificant explanatory variables are removed from the regression equation. To ensure that the final set of explanatory variables is optimal.

#### *2.3. Feature Selection and Random Forest Model*

We used R to establish an RF model to estimate forest biomass. The RF model was a classification and regression algorithm based on decision trees [40]. By establishing and combining multiple decision tree predictions (1000 trees in our study), the average value of all the decision tree prediction results was taken as the final prediction result [41]. The RF model can effectively alleviate the problem of overfitting and is insensitive to the collinearity between variables, so it is suitable for establishing a nonlinear model [42]. RF is increasingly used to perform biomass regression and estimate forest biomass [43,44].

First, subsets of variables were selected as input for the RF prediction using feature selection to ensure that the input variables were highly correlated with biomass. Feature selection refers to the selection of subsets from the original feature set to optimize a certain evaluation criterion so that the model established with the optimal feature subset can achieve a prediction accuracy similar to or better than that of the model established without feature selection. RF provides an increase in the mean-squared error (percentage of IncMSE, where IncMSE indicates the increase in MSE) for each independent indicator, quantifying the increase in the MSE when the indicator is randomly permuted. This error measures the relative importance of each indicator, where a high IncMSE implies that the indicator has a high weight in the model prediction and vice versa [23]. Then, we used the data after feature selection as the independent variable, forest biomass as the dependent variable, and the random forest software package in R to establish an RF model.

#### *2.4. Model Accuracy Evaluation*

To test these models, we assessed the prediction accuracy on randomly selected subsets (20%) of the original dataset retained before the model was developed. To evaluate the advantage of the use of an advanced regression tree model versus more traditional approaches, the performance of the RF model was computed and compared with that of a stepwise multiple linear regression model.

We used the proportion of variance explained (*R*2) and the root mean square error (RMSE) to evaluate the model performance on the complete datasets. In addition, we computed the relative RMSE (RMSE%), the bias and the relative bias (bias%). Bias was calculated as the difference between a population mean of the measurements or test results and an accepted reference or true value, R<sup>2</sup> values were used to judge the model, and RMSE, Bias%, RMSE% reflect the precision of the model [45].

These statistics were calculated as follows:

$$R^2 = 1 - \frac{\sum \left( y\_i - \hat{y}\_i \right)^2}{\sum \left( y\_i - \overline{y}\_i \right)^2} \tag{14}$$

$$\text{RMSE} = \sqrt{\frac{\sum \left( y\_i - \hat{y}\_i \right)^2}{n - 1}} \tag{15}$$

$$\text{RMSE\%} = \frac{\text{RMSE}}{\overline{y}\_i} \times 100\% \tag{16}$$

$$\text{Biasc} = \frac{1}{n} \sum\_{i=1}^{n} (y\_i - \hat{y}\_i) \tag{17}$$

$$\text{Bias\%} = \frac{\text{BIAS}}{\overline{y}\_i} \times 100\% \tag{18}$$

where *yi* is the observed biomass of the plot, *y*ˆ*<sup>i</sup>* is the predicted biomass of the plot, and *yi* is the mean biomass of n plots.

#### **3. Results**

#### *3.1. Univariate Correlation Analysis*

Previous studies have typically analyzed the relationship between a single remote sensing variable and the forest biomass or have used the original band and variables transformed from images for feature selection [46,47]. Through the Pearson correlation coefficient (r), we analyzed the ability of each variable to estimate biomass and obtained the correlation between each variable and the biomass, as shown in Table 3.


**Table 3.** Coefficients of correlation between forest biomass and variables.

Among all variables, forest inventory variables were highly correlated with biomass in three forest types. Different remote sensing variables (OLI data and PALSAR data) showed different degrees of correlation. The shortwave infrared (SWIR) optical band (Band 7) showed the greatest relevant biomass among the Landsat data because it allowed an effective separation between high- and low-biomass data. The importance of the SWIR wavelengths in biomass prediction is consistent with previous studies [40]. In addition, Band 5, Band 4, Band 3, Band 2, the NDVI, and the greenness were also highly correlated with biomass, and the most relevant texture factors were the mean, correlation and second moment. The other Landsat variables had little correlation with biomass. SAR data can penetrate dense forests and obtain the vertical structure information of forests, so the PALSAR HH and HV backscatter coefficients and their derivative variables (sum, difference, ratio) were correlated with forest biomass. In addition, correlations between the forest biomass and HV backscatter coefficients of different forest types were higher than those between the forest biomass and HH backscatter coefficients, which is in line with previous research results [22,48]. All these factors can be considered potential variables for forest biomass estimation.

#### *3.2. Results of Forest Biomass Model Establishment*

#### 3.2.1. Multiple Stepwise Regression Model

To avoid overfitting, the multiple stepwise regression method was used to screen variables and establish a multiple linear model. The results are as follows:

The multiple stepwise regression model of coniferous forests was:

$$\ln(B\_{\mathbb{C}}) = 3.821 + 0.226 \times \text{N}\_{\mathbb{T}} + 0.111 \text{N}\_{\mathbb{Z}} + 0.199 \text{N}\_{\mathbb{S}} - 0.120 \text{X}\_{\mathbb{Z}} - 0.246 \text{X}\_{\mathbb{B}} - 0.089 \text{X}\_{\mathbb{B}} - 0.656 \text{X}\_{\mathbb{B}} + 0.598 \text{X}\_{\mathbb{Z}} - 0.308 \text{X}\_{\mathbb{B}} \tag{19}$$

The multiple stepwise regression model of mixed forest was:

ln(*BM*) = 3.776 + 0.291 × *N*<sup>1</sup> + 0.155*N*<sup>2</sup> + 0.108*N*<sup>3</sup> − 0.153*X*<sup>1</sup> + 0.302*X*<sup>2</sup> − 0.132*X*<sup>3</sup> − 0.05*X*<sup>4</sup> + 0.052*X*<sup>6</sup> + 0.037*X*<sup>13</sup> (20)

The multiple stepwise regression model of broadleaf forests was:

$$\ln(B\_{\rm B}) = 3.810 + 0.127 \times N\_1 + 0.137 N\_2 + 0.181 N\_3 + 0.110 X\_1 + 0.099 X\_6 - 0.088 X\_9 + 0.024 X\_{13} + 0.043 X\_{15} - 0.057 X\_{18} \tag{21}$$

Then, the biomass estimation model was obtained as follows:

$$\begin{aligned} 3.821 + 0.226 \times N\_1 + 0.111 N\_2 + 0.139 N\_3 - 0.120 X\_1 \\ B\_C = e \end{aligned} $$
 
$$B\_C = e^{-0.246X\_6 - 0.089X\_{10} - 0.656X\_{16} + 0.538X\_{17} - 0.308X\_{19}} \tag{22}$$

$$\begin{array}{c} 3.776 + 0.291 \times \text{N}\_1 + 0.155 \text{N}\_2 + 0.108 \text{N}\_3 - 0.153 \text{X}\_1 \\ B\_M = e \\ B\_{M} = e \end{array} + 0.302 \text{X}\_2 - 0.132 \text{X}\_3 - 0.05 \text{X}\_4 + 0.052 \text{X}\_6 + 0.037 \text{X}\_{13} \tag{23}$$
 
$$\text{3.212 + 2.147} \quad \text{N}\_1 + 0.107 \text{Y}\_1 + 0.101 \text{Y}\_2 + 0.241 \text{Y}\_3 \tag{24}$$

$$\begin{aligned} 3.810 + 0.127 \times \text{N}\_1 + 0.137 \text{N}\_2 + 0.181 \text{N}\_3 + 0.110 \text{X}\_1 \\ B\_B = \varepsilon \end{aligned} \quad + 0.039 \text{X}\_6 - 0.088 \text{X}\_9 + 0.024 \text{X}\_{13} + 0.043 \text{X}\_{15} - 0.057 \text{X}\_{18} \tag{24}$$

where *Xi* and *Ni* in each formula correspond to the variables in Table 3.

P values represent the probability that the sample results differ from the original hypothesis. The smaller the P value is, the more significant the results are. Generally speaking, *p* < 0.05 indicates a significant difference, and *p* < 0.01 indicates a very significant difference. Table 4 shows that the P values of the model coefficients of different forest types, which are less than 0.05, some of which are less than 0.01. These results show that the differences in the selected variables are significant.

**Table 4.** P values of the coefficients of models for different forest types.


#### 3.2.2. Random Forest Model

First, feature variables were selected for the variables involved in the modeling. The importance of the variables was ranked according to IncMSE%, and the unimportant variables were eliminated. Generally, the number of final variables is 1/3 of the total number of input variables [49]. Table 5 shows that *d*2*H* and the mean age are two very important variables in the RF model. IncMSE% was more than 20% in the different forest types. NDVI, Band 2 and Band 7 were also important to the model with

regard to optical data. The most influential backscattering coefficient factors were Γ<sup>0</sup> *Hv* and <sup>Γ</sup><sup>0</sup> *HH* <sup>−</sup> <sup>Γ</sup><sup>0</sup> *Hv*. The IncMSE% of these factors were all higher than 10%.


**Table 5.** The IncMSE% of the top five most important variables in the biomass fitting of different forest types in the random forest model.

#### *3.3. Model Precision Evaluation and Comparison of Two Models*

To test the goodness of fit of the model, 20% of the samples were used for validation. We analyzed the scatter diagrams of different forest types in Figure 2 and obtained the accuracy of the linear regression and RF models in Table 6.

**Figure 2.** Graphs of the predicted versus observed values of three forest types for two models. (**a1**) and (**a2**) Coniferous forest; (**b1**) and (**b2**) Mixed forest; (**c1**) and (**c2**): Broadleaf forest.

Figure 2 shows that correlations between the estimated biomass and observed biomass in the coniferous forests, mixed forests, and broadleaf forests were all better for the RF models than for linear regression. In the linear model, as the biomass value increased, the performance of the model decreased, and most of the high-value biomass was underestimated. In particular, the high biomass values of mixed forests were greatly underestimated. RF can improve the performance of the model. When the biomass was less than 100 Mg ha<sup>−</sup>1, the difference between predicted and observed values is

lower than that of the linear model; the error under higher biomass values was slightly larger, and some of the higher values were underestimated.


**Table 6.** Estimation accuracy of different models.

Table 6 shows that the Bias% values were all near 0, and the RMSE% ranged from 27.92% to 32.48% for the MLR model. For the RF model, the Bias ranged from −2.67 to −4.02 and the RMSE% ranged from 22.89%–27.24%. These results showed that the two types of models were relatively stable and could be used to estimate biomass. However, an improvement in performance was found in the RF models for coniferous forest (*R*<sup>2</sup> = 0.66, RMSE = 13.23 Mg ha−1), mixed forest (*R*2= 0.77, RMSE = 11.09 Mg ha−1), and broadleaf forest (*R*<sup>2</sup> = 0.64, RMSE = 11.98 Mg ha−1), in comparison to the linear regression models for coniferous forest (*R*<sup>2</sup> = 0.59, RMSE = 14.15 Mg ha−1), mixed forest *R*<sup>2</sup> = 0.70, RMSE = 14.54 Mg ha<sup>−</sup>1), and broadleaf forest (*R*<sup>2</sup> = 0.53, RMSE = 15.26 Mg ha<sup>−</sup>1). Generally, the RF model was characterized by a high *R*<sup>2</sup> and a low RMSE, indicating a good fitting result.

For the same model with different forest types, the fit of the mixed-forest models was better than that of the models with the other two forest types. The *R*<sup>2</sup> of the linear regression model based on mixed forests was 0.11 and 0.17 higher than the *R*<sup>2</sup> of the linear regression models based on coniferous forests and broadleaf forests, respectively. The *R*<sup>2</sup> of the RF model based on mixed forests was 0.11 and 0.13 higher than the *R*<sup>2</sup> of the RF models based on coniferous forests and broadleaf forests, respectively, and the RMSE was 4.35 and 4.13 Mg ha−<sup>1</sup> lower, respectively. Overall, the model for mixed forests had a high estimation accuracy.

#### *3.4. Results of Biomass for Di*ff*erent Forest Types and Spatial Distribution of the Forest Biomass Density in Beijing*

Based on the model estimation results, the forest biomass and biomass density of coniferous forests, broadleaf forests and mixed forests were estimated, the kriging interpolation was used and biomass density distribution maps of three forest types were obtained. Biomass of different forest types are shown in Table 7 and biomass density distribution in Figure 3.


**Table 7.** Biomass and biomass density of each forest type.

The total forest biomass obtained from the survey data was 74,746.10Mg, and the biomass density was 19.14–195.66 Mg ha<sup>−</sup>1, with an average biomass density of 52.26 Mg ha−1. Among these values, the total biomass of coniferous forest was 35,622.99 Mg, the average biomass density was 53.73 Mg ha<sup>−</sup>1; the total biomass of mixed forest was 13,926.40 Mg, the average biomass density was 51.20 Mg ha<sup>−</sup>1; the total biomass of broadleaf forest was 25,196.80Mg, and the average biomass density was 50.80 Mg ha<sup>−</sup>1.

**Figure 3.** Biomass density distribution in Beijing, China (**a**) Coniferous forest; (**b**) Mixed forest; (**c**) Broadleaf forest; (**d**) all forest sampling plots.

ȱ ȱ

As shown in Figure 3, the biomass distribution of arbor forests was basically consistent with the distribution of forestland in Beijing. The high-biomass area corresponded to dense forestland, and most of these forests were mature and overmature and were mainly distributed in the north and southwestern part of Beijing. The low-biomass area was mainly located in the southeast and central parts of Beijing. Because this area is urban with mostly developed land, the biomass in this area is low. The biomass density in most areas of Beijing was less than 70 Mg ha<sup>−</sup>1.

The distribution of the three types of forest was obviously different. Coniferous forests were mainly distributed to the west and south of Beijing, mixed forests were mainly distributed in the west, and broadleaf forests were mainly distributed in the north and southwest.

This study did not consider shrubs and herbs, so the estimation of biomass can be considered relatively conservative but can also represent the basic situation of biomass in Beijing. At present, China's biomass estimation system is still not perfect. This study provides a feasible method for regional biomass estimation.

#### **4. Discussion**

#### *4.1. Forest Biomass Estimation Model Based on Forest Inventory and Multisource Remote Sensing Data*

In this paper, we propose a novel approach to modeling and mapping the biomass of forests at the regional scale that provides more detailed and accurate information than other approaches, such as estimating using only a single remote sensing data source or forest inventory data.

We combined forest inventory data with multisource remote sensing data (OLI and PALSAR) to estimate forest biomass, capturing almost all forest biomass spatial variability, and producing spatially explicit biomass estimates over regions.

According to the biomass characteristics of different forest types, it is very important to select variables with a high importance to the model [48]. The forest inventory factors selected in this study included not only the DBH and height, which are the two most relevant factors to biomass [38,39,50] but also the mean age and canopy density, which have received increasing attention in recent studies. Many previous studies have demonstrated that these two factors show a good correlation with biomass [40,51], which is consistent with our results. Landsat optical data are sensitive to forest

vegetation, and their spatial resolution is suitable for the sample plot size. Zheng et al. confirmed that the red and NIR bands (Bands 4 and 5, respectively) are effective predictors of biomass [13]. Zheng et al. found that the SWIR band (Band 7) showed a satisfactory estimation ability in forests with a high canopy density [52], Foody et al. found that the NDVI and other vegetation indices are strongly correlated with biomass [53]. These results are consistent with our results. The models established by these factors fit the data relatively well. Multiple-variable PALSAR data have a higher correlation with forest biomass than individual-variable PALSAR HH and HV data because of their ability to detect canopy structure and retrieve forest biomass [54,55]. In addition, the biomass estimation model based on multisource remote sensing data combined with forest inventory data had a higher accuracy than that based on single-source data. Zhao used Landsat TM and ALOS PALSAR data to establish forest biomass models in Zhejiang Province. The *R*<sup>2</sup> of each forest type was below 0.5 [56]. Urbazaev used SAR backscatter, Landsat images and topographic factors to obtain the best *R*2, a value of only 0.62 [57], which was lower than the accuracy in this study. This indicates that our model is suitable for estimating forest biomass in Beijing.

Optical and radar data are an effective supplement to inventory data, provide spatial information for estimating regional forest biomass, and can continuously estimate forest biomass [51]. In the results, the *R*<sup>2</sup> and RMSE of the three forest types were all greater than 0.5 and less than 20 Mg ha<sup>−</sup>1. The models are reliable, but the model accuracy differed among different forest types.

In our results, mixed forests had the highest estimation accuracy, followed by coniferous forests and broadleaf forests, which is inconsistent with previous research.

Previous studies have shown that coniferous forest biomass estimation models have a high accuracy [58]. This inconsistency may be because the three types of modeling factors included in this study are highly sensitive to the structure of mixed forests. Moreover, differences in the study area location, tree species and forest types lead to different model estimation accuracies for different forest types.

However, the *R*<sup>2</sup> values of the model were all less than 0.8, which indicates the present model has less precision than the model established in a previous study [25]. A possible reason is that previous studies have mostly focused on small-scale areas. Our research mainly focuses on meso-scale areas, including a variety of terrain and environmental conditions, causing different environmental factors to have a certain impact on the modeling accuracy.

#### *4.2. Estimation and Spatial Distribution of Forest Biomass in Beijing*

The results of this study suggest that it is possible to produce spatially explicit biomass estimates over regions if adequate inventory data and remote sensing data are available. This meso-scale study was based on a relatively large sample size.

Because the data set used in this research did not contain all sample plot data from the forest inventory in Beijing, it can represent the distribution but not the total amount of biomass, while biomass density can represent the state of forest resources in Beijing well [9]. The 2016 forest biomass density range of Beijing was estimated by this model to be 19.14–195.66 Mg ha−1, and the average biomass density was 52.26 Mg ha<sup>−</sup>1. The forest biomass density in Beijing increased compared with previous studies [24].

Overall, forest resources show a pattern of more forests in mountainous areas, less forests in plains, as well as more forests outside urban areas and less forests in urban areas [59]. This is consistent with the distribution map of biomass density obtained in this study, as shown in Figure 3. However, this study shows that the biomass density in the central city is high, which is inconsistent with the distribution of forest resources in Beijing. This result may be caused by the fact that in recent years, the central city has insisted on the construction of ecological cities and increases in the area of green space, so the forest biomass has increased.

China is vigorously promoting the construction of eco-cities. At the meso-scale, forestry biomass estimation and biomass density mapping can provide decision makers with detailed information on forest resources to strengthen the management of forest resources.

#### *4.3. Comparing Performance between MLR and RF Models*

We also compared the performance of the MLR model and the RF model. We found that the prediction effect of the linear model for extreme biomass values (extremely high or very low) was not completely ideal. Low values were overestimated, and high values were underestimated, which is also common in previous studies [25,60]. The range in biomass for the multiple regression prediction is larger than that for the RF model.

This finding indicates that this model may be applicable only in the estimation of biomass values within a certain range. The use of the RF model had a positive impact on the estimation accuracy of extreme values [23]. The results showed that the RF model captured the complex nonlinear relationship between the optical and SAR data and biomass and compensated for the lack of inventory data, capitalizing on the strengths of both the forest inventory and remote sensing data. Therefore, the fit of the RF model is better than that of the linear regression model.

Linear regression and RF model contain different independent variable factors. Because of the collinearity between variables, some variables will be eliminated in linear model, while RF model can fully consider the fitting problem. The two model types include forest investigation factors such as *d*2*H*, crown density and the mean age, indicating that the FID variables has a greater impact on the estimation of biomass and is not affected by collinearity. The RF models contain more PALSAR variables, and the model accuracy is higher, which indicates that the PALSAR variables are more sensitive to forest biomass.

Certainly, there are still some limitations in our research. Two thirds of Beijing is plain, and one third is mountainous area. Terrain correction is a part of remote sensing image correction under rugged surface, which can offset the influence of terrain to a certain extent, and is helpful to improve the accuracy of biomass estimation. Therefore, the research of terrain correction will be strengthened in future research. In the image preprocessing stage, the inaccuracy of the biomass models based on forest types and age classes and the lack of a consideration for the impact of environmental factors such as topography, soil and hydrology on biomass, which will be strengthened in the future research. Despite these problems, this study aimed to improve the performance of the regional forest biomass model and can provide a reference and support for future plans of relevant forestry departments, which has certain practical significance.

#### **5. Conclusions**

This paper proposed an approach for establishing the forest biomass of different forest type models and calculating forest biomass in Beijing by combining forest inventory data with multisource remote sensing data. The approach can capture all spatial variability and provide a reliable method for estimating forest biomass at a meso-scale with a high efficiency and low cost. In addition, we used this model to predict the forest biomass in Beijing in 2016. Among the three studied forest types, coniferous forest had the highest biomass density. According to the distribution of forest biomass in Beijing, the northern and southwestern parts of Beijing had a high biomass, while the central and eastern parts have a low biomass density. At present, there is no perfect biomass estimation system in China. Therefore, this method can provide a basis for meso-level biomass estimation and a reference for the planning of relevant forestry decision-making departments.

**Author Contributions:** Conceptualization, Z.F. and Y.Z.; Methodology, Z.F.; Validation, Z.F. and Y.Z.; Formal Analysis, Y.Z.; Investigation, Z.F. and Y.Z.; Resources, Z.F.; Data Curation, Z.F. and J.L. (Jincheng Liu); Writing-Original Draft Preparation, Y.Z.; Writing-Review & Editing, Y.Z. and J.L. (Jing Lu); Visualization, Y.Z.; Supervision, Z.F.; Project Administration, Z.F.; and Funding Acquisition, Z.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was jointly supported by the medium long-term project of "Precision Forestry Key Technology and Equipment Research" (No. 2015ZCQ-LX-01) and the National Natural Science Foundation of China (No.U1710123).

**Acknowledgments:** We would like to acknowledge support from Beijing Key Laboratory for Precision Forestry, Beijing Forestry University, as well as all the people who have contributed to this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
