Individual tree detection using template matching of multiple rasters derived from multispectral airborne laser scanning data

ABSTRACT

Multispectral airborne laser scanning (MS-ALS) provides information about 3D structure as well as the intensity of the reflected light and is a promising technique for acquiring forest information. Data from MS-ALS have been used for tree species classification and tree health evaluation. This paper investigates its potential for individual tree detection (ITD) when using intensity as an additional metric. To this end, rasters of height, point density, vegetation ratio, and intensity at three wavelengths were used for template matching to detect individual trees. Optimal combinations of metrics were identified for ITD in plots with different levels of canopy complexity. The F-scores for detection by template matching ranged from 0.94 to 0.73, depending on the choice of template derivation and raster generalization methods. Using intensity and point density as metrics instead of height increased the F-scores by up to 14% for the plots with the most understorey trees.

1. Introduction

Individual tree detection (ITD) can be used to derive forest information from airborne laser scanning (ALS) data. Existing ITD methods are generally based on two-dimensional (2D) or three-dimensional (3D) canopy segmentation. 2D segmentation involves image processing of rasters using a canopy height model (CHM) or a normalized digital surface model (nDSM). Conversely, 3D segmentation extracts information directly from the laser point cloud. Methods used for 2D segmentation include watershed segmentation (Tang et al., 2007; Ene et al., 2012), region growing (Hyyppä et al., 2001; Solberg et al., 2006), and neural networks (Kestur et al. 2018). Methods used for 3D segmentation include k-means clustering (Lindberg et al. 2013), mean shift (Xiao et al. 2019), and adaptive and agglomerative clustering (Lee et al., 2010; Gupta et al., 2010). Existing 2D and 3D methods typically only use the spatial information in the point cloud; the intensity data are rarely exploited.

To compensate for the lack of spectral information provided by traditional one-channel laser scanning, multispectral laser scanning systems (MS-ALS) have been developed in recent years. Such systems include multiple scanners, enabling recording of the intensity of reflected laser pulses at more than one wavelength when scanning (Junttila et al. 2018, 2017). For example, the commercial Optech Titan (Teledyne Optech, Canada) system scans at wavelengths of 532 nm, 1064 nm, and 1550 nm (Wichmann et al. 2015; Yu et al. 2017; Ahokas et al. 2016; Axelsson, Lindberg, and Olsson 2018; Matikainen et al. 2017). MS-ALS has been applied in land use classification and land cover mapping (Matikainen et al. 2017; Morsy et al. 2016; Wichmann et al. 2015), plant physiological and health evaluation (Junttila et al. 2018, 2017; Wallace et al. 2014), and tree species classification (Yu et al. 2017; Axelsson, Lindberg, and Olsson 2018; Budei et al. 2018). However, the applications of MS-ALS reported to date have involved detecting plants’ chlorophyll and moisture content or distinguishing between artificial and natural targets (Gong et al. 2015) rather than segmentation of targets such as individual trees.

Previous studies have used spectral information to improve ITD based on ALS data by applying combined analysis of multispectral images and ALS-derived layers (Leckie et al., 2003; Breidenbach et al., 2010; Ke et al., 2010), or by analysing multispectral images to check ALS segmentations and detect under-segmentation by recognizing different tree species (Heinzel and Koch 2012). The additional information provided by intensity measurements can complement spatial data, making it easier to distinguish tree crown edges. This could be particularly valuable in situations that are challenging for ITD, such as when trees of different species are located in close proximity and many understorey trees are present. While the use of spectral information by combining multispectral images with ALS data has been shown to enhance ITD, the usefulness of MS-ALS in this context remains to be determined.

Template matching was originally developed to detect treetops from aerial images by measuring the similarity between the tree templates and targets (Pollock 1996; Larsen and Rudemo 1998). The technique has since been improved and compared with other methods for the ITD (Larsen et al. 2011) used to separate tree crowns of multiple trees (Leckie et al. 2016). After the introduction of laser scanning in forestry, template matching methods suitable for use with laser data were developed. For example, Pirotti (2010) extracted 65% of all trees present when detecting individual trees from a CHM, while Nyström et al. (2014) detected 38% of windthrown trees from an object height model (OHM) created by subtraction from a digital elevation model (DEM). A similar approach has been used to detect fallen logs by analysis of high-resolution unmanned aerial vehicle (UAV) images (Panagiotidis et al. 2019).

