Individual Tree Diameter Estimation in Small-Scale Forest Inventory Using UAV Laser Scanning

: Unmanned aerial vehicle laser scanning (UAVLS) systems present a relatively new means of remote sensing and are increasingly applied in the ﬁeld of forest ecology and management. However, one of the most essential parameters in forest inventory, tree diameter at breast height (DBH), cannot be directly extracted from aerial point cloud data due to the limitations of scanning angle and canopy obstruction. Therefore, in this study DBH-UAVLS point cloud estimation models were established using a generalized nonlinear mixed-effects (NLME) model. The experiments were conducted using Larix olgensis as the subject species, and a total of 8364 correctly delineated trees from UAVLS data within 118 plots across 11 sites were used for DBH modeling. Both tree- and plot-level metrics were obtained using light detection and ranging (LiDAR) and were used as the models’ independent predictors. The results indicated that the addition of site-level random effects signiﬁcantly improved the model ﬁtting. Compared with nonparametric modeling approaches (random forest and k-nearest neighbors) and uni- or multivariable weighted nonlinear least square regression through leave-one-site-out cross-validation, the NLME model with local calibration achieved the lowest root mean square error (RMSE) values (1.94 cm) and the most stable prediction across different sites. Using the site in a random-effects model improved the transferability of LiDAR-based DBH estimation. The best linear unbiased predictor (BLUP), used to conduct local model calibration, led to an improvement in the models’ performance as the number of ﬁeld measurements increased. The research provides a baseline for unmanned aerial vehicle (UAV) small-scale forest inventories and might be a reasonable alternative for operational forestry.


Introduction
Forests, as one of the essential elements of a terrestrial ecosystem, have a crucial role in terms of regulating fluxes and stores of carbon and water, contributing to biodiversity conservation, and regulating the global climate system [1,2]. To quantitatively assess the amount and map the distribution of forest and its changes, a timely and accurate forest resource inventory is needed [2].
Traditionally, forest inventories depend on the sampling of the ground truth (in situ measurements), where each selected individual's attributes are obtained through tree-by-tree measurements [3,4]. Such inventories are not cost-effective since the field measurements are often labor-intensive and time-consuming, consequently limiting the sampling intensity and number of tree attributes measured [4]. Developments in remote sensing technologies have brought about massive breakthroughs in terms of improving the performance of forest inventory, specifically with respect to the measurement scale and efficiency [5]. Airborne LiDAR is capable of reconstructing the detailed 3D structure of Remote Sens. 2021, 13, 24 3 of 21 though nonparametric approaches involve fewer assumptions and could achieve a higher prediction accuracy, they are not a suitable tool to understand the relationship between LiDAR metrics and forest inventory attributes [36].
One of the inherent problems with using remote sensing data as predictors is the calibration and transferability of the established equations [37,38]. Models generally perform well in regions where the modeling data are located [39]. However, practical inventories are mostly performed using different sensors or imaging parameters, and the structure of the stand-and tree-level attributes can be remarkably varied between sites [38]. Different laser scanning devices, modes, and operational parameters may generate notably distinctive point clouds [40][41][42]. Point cloud-based estimations produced by a particular device or parameter might have a higher accuracy in some sites than in others [39]. Thus, the utilization of established equations often leads to bias and is less accurate; on some occasions, the decrease is trivial and still acceptable [39]. Furthermore, forest structure can also be a source of transferability problems, mainly if there is a considerable difference between the source and the target locations, or if the model fitting data do not properly cover the variations in the data used for model validation [38,39].
In UAVLS-based forest inventories, more attention should be paid to the transferability of the model, since UAVLS acquisition is primarily affected by several environmental factors. Moreover, UAVLS is more suitable for inventorying small forest areas with a high accuracy, and it is practical to use the existing models and the new UAVLS acquisition for predicting certain variables. Our study thus focuses on developing a DBH estimation model for UAVLS-based forest inventories. The aims of the present research were to (1) delineate individual trees and extract various tree-and plot-level metrics derived from UAVLS data, (2) present a nonlinear mixed-effects (NLME) modeling framework for estimating individual tree DBH using UAVLS-derived metrics and random effects, (3) calibrate and assess the established NLME model using the field data from a different region, and (4) compare the accuracies and transferability of the DBH estimation with alternative modeling approaches. This study could provide a baseline and new perspective for DBH modeling using UAVLS for small forest inventories.

