Delineating and Reconstructing 3D Forest Fuel Components and Volumes with Terrestrial Laser Scanning

Abstract

Predictive accuracy in wildland fire behavior is contingent on a thorough understanding of the 3D fuel distribution. However, this task is complicated by the complex nature of fuel forms and the associated constraints in sampling and quantification. In this study, twelve terrestrial laser scanning (TLS) plot scans were sampled within the mountain pine beetle-impacted forests of Jasper National Park, Canada. The TLS point clouds were delineated into eight classes, namely individual-tree stems, branches, foliage, downed woody logs, sapling stems, below-canopy branches, grass layer, and ground-surface points using a transformer-based deep learning classifier. The fine-scale 3D architecture of trees and branches was reconstructed using a quantitative structural model (QSM) based on the multi-class components from the previous step, with volume attributes extracted and analyzed at the branch, tree, and plot levels. The classification accuracy was evaluated by partially validating the results through field measurements of tree height, diameter-at-breast height (DBH), and live crown base height (LCBH). The extraction and reconstruction of 3D wood components enable advanced fuel characterization with high heterogeneity. The existence of ladder trees was found to increase the vertical overlap of volumes between tree branches and below-canopy branches from 8.4% to 10.8%.

1. Introduction

An increasing frequency of wildland fires has been reported across North America [1,2,3] since the 1980s [4] and is projected to continue to increase in coming decades, mainly due to changes in climate (warmer, drier conditions) and anthropogenic drivers [5,6]. The threat of high-intensity fires raises concerns due to reduced environmental resilience and limited fire control resources [7,8]. In recent years, some wildland fire research has been dedicated to the optimization of large-scale wildfire systems, for example, the Canadian Forest Fire Danger Rating System (CFFDRS) [9] and the National Fire Danger Rating System (NFDRS) [10]. Among the system components, climate and weather have been extensively studied and are broadly well-established as general indicators of fire spread [11,12,13,14]. This is primarily because weather variables account for the majority of fire intensity and spread when compared to other factors [15].

The system inclusion of wildland fuel, another fundamental predictor of fire spread and intensity in fire behavior models, is generalized within a geospatial data layer. For example, the fuel inputs for the Canadian Fire Behavior Prediction (FBP) system includes inventory variables such as fuel type, fuel load, and moisture content within relatively broad areas or fuel types [16] but are not indicative of the spatial heterogeneity of vegetation fuels. Despite the general use of fuel types, empirically-based models cannot use spatially varying fuel structures. Physically-based simulation models enable more complex spatiotemporal factorization of fuels such as fuel breaks and shrub/ladder fuel layer mixtures including: (1) 2D raster-based Prometheus [17], FARSITE [18], Burn-P3 [19], and Cell2Fire [20], (2) 3D terrain-based FIRETEC [21] and QUIC-Fire [22], and (3) 3D tree-level WFDS [23,24], FastFuel [25,26], and STANDFIRE [26]. The support for spatially-explicit fuel inputs in modern fire behavior simulators may, in part, be attributed to the availability of new observation technologies represented by multi-spectral and thermal remote sensing systems [27,28].

The increasing availability of LiDAR data can further enhance the opportunity to create high-resolution 3D vegetation and fuel units for more realistic fire prediction and simulation [25] and for a better understanding of the distribution of fuels in classified fuel types [29]. Terrestrial laser scanning (TLS) is particularly useful for measuring the 3D structure of trees, which can be classified into fuel distributions from the ground level. The high density of returns also improves the identification and quantification of understory and ladder fuels, which may be occluded or undetectable using reduced point density airborne lidar systems or from overstory shadowing, which impacts optical remote sensing technologies [30,31,32]. TLS enables fine-scale measurements of micro-scale topographic variations, vegetation dimensions, species, and the ability to quantify live and dead fuel components [33,34]. From these measurements, 3D branch and crown reconstruction, spatial distribution and the connectivity of fuels, and any unexplored fuel-related metrics can be quantified from very high spatial resolution 3D point clouds [34]. For example, García et al. [35] used TLS to extract fuel attributes including canopy cover, canopy base height (CBH), and fuel strata gaps, while Alonso-Rego et al. [36] further extracted canopy fuel load (CFL) and canopy bulk density (CBD) for use in fire behavior models. Observing and modeling such fine-resolution 3D fuels from TLS are useful for quantifying 3D fuel variability within and across fuel types [37,38], which would normally be manually measured in fuel plots, identifying changing fuel conditions in climate mediated or otherwise disturbed or changing vegetation [39,40], and for use in future wildfire behavior simulations [41].

A more in-depth integration between the TLS data and fire simulation models has potential at the point or voxel level [42], but has rarely been realized in the literature. Fire simulators typically rely on artificial 3D scenes with tree geometries interpolated from remotely sensed data or inventory attributes [43]. FUEL3D is one of the 3D scene generators [44,45] that does not yet incorporate true forest heterogeneity such as litter, shrub, and ladder fuels, which may be observed in TLS point clouds [46]. An obvious barrier of applying TLS at the plot scale is the challenge of transferring point data into volumetric fuel components for fire simulation, which requires: (1) the classification of fuel types, and (2) a reduction in the potential for point cloud occlusion.

