Prediction of Individual Tree Diameter and Height to Crown Base Using Nonlinear Simultaneous Regression and Airborne LiDAR Data

Abstract

The forest growth and yield models, which are used as important decision-support tools in forest management, are commonly based on the individual tree characteristics, such as diameter at breast height (DBH), crown ratio, and height to crown base (HCB). Taking direct measurements for DBH and HCB through the ground-based methods is cumbersome and costly. The indirect method of getting such information is possible from remote sensing databases, which can be used to build DBH and HCB prediction models. The DBH and HCB of the same trees are significantly correlated, and so their inherent correlations need to be appropriately accounted for in the DBH and HCB models. However, all the existing DBH and HCB models, including models based on light detection and ranging (LiDAR) have ignored such correlations and thus failed to account for the compatibility of DBH and HCB estimates, in addition to disregarding measurement errors. To address these problems, we developed a compatible simultaneous equation system of DBH and HCB error-in-variable (EIV) models using LiDAR-derived data and ground-measurements for 510 Picea crassifolia Kom trees in northwest China. Four versatile algorithms, such as nonlinear seemingly unrelated regression (NSUR), two-stage least square (2SLS) regression, three-stage least square (3SLS) regression, and full information maximum likelihood (FIML) were evaluated for their estimating efficiencies and precisions for a simultaneous equation system of DBH and HCB EIV models. In addition, two other model structures, namely, nonlinear least squares with HCB estimation not based on the DBH (NLS and NBD) and nonlinear least squares with HCB estimation based on the DBH (NLS and BD) were also developed, and their fitting precisions with a simultaneous equation system compared. The leave-one-out cross-validation method was applied to evaluate all estimating algorithms and their resulting models. We found that only the simultaneous equation system could illustrate the effect of errors associated with the regressors on the response variables (DBH and HCB) and guaranteed the compatibility between the DBH and HCB models at an individual level. In addition, such an established system also effectively accounted for the inherent correlations between DBH with HCB. However, both the NLS and BD model and the NLS and NBD model did not show these properties. The precision of a simultaneous equation system developed using NSUR appeared the best among all the evaluated algorithms. Our equation system does not require the stand-level information as input, but it does require the information of tree height, crown width, and crown projection area, all of which can be readily derived from LiDAR imagery using the delineation algorithms and ground-based DBH measurements. Our results indicate that NSUR is a more reliable and quicker algorithm for developing DBH and HCB models using large scale LiDAR-based datasets. The novelty of this study is that the compatibility problem of the DBH model and the HCB EIV model was properly addressed, and the potential algorithms were compared to choose the most suitable one (NSUR). The presented method and algorithm will be useful for establishing similar compatible equation systems of tree DBH and HCB EIV models for other tree species.

1. Introduction

A tree crown is characterized by crown height, crown width, crown density, leaf area, and crown ratio, and their measurements are useful for forest management and research. The crown ratio is considered a reliable indicator of the vigor and potential growth of a tree [1,2,3,4]. Height to crown base (HCB) is an important tree measure to derive crown ratio and is also regarded as an indicator of log quality. HCB is usually understood as the vertical height from the ground to the bottom of live whorled branch on the bole of a tree [5]. The ground-based measurement of HCB is a time-consuming and labor-intensive process; thus, it is rarely done during field inventory [6,7]. Most researchers have obtained the HCB value by establishing linear or nonlinear HCB models with other variables as predictors, such as DBH, tree height, basal area, basal area larger than a target tree, the sum of basal area of all trees with diameter bigger than a target tree, crown competition factor, climate, and site index [8,9,10,11,12]. Tree diameter at breast height (DBH) is also an important tree attribute that is used as a main predictor in forest growth and yield, taper, and biomass models. In general, the measurement of DBH is very common in ground-based inventory; however, field-inventory data could have a low accuracy, and their measurement needs more time and cost, especially measurements required for extensive forest areas. Therefore, methods of HCB data collection have been transformed from the traditional forest field inventory to modeling and prediction based on remote sensing technology [13,14,15,16].

Light detection and ranging (LiDAR) can accurately determine the geographical position of surface objects by transmitting and receiving laser pulses. Laser pulses travel down the forest canopy, and detailed information on the three-dimensional structures of the forest canopy and understory topography can be obtained [17]. Many tree attributes, such as tree height and crown dimensions [18] can be obtained based on the LiDAR data. The study approaches based on HCB prediction may be divided into two categories: direct and indirect approaches. The direct approaches refer to those derived from HCB with various geometrical shapes of the crown [12,19,20,21,22,23] or predicting HCB according to descriptive statistics of the LiDAR-based data distribution [4,24]. Direct approaches do not require any ground-measured HCB data, which are costly and time-consuming, as they only require point-cloud data processing and analysis including tree detection and the determination of crown base positions. In addition, this approach could also cause considerable uncertainties in determining the base of the first normal green branch as a part of the crown. Therefore, its application is quite limited to estimating HCB. The indirect approach, on the other hand, refers to predicting HCB through the application of statistical modeling [22,25,26,27,28]. This approach requires field-measured HCB data to establish the models for the prediction of HCB. The models for the accurate prediction of individual tree HCB can be built using LiDAR-based information, and so this method has been frequently used in recent years [22,25,26,27,28].

The application of ordinary least square (OLS) regression to estimate the parameters of LiDAR-based DBH and HCB models is not generally preferred, but it is still used [16,29]. This estimation method usually assumes that (i) regressors are random variables with errors, (ii) regressors are fixed variables without errors, and (iii) the associated error is subject to normal distribution with zero mean and constant variance [30]. Any violation of the second assumption leads to the substantially biased estimation of the models [30], which eventually reduces the prediction accuracy.

