Development of a Methodology Based on ALS Data and Diameter Distribution Simulation to Characterize Management Units at Tree Level

Abstract

Characterizing Management Units (MUs) with tree-level data is instrumental for a comprehensive understanding of forest structure and for providing information needed to support forest management decision-making. Airborne Laser Scanning (ALS) data may enhance this characterization. While some studies rely on Individual Tree Detection (ITD) methods using ALS data to estimate tree diameters within stands, these methods often face challenges when the goal is to characterize MUs in dense forests. This study proposes a methodology that simulates diameter distributions from LiDAR data using an Area-Based Approach (ABA) to overcome these limitations. Focusing on maritime pine (Pinus pinaster Ait.) MUs within a forest intervention zone in northern Portugal, the research initially assesses the suitability of two highly flexible Probability Density Functions (PDFs), Johnson’s SB and Weibull, for simulating diameter distribution in maritime pine stands in Portugal using the PINASTER database. The selected PDF is then used in conjunction with ABA to derive the variables needed for parameter recovery, enabling the simulation of diameter distributions within each MU. Monte Carlo Simulation (MCS) is applied to generate a sample list of tree diameters from the simulated distributions. The results indicate that this methodology is appropriate to estimate diameter distributions within maritime pine MUs by using ABA combined with Johnson’s SB and Weibull PDFs.

1. Introduction

In forestry, a management area generally includes a group of forest stands (Management Units, MUs) with varying ages and species compositions that are managed collectively [1]. These forest stands can be classified as either pure (containing a single species) or mixed-species and may consist of even-aged or uneven-aged structures. Given these differences, accurate forest assessment and management depend on the precise quantification of forest information at both the tree and stand levels [2], which is a critical task for sustainable forest management.

Characterizing Management Units at the tree level improves our understanding of forest structure and provides essential data for decision-making. However, detailed forest inventory over large areas presents challenges due to the significant time and resources required. To address this challenge, geotechnologies like Airborne Laser Scanning (ALS) have proven effective for the cost-efficient characterization of large forest areas, capturing the three-dimensional structure of vegetation and facilitating forest mapping, which enhances our understanding of forest dynamics [3,4].

One primary application of ALS data in forestry is the Area-Based Approach (ABA), which combines field inventory data with statistical and spatial generalizations of normalized ALS data [5]. This enables the efficient estimation and mapping of forest inventory attributes like dominant height, diameter, volume, basal area, and tree density at the landscape level [6]. While the ABA approach does not directly measure tree diameters, several stand-level variables can be estimated through predictive models. These variables can then be used to recover Probability Density Functions (PDFs). This approach allows for detailed vertical structure characterization and reproduction of different diameter distributions shapes [7].

Several Probability Density Functions, such as Johnson’s SB, Weibull, Beta, Gamma, Lognormal, and Truncated Weibull, have been explored for simulating diameter distributions in forest stands using field data [8,9,10]. While these functions have proven suitable for diameter distribution simulation, Johnson’s SB and Weibull have demonstrated particular flexibility and adaptability, making them effective for different forest types and conditions. For instance, Johnson’s SB accommodates a wide range of combinations of skewness squared and kurtosis, making it particularly effective for modeling diameter distributions in complex forest stands [8]. Conversely, the Weibull PDF is simpler, requiring fewer parameters and proving suitable for diverse forest types, particularly in managed forests [11].

In Mediterranean forests, Johnson’s SB and Weibull have consistently performed well in simulating diameter distributions, adapting effectively to different forest types. Páscoa (1987) [12] tested the parameter-recovery Weibull of the PDF to simulate the evolution of structure, growth, and yield of maritime pine stands in Portugal. Fonseca et al. [13] demonstrated the effectiveness of Johnson’s SB in simulating maritime pine diameter distributions in Portugal using a parameter recovery approach, highlighting its ability to closely match observed distributions. Similarly, Mateus and Tomé [14] successfully applied Johnson’s SB PDF to simulate diameter distribution in eucalyptus plantations in Portugal, confirming its suitability for these stands. In addition, Palahí et al. [15] compared the Johnson’s SB, Weibull, and Beta PDFs to simulate diameter distributions in several stands with different species compositions, including maritime pine.

Building on this, previous research has explored methods for simulating the diameter distribution of boreal forests using ALS metrics and Probability Density Functions, [16,17,18]. Similar studies have applied this approach to different forest types in the Iberian Peninsula. For instance, Arias-Rodil et al. [19] used the two-parameter recovery of the Weibull PDF to simulate diameter distribution in Pinus radiata stands in Galicia, Spain. Additionally, Cosenza et al. [20] tested the Johnson’s SB and Weibull PDFs for modeling diameter distribution in Eucalyptus globulus and Pinus radiata stands using ALS data and utilized ABA as a method to estimate the variables required for each parameter recovery.