The former step differentiates fine/coarse fuel types of different flammability and fire spread patterns. For example, Chen et al. [47] classified forest fuels from TLS into strata of overstory, elevated, near-surface, and surface using a height slicing method. Loudermilk et al. [37] first used TLS to estimate the fuelbed structure at the submeter scale, finding that simple partitioning of plots into point cloud groups based on k-means significantly improved the fire variable predictions (e.g., temperature) compared to individual plot models. García et al. [35] reported the first quantitative extraction of fuel-oriented variables including the canopy height, canopy cover, CBH, and fuel strata gap to facilitate the assessment of vertical fire spread risk. Similarly, Rowell et al. [42] counted fuel surface voxels and fuelbed porosity to estimate the fuel volume, discovering a strong linear relationship between the TLS-based volume and fuel mass that was higher than the field measurements, particularly within the lowest 10–20 cm stratum with the highest fuel concentration. They also suggested that fuel biomass estimation could be further improved using more sensitive TLS metrics than overall voxel volume occupation. Rowell et al. [43] emphasized the importance of partitioning fine-scale fuelbeds for fire behavior models because fuels combust differently based on changes in relative humidity and ambient temperature. Yrttimaa et al. [33] demonstrated the potential of TLS for mapping coarse woody debris, detecting 68% volume using an automated method, while Muir et al. [48] found it challenging to accurately determine the downed woody debris length from TLS due to grass and tree trunk occlusion. Without calibration, TLS alone might not be the superior tool for estimating the downed woody debris volume. Yrttimaa et al. [33] noted that a maximum of 83% of volume was detectable from TLS with visual interpretation. Extracting ladder fuels from TLS has also been demonstrated as a viable approach. Ladder fuel refers to vegetation such as low-hanging branches, shrubs, and leaning trees that provide a “ladder” for ground fires to climb up into the crown. The presence of ladder fuel increases the risk of more intense, fast-spreading crown fires. Forbes et al. [49] concluded that the use of TLS to quantify ladder fuel density variation across strata was the best predictive model of burn severity among multiple remote sensing approaches. TLS has also been applied to fuel change analysis, with the TLS voxel density significantly relating to shrub fuel bulk density from both pre- and post-fire measurements [49]. One notable constraint of TLS is the limited spectral information from laser returns, meaning that most TLS use is based on geometric information. However, an experiment using two TLS devices with different wavelengths demonstrated significant accuracy (r² = 0.97) in explaining the fuel moisture content using pulse intensity and its standard deviation in a multiple linear regression model [50].

The latter step corrects the bias of occluded wood volume, typically through wood geometrical reconstruction using a quantitative structure model (QSM) or fractal theory [51,52,53,54,55]. To achieve plot-level wood reconstruction, individual trees also need to be isolated and wood/foliage segmented [56]. Such point cloud algorithm developments are spread across separate studies on TLS [57,58,59,60,61,62,63] and are further extended to more portable options like mobile laser scanning (MLS). MLS can either be handheld or backpack-mounted with comparable accuracy [64,65], or vehicle-mounted for expanded coverage with minor accuracy loss [66]. However, these methodological developments have yet to be synthesized in the context of fuel analysis. This gap may be attributed to the relatively recent emergence of many of these algorithms in need of comprehensive benchmarking, validation, and software development at the broader level. Among the benchmarked algorithms, the highest accuracies range from ~90% for tree/ground classification [67], ~80% for individual tree isolation [68,69], ~90% for wood/leaf classification [70,71,72], to >98% for QSM wood volume estimates [73]. However, these have not been investigated for fuel extraction.

The collection of fuel plot data is time consuming, expensive, and does not sample all attributes of the structure of trees. TLS offers an opportunity for the rapid collection and processing of data, which represents many of the 3D structural attributes that are important for characterizing the fuel distribution, change, and their potential use for parameterizing fire models. This study demonstrates a 3D fuel analysis framework for TLS with novel processing examples of wood point extraction, classification, isolation, and reconstruction applied to a variety of pine/fir forests and with mountain pine beetle disturbance at Jasper National Park, Alberta, Canada. The crown, ladder, and surface fuel components were characterized at a high-resolution, 3D scale and their distributions were investigated at the tree- and plot-level, respectively. This study demonstrates the promising application of TLS for detailed fuel–fire assessment and prediction. The major objectives of this study are:

• Delineation of 3D plot components with deep learning classifiers amenable for multi-layer fuel analysis;
• Provision of wood volume distribution at the branch, individual tree, and plot level from the delineated 3D components based on the QSM and scaling methods.

2. Study Area and Data Collection

Our study site is located near the tributary basins of the Athabasca River and Miette River within the eastern main ranges of the Canadian Rockies in Jasper National Park (Figure 1). The elevation of the valley basins ranges between 1050 and 1200 m with a typical annual average temperature of 4.1 °C. The predominant species in the park (and examined in the TLS and field plots) is Pinus contorta (lodgepole pine), categorized as fuel type C-3 in the Fire Behavior Prediction (FBP) system. A decade long attack of mountain pine beetle (MPB) has altered the existing fire regime by dehydrating and often killing many of the pine trees, causing the transition from healthy green trees into the red-phase, where needles have died but remain on branches (and are red in color), and grey-phase, where needles have dropped from trees, leaving remaining grey branches on the trees and causing the accumulation of dry woody debris on the forest floor. This region is characterized as a lightning-shadow area with a low flash density of less than 0.5 per km² per year [74]. The majority of lightning events occur during the night-to-morning period near water bodies, which is atypical of severe lightning commonly associated with local solar heating for thunderstorms [75]. However, studies have indicated that the eastern aspect of terrain features within the elevation range between 1200 and 2000 m could enhance flash density [75,76], consequently increasing the risk of fire occurrence and subsequent spread into adjacent valley basins. The accumulation of fire fuels from mountain pine beetle and drying potential is also of significant concern for the town of Jasper, which is also located at the confluence of the basins (Figure 1). This was exemplified in September 2022 when a fire occurred 15 km north of the town of Jasper, reinforcing the importance of fuel monitoring, prescribed fires, and fire prediction in the area.