The prediction accuracy of the developed HCB and DBH models uses the LiDAR-based tree height, crown width and crown area may not be always satisfactory for a couple reasons. Firstly, LiDAR-based tree height, crown width, and crown area have random or systematic errors caused by LiDAR system configuration and parameter estimation. Any error involved in the variables could increase the residual variance of the model and also lead to invalid statistical tests [31,32]. Secondly, the estimated DBH from a LiDAR-based DBH estimation model contains non-ignorable or inevitable errors [33]. If such erroneous DBH is used as a predictor in a LiDAR-based HCB model, substantial bias would occur due to error transfers [34]. In addition, estimating with a LiDAR-based DBH model and a LiDAR-based HCB model separately or independently using OLS disregards the inherent correlations of HCB with DBH and thus fails to account for the compatibility of the estimated HCB and DBH. Thus, estimating the parameters of both model types independently with OLS may create a remarkable problem, especially in the condition when errors are associated with both the regressors and response variables. An appropriate settlement of this problem is to apply error-in-variable (EIV) modeling, which takes the errors into consideration and can guarantee compatibility between HCB and DBH [29,35,36,37].

Fuller [35] first introduced the theory on the development and application of linear EIV models, and, later on, Carroll et al. [32] applied this concept on the nonlinear EIV modeling in detail. Kangas [31] investigated the effects of EIV on the parameters of the diameter growth model and applied the simulation extrapolation algorithm to adjust the errors in the estimated parameters. Lindely [38] proved that validation data from the same population as the fitting data resulted in predictions that were usually unbiased, even though the regressors were subject to error. Tang and Zhang [36] developed an EIV model to investigate the unbiased parameter estimates. Tang and Wang [39] proposed the two-stage EIV method to estimate the model parameters. In their study, the EIV concept was introduced into forest attribute modeling, which provides a theoretical basis for studying the influence of errors on stand growth and harvest models. Li and Tang [40] compared three methods, namely simulation extrapolation, regression calibration, and EIV to estimate the models and found a better performance with EIV with smaller variances compared to other two methods.

Few studies have been carried out with DBH EIV modeling using remote sensing data. For example, Fu et al. [33] developed an individual tree DBH and above-ground biomass (AGB) EIV model with LiDAR-based tree height and crown projection area as predictors with the application of the two-stage error-in-variable modeling (TSEM) and nonlinear seemingly-unrelated regression (NSUR) to estimate model parameters. Both TSEM and NSUR explain the correlations of DBH with AGB and also effectively explain the errors in DBH on the prediction of AGB. Zhang et al. [29] reported that the DBH EIV model developed with errors associated with both response and regressor variables through the application of the maximum likelihood method was most appropriate. To the authors’ knowledge, no studies have been carried out on developing LiDAR-based HCB EIV models that were attributed to compatibility.

This study thus aimed (a) to develop a compatible simultaneous equation system of DBH and HCB EIV models based on the LiDAR data at the individual tree level for Picea crassifolia Kom forests in northwest China, (b) to evaluate the compatibility of two different nonlinear OLS-based DBH and HCB models with the leave-one-out cross validation method, and (c) to compare various unbiased fitting algorithms including NSUR. To simplify the proposed simultaneous equation system and to guarantee its application in the future, only response variables (HCB and DBH) were assumed as the error-in-variables [39], and predictor variables were regarded as error-out variables [33]. The presented compatible simultaneous equation system of DBH and HCB models will be applicable to other Picea species whose growth and stand conditions are very much similar to the basis of our studied species. This tree species is crucial to the economic and social development of the rural population, as well as regional carbon storage and cycling, and the maintenance of the structures and functions of the forest ecosystems in northwest China. This article is mainly concerned with the methodology employed in this study, which is clearly described in the Methods section; additionally, the major strengths and weaknesses of the methodologies, along with the main findings of the study, are thoroughly discussed while the potential contribution of the study is highlighted.

2. Methods

2.1. Data Collection

The study site is located at the Xishui forest farm of the Su’nan Yuguzu autonomous county, Gansu province (38°29′–38°35′N, 100°12′–100°20′E) (Figure 1a) with Picea crassifolia Kom as the dominate tree species. The climate in this field is a temperate semi-arid zone. It is covered by mountainous forests. Slopes with south-facing aspect are covered by grass, and the slopes with north-facing aspect are covered by natural secondary pure forests with one dominating tree species of Picea crassifolia Kom. The ground is covered by a moss floor, and the average elevation here is around 2993 m. The typical soil type is sandy loam. Along the hill, we established a permanent sample plot (PSP) with 100 m long and 100 m wide in 2008, and the PSP was divided into sixteen sub-plots that were 25 m long and 25 m wide. The PSP designed in this study was very representative of the entire forest of the study area and was mainly used for the carbon flux observation and dynamic monitoring of forest quality.

Airborne LiDAR data were acquired by the LiteMapper 5600 system with laser scanner—Riegl LMS-Q560 by a specification of 50 kH pulse repetition frequency, a 49 HZ scanning frequency, and a 30° maximum scanning angle [41]. The LiDAR data were collected on 23 June 2008, and field-measured data were collected on 1 June through 13 June 2008. The wavelength was 1550 nm, and the pulse length and laser beam divergence were 3.5 ns and 0.5 mrad, respectively. The average flight height was 3699 m, and the average flight speed was 230 km h⁻¹. The scanner’s pulse repetition frequency, scanning frequency, and maximum scanning angle were 50 kHZ, 49 Hz, and 30°, respectively, and the mean density of point cloud was 4.34 m⁻². The spatial distribution of neighbor smoothed 510 Picea crassifolia Kom trees is shown in Figure 1b. Data summary is presented in

Table 1. Descriptive statistics of tree measurements (SD, standard deviation).

Variable	Min.	Max.	Mean	SD
LH (m)	4.62	22.15	13.84	3.20
HCB (m)	0.90	10.20	4.80	3.52
LCW (m)	2.00	7.50	4.20	0.95
DBH (cm)	3.60	81.10	22.57	8.54
LCA (m)	3.19	38.63	12.21	5.47

