Individual Tree Identification in ULS Point Clouds Using a Crown Width Mixed-Effects Model Based on NFI Data

Abstract

A linear mixed-effects model was used to relate crown width to height using an inventory plot as a random effect for trees in Czechia based on data from the National Forest Inventory (NFI). This model was used to estimate window size for a local maximum filter procedure (LMF) to detect individual tree tops in unmanned aerial laser scanning (ULS) point clouds of mixed species forest stands with diverse structures. Random model parameters were estimated for the study site based on several sample trees. Models calibrated with five or more samples achieved significantly better results (mean percentage error; MPE −0.17 for 5 samples) compared to when a fixed-effects model (MPE −0.62) was used. Lower performance was observed in dense stands with trees that were between 5 and 10 m in height. It was concluded that locally calibrated models predicting crown widths from tree heights might serve as a universal point of departure when searching for an optimal window size setting in LMF procedures.

1. Introduction

Due to mass production and increasing affordability, the use of unmanned aerial vehicles (UAV) for forest monitoring is no longer the prerogative of advanced scientific teams. A similar trend has been observed in the software used for processing the data acquired by UAVs into 3D point clouds representing forest trees during laser scanning (light detection and ranging-Lidar) operations or its surface in the case of photogrammetry. Still, subsequent steps are necessary for extracting relevant information from forest point clouds, which tends to be quite fragmented, and it often requires intensive fitting for differing situations.

The detection of individual tree tops from UAV-based photogrammetry point clouds is a frequently discussed procedure, and different methods have been reviewed by several authors, for example [1]. The local maximum filter (LMF) is one of the most commonly used methods. Several studies suggest that using LMF with a variable window size produces more promising results, especially in more diverse conditions [2]. Nevertheless, the correct size of the window, whether fixed or variable, is critical for the results, and an incorrect size can lead to serious errors. Determining the window size is often a trial-and-error process for particular stands or situations. The R-project [3] LidR package [4] allows for a variable LMF window size based on the chosen function of the Z-coordinate of the point clouds in which the tree tops are identified. The tree crown width might serve as a good parameter for setting the window size based on the tree height [2]. Using this method, a model is built using tree height to predict the crown diameter.

Several studies have concluded that the correlation between crown diameter and tree height is rather low; therefore, crown diameter prediction models are usually built using the tree stem diameter and several other independent variables [5,6,7]. Because these variables are not available in our particular case, a linear mixed-effects model with NFI plots treated as a random effect was used. The local calibration [8,9] of this model was expected to provide an ideal means of establishing an appropriate LMF variable window size for a given area of interest. The Czechia National Forest Inventory (NFI) [10] data provide a great opportunity to build such a model at the country level. Unmanned aerial laser scanning (ULS) point clouds together with the field survey results of tree counts on sample plots represent a promising means to assess the performance of local calibration based on different numbers of samples.

Although other individual tree detection methods that use, for example, info on the internal forest structure derived from ULS or airborne laser scanning (ALS) point clouds have started to emerge more recently [11,12], LMF represents one of the most common methods, especially for photogrammetry point clouds that are limited to the forest surface and that lack information on the internal structure. As stated above, using the correct window size in LMF applications is crucial for results. This study tries to confirm, whether local calibrated prediction of crown width might serve as a viable universal approach that can be used to determine the correct window size.

2. Materials and Methods

2.1. Study Site

Our study site, Skorkov (Lat: 50.2143906N, Lon: 14.7160739E; Position WGS 84), which is northeast of Prague, was established in a mixed forest stand comprising mostly flat terrain (188–190 m AMSL) and a species composition dominated by Pinus silvestris L. and Quercus robur L. Different age classes were present in 9.31 ha in rather rectangular arrangements, which is typical for Pinus management in this region.

2.2. Data Acquisition

The ULS point clouds were acquired using the UAV VUX-SYS setup (RIEGL Laser Measurement System GmbH, Horn, Austria) for UAV-borne data acquisition, which consisted of the UAV RiCOPTER, the VUX-1UAV laser scanner, and the AP-20 inertial measurement unit [13]. The laser scanner was set to its maximum pulse frequency of 550 kHz and had a registration of 200 scanning lines per second. The flight was performed at a constant altitude of 90 m aboveground at a constant ground speed of 6 m·s⁻¹. This setup reached an average point density of 200 points·m⁻² for each scanning line. Due to the overlap between the individual scanning lines, the final mean point density was around 600 points·m⁻². All UAV flights complied with the national requirements stated in [14] that were valid at the time of flight (September 2020), namely a maximum flight height below 300 m, a visual line of sight, airspace restrictions, and a minimal distance from people and buildings.

