Modeling Diameter Distributions with Six Probability Density Functions in Pinus halepensis Mill. Plantations Using Low-Density Airborne Laser Scanning Data in Aragón (Northeast Spain)

Abstract

The diameter distributions of trees in 50 temporary sample plots (TSPs) established in Pinus halepensis Mill. stands were recovered from LiDAR metrics by using six probability density functions (PDFs): the Weibull (2P and 3P), Johnson’s SB, beta, generalized beta and gamma-2P functions. The parameters were recovered from the first and the second moments of the distributions (mean and variance, respectively) by using parameter recovery models (PRM). Linear models were used to predict both moments from LiDAR data. In recovering the functions, the location parameters of the distributions were predetermined as the minimum diameter inventoried, and scale parameters were established as the maximum diameters predicted from LiDAR metrics. The Kolmogorov–Smirnov (KS) statistic (D_n), number of acceptances by the KS test, the Cramér von Misses (W²) statistic, bias and mean square error (MSE) were used to evaluate the goodness of fits. The fits for the six recovered functions were compared with the fits to all measured data from 58 TSPs (LiDAR metrics could only be extracted from 50 of the plots). In the fitting phase, the location parameters were fixed at a suitable value determined according to the forestry literature (0.75·d_min). The linear models used to recover the two moments of the distributions and the maximum diameters determined from LiDAR data were accurate, with R² values of 0.750, 0.724 and 0.873 for d_g, d_med and d_max. Reasonable results were obtained with all six recovered functions. The goodness-of-fit statistics indicated that the beta function was the most accurate, followed by the generalized beta function. The Weibull-3P function provided the poorest fits and the Weibull-2P and Johnson’s SB also yielded poor fits to the data.

1. Introduction

LiDAR (Laser Imaging Detection and Ranging) technology allows large areas to be scanned, thus generating continuous information about the entire space. In the last 20 years, LiDAR has been increasingly used in forest inventories at different scales [1,2], because of its ability to provide detailed three-dimensional information on the size and structure of forest cover [3]. The process provides data on the dimensions of the trees and on the structure of the forest cover, and other parameters such as canopy fuel attributes [4]. The variables most closely related to LiDAR data are canopy cover and tree height, which mainly depend on the vertical distribution of the layers forming the tree cover [5,6].

The main advantage of using field data that will then be processed using LiDAR techniques is that a much smaller sampling effort is required (with a good sampling design), regarding both the number of plots and the number of tree height measurements required. LiDAR technology, such as airborne laser scanning (ALS), provides a great deal of information about vegetation height and is usually quite accurate, especially in coniferous plantations. This technology therefore represents the future of forest inventories, especially on a large scale, greatly reducing the costs associated with inventories, which are traditionally the most expensive part of forest management projects.

One of the most common approaches used in ALS forest inventories is the area-based approach (ABA), in which metrics extracted from the normalized height of the LiDAR data cloud (NHD) are used to predict forest variables [7]. Although the ABA does not directly detect tree diameters, it determines the forest stand structure indirectly by using the NHD metrics to estimate probability density functions (PDF) that describe diameter distributions [6].

The empirical diameter distribution (specified by the DBH measurements within the stand) described by a cumulative distribution function (CDF) or a PDF enables prediction of stand growth and planning of various uses and provides data about stand structure. The diameter distribution of the total number of trees in a stand is important because the model parameters and stand variables provided in yield tables are closely related [8]. However, DBH measurements are not always available, and therefore the diameter distribution must be predicted by using stand attributes as explanatory variables, usually under the assumption that they can be described by a specified theoretical PDF [6,9].

However, the parameters of the PDFs can also be estimated from LiDAR metrics. Three approaches have been considered for estimating the diameter distribution from LiDAR data, within the framework of parametric prediction [6]: (a) use of regression analysis and parameter recovery models (PRM) to relate LiDAR metrics directly to PDF parameters or the moments or percentiles of the diameter distribution, (b) modeling the PDF parameters from stand-level variables predicted using area-based LiDAR metrics and (c) predicting the diameter distribution on the basis of recognition of individual trees, which requires high pulse densities. The first approach has recently been used, with good results [10,11], and is the methodology used in the present study.

Several CDFs or PDFs have been used to describe and predict tree diameter frequency in forest stands, including Gram-Charlier [12], normal [13], Birnbaum–Saunders, and gamma [14], beta [15], generalized beta [16], Johnson’s SB [17] and the Weibull [18] functions. To date, the parameters of these distributions have only been recovered from LiDAR metrics by applying PRM to the Weibull distribution [10,19,20] and Johnson’s SB distribution [11]. However, this approach has never been used with distributions such as the beta, generalized beta and gamma-2P functions.

Thus, the main objective of this study was to predict and to compare the diameter distributions in Pinus halepensis Mill. plantations in Aragón (NE Spain) recovered by PRM in combination with the Weibull (2P and 3P), Johnson’s SB, beta, generalized beta and gamma-2P functions, and using height and canopy cover LiDAR metrics from low-density airborne scanning data provided by the Spanish country-wide PNOA (Plan Nacional de Ortografía Aérea), which covers the whole country, so our study can be scalable to other species using the same methodology and type of data.

