Deep Learning Based Retrieval of Forest Aboveground Biomass from Combined LiDAR and Landsat 8 Data

Estimation of forest aboveground biomass (AGB) is crucial for various technical and scientific applications, ranging from regional carbon and bioenergy policies to sustainable forest management. However, passive optical remote sensing, which is the most widely used remote sensing data for retrieving vegetation parameters, is constrained by spectral saturation problems and cloud cover. On the other hand, LiDAR data, which have been extensively used to estimate forest structure attributes, cannot provide sufficient spectral information of vegetation canopies. Thus, this study aimed to develop a novel synergistic approach to estimating biomass by integrating LiDAR data with Landsat 8 imagery through a deep learning-based workflow. First the relationships between biomass and spectral vegetation indices (SVIs) and LiDAR metrics were separately investigated. Next, two groups of combined optical and LiDAR indices (i.e., COLI1 and COLI2) were designed and explored to identify their performances in biomass estimation. Finally, five prediction models, including K-nearest Neighbor, Random Forest, Support Vector Regression, the deep learning model, i.e., Stacked Sparse Autoencoder network (SSAE), and multiple stepwise linear regressions, were individually used to estimate biomass with input variables of different scenarios, i.e., (i) all the COLI1 (ACOLI1), (ii) all the COLI2 (ACOLI2), (iii) ACOLI1 and all the optical (AO) and LiDAR variables (AL), and (iv) ACOLI2, AO and AL. Results showed that univariate models with the combined optical and LiDAR indices as explanatory variables presented better modeling performance than those with either optical or LiDAR data alone, regardless of the combination mode. The SSAE model obtained the best performance compared to the other tested prediction algorithms for the forest biomass estimation. The best predictive accuracy was achieved by the SSAE model with inputs of combined optical and LiDAR variables (i.e., ACOLI1, AO and AL) that yielded an R2 of 0.935, root mean squared error (RMSE) of 15.67 Mg/ha, and relative root mean squared error (RMSEr) of 11.407%. It was concluded that the presented combined indices were simple and effective by integrating LiDAR-derived structure information with Landsat 8 spectral data for estimating forest biomass. Overall, the SSAE model with inputs of Landsat 8 and LiDAR integrated information resulted in accurate estimation of forest biomass. The presented modeling workflow will greatly facilitate future forest biomass estimation and carbon stock assessments.


