Introduction 

Precision forestry is becoming increasingly important in the face of the urgent challenges posed by climate change and the rising call for sustainable forest management (SFM) practices ([22]). This emerging trend depends on leveraging precise data obtained using advanced technologies to guide well-informed decision-making ([15]). In the realm of SFM, a noteworthy potential exists for value creation through implementing improved practices. Besides the ecological advantages related to enhanced productivity, which contribute to carbon sequestration and alleviate pressure on forests, the adoption of this paradigm shift is fundamentally linked with substantial economic and social value ([10]).

One concrete and increasingly prevalent instance of the potential being realized is precision harvesting (PH), which aims to maximize the efficiency of the harvesting operations while minimizing the impact on the environment and overcoming several limitations inherent to conventional logging practices (e.g., lack of precision and inefficient resource utilization due to lack of tree-level data - [37]). Within the scope of PH, the optimal allocation of assortments is a crucial component of the wood products supply chain and the carbon stock projections since different wood products hold different economical values and carbon storage potential ([5], [32]). Stem taper equations, which predict the change in stem form from ground to tip can accurately disaggregate trees into specific products based on certain specifications like log length and diameter and can help reach a better optimization of wood products ([6], [7], [43]).

In the academic literature, a panoply of taper model formulations and methods of parameter estimation is currently in use and the selection of a proper model form is more important than the actual fitting method ([49]). Following the classification of McTague & Weiskittel ([34]), taper models could be grouped into three main categories: (i) parametric models (e.g., simple-taper, variable-exponent, and segmented equations) which include most taper equations, are calibrated by the means of parametric approaches (e.g., nonlinear least squares - NLS - and nonlinear mixed effects - NLME - methods) and can ensure biologically consistent behavior ([35]). The Max and Burkhart model stands as an instance of this class, particularly within the segmented equations category, and has been demonstrated to exhibit high accuracy in predicting stem taper ([17], [23], [33]). (ii) Semiparametric models (e.g., B-splines and P-splines) which offer greater flexibility in the fit without requiring an extensive addition of parameters, may provide a better representation for complex form species ([26]). The B-spline regression model is a particular instance of semi-parametric model that was, first introduced by Kublin et al. ([24]), and has been successfully used to predict stem taper ([14], [19]). (iii) Nonparametric models (e.g., Random Forest - RF - and Artificial Neural Network -ANN) can offer a strong predictive performance without requiring the testing of statistical hypotheses ([1]). In particular, RF, an ensemble learning method in machine learning, has been gaining popularity in modeling taper by providing precise predictions and demonstrating generalizable performance ([36], [46]).

Developing taper equations requires multiple observations on upper-stem diameters and other covariates collected on standing or felled trees ([18]). Precision forestry in general offers substantial opportunities that may be beneficial for the forest sector including the development of decision tools such as taper equations. The advancements in Terrestrial Laser Scanning (TLS) technology present an unprecedented opportunity for acquiring essential data with high precision and non-destructive methods, thereby significantly reducing the time required compared to conventional approaches ([6], [42]). Data acquisition for TLS-based forest missions can be performed following single or multi-scan approaches. The single-scan approach has the simplest data acquisition setting and the fastest speed, however, a major problem with that scan method is that only parts of the trees are covered due to occlusion effects by other objects. On the other hand, the multi-scan approach appears the most accurate technique and has the potential to fully cover trees within the sampling area, since scans are performed from multiple directions to overcome the occlusion issue ([29]). Many studies using the multi-scan mode reported a stem detection rate between 62.1% and 100% depending on the forest structure and the scanning setup ([28], [31]).