Previous studies have demonstrated the effectiveness of template matching in target detection when using spectral and spatial information separately, but its use with combined datasets remains unexplored. In addition, template matching can detect patterns and textures, enabling the use of diverse metrics in ITD (e.g. height, point density, the vegetation ratio, and intensity values for multiple channels). Using template matching, features representing spatial, texture, and spectral information can be fused from intermediate results at the pixel level.

In this study, MS-ALS data were used for ITD by template matching. The method was applied to sites with varying understorey vegetation densities to answer the following questions:

(1) How does the accuracy of ITD based on airborne MS-ALS data compare to that achieved using single-wavelength ALS?

(2) Is height the best metric for ITD?

(3) Are there differences in the tree maps generated based on intensity data from various channels, and can these differences be exploited in ITD?

(4) What metrics should be selected for specific forest structures?

To improve template matching performance, the method was also refined by developing:

(1) A processing flow for raster image generalization, which contained the compensation of local low values that benefited different metrics with different levels.

(2) A method of deriving templates from the target plots instead of field data or manual selection was proposed.

2. Materials

2.1. Study area

The study area is located in Remningstorp (Figure 1), Västra Götaland, Sweden (58°27′18.35″N, 13°39′8.03″E) and covers an area of 1602 ha. The forest in this area is mainly managed spruce (Pinus Sylvestris) and pine (Picea Abies) for wood production.

Figure 1. Study area and plot positions. Red polygon on the right is the laser-scanned area

2.2. Field data

A field inventory was conducted in 2014 in 40 circular plots with diameters of 80 m. Each plot was divided into 16 subplots (20 × 20 m, Figure 3), and the positions of the subplot centres were determined with an RTK-GPS (Trimble Inc., Sunnyvale, CA, USA). The position, diameter at breast height (DBH), and species of each tree in each subplot was recorded. Individual-tree positions were determined using a Postex system (Lämås 2010) and calibrated using TLS data acquired in 2016. The height of a subsample of trees in the plots was also measured and recorded.

In the time between the field inventory and laser scanning, selective cutting and clear cutting were performed in 28 of the total 40 plots; only data for the 12 unchanged plots are considered here. Descriptive statistics for these plots are presented in Table 1. Because detecting understorey trees is challenging, an index was created to quantify the proportion of small trees. Here, small trees are defined as trees with DBH smaller than T, where T is the mean of the largest and smallest measured DBH for the plot. The reason for creating the index is further discussed later related to the performance of the ITD algorithm. For presentational purposes, the plot numbers were ordered according to the proportion of small trees. DBH distributions and cumulative curves are presented in Figure 2.

Table 1. Basic information about the 12 plots, including mean DBH and height, stem density, proportion of pine trees, and proportion of small trees. SD refers to standard deviation

Plot no.	Mean DBH (cm)	SD DBH/Mean DBH	Mean height (m)	SD height/Mean height	Number of stems/ha	Proportion of spruce/pine/deciduous (%)	Proportion of small trees (%)
1	30.2	0.24	21.8	0.16	615	93/1/6	24
2	28.3	0.24	21.2	0.17	704	95/0/5	27
3	24.4	0.24	19.7	0.17	921	68/3/29	28
4	22.6	0.33	18.3	0.25	905	85/0/15	29
5	28.8	0.29	21.2	0.2	740	91/0/9	31
6	27.4	0.22	19.5	0.15	635	49/49/2	32
7	28.2	0.29	20.9	0.21	543	93/0/7	35
8	28.3	0.36	20.7	0.25	820	92/0/8	38
9	21.5	0.38	17.6	0.28	899	77/1/22	44
10	22.2	0.35	17.9	0.25	808	49/0/51	46
11	17.1	0.42	14.7	0.33	1068	86/0/14	60
12	20.2	0.69	14.0	0.51	800	52/40/8	60

Figure 2. The DBH distribution of plots 1 − 12

The precision of the tree positions in each subplot relative to the subplot centres was high when using a Postex ultrasound system (Lämås 2010), but the accuracy of the subplot centres position had potentially larger errors due to the errors in the RTK-GPS measurements below the canopy (Olofsson, Lindberg, and Holmgren 2008). Consequently, all tree positions in the same subplot had an identical shift relative to the ALS data. These shifts were detected and eliminated by visual interpretation of the tree positions relative to the local maxima of the ALS data. A shift vector was added to the coordinates of all trees in each subplot, as shown in Figure 3.