Twelve adjacent TLS plots coincident with the field data were collected at six locations along the valley south and west of the town in July 2021. Two plot centers were located side by side at a distance of 40 m to reduce movement of the TLS and to increase sampling efficiency at each location (Figure 1). A Teledyne Optech Polaris TLS station (1550 nm) was placed at each plot center for a 360° × 120° scan and at four corners 20 m from the plot center to complement the central scan and reduce the occlusion of returns. TLS scans were collected at a rate of two plots per day. Trees within an 11.3 m radius from each plot center were also measured for their species, tree height, diameter-at-breast-height (DBH), live crown base height (LCBH), dead crown base height (DCBH), tree phases under MPB attack, and the approximate degree of dead condition (Dead %). A tree with a Dead % value of 100 characterized a leaf-off tree. A total of 572 trees were surveyed including six major species: lodgepole pine (Pinus contorta), white spruce (Picea glauca), black spruce (Picea mariana), trembling aspen (Populus tremuloides), Douglas fir (Pseudotsuga menziesii), and subalpine fir (Abies lasiocarpa). The selection of twelve plots represented a wide range of fuel conditions with varying stem density, terrain slopes, sub-canopy thickness, and tree health (

Table 1. Summary characteristics of trees within each of the twelve plots (named ‘JP’ for the Jasper plot). DBH, Height, LCBH, and DCBH plot-level averages from individual tree mensuration within 11.3 m from the plot center are described. Subcanopy height and slope were averaged over a grid at 1 m resolution from the georeferenced TLS plot scans.

Plot	Coordinates	Stem Density (ha−1)	Subcanopy Height (std) (m)	Slope (std) (°)	DBH (std) (cm)	Height (std) (m)	LCBH (std) (m)	DCBH (std) (m)
JP1	52.669°, −117.888°	801	1.95 (1.44)	12.1 (7.3)	10.2 (4.6)	8.2 (3.4)	2.9 (2.8)	3.2 (1.5)
JP2	52.669°, −117.888°	1626	2.61 (1.94)	12.0 (9.1)	11.1 (5.2)	9.7 (3.7)	4.1 (3.0)	3.7 (1.5)
JP3	52.669°, −117.888°	776	1.40 (1.14)	3.0 (2.1)	13.7 (4.4)	11.1 (3.5)	5.1 (1.7)	2.3 (0.8)
JP4	52.800°, −118.015°	1426	1.37 (1.02)	3.0 (3.5)	12.7 (5.1)	10.6 (4.0)	5.2 (2.6)	2.3 (0.5)
JP5	52.800°, −118.015°	1301	3.24 (2.50)	9.2 (10.0)	11.9 (7.2)	10.4 (5.8)	2.4 (3.3)	5.7 (2.4)
JP6	52.670°, −117.891°	2051	3.21 (2.28)	11.3 (13.6)	10.1 (7.0)	9.0 (5.8)	2.1 (3.7)	4.0 (2.7)
JP7	52.789°, −118.005°	450	1.19 (0.81)	3.5 (3.1)	16.8 (9.2)	12.7 (6.3)	4.5 (4.3)	2.6 (1.0)
JP8	52.788°, −118.005°	300	1.39 (1.17)	2.5 (2.9)	14.0 (8.6)	9.3 (5.3)	0.8 (0.9)	1.6 (0.4)
JP9	52.788°, −118.005°	1876	2.77 (2.09)	23.4 (12.1)	15.2 (5.8)	15.5 (5.8)	6.8 (4.8)	7.1 (3.0)
JP10	52.870°, −118.253°	1501	2.34 (1.69)	20.8 (7.4)	16.1 (5.5)	16.4 (5.6)	9.4 (4.1)	6.8 (3.3)
JP11	52.724°, −117.929°	1101	2.56 (2.43)	7.2 (7.3)	16.1 (5.8)	14.2 (5.4)	6.0 (4.9)	7.5 (2.2)
JP12	52.725°, −117.930°	1101	2.41 (2.17)	6.9 (8.4)	18.2 (6.1)	16.9 (5.9)	1.5 (2.2)	5.8 (3.5)

). Individual trees were not geolocated due to time constraints. Instead, plot quadrants were separated by four cardinal directions and tree locations within quadrants were visually identified. Figure 1 shows an example ground-level view of the JP9 plot in which ladder fuels and downed woody logs are prevalent in the plot.

For each scan, a point cloud was generated using all laser returns, incorporating an intensity attribute for each point. Point clouds of corner scans from TLS were co-registered manually based on the 3D tie points with CloudCompare (v2.11) [77]. The selection of tie points are iterative trial and error based on the sharp branch edges until no obvious visual displacement is found. A plot consists of 7–10 tie points spread over different corners. The co-registration error was represented by the root-mean-square error of all tie points in the plot from CloudCompare. The registration errors per plot are summarized in Figure 2 with an average error of 1.9 cm. It is worth noting that in situ targets are commonly preferred for achieving optimal registration accuracy. These were excluded from this study due to the inefficiency associated with their setup and their limited discernibility in plots rich with understory vegetation. All registered scans were merged into a single ‘plot’ point cloud, cropped to a circular radius of 25 m, and then resampled to a spatial resolution of 5 mm for further processing. A spatial resolution of 5 mm was specifically selected as it provides a reasonable balance between the ability to discern branches from noise conifer trees and computational efficiency for subsequent analyses. TLS plot scans were georeferenced to point clouds acquired from airborne laser scanning (ALS) collected the same summer (early August, 2021). The georeferencing followed a similar manual registration procedure as that employed for TLS co-registration. Iterative selection of 10–14 tie points was performed on the tree tops and crown edges from both the TLS and ALS datasets until visual overlap of the crown shapes was achieved. The mean error associated with this TLS-to-ALS georeferencing process was determined to be 24.8 cm.

3. Methodology

3.1. Workflow and Preprocessing

The primary workflow decomposed TLS point clouds into eight classes, which are important for quantifying fuels: tree stems, tree branches, tree foliage, seedling stems, downed woody logs, below-canopy branches, ground points, and other non-wood foliage and litters (Figure 3). Within this classification scheme, below-canopy branches refer to non-stem wood points with diameters greater than 5 cm within the saplings, shrubs, and surface layers. Ladder fuels are further differentiated into categories such as leaning trees, below-canopy dead branches, and downed woody logs. Within the scope of this paper, these leaning trees and stems are respectively referred to as “ladder trees” and “ladder stems”. It is worth noting that some dead fallen trees may not make contact with the ground but could possess a moderate leaning angle. From our observations, such instances are relatively infrequent. We identified downed woody debris based on a high inclination angle exceeding 60° as this criterion clearly demarcates the boundary between the leaning tree layer and the downed woody debris layer. Varying classifiers and filters were used to extract the class of interest for each point and the eight classes were combined and used to form a multi-class point cloud per plot. The classifications were programmed using PyTorch (v1.13), Python, and all the remaining processes using MATLAB (R2022b). PC configuration included an Intel^® Core™ i7-12700K 12 × 3.60 GHz, 128 GB RAM, and NVIDIA GeForce RTX 3090 (24 GB).