2.3. Laser Point Cloud Generation

The ULS point cloud was prepared from the data acquired by an UAV Ricopter equipped with a Riegl VUX-SYS1 sensor. The smoothed best trajectory estimates, which were based on the reference station post-processing kinematics with the virtual base station from Trimble VRS Now (Trimble Inc., Sunnyvale, CA, USA), were used and were prepared in the RiPROCESS environment using the RiPRECISSION module [13]. Point cloud classification (for ground and non-ground points) was conducted in the same software environment.

2.4. Sample Plot Surveys

To gather relevant data to evaluate individual tree identification using an LMF procedure, a field survey was conducted. A total of 39 circular plots were established Figure 1, with a variable plot radius related to the average tree height within the main canopy, which is shown in

Table 1. Sample plot size classes based on average tree heights of the main canopy.

Plot Class	Tree Height (m)	Plot Radius (m)	Plot Area (m2)	Minimal Register Limit
1	<3	3.00	28.27	1.5 m height
2	3–10	5.00	78.54	1.5 m height
3	>10	12.62	500.00	7 cm DBH *

* DBH—diameter at breast height (1.3 m height).

.

On each sample plot, all of the trees over the minimal register limit were counted. A randomly chosen tree was also sampled for its height and crown parameters. In the class 1 plots, maximum and minimum crown projection diameters were measured by tape, and crown width was calculated as an average of these two measurements. For the class 2 and 3 plots, a procedure similar to the one used during NFI surveys [10] was used, including measuring the position of the ground projection of at least five points defining the crown extent. The field survey was conducted using the Field-Map software/hardware set (IFER—Monitoring and Mapping Solutions Ltd., Jílové u Prahy, Czechia), which combines a flexible GIS software Field-Map with electronic equipment (Mapstar compass and ForestPro range-finder with inclinometer) for mapping and dendrometric measurements.

2.5. Crown Width Model Data

The NFI data acquired by the Forest Management Institute since 2001 were used to build linear mixed-effects models for the tree crown width predictions of individual trees based on height. The intended use of the model within the LMF procedure was limited to a single predictor variable, height; thus, other frequently employed variables [5,6,7], such as diameter at breast height (DBH), were not employed to construct the crown prediction models. The NFI dataset included data on more than ninety-four thousand trees measured on more than twenty-two thousand plots, as seen in

Table 2. Selected summary of NFI data.

Characteristic	n 1 = 94,066
Plots (n)	22,532
Mean height (m)	25 (4 to 54)
Mean crown width (m)	5.57 (1.35 to 19.33)
Species group
Abies alba Mill.	1202 (1.3%)
Acer spp.	1709 (1.8%)
Alnus spp.	2177 (2.3%)
Betula spp.	2382 (2.5%)
Carpinus betulus L.	890 (0.9%)
Fagus sylvatica L.	6475 (6.9%)
Fraxinus spp.	1340 (1.4%)
Larix decidua Mill.	4300 (4.6%)
Other broadleaved hard spp.	970 (1.0%)
Other broadleaved soft spp.	2672 (2.8%)
Other coniferous spp.	763 (0.8%)
Picea abies (L.) Karst.	48,767 (52%)
Pinus sylvestris L.	13,691 (15%)
Quercus spp.	6728 (7.2%)

¹ Mean (Range); n (%).

. Only trees from the upper tree layer, which was determined as defined by the International Union of Forest Research Organizations and according to observed crown projection area and height measurements, were selected. The NFI field surveys [10] identified tree crowns in the form of a horizontal projection area. For the purposes of this study, the crown width was calculated as a function of crown area by assuming a circular shape centered on the tree stem.

2.6. Model

A simple linear mixed-effects model was built iteratively with different predictor or predicted variable transformations to achieve optimal results. The best performing model was:

$$\ln \left(C{W}_{ki}\right)={\beta }_0+{\beta }_1{H}_{ki}/100+{\alpha }_{0k}+{\alpha }_{1k }{H}_{ki}/100+e_{ki}$$

(1)

where CW_ki is the crown width of tree i in the stand/plot k (m), H_ki is the height of tree i in the stand/plot k (m), β₀ and β₁ are fixed population parameters to be estimated, α_0k and α_1k are the random parameters to be estimated with zero expectations for stand k, and e_ki is the random residual error for tree i on stand/plot k. The H_ki variable had to be scaled, or else the model would not converge. This model is basically the same as the one described by Lappi [8], which used height-DBH data.

The model parameters were first estimated using all of the trees from the NFI dataset with crown width and height measurements (n = 94,066), regardless of species, and the NFI plots were treated as a random parameter. This model was called f1. The model was developed in this manner in an effort to account for parameters other than height, as represented all together by individual inventory plot.

