Comparing Empirical and Semi-Empirical Approaches to Forest Biomass Modelling in Different Biomes Using Airborne Laser Scanner Data

Airborne laser scanner (ALS) data are used operationally to support field inventories and enhance the accuracy of forest biomass estimates. Modelling the relationship between ALS and field data is a fundamental step of such applications and the quality of the model is essential for the final accuracy of the estimates. Different modelling approaches and variable transformations have been advocated in the existing literature, but comparisons are few or non-existent. In the present study, two main approaches to modelling were compared: the empirical and semi-empirical approaches. Evaluation of model performance was conducted using a conventional evaluation criterion, i.e., the mean square deviation (MSD). In addition, a novel evaluation criterion, the model error (ME), was proposed. The ME was constructed by combining a MSD expression and a model-based variance estimate. For the empirical approach, multiple regression models were developed with two alternative transformation strategies: square root transformation of the response, and natural logarithmic transformation of both response and predictors. For the semi-empirical approach, a nonlinear regression of a power model form was chosen. Two alternative predictor variables, mean canopy height and top canopy height, were used separately. Results showed that the semi-empirical approach resulted in the smallest MSD in three of five study sites. The empirical approach resulted in smaller ME in the temperate and boreal biomes, while the semi-empirical approach resulted in smaller ME in the tropical biomes.


