Detection of standing retention trees in boreal forests with airborne laser scanning point clouds and multispectral imagery

In a landscape consisting primarily of intensive forestry interspersed with some protected areas, multifunctional forestry with retention trees can play a crucial role in nature conservation. Accurate mapping of retention trees is important for guiding landscape‐level conservation and forest management and improving landscape connectivity. Sizeable dead and living retention trees play a particularly important ecological role but even their large‐scale inventory is often intensive through field work and/or inaccurate. We aimed to detect and classify retention trees using the novel nationwide Finnish airborne laser scanning (ALS) data (~5 pulses/m2) in conjunction with unrectified colour‐infrared (CIR) aerial imagery. Applying photogrammetric principles, we added spectral information from the CIR imagery to the ALS‐derived point cloud. For a training dataset of 160 retention trees from 19 stands and a geographically separate validation dataset of 79 trees from eight stands, we segmented trees via individual tree detection (ITD), removed most trees belonging to the regenerating vegetation layer, and classified trees into living conifers, living broadleaves and dead trees by linear discriminant analysis. The detection rate via ITD differed considerably for dead and living trees, with 41.7% of all dead and 83.8% of all living trees being detected with relatively low commission error rates. Dead trees with smaller diameters and heights were more likely missed, while grouping caused living tree omission. For classification into living conifers, living broadleaves and dead trees, an overall accuracy of 67.3% was achieved in training and 71.2% in validation using only ALS‐derived metrics. When adding spectral metrics, the overall accuracies were 79.6% and 61.0% for training and validation respectively. Our findings imply that wall‐to‐wall large‐scale high density ALS data can be used to detect retention trees rather accurately—even larger dead trees—and that metrics derived solely from ALS data can accurately classify detected retention trees into living conifers, living broadleaves and dead trees. Considering the ecological value of retention trees, our results are promising and indicate that ALS data of the studied pulse density are a cost‐effective option for large area mapping of retention trees in countries with such data available.