Introduction
Reliable, up-to-date forest aboveground biomass (AGB) mapping is a prerequisite for understanding the relationship between AGB and climate change. Plot-based estimations of forest AGB, while typically of high accuracy, are costly and can only provide quality information for a limited number of stands at the landscape scale [1]. Consequently, numerous remote sensing techniques have been increasingly utilized to assist forest AGB estimation during the last few decades [2]. It is possible to provide updated, consistent, and spatially explicit assessment of forest biomass and its dynamics by using remote sensing images, particularly in large areas with limited accessibility [3,4].
Numerous statistical models have been explored in relating field-measured AGB to remotely sensed variables [5,6]. Multispectral vegetation indices (SVIs), defined with various combinations of visible, near-infrared (NIR) and shortwave reflectance, are the most widely used for retrieving biophysical and biochemical parameters. It is proved that SVIs have strong correlations with vegetation structure characteristics, such as AGB and leaf area index (LAI) [7][8][9]. However, vegetation indices tend to saturate for forests at high biomass levels [10][11][12]. The saturation point varies greatly depending on the source data and the vegetation type and ranges from 15 to 100 Mg/ha for different visible/NIR vegetation indices [13,14]. Additionally, optical remote sensing provides limited information on the vertical distribution of forest structure [15], and it is not always possible to compile a temporally and radiometrically consistent cloud-free datasets over large areas [3]. The past two decades have witnessed a large number of studies using Synthetic Aperture Radar (SAR), which could penetrate clouds and forest canopies with appropriate wavelength and polarization modes, for mapping forest AGB [16][17][18][19]. However, it remains problematic to estimate AGB using SAR backscattering signals due not only to the saturation at high biomass levels but also to the high sensitivity to soil conditions, including surface roughness and soil moisture [19,20].
Recently, light detection and ranging (LiDAR) has proven to be very powerful in estimating forest structure attributes, such as canopy height, leaf area index (LAI), and AGB [21][22][23][24]. The most commonly used LiDAR metrics are height metrics which can be directly measured by LiDAR and provide information related to the vertical structure of individual trees and forest stands. The prevalence of LiDAR data in forest studies is attributed to its penetration by which layered structural echoes through a certain canopy depth can be detected. However, LiDAR data cannot provide sufficient spectral characteristics of vegetation canopies since most LiDAR systems only work at a single wavelength [25]. Although hyperspectral LiDAR systems have emerged to capture spectral and structural information simultaneously [26][27][28], they are mainly tested inside the laboratory. It is, thus, not feasible to apply these new systems to normal data collection in large areas [26]. Therefore, the integration of LiDAR data with optical remote sensing imagery has been identified as the most promising approach to acquire structural and spectral information from forests simultaneously for biomass estimation. Combining both airborne LiDAR and hyperspectral data has shown great capability to map tree species in different forest areas [29,30]. In addition, Graham [31] found that retrieval of canopy LAI in a coniferous and broadleaf mixed forest can be improved by integrating LiDAR with WorldView-2 data.
The parametric models, such as multiple regression, are commonly used to develop relationships between forest attributes and remote sensing predictors [32,33]. In recent years, non-parametric machine learning models have become prevalent. Contrary to the linear regression model, many machine learning techniques (e.g., Random Forest (RF), Support Vector Regression (SVR), K-nearest Neighbor (KNN), and Deep learning (DL)) are able to reveal complicated non-linear patterns [9,34,35]. Additionally, machine learning models are able to address issues associated with data dimensionality [36,37] in fitting models with a large number of predictors. Few studies have integrated LiDAR data with optical remote sensing metrics to estimate forest structure parameters using non-parametric machine learning algorithms, especially the deep learning models. The deep learning algorithms have shown their effectiveness in object detection and image classification [38][39][40][41][42][43]. The deep learning models can automatically extracted invariant and abstract features which have better discrimination than artificial features. As one of the deep learning algorithms, the Stacked Sparse Autoencoder network Remote Sens. 2019, 11, 1459 3 of 17 (SSAE) has been widely used in some fields, especially image classification [43]. Besides, few studies paid attention to the comparisons between deep learning and other machine learning algorithms in predicting forest parameters.
Therefore, in this study, airborne LiDAR data were integrated with Landsat 8 imagery for mapping forest aboveground biomass. The objectives of the study were (i) to evaluate the potentials of optical spectral vegetation indices (SVIs) and LiDAR metrics, respectively, for estimating forest AGB; (ii) to determine if the combined LiDAR and Landsat 8 indices can enhance the AGB estimation; and (iii) to provide a comparison of five predictive models, i.e., K-nearest Neighbor (KNN), Random Forest (RF), Support Vector Regression (SVR), Stacked Sparse Autoencoder network (SSAE), and multiple stepwise linear regressions (MSLR) models, for estimating biomass.