Figure 3. Corrections of the single-tree positions in a representative plot. (a) Measured tree positions. (b) Corrected tree positions. The background rasters represent the output of the smoothed nDSM. The black solid lines indicate the plot’s boundary. The dotted lines show the boundaries of the 20 m × 20 m subplots

2.3. MS-ALS data

The study area was scanned on the 21 July 2016 with an Optech Titan (Teledyne Optech, Canada) commercial multispectral airborne laser scanner. This instrument scans at three laser wavelengths: 1550 nm (channel 1; short-wave infrared), 1064 nm (channel 2; near-infrared), and 532 nm (channel 3; green). To collect point cloud data with high resolution and high intensity, scanning was conducted at an altitude of around 400 m above ground level, generating 10 returns per m² in each channel, and up to four returns for each emitted pulse.

3. Methods

The framework of the ITD method and the study purpose are shown in Figure 4.

Figure 4. Framework of the ITD method tested here and questions considered when analysing the results

3.1. Preprocessing of laser data and raster derivation

The MS-ALS data were cut out with a 5 m buffer around the boundary of each 80 m diameter plot. The height was normalized to the ground and points lower than 2 m above the ground were removed using LAStools (2007–2020, rapidlasso GmbH, Germany) with the default settings. Range correction was done by multiplying the intensity by the squared distance to the scanner for each point.

Six raster images with 0.25 m resolution were derived from the point clouds based on the maximum height (H), point density (D), vegetation ratio (V), and average intensity from three channels (I₁, I₂, and I₃). To minimize the number of missing data cells in the raster images, the point clouds from the three channels were combined when calculating the maximum height, point density and vegetation ratio. Points less than 2 m above ground level were excluded from the point clouds when computing metrics, and the vegetation ratio was calculated as the ratio of the number of points 2 m or more above ground level to the total number of points.

3.2. Individual tree detection

3.2.1. Image generalization

Before template matching, raster cells with missing data and local minima were filled in the raster images and then smoothed to improve performance. Because of the chosen flight and scanning directions, there were missing data and minimal values between adjacent scan lines. Figure 5a depicts the scanning process and shows how the raster images are orientated with respect to the flight and scanning directions, as well as the missing data and minimal value stripes in a representative raster. To prevent the missing data and minimal values from affecting template matching, missing data cells were filled by computing the mean of all data-containing neighbouring cells in a 3 × 3 grid centred on the cell with missing data. Cells with values lower than their transverse two neighbour were replaced by the average of the two neighbours. Finally, the raster images were smoothed using a 3 × 3 mean filter.

Figure 5. The flight and scanning direction of the MS-ALS data used in this study and the generated point cloud. (a) Flight and scanning direction (i) shown on a maximum height raster (ii). (b) The corresponding point cloud seen from different directions

3.2.2. Template derivation

Template raster images with trees were derived directly from the rasters without using field measurements as reference data. According to the principle that a local height maximum has a high probability of being a tree top, a number of local maxima were defined as the centres of template images. For the plot size considered in this study, an n × n grid was created and the positions of the nodes were set as the initial seeds. The seeds were repeatedly moved to the adjacent raster cell with the highest value in a 2 m × 2 m neighbourhood until their positions no longer changed. Seeds whose initial placements corresponded to a gap in the canopy were excluded. Finally, 4 m × 4 m template images were defined based on the metric images of the plot, using the retained seeds as the centre (see Figure 6 for an example). During this process, it was possible for multiple initial seeds to converge onto the same local maximum, so the final number of templates could be lower than the initial one. The number of grid nodes used to create the initial seeds was varied from 1 to 15 × 15 to determine how such variation affected detection accuracy.

Figure 6. Templates (9 × 9) on a height raster image. (a) Initial positions of the nodes. (b) Final positions of the templates. Red boxes delineate the borders of the templates and blue spots indicate their centres. The background rasters are the M_DSM processed after the image generation steps proposed in.Section 3.2.1

3.2.3. Template matching