Study Area and UAVLS Data Acquisition
This research was conducted in Mengjiagang Forest Farm (130 • 32 0 -130 • 52 6 E, 46 • 20 20 -45 • 30 16 N), which is located along the western fringes of Wanda mountain in the northeastern part of Huanan County, Heilongjiang province, China. The terrain mainly features gentle slopes and low-elevation hills, with an average elevation of approximately 250 m a.s.l. [43]. The forest vegetation in the area is dominated by coniferous plantations, predominantly Pinus sylvestris var. mongolica, Pinus koraiensis, Larix olgensis, and Picea asperata.
In the present study, we used larch, one of the most abundant and economically important species in northeast China, for DBH modeling. A total of 11 sites of larch plantations of different age groups, stand densities, and forest conditions were selected to represent all larch stands in the study area ( Figure 1). All the sites were dominantly covered by Larix olgensis plantations with an initial planting spacing of 2 × 1.5 m (3300 stems/ha). The stand densities were adjusted by thinning practice depending on the growth stage and site quality, which were about 3000, 2000, 1000, 700, and 500 stems/ha for young, middle-aged, near-mature, mature, and over-mature forests, respectively.
were selected to represent all larch stands in the study area ( Figure 1). All the si tes were dominantly covered by Larix olgensis plantations with an initial planting spacing of 2 × 1.5 m (3300 stems/ha). The stand densities were adjusted by thin ning practice depending on the growth stage and site quality, which were about 3000, 2000, 1000, 700, and 500 stems/ha for young, middle-aged, near-mature, mat ure, and over-mature forests, respectively. UAVLS data were acquired for these 11 sites from 10 to 12 July 2019, using a RIEGL mini VUX-1UAV LiDAR scanner (www.riegl.com/products/unmanned-sca nning/riegl-minivux-1uav) carried by a Feima D200 UAV platform. The LiDAR sen sor operated at a 100 kHz pulse repetition rate with a maximum scan speed of 10 0 scans per second. The scanner's maximum measurement range was 250 m with a n accuracy of 15 mm. All the flights were designed as crossed transects with 80 m swath overlaps at 80 m altitude and 5.0 m/s speed, conditions which are common ly used in UAVLS data collection [13,15]. The sensor provides a 360° field of view. In interpreting the point cloud data, the scanning angle was specifically set to diff erent sites depending on the trajectory and overlaps, but all were within ±60° to av oid measurement errors caused by excessive angle. The average pulse densities for each site ranged from about 120 to 220 pulses/m 2 , and the final point densities wer e about 150 to 270 pt/m 2 , with up to five echoes. The descriptions of forest charact eristics and UAVLS point densities are listed in Table 1 for each site. UAVLS data were acquired for these 11 sites from 10 to 12 July 2019, using a RIEGL mini VUX-1UAV LiDAR scanner (www.riegl.com/products/unmanned-scanning/rieglminivux-1uav) carried by a Feima D200 UAV platform. The LiDAR sensor operated at a 100 kHz pulse repetition rate with a maximum scan speed of 100 scans per second. The scanner's maximum measurement range was 250 m with an accuracy of 15 mm. All the flights were designed as crossed transects with 80 m swath overlaps at 80 m altitude and 5.0 m/s speed, conditions which are commonly used in UAVLS data collection [13,15]. The sensor provides a 360 • field of view. In interpreting the point cloud data, the scanning angle was specifically set to different sites depending on the trajectory and overlaps, but all were within ±60 • to avoid measurement errors caused by excessive angle. The average pulse densities for each site ranged from about 120 to 220 pulses/m 2 , and the final point densities were about 150 to 270 pt/m 2 , with up to five echoes. The descriptions of forest characteristics and UAVLS point densities are listed in Table 1 for each site.

In Situ Measurements
A total of 118 plots (30 m × 30 m) across 11 experimental sites were established in July 2019. The sample plots were evenly distributed in each site, excluding the edge of the stand and large forest gaps. The distance between plots was larger than twice each site's average tree height. The number of plots for each site is listed in Table 1. DBH and crown width measurements were conducted for all trees with DBH more than 5 cm using a diameter tape. Meanwhile, the height of the first live branch position (at crown base) and the total tree height were measured using a Vertex IV instrument with a height resolution of 0.1 m (Haglöfs, Sweden). The trees' coordinates were recorded by measuring their relative position to the corner of plots. Furthermore, the geographical coordinates of the individual trees and four corners of each plot were determined with a real-time kinetic (RTK) global navigation satellite system (GNSS) (UniStrong G10A, Beijing, China) with a positioning error of approximately 0.1 m. Trees under poor GNSS signals were georeferenced by their relative coordinates, the positioning accuracy of which was estimated to be about 0.3-0.5 m.

LiDAR Metrics Extraction
In this study, both tree-and plot-level metrics were extracted from the UAVLS data. The tree-level metrics include the basic characteristics of each individual and their competitive status within the sample plot. In addition, stand conditions also affect individual stem diameter growth and size distribution. The plot-level metrics were used to describe the crown structural and topography characteristics that were introduced as auxiliary information.