Pinus nigra Arn. is a fast-growing conifer with a wide though fragmented distribution across Europe and Asia Minor, predominantly in mountainous areas. Due to its ecological flexibility, it is one of the most widely used tree species for reforestation worldwide, and it is considered a potential substitute for native coniferous species in Central Europe under future climate change ([48]). Economically, it is one of the most important conifers in Southern Europe, as its wood is highly suitable for general construction, indoor flooring, furniture industry, fuelwood, and paper pulp ([13]). In Italy, Pinus nigra Arn. stands cover an area of 444.785 hectares, representing 4.25% of the total national forest area, with age ranging from 50 to 95 years ([32]). The previously developed yield equations for Pinus nigra Arn. were based on data collected from young trees and have been shown to underestimate the volume of mature trees in the current conditions ([2]). Therefore, the primary aim of this study was to compare the performance of four taper models from tree class categories in predicting diameter over bark (DOB) and total stem volume of Pinus nigra Arn. specifically using data collected from multi-scan approach with TLS. The developed tool will provide forest managers with more accurate estimates of volume and enable better optimization of assortments for industrial purposes.

  Materials and methods 

Study area and sampling protocol

The study was conducted in the Vallombrosa forest (43° 44′ N, 11° 34′ E), a biogenetic reserve located about 50 km east-southeast of Florence, Tuscany, Italy. The native vegetation of this forest is mainly represented by beech (Fagus sylvatica) at higher elevations, oak-hornbeam stands (Quercus spp. mixed with Carpinus betulus L. and Ostrya carpinifolia Scop.), and chestnut (Castanea sativa Mill.) at lower altitudes ([12]). As part of the main reforestation program that occurred after the World War II, the European black pine (Pinus nigra Arn.) was introduced in the extreme western side of this forest ([3]). Covering 11% of the total forest area, Pinus nigra Arn. stands were specifically studied here.

The Pinus nigra Arn. area was stratified as part of the sampling scheme based on two criteria: stand density and total tree height. Each criterion was divided into three levels to ensure optimal representation of forest stand variability, resulting in a total of nine strata. An airborne lidar flight was utilized to achieve this stratification, providing high-resolution data on forest structure and height. The ALS survey was carried out under leaf-on condition using a Eurocopter AS350 B3 equipped with a LiDAR RIEGL LMS-Q680i sensor. The flight height was 1000 m above ground level. Full-waveform LiDAR data were registered and discretized to a point density of 10 points m^-2. Standard procedures for pre-processing ALS data (e.g., outliers and noise removal, classification of ground/non-ground, and computation of height on the ground) were carried out with the LAStools software ([20]) to obtain ALS normalized point clouds. Subsequently, one scan center was randomly placed in each stratum; then, a multi-scan TLS mission was performed on each scan center to collect detailed information on forest characteristics (Fig. 1).

TLS data collection and processing

TLS data were acquired using a FARO Focus 3D × 130 (FARO Technologies Inc., Lake Mary, FL, USA). The instrument uses a phase-shift-based technology with a maximum range of 120 m and acquires data with an azimuth scan angle of 360°. It collects the x, y, and z coordinates, and the intensity of the laser returns with a scan-ranging noise of ± 1 mm.

Using the TLS device, the data acquisition was set to a reasonable scan configuration providing a good trade-off between a sufficient density and the required time for a single scan (resolution of 7.6 mm at 10 m and 1/5 4× overall quality), for a total of about 28 million pulses per scan. With the aim of reducing occlusion due to other obstacles or vegetation and ensuring sufficient coverage, scanning was performed in multiple positions, at least 8 scans per plot subjectively distributed based on the density and structure. Each of the trees in the stand was at least scanned from 3 positions and up to 12 white polystyrene registration spheres (14 cm diameter) were placed in each plot to aid in the digital registration of individual scans.

Individual scans were merged on the plot level using the automatic registration algorithm implemented in Trimble Real Works^® (TRW) software. The software automatically joins overlapping redundant points to create one seamless 3D point cloud suited for the analysis. Details of the operation (settings, criteria, and thresholds) performed by the software are not declared nor accessible. The process attained a very low plot-level registration error, and TRW achieved a high precision scan placement (mean = 2.33 mm, sd = 3.18 mm, max = 8.21 mm). Using the TRW tool, each pine tree stem was manually separated from the rest (i.e., from the ground, the crown, and other non-pine points). Diameter measurements were acquired at approximately 1-meter intervals through the fitting of cylinders to point cloud data (Fig. 2). To determine the volume, Huber’s formula was employed, which entails using the diameter measured at the midpoint of a log segment (eqn. 1):