2. Materials and Methods

2.1. Dataset

The data used in this study were obtained from plantations of Aleppo pine (Pinus halepensis Mill.) in the region of Aragón (NE Spain). P. halepensis is distributed naturally throughout the Mediterranean area of the southeast of the Iberian Peninsula and the Balearic Islands [21]. Plantations of the species represent 15% of those established between 1940 and 1980 in Spain. However, the species is not particularly productive relative to other pines. Scant attention has therefore been given to Aleppo pine, even though at local scales it may be the only economic alternative available and may complement other non-timber harvests [22].

A total of 58 temporary sampling plots (TSPs) of P. halepensis Mill., measured in October 2017, were used for the present study, and LiDAR metrics were able to be extracted from the data corresponding to 50 of these plots (Figure 1). The plot size ranged from 225 to 625 m², to achieve a minimum number of 30 trees per plot. Diameter at breast height (DBH at 1.3 m above the ground) of a total of 1685 trees was measured with callipers, to an accuracy of 0.1 cm. The following field and LiDAR variables were calculated from the inventory data: quadratic mean diameter, mean diameter, maximum diameter, number of trees per hectare, dominant height, basal area and total LiDAR return density within the plots (pulses·m⁻²). Summary statistics are shown in

Table 1. Summary of field and LiDAR data for the 50 plots for which LiDAR metrics were available.

Species	Variable	Mean	Max	Min	SD
Pinus halepensis	dg	17.8	30.6	11.7	5.1
	dmed	17.3	29.6	11.4	4.8
	dmax	25.7	49.6	18.8	7.4
	N	1054	3200	176	588.4
	Ho	10.1	19.1	6.2	3.2
	G	24.1	58.9	6.3	11.9
	LRD	1.070	2.222	0.453	0.421

d_g: quadratic mean diameter; d_med: mean diameter; d_max: maximum diameter; N: number of trees per hectare; Ho: dominant height; G: basal area; LRD: total LiDAR return density within the plots (pulses·m⁻²).

.

2.2. Lidar Metrics

The LiDAR data covering 50 plots in Pinus halepensis stands were acquired in a national survey carried out for the PNOA (Plan Nacional de Ortofotografía de España) project between 16 October and 16 November 2016, under the direction of the Spanish Ministerio de Fomento (Dirección General del Instituto Geográfico Nacional and Centro Nacional de Información Geográfica). The data were obtained with a LEICA ALS80 sensor, at a pulse repetition rate of 45 kHz, scan frequency of 70 Hz and a maximum scan angle of ±50°. A maximum of 4 returns per pulse were registered, reaching an average point density of 1.070 points·m⁻², with a theoretical laser pulse density required for the PNOA project of 0.5 first returns per square meter. Summary statistics of the LiDAR return density per square meter within the plots are shown in Table 1. The LiDAR data were also processed with FUSION software [23]. As reported by the provider (https://pnoa.ign.es/) accessed on 1 April 2020, the vertical accuracy of the LiDAR metrics, given by the RMSE, is ≤0.20 m. The set of metrics from the points laid above 2 m was extracted for each plot (

Table 2. Potential explanatory variables related to height distribution and canopy closure.

Variables Related to Height Distribution (m)	Description
LH_MIN, LH_MAX, LH_MEAN	Minimum, maximum and mean height
LH_MODE, LH_MEDIAN, LH_SD, LH_CV	Mode, median, standard deviation and height’s coefficient of variation
LH_SK, LH_KUR	Skewness and kurtosis
LH_IQ	Interquartile amplitude
LH_AAD	Mean absolute deviation
LH_MAD_MEDIAN, LH_MAD_MODE	Median of the absolute deviations from the overall height median (LH_MAD_MEDIAN) and mode (LH_MAD_MODE)
LH_L1, LH_L2…, LH_L4	L moments
INT_L_SK, INT_L_KUR	Linear combinations of L moments (skewness and kurtosis)
LH_P05,…, LH_P95	Percentiles
LH_P25; LH_P75	First and third quartiles
Variables Related to Canopy Closure (%)	Description
LFCC	Percentage of first returns above 2 m
LFCC_MEAN	Percentage of first returns above LH_MEAN
LFCC_MODE	Percentage of first returns above LH_MODE
LFCC_ALL	Percentage of all returns above 2 m
LFCC_ALL_MEAN	Percentage of all returns above LH_MEAN
LFCC_ALL_MODE	Percentage of all returns above LH_MODE
ALL_MEAN_FIRST	100* all returns above LH_MEAN / total first returns
ALL_FIRST	100* all returns above 2 m / total first returns
R2_COUNT	Number of first returns above 2 m
CANOPY	Canopy relief ratio: (hmean − hmin)/(hmax − hmin)

).

2.3. Diameter Distribution Models and Fitting

2.3.1. The Weibull Function

The three-parameter Weibull PDF has the following expression for a continuous random variable x [24]:

$$f\left(x\right)=\left(\frac{c}{b}\right)·{\left(\frac{x-a}{b}\right)}^{c-1}·e^{-{\left(\frac{x-a}{b}\right)}^c}$$

(1)

where f(x) is the probability density of trees with diameter x, a represents the location, b the scale and c the shape. If a equals zero, the expression corresponds to the two-parameter Weibull function.

The method of moments was used to fit all six distributions on the basis of the relationship between the parameters and the first and second moments of the diameter distribution (arithmetic mean diameter and variance). For the Weibull distribution, it is computed by the following expressions [10,18]:

$$b=\frac{\overline{d}-a}{\mathsf{\Gamma }\left(1+\frac{1}{c}\right)}$$

(2)

$${\sigma }^2=\frac{{\left(\overline{d}-a\right)}^2}{{\mathsf{\Gamma }}^2\left(1+\frac{1}{c}\right)}\left[\mathsf{\Gamma }\left(1+\frac{2}{c}\right)-{\mathsf{\Gamma }}^2\left(1+\frac{1}{c}\right)\right]$$

(3)

where $\overline{d}$ is the arithmetic mean diameter of the distribution, ${\sigma }^2$ is the variance and $\mathsf{\Gamma }(·)$ is the gamma function. Equation (3) was resolved by a bisection iterative procedure [25].

2.3.2. The Beta Function

The general expression of the beta distribution for a random variable x is given by [26]:

$$f\left(x\right)=c·{\left(x-L\right)}^{\propto }·{\left(U-x\right)}^{\gamma }$$

(4)

$$withc=\frac{1}{\frac{{\left(-L+U\right)}^{1+\gamma }·\mathsf{\Gamma }\left(1+\alpha \right)·\mathsf{\Gamma }\left(1+\gamma \right)}{{\left(\frac{1}{-L+U}\right)}^{\alpha }·\mathsf{\Gamma }\left(2+\alpha +\gamma \right)}}$$

(5)

for the interval L ≤ x ≤ U, and x = 0 otherwise, where x is the diameter at breast height and is assumed to be continuous, f(x) is the density associated with diameter x, U and L are respectively the upper and lower limits of the beta distribution, c is the scaling factor, α and γ are respectively the first and the second exponents that determine the shape of the distribution and $\mathsf{\Gamma }(·)$ is the gamma function.

The method of moments for the beta function [9,15,27] is computed by the following equations:

$$\gamma =\frac{\frac{Z}{s_{rel}^2·{\left(Z+1\right)}^2}-1}{Z+1}-1$$

(6)

$$\alpha =Z·\left(\gamma +1\right)-1$$

(7)

$$Z=\frac{x_{rel}}{1-x_{rel}}$$

(8)

$$x_{rel}=\frac{\overline{d}-L}{U-L}$$

(9)

$$s_{rel}^2=\frac{s^2}{{\left(U-L\right)}^2}.$$

(10)

where $\overline{d}$ is the arithmetic mean diameter of the distribution and $s^2$ is the variance.

2.3.3. The Generalized Beta Distribution (GBD)

The general expression of the generalized beta PDF for a random variable x is as follows [16]:

$$f\left(x\right)=\complement ·{\beta }_2^{-\left({\beta }_3+{\beta }_4+1\right)}·{\left(x-{\beta }_1\right)}^{{\beta }_3}·{\left({\beta }_1+{\beta }_2-x\right)}^{{\beta }_4}$$

(11)

$$where\hspace{1em}\complement =\frac{\mathsf{\Gamma }\left({\beta }_3+{\beta }_4+2\right)}{\mathsf{\Gamma }\left({\beta }_3+1\right)·\mathsf{\Gamma }\left({\beta }_4+1\right)}$$

(12)

for the interval $\left({\beta }_1,{\beta }_1+{\beta }_2\right),$ and 0 otherwise, where x is the diameter at breast height and is assumed to be continuous, f(x) is the density associated with diameter x, ${\beta }_1$ and ${\beta }_2$ are respectively the lower and upper limits of the distribution, ${\beta }_3$ and ${\beta }_4$ are exponents that determine the shape of the distribution, $\complement$ is the scaling factor of the function and $\mathsf{\Gamma }(·)$ is the gamma function.

The method of moments for the GBD function is computed by the following expressions [16,28]:

$${\beta }_1=x_{min}\hspace{1em}{\beta }_3=\frac{{\left({\alpha }_1-{\beta }_1\right)}^2\left({\beta }_1+{\beta }_2-{\alpha }_1\right)-{\alpha }_2\left({\beta }_2-{\beta }_1+{\alpha }_1\right)}{{\alpha }_2{\beta }_2}$$

(13)

$${\beta }_2=x_{max}-{\beta }_1\hspace{1em}{\beta }_4=\frac{{\beta }_2\left({\beta }_3+1\right)}{{\alpha }_1-{\beta }_1}-{\beta }_3-2$$

(14)

where x is the tree diameter, ${\mathsf{\alpha }}_1$ is the arithmetic mean diameter and ${\alpha }_2$ is the variance.