According to the literature, many studies rely on Individual Tree Detection (ITD) methods to simulate diameter distributions within a stand [21,22,23,24], as this approach enables the identification of individual trees, with diameters estimated using tree height and crown attributes [25]. However, despite their widespread use, these methods can face significant limitations in dense forests, where overlapping canopies and high tree density often constrain the accurate detection and differentiation of individual trees, potentially leading to inaccuracies in diameter estimates. For instance, Vauhkonnen et al. found that tree density strongly affects the success of single-tree detection algorithms across different forest types [26]. ABA offers an alternative by enabling the wall-to-wall estimation of stand forest attributes through grid cells [27], which combined with the simulation of diameter distributions, facilitates the detailed characterization and mapping of within-stand variability.

Given these considerations, this study proposes a methodology to characterize the Management Units (MUs) at tree level in dense forests. First, we assess the suitability of two flexible Probability Density Functions (PDFs), Johnson’s SB and Weibull, for simulating diameter distributions in maritime pine stands. While Johnson’s SB and Weibull PDF have been tested for maritime pine in Portugal, they have not yet been compared, highlighting a gap that this study aims to address.

Secondly, this study applies the ABA to estimate the forest attributes inputs needed for the parameter recovery of the selected PDF, enabling diameter distribution simulation and tree list generation in each MU. This detailed characterization is crucial for forest management, as it supports the use of individual tree models, which are essential for simulating thinnings or managing complex stands, such as converting even-aged stands to uneven-aged or pure stands to mixed-species stands.

2. Materials and Methods