, and the relationships of HCB with DBH, LiDAR-derived tree height (LH), LiDAR-derived crown width (LCW), and LiDAR-derived crown projection area (LCA) are shown in Figure 2.

The correlation analysis of LiDAR-derived tree attributes and ground-measured tree attributes are shown in Figure 3, which indicates that these LiDAR-derived tree attributes are highly correlated with ground-measure tree attributes. Thus, these LiDAR-derived tree attributes in the sample could be used for our modeling study. Figure 4 presents the Y coordinate value versus the predicted HCB value showing the vertical profile of the LiDAR product.

The point cloud was created from LiDAR waveforms by the data provider [41]. By applying the algorithm of TerraScan 005 (Terrasolid, Helsinki, Finland), the ground points were classified, and this was used to create the digital elevation model (DEM) (Figure 5a) with a 0.25 m resolution. With ground and vegetation points, the digital surface model (DSM) with a 0.25 m resolution was created using the Highest hit z algorithm of TerraScan 005. A canopy height model (CHM) (Figure 5b) with a resolution of 0.25 m and a window size of 3 × 3 m was obtained by subtracting the DSM and DEM [42,43,44]. The pits in the CHM were smoothed by neighbor smoothing algorithm [45]. Using the local maximum method with a window size of 2.0 m to detect the crown top from the CHM, the LH values were estimated as the values of detected crown tops with a prediction accuracy of 0.65. Using the region growing algorithm proposed by Liu et al. [46], the LCW of each tree was estimated to be the average of the horizontal ranges of the identified crown from west to east and north to south [45]. After determining canopy boundary, the LCA was obtained. Ground measurements were done for various tree attributes including individual tree DBH, HCB, crown width, and total tree height (H) of 16 sub-plots for a total of 510 Picea crassifolia Kom trees.

2.2. Base Model

2.2.1. LiDAR–DBH Base Model

Fu et al. [33] developed an exponential LiDAR-based DBH model using LH and LCA as predictors for Picea crassifolia Kom and found a significantly higher prediction accuracy than other three candidate LiDAR-based DBH model forms (linear, Richards, and logistic). Our preliminary analyses exhibited the biggest $R^2$ and the smallest root mean square error (RMSE) of the exponential LiDAR-based DBH model form, indicating its greatest suitability according to our data characteristics, and its prediction accuracy could be further improved by including LH and LCW as predictors (Equation (1)):

$$DBH={\beta }_1\exp (-{\beta }_{2}LH-{\beta }_{3}LCW)+{\epsilon }_{DBH}$$

(1)

where ${\beta }_1,{\beta }_2,{\beta }_3$ are parameters to be estimated, and ${\epsilon }_{DBH}$ is a residual error.

2.2.2. LiDAR–HCB Base Model

Similar to Walters and Hann [40], we used the logistic model as a LiDAR–HCB base model, which had DBH and LCA as predictors in this study.

$$HCB=\frac{LH}{1+\exp ({\gamma }_{1}DBH+{\gamma }_{2}LCA)}+{\epsilon }_{HCB}$$

(2)

where ${\gamma }_1$ and ${\gamma }_2$ are parameters to be estimated, and ${\epsilon }_{HCB}$ is a residual error.

2.3. A Compatible Individual Tree DBH and HCB EIV Equation System

A compatible equation system consisting of tree-based DBH and HCB EIV models (Equation (3)) was built by integrating both the LiDAR–DBH base model (Equation (1)) and the LiDAR–HCB base model (Equation (2)) by following the methods suggested by existing modeling studies [36,37,47].