UAVLS Data Preprocessing
Firstly, the noise points were manually removed from the raw LiDAR point clouds. The cloth simulation filtering was utilized to categorize the remaining points into nonground and ground points [44]. Due to the dense canopy cover, the ground point cloud density is about 10 pt/m 2 on average. Then, the Kriging method with a 0.5 m pixel size was used to interpolate the ground points into digital terrain models (DTMs) [45]. DTM values were subtracted to obtain the normalized height of each point [46,47]. As mentioned in early studies, data pits in the canopy height models (CHMs) disrupt the crowns' integrity and smoothness (Figure a), which negatively affects the individual crown delineation and parameter extraction [47,48]. We therefore applied graph-based progressive morphological filtering (GPMF), a canopy filtering technique, to generate CHMs from the UAVLS data [47]. This algorithm employs an adaptive morphological operation to eliminate depression points from all returns in progressive filtering. All the remaining surface points were then interpolated by triangulated irregular networks. Many studies have generated CHM from UAVLS data with various spatial resolutions, ranging 0.1 to 1 m [10,16,49,50]. Yin and Wang (2019) have recommended that spatial resolution should be finer than one fourth of the crown diameter to correctly delineate crown boundaries and characterize the crown shapes. The grid cell should also not be much smaller than the average pulse spacing [51]. In our study, the CHMs were therefore generated at a resolution of 0.1 m, which could recognize the minimum crown diameter of 0.6 m recorded in the field survey ( Table 1) and was sufficiently supported by the lowest pulse spacing of 0.09 m. As shown in Figure 2b, canopy surfaces could be characterized clearly and with few data pits.

Individual Tree Delineation
In recent years, two main categories of individual tree delineation methods (point cloud-and CHM-based) have been developed for ALS data [9]. Although point cloudbased algorithms are capable of capturing the understory trees under main canopies, a huge amount of computation and complex parameters limit the application's efficiency, especially for high-density UAVLS point clouds. CHM-based individual tree delineation methods were applied in this study since, in a larch plantation, there are only a few understory trees. Zhao et al. (2017) constructed a single-tree automatic detection algorithm called region-based hierarchical cross-section analysis [52]. This particular algorithm used horizontal relationships between the crown within the vertical direction, and the CHM was considered as a mountain-like topographic surface for detecting an individual tree. In order to avoid the influence of shrubs in the segmentation process, the height of the crown for automatic segmentation was limited to 3 m.

Tree-to-Tree Matching
All the sample trees that met the modeling's requirements were screened using a tree-to-tree matching procedure, harmonizing the field-measured data and the segmented trees based on spatial locations and height differences [53]. The segmented trees were assigned as candidates of a reference tree if the horizontal distance from the reference was less than the corresponding crown radius (with an upper limit of 3 m); meanwhile, the height difference from the reference was less than 20% of the top height of the plot [54]. A unique candidate or the closest one among multiple candidates was selected as a match with the reference tree [47]. Segmented trees without a link to references were considered commission errors ( Figure 3). Conversely, reference trees that were not matched to any segmented trees were classified as omission errors ( Figure 3). In total, 8785 trees were correctly matched with the field measurements; 56-100% (mean 76%) of the trees were detected among a total of 118 plots with commission errors of 17-45% (mean 28%). As shown in Figure 3, the relatively high commission errors in our study were mainly due to the trees outside the plot, but their crowns extending into plot boundaries were mistaken as individual trees. After the matching was completed, all irrelevant trees, such as dead trees, miscellaneous trees, and incomplete segmented crowns along plot boundaries, were manually removed from the matching tree datasets. Overall, 8364 trees were selected for further modeling.

Tree-and Plot-Level Metrics Generation
In this study, a total of 100 UAVLS-based tree-and plot-level variables were calculated as potential predictors for future modeling. The laser returns falling within each segment were masked to clip out an individual tree and utilized to determine a set of tree-level metrics. The total height (H) of each individual tree was defined as the maximum height of all LiDAR pulses. The relative root mean square error between the LiDAR-derived and fieldmeasured height was about 7% (Supplementary Materials Figure S1). A 2D convex hull of projected individuals' points on the x-y plane was constructed as a crown polygon [55]. Crown diameter (CD) was calculated by 2 × √ crown area/π, where crown area was the area of the individual convex hull [56]. Moreover, the crown volume and surface area [26,47] and a series of point distribution and density metrics were also applied as candidates to characterize the structure of the detected trees [35,38,57].
The DBH growth of a tree at a particular site is generally influenced by its competitive status. A crown area competition-based measure was evaluated at a certain percentage of crown length to be expanded and calculated using LiDAR data [58]. In this study, the ratio between the crown area and the computed reference height, equal to p% (25%, 50%, 75%, and 100%) of the tree's height and the total crown area, was calculated as a competition index and marked as CC p25 , CC p50 , CC p75 , CC p100 (see Supplementary  Material Table S1 for a detailed description). Other relative dimensions of the crown projection area and tree height derived from the LiDAR data were also generated as potential competitive indicators.
The plot-level metrics were extracted and coded using FUSION and MATLAB 2016a, respectively. The canopy surface height distribution (canopy rugosity (Cnp R ), mean canopy height (Cnp H ), and canopy openness (Cnp O )) was imputed from pit-free CHMs to present the canopy structural heterogeneity [59]. The commonly used point cloud metrics such as height and distribution statistics were also calculated as candidates [35]. Additionally, nine topographic parameters were applied to reflect the variations in topographic conditions [24,30]. A detailed description of all tree-and plot-level metrics, along with the corresponding references in the literature, is given in Supplementary Material Table S1.

