Multi-Platform LiDAR for Non-Destructive Individual Aboveground Biomass Estimation for Changbai Larch ( Larix olgensis Henry) Using a Hierarchical Bayesian Approach

: Individual-tree aboveground biomass (AGB) estimation is vital for precision forestry and still worth exploring using multi-platform LiDAR data for high accuracy and efﬁciency. Based on the unmanned aerial vehicle and terrestrial LiDAR data, this study explores the feasibility of the individual tree AGB estimation of Changbai larch ( Larix olgensis Henry) of eight plots from three different regions in Maoershan Forest Farm of Heilongjiang, China, using nonlinear mixed effect model with hierarchical Bayesian approach. Results showed that the fused LiDAR data estimated the individual tree parameters (i.e., diameter at breast height (DBH), tree height (TH), and crown projection area (CPA)) with high accuracies (all R 2 > 0.9 and relatively low RMSE and rRMSE) using region-based hierarchical cross-section analysis (RHCSA) algorithm. Considering regions as random variables, the nonlinear mixed-effects AGB model with three predictor variables (i.e., DBH, TH, and CPA) performed better than its corresponding nonlinear model. In addition, the hierarchical Bayesian method provided better model-ﬁtting performances and more stable parameter estimates than the classical method (i.e., nonlinear mixed-effect model), especially for small sample sizes (e.g., <50). This methodology (i.e., multi-platform LiDAR data and the hierarchical Bayesian method) provides a potential solution for non-destructive individual-tree AGB modeling with small sample size and high accuracy in both forestry and remote sensing communities.


Introduction
Forest aboveground biomass (AGB) serves as the basis for monitoring and accounting for carbon stock and plays a crucial role in regulating the global carbon balance [1][2][3]. Accurate and efficient AGB estimation is required for improving estimates of terrestrial carbon sources and carbon sinks [4]. AGB estimation at the plot level is typically obtained by aggregating the predicted biomass of individual trees within a plot [5]. Incipiently, destructive sampling of trees was required for measuring individual tree AGB, involving tree felling, component cutting, drying, and weighing [6]. Due to strong correlations between individual tree structure parameters (e.g., diameter at breast height (DBH) and tree height (TH)) and individual tree AGB, species-specific allometric equations for various forest types have been continuously developed and widely applied for estimating individual tree AGB [7][8][9][10]. In recent years, researchers have demonstrated that the shape and size of tree crowns are usually associated with the properties of photosynthesis and nutrient cycling, which affect tree growth [11,12]. Adding crown-related structure parameters (e.g., crown length, crown width, crown volume, and crown projection area) to the AGB model may increase estimation accuracy [13]. framework has been developed in forest biomass estimation in recent years [1,44,46]. The hierarchical Bayesian approach could estimate a comprehensive set of equations and yield realistic assessments of parameter estimation uncertainty [47]. This approach can incorporate random variations, such as regional variations, into the model fitting procedure as the NLME in classical statistical methods [48][49][50][51].
Before sampling, the parameters will be given a prior distribution, which is a crucial component of the Bayesian approach. Two priors were introduced in the Bayesian framework: non-informative and informative priors [52,53]. Because important prior knowledge about the data can be readily included in Bayesian analyses, the Bayesian approach with informative priors tends to perform better than the Bayesian approach with non-informative priors and the classic statistical approach [4,51]. Therefore, selecting appropriate prior distributions for all parameters from external knowledge (e.g., reported parameters from literature or parameter estimation using conventional statistical methods) is critical for improving model accuracy.
Furthermore, the Bayesian technique benefits from estimating with small sample sizes, which overcomes the drawback of classical statistical analysis on the sample size of AGB estimation based on stratified data. For example, Mauricio et al. [44] applied the Bayesian method to estimate individual tree AGB with six sample trees, yielding a similar performance to that using a classic statistical method with 40-60 trees. Dimitris-Zianis et al. [54] found that the model efficiency of the Bayesian approach was superior for AGB estimation using six individual Hungarian oak (Quercus frainetto Ten.) trees. In recent years, some studies attempted to combine the Bayesian method with LiDAR data in forestry biomass applications. For example, Ver Planck et al. [55] demonstrated the applicability and practicality of using LiDAR data and the hierarchical Bayesian method to estimate forest AGB at the stand level. Wang et al. [51] established an independent tree AGB model for Qinghai spruce (Picea crassifolia Kom) using airborne LiDAR data with both hierarchical Bayesian and classical methods. However, few studies were devoted to applying the hierarchical Bayesian method with multi-platform LiDAR data to estimate the individual tree AGB with appropriate sample size, particularly for typical tree species in northeast China.
As commercially valuable timber, larch is widely planted in the mountains of north, northeast, and southwest China because of its straight shape and high resistance to bending and cracking [56]. In temperate regions of China, larch forests account for 6.5% of plantation area and 6.77% of forest stock, dominating the forest ecosystem [11]. Many researchers predicted the AGB of individual larch trees using allometric equations with measured tree variables (e.g., DBH, TH). For example, Wang [57] developed linear component biomass equations for ten tree species, including the Dahurian larch (Larix gmelinii) using Ordinary Least Square (OLS) regression. Dong et al. [58] developed two additive biomass equations for three coniferous plantation species (i.e., Korean pine (Pinus koraiensis Sieb. et Zucc.), larch, and Mongolian pine (Pinus sylvestris var. mongolica)) in northeast China using the measurement data of DBH and tree height, which had excellent fitting performance (R a 2 = 0.958-0.989). However, it is still worth exploring how to obtain individual-tree AGB estimation of larch with low cost, non-destructive samples, and high accuracy using multi-platform LiDAR data and the hierarchical Bayesian method [51].
Therefore, the objective of this study is to investigate the applicability of hierarchical Bayesian models for non-destructive individual tree AGB estimation for a typical larch (Changbai larch (Larix olgensis Henry) belonging to Pinaceae) in northeastern China based on the fusion of UAV and terrestrial LiDAR data (U-T LiDAR data hereafter). Specifically, this study was to: (1) estimate individual tree parameters, including DBH, TH, and crown projection area (CPA), by two individual tree segmentation algorithms (CSP and RHCSA) based on U-T LiDAR data; (2) establish and compare five commonly used AGB models of Changbai larch based on estimated individual tree parameters; (3) establish the hierarchical Bayesian models with varying sample sizes for individual tree AGB estimation and compare their performances to the conventional nonlinear mixed-effects model (NLME) method. This study combines the advantages of the hierarchical Bayesian method and multi-platform LiDAR data to provide the forestry remote sensing community with a solution of non-destructive AGB modeling using small but effective sample sizes.