$$\{\begin{array}{l}db{h}_i={\beta }_1\exp (-{\beta }_{2}L{H}_i-{\beta }_{3}LC{W}_i)\\ hc{b}_i=L{H}_i/(1+\exp ({\gamma }_{1}DB{H}_i+{\gamma }_{2}LC{A}_i)\\ DB{H}_i=db{h}_i+{\epsilon }_{DB{H}_i}\\ HC{B}_i=hc{b}_i+{\epsilon }_{DB{H}_i}\\ {\epsilon }_i={\epsilon }_{DB{H}_i}+{\epsilon }_{HC{B}_i}\end{array}$$

(3)

where $DB{H}_i$ (cm) and $HC{B}_i$ (m) $(i=1,2\dots ,N)$ are the ground-measured diameter at breast height with errors and height to crown base with errors of the i^th tree, respectively; $db{h}_i$ and $hc{b}_i$ are true values (with the assumption of no errors) of $DB{H}_i$ and $HC{B}_i$, respectively; ${\epsilon }_{DB{H}_i}$ and ${\epsilon }_{{HCB}_i}$ represent the errors of $DB{H}_i$ and $HC{B}_i$, respectively. Error ${\epsilon }_i$ is a two-dimensional vector that is assumed to be normally distributed with zero means and variance–covariance matrix $\Sigma$; $L{H}_i$, $LC{W}_i$, and $LC{A}_i$ are the LiDAR-derived tree height (m), crown width (m), and crown projection area (m²) of the i^th tree, respectively. In this simultaneous equation system (Equation (3)), both DBH and HCB are the EIV, while LH, LCW, and LCA are regarded as error-free variables. The other parameters and variables are the same as defined above. The elements in the variance–covariance matrix $\Sigma$ were applied to account for the inherent correlations of DBH with HCB.

It was assumed that the simultaneous equation system (Equation (3)) with an error term ${\epsilon }_{ti}$ (t = DBH and HCB; and i = 1,…, N) that were not correlated among the observations but were contemporaneously correlated across the sub-models. For each observation, we assumed that:

$$\Sigma =\left(\begin{array}{cc}{\sigma }_{DBH\times DBH}& {\sigma }_{DBH\times HCB}\\ {\sigma }_{HCB\times DBH}& {\sigma }_{HCB\times HCB}\end{array}\right)$$

(4)

where ${\sigma }_{DBH\times DBH}$, ${\sigma }_{HCB\times HCB}$, ${\sigma }_{DBH\times HCB}$, and ${\sigma }_{HCB\times DBH}$ are the variance and covariance related elements for both DBH and HCB.

The covariance matrix of the stacked error terms ($\epsilon =({\epsilon }_{DBH}^T,{\epsilon }_{HCB}^T{)}^T({\epsilon }_t={({\epsilon }_{t1},\dots ,{\epsilon }_{tN})}^T$,$t=DBHandHCB$) would be $R=\Sigma \otimes I_N$.

2.4. Parameter Estimation

Four commonly used algorithms, such as NSUR, two-stage least square (2SLS), three-stage least square (3SLS), and full information maximum likelihood (FIML) were applied to estimate the parameters $B=({\mathsf{\beta }}_{\mathrm{DBH}}^T=({\beta }_1,{\beta }_2,{\beta }_3),{\gamma }_{HCB}^T=({\gamma }_1,{\gamma }_2\left)\right)$ in our simultaneous equation system (Equation (3)). We briefly describe these algorithms in a methodological flow chart (Figure 6), and the details are given in the sub-sections below.

(1) NSUR algorithm

The NSUR algorithm considers the disturbance across the two equations as a linkage of the equation system but assumes that disturbances are uncorrelated across the observations; thus, this algorithm is known as a seemingly unrelated regression. The estimation of parameters in the simultaneous equation system (Equation (3)) was done using the NSUR algorithm [34,48,49,50] with the feasible generalized least square regression method described as follows:

Step 1: Two sub-models (Equations (1) and (2)) in the simultaneous equation system (Equation (3)) were fitted with the NSUR algorithm, and the resulting residuals ${\widehat{\epsilon }}_t$ $\left(t=DBHandHCB\right)$ were used for estimating the variance–covariance matrix, $\Sigma$. The residuals of each sub-model were estimated with OLS using following formula:

$${\widehat{\sigma }}_{ij}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\epsilon }_{it}{\epsilon }_{jt}}$$

(5)

The estimated variance–covariance matrix $\Sigma$ $\mathrm{is}$ given by:

$$\widehat{\Sigma }=\left(\begin{array}{l}{\widehat{\sigma }}_{11}{\widehat{\sigma }}_{12}\dots {\widehat{\sigma }}_{1n}\\ {\widehat{\sigma }}_{21}{\widehat{\sigma }}_{22}\dots {\widehat{\sigma }}_{2n}\end{array}\right)$$

(6)

Step 2: Based on the estimated $\widehat{\Sigma }$, a covariance matrix R was defined as $\widehat{R}=\widehat{\Sigma }\otimes I_N$. The parameters in the simultaneous equation system (Equation (3)) were estimated using a feasible generalized least square method.

$$\widehat{B}={(X^T({\widehat{\Sigma }}^{-1}\otimes I_N)X)}^{-1}({\widehat{\Sigma }}^{-1}\otimes I_N)y$$

(7)

(2) 2SLS algorithm

The 2SLS algorithm was applied using the following steps [51]:

Step 1: By composing the reduced function (Equation (8)) for all the error-in variables on the right side of an equation, the estimated error-in variables were obtained using OLS:

$$Y_i={\mathsf{\Pi }}_i{X}_i+{\epsilon }_i$$

(8)

where $X_i$ is the vector of the error-free variables and ${\mathsf{\Pi }}_i$ is the i^th parameter vector for X.

Step 2: The error-in variables on the right side were replaced with estimated error-in variables ${\widehat{Y}}_i$, and the parameters were estimated using OLS.

$$\mathrm{cov}(\widehat{B})={(X^T(diag({\Sigma }^{-1})\otimes I_N)X)}^{-1}$$

(9)

(3) 3SLS algorithm

The 3SLS algorithm considers the correlation of disturbance terms among different equations. It was carried out with following steps [51]:

Step 1: The same as step 1 in 2SLS.

Step 2: The same as step 2 in 2SLS; in addition, the disturbance ${\epsilon }_i$ was estimated.

Step 3: The covariance ${\sigma }_{ij}$ and ${\widehat{\epsilon }}_i$ was estimated in step 2, and then an estimated variance–covariance matrix, $\widehat{\Sigma }$, was obtained. $\mathsf{\Phi }=\Sigma \otimes I_N$ was defined, and parameters were estimated with a feasible generalized least squares regression.

$$\widehat{B}={(X^T({\widehat{\Sigma }}^{-1}\otimes I_N)X)}^{-1}({\widehat{\Sigma }}^{-1}\otimes I_N)y$$

(10)

$$\mathrm{cov}(\widehat{B})={(X^T(diag({\widehat{\Sigma }}^{-1})\otimes V_{\mathrm{N}})X)}^{-1}$$

(11)

where V is a matrix of the instrumental variables and I is an identity matrix.

(4) FIML algorithm

Instead of only making use of the reduced function information, FIML makes full use of all the information by estimating all the parameters in the simultaneous equation system at the same time. There must be equal number of error-in-variables and sub-models in this equation system, which we used. Otherwise, if the number of endogenous variables is more than that of the sub-models, the limited information maximum likelihood method needs to be applied. The FIML maximizes the following conditional log-likelihood function [52]:

$$Q_n(B,\Sigma )=-\frac{M}{2}\ln (2\pi )-\frac{1}{2}\ln (|\Sigma |)-\frac{1}{2}{\displaystyle \sum _{i=1}^m{(Y_i-y_i)}^T{\Sigma }^{-1}(Y_i-y_i)}$$

(12)

where $M$ is the number of sub-models; the other parameters and variables are the same as defined above.

2.5. Other Model Structures for Comparison

2.5.1. Nonlinear Least Squares with HCB Estimation not Based on DBH (NLS and NBD)

The DBH in the LiDAR–HCB base model (Equation (2)) was substituted by the LiDAR–DBH base model (Equation (1)). Therefore, in this case, the HCB estimation was independent of DBH. The HCB based on DBH model is given by:

$$\begin{array}{l}HC{B}_i=L{H}_i/(1+\exp ({\gamma }_{1}DB{H}_i+{\gamma }_{2}LC{A}_i))+{\epsilon }_{HC{B}_i}\\ =L{H}_i/(1+\exp ({\mu }_1\exp (-{\mu }_{2}L{H}_i-{\mu }_{3}LC{W}_i)+{\mu }_{4}LC{A}_i)+{\tilde{\epsilon }}_{HC{B}_i}\end{array}$$

(13)

where ${\mu }_1={\beta }_1{\gamma }_1$, ${\mu }_2={\beta }_2$, ${\mu }_3={\beta }_3$, and ${\mu }_4={\gamma }_2$ are parameters in the model, and the error of the $HC{B}_i$ was changed into:

$${\tilde{\epsilon }}_{HC{B}_i}=f(L{H}_i,LC{W}_i){\epsilon }_{DB{H}_i}+{\epsilon }_{{HCB}_i}$$

(14)

$$f(L{H}_i,LC{W}_i)={\beta }_1{\gamma }_1\exp (-{\beta }_{2}L{H}_i-{\beta }_{3}LC{W}_i)$$

(15)

With this method, $DB{H}_i$ was estimated by the LiDAR–DBH base model (Equation (1)), and $HC{B}_i$ was estimated by the NLS and NBD model (Equation (13)). It should be noted that the inherent correlations of HCB with DBH could not be addressed for this method. In addition, the compatibility between the estimated DBH and HCB could not be achieved.

2.5.2. Nonlinear Least Squares with HCB Estimation Based on DBH (NLS and BD)

The LiDAR–DBH base model (Equation (1)) and the LiDAR–HCB base model (Equation (2)) were fitted separately based on the database by the NLS and BD. This method was applied to quantify the consequences in the HCB estimation by using predicted DBH to take place of an actual value while ignoring its error. The NLS and BD approach could explain the compatibility between DBH and HCB, but it failed to account for the effect of the errors in the estimated DBH on HCB estimation. The estimated values of $DB{H}_i$ and $HC{B}_i$ $(i=1,2\dots ,N)$ are, respectively, given by:

$$D\widehat{B}{H}_i={\widehat{\beta }}_1\exp (-{\widehat{\beta }}_{2}L{H}_i-{\widehat{\beta }}_{3}LC{W}_i)$$

(16)

$$H\widehat{C}{B}_i=L{H}_i/(1+\exp ({\widehat{\gamma }}_{1}DB{H}_i+{\widehat{\gamma }}_{2}LC{A}_i))$$

(17)

where ${\beta }_1,{\beta }_2,{\beta }_3$, ${\gamma }_1$, and ${\gamma }_2$ are the model parameters; the other parameters and variables are the same as defined above.

2.6. Comparison and Evaluation of Models

We only had 510 observations, which was not enough to divide a full data set into fitting and validation sets. As such, we applied the leave-one-out cross validation (LOOCV) method [53,54] for the validation of the models. Each time, one tree from the full dataset was deleted, and the fitting data set was formed by the remaining trees. A fitting data set was used to fit the DBH and HCB models and estimated their parameters. Using the estimated parameter values, the deleted tree’s DBH and HCB were predicted, and commonly used prediction statistics, such as mean bias($\overline{e}$), variance of bias (${\sigma }_e^2$), RMSE, and mean absolute error (MAE) (Equations (18)–(21)) were computed with the difference obtained from the predicted and observed values. Then, we put the tree back in place, deleted another tree, and performed the same model-fitting and prediction processes. This procedure was performed on all the trees in the full data set. We present the LOOCV computational codes with NSUR algorithm as an example in Appendix A, and we used these codes to evaluate the equation system.

Finally, the prediction performance of the simultaneous equation system (Equation (3)) was estimated with each of the six different methods: NSUR [33,47,55], 2SLS, 3SLS, FIML, NLS and BD [33], and NLS and NBD [33] were evaluated by three statistics including mean bias, bias variance, and root mean square error that were calculated with Equations (18)–(21). The model with the smallest $\overline{e}$, ${\sigma }_e^2$, RMSE, and MAE were defined as the final model to predict DBH and HCB. We performed all computations with R software version 3.4.4 [56].

$$\overline{e}={\displaystyle \sum _{i=1}^N{e}_i/N={\displaystyle \sum _{i=1}^N({\mathrm{y}}_i-{\widehat{y}}_i}})/N$$

(18)

$${\sigma }_e^2={\displaystyle \sum _{i=1}^N{(e_i-\overline{e})}^2/(N-1)}$$

(19)

$$RMSE=\sqrt{{\overline{e}}^2+{\sigma }_e^2}$$

(20)

$$\mathrm{MAE}=|{\displaystyle \sum _{i=1}^N{e}_i|/N=|{\displaystyle \sum _{i=1}^N({\mathrm{y}}_i-{\widehat{y}}_i}})|/N$$

(21)

where ${\mathrm{y}}_i$ and ${\widehat{\mathrm{y}}}_i$ are the measured and estimated height to crown base or DBH for the $i^{th}$ observation, N is the number of observations, $\overline{e}$ is the mean bias, ${\sigma }_e^2$ is the variance of bias, $RMSE$ is the root mean square error, and MAE is the mean absolute error.

3. Results

For the DBH model, the RMSE of NSUR was identical to the NLS and BD model and smaller than that of 2SLS, 3SLS, and FIML. For the HCB model, the RMSE of NSUR was smaller than that of the NLS and BD model, 2SLS, 3SLS, and FIML. The MAE of NSUR for the HCB model was the smallest.

3.1. Parameters Estimation

All the parameters in the LiDAR–DBH base model (Equation (1)), the LiDAR–HCB base model (Equation (2)), and the simultaneous equation system (Equation (3)) were estimated with four different methods, namely NSUR, 2SLS, 3SLS, and FIML using all the data. Most of the parameter estimates were significantly different from zero, and their magnitudes and signs could meet biological logics, except for parameters ${\mu }_1,{\mu }_2, \mathrm{and} {\mu }_3$ for NLS and NBD and ${\gamma }_1$ for both 2SLS and 3SLS, which were not significant (p < 0.05) (

Table 2. Parameter estimates of the LiDAR–DBH base model (Equation (1)), the LiDAR–HCB base model (Equation (2)), and the NLS and BD model (Equation (13)), as well as the simultaneous equation system (Equation (3)). The first three models were estimated using ordinary least squares regression, and the last one was estimated using NSUR, 2SLS, 3SLS, and FIML.

Model	Method	Parameters	Estimates	Standard Error	t-Value
LiDAR–DBH base model (Equation (1))	NLS	${\beta }_1$	5.7161	0.3599	15.882
		${\beta }_2$	−0.0567	0.0061	−9.344
		${\beta }_3$	−0.1264	0.0178	−7.108
		${\sigma }^2$	37.22
LiDAR–HCB base model (Equation (2))	NLS	${\gamma }_1$	0.0053	0.0031	1.68
		${\gamma }_2$	0.0367	0.0057	6.4340
		${\sigma }^2$	3.37
HCB based on DBH model (Equation (13))	NLS	${\mu }_1$	−0.0003	0.0547	−0.00
NLS and BD		${\mu }_2$	0.2472	0.1987	−1.24
		${\mu }_3$	0.0800	0.0520	1.54
		${\mu }_4$	0.0283	0.0137	2.07
		${\sigma }^2$	3.33
Simultaneous equation system (Equation (3))	NSUR	${\beta }_1$	5.6992	0.3615	15.77
		${\beta }_2$	−0.0563	0.0061	−9.27
		${\beta }_3$	−0.1283	0.0180	−7.12
		${\gamma }_1$	0.0059	0.0032	1.87
		${\gamma }_2$	0.0361	0.0058	6.26
		${\sigma }_{DBH}^2$	37.22
		${\sigma }_{HCB}^2$	3.32
		${\sigma }_{DH}$	0.608
Simultaneous equation system (Equation (3))	2SLS	${\beta }_1$	5.2496	0.3586	14.64
		${\beta }_2$	−0.0567	0.0064	−8.83
		${\beta }_3$	−0.1451	0.0186	−7.80
		${\gamma }_1$	0.0054	0.0051	1.04
		${\gamma }_2$	0.0370	0.0091	4.08
		${\sigma }_{DBH}^2$	37.93
		${\sigma }_{HCB}^2$	3.32
		${\sigma }_{DH}$	0.551
Simultaneous equation system (Equation (3))	3SLS	${\beta }_1$	5.2641	0.3590	14.66
		${\beta }_2$	−0.0564	0.0064	−8.78
		${\beta }_3$	−0.1457	0.0186	−7.83
		${\gamma }_1$	0.0052	0.0051	1.01
		${\gamma }_2$	0.0373	0.0091	4.12
		${\sigma }_{DBH}^2$	37.88
		${\sigma }_{HCB}^2$	3.32
		${\sigma }_{DH}$	0.553
Simultaneous equation system (Equation (3))	FIML	${\beta }_1$	5.7036	0.3469	16.44
		${\beta }_2$	−0.0561	0.0047	−11.84
		${\beta }_3$	−0.1289	0.0149	−8.66
		${\gamma }_1$	0.0074	0.0042	1.76
		${\gamma }_2$	0.0335	0.0076	4.39
		${\sigma }_{DBH}^2$	37.72
		${\sigma }_{HCB}^2$	3.31
		${\sigma }_{DH}$	0.65

).

The elements of the variance–covariance matrix were significantly different from each other (p < 0.05) in the simultaneous equation system (Equation (3)), which was estimated using NSUR, 2SLS, 3SLS, and FIML, implying that correlation of DBH with HCB was highly significant.

3.2. Model Prediction

The LOOCV was carried out for the LiDAR–DBH base model (Equation (1)), the LiDAR–HCB base model (Equation (2)), the NLS and NBD model (Equation (13)), and the NLS and BD model (Equation (16)) estimated using the nonlinear OLS, as well as the simultaneous equation system (Equation (3)) estimated using NSUR, 2SLS, 3SLS, FIML, and TSEM. The evaluations and comparisons of all these models were carried out using $\overline{e}$, ${\sigma }_e^2$, and $RMSE$ (

Table 3. Prediction statistics of the models: the DBH-based model (Equation (1)), the HCB-based model (Equation (2)) fitted with NLS, the HCB based on DBH model (Equation (13)) fitted with NLS, and the simultaneous equation system (Equation (3)) fitted with the NSUR, 2SLS regression, 3SLS regression, and FIML algorithms.($\overline{\mathrm{e}}$, mean bias;${\mathsf{\sigma }}_{\mathrm{e}}^2$, bias variance; and $RMSE$, root mean square error. All other acronyms are the same as defined in Table 2).

Fitting Method	Variables	$\overline{\mathit{e}}$	${\mathit{\sigma }}_{\mathit{e}}^2$	RMSE	MAE
NLS and NBD	DBH	−0.0426	37.6572	6.1367	3.6998
	HCB	−0.0367	3.3387	1.8276	1.4619
NLS and BD	HCB	−0.4186	3.3530	1.8784	1.4619
NSUR	DBH	−0.0458	37.6577	6.1368	3.7008
	HCB	−0.0276	3.3414	1.8281	1.4606
2SLS	DBH	0.0093	37.7985	6.1481	3.6984
	HCB	−0.0336	3.3421	1.8285	1.4612
3SLS	DBH	0.0026	37.7977	6.1480	3.7004
	HCB	−0.0341	3.3416	1.8283	1.4613
FIML	DBH	−0.0457	37.6598	6.1369	3.7012
	HCB	−0.0195	3.3470	1.8296	1.4613

).

A compatible DBH and HCB EIV equation system fitted with NSUR showed a better prediction ability than those fitted with other alternative methods (Table 3). For the DBH model, the ${\sigma }_{\mathrm{e}}^2$ of NSUR was identical to that of NSL and BD, as well as 0.37%, 0.37%, and 0.18% smaller than that of 2SLS, 3SLS, and FIML, respectively. The RMSE of NSUR was identical to that of NLS and BD, and it was 0.02%, 0.01%, and 0.08% smaller than that of 2SLS, 3SLS, and FIML, respectively. For the HCB model, the ${\sigma }_{\mathrm{e}}^2$ of NSUR was 0.35%, 0.021%, 0.006%, and 0.17% smaller than that of NLS and BD, 2SLS, 3SLS, and FIML, respectively. The RMSE of NSUR was 2.75%, 0.022%, 0.011%, and 0.082% smaller than that of NLS and BD, 2SLS, 3SLS, and FIML, respectively. The MAE of NSUR for HCB was the smallest.

The residuals of six different alternative models and equation systems were calculated based on a full dataset. This analysis indicated that the mean residuals of the NLS method for HCB were higher than other alternative methods, among which NSUR showed the smallest mean residual for HCB (

Table 4. Descriptive statistics of residuals of the LiDAR– DBH base model (Equation (1)), LiDAR–HCB base model (Equation (2)), and model (Equation (13)), and simultaneous equation system (Equation (3)). The first three models were estimated using ordinary least squares regression and last one was estimated using NSUR, 2SLS, 3SLS, and FIML. (SD, standard deviation).

Model	Method	Response Variable	Min. of Residuals	Max. of Residuals	Mean of Residuals	SD of Residuals
LiDAR–DBH base model (Equation (1))	NLS	DBH	−60.0992	21.6027	0.0423	6.1365
LiDAR–HCB base model (Equation (2))	NLS	HCB	−6.1705	4.3487	0.0367	1.8272
HCB based on DBH model, NLS and BD (Equation (13))	NLS	HCB	−6.1708	4.3495	0.0342	1.8300
Simultaneous equation system (Equation (3))	NSUR	DBH	−60.0883	21.5831	0.0458	6.1366
		HCB	−6.2068	4.4099	0.0276	1.8279
Simultaneous equation system (Equation (3))	2SLS	DBH	−60.1132	22.2219	−0.0093	6.1481
		HCB	−6.2193	4.4061	0.0336	1.8282
Simultaneous equation system (Equation (3))	3SLS	DBH	−60.0874	22.1638	−0.0026	6.1480
		HCB	−6.2150	4.3999	0.0341	1.8280
Simultaneous equation system (Equation (3))	FIML	DBH	−60.0819	21.5710	0.0457	6.1368
		HCB	−6.2407	4.4668	0.0295	1.8295

).

The prediction accuracy of the simultaneous equation system (Equation (3)) fitted with all four fitting algorithms appeared almost identical (Figure 7), indicating that each of the fitting algorithms were able to produce almost equally unbiased estimations and prediction accuracies. The prediction accuracy of DBH seemed to be much higher than that of HCB.

The inherent correlations between the ground-measured DBH, model-estimated DBH, ground-measured HCB, and model estimated-HCB were all significantly high (Figure 8). The inherent correlation between DBH and HCB was substantially high. This figure suggested that all models and all fitting algorithms were appropriately suited to our data.

4. Discussion

HCB is an important tree attribute to assess tree productivity and tree vigor. DBH is commonly used to predict HCB model, but DBH estimated with LiDAR-based attributes contains unignorable errors. In addition, the compatibility between DBH and HCB needs to be considered when estimating HCB. In this study, we investigated four algorithms to estimate DBH and HCB in an EIV equation system—NSUR, 2SLS, 3SLS, and FIML—that were compared with two model structures. The prediction accuracy of the four EIV equation system algorithms and two model structures were reflected by RMSE and MAE. The results showed that the impacts of measurement error of DBH on HCB and the compatibility between DBH and HCB were well accounted for by the NSUR algorithm.

HCB is an important indicator for tree vigor and tree stem form, as well as an indispensable measure for retrieving the crown ratio. However, measuring in-situ HCB is quite labor-intensive and costly, especially when conducted for large forest areas. In this situation, an efficient method of obtaining precise HCB is necessary, which can be possible with the HCB prediction model developed from the LiDAR-derived variables, such as tree height, crown projection area, crown width, and ground-measured DBH. The first three variables can be relatively more accurately and easily measured by applying the advanced remote sensing techniques. The HCB can be estimated from the established HCB model, which may also contain DBH as a predictor [11,28]. The DBH estimation model can also be developed using the LiDAR-derived information [33]. The estimation of HCB and DBH from their corresponding prediction models would be substantially biased if separately developed models were used, i.e., DBH model and HCB models developed independently from each other from the same tree data. In order to overcome such a bias, developing a compatible simultaneous equation system is the most appropriate solution. However, this equation system of DBH and HCB models is still unavailable in forest modeling literature. As mentioned in the introduction, other compatible simultaneous equation systems developed through the EIV modeling approach are available, e.g., a system of equations of DBH and individual tree above-ground biomass models [33]. Considering the knowledge gap, we developed the simultaneous equation system of DBH and HCB models using the tree-level predictors (LH, LCW, and LCA), the information of which was derived from the LiDAR imagery. Four different algorithms (NSUR, 2SLS, 3SLS, and FIML) were used to estimate this equation system.

The data used in our study originated from the Picea crassifolia Kom forest, which is crucial to the economic and social benefits to the rural population, as well as regional carbon storage, regional carbon cycling, and the maintenance of the balanced-functions of forest ecosystems in northwest China. Two different model structures (the NLS and NBD model and the NLS and BD model) built by assuming errors associated with all the regressors and response variables were found to be inappropriate because this approach did not account for the inherent correlations of DBH with HCB and all the estimated parameters and variances were biased.

Generally, the structural estimators or fitting algorithms (NSUR, 2SLS, 3SLS, and FIML) should always be preferred to the NLS, as each of them effectively accounted for the errors in variables in an appropriate way. However, surprisingly, we found that NLS could sometimes provide a closer estimation of the structural estimators applied in this study, and it was the same for NLS and NBD. The NLS and NBD model had a smaller bias variance, so it has possibility to produce a smaller RMSE. However, NLS standard errors are, in all the likelihoods, not useful for inference purposes [57]. The prediction accuracy of the NLS and BD model was the worst with the highest ${\sigma }_{\mathrm{e}}^2$ and the biggest RMSE, thus, in this case, the EIV modeling approach clearly displayed the advantage over NLS. In general, individual tree DBH and HCB models based on the LiDAR data and field-measurements contain errors that exist in image capture, image processing, and the extraction of the information processes, and they are therefore very hard to completely avoid [29,33,58].

The NLS and NBD could neither address the compatibility problem of DBH and HCB nor account for their inherent correlations. However, a simultaneous equation system (Equation (3)) can effectively address these issues. Among the four algorithms used in fitting simultaneous equation system (Equation (3)), NSUR and 2SLS are classified into the limited information estimators, while 3SLS and FIML are the full information estimators. The former two estimators can make use of the reduced model information, while the latter two estimators can make use of full information from the model [33,34]. Based on the model validation results with LOOCV, the prediction accuracy of NSUR was slightly better than that of the other algorithms (2SLS, 3SLS, and FIML). This was probably because NSUR has a better ability to address the error transfers caused by DBH in the simultaneous equation system of the DBH and HCB models. Potentially because of this, Parresol [49] applied NSUR to develop the additive tree biomass models in a pioneer modeling study about a simultaneous equation system in forestry. The prediction accuracy of 3SLS was slightly better than 2SLS, confirming the findings of Tang et al. [34], who found that when errors in across equations were correlated, 3SLS outperformed 2SLS, and—when errors involved across equations were uncorrelated—2SLS outperformed 3SLS.

Our HCB equation system developed in this study was based on the most attractive fit statistics of the base model among the five frequently used HCB base candidate models [10,59,60]. The analysis of correlations between the regressors and HCB showed strong connections among LCA, LH, DBH, and HCB. In other words, these tree characteristics strongly influenced HCB variations. Our DBH base model, which replaced LCPA with LCW in the models of Fu et al. [33], showed a better fitting performance with a smaller RMSE. Both the HCB model applied with all the LiDAR-based data (except for DBH data, which were obtained from ground measurement) and the DBH model were developed by LiDAR data, and this enabled the DBH–HCB-compatible EIV models, suggesting the high possibility of the equation system’s application to an extensive forest area. The validation results based on the LOOCV for NSUR, 2SLS, 3SLS, and FIML were almost identical, even though NSUR slightly outperformed others; however, the prediction difference was still insignificant (Table 3). In this study, we only considered DBH and HCB as error-in-variables; however, other regressors may contain various errors including measurement errors, tree crown delineation errors, and errors of parameter estimation. Ignoring all these errors can cause the complex uncertainties while developing models. Future researchers should focus on these issues. Therefore, readers need to be cautious when considering the conclusion of this study.

As mentioned in the introduction section, this study was based on a novel methodology, which resulted in a system of compatible simultaneous equations of DBH and HCB models in which various LiDAR-derived tree attributes were used. The measurement errors of both DBH and HCB were simultaneously taken into consideration to address the problem of compatibility between DBH and HCB models and to account for inherent correlations between these tree variables through a simultaneous modeling approach. The presented equation system of DBH and HCB models can fulfill the gaps of the unavailability of such an HCB EIV model system in forest modeling literature. A compatible simultaneous equation system of the DBH and HCB models developed using the information of the tree-level predictors (LH, LCW, and LCA) derived from LiDAR imagery and ground-based measurements confirmed the accurate prediction of HCB and DBH. Compared to any of the previously developed HCB models using only ground measurements [11] and those based on LiDAR-derived databases [22,25,26,27,28], the presented equation system in this article will be interesting and useful to both researchers and forest managers, as this system is able to accurately predict HCB. Furthermore, the presented modeling approach and algorithm in this article will be useful for establishing similar compatible equation systems of DBH and HCB EIV models for other tree species and other tree variables that have inherent correlations between themselves.

5. Conclusions

This study developed a compatible simultaneous equation system of DBH and HCB EIV models on the basis of LiDAR-derived and ground-measured data of Picea crassifolia Kom trees in northwest China. Four different algorithms—NSUR, 2SLS regression, 3SLS regression, and FIML—were used to estimate the parameters in an equation system. The NLS used for estimating both the LiDAR–DBH base model (Equation (1)) and the LiDAR–HCB base model (Equation (2)) produced biased results, while the other fitting algorithms used for estimating a simultaneous equation system (Equation (3)) produced unbiased results with similar SSE, MSE, and RMSE. Two additional model structures—nonlinear least squares with HCB estimation not based on DBH (NLS and NBD) and nonlinear least squares with HCB estimation based on the DBH (NLS and BD)—were also developed for comparison. All the fitting algorithms and their resulting models were assessed by a leave-one-out cross validation method. This study indicates that only EIV modeling method can effectively account for the effects of errors associated with the regressors on the response variables and can guarantee the compatibility between DBH model and HCB model at the individual tree level. However, neither the NLS and BD model nor the NLS and NBD model exhibited these advantages. Among the various evaluated algorithms and models, NSUR showed a slightly better performance than the others. The results showed that the methodology proposed in this article is a reliable and efficient, and it can estimate individual tree DBH and HCB from LiDAR-based data over the extensive forest area. In addition, the presented simultaneous equation system (Equation (3)) does not need measurement of any stand-level variable, which would require an additional cost. The presented modeling approach and algorithm will be useful for establishing similar compatible equation systems of DBH and HCB EIV models for other tree species and other tree variables that have inherent correlations between themselves.