The flowchart in Figure 1 summarizes the entire process for generating a tree list for each Management Unit. The methodology is divided into three main stages. The first stage involves selecting a Probability Density Function (PDF) for parameter recovery, using a database based on permanent plots and silvicultural trials of maritime pine in Portugal (PINASTER), which is managed by the ForChange research group (https://www.isa.ulisboa.pt/cef/forchange/fctools/ (accessed on 16 January 2024)). This database was used to assess the ability of the highly flexible PDFs, such as Johnson SB and Weibull, in simulating diameter distribution in maritime pine stands. This assessment is conducted through the Kolmogorov–Smirnorv (KS) test and a comparison of the growing stock volume produced with each PDF in relation with the observed data.

The second stage estimates the stand variables required for parameter recovery of the selected PDF using the Area-Based Approach method, where predictive models are developed from forest field inventory and ALS data. The final stage applies the parameter recovery to simulate diameter distribution in each Management Unit and generate a detailed tree list.

2.1. Study Area

This study focuses on 214 maritime pine Management Units (MUs), which cover a total area of 16.15 km² within an aggregated management forest area in northern Portugal, known as Paiva and Entre Douro e Sousa (Figure 2). The maritime pine Management Units are surrounded by eucalyptus MUs, which are the dominant species in the landscape [28]. This region includes parts of the counties of Paredes, Penafiel, Gondomar, Paços de Ferreira, Lousada, Amarante and Marco de Canavezes, situated north of the Douro River, and extends to the border of Castelo de Paiva County and part of Arouca and Cinfães to the south. The area features diverse topography, with altitudes ranging from 67 to 775 m and slopes varying between 5 and 88 degrees.

2.2. Data Used

This study utilized two datasets. The first dataset, known as the PINASTER database, includes permanent plots and silvicultural trials installed in maritime pine stands across Portugal. It was used to assess the suitability of two highly flexible Probability Density Functions (PDFs) for simulating diameter distributions. The second dataset consists of field inventory plots data collected in 2022 within the study area. These data were combined with ALS data to predict stand variables needed for the PDFs parameter recovery in each Management Unit.

The PINASTER database compiles data from 12 silvicultural trials, each with multiple measurements over time. Some measurements were taken when stands were very young, with most trees shorter than 1.30 m. To ensure datasets where most of the trees had a size allowing dbh measurement, we selected only plots where at least 90% of the trees had measured diameters, resulting in 186 plots and a total of 981 plot measurements.

The second dataset was surveyed during a forest inventory conducted by AFVS—Associação Florestal Vale do Sousa, a Forest Owners Association. The survey took place between 5 July and 25 September 2022. A total of 30 maritime pine plots, each covering 500 m², were randomly located and surveyed according to the Portuguese National Forest Inventory Field Manual [29].

Table 1. Description of the field data with the values related to minimum, mean, maximum, and standard deviation (sd) values.

Dataset	Variable 1	Unit	Minimum	Mean	Maximum	sd
	${\mathrm{d}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$	cm	0.50	6.58	23.39	5.43
	${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$	cm	0.46	13.15	34.88	6.69
PINASTER	${\mathrm{d}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$	cm	1.50	21.28	53.00	9.13
	$\mathrm{d}\mathrm{g}$	cm	0.57	13.62	35.23	6.66
	G	m2 ha−1	0.03	22.07	64.07	16.34
	N	stems ha−1	244.00	1428.11	7755.00	904.39
	${\mathrm{d}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$	cm	5.00	10.21	25.00	5.10
	${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$	cm	9.70	21.85	43.66	8.15
Field inventory data	${\mathrm{d}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$	cm	16.00	37.41	60.00	11.94
	$\mathrm{d}\mathrm{g}$	cm	10.03	23.15	45.19	8.39
	G	m2 ha−1	2.97	19.77	52.34	11.20
	N	stems ha−1	120.00	534.66	1860.00	370.53

¹ ${\mathrm{d}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$: minimum diameter at breast height (dbh); ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$: mean dbh; ${\mathrm{d}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$: maximum dbh; $\mathrm{d}\mathrm{g}$: quadratic mean dbh; G: stand basal area; N: tree density.

summarizes both datasets.

2.3. LiDAR Data Characteristics and Pre-Processing

LiDAR data acquisition for the study area took place on 10 June 2022, as part of the FIRERES project [30]. Data collection was conducted using a PATERNAVIA P68C-TC aircraft (RIEGL, Vienna, Austrian) flying at an average altitude of 1600 m above ground level at a speed of 220 km h⁻¹. A RIEGL VQ-1560i LiDAR sensor (TOPCAD Ingeniería, Lugo, Spain) was used operating with a scan angle range of 60° and pulse repetition rate up to 2.0 MHz, resulting in a point density of 5 pts m⁻².

Following data acquisition, point cloud classification was conducted using the lidR package [31] within the R environment [32]. Once the ground returns were identified, they were filtered and triangulated to generate a Digital Terrain Model (DTM) with a 1 m resolution. This DTM was then used to normalize the point cloud, ensuring that all non-ground measurements were adjusted relative to ground elevation.

With the normalized point cloud in place, descriptive statistics (metrics) were calculated, as outlined in

Table 2. Descriptions of the metrics extracted from the normalized ALS data.

Metrics	Description
Zmean, Zmax	Mean and maximum height
Zsd, Zcv	Height standard deviation, height coefficient of variation
Ziq	Height interquartile range
Zskew, Zkurt	Skewness and kurtosis of height distribution
Zsqmean	Quadratic mean height
Zentrpy	Height entropy
Z5, Z10, Z15, Z20, Z25, Z30, Z35, Z40, Z45, Z50, Z55, Z60, Z65, Z70, Z75, Z80, Z85, Z90, Z95, Z98, Z99	Height percentile from 5% to 99%
CRR	$\mathrm{Canopy} \mathrm{relief} \mathrm{ratio} (Hmean-Hmin)/(Hmax-Hmin)$
Para2	Percentage of all returns above 2 m
Param	Percentage all returns above mean/Total all returns
CC	Percentage of first returns above 2 m
Pfram	Percentage of first returns above mean/Total all returns
ADD	Mean absolute deviation
Pzbvzmn	Percentage of returns above mean height (Zmean)
Zpcum1, Zpcum2, Zpcum3, Zpcum4, Zpcum5, Zpcum6, Zpcum7, Zpcum8, Zpcum9	Cumulative percentage (from 10% to 90%) of returns located in lower 10% maximum elevation

. These metrics were computed for both the entire study area, using a grid with a resolution of 20 × 20 m, and for the clipped normalized point cloud that overlapped with the field inventory plots.

2.4. Assessing the Suitability of Johnson’s SB and Weibull PDFs for Diameter Distributions Simulation in Maritime Pine Stands

To evaluate which Probability Density Function better represents diameter distribution for maritime pine, this section compares Johnson’s SB and Weibull distributions. The process involved testing each PDF on the maritime pine plots in the PINASTER database and comparing their adequacy.

For the Johnson’s SB PDF, as described by Johnson (1949) [33], four key parameters characterize the distribution: the location parameter $\epsilon$, the scale parameter $\mathsf{\lambda }$, and shape parameters $\mathsf{\gamma }$ and $\mathsf{\delta }$. These parameters collectively define the distribution’s characteristics, including the lower bound, range, and shape, respectively (Equation (1)).

$$f\left(x\right)={\displaystyle \frac{\mathsf{\delta }\mathsf{\lambda }}{\sqrt{2\mathsf{\pi }}\left(x-\epsilon \right)(\epsilon +\mathsf{\lambda }+x)}}exp\{-{\displaystyle \frac{1}{2}}{[\mathsf{\gamma }+\mathsf{\delta }\mathrm{l}\mathrm{n}(x-\epsilon /\mathsf{\epsilon }+\mathsf{\lambda }-x\left)\right]}^2\}$$

(1)

$$\epsilon <x<\epsilon +\lambda ,\delta >0,-\infty <\gamma <\infty , \lambda >0, and \epsilon \ge 0; f\left(x\right)=0, otherwise.$$

These parameters were estimated using a three-parameter recovery algorithm adapted from Parresol et al. [34], implemented in the R environment (version 4.3.0) [32] with the minpack.lm package [35]. The recovery process applies the Levenberg–Marquardt algorithm to solve complex nonlinear equations, ensuring an optimal fit for each plot in the dataset. The notation $D~SB (\lambda , \epsilon , \gamma , \delta )$indicates that the variable $D$ follows the SB distribution, characterized by these defining parameters.

A transformation can be applied to normalize the distribution; starting from a simple linear transformation in Equation (2), a normally distributed variable in Equation (3) is derived, which connects the transformed variable y to a standard normal distribution. This relationship is then used in Equation (4) to compute the non-centered moments of the distribution of y through an integral involving a logistic function and the normal probability density function. These transformations provided a detailed representation of the diameter distribution shape across the plots.

$$y=f\left(d\right)={\displaystyle \frac{d-\epsilon }{\lambda }}$$

(2)

$$z=\mathsf{\gamma }+\mathsf{\delta }\ln \left({\displaystyle \frac{y}{1-y}}\right)\sim N\left(\mathrm{0,1}\right)$$

(3)

$${\mathsf{\mu }}_r^{\prime }\left(y\right)={\displaystyle \frac{1}{\sqrt{2\mathsf{\pi }}}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left[1+e^{{\displaystyle \frac{\mathsf{\gamma }-z}{\mathsf{\delta }}}}\right]}^{-r}{e}^{-z^2/2}dz$$

(4)

The following equations provide the key parameters of the Johnson’s SB PDF, using stand-level average attributes to estimate the diameter distributions. Equation (5) defines the shape parameter $\gamma$, which adjusts the skewness of the distribution and can be expressed as a function of the other three parameters. Equation (6) describes the mean value $\overline{d}$ of the distribution as a function of the location parameter $\epsilon$ and the first non-centered moment ${\mathsf{\mu }}_1^{\prime }\left(y\right)$, scaled by $\lambda$, which defines the distribution’s range.

To initiate the estimation process, the initial value for the location parameter ε was set to zero, as sampling data was used [36]. The initial values for $\lambda$ and $\delta$ were set to $dma\mathrm{x}$ and 1.2, respectively. These parameters are iteratively solved using Equations (6) and (7), after which the shape parameter $\gamma$ is determined through Equation (5).

$$\mathsf{\gamma }=\mathsf{\delta }\ln \left({\displaystyle \frac{\mathsf{\lambda }}{d_{50\mathrm{\%}}}}-1\right)$$

(5)

$$\overline{d}=\mathsf{\epsilon }+\mathsf{\lambda }{\mathsf{\mu }}_1^{\prime }\left(y\right)$$

(6)

In Equation (7), where $G$ represents the basal area ($m^2\hspace{0.17em}h{a}^{-1}$), calculated by combining $\epsilon , \lambda$ and the first and second non-centered moments of the distribution of $Y$, ${\mathsf{\mu }}_1^{\prime }\left(y\right)$ and ${\mathsf{\mu }}_2^{\prime }\left(y\right)$. This calculation is weighted by $K$, which is determined by $\pi /40,000$, and incorporates tree density (in trees $h{a}^{-1})$.

$$G=KN\left[{\mathsf{\epsilon }}^2+2\mathsf{\epsilon }\mathsf{\lambda }{\mathsf{\mu }}_1^{\prime }\left(y\right)+{\mathsf{\lambda }}^2{\mathsf{\mu }}_2^{\prime }\left(y\right)\right]$$

(7)

The two-parameter Weibull, represented by Equation (8), was applied following the method by Bailey and Dell [37], where $b$ is the scale parameter and $c$ is the shape parameter. This method employed a moment-based two-parameter recovery approach according to Siipilehto and Mehtätalo [38], implemented in the R environment [32] with the lmfor package [39]. The standard input variables used for the recovery process include tree density $N$ (in trees per hectare, $\mathrm{h}{\mathrm{a}}^{-1}$) and basal area $G$ $({\mathrm{m}}^2\hspace{0.17em}\mathrm{h}{\mathrm{a}}^{-1}$) along with either the mean diameter $d_{mean}$ or the median diameter $d_{50\%}$.

$$f\left(x\right)=\left({\displaystyle \frac{c}{b}}\right){\left({\displaystyle \frac{x}{b}}\right)}^{c-1}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left({\displaystyle \frac{x}{b}}\right)}^c\right)$$

(8)

The equations below ensured that the applied stand mean or stand median diameter matched the corresponding characteristic derived from the recovered Weibull distribution [38]. The expected value of the Weibull distribution is given by Equations (9) and (10), where the scale parameter $b$ and the shape parameter $c$ define the distribution:

$$d_{mean}=b\mathsf{\Gamma }\left(1+{\displaystyle \frac{1}{c}}\right)$$

(9)

$$d_{50\%}=b{\left(\ln 2\right)}^{{\displaystyle \frac{1}{c}}}$$

(10)

The recovery is based on solving Equation (11), which is crucial to ensure that the distribution is consistent with the quadratic mean diameter of the stand.

$$dgW\left(\mathrm{c}, \mathrm{b} \left(D,\mathrm{c}\right)\right)-dg=0$$

(11)

where $dg$ is the quadratic mean diameter (cm), and $dgW$ represents the quadratic mean diameter of the Weibull distribution for a given parameter set, using the scale parameter that corresponds to the combination of the shape parameter and the mean or median diameter provided in $D$.

To compare the Johnson’s SB and Weibull Probability Density Functions adjusted using the PINASTER dataset, we applied the two-sample Kolmogorov–Smirnov (KS) test with the dgof R package [40]. This test assesses the D statistic by measuring the largest discrepancy between the observed and simulated distributions, providing a p-value to assess whether the null hypothesis (H0: observed distribution = simulated distributions) should be rejected. A lower KS statistic indicates a better fit, while a higher value suggests greater deviation. A p-value below 0.05 suggests rejection of (H0), indicating significant differences between the distributions, whereas a higher p-value suggests no reasons to reject H0.

Additionally, we compared the estimated growing stock volume from each PDF to the observed values in the PINASTER dataset plots. For this purpose, the diameter distribution estimated for each plot was combined with a height-diameter equation from Tomé et al. [41] and a volume equation from Tomé et al. [42] to estimate the volume for each plot measurement in the original dataset and the volume for the central value of each class of the simulated PDF. The total volume of each plot measurement was used to assess the accuracy by comparing the estimated and observed volumes.

To assess the accuracy of the volumes produced by each tested PDF, we used several statistics metrics: Mean Absolute Error (MAE), mean squared error (MSE), coefficient of determination (R-squared), root mean squared (RMSE%), and bias (%). These metrics were employed to select the most suitable PDF in predicting growing stock volume. Together, these statistics metrics provide a comprehensive assessment of precision, bias, and overall suitability of the tested PDFs in estimating observed volumes (Equations (12)–(16)).

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left|y_i-\widehat{y_i}\right|$$

(12)

$$\mathrm{MSE}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2$$

(13)

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(\widehat{Y_i}-Y_{ii}\right)}^2}{{\sum }_{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}},$$

(14)

$$\mathrm{RMSE} \%={\displaystyle \frac{100}{\overline{y}}}\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{n}}}$$

(15)

$$\mathrm{Bias} \%=100\sum _{i=1}^n{\displaystyle \frac{\left(y_i-\widehat{y_i}\right)}{n\overline{y}}}$$

(16)

where $y_i$is observed value, $\widehat{y_i}$ is the estimated value for the plot $i=1,...,n$, $\overline{y}$is the observed mean value, and n is the number of observations.

2.5. Estimating ABA-Derived Inputs to Fit PDF for Each Management Unit

This section outlines the process of obtaining stand variables through the Area-Based Approach (ABA), which is necessary for the parameter recovery of the selected Probability Density Function (PDFs). This PDF will then be used to simulate the diameter distribution across each Management Unit.

Following the analysis, which demonstrated that the Weibull PDF outperformed the Johnson’s SB PDF (see results section), we estimated the stand variables required for Weibull parameter recovery. These include the median diameter ${(\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}, \mathrm{c}\mathrm{m})$, quadratic mean diameter $(dg, \mathrm{c}\mathrm{m})$, stand density ($N, \mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{h}{\mathrm{a}}^{-1}$), and basal area $\left(G\right)$ $({\mathrm{m}}^2\hspace{0.17em}\mathrm{h}{\mathrm{a}}^{-1}$). In order to guarantee compatibility among estimated G, dg, and N, the basal area $\left(G\right)$ was obtained using the relation $G=dg * N$. Except for $\left(G\right)$, all inputs were estimated using the Area-Based Approach (ABA).

To estimate these variables, stand-level attributes from field inventory plots were linked to ALS metrics through linear regression models, hereafter referred to as “ABA models”. These models were subsequently applied to the grid cells corresponding to the independent variables selected for each model, allowing for a wall-to-wall estimation of each predicted stand variables across the entire MU area. Additionally, the dominant height—variable not needed as an input for PDF parameter recovery—was estimated through ABA to associate it with dominant trees. Finally, the pixel-level estimates within each MU were averaged to be used as the PDF parameter recovery inputs.

The ABA models consisted of a system of linear regression equations with four predictor variables each (Equation (17)).

$$Y={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{X}_1+{\mathsf{\beta }}_2{X}_2+{\mathsf{\beta }}_3{X}_3+{\mathsf{\beta }}_4{X}_4+\mathsf{\epsilon }$$

(17)

where $Y$ represents the dependent variable, ${\beta }_0$ is the intercept, ${\beta }_1, {\beta }_2,{\beta }_3,{\beta }_4$ are the coefficients of each independent variable $X_1,X_2,X_3,X_4$, and ε is the model error.

Several models were initially fitted for each stand variable, and the best models—using one, two, three, and four predictor variables—were selected via an exhaustive search algorithm from the leaps package [43]. This algorithm uses the Ordinary Least Squares approach to fit all possible models with a given number of predictors. Models were ranked based on: (i) lowest Akaike Information Criteria (AIC, Equation (18)), (ii) lowest Bayesian Information Criterion (BIC, Equation (19)), (iii) highest adjusted coefficient of determination (Equation (20)), and (iv) Variance Inflation Factor (VIF) under five. The VIF, calculated with the R car package [44], was used to prevent multicollinearity among the selected ALS metrics.

For model validation, we applied cross-validation using the Prediction Residuals Error Sum of Squares (PRESS, Equation (21)) [45]. Once the best metrics for each model were selected, we used Seemingly Unrelated Regression (SUR, [46]), with the systemfit R package [47] to account for cross-equation errors, ensuring consistent predictions for stand variables across each MU.

$$AIC=2k-2\mathit{ln}\left(L\right)$$

(18)

$$BIC=k ln\left(n\right)-2\mathit{ln}\left(L\right)$$

(19)

$$R_{\mathrm{adjusted}}^2=1-{\displaystyle \frac{\left(1-R^2\right)\cdot \left(n-1\right)}{n-p-1}}$$

(20)

$$\mathrm{PRESS}=\sum _{i=1}^n{\left(y_i-\widehat{y_{\left(i\right)}}\right)}^2$$

(21)

where $k$ is the number of estimated parameters in the model, $L$ is the value of the Likelihood function of the model, evaluated at the maximum likelihood, $y_i$ represents the observed value, and $\widehat{y_{\left(i\right)}}$ is the predicted value for i-th observation when the model is fitted without using the i-th individual observation.

2.6. PDF Parameter Recovery and Tree List Generation for Each Management Unit

This section outlines the Weibull PDF parameter recovery using stand variables obtained through the Area-Based Approach (ABA) in each Management Unit, followed by the generation of a tree list for each Management Unit. After confirming that the Weibull PDF outperformed the Johnson’s SB PDF (see results section), we proceeded with the two-parameter recovery for the Weibull PDF according to Siipilehto and Mehtätalo [38], using the lmfor R package [39] to obtain the shape and scale parameters for each Management Unit.

As mentioned in Section 2.5, all inputs for the two-parameter recovery of the Weibull PDF, except for basal area $\left(G\right)$, were estimated using the ABA. These include the median diameter ${(\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}})$, quadratic mean diameter $(dg, \mathrm{c}\mathrm{m})$, and the tree density ($N, \mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s} \mathrm{h}{\mathrm{a}}^{-1}$). These variables were used as input for the parameter recovery of the Weibull PDF.

Following the parameter recovery in each MU, a Monte Carlo Simulation (MCS) [48] was applied to generate a representative sample list of tree diameters, distributed according to the diameter classes in each MU. From this generated the tree list, dominant trees were identified based on the Burkhart and Tomé [2] definition of dominant height, which averages the height of the 100 thickest trees per hectare. The dominant height, predicted through ABA, was then assigned to these dominant trees, with a random variation added according to the regression standard error of prediction [49]. This methodological approach provides a detailed perspective on forest structure within each Management Unit, enhancing potential applications in forest management and decision-making.

3. Results

3.1. Probability Density Function Assessment

The KS test indicated that Weibull PDF was slightly more efficient in simulating diameter distribution for the PINASTER dataset compared to Johnson’s SB. Although the difference between the two approaches for parameter recovery of the Weibull PDF was relatively small, with similar results, Weibull ${(d}_{mean})$ was not rejected in 14 more plot measurements than Weibull ${(d}_{50\%})$. In contrast, Johnson’s SB was not rejected in 868 plots measurements (88.48%), which is 22 fewer plot measurements than the Weibull ${(d}_{mean})$ (

Table 3. Kolmogorov–Smirnov test (KS) applied to the PDF.

PDF	* KS (%)
Johnson’s SB	88.48 (868)
$\mathrm{Weibull} {(d}_{mean})$	90.72 (890)
$\mathrm{Weibull} {(\mathrm{d}}_{50\mathrm{\%}})$	89.29 (876)

* Kolmogorov–Smirnov test acceptance (%) and the number of plots not rejected by the test in parenthesis.

).

The KS test rejected 11.52% of the plot measurements for Johnson’s SB, 9.28% for Weibull ${(d}_{mean})$, and 10.71% for Weibull${(d}_{50\%})$. A total of 72 plots, mostly associated with very young stands where the maximum diameter was 5 cm, were rejected by the PDFs due to significant divergences. Most of these plots displayed bimodal patterns in the graphical analysis of the Johnson’s SB PDF. However, the other rejected plots showed no major divergences, but their p-values were still below the 0.05 significance level, leading to their rejection.

Regarding the second criteria, the growing stock prediction for Weibull $\left(d_{50\mathrm{\%}}\right)$ demonstrated higher accuracy, with a RMSE% of 13.88%, a lower bias of 7.82%, and high R-squared (0.974), indicating a better explanation of the variability in the estimated volumes. Additionally, the accuracy of Weibull${(d}_{mean})$ was comparable to that of Weibull ($d_{50\mathrm{\%}})$, with a RMSE% of 14.09% and a slightly higher bias of 7.90% (

Table 4. Accuracy assessment of the volume predicted with each PDF.

Statistics	Johnson’s SB	$\mathbf{Weibull} ({\mathit{d}}_{\mathit{m}\mathit{e}\mathit{a}\mathit{n}}$)	$\mathbf{Weibull} {(\mathit{d}}_{50\mathbf{\%}}$)
MAE	17.88	14.30	14.06
MSE	661.56	388.08	376.62
RMSE (%)	18.40	14.09	13.88
R²	0.95	0.97	0.97
Bias (%)	12.45	7.90	7.82

MAE: mean absolute error; MSE: mean squared error; RMSE: root mean squared error; R²: R-squared.

).

However, both Weibull approaches outperformed the Johnson’s SB PDF, which had the highest RMSE% (18.41%), the largest bias (12.45%), and lower R-squared values (0.954). These results underscore the superior predictive capability of the Weibull PDF when compared to Johnson’s SB in estimating growing stock volume for the PINASTER dataset.

3.2. Modeling Forest Attributes Inputs for the PDF Parameters Recovery in Each Management Unit

Weibull $\left(d_{50\mathrm{\%}}\right)$ was selected for parameter recovery in each Management Unit. This choice is justified by the importance of stock volume in the decision-making process for forest management practices, as it directly influences planning and management strategies. Additionally, Weibull $\left(d_{50\mathrm{\%}}\right)$ was also preferred for its robustness in representing the central tendency of diameter distribution, reducing the influence of extreme values that could skew the results.

After modeling and evaluating several models, we noticed that adding one, three, or four predictor variables did not significantly enhance performance. Therefore, we focused on models with two variables and applied the Seemingly Unrelated Regression to ensure consistent predictions for stand variables across each Management Unit. The best-performing system of equations and output statistics from the systemfit package are presented in

Table 5. Fitted equations with their respective accuracy assessment.

Variable 1	Independent Variable	Coefficients	R-Squared	Adjusted R-Squared	RMSE
${\mathrm{d}}_{50\mathrm{\%}}$	Intercept	15.45	0.45	0.41	6.72
	Zmean	1.64
	Pzbvzmn	−0.12
$\mathrm{d}\mathrm{g}$	Intercept	10.99	0.55	0.52	5.78
	Zq25	1.59
	ADD	1.94
N	Intercept	1940.91	0.63	0.60	232.10
	Zq85	−30.84
	Zpcum4	−18.99

¹ ${\mathrm{d}}_{50\mathrm{\%}}$: median diameter; dg: quadratic mean diameter; N: tree density (stems per hectare).

.

The predicted median diameter (${\mathrm{d}}_{50\mathrm{\%}})$ achieved an Adjusted R-squared of 0.41 and an RMSE of 6.72 cm, while the quadratic mean diameter ($\mathrm{d}\mathrm{g})$ model showed an Adjusted R-square of 0.52 with an RMSE of 5.78 cm. The tree density model demonstrated the highest performance, with an Adjusted R-square of 0.60 and an RMSE of 232.10 trees per hectare.

Based on these predictions, the second step of the Area-Based Approach methodology was employed to obtain the necessary inputs for estimating the parameters of the Weibull Probability Density Function for each Management Unit. The predictive equations, developed using stand variables and ALS metrics, were subsequently applied to the raster files corresponding to the selected metrics from the modelling process. These raster files, which cover the entire study area, enabled the generation of wall-to-wall estimates and maps of the predicted forest inventory attributes.

Finally, the mean pixel values within each Management Unit were calculated for each predicted variable, serving as inputs for the parameter recovery of the Weibull Probability Density Function. The results, which include the median diameter $\left(d_{50\mathrm{\%}}\right)$, quadratic mean diameter $\left(dg\right)$, and tree density $\left(N\right)$, are presented in Figure 3.

3.3. Probability Density Function Parameter Recovery and Tree List Generation for Each Management Unit

The two-parameter recovery of the Weibull${(d}_{50\mathrm{\%}})$ PDF was applied using the estimated stand variables for each Management Unit. Once the parameters were recovered, it was possible to graphically identify the shape and scale parameters of each simulated diameter distribution for the Management Units, thereby assessing the vertical structure of the forest in each of the 214 polygons (Figure 4).

Subsequently, a Monte Carlo Simulation (MCS) was conducted to generate a sample list of tree diameters from the simulated distribution within each Management Unit, based on the predicted tree density predicted by ABA.

Following this, we identified the dominant trees according to the definition of dominant height. The predicted dominant height, estimated using ALS metrics and field inventory data, was assigned to the dominant trees, with random variation applied based on the regression standard error [49]. This approach is particularly useful as it enables the calculation of tree height using specific hypsometric equation, although this falls outside the scope of the present study.

4. Discussion

This research conducted an innovative study by assessing the ability to simulate diameter distribution in 214 maritime pine Management Units in Portugal, using ALS and Probability Density Functions (PDFs). Among the highly flexible PDFs tested using the PINASTER dataset, the Weibull function demonstrated the best performance in terms of statistical accuracy (bias and precision) and fit to the observed data when compared to Johnson’s SB. Notably, both Weibull approaches performed similarly; as evidenced by the Kolmogorov test (KS), Weibull $(d_{mean}$) was not rejected in 90.72% of the plot measurements, and Weibull ${(d}_{50\mathrm{\%}})$ in 89.29%, both significantly exceeding the Johnson’s SB function, which was not rejected just in 88.48% of the plot measurements. This suggests that the Weibull PDF effectively captures the structural variation of diameters within the PINASTER dataset, particularly in homogeneous forests, aligning with the findings of Cosenza et al. [20], who also reported superior performance for the Weibull function over Johnson’s SB.

The high acceptance rate of the KS test across both Weibull PDFs evidence their suitability for simulating diameter distribution in Mediterranean forests in Portugal. However, the KS test rejected the same 72 plots in both PDFs, primarily associated with very young stands, with diameter no greater than 5 cm. Based on these results, we recommend against using this methodology for very young stands, as these plots did not show a simulated distribution that closely matched the observed distribution. This limitation suggests that adaptations or alternative approaches might be more suitable for young stands.

The analysis of prediction errors in growing stock volume estimated with each PDF revealed similar results for the two Weibull approaches. $Weibull {(d}_{50\mathrm{\%}})$ achieved an RMSE% of 13.88 and bias of 7.82%, while Weibull $(d_{mean}$) showed an RMSE% of 14.09 and bias of 7.90%. These values were considerably lower than those of Johnson’s SB, which had an RMSE of 18.40% and bias of 12.45%. These results confirm the effectiveness of the Weibull PDF not only in simulating diameter distributions but also in accurately predicting growing stock volume, which is crucial for forest management as it directly influences decision-making. Consequently, we selected the $Weibull {(d}_{50\mathrm{\%}})$ in conjunction with ABA to simulate diameter distribution within the Management Units.

The precision of forest inventory based on ALS data depends on tree density, which influences the choice between the Area-Based Approach (ABA) and Individual Tree Detection (ITC) [50]. We chose the ABA method because the maritime pine Management Units are dense and surrounded by eucalyptus MUs, adding complexity to capturing forest structure. This methodology, in conjunction with the ABA, proved effective for characterizing diameter distribution in the Management Units. The integration of ALS-derived metrics with predictive models enabled the estimation of forest stand variables required for parameter recovery of the Weibull to simulate diameter distributions, highlighting the practical applicability of this approach for characterizing Management Units.

Since Weibull ${(d}_{50\mathrm{\%}})$ was the PDF selected, we used stand variables predicted by the ABA to recover its shape and scale parameters in each Management Unit. Models developed to the predict median diameter ${(d}_{50\mathrm{\%}})$, quadratic mean diameter $\left(dg\right)$, and tree density $\left(N\right)$ of the Management Units (MUs) demonstrated Adjusted R-squared values of 0.41, 0.52, and 0.60, respectively. Although Seemingly Unrelated Regression (SUR) allowed consistent predictions across the different stand variables, a significant portion of the variability remains unexplained. The stand density model performed best, explaining 60% of the variation, though the RMSE of 232.10 trees per hectare highlights the potential for refinement.

These results suggest that models with two predictor variables based on ALS metrics capture part of the variability in the independent variables; however, adding further predictor variables did not increase the explanatory power, as we tested several models by adding one to four predictors. Further testing of additional predictor variables within the PDF parameter prediction method, where parameters are estimated through regressions using stand variables [38], may enhance the model’s ability to capture forest structure in complex stands, rather than using parameter recovery method. Furthermore, the study area, severely impacted by intense wildfires in recent years [51], may have impacted the field data collection carried out in 2022 by reducing the number of trees measured and limiting data collection.

The application of two-parameter recovery on the Weibull DF, using the modeled variables for each MU, enabled the identification of different shapes and scales of the simulated diameter distributions in each Management Unit. For individual tree models, which are required when the objective is to simulate complex forests, namely close-to-nature management, the Monte Carlo Simulation is crucial to obtain a sample list of tree diameter in each MU. Furthermore, the identification of dominant trees, and the correspondent attribution of dominant height to the trees classified as dominant, adding random variation [49], enabled the assignment of predicted heights from ALS metrics and stand variables to each tree via a height-diameter curve.

Beyond the Johnson’s SB and Weibull distributions, other Probability Density Functions such as Beta, Gamma, Lognormal, and Truncated Weibull could be tested for various species to simulate diameter distribution. These alternatives in different contexts may complement Weibull strengths, particularly in MUs with greater structural complexity or heterogeneity. Nevertheless, the focus was placed on Johnson’s SB and the Weibull PDF due to their flexibility in simulating diameter distributions in each MU. Future research could explore the potential of these additional PDFs for applications in varying forest conditions.

Our methodology was essential in providing detailed insights into forest structure and stand-level characterization, both of which are fundamental for simulating stand dynamics based on growth and yield simulators. Consequently, the output from this work enables the use of individual tree models to simulate thinning, conversion from even-aged into uneven-aged stands, or conversion from pure to mixed stands. Characterizing Management Units (MUs) with tree-level data is instrumental for a comprehensive understanding of forest structure and supports informed forest management decision-making.

5. Conclusions

This work assessed the ability of the highly flexible Weibull and Johnson’s SB Probability Density Functions to simulate diameter distribution using the maritime pine PINASTER database from Portugal, which covers 12 trials each with several plots and treatments. The results demonstrated that the Weibull PDF was more effective than Johnson’s SB in simulating diameter distributions for this species. Additionally, this research provided novel insights into the effectiveness of the Weibull PDF in simulating diameter distribution in maritime pine Management Units using ALS data and the Area-Based Approach method.

This study makes a significant contribution to the existing literature by offering an innovative approach to the characterization of forest Management Units at individual tree levels. The developed methodology can be applied across different forest Management Units, providing a practical and efficient solution for estimating tree-level data, which is crucial for the application of individual tree models.

In conclusion, this study highlights the importance of exploring and developing methodologies based on remote sensing technologies, such as ALS, to enhance forest management and conservation. Continued research in this area could open new possibilities for forest management in different ecosystems, contributing to the sustainability and resilience of forested environments.