The subsequent post-classification scheme is illustrated in Figure 4. Individual tree components such as stem, branch, and foliage were clustered. Wood was reconstructed into cylinder-shaped 3D architecture models using a custom QSM model. All other coarse wood components were isolated and reconstructed to estimate the wood volume across each plot.

3.2. Classification of Tree Point and Coarse Wood Components

The major components, namely, tree, stem, and downed woody log points were extracted using a state-of-the-art deep learning classifier SegFormer [78] customized in a 3D voxel version (Figure 5a). The definition of trees in this study did not include saplings with a height < 5 m. SegFormer replaces the classic convolutional neural network (CNN) layers with self-attention-based transformers, surpassing well-recognized semantic classifiers such as DeepLab V3 [79] and Swing Transformer [80] in terms of an overall rating of computation efficiency and accuracy. The SegFormer parameter values remained set to default and its configuration of voxelization, sampling, training, and testing followed the scheme in Xi et al. [71] except for the block size and voxel resolution (

Table 2. Customization of SegFormer for the different wood component classes of interest.

Classifier	Block Size	Voxel Resolution	Attributes	Sample Size (Training)	Sample Size (Testing)
Tree	128 × 128 × 128	10 × 10 × 10 cm3	Binary (0, 1)	4 plots	1 plot
Stem	128 × 128 × 128	4 × 4 × 4 cm3	Binary (0, 1) + Intensity	10 clips (~8 × 8 m2)	2 clips
Downed logs	128 × 128 × 128	2 × 2 × 2 cm3	Binary (0, 1) + Intensity	21 clips (~8 × 8 m2)	4 clips

). A block of voxels was the input unit to the SegFormer. The block size was selected as 128 × 128 × 128, which maximally exploited the GPU RAM capacity. The voxel resolution was customized based on the scale of interest: for example, smaller features such as downed woody logs matched a finer classification resolution (Table 2). The laser intensity was appended to each voxel as an additional layer to enhance the detection of coarse wood (stems and downed woody logs). Four plots with manually delineated tree and non-tree classes were used as the training sample (Figure 5b) with another single plot used as the independent testing sample. It was manually intensive to label wood and non-wood points as the training and testing reference over entire plots. Therefore, a subset of the plot scan with a size of ~8 × 8 m² was clipped. A total of 10 and 21 clips (or subsets) were created from the TLS scans to train the classifiers for the stems (Figure 5c) and downed woody logs and branches (Figure 5d), respectively, along with two and three clips for testing. A training block was randomly located within each random clip and was also rotated horizontally by a random degree. This was considered as a data augmentation process to mitigate overfitting under a small sample size. The accuracy evaluation metric was the mean inter intersection over union accuracy (mIoU). Incorrectly classified points were edited and cleaned for the wood isolation and tree clustering processes.

3.3. Individual Wood Isolation

The classified stem points inevitably contained small isolated fragments detectable using the connected component analysis [82]. Any connected fragment with a length or width < 5 cm was considered as noise and discarded. The remaining points were then clustered into small segments using an l0 cut-pursuit algorithm [83]. The cut-pursuit algorithm has two main parameters: the number of neighbors (K) from the nearest neighbor search as the minimal unit, and the regularizer λ as the overall graph scale [69]. The greater the K and λ set, the coarser the clusters that would be generated. Both were set to 20 in this study, and the resulting clusters were mainly sized between 0.1 and 10 m. Due to stem occlusions in the crown area, these clusters did not encompass the complete stem components and required merging.

Clusters were sectioned at 0.1 m intervals along their principal directions, as determined by the first component of principal component analysis (PCA). These sections captured the general shape of each cluster. A region-growing method was developed to merge clusters belonging to the same stem, representing a cluster by a sequence of its slice centroids. The nearest pair of clusters was tentatively combined and fitted with a 3D third-order polynomial curve. If the curve fitting error fell below a threshold, these clusters were merged. These steps were repeated for the remaining clusters until no suitable fitting was found. In practice, noise and adjacent stems could lead to incorrect curve fitting such as the lower part of one stem merging with the upper part of neighboring stems. To enhance the robustness of the 3D polynomial curve fitting and prevent over-merging across stems, a RANSAC selection scheme was employed, rejecting outlier centroids. The criteria for “goodness of fit” of the points along the curve were represented by the smoothness and consistency of the points surrounding the curve. The smoothness, represented by the absolute median deviation (MAD) between the fitted curve and slice centroids, should be less than 30 cm (corresponding to the maximum stem cylinder radius in this study). When fitting a pair of clusters, the median deviation of the curve fitting for one cluster was compared to that of the other cluster. The resulting ratio, which represents the curve-fitting consistency between the two clusters, was constrained to be less than 4. Depending on the practical situation, this ratio could be set more stringently such as 3. The RANSAC-based region growing merged well-fitted clusters iteratively until smooth and consistent stem curves were segmented from the stem classification points. The same method was applied to segment the downed woody logs into individual components.

3.4. Individual-Tree Clustering

In the study plots, it is common for trees to be densely packed and occasionally intersect with inclined ladder woods. Conventional individual-tree segmentation, based solely on point density, may underestimate fuel levels in areas where trees overlap or cluster together. Our method adapted the simple morphological rule that a tree stem forms the base for branch and foliage attachment. Among the classified trees, points were segmented into individual stem and non-stem points. The latter were transformed to the cluster scale using the cut-pursuit algorithm described above, but at a finer spatial resolution (K = 5, λ = 1). This resolution matched the typical size of branch and foliage clusters, where each cluster defined a graph node $N$. The non-stem clusters were assigned to the nearest stems based on a closeness measure adapted from the method by Wang [58]: (1) any non-stem nodes within 10 cm of the stem points were considered ‘seed nodes’ ($N_s$), while other nodes were considered to be ‘non-seed’ nodes $(N_n$); (2) the shortest path ($L_{\mathrm{m}\mathrm{i}\mathrm{n}})$ to each seed node $N_s$ was calculated for each $N_n$ using Dijkstra’s algorithm [84]; (3) two distances ($D_c$ and $D_p$) along the shortest path $L_{min}$ were calculated for each $N_n$, where $D_c$ was determined by the accumulated distance from the centroid of $N_n$ to the centroid of $N_s$, and $D_p$ by the accumulated 3D point gaps from the point cloud of $N_n$ to the point cloud of $N_s$, respectively; (4) the projected distance $D_m$ between $N_n$ and the medial line of $N_s$ was calculated; (5) a distance score, defined by a weighted average of the three distances $D_c$, $D_p$, and $D_m$, was determined as the closeness measure; (6) an $N_s$ was clustered to the $N_s$ with the minimum score, and both were assigned with the ID of the stem nearest to the $N_s$. The weights of $D_c$, $D_p$, and $D_m,$ defined as 0.4, 0.3, and 0.3, respectively, were customizable but did not require lengthy optimization for satisfactory individual-tree clustering. In practice, individual stems from previous steps could contain errors. A typical error is the misclassification of thick branches in the crown as a stem and isolated as an incorrect stem in the classification. Our amendment was to repeatedly merge two individual trees if the bottom of one tree fell above the top of the other tree with a projected horizontal overlap of >80%. Overall, individual trees were isolated using the shortest-path and conditional merging rules.