where V is the log volume, H is the log length, and s_0.5 is the mid-section sectional area.

Detailed descriptive statistics are provided for the trees inventoried in this study (Tab. 1). Statistics are derived from stem attributes estimated from point cloud data collected by the TLS device. The diameters and total height of the stems were estimated using the method described earlier, while the stem volume was determined by applying the Huber formula to each of the 1m-length cylinders. In total, 219 trees were extracted from the inventoried area and were systematically allocated to 5 folds. Four folds were exclusively designated for model training, while the remaining fold was reserved for rigorous testing during each iteration of the cross-validation procedure.

Used models

The data acquired by TLS was used to fit 4 taper models from 3 different model categories, namely parametric, semi-parametric, and non-parametric. Thereafter, the fitted models were evaluated in terms of their ability to predict diameters at different heights and the total stem volume as well. In the following, the different used taper models, as well as the evaluation metrics, are explicitly introduced.

Max and Burkhart (1976)

The segmented stem profile model developed by Max & Burkhart ([33]) was used to represent the parametric model’s category. Numerous studies have shown that this model is accurate in estimating upper stem diameters of many species ([21], [38], [39]). The model consists of three quadratic functions representing the stem’s lower, middle, and upper parts grafted together at two points (a₁ and a₂) with continuous polynomials and continuous first derivatives at each point. In the present study the taper equation was constrained to pass through diameter at breast height (dbh) while it is by construction constrained for the total height (ht).

Max & Burkhart ([33]) taper equation can be expressed as follows (eqn. 2):

where y_ij = (d_ij / dbh_i)², dbh_i is the dbh of tree i, with i = [1, 2, …, N], N is the number of trees in the sample, d_ij is the i-th upper-stem diameter over bark (dob) at height h_ij on tree i, with j = [1, 2, …, n_i], n_i is the number of diameter measurements for tree i, z_ij = 1 - (h_ij/ht) is the complement of the relative height, a₁ and a₂ are the join points to be estimated from the data, I_k = 1 if z > a_i, and 0 otherwise with k = [1, 2], b_p’s are regression coefficients with p = [1, 2, 3, 4] and ε_ij is the residual error term associated to the prediction of y_ij.

To constrain the taper equation to pass through dbh at the breast height, only one regression coefficient is needed to be changed. In this case, the b₁ parameter was modified as suggested in Cao ([8]) and replaced by the modified b₁* which can be expressed as follows (eqn. 3):

where z_bh = 1 - (1.3/ht). To estimate the parameters of eqn. 2, we used the Levenberg-Marquardt algorithm, which is implemented within the “nls” function in the R stats package ([44]).

Mixed-effects model

The stem taper models may show a high degree of multicollinearity, autocorrelation, and heteroscedasticity, all of which can be addressed with appropriate statistical methods. Many studies reported a reduction in the within-tree autocorrelation when including a tree-level random effect ([16], [9]). Using a mixed model for example, the autocorrelation within the stem can be accounted for at each height measurement level.

In the mixed effect formulation, all parameters of eqn. 2 are expressed as fixed effect (population level coefficients) while some parameters may contain additional random effect (for each tree within the population). In the case of the Max & Burkhart ([33]) taper equation, different combinations of the parameters (b₁, b₂, b₃, b₄, a₁, a₂) could be a potential candidate for random parameters.

Here is an expression for the stem profile including the fixed and random effects (eqn. 4):

where a_1i and a_2i are join points to be estimated from the data for each tree i, I_k = 1 if z > a_ki, and 0 otherwise with k = [1, 2], and b_pi’s are the regression coefficients to be estimated for each tree i, with p = [1, 2, 3, 4].