Introduction
Airborne laser scanning (ALS) has become an important source of auxiliary data for estimating tree heights, diameter distribution, timber volume, and forest biomass [1,2].For forest inventory purposes, the ALS data are often applied within an estimation framework known as the area-based approach.In this method, introduced in Naesset [3] and further described in Naesset [4], models are constructed using spatially coincident observations of ground-based response values (e.g., biomass) measured on field plots and variables derived from ALS data.The models are subsequently used to predict forest attributes for individual cells which tessellate the area of interest (AOI).The cell predictions are finally aggregated to estimates for the entire AOI or regions within the AOI.
Hollaus et al. [5] and Hollaus et al. [6] described three main approaches to modelling the relationship between forest attributes and the remotely sensed data: physical, empirical, and semi-empirical.While physical modelling approaches have had limited practical applications due to difficulties with inverting the models [5], both the empirical and semi-empirical approaches are commonly used.Following the empirical approach, a variety of features are extracted from the ALS echoes and related to field-based observations using statistical methods such as regression analysis, nearest neighbors, neural networks, or ensemble learning to model the relationship (e.g., [7][8][9][10]).Universal for these methods when used in the empirical approach is that the model selection is performed objectively.One of the most common methods, at least the most frequently implemented method in commercial operations in Norway, is ordinary least square regression (OLS).To improve the linear relationship between the response and predictors, transformations of the response or of both the response and predictors are often used.Two common transformation approaches are the square-root transformation (SQRT) of the response, and the logarithmic transformation (LOG) of both response and predictors.Comparisons of SQRT and LOG models generally show small differences in terms of model performance.Naesset [11], for instance, reported only small differences in model performance when modelling aboveground biomass (AGB) in young forests in Norway using SQRT, LOG, and linear model forms.Saarela et al. [12] found that SQRT models performed slightly better compared to linear and LOG models when using ALS data to model timber volumes in Finland.The use of SQRT and LOG transformations requires that model predictions are back-transformed to original scale.Although different methods of back-transformation of LOG models have been suggested (e.g., [13][14][15]), achieving unbiasedness is challenging [16].For back-transformation of SQRT models, the commonly applied correction factor presented by Gregoire et al. [17] results in "negligibly small" biased predictions.Regression methods that render the need for transformation of the response obsolete have also been used e.g., OLS without transformation of the response [11], quadratic polynomials [18], nonlinear least square (NLS) [19], and generalized linear modelling (GLM) [20].These methods have however not reached the same level of use in operational settings as the LOG and SQRT.This could be because of reported advantages with the LOG and SQRT methods, or with difficulties in applying the alternative methods.
Although less used in operational settings, the semi-empirical approach is common in research of the relationship between forest attributes and ALS-derived predictors (e.g., [5,[21][22][23]).Following this approach, the model form and predictors are selected a priori, often based on the theoretical relationship between the response and predictor or predictors.Allometry of tree height and girth is often described by a power model [24,25] and is therefore commonly applied to model AGB, substituting tree height with ALS-derived variables and girth with AGB.This approach was taken by Asner,et al. [23] to model forest AGB on the island of Hawaii.The authors regressed AGB against mean canopy height derived from the ALS data using a NLS modelling procedure.Asner et al. [26] advocated the same approach following the argument of a theoretical relationship between mean canopy height and carbon stored in the AGB.Having only one predictor variable, this approach makes it easy to fit the NLS model using statistical software such as e.g., R [27].Increasing the number of predictors makes it increasingly difficult to fit the model.Although Magnussen et al. [21] used two predictors, other authors have used OLS and LOG transformations to achieve the desired model form (e.g., [22,28]).The variables were however selected subjectively making the approach semi-empirical.Bouvier et al. [22] selected four predictor variables based on their theoretical complimentary properties describing the canopy height, heterogeneity of canopy height, canopy cover, and variation in leaf area density.Tompalski et al. [28] took a similar approach, however selecting the most correlated ALS-derived variable to each of the four model properties.Chen [29] modelled AGB in three different study sites in northwestern USA.He compared two power models, each with one predictor variable, one multiple OLS model with a logarithmic transformed response and three pre-selected predictors, and two nonparametric methods.Results showed that the power model with the mean of all ALS echoes as predictor had the smallest estimated prediction error in two of the study sites, whilst the multiple OLS model form with three variables had the smallest estimated prediction errors in one study site.
Although the empirical approach has been successful and is currently used by commercial companies in many operational applications [30], the method has its limitations.In the presence of a large number of candidate predictors, with many of them being strongly correlated, there is a potential risk of multicollinearity problems that complicates the selection of model predictors [31] (Chapter 7).Furthermore, the empirical models are to some degree often affected by local effects such as geographical region, forest type, ALS acquisition parameters, and forest inventory design [5].In this context, the semi-empirical approach has the advantage of simplifying the modelling step and enabling the re-use and re-calibration of the models with new ALS data [21].
The use of predictive models for supporting forest inventories is traditionally done within either a design-based or a model-based inferential framework.Even though the design-based framework has been predominant, the case for the model-based framework is now being made [32].One of the main advantages of the model-based framework is that it does not, as opposed to design-based, rely on a probabilistic field sample.Instead, the inference is based on the model as a valid model for the distribution of Y i random variables.The sampled values y i are considered as a realization of the model [33] (p.40).The model cannot be observed, but the parameters of the model can be estimated from the sample.Furthermore, the model-based framework provides more flexibility in terms of model transferability and small area estimation [32].The variability of the predictive model can be evaluated by the mean square error (MSE).However, when used within the model-based inferential framework, the error of estimates for predictive models with small MSE can be larger compared to other candidate models that are nearly as good as the best one according to the MSE criterion (e.g., [12]).Evaluation of predictive models within the model-based inferential framework can instead be assessed in terms of (1) model specification to ensure that the model is correctly specified and (2) a model-based variance estimate (e.g., [34]).In order to evaluate the errors of predictive models, we studied a criterion that incorporates both a cross-validated MSE (i.e., mean square difference (MSD)) and a model-based variance estimate.The proposed evaluation criterion is referred to as the model error (ME).
While both the empirical and semi-empirical modelling approaches are used for AGB estimation, to the best of our knowledge, only the study by Hollaus et al. [5] has attempted to compare the two alternatives.The study by Hollaus et al. [5] found that the two approaches resulted in similar accuracy, and advocated the use of the semi-empirical approach based on its simplicity.Because the two approaches are commonly used, and with little scientific basis for choosing one over the other, it is of interest to assess possible differences in performance.It is also of interest to generalize the choice of approach, and the diverse data available for this study provided an opportunity for such generalizations.
The main objective of the present study was therefore to compare the two approaches in a variety of biomes.We used both MSD and ME as model evaluation criteria, and compared the performance of two empirical and two semi-empirical modelling approaches.The analyses were performed on five datasets representing four different forest biomes: tropical moist and -dry, temperate, and boreal forests.

Field Data
Data from five study sites from four different forest biomes were used in the present study (Table 1).These datasets containing 1182 field plots in total allowed for an extensive comparison of both the empirical and semi-empirical approaches.
The first dataset (S1) was collected in Amani Nature Reserve, northeastern Tanzania (5 • 08 S, 38 • 37 E).It represents a tropical moist forest biome.The area receives around 2000 mm rainfall per year, and most of the rain falls in the two wet seasons, April-May and October-November.Daily mean temperatures vary from about 16 to 25 • C. The area is covered by both old-growth natural tropical forest with a multi-layered canopy and previously harvested areas with single layered, light demanding, pioneer species.Forest inventory data such as diameter at breast height (DBH), tree height (H), and tree species were collected on 153 rectangular plots of about 50 m × 20 m in size over a period of four years, 2008-2012 [35].The horizontal area of the plots varies from 639 to 1239 m 2 , with a mean of 914 m 2 (Table 1).The variation in plot size is due to the procedure of plot establishment where the plots were laid out along the terrain slope, without slope correction.A threshold of ≥10 cm was used for recording the DBH.The corner coordinates of each plot were established by means of differential global navigation satellite system (GNSS), using survey-grade receivers.Errors in the x-y coordinates of the plot corners were estimated to an average of 0.57 m based on random errors reported from the post-processing software [36] and empirical experience of the relationship between reported error and the true error documented by Naesset [37].H was measured using Vertex IV hypsometer and a DBH-H model [35] (Equation ( 1)) was developed from the sample trees with height measurements collected on each plot and used to predict the H of all trees.AGB predicted for each individual tree in the sample was obtained using a local AGB model with DBH and H as independent variables [38] (p.43).Further details about the field data can be found in Hansen et al. [35].
The second dataset (S2) was collected in a tropical dry forest during 2011 and the first quarter of 2012 [39].The site is located in Liwale district (9 • 54 S, 37 • 38 E), southeastern Tanzania, and it is characterized by a climate with two rain seasons during November-January and March-May.The mean annual rainfall ranges between 600 and 1000 mm, and daily mean temperatures range between 20 and 30 • C. The forest structure has high diversity both in terms of structure because of varying soil conditions as well as a high species diversity.Dominant species included Brachystegia sp., Julbernadia sp., and Pterocarpus angolensis.Field plots were laid out in a clustered design as part of Tanzania's National Forestry Resources Monitoring and Assessment (NAFORMA) program [40].To avoid problems with spatial auto-correlation, we selected the two plots farthest away from each other from each cluster for the present study.The minimum distance between plots was 905 m, and the average distance within a cluster was 1438 m.The semi-variance of predicted wood volume studied during the design of the NAFORMA sample survey is indicated to level off at a distance of around 250 m for most of the forest types in Tanzania [41].Wood volume and AGB are strongly correlated variables.To test for an effect of the clustering, a log-likelihood ratio test between GLM models with and without a compound symmetry covariance structure that should account for heteroscedasticity was performed.The models were fit using maximum likelihood using the "nlme" [42] package in R [27].The test showed no significant cluster effect.Collection of forest inventory data was performed on concentric circular plots of 1, 5, 10, and 15 m in radius, following the field protocol of NAFORMA [43].Minimum thresholds for DBH for registration on the concentric plots were ≥1, ≥5, ≥10, and ≥20 cm respectively.The center coordinate of each plot was established with the same equipment and procedure as in the first dataset.Errors in the x-y coordinates of the plot corners were estimated to an average of 0.19 m [39].The height of every fifth three was measured using a Suunto hypsometer and DBH-H models were developed for each stratum.These DBH-H models were used to predict the H of all trees.Tree AGB was predicted using allometric models developed for Tanzanian woodlands by Mugasha et al. [44].Further details about the field data can be found in Mauya et al. [39].
The third dataset (S3) was collected in a temperate forest in the Czech Republic during January, July, and August of 2015 [45].The study site was located in Mendel University Training Forest, Křtiny (49 • 17 N, 16 • 44 E).Mean annual rainfall is 600 mm and mean annual temperature is 7.5 • C. Tree species composition can be described as coniferous dominant with Picea abies, Larix decidua, and Pinus sylvestris making up 76%, 4%, and 2% of the trees, respectively.The rest of the trees are Fagus sylvatica (17%) and other species (2%).The canopy is generally closed and with only one layer.Forest inventory data were collected on 50 circular plots with radius 12.62 m.The plots were selected by stratified sampling using an existing forest management plan.Plot centers were determined by GNSS using Topcon Hiper Pro with applied RTK corrections from the CZEPOS (Czech Positioning System) reference network.GNSS measurements were performed at five-second intervals for at least 20 minutes.On each plot, DBH and H of all trees above a minimum DBH ≥7 cm were measured.H was measured using a Vertex laser hypsometer.The age of forest stands was in the range of 60-130 years.Tree specific AGB was calculated using allometric models (Alnus glunitosa, Fagus sylvatica, Picea abies, and Pseudotsuga menziesii [46]; Pinus sylvestris and Abies alba [47]; Larix decidua [48]).
The fourth dataset (S4) was collected in a boreal forest located in Aurskog-Høland municipality (59 • 50 N, 11 • 30 E), southeastern Norway, during the fall of 2006 [49].Conifers (Picea abies and Pinus sylvestris) dominate the forest.Some deciduous species are also present, mainly Betula pubescens, especially in younger stands.A total of 201 circular field plots of 200 m 2 in size were systematically distributed.Plot centers were determined by GNSS using Topcon LegacyE receivers observing pseudo-range and carrier phase of the global positioning system and global navigation satellite system.All trees on the field plots with DBH ≥4 cm were callipered and DBH was recorded in two cm diameter classes.An average of nine trees per plot were sampled for height measurement with a probability proportional to stem basal area.H was measured with a Vertex hypsometer.Volume of each sample tree was estimated by means of species-specific volume models of individual trees dependent on DBH and H [50][51][52].A so-called tariff height for each tree on the plot using tariff height curves [53].For trees with a field measured H and an H estimated from the tariff height curves, separate volumes were estimated using both heights.The ratio between the mean volume predicted from field measurements and the volume for the sample trees was then used to adjust the volumes of all trees.Using the ratio estimated volumes, H were predicted for trees without a field measured H using single tree volume models [50][51][52].Trees with DBH <5 cm and H >1.3 m were counted and their H estimated by means of models presented by Tomter [54].AGB of each component; stump, stem, bark, dead and living branches, and foliage of each tree was calculated using species-specific allometric models developed by Marklund [55] with DBH and field measured or predicted H as independent variables.
The fifth dataset (S5) was collected in Hedmark county (61 as part of the Norwegian national forest inventory.A total of 648 circular field plots of 250 m 2 in size were systematically distributed.Plot centers were determined by GNSS using Topcon LegacyE receivers observing pseudo-range and carrier phase of the global positioning system and global navigation satellite system for a minimum of 15 min.Forest conditions and dominant species in S5 are relatively similar to the forest in S4.The AGB density, however, is considerably lower (Table 1).All trees on the field plots with DBH ≥5 cm were callipered.An average of ten trees per plot were sampled for height measurement with a probability proportional to stem basal area.H was measured with a Vertex hypsometer.AGB for trees with DBH ≥5 cm was estimated using the same procedure as in S4.Trees with DBH <5 cm and H >1.3 m were counted and their H estimated by means of models presented by Tomter [54].The original dataset covered both forest and non-forest areas.In the present study, only plots with registered AGB were used.
For all datasets, the AGB predictions of individual trees were summed for each plot and expanded to per hectare unit (Table 1).

ALS Data Acquisition and Initial Processing
The ALS datasets were acquired at different times and with different sensors and acquisition parameters.The acquisition parameters are summarized in Table 2. Post flight processing of the ALS data was performed using the TerraScan software (Terrasolid Ltd., Helsinki, Finland) [57].A triangulated irregular network (TIN) surface was constructed using the algorithm presented by Axelsson [58].Following the construction of the TIN, the elevation of each ALS echo relative to the TIN was computed, resulting in an elevation above the ground surface for each echo.

ALS-Derived Predictor Variables
The Leica ALS70 sensors provided records of up to five echoes per pulse and the Optech ALTM 3100 sensor four echoes.In the present study, we categorized each echo as "single", "first of many", "intermediate", or "last of many" based on the echo sequence.The "single" and "first of many" echoes were combined into one dataset and denoted as "first" while "single" and "last of many" were combined into another dataset and denoted as "last".The separation of "first" and "last" echoes were based on the assumption that the echoes from a pulse that has already been intercepted by the canopy ("last") provide divergent information from the echoes not previously intercepted ("first"), and is a common procedure for processing ALS data [59].Predictor variables were derived from the vertical distribution of echoes for each plot.The variables used in the empirical approach included percentiles of echo elevation, relative canopy density, and mean (h mean ) and maximum (h max ) elevation.Height percentiles were calculated at 10% intervals (h 10 , h 20 , . . ., h 90 ).These variables were calculated using a canopy threshold of 1.3 m, frequently used to eliminate echoes from undergrowth and rocks.Canopy density was computed by dividing the height from the canopy threshold to the 95th percentile height into ten equal intervals as recommended by Naesset and Gobakken [60].Cumulative proportions of echoes in the ten intervals to total number of echoes were calculated (d 0 , d 1 , . . ., d 9 ).Height percentiles and canopy density variables were calculated from the "first" and "last" data separately.
Studies following a semi-empirical approach have often used only one [5,26] or two [3,21] ALS-derived variables to model AGB or timber volumes.A variable capturing both the height and density of the canopy is often used.By calculating the height of ALS echoes without a canopy threshold, the height variables describe both the height and density of the canopy.Hollaus et al. [5] used a variable derived from "first" echoes only for modelling timber volumes in Austria, whereas others have used both a variable from all echoes [26] and a variable from first echoes [61] to model AGB in tropical biomes.Naesset [3] used the mean height of all echoes together with a canopy density variable to model timber volumes in eastern Norway.Furthermore, Magnussen et al. [21] found that complementing the canopy height of "first" echoes with a measure for the variance in canopy height reduced the RMSE with 11% when modelling timber volumes in eastern Norway.For simplistic reasons, and the ability to compute the heteroscedasticity consistent covariance matrix estimator described in Section 2.4.4,we chose to explore the approach used by e.g., Asner et al. [26] and Asner and Mascaro [61].We therefore derived two variables from the ALS data to be used in the semi-empirical approach.Although the calculations are not identical, we use the terms of Asner et al. [26] and derived the mean canopy height (MCH) from all ALS echoes and the top canopy height (TCH) from ALS echoes in the "first" dataset.

Statistical Modelling
As described in the introduction, there exists a number of statistical methods which have been used in both empirical and semi-empirical approaches.We have chosen two common approaches in operational settings and research i.e., OLS used in the empirical approach (described in Section 2.4.1), and NLS used in the semi-empirical approach (described in Section 2.4.2).

Ordinary Least Square Modelling
Multiple regression models expressed as LOG functions are frequently used for estimating AGB with ALS data in the empirical approach (e.g., [60,62]).Another common transformation used for predictive AGB models is the SQRT (e.g., [11,12]).Both transformations are applied to improve the linear relationship between the response and the predictors and to mitigate heteroscedasticity in the models.Thus, the models used for the empirical approach were formulated as: and sqrt where y i is field values of AGB in plot i, β 0 is the model intercept coefficient, e i is the model residual for plot i, and β is the vector of model parameters associated with the X matrix of ALS predictors.Residuals were assumed to be normally distributed with mean zero and a constant variance e i ∼ N 0, σ 2 i .Selection of predictor variables was performed using a best subset regression procedure implemented in the "leaps" package [63] in R, constrained to include a maximum of three predictors in the models.To avoid overfitting and multicollinearity, the models were selected using the Bayesian information criterion and variance inflation factors were kept below 10.Models without cross validation (see Section 2.4.3) were reported for assessing selected predictors.
Transformation of the response introduce a bias when the predicted biomass is back-transformed to arithmetic scale.To correct for bias in the LOG models, the correction factor for the uniform minimum variance unbiased estimator [14] was applied.Initial bias correction showed that the correction factor presented by Bradu and Mundlak [14] gave significant bias (p > 0.05) at two study sites (S2 and S5).Thus, a correction factor presented by Snowdon [13] was applied for predictions in S2 and S5.SQRT models were back-transformed according to Gregoire et al. [17].In the present study, one LOG (Equation ( 1)) and one SQRT (Equation (2)) model was fit for each study site.These are referred to as OLS LOG and OLS SQRT , respectively.

Nonlinear Regression
Instead of relying on transformations and multiple variables, the choice of modelling technique in the semi-empirical approach is often to use NLS.The use of the NLS technique enables the model to be fitted non-linearly through successive approximations by which only initial starting values for the approximation have to be stated [64] (Chapter 8).With several parameters to be estimated, the selection of starting values becomes difficult.Taking the semi-empirical approach with only one predictor, selection of starting values is simplified.We fit two nonlinear power models using the standard "nls" procedure in R with models formulated as: where a and b are the model parameters to be estimated, and H is the MCH and TCH used in separate models denoted NLS MCH and NLS TCH , respectively.

Model Evaluation Criteria
We first assessed the models based on the mean square deviation (MSD, Equation ( 4)) produced by a 10-fold cross validation procedure with: where K is the total number of folds (10), n is the total number of observations, n k is the number of observations in the k-th fold (k = 1, 2, . . ., K), MSE k is the mean square error of AGB predictions in the k-th fold, MSD is the mean square deviation.Successful use of the MSD as a model criterion assumes a design-based estimation strategy as the criterion is based on the observations and predictions on the sample units.We therefore sought to construct a model evaluation criterion that could be more informative than the MSE in cases where the field sample is not a probability sample.This was done by incorporating a model-based variance estimate based on a 10-fold cross validation.This criterion is referred to as the model error (ME, Equation ( 5)).The ME criterion combines both an expression of the MSD (Equation ( 5), first term) and a model-based variance estimate (Equation ( 5), second term): where DF is the model degrees of freedom and v ar( μ) is a model-based variance estimate calculated from the cross validation as: where X k is a matrix of variables derived from the ALS data in the k-th fold, X k ≈ ∂f X k ; β k /∂β k is a matrix of the approximated mean variables of X k in the k-th fold, and V k is a matrix containing the weights of the observations in the diagonal cells and the error correlations in the off diagonal cells in the k-th fold.

Covariance Matrix Estimators
A central part of the variance estimator in Equation 6is the matrix of estimated covariances between the parameter estimates β of the regression model.If the model errors are homoscedastic (V = σ 2 I), Equation ( 6) simplifies to: In the case of AGB modelling however, the variance of the response is likely to be related to the mean of the response, and thus we often have heteroscedasticity in our regression models.In this case, Long and Ervin [65] recommended using a heteroscedasticity consistent covariance matrix.OLS models with significant (p < 0.05) heteroscedasticity were identified using the Breusch-Pagan test [66].Residual plots were used to visually assess heteroscedasticity in the NLS models.In the presence of heteroscedasticity, heteroscedasticity consistent covariance matrix estimators were used as estimators of var( β).For the OLS models, heteroscedasticity consistent covariance matrix estimators of type HC 3 , presented by MacKinnon and White [67], were computed using the "sandwich" package [68] in R. The HC 3 (Equation ( 8)) estimator is frequently used and recommended for small samples [65].Computation of the HC 3 estimator requires the model projection matrix.Since this is not available for NLS models, a nonlinear heteroscedasticity consistent covariance matrix estimator described by White [69] (p.821) was therefore adopted for a nonlinear heteroscedasticity consistent covariance matrix estimator (NHC), Equation (9).The nonlinear NHC estimator was implemented in R and used to compute the final variance estimate in the presence of heteroscedasticity: where Z ≈ ∂f(X; β)/∂β. is a matrix of the approximated partial derivatives of the NLS model, ê2 i is the estimated residuals squared, and p i . is the diagonal elements of the projection matrix.

Discussion
Both empirical and semi-empirical approaches to modelling are common practice.Even so, to our knowledge, only one study has been published that compares the two approaches for biomass estimation purposes using ALS data.The study by Hollaus et al. [5] found that the performances of an empirical and a semi-empirical approach were similar in terms of the coefficient of determination.Hollaus et al. [5], somewhat unexpectedly, found that the MSE and standard deviation favored the semi-empirical approach.When assessing only MSD in the present study, we made similar observations for two of the five study sites (S2 and S5).However, a model-based variance estimate can be used in addition to the MSD to select a model that not only has a small estimated prediction error (MSD), but also a high estimated stability in the model parameters (slope and intercept), i.e., that model parameters are less influenced by new observations as indicated by the model-based variance estimate.The ME proposed in the present study is a novel combination of the MSD and a model-based variance estimator.Scatter plots of the two components of the ME (Figure 4) can be used to visualize the two components and their respective magnitudes.This can aid in deciding which estimators have the desired properties in terms of small MSD, small model-based variance, or a combination of the two.Using the proposed ME criterion, the empirical approach was found to result in smaller errors in the boreal and temperate biomes.In the tropical biome, the semi-empirical approach resulted in smaller ME.
Comparisons of OLSLOG and OLSSQRT models show that the OLSSQRT models resulted in smaller MSD and ME.In three of the study sites, the difference in MSD was quite small (<8%).In S3 and S5 however, the OLSSQRT models resulted in MSD values that were 48% and 25% smaller than in the OLSLOG models.This is in contrast to Hansen et al. [35] who reported approximately 7% smaller MSD using logarithmic transformation compared to SQRT.However, only the response was transformed in Hansen et al. [35].Other comparative studies have however also found advantages of using SQRT models compared to LOG [12,70].
Comparisons of NLSMCH and NLSTCH using MCH and TCH respectively as predictor variables showed an advantage of using MCH in all study sites except for S3 where NLSTCH resulted in 20% smaller MSD compared to NLSMCH.Although Asner and Mascaro [61] argue that TCH is a more consistent predictor, less affected by sensor effects compared to MCH, our results suggest that using

Results
Two OLS models and two NLS models were fit separately for the five study sites.The empirical approach was assessed by the two OLS models formulated as one LOG and one SQRT model, referred to as OLS LOG and OLS SQRT respectively.Models for the five study sites, fit using all observations, were reported for reference (see Supplementary Table S1).The semi-empirical approach was assessed in terms of the two NLS models with MCH and TCH as predictor variables, referred to as NLS MCH and NLS TCH respectively.Comparison of the model performances showed that the empirical approach resulted in the smallest MSD in three of the five study sites (Figure 1).In S1, S3 and S5, the empirical approach with a SQRT model resulted in the smallest MSD.In S3, the semi-empirical (NLS TCH ) and empirical (OLS SQRT ) approaches were almost equal in terms of MSD, with the OLS SQRT resulting in slightly smaller MSD.In terms of model-based estimated variance (Figure 2), the empirical approach resulted in the smallest variance estimates in four of the five study sites.In S1, the NLS MCH resulted in the smallest variance estimate.Considering the ME, the empirical approach had smaller ME in three of five study sites (Figure 3).In the two tropical study sites (S1 and S2), the NLS MCH model produced a smaller ME estimate compared to the other models.The scatterplot in Figure 4 provides additional information on the ME criterion.For sites S1, S2, S3, and S5, the first term of Equation (4) (i.e., an expression of MSD) was smallest for the NLS models.The second term of Equation (4) (i.e., estimated variance) however was larger, and the OLS models resulted in smaller ME for S3, S4, and S5.
When comparing the two different transformation strategies in the empirical approach, the OLS SQRT models resulted in smaller MSD in all study sites (Figure 1).It also resulted in smaller ME in all study sites, compared to LOG models (Figure 3).For the two NLS modelling approaches, NLS MCH resulted in smaller MSD and ME values in all study sites except S3 where the NLS TCH resulted in the smallest values of MSD and ME.

Discussion
Both empirical and semi-empirical approaches to modelling are common practice.Even so, to our knowledge, only one study has been published that compares the two approaches for biomass estimation purposes using ALS data.The study by Hollaus et al. [5] found that the performances of an empirical and a semi-empirical approach were similar in terms of the coefficient of determination.Hollaus et al. [5], somewhat unexpectedly, found that the MSE and standard deviation favored the semi-empirical approach.When assessing only MSD in the present study, we made similar observations for two of the five study sites (S2 and S5).However, a model-based variance estimate can be used in addition to the MSD to select a model that not only has a small estimated prediction error (MSD), but also a high estimated stability in the model parameters (slope and intercept), i.e., that model parameters are less influenced by new observations as indicated by the model-based variance estimate.The ME proposed in the present study is a novel combination of the MSD and a model-based variance estimator.Scatter plots of the two components of the ME (Figure 4) can be used to visualize the two components and their respective magnitudes.This can aid in deciding which estimators have the desired properties in terms of small MSD, small model-based variance, or a combination of the two.Using the proposed ME criterion, the empirical approach was found to result in smaller errors in the boreal and temperate biomes.In the tropical biome, the semi-empirical approach resulted in smaller ME.
Comparisons of OLS LOG and OLS SQRT models show that the OLS SQRT models resulted in smaller MSD and ME.In three of the study sites, the difference in MSD was quite small (<8%).In S3 and S5 however, the OLS SQRT models resulted in MSD values that were 48% and 25% smaller than in the OLS LOG models.This is in contrast to Hansen et al. [35] who reported approximately 7% smaller MSD using logarithmic transformation compared to SQRT.However, only the response was transformed in Hansen et al. [35].Other comparative studies have however also found advantages of using SQRT models compared to LOG [12,70].
Comparisons of NLS MCH and NLS TCH using MCH and TCH respectively as predictor variables showed an advantage of using MCH in all study sites except for S3 where NLS TCH resulted in 20% smaller MSD compared to NLS MCH .Although Asner and Mascaro [61] argue that TCH is a more consistent predictor, less affected by sensor effects compared to MCH, our results suggest that using TCH instead of MCH comes at a cost of loss in prediction accuracy.A possible explanation could be that the MCH captures variations in the forest canopy structure better, compared to TCH.This is supported by the result for S3 where NLS TCH resulted in smaller MSD and ME compared to NLS MCH .S3 is characterized by low variation in age and AGB (Table 1).
The types of statistical methods used to model the relationship between forest attributes and ALS data have possibly affected the results in the present study.Several other methods for modelling are available, most notably GLM, nearest neighbors techniques, neural networks, and the increasingly popular Random Forest algorithm.To assess these alternative methods is a great undertaking that we considered to be outside the scope of the present study.Instead, we chose to focus on two commonly used modelling methods: OLS and NLS.Furthermore, several of the non-parametric methods mentioned do not have available model-based variance estimators, and would require a different approach to produce a criterion similar to the ME.
Even though the empirical approach resulted in smaller ME in boreal and temperate study sites, the semi-empirical approach could be a viable option based on advantageous properties such as re-use or re-calibration of models in new study sites [21].Models that are relatively unaffected by noise, and that have a robust relationship between the response and predictor variable facilitate re-use or re-calibration.The model-based variance estimate, incorporated in the ME, could be calculated for new study sites, provided available ALS data.Thus, the proposed ME criterion could aid in the decision of modelling approach by providing a means of comparing the performance in terms of estimated variance to the estimated mean square error.

Figure 1 .
Figure 1.Calculated mean square deviation (MSD) for S1, S2, S3, S4, and S5.OLS LOG and OLS SQRT represent the empirical approach to modelling.NLS MCH and NLS TCH represent the semi-empirical approach to modelling.

Figure 2 .
Figure 2. Calculated model-based variance estimates ( v ar( μ)) for S1, S2, S3, S4, and S5.OLS LOG and OLS SQRT represent the empirical approach to modelling.NLS MCH and NLS TCH represent the semi-empirical approach to modelling.

Figure 3 .
Figure 3. Calculated model error (ME) for S1, S2, S3, S4, and S5.OLS LOG and OLS SQRT represent the empirical approach to modelling.NLS MCH and NLS TCH represent the semi-empirical approach to modelling.

Figure 4 .
Figure 4. Scatter plot of the two components of the ME for S1, S2, S3, S4, and S5.FT and ST are the first and second terms of Equation (5), respectively.

Figure 4 .
Figure 4. Scatter plot of the two components of the ME for S1, S2, S3, S4, and S5.FT and ST are the first and second terms of Equation (5), respectively.

Table 1 .
Summary statistics for the five study sites.