Base Model Selection
The most common approach to predicting diameter from LiDAR data is using LiDARderived tree height (H) as a sole predictor or adding delineated crown attributes as other predictors [60,61]. Thus, we first fitted the DBH-H function as the base model. Six candidate Remote Sens. 2021, 13, 24 8 of 21 equation forms were evaluated and the combination of power and exponential function (Equation (1)) was chosen as the base model for future NLME modeling. The detailed selection and evaluation of the candidate base models are presented in Supplementary Material Figure S2. A logical constraint was incorporated to ensure a zero DBH when the tree height equals the breast height (1.3 m).
Here, D ij and H ij represent the DBH and LiDAR-derived tree heights, respectively, for the jth tree in the ith site. a, b, and c are the parameters to be estimated. ε ij is the error term. Due to the residual variance increasing with respect to the prediction, a power-type variance function (Equation (2)) was applied for correcting the variance heterogeneity [62,63]: where σ 2 is a scaling factor for error dispersion and δ is the estimated parameter in model fitting. The model was fitted using all observations (8364 correctly detected individual trees) and nonlinear weighted least-square regressions in the R software (www.r-project.org).

Extension of a Base Model
Then, the selected base model was expanded as a generalized model through the inclusion of various covariate predictors. Besides tree height, DBH is also influenced by the tree's size, competition status, and stand characteristics [24,25,63]. We assessed the influence of other LiDAR-derived parameters on DBH by a two-step covariate selection approach for the generalized DBH estimation modeling. First, each plot data was fitted using the selected base model (Equation (1)) to obtain the corresponding parameter estimation values [31,64,65]. The relationships between model coefficients and the extracted LiDAR metrics (see Section 2.3.4) and their logarithmic transformations were then scrutinized by graphical and correlation analyses [64,66]. As with many studies of LiDAR-derived DBH modeling, crown diameter (CD) was introduced as a predictor to explain the DBH size variation under the same tree height. In addition, the correlation analysis indicated that the competition and site condition had a relatively large impact on the model parameters.
To stabilize the overparameterization and collinearity effects, we selected one competition index (CC p75 ) and one plot-level metric (Cnp R ) as predictors to construct the generalized DBH estimation model. These three variables (CD, CC p75 , and Cnp R ) were applied for the parameterization of a in the base model as a = f CD, CC p75 , Cnp R . The generalized model can be expanded and rewritten as Equation (3): where CD ij is the crown diameter of jth tree in the ith site; Cnp Rij is the canopy rugosity of the plot where the jth tree in the ith site is located; and CC p75ij is the competition index of the jth tree in the ith site. a 0 , a 1 , a 2 , and a 3 are the parameters to be estimated; other parameters were the same as those defined in Equation (1). This expanded generalized model was also fitted by weighted nonlinear least-square regressions and with the variance function in Equation (2).

Nonlinear Mixed-Effects Modeling
The entire dataset contained 11 UAVLS sampling sites; the tree-level attributes were nested in the site-level ones, which could be considered longitudinal data for DBH modeling. We therefore set the sample sites as random effects to account for the DBH variation. For the six parameters (a 0 -a 3 , b, c) in Equation (3), there are 63 combinations to construct a structural model with site-level random effects. Considering the largest logarithm likelihood values (LL) and the smallest Akaike information criterion (AIC) among the converged Remote Sens. 2021, 13, 24 9 of 21 models, the generalized model (Equation (3)) was expanded by adding in a 0 , a 1 , and a 3 to represent the site-level random effects. The final NLME model was: where a 0 -a 3 , b, c are the fixed-effects parameters and u i0 , u i1 , and u i3 are the site-level random-effects parameters. u i is the random-effects parameter vector of the ith site, which is assumed to be normally distributed with zero mean and an unstructured variancecovariance matrix of Ψ, where σ 2 u i0 and σ u i0 ,u i1 represent the variance-covariance components of the site-level random effects [62,67]. ε i is the within-group error, following a normal distribution with an average value vector of zero and a variance-covariance matrix R i [32]. σ 2 is the scaling factor for residual dispersion common to all sites; G i is a diagonal matrix of within-sample site heteroskedasticity variances in which diagonal elements were provided by the variance function Equation (2); Γ i is simplified as an identity matrix, considering no correlation patterns within the same sample site [63]. The NLME models were fitted to all observations using the NLME package [68] in the R environment ( www.r-project.org) by the method of restricted maximum likelihood (REML).

Prediction and Calibration of the NLME Model
In the prediction phase of the mixed-effects model, two different situations-fixedeffects or a combination of fixed and random effects-can be considered for DBH estimation [69,70]. The model without estimated random effects is known as the mean response or uncalibrated prediction. Conversely, the model with estimated random effects from the measurement of a response variable is known as a subject-specific or calibrated prediction [69,71,72].