For each 90 m × 90 m plot, a similarity map with 0.25 m resolution was created using a moving window to illustrate the probability of each pixel being a tree position based on the templates (Figure 7). For a given pixel ${\mathrm{I}}_{i,j}$ in a plot image $\mathrm{I}$ a 4 m × 4 m subimage ${\mathbf{S}}_{i,j}$ centred on ${\mathrm{I}}_{i,j}$ was compared to each of the templates ${\mathbf{T}}_1$, ${\mathbf{T}}_2$ … ${\mathbf{T}}_m$ generated from a fixed number of nodes. The similarity was then calculated using Equations (1)(1) ${\mathbf{D}}_m\left(i,j\right)=\sum _{a=1}^A\sum _{b=1}^B{\left[{\mathbf{S}}_{i,j}\left(a,b\right)-{\mathbf{T}}_m\left(a,b\right)\right]}^2$(1) and (2(2) $\mathbf{P}\left(i,j\right)=1-{\left\{{\mathbf{D}}_m\left(i,j\right)\right\}}_{\mathrm{m}\mathrm{i}\mathrm{n}},m=1,2,\dots ,n$(2) )

(1) ${\mathbf{D}}_m\left(i,j\right)=\sum _{a=1}^A\sum _{b=1}^B{\left[{\mathbf{S}}_{i,j}\left(a,b\right)-{\mathbf{T}}_m\left(a,b\right)\right]}^2$(1)

(2) $\mathbf{P}\left(i,j\right)=1-{\left\{{\mathbf{D}}_m\left(i,j\right)\right\}}_{\mathrm{m}\mathrm{i}\mathrm{n}},m=1,2,\dots ,n$(2)

Figure 7. The similarity calculation process

where ${\mathbf{D}}_m\left(i,j\right)$ is the difference between subimage ${\mathbf{S}}_{i,j}$ and template ${\mathbf{T}}_m$, ${\left\{{\mathbf{D}}_m\left(i,j\right)\right\}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimum difference for all $n$ templates, and $\mathbf{P}\left(i,j\right)$ is the similarity of the pixel $\left(i,j\right)$ to a tree position. $\mathbf{S}$ and $\mathbf{T}$ were normalized against their maximum values, respectively, before calculating differences.

Using the same method, similarity maps (Figure 8) were derived from H, D, V, I₁, I₂, and I₃. To synthesize the similarities for different metrics, new similarity maps were created by averaging different combinations of H, D, V, I₁, I₂ and I₃. The averaged combinations are denoted P with used variables (H, D, V, I₁, I₂, and I₃) marked as subscript, e.g. P_{H, D} represents that the similarity maps from H and D are averaged for the template matching.

Figure 8. Similarly maps derived from different metrics by template matching. (a – e) H, D, V, I₁, I₂ and I_3. Red and blue indicate high and low values, respectively. Histogram equalization was performed to improve the visualization of the results

Finally, tree positions were detected from the similarity maps. The similarity maps were smoothed using a 3 × 3 Gaussian function, and a 3 × 3 moving window was used to detect the local maxima. For comparative purposes, tree positions were also detected directly from the nDSM generated using the preprocessing method described in section 3.1. The resulting tree positions are indicated by the label M_DSM in the results section.

3.3. Validation

To validate the detection accuracy, tree positions derived from ALS (ALS-trees) were matched to the field measurements (field-trees). The ALS-trees were first linked to all possible field-trees within a certain distance, then only the linked pairs with the smallest distances were retained for each tree. For trees with DBH values below 25 cm, the maximum linking distance was 3 m horizontally, while for larger trees the maximum distance was 12 × DBH (Olofsson, Lindberg, and Holmgren 2008). After matching, only trees within 35 m of the plot centres were used to assess the detection accuracy to reduce the statistical impact of trees close to the plot boundary. Precision, recall, and F-scores were calculated using Equations (3)(3) $P=\frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{d}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{b}\mathrm{y}\mathrm{A}\mathrm{L}\mathrm{S}}$(3) – (5(5) $F=\frac{2\times P\times R}{P+R}$(5) ):

(3) $P=\frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{d}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{b}\mathrm{y}\mathrm{A}\mathrm{L}\mathrm{S}}$(3)

(4) $R=\frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{y}}$(4)

(5) $F=\frac{2\times P\times R}{P+R}$(5)

In this study, F-scores were used for parameter comparison and selection in Section 4.1 and 4.2 because it is a comprehensive evaluation of the commission and omission effects. In Section 4.3 where presenting the detection accuracy, detection rate (R_D, same value as R in Eq. 4)(4) $R=\frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}}{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{s}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{y}}$(4) was used for better cross-interpretation of results from different studies. Specially, the detection rate of the small trees (R_D,s) was also presented to show the performance of understorey detection.

4. Results