| INTRODUC TI ON
As part of ecological considerations in forestry, the practice of leaving retention trees during clearcutting operations has been heavily employed in recent decades across the globe . Retention trees are trees that have been and will be left unharvested during felling operations, either standing or fallen, whole or broken. These trees have proven valuable for several groups of species (meta-analyses by Fedrowitz et al., 2014;Mori & Kitagawa, 2014;Rosenvald & Lõhmus, 2008), being able to maintain some of the biodiversity present in the stand prior to final felling, act as lifeboats, and potentially serve as stepping stones for dispersal in otherwise heavily fragmented landscapes.
The effectiveness of retention trees is strongly influenced by the types of trees left and their positioning in the stand and the wider landscape. Key tree-level attributes that affect the value of a retention tree include tree species, size and condition, that is, dead, weakened or living, and on a within-tree level, the presence of tree-related microhabitats (as defined by Larrieu et al., 2018). Of particular importance are large trees rich in microhabitats (Bütler et al., 2013), dead trees that play a key role in the protection of saproxylic species, including rare and red-listed species (Gustafsson et al., 2020), and certain tree species that are otherwise considered ecologically valuable (Koch Widerberg et al., 2018;Lundström et al., 2013;Murray et al., 2021).
Besides tree-level attributes, the positioning of retention trees plays a crucial role in the effectiveness of retention trees, both within a stand (e.g. grouped vs. single tree retention and vicinity to intact stands, Franklin et al., 2019;Gibbons et al., 2008) and in the wider landscape (Mori et al., 2017;Vandekerkhove et al., 2013).
Due to their ecological value, retention trees should ideally be accurately mapped, and such mapping can in turn guide landscape-level conservation and forest management and improve landscape connectivity. Inventories of retention trees through traditional field work means are, however, time-consuming and costly. Furthermore, inventories of retention trees should preferably be repeated frequently due to high tree fall rates (Rosenvald et al., 2019). Remote sensing may present a cost-and time-effective solution to map and characterize retention trees. Since retention forestry is applied in multiple countries worldwide Lindenmayer et al., 2012;Martínez Pastur et al., 2020), this methodology is also considered globally relevant.
Mapping of retention trees and identifying at least some tree-level attributes through remote sensing may also assist in expanding ecological research on retention trees to the landscape scale; an area of research that has not been explored well thus far (Gustafsson et al., 2020).
To map retention trees most efficiently, trees should be recognized individually via individual tree detection (ITD; Hyyppä et al., 2001) and separated into various classes based on tree-level attributes. In general, individual living trees can nowadays be reliably and cost-effectively detected and their crowns delineated on large scales. Different remotely sensed data can be used for this purpose; for example, airborne laser scanning (ALS) data with sufficient point density (for large living trees 2 points/m 2 is sufficient, Kaartinen et al., 2012), optical images taken with uncrewed aerial vehicles (UAVs) or aircrafts that can be converted into photogrammetric point clouds (Hardenbol et al., 2021;Nevalainen et al., 2017), or combinations of the above. Conversely, detecting individual standing dead trees generally presents a much more difficult task requiring relatively costly data. Nevertheless, promising results have been obtained with ALS data at a relatively low point density of 6.7 points/m 2 (Wing et al., 2015) and with aerial imagery at a ground sampling distance of 20 cm (Zielewska-Büttner et al., 2020).
Considering retention trees, their detection, compared with trees in a dense forest, might be facilitated by their positioning as they are often few and far between. Furthermore, most retention trees, if they do not fall, will be taller than the main canopy for many decades following final felling. As an alternative to the ITD approach, an areabased approach (ABA) could provide similar stand-level information such as retention tree volume. However, ABA is not well suited for predicting the actual numbers and spatial locations of the retention trees, as it only provides averaged predictions at grid cell level.
Besides detecting retention trees, separating the individual trees into several classes with different tree-level attributes from remotely sensed data would provide further beneficial information.
Herein it is important to consider that tree species, condition and size are especially important. Tree species classification has garnered a lot of attention and has been successfully implemented with both passive and active remotely sensed data, or combinations of both (reviewed by Fassnacht et al., 2016). However, the focus has often been on separating broader classes like conifers from broadleaves, especially with information derived solely from ALS data which are insufficient for detailed species separation (Ørka et al., 2009;cf. Persson et al., 2022). The separation of living and dead trees has worked particularly well with high-spatial-accuracy optical data like hyperspectral imagery with the aim of providing tree-level health assessments (Näsi et al., 2018) but spectral values from colour-infrared (CIR) imagery fused with ALS data (Zielewska-Büttner et al., 2018) and even just intensity values from ALS data have shown the capability to distinguish living and dead trees (Kim ALS data of the studied pulse density are a cost-effective option for large area mapping of retention trees in countries with such data available.

K E Y W O R D S
coarse woody debris, forest management, individual tree detection, landscape connectivity, multifunctional forestry Wing et al., 2015). Finally, detailed size measurements like height, diameter at breast height (DBH) and volume have been provided with several methods although with method-specific limitations (Hyyppä et al., 2001;Wang et al., 2019).
As indicated above, ITD has been successfully employed, both for living and dead trees. However, studies aimed at dead tree detection have been using demanding high-resolution ALS or optical data and such data are not readily available and are costly to obtain.
In various countries, including Finland, there are efforts underway to collect high density (~5 pulses/m 2 ) nationwide ALS data. It is important to research whether such data are sufficient for a multitude of applications as it would save on valuable resources. As mentioned above, ALS data with a spatial coverage of around 2 points/m 2 have proven sufficient for ITD of living trees (Kaartinen et al., 2012), but this requirement may be considerably higher for dead trees and it is currently uncertain whether 5 pulses/m 2 is sufficient. However, while ALS-derived values such as intensity have proven sufficient for the separation of certain tree species (Ørka et al., 2009) and tree conditions (Wing et al., 2015), more accurate classifications of tree species and tree conditions could be obtained when fusing ALS data with spectral information (Kukkonen et al., 2018;Zielewska-Büttner et al., 2018). In some Nordic countries, including Finland, it is standard practice to collect multispectral aerial imagery in conjunction with ALS data for the purpose of species-specific forest inventories  and updating of topographic maps. As such, it is a cost-effective option to fuse ALS data with this spectral imagery and explore the capabilities of the fused data.
Our aim in this study was to determine the suitability of nationwide ALS data with a pulse density of 5 pulses/m 2 and the combination of these ALS data with colour-infrared (CIR) imagery (near infrared [NIR], red, and green bands) to detect retention trees, delineate tree crowns, and separate detected retention trees into living and dead (specifically living conifers, living broadleaves, and dead trees). To achieve these objectives, we used ALS-derived point clouds fused with unrectified CIR imagery via photogrammetry and field inventory data of retention trees. We then used an ITD algorithm and a watershed segmentation algorithm to delineate the crowns of detected trees. These individual tree segments were classified into living conifers, living broadleaves and dead trees based on intensity and relative height variables when only using ALS data and, additionally, spectral variables when using the fused ALS and CIR data which the returns within each segment contained.

