Predicting Height to Crown Base of Larix olgensis in Northeast China Using UAV-LiDAR Data and Nonlinear Mixed Effects Models

: As a core content of forest management, the height to crown base (HCB) model can provide a theoretical basis for the study of forest growth and yield. In this study, 8364 trees of Larix olgensis within 118 sample plots from 11 sites were measured to establish a two-level nonlinear mixed effect (NLME) HCB model. All predictors were derived from an unmanned aerial vehicle light detection and ranging (UAV-LiDAR) laser scanning system, which is reliable for extensive forest measurement. The effects of the different individual trees, stand factors, and their combinations on the HCB were analyzed, and the leave-one-site-out cross-validation was utilized for model validation. The results showed that the NLME model signiﬁcantly improved the prediction accuracy compared to the base model, with a mean absolute error and relative mean absolute error of 0.89% and 9.71%, respectively. In addition, both site-level and plot-level sampling strategies were simulated for NLME model calibration. According to different prediction scale and accuracy requirements, selecting 15 trees randomly per site or selecting the three largest trees and three medium-size trees per plot was considered the most favorable option, especially when both investigations cost and the model’s accuracy are primarily considered. The newly established HCB model will provide valuable tools to effectively utilize the UAV-LiDAR data for facilitating decision making in larch plantations management.