4.1. Generalization of raster images

Figure 9 shows representative rasters generated during each of the generalization steps based on data for Plot 1. For height rasters, replacing missing data values had a clear effect in that it eliminated striped texture from the images. However, the striped texture in the density and vegetation ratio rasters only disappeared after replacing low values. In the case of the intensity rasters (particularly those for channels 1 and 2), it was difficult to distinguish tree crowns in the original images, and all generalization steps improved the detection of treetops and crown boundaries.

Figure 9. The effects of the processing steps on rasters based on different metrics within a 25 × 25 m area

The quantitative description of how each generalization step improved the detection performance was conducted using F-scores (Figure 10). Consistent with Figure 9, different metrics showed significant different improvements within each step. For M_DSM and H, F-scores increased most obviously when only replacing default values, while replacing the low values and mean filter did not contribute much. However, for density, only replacing default values showed little improvements compared to the raw data, while replacing low values was the most crucial step for all plots. For intensity, especially I₂ and I₃, F-scores improved significantly after every step, illustrating the necessity of conducting every step in the generalization.

Figure 10. F-scores for all plots based on individual metrics: M_DSM, H, V, D, I₁, I₂ and I_3.

4.2. Template derivation

The numbers of templates generated using different numbers of nodes, and the corresponding precision values, are given in Table A1 and Table A2 in Appendix. Higher numbers of nodes resulted in the generation of more templates and increased the probability of fake templates emerging (i.e. templates that could not be matched to field-measured trees). However, such fake templates did not reduce the detection accuracy, as shown in Figure 11. For example, even when one in three templates was fake for Plot 10, the F-scores were not appreciably lower than in cases without fake templates. However, additional templates did not improve the detection accuracy; broadly speaking, the detection was insensitive to the number of templates.

Figure 11. F-scores for detection using different numbers of templates

4.3. Detection accuracy

The detection accuracy values presented in this section are based on rasters obtained after applying all generalization steps and a 9 × 9 grid for template derivation. All results for individual and combined metrics for the 12 plots are presented in Figure 12. For all plots, template matching from H improved accuracy relative to detecting local maxima directly from the generalized nDSM. The effects of combining metrics on accuracy differed between plots. For plots 1–7 (see blue lines in Figure 12), the F-scores were highest using combinations containing H and particularly low for combinations including D. For plots 11–12 (green lines), performance was best for combinations including D but not H. Plots 8–10 (red lines) behaved similarly to plots 1–7 but yielded better performance when template matching was done using I instead of H or D.

Figure 12. F-scores using different metrics/combined metrics for 12 plots

The different performances of H, D and I for plots with different structures could be explained by the raster textures (Figures 13 and 14). Figure 13 presents an example of a small tree standing beside a dominant tree in Plot 12, which had 60% small trees. The H raster (Figure 13b) shows significantly lower values for the small tree compared to the surroundings, leading to difficulty of detection. However, in the rasters of density and intensity, the pixel values of small trees were influenced only slightly by the dominant tree. This example illustrates that the rasters of density and intensity provided more possibilities for detecting small trees beside dominant trees.

Figure 13. An example for a small tree standing beside a dominant tree in Plot 12. (a) 3D point cloud from laser data. (b) rasters of (i) H, (ii) D, (iii) I₂ and (iv) I₃. Red asterisks and blue circles indicate field-measured trees and detected trees, respectively; blue lines represent the links between them. Yellow cycle makes out the small tree in different rasters

Figure 14. An example for several understorey trees standing under a dominant tree in Plot 12. (a) 3D point cloud from laser data. (b) rasters of (i) H, (ii) D, (iii) I₂ and (iv) I₃. Red asterisks and blue circles indicate field-measured trees and detected trees, respectively; blue lines represent the links between them

For the cases of understorey trees, Figure 14 shows an example of how D and I improved the detection rate. The returns from the understorey trees changed D and I, causing more fragmentized texture compared to H, improving the probability of detecting understorey trees. Instead of only focusing on the top layer as H, D and I were more sensitive to the point cloud from different vertical layers.

Figure 15 shows the similarity maps for the height (H), intensity of the 1064 nm channel (I₂), and point density (D) for Plot 12. The similarity map for H looks similar to the nDSM, only showing the tallest layer of the canopy. The similarity map for I₂ exhibits higher values connecting the tree crowns around small trees, and the similarity map for D exhibits greater variation in the regions containing small trees, resulting in more local maxima. Although D provided more information that could be used to detect understorey trees, it was less easy to generalize than H, as shown in Figure 11; the D rasters had too much texture inside the tree crowns (Figure 9), which also might lead to over-segmentation.