| Study area
Our study area is in the municipalities of Liperi, Outokumpu and

| Field inventory
Two field inventories were performed to gather data for a training dataset (collected in May 2021) and a validation dataset (collected in October 2021) for which no permits were required. To include a sufficient number of standing dead trees in both datasets, we mostly inventoried stands that had at least one standing dead tree. We define retention trees as trees that stuck out at least 3 m from the regenerating vegetation layer. We included stands that underwent final harvesting at any time in the past and the regenerating vegetation layer thus ranged from 0 to 13 m (recent to old final felling).
All standing retention trees above 9 cm in DBH were inventoried.
We selected 9 cm DBH as a cut-off point because it is reasonable to assume that thinner dead trees are likely missed frequently by the laser pulses. For suitable retention trees, we recorded tree species, conditions (dead or living), positions, heights, diameters, and only for dead trees, stem and crown conditions. Tree species that were recorded were Norway spruce, Scots pine, silver birch, downy birch, European aspen, rowan, goat willow and grey alder. Note, however, that we grouped broadleaves and conifers to obtain good samples sizes for the classifica- in two distinct ways. If a tree was intact, DBH was measured at 1.3 m height above-ground. If a tree was broken, diameter was estimated at the middle height of the remaining standing tree column. For dead trees only, we recorded if the stem was intact or broken (either naturally or artificially), and if the stem was intact, whether the crown was entirely present, partially present or absent. In total, we collected field data on 160 retention trees from 19 stands for the training dataset and on 79 retention trees from eight stands for the validation dataset (Table 1).

| Remotely sensed data
Two remotely sensed datasets were used in this research, both covering our entire study area. ALS data with an average pulse density of 5 pulses/m 2 and unrectified colour-infrared aerial (CIR) imagery.
The ALS data were acquired on 12-14 June 2020 and their col- This scanner records up to four returns (echoes) per pulse, which were later categorized as first of many, intermediate, last of many and single returns. The multipulse mode was turned on. In this scanner, the divergence (1/e 2 ) of the laser beam (λ = 1064 nm) is 0.25 mrad, which leads to a footprint of about 0.5 m.
The CIR imagery was collected on 14 June 2020 with an UltraCam Eagle Mark 1 camera (Vexcel Imaging) at an average altitude of 7821.7 m above ellipsoid. This altitude resulted in a nominal ground sampling distance of 0.4 m. The system was equipped with a 100.5 mm focal length lens and a 4360 × 6670-pixel (15.6 μm pixel size) sensor for colour images (red, green, blue and NIR bands) and a sensor for panchromatic images that was irrelevant in our current

TA B L E 1
The reference trees of the training (19 stands) and validation (eight stands) datasets for the remotely sensed detection of retention trees.

| Processing of remotely sensed data
Applying photogrammetric principles, we added spectral information (red, green and NIR bands) from the unrectified multispectral imagery to the ALS-derived point cloud (Packalén, 2009). This was done through collinearity equations with knowledge of camera locations and orientations at the time of exposure of each image and internal camera parameters. The per channel digital number value for each pixel was averaged from all the images where a 3D point was observed. The result of this process was a point cloud that included both XYZ coordinates and spectral information for each return.

