Deriving Tree Size Distributions of Tropical Forests from Lidar

Abstract

Remote sensing is an important tool to monitor forests to rapidly detect changes due to global change and other threats. Here, we present a novel methodology to infer the tree size distribution from light detection and ranging (lidar) measurements. Our approach is based on a theoretical leaf–tree matrix derived from allometric relations of trees. Using the leaf–tree matrix, we compute the tree size distribution that fit to the observed leaf area density profile via lidar. To validate our approach, we analyzed the stem diameter distribution of a tropical forest in Panama and compared lidar-derived data with data from forest inventories at different spatial scales (0.04 ha to 50 ha). Our estimates had a high accuracy at scales above 1 ha (1 ha: root mean square error (RMSE) 67.6 trees ha⁻¹/normalized RMSE 18.8%/R² 0.76; 50 ha: 22.8 trees ha⁻¹/6.2%/0.89). Estimates for smaller scales (1-ha to 0.04-ha) were reliably for forests with low height, dense canopy or low tree height heterogeneity. Estimates for the basal area were accurate at the 1-ha scale (RMSE 4.7 tree ha⁻¹, bias 0.8 m² ha⁻¹) but less accurate at smaller scales. Our methodology, further tested at additional sites, provides a useful approach to determine the tree size distribution of forests by integrating information on tree allometries.

1. Introduction

Forests shape the Earth’s ecosystem, covering more than 30% of land area worldwide and storing about 45% of the terrestrial carbon [1,2,3]. They represent an important habitat for biodiversity and are relevant for economy and society [4,5]. At the same time, they are increasingly affected by climate warming and weather extremes [6,7,8], logging and fire [5], deforestation and fragmentation [5,9,10,11,12].

To conserve forests, their current status and future development must be monitored. In the past, the state of a forest has been determined using field measurements of sample plots (e.g., 1 ha). Recent studies, in turn, gathered large-scale plot data (e.g., 50 ha), regularly collecting detailed information such as the stem diameter or the species and location of each tree [13]. Over recent decades, remote sensing has become an increasingly used tool to monitor the state of forests [14,15].

Remote sensing techniques like light detection and ranging (lidar) (terrestrial, airborne or spaceborne, e.g., Ice, Cloud, and land Elevation Satellite (ICESat) [16], Global Ecosystem Dynamics Investigation (GEDI) [17]), radar (e.g., TerraSAR-X-Add-on for Digital Elevation Measurements (TanDEM-X) [18], Advanced Land Observing Satellite (ALOS), Phased Array type L-band Synthetic Aperture Radar (PALSAR) [19]) and optical reflectance (Landsat [20], Sentinel [21] or Moderate-resolution Imaging Spectroradiometer (MODIS) [22]) provide information on forests at different spatial and temporal resolutions. Earth observations have been used to gain insight on (a) forest cover, (b) forest biomass, and (c) forest structure. For example, global forest cover has been explored from Landsat data at the 30-m resolution [3] and has been used for the analysis of deforestation, forest fragmentation [9], and their respective drivers [10].

Often, remote sensing observations have been analyzed by computing correlations with forest attributes (e.g., lidar metrics correlated with biomass or basal area) [23,24,25,26,27,28]. Such approaches have a high predictive power, but typically focus on predicting single forest attributes (like biomass) and require field-data-intensive calibration of site-specific parameters in the statistical relationships. Studies also combined statistical methods with information on tree geometry [29,30,31]. A good example for an approach that integrates tree geometry is Spriggs et al. [30], who combined allometric relations with a Bayesian optimization approach.

Recent studies have followed the approaches of tree segmentation [32,33,34,35,36] by detecting single tree crowns from airborne lidar to estimate tree size attributes and frequencies. Those techniques interpret lidar point clouds and allow the detection of tree geometry relations (in particular for large trees). Similar approaches also use information from terrestrial lidar [37] with its strength to detect detailed tree leaves and branches, particularly in the lower canopy (often constrained to small forest plots). Both approaches require high-resolution lidar point clouds to have sufficient information on the trees.

Here, we present a novel approach for inferring forest structure from lidar measurements that is complementary to the recent approaches. Our approach integrates the available information on tree geometry in a theoretical model. Through its combination with the measured vertical lidar profile (e.g., derived from waveforms or lidar point clouds), we calculated how many trees of a certain size could occur in a forest to match the measured profile (Figure 1). This allowed us to estimate the tree size distribution of a forest (i.e., number of trees in different stem diameter classes) based on lidar measurements.

We tested our approach on a 50-ha forest plot in Panama. After determining the leaf–tree matrix, we derived the stem diameter distribution from lidar by applying classical methods of mathematics (matrix inversion). We validated the estimated stem diameter distribution by using field inventory data and demonstrated the accuracy of our approach at different spatial scales (50, 25, 5, 1, 0.25, and 0.04 ha). We analyzed how reliably we can estimate tree numbers per stem diameter class. In addition, we derived estimates for the total stem number and basal area from the lidar measurements.

2. Materials and Methods

2.1. Study Site