Study Area and Sampling
The study area is located in the Maoershan Experimental Forest Farm, Shangzhi City, Heilongjiang Province, China, from 127 • 29 E to 127 • 44 E and 45 • 14 N to 45 • 29 N (Figure 1a). The slope ranges from 5 • to 25 • , and the terrain is mountainous, rising from south to north with an average elevation of about 300 m. This region has a temperate continental monsoon climate. Maoershan is a typical natural secondary forest in northeastern China surrounded by various broadleaved trees, such as white birch (Betula platyphylla Suk.), Mongolia oak (Quercus mongolica Fisch. ex Ledeb.), and Korean aspen (Populus davidiana), and a few coniferous trees, such as Changbai larch (Larix olgensis Henry), Mongolian pine (Pinus sylvestris var. mongolica Litv.), and Korean pine (Pinus koraiensis Sieb. et Zucc.).
The eight sample plots of 0.09 ha (30 × 30 m) were selected from three larch plantations regions according to different site conditions and forest stages (A1: middle-age forest; A2: near-mature forest; A3: mature forest) (Figure 1b). The normalized UAV-LiDAR and TLS point data of the eight sample plots were shown in Figure 1c.
Changbai larch based on estimated individual tree parameters; (3) establish the hierarchical Bayesian models with varying sample sizes for individual tree AGB estimation and compare their performances to the conventional nonlinear mixed-effects model (NLME) method. This study combines the advantages of the hierarchical Bayesian method and multi-platform LiDAR data to provide the forestry remote sensing community with a solution of non-destructive AGB modeling using small but effective sample sizes.