Based on the field surveys of the study site, a limited number of tree species are present, with three species representing the majority—P. sylvestris, Q. robur, and Quercus rubra L. Therefore, it was decided to create a subset of the NFI dataset to limit it to only the P. sylvestris and Quercus species, and this subset was then used to construct model f2.

2.7. Model Calibration

The use of a linear mixed-effects model allows for calibration within particular stands. Based on the sample data of crown widths and tree heights obtained at the stand of interest, the values of random parameters from equation Equation (1) can be estimated [9] according to Equation (2).

$$\hat{a}=(Z^{\prime }{\hat{R}}^{-1} Z+{\hat{D}}^{-1}{)}^{-1}{Z}^{\prime }{\hat{R}}^{-1}\left(y-\hat{\mu }\right)$$

(2)

where $\hat{a}$ is the matrix of estimated random parameters; Z is the design matrix associated with random parameters; $\hat{D}$ is an estimate of $D=var\left[\begin{array}{c}{\propto }_{0k}\\ {\propto }_{1k}\end{array}\right]=\left[\begin{array}{cc}{\tau }_{0k}& {\rho }_{01k}\\ {\rho }_{01k}& {\tau }_{1k}\end{array}\right]$ the variance–covariance matrix for the random parameters; $\hat{R}={\hat{\sigma }}_{e_{ki}}^2\times I$ is the estimated variance–covariance matrix for residual errors of individual trees where ${\hat{\sigma }}_{e_{ki}}^2$ is the square of the residual standard error; $\hat{\mu }$ is a vector of predicted values from the fixed effect model; and y is a vector of observations of the dependent variable from the particular stand.

To evaluate this calibration method on our model, one random plot from the NFI dataset was iteratively taken out as a calibration plot. The LME model (Equation (1)) was fitted to the NFI dataset without the calibration plot. Calibration was performed repeatedly based on the height and crown width data of 1, 3, 5, and 10 random trees from the calibration plot according to Equation (2).

The crown width was then estimated based on both the marginal model (without random parameters) and the calibrated model with estimated random parameters. Both estimates were compared to field-measured crown width values from the calibration plot using the normalized root mean squared error (nRMSE; Equation (3)). This procedure was repeated one hundred times to assess the performance of the calibrated models.

$$nRMSE=\frac{\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\hat{y}}_i-y_i\right)}^2}{n}}}{\overline{y}}$$

(3)

where n is the number of observations, y is the observed value, $\overline{y}$ is mean of observed values, and $\hat{y}$ is the estimated value.

2.8. Study Site Calibration

For calibration of global models explained in the previous section, data that were obtained for the crown widths and heights in sample plots were used. A total of 39 trees were measured for this calibration. To evaluate the influence of the different number of calibration measurements both models were calibrated with randomly chosen 3, 5, 10, and all 39 calibration samples. Together with the uncalibrated marginal models, a total of 10 formulas were prepared for use during the LMF procedure. Model function names are identified by the model used (f1; f2) and the number of sample trees used to calibrate the model (e.g., cal3 = 3 sample trees used); for example, Model f2 calibrated with 10 sample trees is identified as f2_cal10. The details of both models, which are necessary for local calibration, are presented in Appendix A.

2.9. Trees Identification

The ULS point cloud was processed through the following pipeline performed using the R statistical package [3] with the LidR [4] and ForestTools [15] packages:

• Read ULS point cloud data.
• Perform normalization of point cloud with ground points classified by the cloth simulation filter [16].

Then, for each of the 10 model functions:

3.

Run the LMF procedure to identify the local maximum within a circular window. The diameter of the window is defined by the crown width predicted according to model function using Z coordinate as the height variable.

4.

Intersect the resulting tree tops with the sample plots.

5.

Compare the number of trees identified by the LMF procedure and the number of trees from our field survey. Comparisons were made based on the percentage error (PE) determined using Equation (4) and the mean percentage error (MPE) determined using Equation (5).

$$PE=\frac{{\hat{y}}_i-y_i}{y_i}$$

(4)

$$MPE=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\frac{{\hat{y}}_i-y_i}{y_i}$$

(5)

where n is the number of sample plots, y is the observed number of trees at sample plot i, and $\hat{y}$ is the estimated number of trees at sample plot i.

An overview of all of the processes described above part by part, leading to final results, is provided in Figure 2. The model calibration test is not included in this flow chart.

All of the analyses were conducted in R-Studio version 1.4.1103 [17] with R version 4.0.3 [3] and using the following libraries: LidR [4], ForestTools [15], sjPlot [18], JTools [19], Readr [20], GGplot2 [21], and gtsummary [22].