•
Prediction of mean response: For the prediction of mean response, there is no need to conduct new in situ field measurements of response variables. The prediction only utilizes fixed values of the mixed-effects model [63,64].

•
Prediction with local calibration: For subject-specific prediction, the tree-level attributes measured from validation sites were used to predict the site effects in the model calibration process. The best linear unbiased predictions (BLUPs) method was then used to calculate the parameter of random effects [73]. A vector of random effects parameters of sampled plot i was calculated with Equation (5): whereû i is a vector of the random-effects of the ith sampled site,Ψ is the estimated variance-covariance matrix for the random effects,R i is the variance-covariance matrix of within-group errors, Z i is the design matrix of the partial derivatives of the nonlinear function corresponding to the random parameters, and Z T i represents the transposition of Z i . e i is the dimensional error terms of the ith sampled site predicted by the fixed effects parameters of the mixed-effects model.
Generally, the more measured trees are used for estimating the random-effects parameters, the higher the prediction accuracy [70]. Thus, we had two strategies for resampling calibration, considering both measurement costs and estimation accuracy: (1) Random selection of 1-50 individual sample trees across a validation site.
(2) Random selection of 1-5 square subsample plots with various sizes (length of 5-30 m) within a validation site. Furthermore, all trees located in the subsample plots were selected for calibration.

Benchmarking with Nonparametric Models
Nonparametric modeling methods have been used to estimate forest attributes from LiDAR-derived data which involve fewer assumptions and could achieve higher prediction accuracy [36,74]. In the present study, the two most common nonparametric modeling methods (random forest and k-Nearest Neighbors) were applied for further comparison with the NLME models.

Random Forest
Random forest (RF) is a technique that creates a set of decision trees and then aggregates the results for classification and regression [75]. Each tree is generated independently using bootstrap samples from the training dataset called "in-bag samples" (usually twothirds of the data). Meanwhile, the remaining "out-of-bag" (OOB) samples are used for internal cross-validation. The relative importance of each metric was ranked by quantifying the mean square error increase when each variable of the OOB samples is randomly permuted [75]. Depending on the variable importance ranking, an iterative backward elimination procedure is used for stepwise variable selection [76]. All the variables were first added in an RF, and the less important variable was eliminated; the variable importance was then recalculated using the remaining variables. This procedure was iteratively repeated until a given number of variables was obtained. Herein, 15 predictors were selected from all extracted variables for RF modeling since some additional variables could not significantly decrease the OOB error (mean square error of the OOB samples); see Supplementary Materials Figure S3. Two parameters in the RF modeling, n tree and m try , were set as 400 and 1/3, respectively.

k-Nearest Neighbors
The k-Nearest Neighbors (k-NN) algorithm is a nonparametric estimation method based on the statistical difference between the predictor values and the reference samples. The nearest neighbors used for prediction have been widely applied and discussed for forest applications [77,78]. Variable selection strategies and parameter optimization methods for nearest neighbor imputation have been discussed in detail in [76,78]. In this study, the variables obtained from the RF method were also applied for the k-NN method, as in [76]. For the k-NN computation, the most similar neighbor (MSN) distance metric was used with a canonical correlation analysis based weighting matrix to choose the most similar neighbors [79]. The neighbor number was set to five.

Model Evaluation and Validation
The extrapolation and transferability of the DBH estimation (via the base, generalized, NLME, RF, and k-NN models) were assessed using the observed data from independent regions (outside the scope of the modeling data). The leave-one-out cross-validation (LOOCV) was employed to avoid overestimations [80,81]. In particular, we adapted the site-level LOOCV method (named leave-one-site-out cross-validations) instead of the commonly used tree-and plot-level. It was run by iteratively leaving out one site (N − 1) from the full dataset, aiming to simulate the estimations' expansibility bias across different sites [39]. The mean error (BIAS, in cm) and root mean square error (RMSE, in cm) were computed using the predicted and observed DBH in the site-level LOOCV as follows: where D ij is the observed DBH value of the jth tree in the ith sample site;D ij,−i is the predicted value of the model, which was fitted using all observations without the ith sample site; m is the number of sample sites; and n i is the number of observations in the ith sample sites. Simultaneously, relative BIAS and RMSE (BIAS% and RMSE%) were also applied for evaluation. Both estimation types (with and without local calibration) were included for validating the NLME prediction. The single-or multivariable (base and generalized model in the NLME modeling) weighted nonlinear least square regression and nonparametric methods (RF and k-NN) were compared with the NLME model for DBH estimation.