3.5. Branch Classification

A similar deep learning model for branch classification is also feasible, however, creating sufficient branch references to train deep learning from fuel-rich plots is laborious and impractical. Furthermore, human-interpreted thin branch samples can be erroneous and difficult to identify within a dense and crowded canopy. This study proposed an alternative solution to deep learning, circumventing the dilemma of reference making while potentially achieving a similar level of accuracy compared to the branch filtering accuracy of 0.54 reported by Xi et al. [85].

Ideally, branches are linear geometrical features detectable using a linearity-based filter. The linearity metric can be derived from PCA components [71]. The challenge is the weak signal of smaller diameter branch points surrounded by foliage and shrubs. The point clouds of larger diameter branches might also become fragmented and indiscernible due to the high occlusion of points by objects closer to the TLS. In this situation, the geometric identity of a branch can diminish, but its spectral identity remains visible, characterized by the intensity of laser returns [86]. The intensity alone, however, is not a robust indicator of branches when there are other environmental factors that influence the return intensity (e.g., local surface ruggedness) [87]. Therefore, a weighted-linearity-intensity index (WLI) combining the return intensity and the linearity of the branch was adapted in this study (Equation (1)):

$$\mathrm{W}\mathrm{L}\mathrm{I}=w\times \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{d}\left(\mathrm{L}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{y}\right)+\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{m}\mathrm{o}\mathrm{i}\mathrm{d}\left(\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y}\right),$$

(1)

Here, a weight parameter $w$ was used to leverage the importance of either term, which was set to 1.0 for simplicity. Approximate branch points were classified using a WLI > 0.75.

Despite this initial filter, the classified branch points were still noisy. A Frangi vesselness (FV) filter, which is used in multi-scale vessel detection from 3D X-ray medical images [88], was adopted to enhance the branch pattern. In this study, the filtered point clouds were voxelized at a resolution of 1.5 cm. Five levels of scales (2, 4, 6, 8, and 12 times voxel resolution) were input to the FV filter corresponding to the range of branch thicknesses between 3 and 18 cm. A 3D voxel morphological operator was then used to remove the isolated blank voxels surrounding them. Branch points were then filtered within each individual tree. Ten pine trees within the twelve plots were manually cleaned and delineated based on the results of branch filtering as the validation reference. It is important to note that while the WLI, noise filtering, and QSM branch delineation collectively contributed to mitigating the noise issues associated with branch identification solely based on intensity, they did not fully address the attenuation of intensity signals caused by the distance and angle incidence. While multi-scan co-registration alleviated the attenuation effect, the inclusion of intensity calibration is expected to further optimize the accuracy of branch detection.

3.6. Tree Architecture Reconstruction

Tree architecture reconstruction, also known as QSM, not only includes the 3D geometric modeling of branches and stems using cylindrical proxies, but also the establishment of the topologic relationship between wood components. Each individual branch requires isolation and to be subsequently arranged in a hierarchical order, beginning with the main branch connected to the stem, followed by the smaller branches attached to it, and continuing in this manner. This is achieved in a similar manner to the tree clustering approach (Section 3.4). Here, branch points are first clustered using the cut-pursuit algorithm. Branch clusters near stem points are considered as seed nodes. The shortest paths from the remaining nodes to the seed nodes are generated with a maximum point gap between two clusters of 3 cm. The shortest path with the greatest accumulated distance is classified as an individual branch, and all nodes along this path become branch nodes. As an example, Figure 6 demonstrates that the parent node of a branch is the stem node closest to the first node of the path. Likewise, nodes along the path with the next greatest accumulated distance form the second branch. The second branch discards any nodes already belonging to other branches, where the parent node is the joint node (i.e., the last node of the discarded part) of the previous branch. The branch generation process is repeated until no close clusters (gap < 3 cm) are found. The remaining isolated clusters with a gap between 3 and 30 cm are further attached to the existing branch architecture using the same shortest path rule, except that the seed nodes are closest to the existing branches instead of to the stem. Remaining isolated clusters are discarded. All individual branches are fitted with a sequence of circles [89] every 10 cm (customizable) using least-square minimization, which determines the branch direction and radius at each segment [85]. The branch radii occasionally encounter large values due to the inclusion of locally noisy points included in the circle fitting algorithm. To address this issue, the impact of outlier points on the branch radius was reduced using a moving-window median filter along each branch, employing a maximum cut-off value based on the median of all branch radii and a minimum cut-off value of 2.5 mm for the branch radius. This approach mitigated the problem of volume overestimation for small branches with diameters < 5 cm as reported by Demol et al. [90].

3.7. Fuel Attribute Extraction and In Situ Validation