To define the statistical properties of the terms in eqn. 4, its vectorial formulation presented in the following offers a more simplified notation and will be considered (eqn. 5):

where Y_i is a vector of predictions, X_i is a matrix of predictors, A_i is a matrix for the fixed effects, B_i is a matrix for random effects, B is a vector of fixed parameters, U_i is the vector for random effects and ε_i is a vector of random errors, with the assumptions that ε_i = N(0,R) and U_i = N(0,D) where R is the diagonal variance-covariance matrix of ε_i and D is a diagonal variance-covariance matrix of U_i.

Using the “saemix” package in R ([11]), the stochastic approximation expectation maximization algorithm was used to estimate the parameters of eqn. 4.

The volume equation derived from the integration of the Max and Burkhart taper equation can be expressed as follows (eqn. 6):

where z_y = 1 - (h_y/ht) and z_x = 1 - (h_x/ht), h_x and h_y being respectively the lower height of interest (m) and the upper height of interest, K is a conversion factor (= π/40.000), I_k = 1 if z_y > a_i, and 0 otherwise, J_k = 1 if z_x > a_i, and 0 otherwise, with k = [1, 2].

B-splines

Semiparametric models are a class of statistical models that provide a flexible compromise between the rigidity of fully parametric models and the complexity of fully non-parametric models. By avoiding rigid assumptions about functional form, semi-parametric models offer the potential to capture the underlying relationship between variables more accurately, particularly for unknown or complex functional forms, while parametric approaches can be used for estimating model coefficients. Basis-splines (B-splines) are a type of semiparametric model that possess favorable numerical properties. To construct a B-spline model, a fixed number of splines, also called polynomials, are connected at specific points known as knots while ensuring the condition of continuity of their second derivatives across the knots.

To formalize, here we consider that the population mean diameter for a given relative height (h_r) may be approximated by a B-spline function of degree p and can be expressed as follows (eqn. 7):

where f(hr) is the response of the model which renders the diameter value d_i at a given h_r, β_l’s are the parameters of the model, d₁ is the number of parameters in the B-spline function with the condition that d₁ ≤ k₁ + 1 + p; k₁ corresponds to the number of internal knots, B_l,p⁽¹⁾ (hr) is the B-spline basis function computed based on the “de Boor recurrence relation”.

Due to the nonlinear nature of stem tapering and the hierarchy of data acquired from sampled trees, semiparametric models can be calibrated more efficiently using mixed approaches, which can also consider the random deviation of an individual tree from the population average. Linear mixed effect models can be used for the calibration if knot values are fixed a priori, which has the advantage of being less computationally expensive and more numerically stable than nonlinear mixed models.

The B-spline based taper equation considering both fixed and random effects can be written as (eqn. 8):

where f(hr_ij) corresponds to the population mean response, g(hr_ij) represents the tree-specific deviation from f(hr_ij), ε_ij is the tree-specific residual error, d_ij corresponds to the j-th diameter from the i-th sample tree, hr_ij corresponds to the j-th relative height from the i-th sample tree, β_l’s are the fixed regression coefficient of the model, θ_i,l’s are the random effect parameters to be estimated for each tree, d₂ is the number of random parameters in the g function with the condition that d₂ ≤ k₂ + 1 + p; k₂ corresponds to the number of internal knots in the random effect part of the model, B_l,p⁽²⁾(hr_ij) is the B-spline basis function for the random effect part of the model.

The vectorial formulation of eqn. 8 can be expressed as follows (eqn. 9):

where Y_i is the vector of diameter predictions at the level of the i-th tree, X_i is the matrix of B-spine basis values for the fixed effect part of the model, Z_i is the matrix of B-spine basis values for the random effect part of the model, β is a vector of fixed parameters, θ_i is the vector of random effect parameters, ε_i is the vector of residual error with the assumption that ε_i = N(0, R) and θ_i = N(0, G) where R is the diagonal variance-covariance matrix of ε_i and G is a positive non-diagonal matrix of θ_i. For more details on this method, see Kublin et al. ([24]).