2.3.4. The Johnson’s SB Function

The model of the Johnson’s SB PDF has the following expression for a continuous random variable x [29]:

$$f\left(x\right)=\frac{\delta }{\sqrt{2\pi }}·\frac{\lambda }{\left(\epsilon +\lambda -x\right)\left(x-\epsilon \right)}·e^{-\frac{1}{2}{\left[\gamma +\delta ·ln\left(\frac{x-\epsilon }{\epsilon +\lambda -x}\right)\right]}^2}$$

(15)

where f(x) is the probability density associated with diameter x, ε < x < ε + λ, −∞ < ε < ∞, −∞ < γ < ∞, λ > 0, and δ > 0.

The model is characterized by the parameters ε (location), λ (scale), γ (asymmetry) and δ (kurtosis). The method of moments used in [11,30] is computed as follows:

$$\delta \cong \frac{\mu \left(1-\mu \right)}{Sd\left(x\right)}+\frac{Sd\left(x\right)}{4}\left[\frac{1}{\mu \left(1-\mu \right)}-8\right]$$

(16)

$$\gamma \cong \delta ln\left(\frac{1-\mu }{\mu }\right)+\left(\frac{0.5-\mu }{\delta }\right)$$

(17)

$$\mu =\frac{\overline{d}-\epsilon }{\lambda }$$

(18)

$$Sd\left(x\right)=\frac{{\sigma }_x}{\lambda }$$

(19)

where $\overline{d}$ is the arithmetic mean of the plot diameters, Sd(x) is the modified standard (δ) deviation and ${\sigma }_x$ is the plot diameter standard deviation.

2.3.5. The Gamma Function

The model of the gamma PDF has the following expression for a continuous random variable x [14,31]:

$$f\left(x\right)=\frac{{\left(x-\gamma \right)}^{\alpha -1}}{{\beta }^{\alpha }\mathsf{\Gamma }\left(\alpha \right)}·e^{\left[-\left(\frac{x-\gamma }{\beta }\right)\right]}$$

(20)

with x > $\gamma$, $\alpha$ > 0 and $\beta$ > 0, where $\alpha$ is the shape parameter, $\beta$ is an inverse scale parameter, $\gamma$ is the location parameter ($\gamma =0$ for the two-parameter gamma distribution) and $\mathsf{\Gamma }(·)$ is the gamma function.

The method of moments for the gamma-2P function is computed by the following expressions [28]:

$$\alpha ={\left(\frac{\overline{d}}{s}\right)}^2$$

(21)

$$\beta =\frac{s^2}{\overline{d}}$$

(22)

where $\overline{d}$ is the mean diameter and s is the standard deviation.

For fitting the distributions to the observed data (phase 1), a value of 0.75·d_min was used for the location parameters, as [32] achieved good results with this value compared to other constraints. For the Weibull-2P and gamma-2P functions, the value of the location parameters was zero. Scale parameters and the upper limit of the four-parameter functions (Johnson’s SB, beta and generalized beta) were established as the maximum diameter of the distributions (d_max) and the upper limit of the largest diameter class for the beta function [27,28] to improve parameters’ convergence. The size of the diameter classes was assumed to equal 1 cm.

2.4. Goodness of Fit Evaluation

The consistency of the functions was evaluated using the Kolmogorov–Smirnov (KS) (D_n) and Cramér von Mises (W²) statistics, bias, mean square error (MSE) and number of acceptances by the KS test. For a given cumulative distribution function, F(x), $D_n=su{p}_x|F_n\left(x\right)-F_0\left(x\right)|$, where sup_x is the supremum of the set of distances, calculated as follows [33]:

$$D_n=max\left\{ma{x}_{1\le i\le n_i}\left[F_n\left(x_i\right)-F_0\left(x_j\right)\right],ma{x}_{1\le i\le n_i}\left[F_0\left(x_j\right)-F_n\left(x_{i-1}\right)\right]\right\}$$

(23)

where the cumulative observed and estimated frequencies, i.e., F_n(x_i) and F₀(x_j), are compared.

The recovered distributions were also assessed using the two-sample Kolmogorov–Smirnov test (KS) under the null hypothesis that the plot data originate from the recovered distribution. However, as the parameters were estimated empirically, the theoretical distribution for each plot is unknown, and the KS test should therefore be conducted using a Monte Carlo simulation [34]. The procedure is explained in further detail in [10,11].

The Cramér von Mises statistic (W²) is a measure of the square of the distance between the empirical and the cumulative theoretical distribution [35]:

$$W^2={\displaystyle \sum }_{i=1}^n{\left\{\widehat{F}\left(x_i\right)-\frac{\left(i-0.5\right)}{n}\right\}}^2+\frac{1}{12n}$$

(24)

where $\widehat{F}$($x_i)$ is the cumulative theoretical distribution in diameter class i, and n is the number of diameter classes.

The bias and the mean square error (MSE) were also used as goodness-of-fit measures and were expressed as follows:

$$Bias=\frac{{{\displaystyle \sum }}_{i=1}^N{Y}_i-{\widehat{Y}}_i}{N}$$