The above-described classifications of specific plot components such as stems and branch points are not mutually exclusive. Some points can be associated with more than one class. Consequently, an additional step was needed to synthesize the classes and prevent ambiguous one-to-many class mappings. Specifically, the tree, within-tree stem, branch, and foliage classes were labelled with class 1, 2, and 3, respectively. Non-tree points mainly included subcanopy shrubs, sapling layers, and surface grass. From the viewpoint of fuel analysis, these points were treated as a whole without the individual segmentation used for individual trees. Point clouds of downed woody logs within the non-tree class were assigned class 4. Stems and branches detected within the non-tree class represented sapling stems (class 5) and branches of seedlings and shrubs (class 6), respectively. The elevations of the lowest points per 3 × 3 m² area of the non-tree class were extrapolated into a plot-size grid as a digital elevation model (DEM). Any points within the 20 cm buffer of the DEM were filtered as ground points including litter and surface grass (class 7). The remaining points presented tall grass and any below-canopy foliage (class 0). The class labeling followed a priority order to avoid ambiguity: class 1 > class 2 > class 3 > class 5 > class 4 > class 7 > class 6 > class 0.

This multi-class plot delineation enabled the fine-scale attribute extraction of separate fuel components. The stem tilt angle (TA) was quantified using the third coefficient of the first principal component of the stem point clouds from the PCA analysis. Ladder stems (including stems from shrubs and saplings) were filtered from individual stems with a TA > 10°. The ladder stems with any attached branches constituted ladder trees. The live crown base height (LCBH), which was measured in the field, was approximated based on the height of a point with the maximum curvature along the height-specific voxel volume profile, in other words, the knee point or inflection point similar to Popescu and Zhao [91]. Branch and stem volumes were determined as the total volume of 3D fitted cylinders from the QSM model described in Section 3.6. At the plot level, it is difficult to reconstruct complete geometries of downed woody logs and branches near the ground and in the duff layers of the soil due to the high occlusion of returns and complex shapes. The volumes of these components were inferred based on the linear relationship between the voxel volumes of the stems and branches, namely, SVOL and BVOL, and their QSM volumes of the counterparts, namely, QSVOL and QBVOL.

In situ measurement of tree heights, DBH, and LCBH were used to partly validate the accuracy of the diameter and length extraction from the delineated plot components using TLS. RMSE and Pearson’s r were used as the accuracy metrics. DCBH, while measured, was not used here due to the interpretation of what constitutes the height of dead branches vs. live ones (as a proportion), and therefore there is measurement uncertainty due to the vision occlusion from the under-canopy and shrub layers. Measured trees from each field quadrant (separated by cardinal coordinates) were manually linked to the isolated TLS trees based on tree orientation, DBH, LCBH, species, and tree height.

4. Results

4.1. Classification, Segmentation, and Filtering Accuracy

The SegFormer tree classification is visualized in Figure 7 with trees in red. Large variations in the tree distribution pattern and density were observed in the classified tree points due to time since post mountain pine beetle attack (and phase). JP8 had the visually sparsest tree distribution and JP9 the densest, also observed in the stem density survey in Table 1. A species mixture of mature pine, fir, and aspen was also observed in JP5 and JP6, and fir saplings were visible in JP7, JP8, JP11, and JP12. The crown apices of fir saplings were prone to misclassification as trees (e.g., in JP 8 and JP11) due to their intermediate height. Note that the intent of tree classification might be identical to the ground removal scenario [67] where algorithms such as the cloth simulation filter (CSF) [92] abound. The subtle difference is that a tree classifier needs to identify tree points from the subcanopy complexes instead of simple height slicing using the ground DEM buffer (e.g., CSF). The stem bottoms within the DEM buffer need to be preserved as the stem class. Deep learning addressed this complexity, yielding a high mIoU accuracy of 0.91 from four testing plots (

Table 3. Average testing accuracy and total processing time of the twelve plots for each main step. The accuracy of QSM wood reconstruction was not investigated directly.

Process	Accuracy (mIoU)	Processing Time (Minutes)
Tree point classification	0.91	38
Stem point classification	0.94	48
Downed woody log point classification	0.81	70
Stem and downed woody log isolation	0.73	61
Individual-tree clustering	0.86	259
Branch filtering	0.62	1193
QSM wood reconstruction	-	1534

).

The stem point classification was contrasted with the reference in the example JP9 clip in Figure 8. The SegFormer could identify most stem points with slight over-detection: major stems were thickened, and additional linear shapes in the crown were misclassified as stems. The former edge issue could be caused by the setup of a slightly coarse voxel resolution of 4 cm and the latter inscrutable noise is a common issue of deep learning [93]. The challenge was to detect sparse stem signals within dense, needle-leaf fir crowns where existing studies either excluded the classification of fir species or demonstrated insufficient accuracy [94], however, for wildland fire, fir is an important species because it is highly flammable. These issues were mitigated in Figure 8, where the mIoU of the fir stem classification was 0.89, and the stem classification over all testing samples was 0.94 on average (Table 3). The example classification results of downed woody logs are also visualized in Figure 8. SegFormer tended to detect outlying small, round-shaped pieces of wood in addition to the reference ones. Using the manual reference as the baseline, the mIoU of woody log classification was 0.81. Figure 8 also shows the example of stem and downed woody log classification over JP1. From the birds-eye view, not all coarse woody log points from the reference were captured by SegFormer, especially among the oblique stems and highly-occluded thin ground logs. This issue partially arose because the manually created training samples did not represent sufficient instances of oblique or highly-occluded woody logs, and the deep learning model was not adequately trained to handle such intricate situations. While increasing the representation of these complex instances in the training dataset could potentially address this issue, such an approach would require further manual editing, a resource-intensive task that was beyond the scope of this study.

The example results of cut-pursuit clustering and individual wood isolation are illustrated in Figure 9 using a JP2 clip as an example. The classified stem points were manually cleaned before the clustering and isolation. Small stem clusters with linear shapes were generated by cut-pursuit and further merged into complete wood segments using the RANSAC-based region growing method. Some bending and occlusion (Figure 9c) could be overcome due to the polynomial fitting of the RANSAC method. The intersection of two or more stems or downed woody logs, a typical barrier for correct segmentation, can also be addressed to a certain extent using our method. Issues with the segmentation included the failure of merging several small irregular or highly curved clusters. Plot-level tree clustering based on individual stems is illustrated in Figure 9d–f using JP1 as a second example in this figure. The classified stems were uncleaned, resulting in an issue of over-isolation (Figure 9e), where more stems were isolated than actually exist. This could be mitigated by the use of overlap-based merging during tree clustering (Section 3.3) (Figure 9f). Many trees that were intersected in the crown could be separated with well-delineated boundaries due to the morphology-based shortest path rule (Figure 9f). Accuracy assessment did not include incomplete and sparse trees near the edge of the plot caused by the partial scan coverage. The mIoU of stem isolation and individual tree clustering was 0.73 and 0.86, respectively.