3. Results

3.1. Models

The results of models f1 and f2 are visible in

Table 3. Linear mixed-effects models f1 and f2 results.

	Model f1 ln(CW)			Model f2 ln(CW)
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.9692	0.9552–0.9831	<0.001	1.0471	1.0140–1.0802	<0.001
H/100	2.9192	2.8671–2.9714	<0.001	2.7494	2.6134–2.8854	<0.001
Random Effects
σ2	0.0579			0.0588
τ0k	0.3248			0.3792
τ1k	3.0351			4.7034
ρ01k	−0.8865			−0.8578
ICC	0.6090			0.6663
N k	22,532			7009
Observations	94,066			20,419
Marginal R2/Conditional R2	0.226/0.697			0.129/0.709

CW—crown width. CI—confidence interval. H—tree height. σ², τ_0k, τ_1k, ρ_01k explained in Equation (2). ICC—Intra-class correlation. N_k—number of plots. Marginal and conditional R-squared statistics, based on [23].

as well as in Appendix A.

3.2. Model Calibration Test

Using the local calibration procedure, a single plot from the NFI dataset (randomly selected from plots with more than 15 trees measured) was used as a hypothetical new plot. The model was calibrated repeatedly so that the model could be fit to the rest of the data set using 1, 3, 5, or 10 trees from the calibration plot. For each of these calibrated models, the crown width was estimated and compared to the original values from the NFI dataset Figure 3.

This procedure was repeated one hundred times. The overall performance of the calibration process is evident in Figure 4, which shows that the model has a better fit when calibrated with more sample trees, although differences when using 3, 5, and 10 sample trees were relatively small. Using this procedure, it is evident that the fixed-effects model could be improved by the calibration process.

A Kruskal–Wallis test was conducted on the nRMSE values based on the number of samples used for calibration and demonstrated a significant difference between groups (p-value = 1.219 × 10⁻¹³). Subsequently, Dunn’s test confirmed the optimal use of at least three samples for calibration; further increases in the sample numbers did not contribute any significant improvements, as seen in

Table 4. Dunn’s test p-values for nRMSE grouped by number of samples used for calibration.

Samples (n)	0	1	3	5
1	0.001
3	0.000	0.033
5	0.000	0.003	0.443
10	0.000	0.001	0.219	0.608

p-values lower than 0.05 are bold.

.

3.3. Study Site Model Calibration

A similar procedure was used for the local calibration of model f1 and f2 based on the data obtained from the sample plots at the study site. Again, several versions with 3, 5, 10, and the complete set of field-measured tree samples (39) were assessed. Together with the uncalibrated fixed-effects model, 10 formulas were produced overall to model crown widths based on tree heights, as seen in Figure 5 and Appendix A.

3.4. Tree Identification

The MPE between the field survey-observed and LMF-identified number of trees on the sample plots was −0.23 on average, varying from −0.62 to 0.04 (standard deviation of 0.21). The nRMSE had values between 0.50 and 0.78. The MPE distribution determined by the function used for the LMF procedure is described in Figure 6a, which indicates better results for calibrated functions overall, with the surprising exception of the fixed-effects only model for Pinus and Quercus spp., f2mar.

To better asses the influence of tree height, sample plots were divided into four consecutive categories according to max height of the trees observed in the plots. Figure 6b and

Table 5. Mean percentage error between LMF-identified and field survey-observed number of trees on sample plots.

Function	nRMSE	Tree Count MPE by Height Category
		All	0 to 5 m	5 to 10 m	10 to 20 m	20 to 30 m
f1mar	0.78	−0.62	−0.77	−0.77	−0.68	−0.40
f1_3cal	0.68	−0.48	−0.60	−0.67	−0.57	−0.23
f1_5cal	0.53	−0.17	−0.19	−0.49	−0.21	0.05
f1_10cal	0.53	−0.18	−0.17	−0.49	−0.23	0.02
f1_39cal	0.51	−0.08	0.02	−0.42	−0.15	0.09
f2mar	0.61	0.03	−0.38	−0.52	−0.06	0.69
f2_3cal	0.61	−0.39	−0.42	−0.58	−0.45	−0.22
f2_5cal	0.54	−0.23	−0.17	−0.49	−0.29	−0.09
f2_10cal	0.52	−0.09	−0.00	−0.43	−0.17	0.09
f2_39cal	0.50	−0.08	0.07	−0.41	−0.17	0.07