According to the suggestions by Kublin et al. ([24]). The parameters of eqn. 8 were estimated considering knots [0.0, 0.1, 0.75, 1.0] corresponding to the B-spline basis function B⁽¹⁾(hr) for the population mean function and [0.0, 0.1, 1.0] for B⁽²⁾(hr) representing the deviation from the population average. Hence, the regression model depends only on relative height. Measurements of dbh and further diameter measurements of a specific tree are merely used to calibrate the tree specific deviation from the average taper curve by estimating the corresponding random effects. To ensure diameter estimates of zero at tree top, the last spline in both B⁽¹⁾(hr) and B⁽²⁾(hr) are omitted (see Fig. S1 in Supplementary material for a graphical visualization and [24]).

The R package “TapeR” was used to calibrate the model by the Restricted Maximum Likelihood (REML) method (default in the internally used R package “nlme” - [41]) providing the dbh values for each sampling tree of training data. To estimate stem volumes, the implemented numerical integration method included in the “TapeR” package was used ([25]).

Random Forest

Nonparametric methods have been recently used for taper predictions providing reasonable performances ([36], [50]). The present study used Random Forest to estimate the diameter at different heights of the stem. This algorithm works by training many decision trees on random subsets of the features, then averaging out their predictions. Regarding the structure of the model, the diameter was set as the dependent variable, which were predicted by four independent predictors namely dbh, h, ht and h/ht (abbreviations as above). For fine-tuning the model, the grid search method was used to evaluate all possible combinations of hyperparameters (viz., the number of decision trees; the number of features to consider when looking for the best split; the minimum number of samples required to split an internal node) within the search space as well as the cross-validation method to prevent overfitting. The Python library “scikit-learn” (ver. 0.24.2) which includes an implementation of the Random Forest algorithm was used here ([40]).

Since nonparametric methods are not algebraically integrable, including the random forest algorithm, the Monte Carlo integration, which is a numerical integration approach, was self-implemented and used to predict the stem volume for each sample tree (Fig. 3). This technique relies on a random choice of points at which the integrand is evaluated. Monte Carlo integration could be expressed as follows (eqn. 10):

where h_a and h_b being respectively the lower and the upper heights of interest (m), K is the conversion factor (= π/40.000), h_i corresponds to the i-th height on the tree generated randomly from a uniform distribution U(h_a, h_b), i = [1, 2, …, N], and N is the number of samples generated with the distribution U.

Evaluation of the methods

Taper models should yield unbiased estimates with minimal variance in both diameter outside bark (dob) and stem volume. A single overall measure of bias or residual error of estimation calculated for dob and stem volume does not provide an exclusive indicator of the goodness in evaluating several taper equations for a set of trees of a given species. This holds true when calculating bias metrics at different parts of stems and for different trees can result in a zero or close-to-zero average given that bias can be positive or negative. Alternatively, the standard error which can report the variability of the biases is considered of better reliability. Moreover, when dealing with small sample size, the residual variance can be influenced by the degrees of freedom, making the correlation index, also known as the fitness index, a useful metric to address this concern. All those metrics were used to rank the fitted models.

The most reliable measure of mean bias (B), overall residual variance (SSE) and the fitness index (FI) can be expressed as follows (eqn. 11 - eqn. 13):

where Y_i is the actual observation, Ŷ_i is the predicted value of the actual observation, bar{Y}_i is the mean of the actual observations, n is the number of observation and k corresponds to the number of estimated parameters.

As part of our evaluation process, we employed a grouping methodology to assess the performance of the taper models in conjunction with the overall assessment presented earlier. For dob, biases and standard errors were calculated for each 10% step of the height along the stem, while for the total stem volume (from ground to top), the biases and standard errors were calculated by DBH classes. This allowed us to gain a comprehensive understanding of the accuracy of the taper equations across various tree dimensions as well as across different parts within the trees.

  Results 