Introduction
The forest biome is valuable for providing abundant wood resources, protecting wildlife habitats, storing a high amount of carbon, regulating micro-and macro-climates, and possessing other numerous ecological functions. It also plays a vital role in the terrestrial ecosystem, in which its fluctuation highly affects the terrestrial biosphere and other surface processes [1,2]. Traditionally, forest resource inventory relies on field sample measurements by collecting and summarizing the tree-level attributes within a designated sampling area through tree-by-tree measurements. Ground-based forest inventory has a high value of precision, consequently needing an immense amount of manpower, material resources, and time [3]. Hence, delivering a prompt forest inventory becomes a problematical obstacle, specifically for the inventory of an extensive forest area [4].
The continuous development of remote sensing technology has brought a promising solution to punctuality and the high-spatial limitation, providing a breakthrough for highly efficient forest inventory [5,6]. Optical remote sensing data has achieved effective results in large spatial extents for stand age identification, volume estimation, and biomass mapping [7,8]. Nevertheless, the optical sensors carried by this passive remote sensing platform may have serious signal saturation dilemmas, leading to the deviation of some forestry parameter extraction results [9]. As a burgeoning active remote sensing technology, airborne light detection and ranging (LiDAR) can provide precise measurements of vertical forest structure. Combined with the ground-based sample data, LiDAR data can obtain more accurate distribution characteristics of forest types and a wider range of three-dimensional spatial structures. There is great potential in estimating forest inventory parameters by analyzing the forest structure dynamics through multi-temporal remote sensing data series [10,11]. Various studies have demonstrated that airborne LiDAR has a special benefit in forest fire behavior prediction, canopy structure extraction, volume and biomass estimation, and other forest attributes and ecosystem structure measurement, in which the accuracy is beyond the reach of its traditional counterparts [12][13][14][15].
For the past few years, with the rapid development of unmanned aerial vehicle light detection and ranging (UAV-LiDAR) remote sensing technology, the research on estimating stand structure parameters by UAV-LiDAR is increasing gradually [16,17]. Compared to the conventional airborne data, the UAV-LiDAR has quite a few advantages, such as lower data acquisition cost, higher pulse sampling density, simpler operational procedure, and higher flight route flexibility [18,19]. Hence, it is gradually becoming a powerful tool for 3D forest mapping and convincingly appears as a low-cost remote sensing alternative to airborne and satellite platforms [20,21]. The high-density point clouds UAV-LiDAR system can obtain relatively high precision for tree height and crown structure measurement [22]. The high spatial resolution of UAV-LiDAR provides high-precision single-tree segmentation which leads to stronger data availability and a more accurate algorithm [23,24]. At present, the extraction of tree height and crown information is crucial for the single-tree-based forest inventory [25,26].
Crown is a crucial component of an individual tree that facilitates a material exchange and energy transformation between the forest and environment, greatly determining trees' vegetative space, such as sunlight, photosynthesis, and water utilization [27]. The height to crown base (HCB) of a standing tree is defined as the distance from the first living branch at the crown's base to the ground position [28]. HCB reflects the crown's vertical structure and is closely related to the number of foliage [29]. As an essential indicator of an individual tree's crown characteristics, HCB not only determines the yield, productivity, and growth vigor of an individual tree [30] but also reflects the competition status of a tree within the stand along with other non-natural factors [31][32][33]. HCB is also a decisive variable in the stand growth and yield model and plays a pivotal role in the crown width model [34,35], crown shape model [36], biomass model [37], and fire behavior model [38]. However, in the actual measurement of HCB, especially in stands with large canopy density and complex forest conditions, the measurement accuracy is poor, and thus the efficiency is low [39,40]. Therefore, using LiDAR data to accurately and efficiently predict HCB in a large-scale inventory project is particularly momentous for effective forest management.
In recent years, many studies have used LiDAR point cloud data to predict HCB, in which the method is mainly divided into direct and indirect methods. The direct method is to directly derive HCB from LiDAR data using a canopy, canopy approximation, or percentile ranking based on a polygon or voxel [41][42][43], and the predicted distribution of HCB is based on the descriptive statistics of LiDAR data [44]. However, the prediction efficiency largely depends on the penetration of laser sensors and forest conditions; hence, it is difficult to be transformed for other applications. On the other hand, the indirect method utilizes a model which can predict HCB from other tree and/or stand-level variables [35,38,[45][46][47][48]. The model can accurately predict HCB in other forests with similar species and stand and site conditions. The majority of the developed HCB models use ordinary least squares (OLS) for parameter estimation. However, the data used in modeling is often hierarchically structured (multiple measurements within the same subject, such as plots within sites). Hence, these measurements will be correlated, yielding a relatively high prediction error in the resulting Remote Sens. 2021, 13, 1834 3 of 20 model's deviation. The mixed-effect modeling method can effectively solve this problem and greatly improve the prediction accuracy of the model. A mixed-effect model requires both fixed and random parameters to be simultaneously estimated, allowing variation from each level to be modeled [49]. The fixed parameters describe the variation between covariates and experimental treatment, while the random components describe the data's correlation and heterogeneity [50,51]. Maltamo et al. (2018) utilized the linear mixed effect (LME) approach to develop an HCB model using several predictors, such as the LiDAR-derived diameter at breast height (DBH) and tree height percentile data, along with the field-measured total tree height data [52]. Yang et al. (2020) considered the potential linear and nonlinear relationship between HCB and other specific predictors to establish and compare the applicability of the LME and nonlinear mixed effect model (NLME) using the LiDAR-derived tree height and crown width along with the field-measured DBH as predictors [46]. Nevertheless, few studies considered the nested data structures when UAV-LiDAR is applied for data extraction. Furthermore, the field-measured data is still necessary for some of the model's predicted variables.
Korean larch (Larix olgensis) has high economic and ecological value, such as fast growth, excellent wood properties, and strong resistance to insects and diseases, which are widely distributed in Northern China, North Korea, Japan, and Russia [53,54]. As one of the main afforestation and reforestation tree species, the Korean larch plantation covers an area of 3.16 million hectares (accounting for 5.54% of the total plantation area in China), and its volume reaches 23,700 m 3 (accounting for 7.01% of the total plantation volume in China) [55,56]. Timely and effective data acquisition of a series of a tree or stand attributes is crucial to the rational management of larch forests.
Therefore, this study aims to develop a two-level nonlinear mixed effect model for Larix olgensis to explain the difference in each variable's influence on the HCB at both site and plot level, in which all predictors were extracted from LiDAR data. The specific objectives of this study are to (1) extract tree-and stand-level attributes from UAV-LiDAR point cloud data; (2) propose a two-level nonlinear mixed effect model (NLME) framework based on UAV-LiDAR-derived metrics; (3) calibrate the established NLME model using the field-measured data within diverse test areas; and (4) assess the models' predictability using the site-level leave-one-out cross-validation method. This study is expected to improve forestry investigation efficiency, economize survey cost, and provide guidance for future research and forest management decision-making.

Study Area
This study was conducted in the Mengjiagang Forest Farm (130 • 32 42 −130 • 52 36 E, 46 • 20 16 −45 • 30 50 N), located in the northeastern area of the Huanan County, Heilongjiang province, China. The forest farm is specifically situated in the western foothills of Wanda mountain. The slope is relatively gentle, most of which is between 10 • and 20 • . The terrain is higher in the northeast and lower in the southwest, with a maximum, minimum, and average altitude of 575 m, 170 m, and 250 m above sea level, respectively [57]. This area pertains to a temperate continental monsoon climate, and the soil type is mostly dark brown forest soil. The forests are primarily dominated by artificially grown coniferous trees, including Larix olgensis, Pinus sylvestris var. mongolica, and Pinus koraiensis.

Field Measurement Data
The field survey was carried out in July 2019; a total of 118 square sample plots (30 m × 30 m) with normal growth were established in 11 sites (Figure 1) for long-term growth monitoring (the descriptive statistics of field-measured data are shown in Table 1). The initial planting density of all sites was 3300 stems/ha (spacing 2 m × 1.5 m). The young forest was treated with crown thinning (interval 3−5 years) to adjust the stand's composition and density. The middle-aged forest and near-mature forest were treated with growth tending (interval 6−10 years) to promote the growth of trees. The intensity of crown thinning is 25−45%, or the accumulation intensity is 15−30%; the intensity of growth tending is 15−30%, or the accumulation intensity is 10−20% [58]. All trees in sample plots with DBH greater than 5 cm were measured by diameter tape. The four-directional canopy width was measured by steel tape, and the values were averaged to obtain the crown width. The total tree height and HCB were measured by Vertex IV Ultrasonic Hypsometer made by Haglöf Sweden. The relative coordinates of the sampled trees were measured according to their relative position to the corner of the plots. In addition, the geographic coordinates of the individual trees and four corners of each plot were thoroughly measured with a real-time kinetic (RTK) global navigation satellite system (GNSS) (UniStrong G10A, China), except for trees having poor GNSS signals which were georeferenced by their relative coordinates.
according to their relative position to the corner of the plots. In addition, the geographic coordinates of the individual trees and four corners of each plot were thoroughly measured with a real-time kinetic (RTK) global navigation satellite system (GNSS) (UniStrong G10A, China), except for trees having poor GNSS signals which were georeferenced by their relative coordinates.

Unmanned Aerial Vehicle Laser Scanning Data Acquisition
The RIEGL mini VUX-1UAV LiDAR scanner mounted on the Feima D200 UAV platform was used to obtain UAV-LiDAR data from the 11 sites (from 10-12 July 2019). The working pulse repetition frequency, maximum measuring range, and maximum scanning speed of the scanner were 100 kHz, 250 m, and 100 scans per second, respectively. The scanning angle was controlled within ±60 • to reduce the measurement error caused by an immense angle. The flight speed was 5.0 m/s, maneuvering around 80 m of above-ground altitude, and the air route was a cross-transect with overlapping strips of 80 m. The average point cloud density results for each site were between 150 and 270 pt./m 2 .

UAV-LiDAR Metrics Extraction
The UAV-LiDAR data processing and metrics extraction can be found in Figure 2. The noise points were masked off manually from the raw LiDAR data, and the remaining points were divided into ground points and non-ground points by using cloth-simulated filtering [59]. To generate digital terrain models (DTMs), the ground points were interpolated by a Kriging spatial interposition method with a 0.5 m pixel size [60]. Then, the normalized height of each point was obtained by subtracting the DTM value from the elevation of all points [61]. We applied a canopy filtering method called graph-based progressive morphological filtering (GPMF), which was used to obtain the canopy height model (CHM) from LiDAR data. GPMF might prevent data pits obtained by traditional methods from damaging the integrity and smoothness of the tree canopy, leading to a large error in the extraction of tree parameters [62]. The resolution of CHM was set to 0.1 m. The individual tree canopies were detected automatically from the CHM using the region-based hierarchical cross-section analysis (RHCSA) algorithm [63]. This algorithm treats CHM as a mountain terrain and uses the crown's horizontal relation in the vertical direction to detect a single tree.

Base HCB Development
Following previous field-based HCB modeling studies [45], we tested three candidates of HCB models [65,69] for model development (logistic and exponential form). Instead of DBH, we applied as the primary predictor ( = ) since DBH mostly cannot be directly extracted from aerial point cloud data due to the canopy obstruction. The following exponential model was originally formulated by Wykoff, et al. (1982) [65] and was found to be the most suitable for our dataset: Tree-to-tree matching generation was performed between field measurements and segmented trees according to spatial position and height difference [2]. If the segmented tree was located within a circular buffer (the corresponding crown radius) around the reference point and the height difference was less than 20% of the plot's highest height, the segmented tree was designated as the reference tree candidate [64]; then, a unique candidate or nearest individual in multiple candidates was selected as a different candidate to match the reference tree. Eight thousand seven hundred eighty-five (8785) trees in 118 plots were matched correctly with the field data, and the detection rate was 56-100% (mean 76%). After the automatic matching generation, the dead trees and other irrelevant data were removed from the dataset. Finally, the remaining 8364 trees were correctly matched with the field-measured data and utilized for subsequent modeling.
Previous studies have shown that HCB is significantly affected by tree size, competition, and stand characteristics [31,65,66]. In this study, those three categories of predictors were generated from UAV-LiDAR data ( Table 2) and used to develop a generalized HCB estimation model (Section 2.5). The indicators of tree size-including LiDAR-derived tree height and crown width-were defined as the maximum height of all LiDAR pulses and calculated by 2 × √ crownarea/π, where crownarea was the area of the convex hulls of delineated crowns [15]. Furthermore, distance-independent competition indices were calculated from LiDAR data [67]. The LiDAR-derived relative dimensions of tree height and crown projection area were firstly introduced, including the ratios of a target tree's height to the maximum and mean tree height (RH max and RH mean ) and the ratios of a target tree's crown area to the maximum, mean, and total crown area (RCA max , RCA mean and RCA total ), respectively. Secondly, a measure of competition based on crown areas evaluated at a certain percentage of crown length was expanded and calculated using LiDAR data [68]. The ratio between the subjected and the total crown areas computed at a reference height equal to p% of the height of the subject tree (h p ) was calculated as a competition index. Finally, the plot-level metrics were generated to characterize the stand conditions, including height percentiles (H P 5 , H P 10 , . . . , H P 99 ), variance, standard deviation and coefficient of variation of height (H P var , H P std , H P cv ), skewness and kurtosis of height (H P skw , H P kur ), and the proportion of points above the corresponding percentiles (H P 5 , H P 10 , . . . , H P 99 ) to the total number of points within a plot (D P 5 , D P 10 , . . . , D P 99 ). The candidate variables of this study are given in Table 2. the ratio of the crown area above p% relative height of the target tree to the sum of all crown areas above this height in the sample-plot RH mean the ratio of the total height of the target tree to the mean total height in the sample-plot RH max the ratio of the total height of the target tree to the maximum total height in the sample-plot RCA mean the ratio of the crown width of the target tree to the mean crown width in the sample-plot RCA max the ratio of the crown width of the target tree to the maximum crown width in the sample-plot RCA total the ratio of the crown width of the target tree to the total crown width in the sample-plot Following previous field-based HCB modeling studies [45], we tested three candidates of HCB models [65,69] for model development (logistic and exponential form). Instead of DBH, we applied CW as the primary predictor (X = CW) since DBH mostly cannot be directly extracted from aerial point cloud data due to the canopy obstruction. The following exponential model was originally formulated by Wykoff, et al. (1982) [65] and was found to be the most suitable for our dataset: where HCB and H is height to crown base and total tree height, respectively; X is vectors of stand or tree variables; and β is the estimated parameter vector.
H and CW are closely related to HCB, affecting photosynthesis and interspecific competition of trees, and both can be directly extracted from LiDAR data; hence, they are often used as a predictive variable for the HCB model [4]. As an effort to improve the fitting effect and prediction accuracy of the model, we introduced covariates to reflect the stand quality and competition factors in constructing the based model (Table 2), the predictors could be extended as: where CW, Stand., and Comp. is the crown width, stand quality factors, and competition factors, respectively. All covariates were significant and had consistent correlation, which were selected using optimal subset regression. To ensure the simplicity of the model and prevent excessive parameterization and collinearity, one covariate within each variable group was introduced into the model.
The models were fitted to the whole data using nonlinear least-squares regression in R4.0.3 software. Several statistical criteria were used to select the best fitting performance model, including the coefficient of determination (R 2 ), the mean difference (Bias), root mean square error (RMSE), and Akaike information criterion (AIC) as in Dong, et al. (2014) [70].

Two-Level NLME HCB Model
A two-level NLME HCB model was further introduced to consider the random interference of both site-and plot-level variation. The model expression is given below [71]: where the indices i, j, k are the site-level, the plot-level within the site-level, and the observation of an individual tree, respectively; HCB jik is the height to the crown base of the k th tree on the j th plot within the i th site; M is the number of site; M i is the number of the sample plots within the i th site; n ij is the number of trees on the j th plot in the i th site; f (.) is a real-valued and differentiable function of a plot-specific parameter vector ϕ ijk and a covariate vector x ijk ; and ε ijk represents the within-group error with zero expected value and follows a normal distribution, and has R ij as the positive-definite variance-covariance structure. Furthermore, ϕ ijk can be expressed as: where β is the p-dimensional fixed-effect parameter vector; u i and u ij are the site-and plot-level random-effects parameter, which assumed to obey the normal distribution with the expectation of zero value and the variance-covariance of ψ 1 and ψ 2 , respectively; and A ijk , B i,jk , and B ijk are respectively the design matrices corresponding to β, u i , and u ij . ε ijk , u i , and u ij are mutually independent. The most important step for applying a mixed-effects model is to determine which parameters are categorized as fixed effects and random-effects parameters. In this study, all parameter combinations were simulated as mixed parameters with Akaike's information (AIC), Schawarz's Bayesian information criterion (BIC), and likelihood (LL) as the main criteria to evaluate the fitting performance.
Further analysis was carried out by selecting a mixed-effect parameter combination having the smallest AIC, BIC, and LL. We performed a likelihood-ratio test (LRT) to avoid overparameterization [72].
The predictions residuals of the NLME HCB model with both site-and plot-level random effects were analyzed for potential spatial autocorrelation and heteroscedasticity. The preliminary analysis showed that heteroscedasticity was detected, but there was no spatial autocorrelation between the observed values. To solve this problem, the following error term variance-covariance matrix structure was analyzed and applied: where R ij is the variance-covariance matrix of the error term ε ij in the j th plot within the i th site (i = 1, . . ., M, j = 1, . . . , M i ); σ 2 represents the scaling factor of the sample plot error dispersion; G 0.5 ij is a n ij × n ij dimensional diagonal matrix which accounts for the heteroscedasticity of data in the sample plot; Γ ij is a n ij × n ij dimensional matrix explaining within plot autocorrelation structure of errors. The Γ ij was supposed as an identity matrix since there was not any spatial autocorrelation detected. Therefore, only the effect of heteroscedasticity needed to be considered on the model. Three commonly used variance stability functions were evaluated and compared, including exponential function, power function, and constant plus power function. The results showed that the power variance function with H as the independent variable (Equation (6)) effectively explained the heteroscedasticity in our data.
var ε ijk = σ 2 H 2γ ijk (6) where H ijk is the tree height derived from LiDAR data of the k th tree within the j th plot and the i th site; and γ is the estimated parameter.

Prediction and Calibration
Two situations-using fixed-effects only and a mixture of fixed-and random-effectscan be considered when using a two-level NLME model to predict HCB. The fixed-effect only model can also be called a population average or uncalibrated response model. In contrast, the model that contains random effects is typically called a localized or subjectspecific model, in which the localizing process is mostly known as model calibration [73]. The uncalibrated model nullifies the random effect parameter and does not need any prior information. The subject-specific models were calibrated by predicting the specific plot and site effects using several sampled trees' measured attributes from the validation site. The values of random effect parameters were determined by the best linear unbiased prediction (BLUPs) [74]. The expression is as follows: whereû i = û i ,û i1 ,û i2 , . . . ,û iM i T is a q 1 + Mq 2 -dimensional vector of the estimated random-effects parameters for the i th site;û i is the q 1 -dimensional vector of the estimated value at site level;û i1 toû iM i are the estimated q 2 -dimensional vectors of the randomeffects parameters at the sample plot level;R i is the variance-covariance matrix of withingroup errors; Z i is the design matrix of the partial derivatives of the nonlinear function corresponding to the random parameters; andψ is an (q 1 + Mq 2 ) × (q 1 + Mq 2 ) block diagonal matrix consisting ofψ 1 andψ 2 , the two estimated variance-covariance matrices for the random effects parameters u i and u ij ; and e i is the error terms of the predicted by the fixed effects parameters of the mixed-effects models. The maximum likelihood method of NLME function in NLME package [75] in the software R4.0.3 was employed to estimate the model parameters.

Model Assessment
The prediction applicability of the HCB models including the base model, uncalibrated model, and calibrated NLME model were evaluated by using independent observation data. In this study, all data was utilized for model development, and the leave-one-siteout cross-validation (LOOCV) was used to test the independence and adaptability of the model. Base model, uncalibrated model, and calibrated NLME model were compared using the average statistics of cross-validation within a sample plot. The performance of the model was evaluated by calculating four model validation statistics (Bias, Bias%, MAE, and MAE%) as follows: where HCB t andĤ CB t is the t th observed and predicted height to crown base (t = 1, . . . , N); N is number of the observations; and mean(HCB) is the mean value of HCB observations.

Comparison of Different Sampling Strategies
The mixed-effect model's calibration was calculated using the field-measured HCB of some multiple sample trees as the prior information to predict the specific random-effects parameters. Therefore, we proposed multilevel prediction (both site-and plot-level), which is convenient for application across different scales yet still provides a high accuracy [76].

Site-Level Calibration
The site-level calibration was achieved by setting the plot-level random parameter to zero, which can be completed without utilizing sample plots. Hence it is more suitable for large-scale and efficient prediction. Generally, the larger the calibration sample size, the more accurate the calibration result. However, sampling a large number of sample trees only to calculate the random-effect parameters is impractical. Hence, we proposed a simpler sampling strategy: Selecting 1−50 trees randomly per site. The simulation was repeated 1000 times to calculate the average results, preventing the prediction from being biased.

Plot-Level Calibration
Compared with the site-level calibration, plot-level calibration has higher prediction accuracy since it considers the random effects of both site-level and plot-level nested within the site. Hence, it is more suitable for small-scale and high-precision prediction. Considering the measurement cost and potential error of UAV-LiDAR in extracting forest structure parameters, eight sampling strategies with different subsampling schemes were proposed as follows: Type I: Selecting l trees randomly per plot (l: 2, 3, . . . , 18). Type II: Selecting l largest trees per plot (l: 2, 3, . . . , 18).
Bias% and MAE% were used to assess the prediction accuracy under different sampling strategies and sizes.

Base Model Development
The optimal subset method was utilized to select covariates and H, CW, CCp 75 , and H P 99 were used as predictive variables to expand the basic model since they have the largest value of R 2 and the smallest value of RMSE, Bias, and AIC. The CCp 75 -the ratio of the target tree's crown area to the sum of all crown areas above their 75% relative tree's total height within the sample plot-can commendably reflect the competitive situation [77]. Meanwhile, the elevation of the 99%-point cloud of the subject plot (H P 99 ) was used to describe the stand variation's effect on HCB. The final multivariate model is as follows: where HCB ijk , H ijk , and CW ijk are height to crown base, total tree height, and crown width of the k th tree in the j th plot within the i th site, respectively; β 0 , β 1 , β 2 , and β 3 are model parameters to be determined; and ε ijk is an error term. The parameter estimates of the base model were presented in Section 3.3. All model parameters are significant, and the model fitting has substantially improved after adding covariates ( Table 3). The basic model explained more variations of HCB after extending the site quality and competition indices, which could be useful for the further construction of the two-level mixed effect model. Table 3. Parameter estimates and fitting statistics for the base and NLME model.

Two-Level Nonlinear Mixed-Effects HCB Model
The influences of site and plot variation on HCB were considered in model (12). There were four parameters (β 0 −β 3 ) considered in the base model, yielding 15 different combinations of random-effects parameters. However, only nine combinations of parameter estimates were successfully converged. The model had the smallest AIC value (25,565) and the highest fitting accuracy (R 2 = 0.9424 and RMSE = 1.1282) when β 0 and β 1 were included as the random effect parameters. All parameter estimates can be found in Section 3.3. The specific form is as follows: where: where u 0i and u 1i are the random effects caused by the i th site on β 0 and β 1 , respectively; u 0ij and u 1ij are the random effects caused by the j th sample plot nested in the i th site on β 0 and β 1 , respectively; and other symbols have been described in previous sections. Figure 3 shows the standardized residuals distribution of the base model and calibrated NLME model. The weighted power function significantly stabilized the heteroscedasticity. The values of the power variance functions are listed in Table 3.
where , , and are height to crown base, total tree height, and cro width of the tree in the plot within the site, respectively; , , , and are model parameters to be determined; and is an error term. The parameter e mates of the base model were presented in Section 3.3.
All model parameters are significant, and the model fitting has substantially i proved after adding covariates ( Table 3). The basic model explained more variations HCB after extending the site quality and competition indices, which could be useful the further construction of the two-level mixed effect model.

Two-Level Nonlinear Mixed-Effects HCB Model
The influences of site and plot variation on HCB were considered in model (12). Th were four parameters ( − ) considered in the base model, yielding 15 different com nations of random-effects parameters. However, only nine combinations of parameter timates were successfully converged. The model had the smallest AIC value (25,565) a the highest fitting accuracy (R 2 = 0.9424 and RMSE = 1.1282) when and were cluded as the random effect parameters. All parameter estimates can be found in Sect 3.3. The specific form is as follows: where: and are the random effects caused by the site on and , resp tively; and are the random effects caused by the sample plot nested in site on and , respectively; and other symbols have been described in previo sections. Figure 3 shows the standardized residuals distribution of the base model a calibrated NLME model. The weighted power function significantly stabilized the het oscedasticity. The values of the power variance functions are listed in Table 3.

Parameter Estimates
All parameters estimated in the NLME HCB model were statistically significant ( 0.05). The specific parameter estimation values and fitting statistics of each model

Parameter Estimates
All parameters estimated in the NLME HCB model were statistically significant (p < 0.05). The specific parameter estimation values and fitting statistics of each model are shown in Table 3. The likelihood ratio test (LRT) of the base nonlinear model (12) and the NLME model (13) showed that the mixed effect model was statistically significant (p < 0.0001), indicating that the plot variation had a significant random effect on HCB.

Model Assessment
The prediction ability of the base model, NLME model with fixed parameters only (uncalibrated model), and NLME model with fixed and random parameters (calibrated model) was assessed by the leave-one-site-out cross-validation method, and the results were compared and are demonstrated in Table 4. The calibrated NLME model had the best adaptability and stability performance compared to the two others, with the MAE and MAE% as low as 0.89% and 9.71%, respectively. The base model's performance was worse than the NLME model but slightly better than the uncalibrated model. Figure 4 visualized that the prediction calculated by the calibrated NLME model is more consistent along the trend line Y = X than the base and uncalibrated model. In addition, the base and uncalibrated models show an apparent overestimation in predicting the individual trees with lower HCB, which has been effectively solved by the calibrated NLME model (Figure 4).  Moreover, the prediction accuracy of the three models was scrutinized across 15 diameter classes ( Figure 5). The upper limit exclusion method was used for diameter class integration with 2 cm equal difference (i.e. [5,7), [7,9), ...). The trees with DBH more than 33 cm are classified into 34 diameter classes. The uncalibrated model shows the worst accuracy among all methods with the MAE value ranging from 0.69 to 1.88. Meanwhile, the calibrated NLME model has the highest prediction accuracy with MAE ranging from 0.61 to 1.47, having approximately a 0.09 to 0.51 decrease compared with the uncalibrated model. The MAE% values were the worst for the low diameter classes and were found to be decreased with the increasing DBH. The calibrated NLME model presented the best prediction accuracies across different diameter classes. It is worth mentioning that the base model is generally underestimated in the lower and larger diameter classes but overestimated in the medium diameter classes ( Figure 5). Underestimation and overestimation lead to the "canceling out" of positive and negative deviation, which also explains why the base model has a lower bias (Table 4). Figure 5 also shows that the error range and prediction accuracy of the calibrated NLME model are relatively stable, while the performance of the uncalibrated model is often worse than that of the base model. The calibrated NLME model provides a clear visualization of the whole population's average responses and the changes between different levels. The estimation was consistent across different diameter classes, proving the mixed effect model's robustness in predicting the property's change of variables. Moreover, the prediction accuracy of the three models was scrutinized across 15 diameter classes ( Figure 5). The upper limit exclusion method was used for diameter class integration with 2 cm equal difference (i.e., [5,7), [7,9), ...). The trees with DBH more than 33 cm are classified into 34 diameter classes. The uncalibrated model shows the worst accuracy among all methods with the MAE value ranging from 0.69 to 1.88. Meanwhile, the calibrated NLME model has the highest prediction accuracy with MAE ranging from 0.61 to 1.47, having approximately a 0.09 to 0.51 decrease compared with the uncalibrated model. The MAE% values were the worst for the low diameter classes and were found to be decreased with the increasing DBH. The calibrated NLME model presented the best prediction accuracies across different diameter classes. It is worth mentioning that the base model is generally underestimated in the lower and larger diameter classes but overestimated in the medium diameter classes ( Figure 5). Underestimation and overestimation lead to the "canceling out" of positive and negative deviation, which also explains why the base model has a lower bias (Table 4). Figure 5 also shows that the error range and prediction accuracy of the calibrated NLME model are relatively stable, while the performance of the uncalibrated model is often worse than that of the base model. The calibrated NLME model provides a clear visualization of the whole population's average responses and the changes between different levels. The estimation was consistent across different diameter classes, proving the mixed effect model's robustness in predicting the property's change of variables.
base model is generally underestimated in the lower and larger diameter classes but overestimated in the medium diameter classes ( Figure 5). Underestimation and overestimation lead to the "canceling out" of positive and negative deviation, which also explains why the base model has a lower bias (Table 4). Figure 5 also shows that the error range and prediction accuracy of the calibrated NLME model are relatively stable, while the performance of the uncalibrated model is often worse than that of the base model. The calibrated NLME model provides a clear visualization of the whole population's average responses and the changes between different levels. The estimation was consistent across different diameter classes, proving the mixed effect model's robustness in predicting the property's change of variables.

Site-Level Calibration
We randomly selected 1−50 trees from each site for site-level local calibration using the BLUPs theory, which is repeated continuously for 1000 times to calculate the average of the error statistics (MAE% and Bias%). The results of the calibration response mode are shown in Figure 6. For site-level calibration, only when the subsample is more than five, the prediction effect is better than that of the uncalibrated model. MAE% gradually decreased with the increase of sampling number. The overall trend of Bias% is similar to MAE%, although there are slight fluctuations. When the subsample size is 15, MAE% was reduced to 11.5%. Increasing the number of samples was not significant to improve the model's accuracy but it will increase the cost of measurement. Therefore, there is always be a tradeoff between the cost and accuracy, in which selecting 15 trees randomly for site-level calibration might be recommended as the best compromise between the two.

Site-Level Calibration
We randomly selected 1−50 trees from each site for site-level local calibration using the BLUPs theory, which is repeated continuously for 1000 times to calculate the average of the error statistics (MAE% and Bias%). The results of the calibration response mode are shown in Figure 6. For site-level calibration, only when the subsample is more than five, the prediction effect is better than that of the uncalibrated model. MAE% gradually decreased with the increase of sampling number. The overall trend of Bias% is similar to MAE%, although there are slight fluctuations. When the subsample size is 15, MAE% was reduced to 11.5%. Increasing the number of samples was not significant to improve the model's accuracy but it will increase the cost of measurement. Therefore, there is always be a tradeoff between the cost and accuracy, in which selecting 15 trees randomly for sitelevel calibration might be recommended as the best compromise between the two.

Plot-Level Calibration
Eight BLUPs strategies in local plot-level calibration were proposed to improve the model's accuracy. The calibration prediction results were visualized by two error statistics (Bias% and MAE%) and shown in Figure 7. Apart from type III, all sampling strategies had similar trends, in which the calibration performance improves with the increase of sample size, and the subject-specific predictions obtained higher accuracy than uncali-

Plot-Level Calibration
Eight BLUPs strategies in local plot-level calibration were proposed to improve the model's accuracy. The calibration prediction results were visualized by two error statistics (Bias% and MAE%) and shown in Figure 7. Apart from type III, all sampling strategies had similar trends, in which the calibration performance improves with the increase of sample size, and the subject-specific predictions obtained higher accuracy than uncalibrated predictions, even though only a small number of trees (i.e., two) are used as the basis for the random effect predictions. The type VI sampling strategy obtained the smallest MAE% when the prediction pre-measured sub-sample size was less than eight trees, while similar prediction performance was delivered by both type VI and type IV when the sub-sample size was more than eight trees. The largest MAE% was always acquired by the type III sampling strategy. For the type VI, more stable prediction was obtained after including seven sampling trees, even the results were still outperformed by other sampling strategies. In addition, a paired t-test was used to compare the calibrated NLME predictions using between 6 trees and 18 trees for calculating the random parameters, and the result was found to be not statistically significant. Hence, using six sampling trees in the HCB model calibration might present the most optimum results when the time and cost efficiency is considered. A relatively high prediction accuracy with a low measurement cost might be achieved by sampling the three largest trees and three medium-size trees (type VI) from each sample plot.

Discussion
As an essential individual tree variable, height to crown base (HCB) can effectively reflect the crown size, subsequently affecting the forest's vitality, wood quality, and commercial value of trees [30]. Hence, it is very important to obtain accurate and efficient data of HCB, which has been conveniently facilitated by utilizing remote sensing data. However, due to the tree canopy occlusion and point cloud density limitation, acquiring a stable method to accurately and directly extract the HCB from LiDAR data remains challenging. Therefore, this study was intended to develop and provide an accurate HCB model using a generalized NLME method.
In this study, three commonly used HCB models were evaluated to construct the HCB model based on LiDAR data [65,69]. Tree height ( ) and crown width ( ) were used as the primary predictors since they account for the most variations on HCB and are relatively easy to be accurately extracted [46]. Furthermore, stand quality and interspecific competition were added to the model as potential predictors. In previous forest modeling studies, dominant height ( ) was the most commonly used variable to represent stand quality [78,79]. However, the LiDAR-based relative percentile height can be directly extracted from each sample plot to replace and was not affected by the calculation method; hence, it is more suitable for our model. All combinations of stand and competition metrics were analyzed and used to finally determine both and , which are novel to LiDAR-derived variables. The RMSE decreased by 38.46% after these two covari-

Discussion
As an essential individual tree variable, height to crown base (HCB) can effectively reflect the crown size, subsequently affecting the forest's vitality, wood quality, and commercial value of trees [30]. Hence, it is very important to obtain accurate and efficient data of HCB, which has been conveniently facilitated by utilizing remote sensing data. However, due to the tree canopy occlusion and point cloud density limitation, acquiring a stable method to accurately and directly extract the HCB from LiDAR data remains challenging. Therefore, this study was intended to develop and provide an accurate HCB model using a generalized NLME method.
In this study, three commonly used HCB models were evaluated to construct the HCB model based on LiDAR data [65,69]. Tree height (H) and crown width (CW) were used as the primary predictors since they account for the most variations on HCB and are relatively easy to be accurately extracted [46]. Furthermore, stand quality and interspecific competition were added to the model as potential predictors. In previous forest modeling studies, dominant height (H d ) was the most commonly used variable to represent stand quality [78,79]. However, the LiDAR-based relative percentile height can be directly ex-tracted from each sample plot to replace H d and was not affected by the calculation method; hence, it is more suitable for our model. All combinations of stand and competition metrics were analyzed and used to finally determine both H P 99 and CCp 75 , which are novel to LiDAR-derived variables. The RMSE decreased by 38.46% after these two covariates were introduced into the model, showing that the covariates extension was effective and could be further used for the random effects. The negative parameters of both stand and competition index indicated that the HCB was positively correlated with the stand quality and competition size. This is ecologically reasonable since the better the stand conditions, the higher the tree height, consequently promoting the crown's overall upward movement and increasing the HCB. The competition index expressed by CCp 75 was constantly recognized in the description of the crown competition. The larger the CCp 75 of the targeted trees, the weaker the competition ability, and the more difficult for the trees to acquire sunlight, water, and soil nutrients, resulting in a shorter and smaller crown [80,81]. The proposed multivariable generalized base model has high adaptability and was found to be effective for detecting the variation of HCB with the same height and crown width. The selected experimental sites have similar geographical conditions; hence, we did not consider the influence of altitude, terrain, and slope on the developed model.
The nonlinear mixed effect (NLME) model can reflect the potential variations among different levels and has been widely utilized to construct the HCB model [45,47,66]. In this study, the data was sampled from an extensive and widely distributed area. Hence, applying a two-level nested random effect model will be appropriate to describe the specific effect of the plot and site variations on HCB. Generally, the mixed effect model's convergence becomes harder to be achieved as the number of random effect parameters increases, especially when the model contains more than two parameters [82]. In this study, all possible combinations of random effect parameters were considered. The HCB model obtained the smallest AIC and the highest R 2 after the random effects were added into the intercept β 0 and the regression coefficient β 1 of the crown width. The calibrated model's accuracy was significantly higher than that of the generalized base model, which was indicated by LRT. In addition, three functions were compared for selecting the appropriate weighting function, in which the power function was finally chosen due to its simplicity and convenient application.
The leave-one-site-out cross-validation method was applied to assess the independence and applicability of the models. As shown in Table 4, the introduction of random effects has greatly improved the model's prediction accuracy, even after it was scrutinized throughout various diameter classes ( Figure 5). The relative bias of the calibrated NLME model was closer to zero and relatively more stable than the other two models. However, the model generally produced a larger MAE% and Bias% in small diameter classes (i.e., 6−8 diameter classes). The possible reason is that the younger trees often have no crown overlap between individuals, leading to less competition among trees. Thus, the lower canopy's living branches have more opportunities to absorb sunlight for photosynthesis, which contributes to enhancing their growth activity. Furthermore, HCB is usually more affected by non-natural factors (e.g., tending and pruning activity) [33]; hence, various management measures should be considered to accurately predict the HCB of young trees.
Compared with the recently developed HCB models which are only based on fieldmeasured attributes [45,81], our LiDAR-based models might provide a more highly efficient prediction that requires less effort to obtain the input information. Several researches [41,42,83] extracted HCB directly from LiDAR data. However, the algorithms are not applicable for large-scale estimation. In addition, some studies [46,52] include field-measured variables (i.e., DBH) in building LiDAR-based HCB models that was unnecessary in our study, which improves the flexibility of the model.
We proposed both site-and plot-level model calibration to improve the flexibility of model application throughout different scale and accuracy requirements. Site-level calibration does not need additional plot information; hence, it is convenient for efficient prediction in large areas. Therefore, we proposed a simple sampling strategy (1−50 trees randomly selected per site). The model's prediction accuracy will improve as the number of sampled trees per plot increases [83]. The results showed that randomly selecting 15 trees from each site may be a good option when both precision and cost are considered. In addition, for plot-level calibration, eight complex sampling strategies and various numbers of subsampled trees were assessed by leave-one-site-out cross-validation, aiming to ascertain the best calibration scheme for predicting the HCB of new stands. The results showed that using the smallest trees for model calibration (type III) obtained the worst result, consistent with other studies [45,46,51]. This may be attributed to the fact that some fixed-effect variables in the NLME HCB model have reflected the variation of small-size trees; hence, the type III sampling strategy did not provide any additional information for calibration. On the contrary, using the largest trees and medium-size trees (type VI sampling strategy) obtained the best calibration results. It is worth recalling that type IV was worse than type VI when the sample size was small, but the difference between them gradually decreased as the sampling size increased. However, this correction strategy is necessary to obtain the DBH of each tree in the sample plot, which undoubtedly increases the workload of the field measurement. The result of random sampling (type I) is also given, which is only slightly worse than type VI (Figure 7). It can be used when the field-measured DBH is not available. Generally, four to nine sample trees are used for mixed effect model plot-level calibration to ensure a tradeoff between the model predictability and inventory cost [79,84]. This study suggests that using six trees for NLME model calibration is appropriate when prediction accuracy and measurement cost are considered. Adding more sample trees in model calibration seems inefficient since it brings an insignificant improvement to the model performance and yet causes a remarkable increase in inventory cost.
This study only used the correctly matched individual trees for constructing and evaluating the model. However, the commission and omission errors cannot be ignored since they might affect the method's actual application. The segmentation errors often increase with stand density, which mainly miss detection of suppressed trees, especially in young forests [85]. Therefore, a careful application should be conducted when applying this method in high-density young forests. In the future, a segmentation algorithm suitable for the high-density forest will substantially improve the application of this model.

Conclusions
Larix olgensis is one of the main afforestation tree species in Northeast China, which has fast growth and provides high timber output. In this paper, a generalized nonlinear mixed effect HCB model of larch plantation was established using UAV-LiDAR data. The newly developed model could complement the missing HCB data in forest inventory and reduce the field workload without ignoring the prediction accuracy. The tree-and plot-level LiDAR-derived metrics (i.e., H, CW, CCp 75 and H P 99 ) were introduced into the base model, which significantly improved the model's ability to explain the HCB variation. The introduction of two-level random effects greatly improved the application ability and prediction accuracy of the model. The newly developed model not only has an important significance for HCB prediction, but also supplies a more flexible and convenient method for forest application based on UAVs.