Model Fitting
The parameter estimates and fitting performances of the one-variable base model (Equation (1)), multivariable generalized model (Equation (3)) and NLME model (Equation (4)) are presented in Table 2. The base model described about 81% of variation while applying tree height as a sole predictor. After adding covariates (CD, CC p75 , and Cnp R ) into the univariable base model, there was a considerable enhancement in model fitting; the RMSE decreased by about 30% (from 2.6397 to 1.8926) and the R 2 a increased by about 10% (from 0.8105 to 0.9026). Meanwhile, both the AIC and LL in Equation (3) had a 14% decrease and increase, respectively. These results show that the fitting performance of Equation (3) is more generalized than that of the base model. Furthermore, there was a further improvement in the model fitting after introducing site random effect parameters of u 0 , u 1 , and u 5 ; the NLME model achieved a higher R 2 a and LL and a lower RMSE and AIC than the generalized model. In addition, the result of LRT between the generalized and NLME models was statistically significant (p < 0.0001), which implied that there were significant site-level random effects on the variation in DBH.

Different Calibration for NLME Model
We applied two strategies (subplots' and trees' random sampling) to locally calibrate the NLME model and calculate the site-level random effects for predicting DBH. Both of the calibration methods were repeated continuously 1000 times. The average RMSE is presented in Figure 4. For each method, the subject-specific prediction with local calibrated random effects could achieve lower RMSE than the uncalibrated NLME model (mean response prediction by the fixed parameters presented in Table 2). The RMSE values decreased with the increasing number of sampling blocks and the width of blocks. In addition, the RMSE also decreased as the number of trees increased. Overall, the prediction accuracy has a positive correlation with the sampling number, indicating that more trees being used for estimating the random-effects parameters will yield a higher prediction accuracy. Considering the time and cost of the in situ measurements, we herein applied the feasible scheme of sampling 20 trees at a particular calibration site as a baseline for NLME prediction to compare with other methods.

Comparison of Prediction
The results of the six DBH estimators are shown in Table 3 and Figure 5. The base univariable model with the LiDAR-derived tree height data as the only predictor had a relatively poorer performance than the other multivariable equations. This indicated that using tree height as a sole predictor is insufficient to explain the variation in DBH. The NLME prediction calibrated with 20 sampling trees exhibited the lowest BIAS (0.02 cm) and RMSE (1.86 cm). Compared to the generalized model and the uncalibrated NLME model, the insertion of site-level random effects led to an improvement in the prediction accuracy. The generalized and uncalibrated model performed slightly better than the two nonparametric models (RF and k-NN). The accuracy assessments based on leave-onesite-out cross-validation showed the transferability of the multivariable parametric model across different sites. Figure 5 presents the relationship between the predicted and the fieldmeasured diameter at breast height (DBH) calculated with six different methods, which shows that the calibrated NLME method has more compact and uniformly distributed points on both sides of the y = x trend line.  The predictions across 11 different sites were separately assessed, and the RMSE values are given in Table 4. The univariable base model always had the worst accuracy at every site, while the multivariable parametric and nonparametric models achieved more reliable results. Meanwhile, the NLME model with local calibration always achieved the highest prediction accuracy, with the RMSE values ranging from 1.64 to 2.37. Using 20 sample trees per site for local calibration might decrease the RMSE value by 0.02-0.41 with respect to the uncalibrated prediction. The six models' residuals were plotted across 11 sites in Figure 6. The error range of the calibrated NLME model was relatively stable compared with other methods across all 11 sites, which indicated a more stable prediction accuracy and better transferability between different sites. The calibrated NLME model indicates the robustness of the mixed-effects modeling technique, which is suitable when the predictor properties are highly varied. The mixed-effects model comprises both fixedand random-effects parameters, which can express not only the mean response of the whole population but also the variation between individuals.  In addition, the prediction results of the six approaches were assessed across different forest structures (age groups). The RMSE and RMSE% values are plotted in Figure 7. The RMSE% values were the worst in young stands and decreased with the increase in age for all methods. The univariable base model always exhibited the worst performance, while the multivariable parametric and nonparametric models had substantial improvements in the RMSE and RMSE% values in each age group. Furthermore, the calibrated NLME model almost presented the best prediction accuracies across different age groups. The calibrated NLME model decreased the RMSE and RMSE% values by 1-8% compared with the uncalibrated prediction and achieved the smallest improvement in young stands.