Study Area
The study site is located in the northeast part of Conghua (23 •  algorithms, the Stacked Sparse Autoencoder network (SSAE) has been widely used in some fields, especially image classification [43]. Besides, few studies paid attention to the comparisons between deep learning and other machine learning algorithms in predicting forest parameters. Therefore, in this study, airborne LiDAR data were integrated with Landsat 8 imagery for mapping forest aboveground biomass. The objectives of the study were (i) to evaluate the potentials of optical spectral vegetation indices (SVIs) and LiDAR metrics, respectively, for estimating forest AGB; (ii) to determine if the combined LiDAR and Landsat 8 indices can enhance the AGB estimation; and (iii) to provide a comparison of five predictive models, i.e., K-nearest Neighbor (KNN), Random Forest (RF), Support Vector Regression (SVR), Stacked Sparse Autoencoder network (SSAE), and multiple stepwise linear regressions (MSLR) models, for estimating biomass.

Field Data
Aboveground biomass values of 236 inventory plots distributed across the study area were collected in 2013. The locations of plots were selected according to the subjective sampling evaluations, and the size of each plot was 30 m × 30 m. The tree height (H) was measured using a laser hypsometer, and the tree diameter at breast height (DBH) was measured using a tape. The center of each plot was correctly determined using a GPS (Garmin MAP 60CS, accuracy ±3 m). Additionally, the species and type (i.e., evergreen or deciduous) of each tree were recorded. Trees with DBH more than or equal to 3 cm were measured in the survey. Subsequently, a plot would be classified into a coniferous (deciduous) plot, when it consisted of over 75% coniferous (deciduous) trees [44]. Otherwise, we identified the plot as the mixed forest [45]. The observed plots (n = 236) were randomly split into training (n = 177) and validation (n = 59) datasets at a ratio of 3:1.
In general, AGB of each individual tree in the plot can be derived by using species-specific allometric equations with inputs of the DBH and H [20,46,47]. However, no allometric equations are available for the study area. Accordingly, the methodology put forward by Fang et al. [48] was employed to calculate the biomass at the plot level. In this method, the relationship between total volumes and total AGB of each plot is shown in Equation (1). The total volumes could be obtained by adding up the volumes of all individual trees in the plot. The single tree volume was estimated based on a volume table with inputs of H and DBH.
where TAGB and TV are separately the total AGB and total volumes in a plot, and the plot AGB needs to be finally computed using TAGB at a megagrams per hectare conversion unit. a and b are coefficients from [48] and they are different for diverse forest types. The biomass values were in the range from 17.464 to 313.918 Mg/ha with a mean value of 125.745 Mg/ha and a standard deviation of 71.13 Mg/ha.

LiDAR Data Acquisition and Preprocessing
Airborne observations for the study area were conducted in June 2013 using a Leica ALS50-II laser scanning system. The system recorded both first and last return data for each laser pulse. The pulse frequency was 52.9 kHz and the flying speed was 80 knots, producing an average density of 0.8 pulses/m 2 . The LiDAR data was preprocessed using the Terrascan software (v4.006-Terrasolid, Helsinki, Finland). Firstly, the points were removed if they had only few neighbors and/or their elevations were higher than the median elevation of surrounding points. The points were then classified as ground and non-ground returns. The ground returns were interpolated to produce a digital elevation model (DEM), and the first returns were interpolated to derive a digital surface model (DSM) with a resolution of 1 m. Finally, a Canopy Height Model (CHM) was generated by subtracting the DEM from the DSM. According to the minimum and maximum height of field-measured trees within the area, the CHM pixels with values ranging from 2 m to 35 m were extracted to ensure the understory vegetation and objects exceeding the tree height were excluded [21,49,50]. Figure 2 shows the workflow of the biomass estimation processes, including the variables calculation, model calibrations using different data scenarios, accuracy assessment, and wall-to-wall biomass mapping using the calibrated model. subtracting the DEM from the DSM. According to the minimum and maximum height of field-measured trees within the area, the CHM pixels with values ranging from 2 m to 35 m were extracted to ensure the understory vegetation and objects exceeding the tree height were excluded [21,49,50]. Figure 2 shows the workflow of the biomass estimation processes, including the variables calculation, model calibrations using different data scenarios, accuracy assessment, and wall-to-wall biomass mapping using the calibrated model.

Landsat 8 OLI and LiDAR Variables
According to previous studies [7,9,51], seven SVIs were selected, including the Normalized Difference Vegetation Index (NDVI), Simple Ratio Vegetation Index (SR), Enhanced Vegetation Index  Table 1). A total of 16 LiDAR variables were computed at the plot level. The canopy cover (Cov) variable was derived using the proportion of the number of pulses returned from the canopy to all the returns. Six classes of canopy height metrics calculating the height distribution and height summary statistics of the canopy were extracted, including the mean (H mean ), maximum (H max ), standard deviation (H std ), variance (H var ), coefficient of variation (H cv ), and percentiles with interval of 10% (H p : p 10, p 20 , . . . , p 90 ). In addition, canopy relief ratio (CRR) which is a quantitative index of the canopy relative shape, describing the ratio of all returns above the mean heights was extimated (Equation (2)) [52].