| Individual tree detection and computation of predictor variables
The point cloud elevations were normalized into above-ground level with the publicly available ALS-based 2-metre resolution elevation model provided by the National Land Survey of Finland. With normalized elevations, we first removed some clear outliers/noise with high height values. Next, canopy height models (CHMs) with a 50 cm resolution were constructed from the normalized point clouds using the LAStools software (Isenburg, 2020) and the pit-free method presented by Khosravipour et al. (2014). This method avoids the generation of no-data areas inside the CHMs by stacking a set of temporary CHMs constructed with different height thresholds.
The ITD was performed with the R software version 3.5.3 (R Core Team, 2019) using appropriate functions from the libraries LidR (Roussel et al., 2020), rLiDAR (Silva et al., 2017) and ForestTools (Plowright & Roussel, 2018 containing a field-measured tree, the segment was linked with the largest field-measured tree within it. Note that some trees had fallen between the ALS data inventory and our field inventory (n = 14).
Such trees showed as commission trees, but we decided to remove them from the results since they were clear artefacts.
Returns within segments with maximum heights above what we  (Yu et al., 1999).
Finally, each predictor was normalized to the standard unit interval which made the model interpretation easier.

| Statistical analyses
We first investigated to what degree tree-level variables (diameter, height, stem and crown conditions and tree species) affected the detection rate for dead trees. We used logistic regression with a stepwise variable selection to test the effects of these variables on detection after the assumptions of the test were validated. This analysis was only performed for dead trees since visual assessment revealed that the few living trees that were left undetected were in dense retention groups whereby some trees were not individually segmented.
The variables extracted from the segments that were linked with the field-measured trees from the training dataset were imported into R for classification into living conifers, living broadleaves and dead trees with linear discriminant analysis (LDA). The LDA model was trained using the R-package MASS (Ripley et al., 2020). The number of predictor variables in the LDA was set at five, which is a compromise between the accuracy of predictions and avoidance of overfitting. We implemented the variable selection based on a heuristic optimization algorithm (Packalén et al., 2012) that is a modification of the original simulated annealing method (Kirkpatrick et al., 1983). This algorithm tested many different predictor combinations attempting to improve from the best currently found predictor set. To avoid convergence on local optima in the variable space, worse solutions also had a probability to be accepted. This probability of accepting worse solutions decreased linearly as the number of iterations increased and was controlled by the cooling parameter temperature. We iterated the algorithm 10,000 times with the initial temperature set at 0.2 (Packalén et al., 2012). The variable combination that provided the highest classification accuracy was selected for reporting. The accuracy was assessed based on the criterium of the greatest overall accuracy value. This provided the optimal variables for classification of all three categories. After an optimal classification model was produced from the training dataset, we tested this model on the validation dataset.
We calculated accuracy assessments of the optimal model based on confusion matrices, overall/user's/producer's accuracies, Cohen's kappa coefficients and macro-averaged F1-scores (mean of per-class Meanwhile, the overall number of commission cases was 38 trees with 22 of these resulting from living trees being split into two or more segments (Figure 2 tree detection was significantly affected by diameter (p = 0.01) and height (p < 0.001; Figure 3). It was found that, holding height con-

| Classification into dead and living trees
Classifying the detected retention trees into living conifers, living broadleaves and dead trees proved relatively successful. When using all predictors (spectral, height and intensity related) from the training dataset with leave-one-out cross-validation, an overall accuracy of 79.6% (Cohen's kappa: 0.69, macro-averaged F1-score: 0.79) was achieved (Table 3a). Testing the model derived from this training data on the validation dataset resulted in an overall accuracy of 61.0% (Cohen's kappa: 0.43, macro-averaged F1-score: 0.60; Table 3b).
When only employing ALS-derived predictors (height and intensity related), the classification accuracies changed considerably.
On the training dataset with leave-one-out cross-validation, an overall accuracy of 67.3% (Cohen's kappa: 0.50, macro-averaged F1score: 0.68) was achieved ( In both cases, the first discriminant separated dead trees from living trees (Figure 4). With all predictors, the second discriminant managed to separate living conifers and broadleaves rather reliably (Figure 4a), while this was not the case with only ALS-derived F I G U R E 2 A detailed view of part of a stand with two living conifers and a few dead trees and the segments resulting from the watershed segmentation algorithm. On top, the segments are laid over the canopy height model created from the airborne laser scanning data. Below, the same segments are (in the same colour and just for visualization) laid over a cross-section of the point cloud which portrays fused airborne laser scanning data and colour-infrared imagery. This figure clearly shows how dead trees can be left undetected (omission) with no returns above 2 m and living trees can be split into multiple segments when returns show more tops resulting in commission trees.
predictors where the second discriminant was of minor importance ( Figure 4b). With our validation data, the best possible separation of the three classes was achieved with tree shape and intensity characteristics (Table 4). The most important predictor separating dead and living trees was maximum crown radius when using either all predictors or only ALS-derived predictors. Since dead trees often lack a crown, it is expected that crown radius is much smaller for dead trees. With only ALS-derived predictors, however, crown radius was also instrumental in separating living broadleaves and living conifers, while intensity and height (characterizing crown shape) were also important in this separation.

| Methodology effectiveness
We found that retention trees can be detected at adequate rates and subsequently classified accurately into living conifers, living broadleaves and dead trees. This methodology presents a foundation to improve both future forest management for biodiversity and even ecological research. Directly comparing our findings with previous studies is difficult considering that we are dealing with different forest environments than previously dealt with in literature and unique remotely sensed data. Nevertheless, based on multiple previous studies, at an approximate point density of 5 points/m 2 , living trees are expected to be detectable at high rates (Kaartinen et al., 2012) and this is verified by our findings for living retention trees. For standing dead trees, Wing et al. (2015) suggested that a point density of at least 4 points/m 2 (first and single returns only) should result in a satisfactory detection rate for large trees. Our results corroborate this suggestion, and like the study by Wing et al. (2015), we acquired a high detection rate for large dead trees and detection became progressively poorer for smaller trees. Considering, however, that dead trees with a diameter below 35 cm were well represented in our data, it meant that many trees were left undetected.
This presents a clear argument for higher pulse density, especially when bearing in mind that canopy cover in this study is low, which improves detection rate (Wing et al., 2015). Yet, from an ecological perspective, while even small pieces of coarse woody debris are ecologically valuable (Suominen et al., 2018), large dead trees are certainly considered more important for biodiversity (Stokland et al., 2012). As such, accurate detection of these large dead trees is more urgent, and our methodology proves suitable for their detection. Another consequence of the employed pulse density on dead tree detection comes in the form of height underestimation. While we only inventoried retention trees that stuck out from the regenerating vegetation layer by 3 m in the field, some dead tree heights based on the ALS data were underestimated. Our regenerating vegetation layer rules thus removed actual detected dead trees.
F I G U R E 3 Dead tree diameter and height distributions and the detection rates along these distributions.
Classifying the detected retention trees based on tree-level attributes is an important next step following detection as different classes provide distinct biodiversity values. In our methodology, we limited the classification to living conifers, living broadleaves and dead trees, but this could be expanded upon, especially with more costly data. Regardless, our classification proved successful, and while we would expect that using CIR-derived information in addition to the ALS data would improve our aimed at classification While specific to our collected datasets, CIR-derived information might therefore not be essential to classify retention trees into living conifers, living broadleaves and dead trees.

| Additional limitations and future perspectives
One of the main limitations in our study is that several ecologically important tree-level attributes of retention trees were not distinguished with our methodology due to the limitations of the employed data. A more accurate distinction between retention trees, especially concerning tree species, would be useful considering the ecological value of different tree species. For example, European TA B L E 3 Confusion matrices performed to separate the three tree type classes for (a) the training data with all predictors and leaveone-out cross-validation, (b) the validation data with all predictors, (c) the training data with only ALS-derived predictors and leave-one-out cross-validation and (d) the validation data with only ALS-derived predictors aspen in European boreal forests is highly valuable for biodiversity as a retention tree (Lundström et al., 2013) which in our case is merged with other broadleaves. By using other remote sensing methods, data may be collected that provide opportunities for more detailed species classifications of living retention trees. This includes ALS data with a higher pulse density that would provide more detailed structural characteristics (Li et al., 2013), or multispectral and hyperspectral imagery collected at low altitudes that provide detailed spectral characteristics (Briechle et al., 2021;Hardenbol et al., 2021;Viinikka et al., 2020). Separating dead retention trees into the various tree species through remote sensing is more data intensive and has been accomplished, although at less detailed classifications than for living trees (Kamińska et al., 2018). Moreover, habitat trees containing a wealth of microhabitats may also be detected with remote sensing eventually, but thus far even terrestrial laser scanning with a very high pulse density has been poor in detecting microhabitats directly (Frey et al., 2020;cf. Rehush et al., 2018). Nevertheless, inferring habitat trees based on tree-level attributes presents a rather suitable alternative solution (Zlonis et al., 2022). Finally, fallen retention trees were also ignored in our analyses despite their ecological importance (Suominen et al., 2015). Fallen trees have been detected successfully in certain habitats with higher pulse density LiDAR data (Heinaro et al., 2021;Tanhuanpää et al., 2015) providing further grounds for such data.

F I G U R E 4
Class separation by two discriminant axes for classification into living conifers, living broadleaves and dead trees for (a) the training data with leave-one-out cross-validation and the validation data with all predictors (ALS and multispectral derived), and (b) the training data with leave-one-out cross-validation and the validation data with only ALS-derived predictors. The depicted points show the field-observed class of each tree. The circles shown represent the Euclidean distance from the class means at a level of 0.95 and 2.5.  Albeit less cost-effective for operational purposes, alternative data sources to those used in this study can certainly provide better detection and classification performances. These sources include conventional linear-mode airborne LiDAR data with a higher pulse density from any airborne platform and even single-photon LiDAR data (e.g. Räty et al., 2022). Balancing efficiency and cost are of the essence, however, especially if we consider application of our findings to large areas and other parts of the world, hence why we decided to employ nationwide data with a lower pulse density. It is worth noting that although the pulse densities for nation-and statewide ALS inventories are often around 0.5-5 pulses/m 2 to limit costs, these densities are still good enough for many applications, such as creating digital terrain models and conducting forest management inventories.
Overall, our findings indicate that the nationwide 5 pulses/ m 2 ALS data in Finland are adequate for retention tree detection, although higher pulse densities would be beneficial, especially for dead trees. Based on our methodology, creating large-scale maps of retention trees with a reasonable accuracy should now be achievable, and be able to guide landscape-level conservation and forest management. Additionally, our methodology could be applied for other purposes such as mapping carbon stocks more accurately, as more retained trees increase carbon sequestration (Santaniello et al., 2017). Note, however, that an accurate automatic methodology to obtain stand outlines will be required for the creation of such maps which often utilizes fused ALS data and CIR imagery (Dechesne et al., 2017;Leppänen et al., 2008). With the advent of such maps, research on retention trees might also shift gear. Research on retention trees has primarily been aimed at local stand scales while the landscape scale has been widely ignored (Gustafsson et al., 2020;Rosenvald & Lõhmus, 2008). Our methodology could facilitate research on larger scales by reducing constraints on time and costs of the associated larger data acquisition. For example, with assistance from our methodology, we might uncover the role of retention trees in altering landscape connectivity and species dispersal. Multiple publications have in recent decades claimed that retention trees should improve landscape connectivity and facilitate species dispersal (e.g. Vandekerkhove et al., 2013) but conclusive proof was (Rosenvald & Lõhmus, 2008) and is, based on our knowledge, still lacking. Another option that

ACK N OWLED G EM ENT
We thank Aki Suvanto for his assistance with receiving all the required remotely sensed data for our analyses and providing clarifications.

FU N D I N G I N FO R M ATI O N
This study was financially supported by the Academy of Finland

CO N FLI C T O F I NTE R E S T
None of the authors have any conflict of interest to report.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/2041-210X.13995.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available from the Dryad Digital Repository https://doi. org/10.5061/dryad.fqz61 2jw3 (Hardenbol et al., 2022).