(25)

$$MSE=\frac{{{\displaystyle \sum }}_{i=1}^N{\left(Y_i-{\widehat{Y}}_i\right)}^2}{N}$$

(26)

where Y_i is the observed relative frequency of trees in each diameter class, $\widehat{Y_i}$ is the theoretical value predicted by the model and N is the number of diameter classes. The bias and MSE were calculated for each fit as the mean relative frequency of trees.

2.5. Recovering the Parameters of the Distributions from LiDAR Metrics

The function parameters were recovered from the first two moments of the distributions (mean diameter and variance) by using parameter recovery models (PRMs). The variance $\left({\sigma }_d^2\right)$ is related to the quadratic mean diameter (d_g) and the mean diameter (d_m), as follows:

$${\sigma }_d^2=d_g^2-d_m^2$$

(27)

Within the LiDAR data framework, the stand level attributes usually used for the recovery procedure (e.g., t, N, H, S) can be replaced by LiDAR metrics as explanatory variables to estimate d_g and d_m [10]. However, in the four-parameter functions (Johnson’s SB, beta and generalized beta), the range and the upper limit of the distributions can also be estimated. In this case, both values were used as the maximum diameter (d_max), which was then related to LiDAR metrics. Location parameters of the recovered functions were predetermined as the minimum diameter inventory (7.5 cm) for the Weibull-3P, Johnson’s SB, beta and generalized beta functions, and as zero for the Weibull-2P and gamma-2P functions.

We used linear models to establish the empirical relationships between d_max, d_g and LiDAR metrics:

$$d_{max}={\alpha }_0+{\alpha }_1{X}_1+{\alpha }_2{X}_2+\dots +{\alpha }_m{X}_m+\epsilon$$

(28)