Taper prediction

The model by Max & Burkhart ([33]) was fit using both fixed and mixed-effects approaches. Different mixed-effects combinations of parameters were tested (Tab. S1 in Supplementary material), and the combination involving a₁, a₂, b₁, b₂ and b₄ produced the lowest values of Akaike’s information criterion (AIC) and Bayesian information criterion (BIC), resulting in the following parameter values (Tab. 2) and model form (eqn. 14):

The mixed B-splines model was fit using the linear mixed-effects approach implemented in the TapeR package and reached convergence. In order to constrain the model predictions to equal zero at the tree top, d₁ which is the number of fixed parameters was set to 5, while d₂ which corresponds to the number of random parameters was set to 4, omitting the topmost spline, respectively (see Fig. S1 in Supplementary material). The results of B-splines model fitting are represented in Tab. 3.

Regarding the RF model, after rigorous fine-tuning, the model reached a state of convergence as the number of decision trees stabilized at an optimal value of 267. In addition, the hyperparameter tuning process led to the following optimized configuration: number of features to best split = 3; minimum samples to split internal nodes = 4.

For each taper equation, overall statistics of fit (Bias, SEE and FI) for the entire stem were calculated and are represented in Tab. 4 for dob by model. Among the four models under investigation, discernible trends in bias emerged. The Max and Burkhart models exhibited a noteworthy positive bias, suggestive of their propensity to overestimate diameters. In contrast, the B-spline and Random Forest models displayed a discernible negative bias, indicating a tendency to underestimate diameters. Notably, the Max and Burkhart fixed effect model demonstrated a particularly pronounced inclination towards overestimation. The results suggest that mixed-effect models outperformed models considering only fixed effects. In particular, the Max and Burkhart and B-splines mixed-effects models showed relatively comparable SEE values of 0.78 cm and 0.96 cm, respectively. In contrast, the Max and Burkhart fixed-effect and random forest models had comparable SEE values of 1.90 cm and 2.26 cm, respectively. This comparison holds true for FI as well, although to a lesser extent than for SEE.

To assess the performance of taper models at different positions throughout the stem, the statistics of fit, derived from the 5-fold cross-validation results, were analyzed by relative height classes. These statistics are documented in Tab. S2 (Supplementary material) and graphically depicted in Fig. 1. To comprehensively assess the taper of all models, the average bias, SEE, and FI were calculated for each model by relative height class along the stem. The results show that models including mixed effects outperformed those considering only fixed effect roughly along all the section heights. Noticeably, a significant increase of SEE was observed for all the models from a relative height class of (70-80%) and was shown to coincide with the crown base height (CBH) for most of the trees (Fig. 4).

Bias measurements indicate that the Max and Burkhart fixed-effect model tends to consistently overestimate upper stem diameters throughout the entire stem. On the other hand, the random forest model tends to underestimate upper stem diameters. B-spline-based predictions tend to overestimate upper stem diameters in the lower and upper portions while underestimating in the middle portion. In contrast, predictions from the Max and Burkhart mixed-effects models exhibit an opposing bias pattern in the lower and middle sections but demonstrate similar behavior in the upper portion of the stem.

Volume prediction

The predictions of volume (over bark) by four taper models were compared with each other. Overall statistics of fit (Bias, SEE, and FI) for the total stem volume are represented in Tab. 5. The results show that mixed-effects models outperformed models considering only fixed effect. Specifically, Max and Burkhart and B-splines models had a comparable SEE values of respectively 0.021 and 0.034 m³, while the Max and Burkhart fixed effect and the Random Forest models had comparable SEE values of respectively 0.221 and 0.273 m³. In addition, the overall results show that all the four models tend to generally overestimate the volume as the bias resulted in positive values for all the candidate models. Notably, the Max and Burkhart models (both fixed and mixed effects) exhibited the smallest overprediction, with an absolute bias value of less than 0.01 m³.