Figure 15. Similarity maps for (a) H, (b) I₂ and (c) D in Plot 12. Red and blue indicate high and low values, respectively; black points indicate the positions of trees according to field measurements

When combining the similarity maps from different metrics, the performance of detection improved with the average effects among different metrics. Figure 16 depicts the linkage of ALS-trees and field-trees for Plot 12; the yellow boxes highlight regions with notable differences. It is clear that P_{I2, I3} was better at detecting small trees standing at the edges of big tree crowns. However, for clusters of small trees under big tree crowns (highlighted using green boxes), intensity was insufficient for detection. In such cases, the density and vegetation ratio were excellent for detecting clusters of small trees.

Figure 16. Linkage between trees detected by template matching and field-measured trees using the (a) H, (b) P_{I2, I3} and (c) P_{D, V} metrics for Plot 12. Red asterisks and blue circles indicate field-measured trees and detected trees, respectively; blue lines represent the links between them. The background rasters are the nDSM processed after the image generation steps proposed in.Section 3.2.1

To identify metric combinations suitable for use in plots with different forest structures, we considered the combinations with the highest and second highest accuracy for each plot (collectively referred to as top accuracies). We then selected the metric combinations with the most top accuracy results for plots with similar proportions of small trees. Table 2 summarizes the accuracies achieved with different combinations and the optimal combinations selected for different plot types. When detailed information on the plot conditions is available, the metric combinations listed under Option 1 in Table 2 are recommended; this option provides five different metric selections for five kinds of plot conditions. Option two represents a simplified solution for three levels of canopy complexity. For comparative purposes, Table 2 also shows the accuracy achieved using local maxima detection, which is a commonly used ITD method.

Table 2. Summary of the plot attributes, optimal metric combinations for single tree detection in each plot, the accuracy achieved using these combinations, and the accuracy achieved by direct detection based on local maxima for each plot

Plot no.	Proportion of small trees (%)	Template matching								Local maxima detection
		Option 1				Option 2
		Metrics	F-scores	RD	RD,s	Metrics	F-scores	RD	RD,s	F-scores	RD	RD,s
1	24	H	0.94	0.90	0.82	PH, I2, I3	0.93	0.90	0.82	0.91	0.87	0.72
2	27	H	0.93	0.89	0.68	PH, I2, I3	0.92	0.89	0.72	0.91	0.86	0.61
3	28	H	0.88	0.83	0.67	PH, I2, I3	0.88	0.86	0.73	0.87	0.81	0.64
4	29	H	0.83	0.73	0.42	PH, I2, I3	0.83	0.72	0.42	0.83	0.73	0.39
5	31	H	0.93	0.89	0.69	PH, I2, I3	0.93	0.89	0.75	0.93	0.89	0.72
6	32	PH, I2, I3	0.90	0.88	0.83	PH, I2, I3	0.90	0.88	0.83	0.89	0.87	0.84
7	35	PH, I2, I3	0.91	0.89	0.77	PH, I2, I3	0.91	0.89	0.77	0.87	0.90	0.75
8	38	PV, I2, I3	0.85	0.79	0.60	PH, I2, I3	0.84	0.75	0.52	0.83	0.72	0.47
9	44	PV, I2, I3	0.85	0.77	0.65	PI2, I3	0.85	0.81	0.72	0.80	0.73	0.52
10	46	PV, I2, I3	0.85	0.84	0.74	PI2, I3	0.84	0.87	0.81	0.84	0.73	0.56
11	60	I3	0.78	0.69	0.59	I2	0.77	0.70	0.61	0.74	0.61	0.46
12	60	PD, V	0.73	0.68	0.56	I2	0.69	0.59	0.46	0.59	0.44	0.24

5. Discussion

In this study, the individual tree detection (ITD) performance of MS-ALS was tested in plots with varying forest structures. ITD performance was evaluated on the basis of six metrics (height, point density, vegetation ratio, and intensity from three channels) and various combinations of these metrics. Height is the metric most commonly used for ITD based on ALS data (Eysn et al. 2015; Sačkov, Kulla, and Bucha 2019; Aubry-Kientz et al. 2019) even though its limitations with respect to detecting understorey trees have been widely noted. The results presented here show that this deficiency can be overcome by using additional metrics. The inclusion of intensity substantially improved ITD performance for some plots, demonstrating the advantages of MS-ALS over single-wavelength ALS.