Discussion
In recent years, miniaturized laser scanning sensors on unmanned aerial vehicle platforms have gained an outstanding reputation in practical forest inventory [10,82,83]. Using remotely sensed data for individual tree DBH estimation is an appealing prospect, providing a new way to quantify the stand's volume, biomass, and carbon stock and to reconstruct other forest attributes (e.g., stem diameter distributions [18][19][20]). To address the issue of DBH prediction from UAVLS data, herein we utilized a multi-echo RIEGL UVX1-mini sensor to scan 11 plantation sites of larch, one of the most abundant and economically important species in northeast China. We applied the NLME modeling approach to establish the generalized individual tree DBH equation. The accuracy of the developed equation was evaluated and compared with the weighted nonlinear least squares model and nonparametric regressions through leave-one-site-out cross-validation.
LiDAR-derived tree height data are often utilized to estimate DBH because of the robustness of the LiDAR-derived height measurements and the strong relationship be-tween DBH and tree total height [21,22,63]. In this study, we also used tree height as a basic predictor for NLME modeling. The covariates and random effects were then introduced for model generalization. The results show that the base univariable model performed less well than the multivariable model. This finding corroborates previous research conducted by Jucker et al. (2017), which reported difficulties in using LiDAR-derived tree height data as a sole predictor to reflect the DBH variation [22]. Although some studies have achieved relatively a good performance and high accuracy, the diameter ranges and species information should be taken into account in developing the model for a proper comparison [21,60]. On most occasions, each species has a maximum growth limit of the tree height, and the relationship between DBH and tree height across various ranges of tree size is generally nonlinear [27,63]. The inclusion of crown attributes is essential in order to differentiate trees of a similar height with substantially different diameter sizes. Several crown attributes, such as crown project area, crown surface area, and crown volume, have been reported to have a significant influence on increasing the prediction accuracy of DBH [28,29,63]. In the present study, adding these crown attributes to the model did not lead to any significant improvement compared to using crown diameter.
The diameter growth of individual trees is also affected by the growing environment and stand conditions [30,31]. Hence, our study introduced the LiDAR-derived competition index (CC p75 ) and plot-level metrics (Cnp R ) into DBH modeling. CC p75 represents the competition status, which has a negative effect on DBH growth. Cnp R is the canopy rugosity within each plot, which could be considered as the explanation of the stand density and reflects the variation in the canopy surface's height; the smaller the overlaps between crowns, the greater the variation. The proposed multivariable generalized model shows a high flexibility, which could further explain the DBH variation under similar height and crown sizes. In previous studies, Paris and Bruzzone (2019) also took topography metrics to improve the DBH estimation [24]. However, due to the small terrain difference in our study area, no particular topographical metrics were incorporated in the generalized model. The features selected for stepwise selection using the random forest nonparametric are the height, crown, and competition metrics of the individual tree-and plot-level canopy height metrics (Supplementary Material Figure S3), which emphasizes the importance of the abovementioned variables in describing the DBH variation. With respect to the comparison of predictions among different age groups, all the methods performed the worst in the young stand. This might be because height measurements were taken as the main predictor for DBH estimation, but DBH variation may not be explained well by tree height variation in young stands compared with other stands (see the RMSE% values of the base model in Figure 7). The DBH of larch grows fast at a young stage, so factors reflecting the tree's vigor (such leaf area, site index) may be helpful and should be tried in studies focusing on young stands.
A lot of previous studies have applied nonparametric models to predict DBH utilizing LiDAR-derived variables [34,35]. In our studies, the benchmark from nonparametric models (RF and k-NN) had a relatively reliable performance (see Table 4). However, they still showed slightly higher RMSE values than both uncalibrated and calibrated NLME models in site-level LOOCV, which confirmed the advantage of the parametric model for extrapolating data outside the coverage of the model fitting data [36].
The inclusion of random-effects parameters led to a further improvement in both the model fitting and prediction. Specifically, the calibrated NLME prediction with a small number of resampling field measurements led to more stability and a better accuracy than other methods across all UAVLS-inventoried areas. Although using more variables can improve the generalization of the model, the site environmental variation and the uncertainty in the data acquirement process reduce the model transferability [39,57]. Many studies have attempted to analyze the transferability of airborne LiDAR-based models. Breidenbach et al. (2008) employed mixed-effects models across separate datasets. They found that the mixed-effects model (with both fixed and random effects) was able to more precisely predict DBH using stand attributes from two different inventory areas than the Remote Sens. 2021, 13, 24 16 of 21 fixed-effects model [84]. Korhonen et al. (2019) applied an established mixed-effects model to other inventory areas outside the scope of the data used in the model's development process [57]. Their results revealed that measuring a small number of calibration trees could decrease systematic errors, which increases the model's transferability. The mixed-effects model includes both fixed and random effects parameters, and so can be adapted to a specific site. It can reflect not only the general trend of the sites but also the variation between individuals. The mixed-effect model's transferability can be utilized to calibrate the models for further purposes, providing a new way to improve the model's portability and ductility, especially for UAV-based small-scale forest applications.
Different from previous studies, which mainly applied plot-level mixed-effects models [57,63,85], in the current study we introduced a site-level NLME model to explain the site variability in the tree growth using both LiDAR point clouds and field data. Even after applying similar scanning and flight parameters, different UAVLS devices may produce remarkably different point clouds due to the various forest conditions [6,41]. For smallarea forest inventories using UAV platforms, high-precision point clouds could obtain the individual tree information at a 1 km 2 coverage per scan, breaking the spatial limitation of traditional forest survey sample plots [14,86]. After establishing the models, site-level random effects can be directly and easily used to calibrate the entire UAV flight site, without needing to establish plot-level mixed-effects models for independent calibration. Although a higher prediction accuracy can be obtained by using a smaller scale of mixed effects, a shift from the sample plot survey to focus on the whole stand is needed when UAV is used for small-area forest surveys. The approach proposed in this study might provide a more affordable option for operational forestry.
Using an optimal number of sample trees in mixed-effects model calibration will deliver a relatively high prediction accuracy and appears to be a more efficient strategy for forest management [37,85]. Both of the sampling strategies (based on subplots and individual trees) behaved logically in relation to the amount of calibration information (Figure 4), as corroborated by previous studies on forest modeling [63,65,66]. In the field measurements, plot delineation often required extra labor compared to random tree sampling. As a result, when comparing the prediction accuracy (Section 3.2.2) we only selected 20 sample trees for local calibration, which is more realistic in the actual application. From the prediction results among different age groups, the calibration performance exhibits the lowest improvements in the young stand because the higher stand density in young forests leads to a lower calibration quantity. Proportional resampling methods may obtain more consistent improvements in forests with different densities, but increase the workload of high-density stands. Extending the sample quantity will improve the model's accuracy in a linear correlation with the increasing inventory cost; thus, the optimum practical calibration approach should be determined based on a compromise between accuracy and efficiency [85].
From the perspective of practical applications, another problem that needs to be addressed is the impact of individual tree segmentation errors. It seems that there were relatively high commission errors in our study for larch plantations compared with others; this is mainly due to the fact that the segmentation methods was sensitive to the outer trees with their crowns extending into the plot boundaries, as shown in Figure 3. However, this indeed has few impacts on the subsequent parameter extraction and DBH modeling in our study, since we only used matched trees and the incomplete segmented crowns along plot boundaries were also excluded for model development. In practical application in forest inventory, commission and omission errors may lead to potential errors in the applications of established DBH estimation models [20]. In particular, segmentation errors and DBH estimation bias may cause error transfer from individual tree to plot-level forest parameter estimation such as stem diameter distribution and forest biomass estimation [18,19,87]. An edge-tree correction could have a significant contribution for the following individual as well as plot-level application of the developed models, as in [88] and [89]. On the other hand, many studies have revealed that stand densities and dominant positions may strongly affect the performance of individual tree delineation. More segmentation errors brought about greater challenges in DBH estimation for small trees and high-density young forests [90,91]. Together with the relatively poor prediction accuracy in young stands (Figure 7), it is necessary to further explore the proposed method for application in young trees. In addition, our study only proposed a species-specific DBH estimation modeling method. In practical applications, it is often necessary to classify tree species in the first step and then choose the modeling approach depending on the tree species, as in [55].
With the development of LiDAR technology, many studies have attempted to directly extract DBH from high-density UAV-LiDAR point clouds. The previously published study of Wieser et al. (2017) focused on utilizing UAVLS point cloud data to estimate DBH from manually delineated tree stems during the leaf-off season in a deciduous forest [92]. Kuželka et al. (2020) applied an automatic diameter measurement and tree stem detection procedure in mixed Norway spruce and Scots pine forests with a stand density of less than 500 stems/ha; they achieved a relative RMSE value of 19% for DBH estimation [93]. Liang et al. (2019) compared the results of manual and automatic DBH measurements using UAVLS. The relative RMSE of automated and manual measurements was around 15-30% with different stand complexities [86]. These studies implied that DBH direct measurements lacked robustness, as the performance was affected by the penetration of laser sensors and forest conditions [86]. Overall, we are looking forward to using direct extracted UAVLS point cloud DBH for calibration in future research, minimizing the burden of field measurement and model application.