To assess the performance of stem volume prediction for different tree dimensions, the statistics of fit (SEE, bias, and FI) were analyzed by DBH class categories. These statistics are documented in Tab. S3 (Supplementary material) and graphically depicted in Fig. 6. In general, results show that prediction errors are larger in large diameter trees for both Max and Burkhart fixed effect and Random Forest models, while the prediction error tend to be constant over all the DBH classes for Max and Burkhart mixed effects and B-splines models. The results also demonstrate a strong correlation between volume predictions and volume observations (FI > ~ 0.5) for all four models across all DBH classes, except for the DBH class category (55-65 cm) and (> 65 cm), where both the M&B fixed effect and Random Forest models exhibit limitations (FI ~ 0.42 and 0.21, SEE ~ 0.39 m³ and 0.71 m³, respectively).

Model ranking

A taper-estimating system should provide an unbiased estimation with minimum variance of both diameter outside bark (dob) and total stem volume outside bark. Using the results from Tab. 4 and Tab. 5, it is possible to rank the four estimating systems. The rank sums were created by ranking the performance by dob and stem volume which were equally weighted here in this study. The rank sums reported in Tab. 6 were generated in such a way that each estimating system was assigned a rank separately for every relative height for dob and for every DBH class for the volume. These ranks were also summed for standard errors and biases. As a result, Max and Burkhart mixed-effects system ranked first in most of all the ranking categories. B-splines based system ranked second in estimating diameter and volume by standard error of estimates while random forest system ranked third and Max and Burkhart fixed effect ranked fourth with regards to the same previously discussed estimates.

  Discussion 

Taper equations are an indispensable tool in the field of forest management and industry. However, the development of an accurate taper model present inherent challenges, mainly due to the need for comprehensive customization to account for the characteristics exhibited by each species and site-specific environmental conditions, and the availability of high-quality data required for the model development. In the present study, four taper models from three different class categories were compared in terms of their capacity to predict diameter over bark (dob) and total stem volume of Pinus nigra Arn. Across models, the Max and Burkhart and B-splines mixed-effects models demonstrated superior accuracy and reliability in predicting both dob and total stem volume, with a small advance for the Max and Burkhart model.

The Terrestrial Laser Scanner (TLS) has previously been employed as an effective tool for measuring trees in forest environments and generating taper equations ([45], [47]). The case study conducted on black pine forests, which is presented here, demonstrated the feasibility of using TLS for surveys and highlighted the high quality of the acquired data. In the Italian context, where tree measurements for quantitative purposes have significantly decreased in recent decades, technological tools like TLS represent a pivotal advancement for rejuvenating our understanding of the Italian forest sector.

In the present study, biases and standard errors of estimate were evaluated for diameter outside bark (dob) and total stem volume. These biases and standard errors were also evaluated for different heights within the trees and across different tree sizes. The overall goodness-of-fit statistics (Bias, SEE, and FI) were calculated and reported in Tab. 4, Tab. 5 and Fig. 5, Fig. 6 for the four tested models. The results indicate that the Max and Burkhart and B-splines mixed-effects models explained more than 99% of the total variation in predicting the upper-stem diameter and total stem volume. The Max and Burkhart fixed-effect model explained more than 99% of this total variation, while the Random Forest model explained more than 96% for the diameter variation and 96% for the total volume variation.

The superiority of models that include random effects over those that only consider fixed effects has been demonstrated in various previous studies, and our findings are in agreement with this. For example, a study conducted by Leites & Robinson ([27]) showed that the inclusion of random effects considerably improved the fit of the Max and Burkhart taper equation for Pinus taeda. In a different study by Bronisz & Zasada ([4]), the Kozak’s taper equation fitted using a mixed-effects approach provided better results for the upper diameter and total stem volume for Pinus sylvestris. In our study, the Max & Burkhart ([33]) model was fitted using both fixed and mixed effects approaches. Notably, the error was reduced by 52% for the upper-stem diameter and by 89% for the total stem volume when random effects were additionally considered.