Combined Optical and LiDAR Index (COLI)
As described in the Section 1, optical data can provide spectral information about vegetation canopy, while LiDAR data measures the vertical structure of the forest. Thus, the synergistic utilization of airborne LiDAR and optical remote sensing data is highly valuable for estimating forest biomass. Inspired by the SVIs of optical remote sensing, two types of new indices incorporating optical and LiDAR information (COLI1 and COLI2) were established by integrating the best-performing LiDAR variable (i.e., the LiDAR variable achieving the best AGB estimation accuracy) with each SVI. COLI1 and COLI2 can be written as: where BLV is the best-performing LiDAR variable and SVI i is the optical spectral vegetation index. Prior to calculating COLI, the BLV and SVI i were normalized to allow direct comparisons between optical and LiDAR variables, e.g., with different scales and dynamic ranges.

Regression Algorithms
In the study, five prediction techniques including KNN, RF, SVR, SSAE, and MSLR were used to estimate biomass. The RF technique is based on an ensemble of binary regression trees that are fitted to randomly selected subsets of the training data. In addition to resampling the observations to obtain multiple trees, the RF algorithm also selects a random subset of predictors in tree construction, which is particularly useful when a great number of possibly redundant predictors are available. It has been reported to be an efficient prediction approach for estimating vegetation attributes, especially when the number of predictors is very large [35,53,54]. In [35], the RF was compared with other prediction methods in regard to their predictive power of forest AGB, and RF showed a better estimation accuracy than other algorithms. SVR identifies optimum hyperplanes using kernel functions to separate groups of input data with similar responses to predict a target variable. The SVR model has been successfully applied in biomass mapping and other remote sensing applications [55,56]. In KNN method, the distance from an unknown pixel of the target dataset to each of the known pixels in the reference dataset is computed. The unknown target pixel is assigned a weighted mean of the k most similar neighbours, which is computed based on their distance to the target pixel in the feature space [57]. Tian et al. [34] found the performance of KNN outperformed those of other algorithms for estimating forest biomass using Landsat Thematic Mapper (TM) data. SSAE is a kind of deep learning method that can automatically learn useful features layer by layer. The SSAE deep learning algorithm has been successfully used in several fields, such as image classification [39,43]. Nonetheless, as far as we know, this model was rarely applied to estimate forest parameters. Therefore, the SSAE algorithm was introduced in detail in the following. When SSAE model builds a deep neural network to extract deep features, it utilizes a hierarchical training strategy. The model comprises some sparse autoencoder networks (SAEs), and each SAE consists of two parts: an encoder and a decoder. After each SAE is trained, its decoded layer is removed and then an SSAE is established using the encoder parameters of all SAEs in a layer-by-layer manner. Afterwards, this network is linked to a regression model to implement the prediction work. That is to say, an SSAE model is comprised of the trained SAEs and the regression model. The MSLR was taken into consideration in order to compare the prediction performances of the more prevalent and complex machine learning models with that of the frequently-used parametric model.
We employed the R statistical package [58] to implement RF, SVR, KNN, and MSLR models. The Deep Learning toolbox, an open source library (https://github.com/rasmusbergpalm/ DeepLearnToolbox#deeplearntoolbox), was modified to perform SSAE algorithm. Table 2 showed the three experiments (Experiments 1-3) designed with different data scenarios. Experiment 1: univariate linear or nonlinear regression models were developed based on a single SVI or LiDAR variable; Experiment 2: each COLI was used to establish the univariate linear or nonlinear regression model; Experiment 3: four data combinations were individually fitted to all the five prediction models described in Section 3.2. The four data combinations included (i) all the COLI1 (ACOLI1), (ii) all the COLI2 (ACOLI2), (iii) ACOLI1 and all the optical (AO) and LiDAR variables (AL), and (iv) ACOLI2, AO and AL. Individual model predictions were validated against the same independent dataset (n = 59). Three indicators were calculated, including coefficient of determination (R 2 ), root mean squared error (RMSE), and relative root mean squared error (RMSE r ), for accuracy assessment.

Relationships between SVIs or LiDAR Metrics and AGB (Experiment 1)
The statistical regression analysis was used to model relationships between SVIs and biomass. Accuracy assessment of each modeling case is shown in Table 3. The best statistical model was found by using the OSAVI as the input variable (R 2 = 0.594, RMSE = 37.097 Mg/ha, RMSE r = 27.004%), followed by using the NDVI (R 2 = 0.563, RMSE = 37.702 Mg/ha, RMSE r = 27.445%). The EVI fitted model yielded the lowest R 2 and highest RMSE and RMSE r values (0.209, 52.012 Mg/ha and 37.861%, respectively). Overall, in the experiment 1, statistical models with the single SVI did not show satisfactory performance with the R 2 < 0.6 and RMSE > 37 Mg/ha. Table 4 shows moderate to good modeling performances for AGB estimation by using LiDAR metrics. The variables of p 60 and Cov had the highest (R 2 = 0.696, RMSE = 31.326 Mg/ha, RMSE r = 22.803%) and lowest correlations (R 2 = 0.463, RMSE = 47.298 Mg/ha, RMSE r = 34.43%) with biomass, respectively. The model with the input variable of H mean yielded the next best performance with an R 2 of 0.684 and RMSE of 32.127 Mg/ha. It was noted that statistical models fitted by LiDAR metrics generally resulted in higher R 2 and smaller RMSE values than those fitted by the optical variables. The differences in R 2 and RMSE between the best LiDAR variable based model and the best optical variable based model were approximately 0.1 and 5.77 Mg/ha, respectively.

Relationships between COLIs and AGB (Experiment 2)
As shown in Table 4, the variable of H mean had a relatively good correlation with biomass, following the variable of p 60 . Previous studies have used the LiDAR height metrics H mean as effective prior knowledge in the sampling design [59,60] because H mean can be treated as a direct indicator of the forest structure and growth status. Therefore, two types of the combined optical and LiDAR indices (i.e., COLI1 and COLI2) were generated by integrating H mean with each of the seven SVIs according to Equations (3) and (4), respectively. More details on COLI1 and COLI2 were provided in Section 3.1.2. Table 5 shows AGB predictions could be improved by statistical models that have both H mean and one of the SVIs as input variables, compared to those using either optical or LiDAR data alone. Among all the models in COLI1, the statistical model with input variable of OSAVI×H mean had the best modeling performance with an R 2 of 0.788 and RMSE of 24.787 Mg/ha. In COLI2, the best statistical model (R 2 = 0.763, RMSE = 26.568 Mg/ha) was achieved by using OSAVI_H mean as well. Additionally, it was observed that statistical models in COLI1 generally had slightly better perfomance than those in COLI2. The differences in R 2 and RMSE between the model using the best COLI1 (OSAVI×H mean ) and the model using the best COLI2 (OSAVI_H mean ) reached about 0.03 and 1.78 Mg/ha, respectively. Figure 3 shows AGB comparisons between field observations and model predictions using different data scenarions: (a) OSAVI, (b) H mean , (c) OSAVI×H mean , and (d) OSAVI_H mean . Obviously, the models with combined indices showed fewer deviations in the slopes of the fitted trend-lines (the red lines) away from the 1:1 line (the black lines) than those using either the optical-only or LiDAR-only data.

Multivariate Model Performance (Experiment 3)
The MSLR, KNN, SVR, RF and SSAE algorithms were individually adopted to estimate biomass using ACOLI1, ACOLI2, the combination of ACOLI1, AO and AL, and the combination of ACOLI2, AO and AL. Table 6 shows that the SSAE model has the best modeling performance, followed by the RF model. The MSLR model, in most cases, produced the most error in AGB estimate with the lowest R 2 and highest RMSE values. Figure 4 shows comparisons between field surveys and predicted AGB values derived from the five prediction models calibrated by variables of ACOLI1, AO and AL. Obvious differences in modeling performances among the five prediction algorithms could be identified. The SSAE model showed higher agreements between the model estimates and field observations than other models. Additionally, results showed that the prediction accuracy was improved with increased input variables, i.e., the models with input variables of ACOLI1 (ACOLI2), AO and AL achieved better performances than those with ACOLI1 (ACOLI2). Furthermore, slightly better prediction accuracies were observed for the models calibrated by ACOLI1, AO and AL integrated variables compared to those by ACOLI2, AO and AL. The best performing model, SSAE with ACOLI1, AO and AL integrated variables, reached a R 2 of 0.935, a RMSE of 15.67 Mg/ha and a RMSEr of 11.407%.

Multivariate Model Performance (Experiment 3)
The MSLR, KNN, SVR, RF and SSAE algorithms were individually adopted to estimate biomass using ACOLI1, ACOLI2, the combination of ACOLI1, AO and AL, and the combination of ACOLI2, AO and AL. Table 6 shows that the SSAE model has the best modeling performance, followed by the RF model. The MSLR model, in most cases, produced the most error in AGB estimate with the lowest R 2 and highest RMSE values. Figure 4 shows comparisons between field surveys and predicted AGB values derived from the five prediction models calibrated by variables of ACOLI1, AO and AL. Obvious differences in modeling performances among the five prediction algorithms could be identified. The SSAE model showed higher agreements between the model estimates and field observations than other models. Additionally, results showed that the prediction accuracy was improved with increased input variables, i.e., the models with input variables of ACOLI1 (ACOLI2), AO and AL achieved better performances than those with ACOLI1 (ACOLI2). Furthermore, slightly better prediction accuracies were observed for the models calibrated by ACOLI1, AO and AL integrated variables compared to those by ACOLI2, AO and AL. The best performing model, SSAE with ACOLI1, AO and AL integrated variables, reached a R 2 of 0.935, a RMSE of 15.67 Mg/ha and a RMSE r of 11.407%.

Wall-to-Wall Predictions
The spatial distribution of AGB was mapped by applying the SSAE model (calibrated by ACOLI1, AO and AL integrated variables) to the whole study area ( Figure 5). The predicted values varied from 24.64 Mg/ha to 299.206 Mg/ha, with a mean value of 113.456 Mg/ha. We have masked out the cloud and river areas. The non-forested regions were considered to have low biomass levels. In view of the raw LiDAR information located in the black rectangle of Figure 5 was missing, this part did not participate in the AGB retrieval.

Wall-to-Wall Predictions
The spatial distribution of AGB was mapped by applying the SSAE model (calibrated by ACOLI1, AO and AL integrated variables) to the whole study area ( Figure 5). The predicted values varied from 24.64 Mg/ha to 299.206 Mg/ha, with a mean value of 113.456 Mg/ha. We have masked out the cloud and river areas. The non-forested regions were considered to have low biomass levels. In view of the raw LiDAR information located in the black rectangle of Figure 5 was missing, this part did not participate in the AGB retrieval.

Biomass Estimation Using COLI
In the study, the combined optical and LiDAR indices, i.e., COLI1 and COLI2, were developed

Biomass Estimation Using COLI
In the study, the combined optical and LiDAR indices, i.e., COLI1 and COLI2, were developed by integrating the best-performing LiDAR predictor (i.e., H mean ) with each of the SVIs derived from the Landsat 8 OLI image. Both COLIs (i.e., COLI1 and COLI2) improved the biomass estimation accuracy compared to either the optical-only or LiDAR-only variable. This was in accordance with earlier studies [32,35,61,62]. Most of these studies used other approaches to combine optical and LiDAR data for estimating vegetation parameters. Usually, the common approach is the loosely coupled design [32,35,[61][62][63][64]. For example, Kulawardhana et al. [32] found that the multiple regression models with input variables from the fusion of LiDAR metrics and multispectral vegetation indices had better performances than those calibrated by either the sole LiDAR or multispectral data for vegetation biomass predictions. In Li et al. [62], the statistical model using combined LiDAR and GF-1 data achieved a satisfactory performance for the maize AGB estimation. In this study, new indices were designed, and this approach presented a novel synergistic way of effectively estimating forest biomass. The satisfactory performances of COLIs can be attributed to the ability of multispectral images to provide surface information about canopy density and cover and the ability of LiDAR data to measure forests information about the branches and stems of the trees. Additionally, although COLI1 fitted models generally achieved smaller RMSE and higher R 2 values than COLI2 fitted models, differences in model predictions between them were slight. It was revealed that OSAVI×H mean and OSAVI_H mean could be used for accurate estimates of forest biomass, along with the fact that the other combined indices were also valuable (Table 5). Overall, it was concluded that the proposed COLIs could improve forest biomass predictions since weaknesses and limitations of one predictor might be compensated through the other one. However, the COLIs should experience further investigations and validations before the extensive application.

Biomass Estimation Using Different Prediction Models
To the best of our knowledge, the applications of deep learning models in estimating vegetation parameters are much less common, although they have been investigated by previous studies for many other fields, such as object detection and image classification [38][39][40][41][42][43]. The estimation performance of SSAE deep learning technique outperformed the other four models. This may be attributed to the fact that SSAE can automatically learn deep spatial sparse features of data layer by layer. In SSAE model, the input variables were mapped into another feature space, where the AGB can be accurately predicted. The overcomplete sparse features were derived due to the sparsity in the hidden layer, and the deep learning features were more representative because of stacked layers of neural networks. In future research, the SSAE modeling should be applied separately for different forest types to improve AGB estimation accuracy. Additionally, the RF algorithm achieved more outstanding predictive accuracy than other models besides SSAE, and its satisfactory performance has been demonstrated in previous studies [35,53,54]. The relatively outstanding performance of the RF model depends on its distinctive regression technique and conceptual design which results in the robustness and flexibility with respect to noise and outliers [65]. It should be noted that RF may not work effectively using a small amount of samples on account of the employment of a soft linear boundaries combination [66]. In the study, the prediction accuracy of SVR was worse than that of RF. This result may be attributed to the "noise" variables which might bias the optimal hyperplane in the SVR and thus weaken the model performance. In all the data scenarios, the results show that the SVR achieved a higher accuracy than the KNN. SVR was more powerful at dealing with high-dimensional and nonlinear problems than KNN, which has been suggested by previous studies [35,53]. Besides, the relatively poor performance of KNN might be related to the lack of an effective procedure to weight the prediction strength of the variables in the implementation process of KNN. In the majority of cases, MLSR performed slightly worse than the other tested methods possibly due to its inherent limitations, e.g., the modeling accuracy no longer Remote Sens. 2019, 11, 1459 13 of 17 changes when additional variables are added [67] because not all variables are linearly correlated with biomass.

Overall Performance of the Biomass Estimation
The biomass prediction result of the SSAE model with inputs of ACOLI1, AO and AL integrated variables (R 2 = 0.935, RMSE = 15.67 Mg/ha, RMSE r = 11.407%) in our study showed a high accuracy comparable to or better than other similar studies. Li et al. [62] integrated airborne LiDAR with satellite GF-1 data for estimating maize AGB during peak growing season with an R 2 of 0.69 and RMSE r of 39%. Kulawardhana et al. [32] evaluated the performance of data fusion based on LiDAR metrics and multispectral vegetation indices using multiple regression models for vegetation AGB estimation with an RMSEr of 25.9%. Fassnacht et al. [35] investigated the potential of hyperspectral and airborne LiDAR integrated dataset for forest biomass estimation with a mean R 2 of 0.72 and RMSE r of 28.46%. In Laurin et al. [61], the Partial Least Square Regression (PLSR) with inputs of the combined airborne LiDAR metrics and hyperspectral bands achieved the best performance (R 2 = 0.7, RMSE r = 35.83%) for the forest biomass estimation.
There is the possibility of improving model accuracy in some effective ways. For example, the methods to select the predictor variables should be introduced, and separate estimates for different forest types should be implemented to further improve the AGB prediction performance in future studies.

Conclusions
In this research, univariate regression models using combined optical and LiDAR indices (i.e., COLI1 and COLI2) significantly improved the biomass estimation accuracy compared to those using either data type alone. Therefore, the proposed COLIs provide an avenue to synergistically use optical and LiDAR data for mapping biomass and other vegetation parameters. However, these COLIs should be expanded to other forest types for further validations and wide applications. In addition, it should be noted that the presented results are the first step towards the integration of remotely sensed spatial and spectral information for a precise and non-destructive estimation of forest biomass. Other data fusion methods may further increase the prediction power. The deep learning model (i.e., SSAE), which has previously been rarely explored in mapping forest attributes, was found to be superior over other prediction models for estimating forest AGB. This conclusion may facilitate further applications of deep learning models for mapping vegetation structure parameters. The SSAE prediction model calibrated by variables of ACOLI1, AO and AL yielded the best biomass estimation accuracy. Overall, the deep learning model with inputs of structural and spectral integrated information can become a powerful tool for applications in precision forest biomass monitoring.