In this study, the best metric differed between the plots. Therefore, rules for metric selection were proposed based on an analysis of the plots’ forest structure. Two main groups of plots were identified: ITD based on H yielded the best performance for plots 1–7, which had a typical sigmoidal cumulative DBH distribution curve (shown in Figure 2), whereas for plots with more small trees, ITD based on intensity or density yielded better results. To provide quantitative guidance when selecting metrics for ITD, an index was defined based on the proportion of small trees in a plot. When this index was low, height was sufficient to identify most trees. Conversely, in plots with high values of this index, intensity should be included to improve the detection of small trees beside or beneath the crowns of large trees, and height should have less weight. When this index reached 60%, performance was improved by including the point density metric because of the high point density associated with small tree groups. The index used in this study was derived from the DBH, but indicators reflecting the vertical complexity of the canopy structure could also be derived from the ALS point cloud (Coops et al. 2016; Wilkes et al. 2016). Future studies could use such indices as criteria for selecting raster metrics for ITD.

The template matching method developed in this study involves generating templates from the target data. In other studies on ITD using template matching, templates were created using simulation models (Gulbe and Mednieks 2013; Holmgren and Lindberg 2013; Nyström et al. 2014; Pirotti 2010) or by selecting sample trees manually in the image (Gulbe and Mednieks 2013). Gulbe and Mednieks reported detection rates of 56% and 71%, respectively, using these two types of templates. In this study, the templates were automatically derived from the height rasters. Local height maxima have a rather high probability of corresponding to tree positions, as shown in Table A2. This approach enabled fully automated detection without requiring manual interaction. Another benefit of the consistent generalization strategy was that crown shape deviations in both templates and the target plot data were eliminated to the same degree, reducing sensitivity to template shapes and textures. Figure 11 shows that using generalized rasters made the detection accuracy sufficiently stable with respect to variation in the number of templates, eliminating the need to determine the optimal number of templates.

Generalization of rasters has been a common preprocessing step in ITD, while the step of compensation of local low values used in this study was newly proposed for template matching for ITD. Although the local low values contained certain texture information of individual crowns, the additional texture within crowns might disturb the detection of the crown edges, which were more crucial for ITD. In addition, the random local low values also decreased the generalized level of templates. In this study, similarity maps based on different metrics were averaged by applying the same weighting to all the included metrics. In the future, it may be desirable to optimize the weightings of the metrics based on the forest structure and to thoroughly assess the similarity between target and template rasters for height, point density, and intensity distribution.

6. Conclusion

In this study, ITD from MS-ALS data was conducted, resulting in a detection rate of 0.90 to 0.68 for different forest structures. The performance of ITD was compared for different metrics, including height, point density, vegetation ratio, and intensity from three channels. For homogeneous forest stands, height remained the best and most easily used metric for ITD, but metrics such as intensity and density yielded better performance than height for forests with complex structures. For the plot with the most understorey trees in this study, ITD based on point density and vegetation ratio yielded 0.32 higher detection rate of small trees than local maxima detection based on height.

In addition, a framework for ITD based on a template matching algorithm was developed. Deriving templates from the target plot instead of field data or manual selection was shown to be feasible. The image generalization flow presented here increased the robustness among different metrics. The compensation for the local low values was shown to be crucial for generalizing density and intensity metrics.

This study has confirmed that the additional spectral information provided by MS-ALS compared to single-wavelength ALS can significantly improve the detection of understorey trees. Using intensity metrics from laser data can be a basis for other ITD algorithms. Future studies can be conducted on forests with even more layers to explore how the intensity metrics respond to the canopy layers. Additionally, more ITD methods can be developed using the point density distribution, which showed annular textures for individual trees in this study. The pixels with abnormal increasing density in the annular textures have high potential to be used for detecting understorey trees.

Contributions

All authors conceived and designed the study, analyzed the results and wrote the manuscript; Langning Huo improved and implemented the algorithm.

Acknowledgements

The data from field inventory and the laser scanning were supported by the Hildur and Sven Wingquist Foundation for Forest Science Research. We would like to thank Dr Johan Holmgren, who contributed to the study’s design.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

This work was supported by the Hildur and Sven Wingquist Foundation.