The branch processing and classification results are shown in Figure 10. The WLI index, which integrates the intensity (Figure 10a) and linearity (Figure 10b), filtered out dense foliage (Figure 10c, without stems), and the FV filter smoothened and enhanced linear shapes (Figure 10d). Figure 10c also reveals that a fraction of the twig points was also filtered out, which can underestimate the fuels available for fire as these burn more readily than larger diameter branches. There was a tradeoff between the partial loss of branch points and the need to exclude foliage and noise points from the classification. The plots in our sites consisted of mostly mature, high foliage trees, which also contributed to a high level of return noise. Therefore, a high threshold of 0.75 was set for the WLI, leading to the filtering of branches from fine branch material and needles. Using our testing reference, branch filtering had an mIoU of 0.62.

Following the synthesis of classifications into eight classes, the fine-scale composition of a plot could be represented, as demonstrated in detail for a small area of the JP7 clip (Figure 11a) and across the larger plot using JP1 (Figure 11b). It is evident that foliage in the dark green color (class 3) was primarily distributed among fir trees, short pine trees, and the lower portions of tall pine trees in JP1. Occlusion occurred mostly near the plot edges and on the ground near the corner TLS stations. Understory grass and herbaceous layers indicated by light green (class 0) tended to have an overstory of large crown diameter trees or tree clusters, which enhance shading and reduce water loss, while large spaces between trees had shorter/less dense grasses and was primarily represented by ground returns. The horizontal orientations of ladder trees appeared randomly distributed.

The complete QSM reconstruction at the tree- and plot levels is presented in Figure 12. Results of the branch clustering based on the shortest path rule, skeletonization, and cylinder-fitting are illustrated in Figure 12a–c, respectively, based on the same point cloud in Figure 10 (JP7) and the 3D wireframes from QSM reconstruction (Figure 12d,e) based on the point clouds in Figure 11 (also for JP7). Here, narrow branch diameter branch structures were also captured by the 3D reconstruction. A notable issue was the exaggeration (or artificial increase) in the reconstructed stem diameters near the apex of the tree crown due to the noise points during the circle fitting.

4.2. Field Validation

A total of 269 trees with field measurements were matched to isolated TLS trees within the twelve plots. From these, a subset of 84 trees had live crowns, which was used to calculate LCBH. The r² between the field measurements and TLS extraction was 0.95 for tree height, 0.96 for DBH, and 0.60 for LCBH, respectively (Figure 13). RMSE (RMSE%) was 1.3 m (9.0%) for tree height, 1.3 cm (8.7%) for DBH, and 2.6 m (39.6%) for LCBH, respectively. The bias of TLS height, DBH, and LCBH was −0.53 m (−3.8%), 0.54 cm (3.5%), and −0.35 m (−5.3%), respectively, compared to the field measurement. TLS tended to underestimate the LCBH. The differences are algorithm- and species-specific, and are not consistent among the different studies that have used TLS [95,96].

When comparing the QSM approach for volume calculation, the wood volume determined by counting voxels was subject to the choice of voxel size and the completeness of the scanned surface. Utilizing a voxel size equal to the point cloud resolution of 5 mm led to the underestimation of individual tree wood volume by 81% (branch)–92% (stem), respectively, compared to the wood volume derived from QSM reconstruction. Increasing the voxel size might mitigate the underestimation but would introduce complexity in adjusting the voxel size for different tree and plot components. A linear relationship between these two volumes was observed excluding minor extreme outliers (Figure 14). The r² was 0.54 between SVOL and QSVOL from all individual trees, smaller than the r² of 0.95 between BVOL and QBVOL. The diameter of stems was mostly greater than 5 cm, and simple voxelization of stem surface points did not account for the volumes within the stem body. Conversely, branch diameters were smaller with many contained by voxels, causing better agreement between the voxel and QSM volumes.

4.3. Plot-Level Attribute Analysis and Volume Calibration

At the plot level, the downed woody volume, sapling stem volume, and ladder stem volume were scaled from the voxel volumes based on the ratio of QSVOL to SVOL, and the below-canopy branch volume and ladder tree branch volume based on the ratio of QBVOL to BVOL. The result component volumes and other plot-level attributes are provided with a summary chart in Figure 15. Except for JP12, which was heavily affected by MPB, JP9 and JP10 had the highest concentration of stem to the total plot volume (58.3% and 53.2%, respectively) in the steepest slope (mean slope angle = ~20°), and the lowest branch volume proportions to the total tree volumes (19.5% and 23.3%, respectively). This aligned with the visual appearance in Figure 7 that trees in JP9 and JP10 were dense and straight but with narrow individual tree crowns. JP5 and JP6 presented another distinct plot configuration: mixed species (pine% = 65–75%, fir% = 13–16%) and a thick sapling layer (~4% of the total plot volume compared to the average of 2% among other plots). However, the concentration of ladder tree volume in both plots remained at an average level compared to other plots. In contrast, JP7 and JP8 had only four ladder trees and few saplings. JP8 had the lowest wood volume concentration (60.7%). JP1 and JP2 had the highest ladder stem volume (4–6% compared to the average of 2% among other plots). The ladder stem volume was strongly associated with stem density (r = 0.85), existing straight stem volume (r = 0.49), and DEM slope (r = 0.81), indicating a positive effect from competition and gravity factors. The vertical distributions of each fuel component with slope correction are visualized as stacked area charts in Figure 16. It can be seen that the plot-level composition varied significantly with geographic locations but remained similar between neighboring plots. The stem volume distribution mimicked a Weibull curve, descending at a faster rate toward the tree top than the branch volume distribution. The vertical distribution of foliage volumes was similar to that of branch volumes. Ladder fuels and saplings contributed to a negligible proportion at the plot level in terms of tree volumes in Figure 16. In most plots (JP2, 3, 7, 8, 9, 10, 12), a boundary existed between the tree branch layer and the below-canopy branch layer. The presence of a ladder tree layer between them created a noticeable connection between the two potentially combustible layers. The average overlapping area of the vertical volume profile excluding the ladder tree layer was 8.4% relative to the union area of the tree layer (stem, branch, foliage) and the below-canopy layer (all other components). Including the ladder layer, the overlapping degree increased to 10.8%.