reveal that the unexpectedly good performance demonstrated by the f2mar function is the result of overestimating the number of the tallest trees on one hand and underestimating the number of the smallest trees on the other. Table 5 also displays the increased accuracy of the calibrated models—with 5- and 10-sample calibrations showing comparable performance to the full 39-sample calibration. According to nRMSE, the best performing model, f2_39cal, is closely followed by all models when 5 or more samples are used for calibration. In terms of height categories, the most underestimated category was the number of trees that were 5 to 10 m in height.

A Kruskal–Wallis test suggested significant differences in the tree count error determined by the model functions used for LMF (p-value = 9.302 × 10⁻⁶). A subsequent Dunn’s test revealed a significant improvement in the results when at least five samples were used for calibration (

Table 6. Dunn’s test p-values for tree count error grouped by function used for LMF.

	F1mar	f1_3cal	f1_5cal	f1_10cal	f1_39cal	f2mar	f2_3cal	f2_5cal	f2_10cal
f1_3cal	0.439
f1_5cal	0.006	0.071
f1_10cal	0.005	0.067	0.998
f1_39cal	0.001	0.014	0.665	0.683
f2mar	0.001	0.009	0.537	0.554	0.905
f2_3cal	0.175	0.706	0.228	0.218	0.063	0.037
f2_5cal	0.013	0.167	0.810	0.791	0.451	0.347	0.445
f2_10cal	0.001	0.018	0.692	0.710	0.961	0.837	0.068	0.511
f2_39cal	0.001	0.014	0.646	0.665	0.966	0.942	0.056	0.435	0.925

p-values lower than 0.05 are bold.

). Using model f2, which was only fitted to data from the tree species present in the study site, did not bring expected improvements to any of the models.

4. Discussion

Iterative local calibration tests of a global LME model revealed significant improvements in model performance, even when only a single calibration sample was used. Three samples were interpreted to be the optimal amount because further increases in the number of samples did not create significant improvements. Several studies have similarly recommended only a limited number of samples for the calibration of mixed-effect height–diameter models [24,25,26].

The local maximum filter belongs to one of the most commonly used approaches for individual tree top detection in point clouds [1]. Determining the appropriate window size for LMF is often a trial-and-error process for particular stands or situations. The presented study tries to improve this situation by using a tree crown width that is to set the LMF window size.

Most reported approaches use other predictor variables, most often DBH, for crown projection rather than height [5,27,28] and thus disqualifies using those models for local maximum filter techniques to detect individual tree tops.

Our suggested approach using a linear mixed-effects model to predict the tree crown width based on tree height using a few samples (results at our study site suggested at least five) that were obtained at the study site for local calibration might serve as a viable option. Our presented model reached a conditional R-square value of 0.70, which seems more promising than that obtained for a similar model built on a much smaller data sample, which had an R-square value of 0.51 [2].

With the proposed LMF individual tree-detection method based on the crown width-determined window size using ULS point cloud data, an nRMSE of around 0.50 was achieved. The MPE metric, ranging between −0.62 and +0.03, revealed that in most cases, our approach led to a lower number of trees identified compared to the field survey. This can be attributed to the problematic distinction of single trees in Pinus and Quercus forest types, in our study site with a diverse tree structure. Not all of the trees have grown into the canopy layer, and the distinction of those trees is rather difficult, especially with broadleaved species. This is also supported by the poor performance of trees between 5 and 10 m. Somewhat similar tree detection results ranging from 50% to 140% using different ITD on ULS point cloud methods were reported by Wang et al. [12]. A study by Grznárová et al. [29] produced a detection rate of 95% for coniferous forests and a rate of 71% in broadleaved forests. Nevalainen et al. [30] achieved tree identification rates between 64 and 97%. Results with 38 to 85% of trees undetected across plots were also reported by Jeronimo et al. [31], who also reported that smaller trees were unable to be detected with the same level of success. Throughout the cited studies, it is clearly visible that ITD precision decreases with increased forest structure complexity.

5. Conclusions

The idea of using crown width to set the LMF window size seems to be logically justified. Still, correctly estimating the crown width based on tree height alone poses a significant challenge, as other influential factors, such as tree species, genetic variability, DBH, stand and site characteristics, etc. cannot be directly addressed in LMF. The locally calibrated mixed-effects model built on an extensive NFI dataset represents an approach that strives to cover all these other factors using the plot itself as a random effect.

Our individual tree top detection results only reached limited accuracy in terms of the number of trees identified in the sample plots. On the other hand, these results are comparable to those obtained in other research and suggest that the presented approach using a local maximum filter with a variable window size governed by locally calibrated models predicting crown width from tree height might serve as a universal point of departure when searching for optimal window size settings in LMF procedures and can provide reasonable accuracy, even in more complex forest structure.