Study Area and Sampling
The study area is located in the Maoershan Experimental Forest Farm, Shangzhi City, Heilongjiang Province, China, from 127°29′E to 127°44′E and 45°14′N to 45°29′N ( Figure  1a). The slope ranges from 5° to 25°, and the terrain is mountainous, rising from south to north with an average elevation of about 300 m. This region has a temperate continental monsoon climate. Maoershan is a typical natural secondary forest in northeastern China surrounded by various broadleaved trees, such as white birch (Betula platyphylla Suk.

LiDAR Data and Preprocessing
UAV-LiDAR data used in this study were acquired in August 2020. The UAV-borne LiDAR equipment was RIEGL mini VUX-1UAV-LiDAR scanner (Horn, Austria, https://www.riegl.com/ products/unmanned-scanning/riegl-minivux-1uav, accessed on 14 June 2022) carried by a Feima D200 UAV platform (Shenzhen, China, https://www.feimarobotics.com/en/product Detail D200, accessed on 14 June 2022). All the flights were designed as crossed transects with 50% swath overlaps at 80 m altitude and 5.0 m/s speed. Raw UAV-LiDAR data were denoised and then classified into nonground points and ground points using the improved progressive triangulated irregular network (TIN) densification (IPTD) filtering algorithm [59] in the Green Valley
TLS data of eight sample plots were acquired in September 2020 using a Riegl VZ-400i (RIEGL, Horn, Austria, https://www.riegl.com, accessed on 14 June 2022). In order to ensure the spatial coverage of each station, 11 to 13 scanning stations were set up for each plot. TLS data were registered to WGS 84-UTM zone 52N projection coordinate system using the real-time kinematic (RTK) global positioning system receiver. The point cloud data of multiple stations were co-registered and integrated into plot-level point cloud data using the Riegl RiSCAN PRO software package (v 2.7.1). Due to a large amount of TLS point cloud data, the point cloud was thinned by an octree-based algorithm [60] in terms of speed and computational complexity, following noise elimination. Subsequently, the thinned point clouds were filtered into ground and non-ground points using IPTD filtering algorithm in LiDAR360 software. The main parameters of the two LiDAR data used in this study are shown in the Table A1 of Appendix A.

Field Inventory Data
In this study, all the trees with a DBH equal or greater than 5 cm in eight plots were recorded. In total, four individual tree parameters of 370 Changbai larch trees in eight plots were measured and recorded, including DBH (cm), tree height (TH) (m), tree species, and location. There were 123, 98, and 149 trees in middle-age, near-mature, and mature forest, respectively. DBH was measured using a perimeter ruler; tree height was measured using the Vertex IV ultrasound instrument system, and the location of each tree was recorded using RTK with positional error estimated to be within 5 cm. The descriptive statistics of the main variables are presented in Table 1.

Methods
In this study, the preprocessed UAV-LiDAR and TLS data were registered and fused to create the U-T LiDAR data for comprehensively describing individual tree structure. Firstly, based on the U-T LiDAR data, two algorithms (CSP and RHCSA) were used to delineate individual trees and obtain the optimal LiDAR-derived individual tree parameters (i.e., DBH, TH, and CPA) based on field inventory data (objective 1). Secondly, five widely used AGB model forms (model I-V) were selected and compared by NLS and Bayesian approach based on the reference AGB of individual tree, and the corresponding mixed-effects LiDAR-AGB model was developed with random effects for different regions to further improve the model fitting (objective 2). Then, the hierarchical Bayesian and NLME methods were compared to estimate individual tree AGB with five sample sizes (i.e., 10%, 25%, 50%, 75%, and 100%) in order to explore the appropriate sample size that could balance between the estimation accuracy and cost (objective 3). A flowchart of this study is depicted in Figure 2.

Estimation of Individual Tree Parameters
The UAV-LiDAR and TLS data for each sample plot were coarsely registered based on the projected coordinate system. Then, the iterative closest point (ICP) algorithm, a point-based matching method based on minimizing the cumulative distance between two LiDAR data [26,61], was used to conduct a fine co-registration of the two LiDAR data in this study. In this study, UAV-LiDAR was set as the reference point cloud for registering TLS data. After the coarse and fine registration (Figure 3), the UAV-LiDAR and TLS data for each plot were fused in LiDAR360 software.

Estimation of Individual Tree Parameters
The UAV-LiDAR and TLS data for each sample plot were coarsely registered based on the projected coordinate system. Then, the iterative closest point (ICP) algorithm, a point-based matching method based on minimizing the cumulative distance between two LiDAR data [26,61], was used to conduct a fine co-registration of the two LiDAR data in this study. In this study, UAV-LiDAR was set as the reference point cloud for registering TLS data. After the coarse and fine registration (Figure 3), the UAV-LiDAR and TLS data for each plot were fused in LiDAR360 software. This study used two segmentation algorithms (CSP and RHCSA) to segment individual trees based on the normalized fused LiDAR (U-T LiDAR) data. The CSP algorithm proposed by Tao et al. [34] finds seed points of individual trees by recognizing tree bases and then labels other points by finding the shortest path to the seed points, which has been directly applied to LiDAR point data [23,62]. In addition, the RHCSA algorithm was This study used two segmentation algorithms (CSP and RHCSA) to segment individual trees based on the normalized fused LiDAR (U-T LiDAR) data. The CSP algorithm proposed by Tao et al. [34] finds seed points of individual trees by recognizing tree bases and then labels other points by finding the shortest path to the seed points, which has been directly applied to LiDAR point data [23,62]. In addition, the RHCSA algorithm was applied to segment individual trees based on the canopy height model (CHM) derived from the U-T LiDAR data. RHCSA is a one-step individual tree crown delineation (ITCD) algorithm, which segments individual trees with a few user-defined parameters [15]. The appropriate cell size of CHM could reduce the errors of segmentation and ensure that sufficient detail is maintained in the height model [63]. This study employed kriging interpolation with a pixel size of 0.25 m through trial and error. The smoothed CHM eliminated the majority of noise and empty sinks, allowing for an accurate interpretation of tree crowns [64].
The local maxima of height within each tree crown and the projection area of delineated crowns from the segmentation algorithms were defined as TH and CPA, respectively. The DBH of individual trees was estimated using a nonlinear least-squares circle fitting algorithm [65,66]. The position of the treetop was considered to be the tree location. All of the individual tree parameters (i.e., tree location, DBH, TH, and CPA) were estimated using LiDAR360 5.0 and ArcGIS 10.4 (ESRI, Redlands, CA, USA) software.
The accuracy assessment of two segmentation algorithms was evaluated based on the three indices, including recall (r), precision (p), and F-score (F) [30]. The three indices were calculated using the following equations: where TP (true positive) denotes the number of trees correctly detected, that is, 1:1 matched trees; FN (false negative) denotes the number of trees that were not detected; and FP (false positive) denotes the number of trees falsely detected. r indicates the tree segmentation completeness, p indicates the correctness of the detected trees, and F is the overall accuracy considering both commission and omission errors. In this study, 1:1 matched trees were defined as the trees located within 3 m of the reference tree position and had the minimum differences of DBH from references within the 20% of the average plot value [67]. Based on 1:1 matched trees, the accuracy of DBH and TH were evaluated using the field measurements. Since two-dimensional parameter CPA is difficult to measure in the field and no allometric equation is available to estimate CPA for Changbai larch, the accuracy of CPA was evaluated using the CHM-based manual delineation. The accuracies of estimated individual tree parameters were assessed by the coefficient of determination (R 2 ) of the regression between estimated and measured parameters, root-mean-square error (RMSE), and relative root-mean-square error (rRMSE) of estimated parameters as follows. where y i is the reference parameters of individual trees,ŷ i is the estimated parameters, y i is the mean value of the reference parameters,ŷ i is the mean value of estimated parameters, and n is the number of samples.

Establishment of Individual-Tree AGB Model Based on U-T LiDAR
Five common candidate models established with the three optimal individual tree parameters (DBH, TH, and CPA) were applied as the basic U-T LiDAR-AGB models (baseline) in this study [21,68]. The five AGB models are listed as follows.
Mode II : Model IV : where DBH denotes the diameter at breast height (cm), TH is tree height (m), CPA presents tree crown projection (m 2 ), AGB is reference AGB, and ε is an error term. The reference AGB of individual tree in this study was calculated using the additive biomass equations of larch plantations proposed by Dong et al. [69], shown in Table A2 of Appendix A.
Due to the stratified structure of the data (i.e., sample trees in three regions of different forest stages, including A1: middle-age forest, A2: near-mature forest, and A3: mature forest), the NLME approach [70] was used for the basic LiDAR-AGB model in this study. All parameter combinations were simulated as mixed parameters, with Akaike's information criterion (AIC), Bayesian information criterion (BIC), and log-likelihood (LL) serving as the primary evaluation criteria for the fitting performance. The specific mixed-effects U-T LiDAR-AGB model form is established based on the basic U-T LiDAR model, with the mixed-effects model parameters estimated using the maximum likelihood method (ML) with the nlme function of the nlme library in R software 4.0.3 (New York, NY, USA, https://mran.microsoft.com/, accessed on 14 June 2022).

Establishment of Hierarchical Bayesian Model
Based on the advantage of the Bayesian method for small sample estimation and the stratified structure of the data, hierarchical Bayesian models were established. However, to ensure an adequate sample size, this study selected five sample sizes (samples 1-5) using stratified random sampling with a proportional allocation. Sample 1-5 was stratified randomly selected using a proportion of 10%, 25%, 50%, 75%, and 100% in each forest stage (i.e., A1: middle-age forest, A2: near-mature forest, and A3: mature forest), respectively.
The Bayesian approach is a statistical framework that combines new evidence (data) with prior distributions of parameter values to derive new probability for various parameter values [47,71]. It is also suitable for multilevel analysis [72], where the regions (i.e., A1, A2, and A3) are considered as random variables. The AGB data can be used to estimate parameter θ for each region. Let y = y 1j , . . . , y ij represents the AGB data vector, y ij is the AGB of the ith tree in the region j; let θ = (θ 1 , θ 2 , θ 3 . . .) represent the vector of parameters to be estimated. Then π(θ|λ) is determined, where λ is a hyperparameter vector [73]. The inference parameter θ is based on its posterior distribution: The prior distribution of π(θ|λ) in this study is obtained from the parameters estimated by NLME [73]. The posterior distribution is used for Bayesian statistical inference, assuming θ is known, f (y|θ) provides the distribution of y, which is considered a likelihood function when viewed as a function of the parameters.
For the Bayesian estimation, a burn-in period of 30,000 steps and 330,000 iterations were used to estimate parameters. The thinning parameter was set to 3 to reduce the correlation between neighboring iterations. In this study, Bayesian models were developed using the MCMC procedure in SAS 9.4, and the average running time for each model was approximately three minutes.

Model Evaluation
The best AGB model was selected by both the classical method and Bayesian method from the smallest AIC, BIC, and deviance information criterion (DIC) when pooling the data. Moreover, the stationarity test of Heidelberger-Welch Diagnostics [74,75] was conducted to test whether the model converged in this study. The DIC is calculated by: where D is the posterior mean of the deviance (−2 × Log likelihood of the given data and parameters), and p D is the model complexity, which is summarized by the effective number of parameters. The Fit Index (FI) [76,77] and the root mean square error (RMSE) were applied to compare the Bayesian method with the classical method on different sample sizes. Larger FI and smaller RMSE values indicate a better model fitting. The FI are calculated by: where y i ,ŷ i and y i represent the observed, estimated, and mean values of AGB, respectively.

Estimation of Individual Tree Parameters
The average F (F-score) value of eight sample plots for the two segmentation algorithms was above 0.90 ( Table 2). The CPS algorithm was slightly better than the RHCSA algorithm for individual tree segmentation (0.92 vs. 0.90). A total of 337 sample trees were correctly segmented and matched using the CSP algorithm from 370 reference trees. The RHCSA algorithm, based on CHM for individual tree segmentation, had weaknesses in detecting small trees: the number of correctly detected trees (TP) of RHCSA was slightly fewer than the CSP algorithm (327 vs. 337 out of 370). Note: * represented the average r, p, and F value of eight sample plots; TP, FP, and FN were the total value of true positive (1:1 matched trees), false positive, and false negative for eight sample plots, respectively. CSP: comparative shortest-path algorithm, RHCSA: region hierarchical cross-sectional analysis algorithm. Table 3 presents that the R 2 of individual tree parameters are all greater than 0.9, except the CPA estimated by CSP. It is because that lots of points from shrubs near the ground were misclassified as tree points by the CSP algorithm, and the CSP algorithm calculated CPA based on the average crown diameter estimated by the tree points, therefore resulting in larger estimated canopy diameters and CPA. The RMSE and rRMSE values of TH and CPA estimated by RHCSA were lower than those estimated by CSP. In addition, the accuracies of DBH estimated by CSP and RHCSA were very similar (R 2 : 0.983 vs. 0.990, RMSE: 1.017 vs. 1.024; rRMSE: 4.9 vs. 4.8). To ensure the accuracy of individual tree AGB estimates, we applied the individual tree parameters of 327 sample trees based on the RHCSA algorithm as the predictor variables for the subsequent individual tree AGB estimation.  Table 4 shows the fitness statistics of the five widely used AGB models using both classical and Bayesian approaches in this study. Model V showed the best evaluation statistics with the smallest AIC (3039.067), BIC (3058.017), and DIC (3032.190) values among the five models, while Model IV had the worst model performance and even failed the Heidelberger-Welch Diagnostics test for the stationarity. Model V (Equation (11)) with three predictor variables (i.e., DBH, TH, and CPA) was selected as the basic U-T LiDAR-AGB model due to the best model performance (Table 4). According to AIC, BIC and LL values of all combinations of parameter and random effect (i.e., A1, A2, A3) (see the Table A3 of Appendix A), mixed-effects U-T LiDAR-AGB models were established by adding a random effect parameter to a 4 based on model V as Equation (15). The model evaluation statistics of all parameter combinations are shown in Table A1.

Establishment of Individual-Tree AGB Model Based on U-T LiDAR
where a 1 -a 4 are the parameters of the model, a i is the random-effect parameter of the ith region; AGB ij , DBH ij , TH ij , CPA ij and ε ij are the individual tree AGB (kg), DBH (cm), TH (m), CPA (m 2 ), and random error term of tree j from the ith region, respectively. Table 5 shows the fitting goodness of the basic U-T LiDAR-AGB model (model V) using nonlinear least squares (NLS) and NLME methods. All model parameters were significant, and the NLME method fitting substantially improved. The FI of the mixed-effects model was higher than that of the basic model (0.981 vs. 0.979), and RMSE, AIC, and BIC values decreased. However, there was no significant difference in the standard deviations of the parameters between the two methods. Both models showed an ideal fitting effect.

Establishment of Hierarchical Bayesian Models with Different Sample Sizes
According to stratified random sampling with proportional allocation, five sample sizes decreased from 327 (100%) to 34 (10%) trees were applied for hierarchical Bayesian modeling in this study (see the Table A4 of Appendix A). The basic statistics of individual tree parameters (DBH, TH, CPA, and AGB) for each sample size are summarized in Table A4.
The hierarchical Bayesian models were compared to the corresponding NLME models using different sample sizes. The prior parameter distribution of Bayesian estimation was obtained from the parameters of NLME estimation. The posterior probability distributions of the parameters estimated by the Bayesian theory for U-T LiDAR-AGB models are summarized in Table 6. Although the parameter estimates using the two approaches were very similar, the standard errors of the estimates obtained by the hierarchical Bayesian method in all sample sizes were smaller than those obtained by the classical NLME method. This implies that the Bayesian method provided more stable estimates than the classical method. Figure 4 compares the parameters estimated by the hierarchical Bayesian and NLME methods with the 95% confidence intervals based on the five sample sizes. The hierarchical Bayesian method produced estimates with a narrower 95% confidence interval than the NLME method for all sample sizes. The parameter estimates were more stable for smaller sample sizes, especially the sample size of 82, than for larger sample sizes (e.g., 327) when the hierarchical Bayesian method was applied. Note: The values in parentheses presented the standard errors or standard deviation of the parameters estimates obtained by the hierarchical Bayesian approach or NLME approach, respectively.  Table 6 also confirms that the difference between FI and RMSE between the two methodologies became obvious with decreasing sample sizes. For example, when pooling the data (n = 327), the FI and RMSE of the Bayesian method were very similar to that of the classical method (0.980 vs. 0.981; 24.863 vs. 24.164). When the sample size was reduced to 34, however, the performance of the Bayesian method was much better than that of the classical method: the RMSE decreased by 31.4%.

Discussion
Although it is a widely used method to estimate individual tree AGB using an allometric equation based on field inventory data (e.g., DBH, TH), its efficiency and accuracy are still unsatisfactory, particularly for large-scale forest inventory [78]. This study provided a potential solution for establishing an efficient model using combined UAV and terrestrial LiDAR data with an appropriate sample size based on the hierarchical Bayesian method.  Table 6 also confirms that the difference between FI and RMSE between the two methodologies became obvious with decreasing sample sizes. For example, when pooling the data (n = 327), the FI and RMSE of the Bayesian method were very similar to that of the classical method (0.980 vs. 0.981; 24.863 vs. 24.164). When the sample size was reduced to 34, however, the performance of the Bayesian method was much better than that of the classical method: the RMSE decreased by 31.4%.

Discussion
Although it is a widely used method to estimate individual tree AGB using an allometric equation based on field inventory data (e.g., DBH, TH), its efficiency and accuracy are still unsatisfactory, particularly for large-scale forest inventory [78]. This study provided a potential solution for establishing an efficient model using combined UAV and terrestrial LiDAR data with an appropriate sample size based on the hierarchical Bayesian method.

Individual Tree Parameters Estimation Using U-T LiDAR Data
The near-surface LiDAR is a non-destructive technology widely used in forest inventory and AGB estimation in recent years [19,21,51]. However, single platform/single-scan LiDAR data have limitations; the U-T LiDAR data combines the rich point cloud beneath the canopy obtained from TLS with middle and upper canopy information from UAV-LiDAR. U-T LiDAR data could accurately capture tree trunk parameters (e.g., DBH), TH, and crown structural parameters (e.g., CPA) with a broader viewing angle than field measurements, resulting in a satisfactory AGB estimation accuracy for individual trees.
The accuracy of individual tree segmentation may significantly impact the acquisition of individual tree parameters [26]. Based on the U-T LiDAR data, CHM-based and point cloudbased segmentation algorithms (RHCSA and CSP) were used to segment individual trees in this study. The average plot-level F-scores (0.90 and 0.92 using RHCSA and CSP, respectively) were similar but slightly higher than the other studies. For example, Li et al. [30] reported an average F-score of 0.90 for mixed conifer forests based on the airborne LiDAR data using the PCS algorithm. Tao et al. [34] reported an average F-score of 0.87 for three forest types (i.e., conifers, broadleaf, and mixed) using the CSP algorithm based on terrestrial and mobile LiDAR data. Similar to the previous studies, compared to the PCS algorithm, the CHM-based methods (RHCSA in this study) had weaknesses in detecting small trees, especially for high canopy closure forests [79]. A total of 91.1% and 88.4% of sample trees were correctly segmented and matched by the CSP and the RHCSA algorithm. However, the CPA estimated using the CSP algorithm had much lower accuracy than that of the RHCSA algorithm. This could be because the CSP is a bottom-up algorithm; points from shrubs near the ground could be misclassified as sample tree points, resulting in trees with larger estimated canopy diameters. Furthermore, the CSP algorithm computed the circle area as CPA using the average crown diameter derived from the point cloud. In contrast, the RHCSA algorithm delineated entire tree crowns using CHM from a vertical perspective, making it easier to visualize and estimate CPA than CSP. Due to the similar reasons, the tree height estimated using CSP had a slightly lower accuracy (R 2 = 0.923, RMSE = 1.494 m, rRMSE = 8.3%) than that of RHCSA (R 2 = 0.934, RMSE = 1.247 m, rRMSE = 7.3%). Little difference was observed in the accuracy of DBH extracted by the two algorithms. The results demonstrated the potential of applying multi-platform LiDAR data to derive individual tree parameters for individual tree AGB estimation in lieu of field measurements to some extent.

The Hierarchical Bayesian Method in AGB Estimation
With the widespread use of the MCMC algorithm and the advancement of computing technology, Bayesian-based applications in forestry have become increasingly popular in recent years [41,46,55]. Although the Bayesian and classical methods in this study showed a similar trend in the basic LiDAR-AGB model selection (Model V was the best, followed by Model II, and Model IV was the worst), the Bayesian method is more rigorous ( Table 4). The findings were consistent with that of Fu et al. [21]. However, Wang et al. [51] selected a basic LiDAR-AGB model with DBH and CPA as predictors (model III in this study) due to the LiDAR data obtained from a different platform (i.e., airborne) and the method of estimating individual tree parameters.
The mixed-effects and hierarchical Bayesian models are versatile and applicable for constructing regional biomass models [1]. Based on the stratified data structure, the NLME method of the frequentist paradigm and the hierarchical Bayesian method of the Bayesian paradigm were implemented and compared using regions as the random effect. The hierarchical Bayesian model assumes that parameters are defined from prior distributions. In previous studies [4,80], non-informative priors with enormous or infinite variance were commonly used, and the parameter estimation was always identical to the frequentist paradigm. In this study, the prior distribution of the hierarchical Bayesian method was obtained from the parameters estimated by the classical statistical approach (NLME) according to Zhang et al. [81] and Wang et al. [51]. Alternatively, parameter prior information of some studies was derived from external knowledge (reported parameters from the literature) (e.g., [1]). The impact of prior information on model fitting should be investigated in future research. In addition, the parameters estimated by the hierarchical Bayesian method were more stable than the classical method, with smaller standard errors/deviations and narrower 95% intervals (Table 6 and Figure 4). That is consistent with the results of Wang et al. [51]. Based on the evaluation statistics of FI and RMSE, with the decreasing sample size, the Bayesian method showed obvious advantages (Table 6). For the sample size of 34, the RMSE was a 31.4% reduction compared to the classical method. This is in agreement with the findings of Wang and Zapata-Cuartas [44,51].
This study found that the LiDAR-AGB model with the U-T LiDAR-derived individual tree parameters as predictor variables could achieve better prediction results using the hierarchical Bayesian method, particularly when the sample size was small (34 trees in this study). The parameter estimations were more stable, and the model fitting was better than the classical method. Similarly, Wang et al. [51] estimated the individual tree AGB based on 39 sample trees of Qinghai spruce (Picea crassifolia Kom) using the layered Bayesian method, achieving a result of R 2 = 0.8611. Thus, there is great potential for using the hierarchical Bayesian method with small sample sizes (less than 50) to estimate the individual tree AGB.

Limitations
There are still limitations that need to be discussed. First, it should be recognized that field-estimated AGB is commonly referred to as reference data in remote sensing studies (e.g., [51,55]) in order to avoid destructive biomass sampling in the field and improve efficiency. However, the accuracy of field estimates is unknown, indicating that the 'ground truth' of the individual tree-level AGB may have uncertainty [82]. Although both errors due to tree measurement and the choice of an allometric model are sources of AGB estimation uncertainty, several pieces of research have revealed that the most significant sources of error are currently related to the choice of an allometric model [83]. More work should be devoted to considering these uncertainties derived from remotely sensed data. Second, the efficacy of TLS data in forest stands with dense canopies is inferior to that of BLS due to the inconvenient transport of TLS equipment to sample plots or the difficulty in changing stations in dense forests. In this study, for instance, about ten scanning stations were set up for each sample plot (30 × 30 m), requiring two personnel to operate for 2-3 h to complete the TLS data collection of each plot, with an average stand density of 810 trees/ha. However, BLS may have relatively low accuracy for individual tree parameters due to the lack of points for the middle and upper parts of trees [27]. In practice, the number of TLS scanning stations should be adjusted according to study objectives and stand density to balance efficiency and accuracy. Compared to conventional measurements, combining ULS and TLS data provides a comprehensive, efficient, and cost-effective method for obtaining individual tree properties. Moreover, although the hierarchical Bayesian approach has improved the stability of model parameters compared with the classic NLME method, it does not show significant advantages in model parameter estimation, especially for large sample sizes. Therefore, the hierarchical Bayesian method should be applied according to the study objective and sample size and is still worth further exploration in forestry.

Conclusions
This study investigated the applicability of fusing UAV and terrestrial LiDAR data and the feasibility of hierarchical Bayesian for non-destructive individual tree AGB estimation of Changbai larch (Larix olgensis Henry) in northeastern China. When taking the full advantages of UAV and terrestrial platforms, the U-T LiDAR data estimated the individual tree parameters (i.e., DBH, TH, and CPA) of Changbai larch with high accuracies (all R 2 > 0.9) using RHCSA segmentation algorithm. Considering regions as random variables, the nonlinear mixed-effects U-T LiDAR-AGB model with three predictor variables (i.e., DBH, TH, and CPA) performed better than its corresponding nonlinear model. In addition, the hierarchical Bayesian method provided better performances and more stable parameter estimates than the classical method, especially for small sample sizes (e.g., <50). This study implied that the fused multi-platform LiDAR data combined with the hierarchical Bayesian method have the potential to improve the accuracy of individual tree AGB estimation, which provides a non-destructive approach for individual tree AGB modeling with a small sample size. This methodology of this study (i.e., multi-platform LiDAR data and the hierarchical Bayesian method) could provide a scientific basis for non-destructive individual tree AGB modeling with a small sample size and high accuracy in the future.      Note: DBH presented LiDAR-derived diameter at breast height, TH presented LiDAR-derived tree height, CPA presented LiDAR-derived crown projection area.