5. Discussion

The study of forest component classification from TLS has gained prominence in the research community over the last decade. Among the available methods, many geometric algorithms have demonstrated excellence in stem detection [97,98,99,100] and wood and foliage classification [70,101,102,103], which are alternatives to the emerging deep learning solutions [71,85,104]. However, the limitations of this method are also obvious, whereby rule-based algorithms are too specified to be adaptive to complex environments, which include standing and ladder trees, saplings, and downed woody logs. In comparison, deep learning classifiers can exhibit flexibility without sacrificing accuracy. In this study, SegFormer was employed as a binary classifier (e.g., tree vs. non-tree, and stem vs. non-stem). Altering SegFormer to function as a multi-class classifier, which simultaneously classifies stems and downed woody logs together, resulted in a reduced accuracy. The mIoU decreased by 0.05 (6%) to 0.82 compared to the average mIoU of 0.87 when classifying these components separately. This finding coincides with a recent study by Li et al. [105] where binary CNN-based classifiers exceeded the multi-class situation in wetland vegetation classification with an mIoU increase by 22%. However, the processing time would significantly increase if each class is trained individually.

In addition, deep learning was not used for branch classification due to the low efficiency of branch sample delineation from dense conifers. The detection of fine branches in spruce remains an issue. The geometrical algorithm developed by Hui et al. [70] has demonstrated overall superiority in branch filtering, with the exception of narrow spruce branches. The filter in our study could detect twigs but also contains cones and some noise as a function of the point cloud distribution. Several tree species such as the subalpine fir and Douglas fir were not widely sampled within the training data but coexisted with other conifers found in the plots, adding to the uncertainty of the stem classification. In addition, one stem per tree was the underlying assumption of our stem isolation algorithm, which could be violated by the “co-dominant stem” strategy where two leader stems grow from single trees to withstand environmental stressors [106]. This phenomenon was found occasionally within plots and could potentially result in the underestimation of the true stem volume. It was also observed in the field that dead pine trees in grey phase following mountain pine beetle attack retain some mosses and lichens on the branches. Lichens could be mistaken for foliage from TLS without field validation [107] and may increase return classification biases to increase the branch volume.

The result RMSE% of tree height and DBH in Figure 14 fell within the average range of 13–30% and 5–15%, respectively, a common range from a benchmarking study on eighteen contemporary methods in 2018 [108]. Tree heights were underestimated by TLS (Figure 14), as also observed in other studies [95,108]. The near-zero biases of the TLS DBH aligned with studies [108] using similar extraction methods. Compared to the height and DBH, the TLS LCBH had weaker correspondence with the field measurement. Causes could be an unclear definition of LCBH from an abnormal crown shape, especially the unhealthy trees [109], and the fallible measurement of crown bases from heavily overlapped crowns [110]. From the literature, the former cause was insignificant regarding the strong r² (>0.9) between the human-interpreted and the automatically extracted crown base height (CBH) from TLS by Seidel et al. [111]. The field measurement error of CBH was considered the main issue [112]. The literature with TLS CBH validation is rare due to the strenuous efforts of obtaining CBH from the field [112]. Among the available studies, Jung et al. [110] reported a high r² of 0.82 but had a limited validation sample size of fifteen trees. Fernández-Sarría et al. [113] validated the LCBH accuracy with thirteen field trees and reported an r² of 0.52, in line with this study. A stand-level assessment of CBH accuracy was conducted by Pyörälä et al. [114] with a relatively lower validation r² of 0.20. All of these studies demonstrated a significantly stronger r² from the height and DBH validation compared to CBH, indicating a reliability issue of field mensuration to determine CBH.

6. Conclusions

The capability of TLS has facilitated the enhancement of our understanding of fuels within forests, from landscape- and stand-based scales to a more nuanced focus on individual trees and branches. While TLS is gaining wider application in fire fuel studies as a replacement for inventory mensuration, it is primarily based on voxel profile information, without the in-depth incorporation of advanced data processing and AI techniques. An analysis of the fine-scale forms of trees and branches based on AI methods provides a more detailed understanding of the intricate nature of forest fuels and their connectivity, which may be further assessed along with local climate including the drying potential. This study broadens the extent of existing TLS investigations to characterize fuel components at the plot scale based on deep learning, incorporating ladder trees, fallen woody logs, subcanopy branches, and surface layers. Twelve TLS plot scans were collected from a C3 fuel type with lodgepole pine that were, in many plots, at various stages of mountain pine beetle attack and mortality. The TLS extraction of tree height, DBH, and LCBH matched the in situ mensuration, while TLS provided additional detail of the volume distribution of each plot component along the vertical height. This scheme of component division could help refine the fuel load quantification among materials with different flammability and also indicate the vertical continuity of different fuel components. The existence of a ladder tree layer enhanced the overlap degree between the two combustible layers of individual-tree branches and below-canopy branches from 8.4% to 10.8%. In addition, a significant underestimation of wood volume based on conventional voxel-based counting at fine resolution was found compared to the more reasonable QSM method due to scan occlusion. It is suggested that QSM is applied to calibrate voxel-based wood volumes for fine-scale fuel analysis.

The differentiation of forest components serves as a foundational dataset for a range of fire-related studies. With respect to fire fuels, the grass and litter layer is associated with fine fuels, coarse fuels represented by individual-tree stems and downed woody logs, and intermediate fuel loading embodied by branches. The fraction, length (or height), and diameters of individual fuel components provide essential input in the allometric equations of fuel load [115]. Overall, the parameters of the wood components and terrain hold significance for the planning (e.g., prescribed burns), the potential for ignition, and the prediction of the fire rate of spread. The proof-of-concept in this study illustrates the bright future of TLS in accurately quantifying the fuel physical properties for wildland fire analysis.