We applied our approach to a tropical forest located on Barro Colorado Island in Panama (BCI, 9°9′ N, 79°51′ W) [38,39,40,41], where a 50-ha plot of old-growth tropical moist forest has been monitored every five years since 1980. The study site is located in a natural reserve area in the center of the island. In each census, each tree was tagged and mapped, its species was identified and its stem diameter was measured at breast height. In our study, we used the census of the year 2010, considering all trees of stem diameter ≥ 1 cm. The analyzed forest census included 244,269 trees and 301 tree species (1.5% of the originally recorded trees were disregarded because they were dead or lost, see Supplementary S1). The forest inventory data can be requested from the ForestGeo Global Earth Observatory Network [42] or downloaded for free from Dryad [43].

To relate the height and crown dimension of a tree to its stem diameter, we used three allometric equations:

(a)

tree height h (m) is related to the stem diameter d (m) by an asymptotic function [24,44]

h = (57.4 × d)/(0.43 + d),

(1)

(b)

tree crown radius cr (m) is related to the stem diameter d (m) by a power law [23,45]

cr = 9.08 × d^0.68

(2)

and

(c)

vertical tree crown length cl (m) is linearly related to the tree height by [23]

cl = 0.4 × h.

(3)

The tree crown radius, length and the assumed shape were used to calculate the tree crown volume, which is multiplied with a density factor of ρ = 0.44 m²/m³ to estimate the tree leaf area [23]. The tree height and crown length were used to estimate the leaf area of a tree at different heights within the forest canopy (we assumed that trees’ leaves are homogenously distributed within their crown).

Tree allometries have been derived from independent field studies [23,24,44,45]. Here, we averaged tree allometries across occurring tree species, which has been shown to be a valid approach when studying the forest structure [46]. In the standard case, we assumed that tree crowns have an ellipsoidal shape, but we tested different crown shapes (spheres and cylinders) as well. We further tested the sensitivity of different allometries for tree height (non-asymptotic power law of tree height [45])

h = 43.4 × d^0.6

(4)

and of leaf density within crowns (factor of ρ = 1 m²/m³).

Airborne lidar data were available for the 50-ha forest plot for the year 2009 [47]. The point densities ranged from 0 to 103 m⁻² (in flight swath overlaps) with a median of 19 m⁻². The point cloud was terrain-normalized and thinned using random subsampling to obtain a near homogeneous point density of four returns per m². The latter was achieved by iterative subsampling to different point densities and inspection of the resulting density rasters until no further density differences from the flight pattern were visible. The thinned point cloud contained 19% of the returns from the original point cloud. For details on the instrumentation and processing of the lidar data, see Lobo and Dalling [47] (Materials and Methods). The lidar dataset is publicly available and can be requested from J.W. Dalling (see statement on data accessibility in [47]).

2.2. Derivation of Leaf Area Profiles from lidar

The lidar data (i.e., lidar profile) includes information on the site-specific forest structure as lidar signals get reflected by the leaves and branches of trees. We assumed that the point density PD (1/m²) of the lidar profile in a specific height layer (h_(i−1),h_i) of constant width Δh = 1 m is a result of (a) the tree crown density (LAD) at which lidar signals could get reflected and (b) the probability (W) that the lidar signal could penetrate into the respective height layer i (i.e., has not already been reflected in the layers above):

PD_i = l × LAD_i × W_i

(5)

The parameter l represents a density factor (here, l = 1 m²/m³) that combines different aspects: the density of lidar shots sent to the forest canopy, the fraction of tree leaves (relative to tree branches) at which signals can be reflected and the density of successful lidar signals returned at the top of the forest (due to gaps in the canopy).

The observed lidar point density PD_i in (m^-²) is then defined as

$$P{D}_i= \left(1/A\right) \times {{\displaystyle \sum }}_{k=1}^K{\delta }_{ik}$$

(6)

with the plot area A (m²), the total number of lidar points K and