It is also noteworthy that the evaluation of the proposed models exhibited a higher Standard Error of Estimation (SEE) in predicting the upper-stem diameter at a relative height of [0.7, 0.8] (Tab. S2 - Supplementary material). The analysis of the distribution of the tree’s relative crown base height reveals that 50% of the observations are located between a value of 0.7 and 0.8, and the maximum relative crown base height was found to be approximately equal to 0.8 (Fig. 4). This implies that all observations above a relative height of 0.8 were necessarily measured inside the crown. This last result may partly explain the notable increase in error in predicting the upper-stem diameter at a relative height of 0.7. In fact, many studies have reported the inability of Terrestrial Laser Scanners to collect reliable data beyond the base of the crown due to interference from branches and leaves as well as the mutual occlusion between canopies. Practically, when approaching higher parts of the stem, occlusion starts to deteriorate the quality of cylinder fitting ([30], [42]). Additionally, as the distance from the scanner increases, spatial resolution decreases, which is a well-known effect in scanning that particularly impact the upper parts of trees ([29]). The interference from the crown, compounded with the increased distance, may explain the models’ inability to detect a well-defined taper pattern above the crown base.

In this study, Monte Carlo numerical integration was employed for estimating the volume of tree stems using the non-differentiable and complex Random Forest model. This decision was justified by the model’s inherent characteristics mentioned earlier, which render traditional integration methods less suitable. It is imperative to underscore the significance of comparing our chosen method with alternative numerical integration techniques to gain insights into result robustness and accuracy, potentially informing future research. Notably, the RF model, unlike the parametric and semi-parametric models used in this study, typically demands a substantial amount of training data. In our analysis, the RF model obtained the lowest score in terms of SEE (dob and volume). This highlights the potential impact of our relatively small dataset (219 trees) on RF performance. In addition, the use of various integration methods (Huber, Monte Carlo, and analytical integration) in our study may introduce minor biases affecting our results. Further research is needed to assess and potentially mitigate these biases, as their combined impact on the findings remains uncertain.

Considering the results of this study, the Max & Burkhart ([33]) taper model, calibrated by means of a mixed-effects approach, and its compatible volume integration could be considered for operational use in estimating diameter outside bark at different heights, as well as total stem volume, respectively, with an overall precision of 0.781 cm and 0.021 m³ in Pinus nigra Arn.

  Conclusion 

This study explored the application of precise and non-destructive techniques, such as Terrestrial Laser Scanning, for collecting data to develop taper equations specific to Pinus nigra Arn., in line with the principles of precision forestry. The collected data were used to fit four taper models from three different class categories, evaluated for their ability to predict diameter at various heights and total stem volume. The Max & Burkhart ([33]) taper model, calibrated with a mixed-effects approach, proved superior in estimating diameter and volume compared to other models. However, a common limitation among all models was their inability to accurately predict diameter beyond a relative height of 0.7 due to crown interference. The selected model, within the scope of this study, can be operationally utilized by forest managers to enhance predictions of diameter and stem volume, aiding in optimizing wood allocation to different harvested products. Future research perspectives may focus on refining algorithms and methods for extracting diameter at the upper parts of the stem.

  Data, code availability and setup 

The whole dataset is made available at ⇒ https:/­/­doi.­org/­10.­5281/­zenodo.­8414408. The study was performed on a laptop with Win10 22H2 operating system, an Intel Core i5 7200U CPU with maximal capacity of 2.70 GHz, and 8 Gb of RAM.

  Authors contribution 

NP and CV contributed equally to the manuscript.

  References

  Supplementary Material

Tab. S1 - Model selection statistics for evaluating the inclusion of random parameters in the Max and Burkhart taper model.
Fig. S1 - Population mean-fixed effects B-spline basis.
Tab. S2 - Bias (cm), standard error (cm), and fitness index of the estimate of dob by relative height for four models.
Tab. S3 - Bias (m³), standard error (m³), and fitness index of the estimate of stem volume by DBH class for four models.