$$d_g={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\dots +{\beta }_m{X}_m+\epsilon$$

(29)

where d_max and d_g are the dependent variables, X₁, X₂, …, X_m represents the independent, potentially explanatory variables related to the LiDAR-derived height distribution and canopy closure, ${\alpha }_0$, ${\alpha }_1$, …${\alpha }_n$ and ${\beta }_0$, ${\beta }_1$, …${\beta }_n$ are the parameters to be estimated in the fitting process and $\epsilon$ is the additive error term, which is assumed to be independent and normally and identically distributed, with zero mean.

For a given stand, d_m is always smaller than or equal to d_g, and we therefore used the following model expression to take this restriction into account [9,10,11,36]:

$$d_m=d_g-\exp \left({\delta }_0+{\delta }_1{X}_1+{\delta }_2{X}_2+\dots +{\delta }_m{X}_m\right)+\epsilon$$

(30)

where X₁, X₂, …, X_m are the potential explanatory variables related to the LiDAR-derived height distribution and canopy closure, ${\delta }_0$,${\delta }_1$,…${\delta }_m$ are the parameters to be estimated in the fitting process and $\epsilon$ is the additive error term.

Equation (30) was linearized by applying a natural logarithmic transformation to facilitate selection of the independent variables. Finally, Equation (29) (when d_g was the dependent variable) and Equation (30) were fitted simultaneously with “seemingly unrelated regression” (SUR) to prevent cross-correlation between error components. Goodness of fits were evaluated with the coefficient of determination (R²) and the root mean square error (RMSE).

3. Results

The mean, standard deviation, minimum and maximum values of the parameters estimated in step 1 (fitting the distributions to observed data) and step 2 (recovering the parameters of the distributions from LiDAR metrics) are shown in

Table 3. Descriptive statistics for the parameters of the functions in the fitting and recovery steps.

Function	Step	Param	Mean	SD	Min	Max
Weibull-2P	Fitting	b	18.757	5.072	12.415	32.084
		c	4.905	1.125	2.163	8.020
	Recovery	b	18.918	4.614	13.699	33.740
		c	4.873	0.825	2.215	6.577
Weibull-3P	Fitting	a	7.184	1.750	5.625	14.100
		b	11.222	4.005	6.312	24.255
		c	2.598	0.460	1.535	3.560
	Recovery	a	7.500	-	7.500	7.500
		b	11.008	4.637	5.513	25.319
		c	2.489	0.619	1.245	3.646
beta	Fitting	c	0.006	0.015	1.47 × 10–6	0.107
		L	7.184	1.750	5.625	14.100
		U	25.739	7.403	16.700	49.600
		α	1.191	0.633	0.250	2.574
		γ	1.015	0.720	0.166	3.178
	Recovery	c	0.010	0.015	4.34 × 10–7	0.069
		L	7.500	-	7.500	7.500
		U	26.100	7.229	18.000	55.000
		α	0.972	0.690	0.017	3.023
		γ	0.734	0.488	0.046	2.003
Generalized beta	Fitting	$\complement$	471.180	925.744	1.251	5827.044
		B1	7.184	1.750	5.625	14.100
		B2	25.739	7.403	16.700	49.600
		B3	2.253	1.034	0.543	4.885
		B4	4.285	1.806	0.358	9.363
	Recovery	$\complement$	365.657	868.723	1.428	5125.128
		B1	7.500	-	7.500	7.500
		B2	26.104	7.222	18.286	55.318
		B3	2.057	1.100	0.435	4.647
		B4	4.116	1.015	0.974	6.354
Johnson’s SB	Fitting	ε	7.184	1.750	5.625	14.100
		λ	25.739	7.403	16.700	49.600
		γ	0.727	0.344	0.027	1.379
		δ	1.327	0.235	0.630	1.879
	Recovery	ε	7.500	-	7.500	7.500
		λ	26.104	7.222	18.286	55.318
		γ	0.801	0.409	−0.152	1.457
		δ	1.295	0.190	0.698	1.724
gamma-2P	Fitting	α	19.906	8.426	4.636	48.331
		β	1.062	0.815	0.382	5.329
	Recovery	α	18.594	5.410	4.398	31.558
		β	1.089	0.802	0.562	5.093

Param: parameter; SD: standard deviation; Min: minimum; Max: maximum.

.

The mean values of the parameters were reasonable in both cases, with no large differences in the values produced by each method.

Table 4. Mean values of statistics used for comparison of the functions in the fitting step.

Function	Dn	W2	Bias	MSE
Weibull-2P	0.146922	0.081863	0.003424	0.002009
Weibull-3P	0.150761	0.046157	0.002637	0.001904
beta	0.138560	0.062044	0.001643	0.001851
Generalized beta	0.153656	0.048215	0.002563	0.001879
Johnson’s SB	0.162578	0.043088	0.002276	0.001892
gamma-2P	0.164750	0.054360	0.002736	0.001996

D_n: Kolmogorov–Smirnov statistic; W²: Cramér von Mises statistic; MSE: mean squared error.

shows the mean value of the statistics used to test the goodness-of-fits to the observed distributions: Kolmogorov–Smirnov (D_n), Cramér von Mises (W²), bias and mean squared error (MSE).

In the fitting phase, the lowest value of D_n was yielded by the beta function (0.138560), followed by the Weibull-2P function (0.146922) and the Weibull-3P function (0.150761). The highest value was yielded by the gamma-2P function (0.164750), followed by the Johnson’s SB function (0.162578) and the generalized beta function (0.153656). Regarding the Cramér von Mises statistic (W²), the lowest value corresponded to the Johnson’s SB function (0.043088) and the highest to the Weibull-2P function (0.081863). The order of the functions in terms of W² was Johnson’s SB < Weibull-3P < generalized beta < gamma-2P < beta < Weibull-2P. The smallest MSE value corresponded to the beta function (0.001851), followed by the generalized beta function (0.001879). The highest MSE value was yielded by the Weibull-2P function (0.002009). The order for the six distributions studied was beta < generalized beta < Johnson’s SB < Weibull-3P < gamma-2P < Weibull-2P. The bias may be less important for comparison of the results because errors with different signs can compensate each other in the total mean value. Thus, considering the values of D_n, W² and MSE for the six distributions studied, the beta function provided the most accurate fits to the observed data and the Weibull-2P and the Gamma-2P functions yielded the poorest fits.

The parameter estimates and goodness-of-fit statistics (R², RMSE and % RMSE) of the simultaneous fitting of Equations (29) and (30) used to estimate the mean diameter (d_m) and the quadratic mean diameter (d_g) from LiDAR data and for Equation (28) used to estimate d_max are shown in

Table 5. Parameter estimates and goodness-of-fit statistics for the simultaneous fitting of Equations (29) and (30) used to estimate d_med and d_g from LiDAR data and for Equation (28) used to estimate d_max.

Equation	Dep var	Independent Variable	Param	Param Estim	$\mathbf{P}>|\mathit{t}|$	R2	RMSE	RMSE%
(28)	dmax	Intercept	${\alpha }_0$	6.733	<0.0001	0.873	2.758	10.56
		LH_AAD	${\alpha }_1$	8.143	<0.0001
		LH_P95	${\alpha }_2$	7.728	<0.0001
(29)	dg	Intercept	${\beta }_0$	5.329	<0.0001	0.750	2.542	14.24
		LH_P90	${\beta }_1$	1.335	<0.0001
(30)	dmed	Intercept	${\delta }_0$	−1.789	<0.0001	0.724	2.549	14.72
		LH_MAD_MEDIAN	${\delta }_1$	1.049	<0.0001

LH_AAD: Mean absolute deviation for height; LH_P95: 95% height percentile; LH_P90: 90% height percentile; LH_MAD_MEDIAN: Median of the absolute deviations from the overall median height.

.

The maximum diameter of the distributions (d_max) was considered the scale parameter for the Weibull-3P, Johnson’s SB and generalized beta functions and for establishing the upper limit of the beta function. The maximum diameter was predicted by the following independent variables: mean absolute deviation for height (LH_AAD) and 95% height percentile (LH_P95). The linear model yielded an accuracy of R² = 0.873 and RMSE = 2.758. Models [29,30] used for simultaneous fitting of quadratic mean diameter (d_g) and the mean diameter (d_med) yielded R² values of 0.750 and 0.724 and RMSE values of 2.542 and 2.549, respectively. Good results were obtained for these key variables in terms of recovering the diameter distributions with the six functions.

For recovering the distributions from LiDAR data (

Table 6. Mean values of statistics used for comparison of functions in the recovery step.

Function	Dn	W2	Bias	MSE	KS Acceptance (%)
Weibull-2P	0.282500	0.498277	0.004303	0.002732	28 (56%)
Weibull-3P	0.286170	0.497091	0.003619	0.002746	26 (52%)
beta	0.254078	0.381916	0.003744	0.002649	35 (70%)
Generalized beta	0.264840	0.389917	0.003673	0.002711	33 (66%)
Johnson’s SB	0.273945	0.411657	0.003578	0.002851	32 (64%)
Gamma-2P	0.272107	0.398073	0.004047	0.002792	34 (68%)

D_n: Kolmogorov–Smirnov statistic; W²: Cramer von Mises statistic; MSE: mean squared error.

), the smallest value of D_n was obtained by the beta function (0.254078), as in the fitting step, and then by the generalized beta function (0.264840) and the gamma-2P function (0.272107). The highest value was yielded by the Weibull-3P function (0.286170), followed by the Weibull-2P function (0.282500) and the Johnson’s SB function (0.273945).

The lowest value of the Cramér von Mises statistic (W²) was also yielded by the beta function (0.381916) and the highest by the Weibull-2P function (0.498277). The order of the distributions for the W² value was beta < generalized beta < gamma-2P < Johnson’s SB < Weibull-3P < Weibull-2P. As for D_n and W², the smallest value for the MSE corresponded to the beta distribution (0.002649), followed by the generalized beta function (0.002711). The highest value of the MSE corresponded to the gamma-2P function (0.002792). The order in terms of MSE for the six distributions recovered was beta < generalized beta < Weibull-2P < Weibull-3P < gamma-2P < Johnson’s SB. As for the other statistics, bias was also higher for recovery of the distributions from LiDAR metrics than for fitting to the observed data. Results for the KS test are also consistent with the other statistics, with the following order for plots accepted: beta function (35 plots: 70%), gamma-2P function (34 plots: 68%), generalized beta function (33 plots: 66%), Johnson’s SB function (32 plots: 64%), Weibull-2P function (28 plots: 56%) and Weibull-3P function (26 plots: 52%). Thus, considering the values of D_n, W², MSE and KS test for the six distributions studied, the beta function was the most accurate for recovering the parameters from LiDAR data, followed by the generalized beta function. The Weibull-3P, Weibull-2P and the Johnson’s SB functions yielded poorer results.

The mean values of the MSE in each diameter class for the fits to the observed data and for recovering from LiDAR data for the best (beta) and the poorest (Johnson’s SB) functions are shown in Figure 2a, while those for the generalized beta vs. Weibull-3P and for the gamma-2P vs. Weibull-2P functions are also shown in Figure 2b,c. Figure 2a shows that the MSE for the fits to the observed data was almost the same for both functions (beta and Johnson’s SB) over the diameter range studied, decreasing when the diameters increased. However, significant differences were observed for recovery of the functions. The recovered distributions yielded much larger values of MSE, up to 25 cm, after which both values (fitted and recovered) were more similar. The recovered Johnson’s SB yielded the largest values up to 19 cm, while the beta function was less accurate at between 19 and 24 cm. The difference between fitting and recovery of the smallest diameter classes increased with the number of data points.

Figure 2b shows the same trends for the range studied, with similar behavior obtained for fitting and recovery of both the generalized and Weibull-3P distributions. Figure 2c shows similar results produced by the gamma-2P and Weibull-3P functions.

Graphical assessment of the observed, described (fitted) and recovered distributions for six plots fits are shown in Figure 3. Pairwise comparisons of the functions were made: beta vs. Johnson’s SB (Figure 3a), generalized beta vs. Weibull-3P (Figure 3b) and gamma-2P vs. Weibull-2P (Figure 3c). The predictions were reasonable and followed the observed diameter distributions, especially for plots 1, 3, 4 and 6.

4. Discussion

We used a total of 58 TSPs from P. halepensis Mill. plantations for the present study, and LiDAR metrics were able be extracted from 50 of these plots. The number of plots is greater than in other studies of diameter distributions recovered from LiDAR data in the Northwest Iberian Peninsula [10,11]. These studies recovered the Weibull-2P and the Johnson’s SB distributions from LiDAR metrics. The Weibull function (2P and 3P) is the most commonly used distribution in this type of study [20,37,38]. However, in the present study, we used, for the first time within the framework of LiDAR-based research, three functions commonly used to describe and predict diameter distributions in forest stands, i.e., the beta [26], generalized beta [16] and gamma-2P [14] functions. In some cases, the results obtained were better than those obtained with the more usual Weibull (2P and 3P) and Johnson’s SB functions, for both fitting to the observed data and recovering with LiDAR metrics.

For the fitting phase, the results obtained for the main statistics used (D_n, W² and MSE) are similar to those obtained in [28] for Eucalyptus globulus stands in NW Spain. However, more accurate results were obtained in the same study for Pinus radiata stands by using the same functions. In a previous study, similar results were also obtained for Pinus sylvestris and better than those reported for Pinus pinaster stands in NW Spain [32].

The models used to recover the parameters of the distributions from moments of the distributions were accurate. Explanatory variables were LH_P90 (height percentile of 90%) for d_g, and LH_MAD_MEDIAN (median of the absolute deviations from the overall height median) for d_med. In other studies, height percentiles were also used as independent variables for estimating d_g (75% percentile and number of LiDAR last returns above a height of 1 m) and for d_med (1% percentile) in the models obtained in [10] for recovering the Weibull-2P distribution in 25 plantation plots of Pinus radiata in Northwest Spain, with R² for the d_g and d_med models of 0.80 and 0.77, respectively. For Pinus halepensis, we obtained similar R² values (0.75 and 0.72 for d_g and d_med, respectively). The RMSE values (2.54 cm for d_g and 2.55 cm for d_med) are consistent with the values reported in international forestry literature. For example, for German forests dominated by Picea abies (L.) Karst., the authors of [39] used data from a 0.44 pulse m⁻² LiDAR flight and reported an RMSE of 2.44 cm for d_med, while the authors of [40] studied a broad range of forest types (coniferous and hardwoods) and conditions across Ontario by using an artificially reduced LiDAR database of 0.5 pulses m⁻², reporting RMSE values ranging from 0.76 to 4.3 cm for d_g. The authors of [10] reported RMSE values of 3.42 for d_g and 3.63 for d_med, while the authors of [11] used exponential models instead of linear models to estimate the moments of the Johnson’s SB and the Weibull-2P functions, obtaining an R² of 0.82 and 0.86 for the d_med of Pinus radiata and Eucalyptus globulus respectively, and 0.84 and 0.89 for d_g of the same species. We compared the use of linear and exponential models for obtaining both variables and found that the results were similar, with the exponential model even including more independent variables. Thus, the linear models were considered more suitable.

Use of the four-parameter Johnson’s SB, beta and generalized beta distributions requires knowledge about the location parameter of the three distributions, the scale for the Johnson’s SB and the upper limit of the distributions for the beta and generalized beta functions. The three-parameter Weibull-3P requires knowledge about the location parameter, while in the Weibull-2P and gamma-2P functions, the value of the location parameter is zero. In functions with a location parameter, the value considered in the recovery phase with LiDAR metrics was the minimum diameter inventory (7.5 cm). This avoided having to model this parameter with LiDAR metrics, which can produce inaccurate results if stand variables are used in the model [9]. Thus, in all cases, the location parameter was considered equal to 7.5 cm for recovering the distributions from LiDAR data, while for fitting to the observed data, it was assumed to be 0.75·d_min, producing a similar mean value considering all plots of 7.18 cm (Table 3). The authors of [11] considered the location parameter for the Johnson’s SB to equal zero, thus avoiding modeling it with LiDAR metrics, while the scale parameter considered was also the maximum diameter (d_max).

The model obtained for d_max included, as LiDAR explanatory variables, LH_AAD (mean absolute deviation for height) and LH_P95 (95% height percentile), with R² = 0.87 and % RMSE = 10.56. The authors of [11] obtained similar values for Pinus radiata (R² = 0.93 and % RMSE = 8%) and Eucalyptus globulus (R² = 0.83 and % RMSE = 12%). The results for this variable modeled from LiDAR data are more accurate than those obtained in [27] and [9] with stand variables. The statistics obtained for recovery of the six functions were also reasonable. For example, the authors of [32] reported D_n values of 0.1924 for the Weibull-2P fitted by Maximum Likelihood, 0.2285 for the Johnson’s SB fitted by conditional maximum likelihood (CML) and 0.1812 for the beta fitted by moments to observed distributions of Pinus pinaster in Northwest Spain. The authors of [33] obtained a D_n value of 0.193 in the fits by maximum likelihood of the Weibull-3P to observed data of Pinus taeda plantations in the USA. The value obtained for the recovered beta from LiDAR metrics (0.2541) is close to these values. The number of KS acceptances are consistent with those found in previous studies. The percentage of acceptance (between 52% and 70%) was similar for Eucalyptus globulus in northwest Portugal [11] and higher than in Pinus radiata stands in northwest Spain [10,11].

5. Conclusions

Scant attention has been given to Aleppo pine plantations in Spain. However, we recovered the parameters of six probability density functions, i.e., the Weibull (2P and 3P), Johnson’s SB, beta, generalized beta and gamma-2P functions, from LiDAR metrics in 50 Pinus halepensis stands in Aragón (northeast Spain). The beta, generalized beta and gamma- 2P functions are of novel application in the field of LiDAR-based research. The results obtained were reasonable compared with previous studies and they showed that the beta and generalized beta functions were more accurate than the Johnson’s SB and Weibull (2P and 3P) functions, which are more commonly used in this type of study. None of the distributions predicted anomalous results.