$${\delta }_{ik}=\{\begin{array}{cc}1& h_{i-1}{\le z}_k{<h}_i\\ 0& \mathrm{else}\end{array}$$

(7)

The sum of the function δ_ik yields the number of lidar points with the return height z_k (m) in the respective height layer i. To compute the probability W, we used the Beer–Lambert law of light transmission and extinction [48].

$$W_i = exp\left(-k\times \Delta h\times {{\displaystyle \sum }}_{j=i+1}^{n}LA{D}_j\right)$$

(8)

where ΣLAD_j is the cumulative leaf area density (m²/m³) from the top of the forest (j = n) down to the height layer j = i+1. The parameter k in Equation (8) represents the average light extinction coefficient (k = 0.2 for near-infrared or infrared signals with wavelengths of 600–1400 µm) [48].

Based on Equation (5), we can then recursively estimate the forest’s leaf area profile (LAD^lidar) via

$$LA{D}_i^{lidar}=\frac{P{D}_i}{l\times W_i}$$

(9)

starting with a height layer above the forest canopy (i.e., LAD_n = 0 results in W_n = 1). We did not include ground returns in our approach and considered only height layers ≥ 3 m. The values of the parameters k and l are relevant for the derived leaf area profile, but had only minor influence on the estimation of the stem diameter distribution in our study (Figure A1).

2.3. ‘Leaf–Tree Matrix’ of the Tree Geometry Model

To estimate the number of trees for different stem diameter classes (here, assumed to be unknown) from the lidar-derived leaf area profile, we developed a theoretical model that integrates the available information on tree geometry and allometry.

For this, we introduced the ‘leaf–tree matrix’. The matrix summarizes the consequences of the assumed tree allometries (see Section 2.1) in a virtual forest in which exactly one tree is present for each stem diameter class. Figure 2 visualizes this matrix, which contains for each stem diameter class (columns) and height layer (colored rows) the corresponding leaf area per tree. White cells indicate that trees in the corresponding stem diameter class do not have leaves in the respective height layer. Height layers are defined in 1-m width, and stem diameter classes (of increasing width) correspond to the respective height layers (based on the assumed stem diameter-height relation, Section 2.1). Each tree has a maximum height (upper boundary of the height layer, based on the tree height allometry and diagonal entries of the matrix). Crowns reach only into the lower height layers (according to the tree crown allometry). Leaf area was assumed to be proportional to the tree crown volume (based on the tree crown allometry) and distributed uniformly among the crown-covered height layers (Section 2.1). See Figure A2 for the leaf–tree matrix with different assumed crown shapes and model assumptions.

As a consequence, the vertical leaf area distribution for a forest with only one tree per stem diameter class can be derived by summing up all column entries (i.e., stem diameter classes) per row (i.e., height layers) in the leaf–tree matrix.

We denote the leaf–tree matrix by F = (f_ij) with i,j = 1, …, n and the number of trees for each stem diameter class i by the vector N = (N_i)_i=1,...,n (with n being the number of stem diameter classes, here n = 55). Multiplying the leaf–tree matrix F with the vector N results in the leaf area profile of a forest plot (LAD in m²/m³), in which we normalize by the plot area A (m²) and height layer width Δh (m)

$$LAD= \left(F \times N\right)/\left(A\times \Delta h\right)$$

(10)

We aim to develop an approach to infer the stem diameter distribution based on lidar measurements. Therefore, we assumed that sufficient information about forest structure was included in the lidar-derived leaf area profile (LAD^lidar of Equation (9)), which was used for the LAD (of Equation (10)). Rearranging Equation (10) yields the stem diameter distribution (unknown vector N).

2.4. Linear Equation Solving to Derive Forest Structure from lidar Profiles

Equation (10) represents a system of n linear equations with n unknowns. To rearrange this equation for inferring the unknown stem diameter distribution (vector N), the inverse of the leaf–tree matrix F is required. Since, by construction, the height of trees in the n-th stem diameter class never exceeds the n-th height layer, all above-diagonal matrix entries are zero. This makes the leaf–tree matrix F a lower triangular matrix, which is invertible (n × n matrix) [49]. Rearranging Equation (10) for the vector N then results in

$$N^{lidar} = F^{-1} \times \left(LA{D}^{lidar} \times A\times \Delta h\right)$$

(11)

Though Equation (11) yields the number of trees (N) for each stem diameter class, some of the computed values can be negative. Therefore, we developed an iterative numerical backward approach for solving the linear equation system of Equation (10) (using LAD^lidar to represent LAD). See Appendix A

Table A1. Results of different calculations (direct: calculation of stem numbers N from Equation 11, numerical backward: numerical solving of Equation 10 for stem numbers N, see methods and Appendix B) to estimate the stem diameter distribution from lidar measurements (50-ha scale)

Method	Direct Calculation	Numerical Backward Solving
Regression slope	1.22	1.24
Coefficient of determination R²	0.88	0.89
Number of stem classes with negative predictions (%)	26.9	0
RMSE (trees/ha)	29.6	22.8
nRMSE (%)	8.1	6.2

for a comparison of results between the direct and numerical calculation and Appendix B for details on the algorithm of numerical calculation.

2.5. Analysis and Statistics of Results

We evaluated the accuracy of our approach by comparing the estimated stem diameter distribution with forest inventory data (i.e., stem numbers per ha and stem diameter class). To compare the stem diameter distribution observed in the field $\left(N^{field}\right)$ with that estimated from lidar data and our tree geometry model $\left(N^{lidar}\right)$, we performed a logarithmic regression analysis, setting $\ln \left(N^{lidar}\right)=I+s\times \ln \left(N^{field}\right)$ with the intercept I and the regression slope s. We analyzed the regression slope s and the coefficient of determination R² for stem diameter classes with stem numbers larger than zero only (in both inventory and estimates). Optimal values of R² = 1 and s = 1 would reflect the highest achievable accuracy of the lidar -derived tree numbers.

We further calculated the root mean square error RMSE (in trees per ha)

$$RMSE=\frac{1}{A_{ha}}\times \sqrt{\frac{{{\displaystyle \sum }}_{j=1}^n{\left(N_j^{lidar}-N_j^{field}\right)}^2}{n}}$$

(12)

and its normalized counterpart nRMSE (in %, using the range of censused stem numbers across all stem diameter classes to normalize)

$$nRMSE=A_{ha}\times \frac{RMSE}{\left(\underset{j}{\max }{N}_j^{field}-\underset{j}{\min }{N}_j^{field}\right)}$$

(13)

with $A_{ha}$ being the plot area in ha. We calculated both RMSE and nRMSE for (non-logarithmic) stem numbers aggregated to 10 cm stem diameter classes and for trees with stem diameter d ≥ 10 cm.

We further estimated the total tree density (N) and basal area (BA) of the forest based on the lidar-derived stem diameter distribution (using stem diameter classes of the leaf–tree matrix and class mid values for basal area estimates) and compared them with the field inventory. Again, we focused on classes with stem diameter above 10 cm. In addition, we compared the tree density and basal area estimates also for each 10-cm aggregated stem diameter class (j) with the inventory

$${\Delta }_j^N=\frac{1}{A_{ha}}\times \left|N_j^{lidar}-N_j^{field}\right|$$

(14)

$${\Delta }_j^{BA}=\left|B{A}_j^{lidar}-B{A}_j^{field}\right|$$

(15)

To assess the accuracy of our approach on different spatial scales, we varied the size of the analyzed plot. First, we applied our approach at the scale of 50 ha (entire forest plot). In a second step, we divided the forest plot into subplots of equal sizes of, respectively, 25, 5, 1, 0.25, and 0.04 ha. Subplot sizes of 25 × 25 m or lower are especially interesting as they correspond to the size of the footprint of satellite-based lidar observations (GEDI [17]). For a specific plot size, we calculated, for each subplot, the nine statistical measures (regression slope, R², RMSE, nRMSE, tree density, and basal area in total and absolute differences per 10-cm stem diameter class). For each spatial scale, we calculated the arithmetic mean, standard deviation, and range of the statistical measures across the respective subplots. For the mean values of R², regression slope, RMSE, and nRMSE, we used linear regression to identify scaling relations expressing these statistical quantities as functions of the plot size. To compare the estimated tree density and basal area with inventory data, we used linear regression, computed the bias (average difference of field from lidar-derived attribute), RMSE and nRMSE (normalized by the mean of the observed attribute).

We further correlated the R² values we obtained for the estimated stem diameter distribution (per subplot) (a) with the other obtained statistical measures (RMSE, nRMSE, and regression slope for the stem diameter distribution), (b) with characteristics of the lidar-derived leaf area profile (e.g., median profile height), and (c) forest attributes (e.g., tree density and basal area). For a full list of attributes used for the correlations, see

Table 1. List of attributes used for the correlations with R² (from the comparison of the estimated stem diameter distribution with field data). Columns denote the data used, a description of the attribute, its unit and the different statistical calculations applied.

Data Source	Attribute	Unit	Calculations
lidar data	Profile height	m	Maximum, Median, Variance
	Leaf area density	m²/m³	Maximum, Median, Variance
	Number of lidar returns	1/m2	-
	profile height weighted by leaf area	m	Median (WMPH 1), Variance (WVPH 2)
Forest inventory	Basal area 3,4	m²/ha	Sum
	Tree density 3,4	1/m2	Sum
	Tree height 3	m	Median, Standard deviation

¹ Median height of the vertical leaf profile derived from lidar (weighted by the derived leaf area density LAD^lidar per height layer). ² Variance of height of the vertical leaf profile derived from lidar (weighted by the derived leaf area density LAD^lidar per height layer). ³ Including trees of stem diameter d ≥ 1 cm. ⁴ Including trees of stem diameter d ≥ 10 cm.

. We evaluated each correlation by calculating Spearman’s rank correlation coefficient r.

3. Results

3.1. Stem Diameter Distribution Derived at the 50-ha Scale

The stem diameter distribution of the 50-ha forest plot shows a typical heavy-tailed pattern, as often observed in old-growth forests, with few large trees and numerous small ones (Figure 3a, linear appearance on log–log axes). To test our approach at the 50-ha scale, we derived the stem diameter distribution from lidar and compared the estimated tree numbers per stem diameter class with the observed numbers from the inventory (Figure 3). The stem diameter distribution was estimated from lidar with a high accuracy (RMSE = 22.8 tree per ha, nRMSE = 6.2%, R² = 0.89, Figure 3, Table A1).

The number of mid-sized and large trees was estimated accurately (10 cm < stem diameter ≤ 50 cm with RMSE = 39.4 trees per ha, stem diameter > 50 cm with RMSE = 1.2 trees per ha) while trees smaller than 10 cm in stem diameter were overestimated (RMSE = 590.5 trees per ha, Figure A3a). Based on the lidar-derived stem diameter distribution, the estimated tree density of 516.7 trees per ha (stem diameter d ≥ 10 cm) slightly overestimated the censused density of 447.3 trees per ha (bias of −69.4 trees per ha). The basal area estimates showed in turn, a slight underestimation (inventory: 30.1 m²/ha, lidar: 23.6 m²/ha, bias: 6.5 m²/ha, Figure A3a). A sensitivity analysis revealed only minor influences of differently assumed crown shapes, tree height allometry or the leaf density within crowns (

Table 2. Sensitivity analysis of different model assumptions (see also Figure A2). The changed model assumptions for each scenario are written in bold.

Sensitivity Scenario	1	2	3	4
Crown shape	Sphere	Cylinder	Ellipsoid	Ellipsoid
Height allometry	Asymptotic	Asymptotic	Asymptotic	Power law
Crown length	Crown radius	0.4 × height	0.4 × height	0.4 × height
Crown leaf density	0.44 m²/m³	0.44 m²/m³	1 m²/m³	0.44 m²/m³
Regression slope	1.22	1.29	1.21	1.1
R²	0.89	0.98	0.92	0.89
RMSE (trees/ha)	32.0	4.5	15.0	15.3
nRMSE (%)	8.8	1.2	4.1	3.8

, Figure A2).

3.2. Small-Scale Derivations of Stem Diameter Distributions

We tested our approach also for smaller plot sizes (25, 5, 1, 0.25, and 0.04 ha,

Table 3. Comparison of forest attributes derived from lidar and from the forest inventory for different plot sizes. Results are shown for the regression analysis (slope and R²), RMSE and nRMSE (overall and for three different categories of small, mid-sized and large trees grouped according to stem diameter d) and aggregated forest attributes (tree density and basal area, for stem diameter d ≥ 10 cm). For each plot size, the mean ± standard deviation (and in brackets, minimum and maximum) of the respective attribute is given. The plot size is given in ha with the number of plots and plot dimension in m × m. Forest attributes are compared between lidar and inventory also in terms of the bias, RMSE, and nRMSE (%) (see methods for details). See Table A2 for additionally analyzed plot sizes.

Plot Size (ha)		50	1	0.25	0.04
Side length (m) in x-direction		1000	100	50	20
Side length (m) in y-direction		500	100	50	20
Number of plots		1	50	200	1250
Overall quality	Regression slope	1.24	0.9 ± 0.14 (0.66, 1.23)	0.76 ± 0.2 (0.11, 1.3)	0.55 ± 0.38 (−0.94, 2.17)
	R²	0.89	0.76 ± 0.13 (0.38, 0.95)	0.67 ± 0.17 (0.13, 0.94)	0.44 ± 0.27 (0.00005, 1)
	RMSE (trees/ha)	22.8	67.6 ± 45.1 (15.2, 187.4)	118.1 ± 77.6 (13.4, 385.6)	219.2 ± 160.4 (17.7, 1074.6)
	nRMSE (%)	6.2	18.8 ± 13.2 (4.2, 55.0)	33.2 ± 21.7 (3.0, 105.7)	70.9 ± 62.0 (3.2, 590.1)
RMSE per size group (trees/ha)	Small trees 1	590.5	719.5 ± 379.6 (188.6, 1844.9)	796.8 ± 492.6 (153.6, 3028.1)	902.2 ± 751.6 (128.4, 5423.6)
	Mid-sized trees 2	39.4	77.1 ± 30.5 (28.0, 175.3)	112.9 ± 43.1 (39.9, 267.6)	153.9 ± 73.9 (33.5, 606.7)
	Large trees 3	1.2	5.1 ± 2.5 (2.2, 15.9)	9.1 ± 5.5 (2.4, 36.1)	26.2 ± 21.2 (0, 230.3)
Tree density 4 (trees/ha)	Inventory-based	447.3	447.3 ± 46.1 (353, 597)	447.3 ± 60.4 (300, 648)	447.3 ± 115.7 (175, 1025)
	lidar-derived	516.7	615.5 ± 159.3 (364, 1320)	758.9 ± 249.2 (316, 1864)	992.6 ± 402.5 (150, 3125)
	bias	−69.4	−168.2	−311.6	−545.3
	RMSE	69.4	228.4	394.4	677.8
	nRMSE (%)	15.5	51.1	88.2	151.5
Basal area 4 (m²/ha)	Inventory-based	30.1	30.1 ± 5.1 (20.4, 45.8)	30.1 ± 8.3 (15.5, 67.9)	30.1 ± 20.9 (3.5, 206.1)
	lidar-derived	23.6	29.3 ± 6.9 (18.8, 58.1)	38.1 ± 15.0 (9.5, 163.4)	59.0 ± 45.0 (4.7, 762.5)
	bias	6.5	0.8	−8.0	−28.9
	RMSE	6.5	4.7	14.3	49.8
	nRMSE (%)	21.6	15.7	47.4	165.3

¹ Stem diameter classes d ≤ 10 cm. ² Stem diameter classes 10 cm < d ≤ 50 cm. ³ Stem diameter classes d > 50 cm. ⁴ Stem diameter classes d ≥ 10 cm.

,

Table A2. Comparison of the forest attributes derived from lidar and from the forest inventory for different plot sizes. Results are shown for the regression analysis (slope and R²), RMSE and nRMSE (overall and for three different categories of small, mid-sized and large trees grouped according to the stem diameter d) and aggregated forest attributes (tree density and basal area, for stem diameter d ≥ 10 cm). For each plot size, the mean ± standard deviation (and in brackets, the minimum and maximum) of the respective attribute is given. The plot size is given in ha with the number of plots and plot dimension in m × m. Forest attributes are compared between lidar and inventory also in terms of the bias, RMSE, and nRMSE (%) (see methods for details).

Plot Size (ha)		25	5
Side length (m) in x-direction		500	200
Side length (m) in y-direction		500	250
Number of plots		2	10
Overall quality	Regression slope	1.17 ± 0.05 (1.14, 1.21)	1.0 ± 0.1 (0.8, 1.1)
	R²	0.92 ± 0.005 (0.91, 0.92)	0.84 ± 0.09 (0.65, 0.94)
	RMSE (trees/ha)	25.4 ± 2.3 (23.8, 27.0)	34.1 ± 11.5 (17.1, 50.8)
	nRMSE (%)	7.0 ± 1.1 (6.2, 7.7)	9.4 ± 3.3 (4.9, 14.1)
RMSE per size group	Small trees 1	596.2 ± 192.2 (460.3, 732.1)	622.6 ± 240.2 (402.2, 1191.5)
	Mid-sized trees 2	36.0 ± 1.4 (35.0, 37.0)	45.9 ± 7.3 (36.8, 55.0)
	Large trees 3	2.1 ± 0.6 (1.6, 2.5)	3.3 ± 1.4 (1.5, 5.1)
Tree density 4 (trees/ha)	Inventory-based	447.3 ± 3.0 (445.2, 449.4)	447.3 ± 25.1 (420.0, 489.4)
	lidar-derived	562.3 ± 29.4 (541.5, 583.1)	561.4 ± 65.4 (476.2, 715.4)
	bias	−115.0	−114.1
	RMSE	116.5	128.3
	nRMSE (%)	26.1	28.7
Basal area 4 (m²/ha)	Inventory-based	30.1 ± 2.4 (28.4, 31.8)	30.1 ± 2.8 (25.5, 34.7)
	lidar-derived	24.3 ± 1.3 (23.4, 25.2)	25.3 ± 3.4 (19.6, 31.8)
	bias	5.8	4.7
	RMSE	5.8	4.9
	nRMSE (%)	19.4	16.4

¹ Stem diameter classes d ≤ 10 cm. ² Stem diameter classes 10 cm < d ≤ 50 cm. ³ Stem diameter classes d > 50 cm. ⁴ Stem diameter classes d ≥ 10 cm.

). For most plot sizes (e.g., at the 1-ha and 0.25-ha scale), the stem diameter distributions were estimated well from the lidar data (Figure 4a–f, e.g., average R² values of 0.76 and 0.67, average RMSE values of 67.6 and 118.1 trees per ha, Figure A4). Few plots at the 0.25-ha scale showed lower R² (down to 0.13), higher RMSE values (up to 385.6 trees per ha) or higher nRMSE values (up to 105.7%) (Figure A4). At the 0.04-ha scale, more plots with lower accuracy of the estimated stem diameter distribution occurred; nevertheless, we obtained a mean R² of 0.44 (mean RMSE of 219.2 trees per ha, mean nRMSE of 70.9%) (Figure 4i, Table 3, Table A2). Similar results can also be observed for the regression slope (Figure A4). Interestingly, the R² values scaled with plot size (Figure 5, similarly also RMSE, nRMSE, and regression slope) and correlated non-linearly with the other statistical quantities (

Table A3. Correlation coefficients (Spearman’s r) of R² values (of the compared lidar-derived and inventoried stem diameter distribution) correlated with other statistical measures, lidar data measures, and forest structure attributes for different plot sizes (1, 0.25, and 0.04 ha). See methods and Table 1 for details.

Category	Measure	Spatial Scale
		1 ha	0.25 ha	0.04 ha
Statistical measures	Regression slope	0.63	0.54	0.73
	RMSE	0.09	−0.15	−0.20
	nRMSE (%)	0.05	−0.19	−0.18
lidar data	Profile height—median	−0.06	−0.09	−0.06
	Profile height—maximum	−0.03	−0.10	−0.05
	Profile height—variance	−0.05	−0.09	−0.05
	Leaf area density—median	−0.13	−0.03	−0.01
	Leaf area density—maximum	0.17	−0.03	−0.03
	Leaf area density—variance	0.21	0	−0.04
	WMPH	−0.31	−0.25	−0.19
	WVPH	−0.11	0.03	0.08
	Number of lidar returns	−0.12	0.09	−0.04
Forest inventory	Basal area (m²/ha)	−0.15	−0.14	−0.09
	(for stem diameter d ≥ 10 cm)	−0.15	−0.16	−0.10
	Tree density (ha−1)	0.12	0.27	0.21
	(for stem diameter d ≥ 10 cm)	0.17	−0.02	−0.06
	Standard deviation of tree height (m)	−0.16	−0.28	−0.19
	Median tree height (m)	0.09	−0.10	−0.09

; however, the correlations of R² with RMSE or nRMSE were less strong, r = 0.09 and 0.05 at 1-ha scale).

Estimates on the tree density and basal area showed good results for the 1-ha scale with nRMSEs of 51.1% and 15.7%, respectively (bias of −168.2 trees/ha and 0.8 m²/ha, Figure 6, Table 3). The estimation capability for the tree density and basal area per stem diameter class did not change with decreasing plot size (mostly for small-sized classes) but showed larger variations for smaller plots (Figure A3c,d). Interestingly, the forest basal area was best estimated at the 1-ha scale (Figure 6, Figure A5, Table 3, underestimated for larger plot sizes, Table A2), while the tree density was best estimated at the 50-ha scale (Figure 6, Table 3). Both attributes were overestimated by lidar at smaller spatial scales (Table 3, Table A2).

To understand how well stem numbers can be derived from lidar using our approach (e.g., in case no inventories are available), we correlated the R² values of comparing inventory- and lidar-derived stem diameter distributions with other measures characterizing each plot for the spatial scales of 1, 0.25, and 0.04 ha. Those measures describe either the lidar profile or the forest structure (Table 1). We found correlations for the median height of lidar returns WMPH (weighted median profile height, weighted by the derived leaf area density LAD^lidar, r = −0.19 to −0.31, Table A3) and total tree density (stem diameter d ≥ 1 cm, r = 0.12 to 0.27, Table A3, Figure 7).

The WMPH became more important with increasing spatial scale (as well as variance of leaf area density), while the total tree density showed less strong correlation at larger scales. Interestingly, the tree density, only for stem diameter d ≥ 10 cm, revealed an opposite trend with stronger correlation at the 1-ha scale than at smaller scales (Table A3). This suggests that forest plots characterized by a low and dense canopy could enable a more reliable estimation of the stem diameter distribution from lidar data. In addition, the forest heterogeneity (standard deviation of tree height) also plays a relevant role (r = −0.16 to −0.28, Figure 7c, Table A3). Surprisingly, a reliable estimation of the stem diameter distributions (in terms of high R² values) showed weak correlations with the number of lidar returns, median tree height, and forest basal area (Table A3).

4. Discussion

We presented here a novel approach for estimating the tree size distributions from lidar measurements. Our approach is solely based on the vertical profile of lidar data and infers the number of trees per stem diameter class by integrating the available tree size allometries. On the example of a tropical forest in Panama, we demonstrated that the presented method was able to estimate the stem diameter distribution with high accuracy not only at the 50-ha scale (RMSE of 22.8 trees per ha, nRMSE of 6.2%, R² of 0.89), but also at the 1-ha scale (average RMSE of 67.6 trees per ha, mean R² of 0.76). The basal area was estimated well at the 1-ha scale (bias 0.8 m² ha⁻¹), but the bias increased at smaller scales (−8.0 m² ha⁻¹ at 0.25-ha scale). We identified the most decisive factors for a good estimation of stem diameter distributions at smaller spatial scales as (a) the mean profile height (WMPH), (b) tree density, and (c) forest heterogeneity (in terms of tree height).

The change from a still good estimation of stem diameter distributions at larger spatial scales (50 ha to 0.25 ha) toward slightly biased estimates at the 0.04-ha scale probably could be related to the crown overlapping of large trees and edge effects (crowns from trees of neighboring plots). This might also affect the identified correlations between the goodness of fit and forest attributes, which are not strong (maximum values of 0.3 or −0.3, respectively). Forest inventory data including explicit tree locations allowed us to calculate the overlapping crown parts from neighboring plots and, thus, enabled us to quantify and to account for such effects. Furthermore, we demonstrated our approach on the example of an old-growth tropical forest site (at BCI, Panama). The results of small-scale estimations and the identified correlations could differ for forests of other biomes (like temperate or boreal forests) or for disturbed or managed sites (e.g., logged forests or forest plantations).

4.1. Strengths and Limitations of the Presented Approach

Our approach complements recent methods [23,24,25,26,27,28,29,30,31,32,33,34,35,36] as it (i) integrates already available information on tree geometry (captured in the leaf–tree matrix) and (ii) estimates the full range of size classes in tree size distributions. Tree size distributions can be used to derive different forest attributes (e.g., tree density or forest basal area) and allow us to estimate a forest’s successional state or disturbances [50]. By this, our approach extends statistical methods that often focus on single forest attributes. In comparison to approaches that also integrate tree geometry with statistics, our approach obtained, for almost 85% of the analyzed plots at the 0.25-ha scale, an R² > 0.5 (e.g., Spriggs et al. [30] reached 73% of 0.25-ha plots with an R² > 0.5 in temperate forests). The presented approach complements tree segmentation methods. Such methods show their strength in predicting the size and geometry of large trees, while our approach effectively predicts small- and mid-sized trees in the forest, also at larger spatial scales.

Our approach is based on information about tree allometries, which are available for a wide range of different forest types and biomes (e.g., [44,51], information can also be found in national forest inventories). Tree allometries are assumed a priori in our approach and can differ across regions and forest types. This approach may not be suitable especially for large trees, leaning trees or trees with a complex geometry. We identified only minor influences of these model assumptions on the estimated stem diameter distribution for our study site (Figure A2, Table 2). We demonstrated a good estimation of the forest basal area and total tree density based on lidar-derived tree numbers for stem diameters larger than 10 cm.

To improve the allometries, crown size relations can also be derived by tree crown detection algorithms based on lidar point clouds [32,33,52]. For example, Ferraz et al. [33] reported tree and crown allometries based on detected single trees and their respective uncertainties. This approach requires high-resolution lidar point clouds, particularly for correctly detecting small- and mid-sized trees (with small crowns) in the understory where the lidar point density is normally lower. While our approach shows its strength in estimating tree numbers of small- and mid-sized trees (based on the vertical lidar profile), tree crown detection approaches are especially interesting for applications at smaller spatial scales. Methods based on crown delineation and classification approaches or using terrestrial lidar measurements will, in the future, provide further knowledge on the crown shape of single trees (e.g., National Ecological Observatory Network (NEON)) [53,54].

Our approach includes only average tree size allometries across occurring tree species. We already demonstrated, in a previous study [46], that forest size structure can be explained by using average tree allometry. However, forest sites where crown shapes, in particular, (and the amount of leaf area per crown volume) differ more between tree species will be more difficult to analyze and will require a closer look. In this case, our method could be extended toward a linear combination of multiple leaf–tree matrices (one for each dominant tree species) if information on the species composition is available (e.g., by optical sensors) [55,56]. A recent study [31] demonstrated how stem diameter distributions of different plant functional types can be derived based on site-specific calibration.

Our method is based on vertical lidar profiles only (instead of the full lidar point cloud) and, thus, appears promising for applications to other large footprint lidar data. Previous studies demonstrated the comparability of profiles derived from lidar point clouds and airborne large-footprint waveforms (Land, Vegetation, and Ice Sensor (LVIS)) [57,58]. Two parameters of our methodology (the extinction parameter k and density factor l) for deriving the vertical leaf profile from lidar profiles, could depend on the used remote sensing technique.

While the extinction coefficient k reflects the light extinction in the forest canopy, which depends on the wavelength (e.g., optical vs. near infrared) [48], the density factor l is more phenomenological and might correlate with the used tree allometries. Terrestrial lidar could contribute to the estimation of l for forest plots by reconstructing the detailed tree crown architecture of branches and leaves. We identified only minor influences of both parameters (k and l) on the estimated stem diameter distribution for our study site (Figure A1).

Other methods to derive leaf density profiles from lidar data include similar parameters in their approaches [59,60,61,62]. For example, the method used by Stark et al. [60] and Harding et al. [61] is based on the (logarithmic) ratio of point densities of two adjacent height layers [63] and includes an extinction coefficient. Tang et al. [59] combined the MacArthur–Horn-method with a gap probability approach [64]. The approach of Detto et al. [62], in contrast, is based on a stochastic radiative transfer model and accounts for multiple lidar returns. Methods, like small-footprint full-waveform lidar analyses, could improve estimates on leaf density profiles in the future [65].

4.2. Future Applications and Challenges

Airborne campaigns of lidar measurements have the advantage of providing forest measurements for hundreds of hectares, but are limited to the regional scale and are not able to cover continents. Satellite missions can provide full-waveform lidar data but only provide samples of space (e.g., GEDI has a footprint of 25 m diameter and spacings of 60 m along track and 600 m across track) [17]. Assuming that such spaceborne-generated profiles are comparable with profiles derived from point clouds or airborne-generated waveforms [57], the presented approach could be tested for different remote sensing data and for forest sites at which inventory and spaceborne lidar data is available.

This could be particularly interesting for exploring the site-dependency of tree allometries (and the parameters k and l) in our approach and testing its applicability for a wide range of forest types (especially regions characterized by a heterogeneous forest structure). As our results showed lower quality at smaller spatial scales (of e.g., 0.04 ha), questions on sampling efforts also arise. Do the results improve if we use several plots distant to each other (equivalent to spaceborne lidar footprints)? More specifically, it would be interesting to explore if we can aggregate discontinuously sampled lidar shots over an area to estimate the stem diameter distribution. How many lidar shots are required to be sampled to obtain good estimates of a specific forest site?

Recent studies (based on individual-based forest growth models) [31,66,67] demonstrated how remote sensing information can be combined with forest models to provide additional information about forests. Forest models simulate the long-term dynamics of forests (based on single trees) and can estimate forest successional states by matching simulated and remotely sensed vertical leaf profiles [66]. Our approach—which is based on related methodological assumptions—thus, also appears promising for estimating the successional state (caused, e.g., by disturbances).

As we demonstrated that stem diameter distributions in forest plots with a low dense canopy and low tree height differences can be estimated reliably, we expect that disturbances that occurred previously could be well identified. Challenges for spaceborne lidar could make such an analysis difficult as cloud cover can limit observations and repeated shots will not be located at the same position as before. The latter, however, would be interesting for the tracking of previously logged or disturbed forests. The combination of lidar measurements with other remote sensing techniques (e.g., radar data from the TanDEM-X satellite), thereby, seems promising to allow for an interpolation between lidar shots for large-scale analyses [68].

lidar measurements covering continents and combined with maps of forest cover (e.g., [3]) could enable the calculation of the forest structure and stem diameter distributions at large scale by using the suggested approach if the tree crown allometry is known. Various forest attributes can then also be derived from the estimated lidar-based stem diameter distributions (e.g., forest height, biomass, or basal area) and compared not only with local inventories but also with already existing information, for example, for the Amazon in terms of the basal area or forest biomass [69]. The approach proposed here can, thereby, contribute to the in-depth mapping and monitoring of forests and for supporting sustainable management, conservation, and the protection of forests.

5. Conclusions

We presented a novel approach for predicting the stem diameter distribution of forests based on lidar measurements. Our approach is profile-based and structure-based. This method requires only a vertical lidar profile (instead of point clouds) and integrates the available tree size allometries to infer the number of trees from the profile. By this, our approach complements previous methods based on statistical correlations between lidar and forest attributes or individual tree crown detection methods based on high-resolution lidar point clouds.

We demonstrated, in this study, a test of our approach with good accuracy for an old-growth tropical forest. Further tests of our method for other forest sites (at which tree allometries, lidar data and, for comparison, inventories are available) could help to comprehensively understand the impact of forest stand characteristics and plot size on the accuracy of the estimated stem diameter distributions from lidar.