Conclusions
This study introduces a framework for model-based individual tree diameter at breast height (DBH) estimation in UAVLS forest inventories. The tree-and plot-level LiDARderived metrics (H, CD, CC p75 , and Cnp R ) were found to be statistically significant, and hence were included in the DBH model development process. Although the selected variables improved the model's generalization ability, there are still several unexplained sources of variation, which can be further described by adding the site-level random effects. The site-calibrated NLME model showed a more stable performance across different sites and achieved a higher prediction accuracy than other approaches, such as uncalibrated NLME, uni-or multivariable weighted nonlinear least square regression, and nonparametric regressions (random forest and k-nearest neighbors). Furthermore, the calibration led to logical behavior with respect to the amount of calibration information. The practitioners could realistically choose the sample sizes and calibration method according to the tradeoff between the accuracy requirements and field-measurement costs. Utilizing site-level random effects improved the transferability of the LiDAR-based DBH estimation model, leading to a breakthrough in how we interact with forests in the future. The mixed-effect modeling approach is a flexible method and provides a foundation for UAVLS-based inventories in small-scale forests.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-4 292/13/1/24/s1. Figure S1: A linear fit between field-measured and LiDAR-derived tree height, Figure S2: Computation of the competition index based on crown cross-sectional areas calculated at a reference height equal to a certain percentage of the height of the subject tree, Figure S3: Out-of-bag (OOB) error with variables being removed by a backward stepwise variable selection of random forests. The dash line represents the number of variables equal to 85, Table S1: Summary of the treeand plot-level metrics derived from unmanned aerial vehicle light detection and ranging (UAVLS) point clouds that were used for the DBH estimation, Table S2: The basic equation forms considered for the base model selection.