Automatic Estimation of Tree Position and Stem Diameter Using a Moving Terrestrial Laser Scanner

Airborne laser scanning is now widely used for forest inventories. An essential part of inventory is a collection of field reference data including measurements of tree stem diameter at breast height (DBH). Traditionally this is acquired through manual measurements. The recent development of terrestrial laser scanning (TLS) systems in terms of capacity and weight have made these systems attractive tools for extracting DBH. Multiple TLS scans are often merged into a single point cloud before the information extraction. This technique requires good position and orientation accuracy for each scan location. In this study, we propose a novel method that can operate under a relatively coarse positioning and orientation solution. The method divides the laser measurements into limited time intervals determined by the laser scan rotation. Tree positions and DBH are then automatically extracted from each laser scan rotation. To improve tree identification, the estimated center points are subsequently processed by an iterative closest point algorithm. In a small reference data set from a single field plot consisting of 18 trees, it was found that 14 were automatically identified by this method. The estimated DBH had a mean differences of 0.9 cm and a root mean squared error of 1.5 cm. The proposed method enables fast and efficient data acquisition and a 250 m2 field plot was measured within 30 s.


Introduction
Forest inventories are of paramount importance for sustainable and effective management of forest resources.Airborne laser scanning (ALS) is now widely used in forest inventories [1].An essential part of forest inventories based on ALS is field reference data collection.Field data are typically acquired by manual measurements on circular sample plots of 200-400 m 2 in size distributed over the landscape in question according to statistically rigorous sampling principles [2].The use of terrestrial laser scanning (TLS) to acquire field reference data in these ALS-based forest inventories has been proposed.In a recent review study by Liang et al. [3] it is acknowledged that there is a potential for utilization of TLS in forest inventories.A TLS system uses a laser to measure distances in a regular pattern around the scanner.This information is used to create a dense point cloud covering the surroundings of the scanner position.A setup where the scanner position is fixed is often referred to as a static TLS.If the scanner moves during data capture, the system can be referred to as a kinematic TLS.
There are several studies providing methods for highly detailed description of single trees from single or multiple scanning positions in the field, often including all branches and leaves [4,5].Tree architecture modeling has been used for biomass estimation [6][7][8] and forest structure description [9] at the plot or single tree level.
The use of TLS for the acquisition of specific biophysical tree parameters has been investigated in several studies, including estimation on diameter at breast height (stem diameter at 1.3 m above ground level; DBH) and tree stem volume [10][11][12][13][14][15].A full or partial reconstruction of individual trees has also been investigated [5,[16][17][18][19].Direct estimation of single tree biomass from TLS data has been described in [6,8].Others have estimated biomass from TLS data through tree stem reconstruction [20] or by using related TLS-estimates with existing allometric models [15,21].
Yang et al. [9] used multiple scans in individual field plots and found good agreement between field reference DBH, tree heights, and number of trees and corresponding properties derived from the TLS data.Positioning field plots using global navigation satellite system (GNSS) in forests can be affected by unfavorable conditions for reception of satellite signals, and reliably obtaining an accurate geo-referenced position typically requires a survey-grade GNSS receiver [22].As an alternative to GNSS positioning of field plots and trees in the field, Hauglin et al. [23] demonstrated that the ALS data that exist in many forest inventory projects can be used together with TLS data for positioning purposes by matching ALS and TLS data and exploiting the inherent absolute positions of the ALS data.Both the use of multiple TLS scans and the use of a single scan-the latter typically acquired from the plot center-have been proposed.By using only a single scan, the data acquisition and post-processing procedures are greatly simplified and the acquisition time is reduced [24][25][26].The main challenge in using a single scan setup is to handle the tree shadows that will occur in the TLS data.This is acknowledged in previous studies [12,27].Depending on the tree density and factors such as the presence of undergrowth and low branches, a number of trees will be occluded from the view of the scanner and hence cannot be detected.Automated detection algorithms are in practice associated with omission and commission errors, which means that some trees visible from the scanner will go undetected and that some patterns in the data occasionally will cause the algorithm to detect trees that are actually not present.For these reasons, some of the trees within a field plot will often not be detected from a single TLS scan.Liang et al. [28] reviewed several previous studies and reported that 10-32% of the trees were reported to be occluded from the plot center, depending on tree density and plot size.In the same studies, an additional 4-33% of the visible trees were not detected, depending on the algorithm that was applied and other factors.Tree occlusion and handling of omission errors when using multiple scanning positions were the main objectives of the studies [29,30].Astrup et al. [31] tested several approaches to estimating a corrected stand volume, taking into account the omitted trees.They found that without a correction the volume was substantially underestimated using single-scan TLS data.When applying a statistically founded correction, they obtained results which were closer to the volume computed from field measurements, but differences could still be observed.This suggests that it could be beneficial to investigate other methods to handle or avoid tree occlusion.
Time consumption is an important factor when considering the use of TLS on forest field plots.Bauwens et al. [32] noted that a plot typically can be measured within 10 min with a single static scan setup.The use of multiple scan locations in order to reduce occlusion will substantially increase the time needed for each plot.
As an alternative to the use of multiple static scanner positions, the problem of occlusion can be eliminated by moving the scanner while scanning, typically referred to as mobile laser scanning.A few studies have described the application of this technique in forests.Forsman et al. [33] used a laser scanner and cameras mounted on a car to measure trees along a forest road.Bauwens et al. [32] used a wearable rig consisting of a TLS scanner and additional sensors such as a GNSS receiver and an inertial measurement unit (IMU), which can be referred to as a kinematic TLS system.Walking with this rig through the forest produced a point cloud from which positions and other single tree properties could be derived.Bauwens et al. [32] and Ryding et al. [34] used a hand-held laser scanner with an integrated IMU to produce similar point clouds by walking through a forest area with the scanner.In these and other studies, point cloud data are typically captured from multiple positions and merged prior to a feature extraction such as detection of tree stems.Accurate position and orientation are needed in order to merge multiple scans, and lack of accuracy will produce noisy point clouds due to errors when aligning data captured from different positions.In the present study, we propose a novel approach by which point data from as little as a single scanner rotation of a moving laser scanner is used for tree stem detection.In the current set-up of the method, a scanner rotation typically would take 0.05 s.This approach ensures that a spatially consistent point cloud is used for tree stem detection, and could form the basis for simultaneous localization and mapping (SLAM) using the detected trees.Thus, the aim of the current study was to develop and implement the proposed method and test it in a real forest environment by assessing the obtained accuracy of tree detection, tree position, and estimated DBH.

Study Area and Field Reference Data Collection
The study area is located in Gran municipality in southeastern Norway (60 • 27' N 10 • 33' E, 500 m above sea level).The forest in Gran is boreal and dominated by Norway spruce (Picea abies (L.) Karst.) and Scots pine (Pinus sylvestris L.).An existing circular sample plot of 250 m 2 in a pine-dominated forest stand was subjectively selected from a dataset collected as part of another research project in the area.The plot was measured in field in August 2015.At the plot, all trees with a DBH >4 cm were accurately positioned, using an SOKKIA SET5 total station with additional measurements to at least two known points.The position of the two known points was obtained by post-processed GNSS baselines with base station data obtained from the Norwegian Mapping Authority.The receiver logged satellite data for more than 30 min.For all positioned trees, the DBH and species were also registered.The individual tree stems were calipered for diameter only once (a single measurement).In Table 1 a summary of the trees in the sample plot is given.

ALS Data Acquisition
The present study took advantage of ALS data acquired as part of an operational stand-based forest inventory.The ALS data were acquired on June 2015 with a Leica ALS70 sensor mounted on a fixed-wing aircraft, and the mean point density of the acquired ALS data was five points per m 2 .The ALS echoes were classified into "ground" and "non-ground" points by the data provider using the TerraScan software (Version 015, Terrasolid: Kanavaranta, Finland) [35].From the ALS points classified as ground points, a triangular irregular network (TIN) was created to represent the terrain surface height.The vertical accuracy of the surface model was expected to be approximately 20-30 cm [36].Only the ALS echoes classified as ground points were used in the present study.

TLS Data Acquisition
TLS data were collected on the sample plot in April 2016 using an instrumentation integrating three primary components.The three hardware components used are illustrated in Figure 1.The VLP16 from Velodyne [37] was used as the terrestrial laser scanner.This instrument can be referred to as a low-cost laser scanner compared to traditional TLS systems.All laser measurements were timestamped with GPS time during data capture.The Velodyne VLP16 laser scanner has 16 individual laser beams rotating around the z-axis.The rotation speed is 5-20 rotations per second, referred to as scan rate.The vertical angle distance between each beam is 2 degrees.This gives a total field of view of 30 × 360 degrees.In this project, the pulse repetition rate was 300 kHz and the scan rate was 20 Hz.The data were stored and exported with VeloView v3.1.1.
Remote Sens. 2017, 9, 350 4 of 15 field of view of 30 × 360 degrees.In this project, the pulse repetition rate was 300 kHz and the scan rate was 20 Hz.The data were stored and exported with VeloView v3.1.1.The Applanix APX-15 UAV [38] system was used to capture GNSS and IMU data.Position and orientation were post processed with Applanix POSPAC UAV v7.2.The algorithm used an inertial aided differential GNSS technique [39] with a baseline length up to 30.5 km.The GNSS reference station was obtained from the Norwegian Mapping Authority.Figure 2a shows how the position dilution of precision (Pdop) varied during data capture.The average estimated root mean squared error (RMSE) for the position in the horizontal direction was 0.22 m and 0.32 m in the vertical direction.The estimated RMSE for the orientation in roll and pitch was 0.06 degrees, and 1.30 degrees in heading.Some studies have shown that the GNSS solution can be improved by using the digital elevation model obtained from high-resolution ALS data.The techniques are often referred to as terrain aided GNSS [40], but was not tested in this study.The Applanix APX-15 UAV [38] system was used to capture GNSS and IMU data.Position and orientation were post processed with Applanix POSPAC UAV v7.2.The algorithm used an inertial aided differential GNSS technique [39] with a baseline length up to 30.5 km.The GNSS reference station was obtained from the Norwegian Mapping Authority.Figure 2a shows how the position dilution of precision (Pdop) varied during data capture.The average estimated root mean squared error (RMSE) for the position in the horizontal direction was 0.22 m and 0.32 m in the vertical direction.The estimated RMSE for the orientation in roll and pitch was 0.06 degrees, and 1.30 degrees in heading.Some studies have shown that the GNSS solution can be improved by using the digital elevation model obtained from high-resolution ALS data.The techniques are often referred to as terrain aided GNSS [40], but was not tested in this study.The Applanix APX-15 UAV [38] system was used to capture GNSS and IMU data.Position and orientation were post processed with Applanix POSPAC UAV v7.2.The algorithm used an inertial aided differential GNSS technique [39] with a baseline length up to 30.5 km.The GNSS reference station was obtained from the Norwegian Mapping Authority.Figure 2a shows how the position dilution of precision (Pdop) varied during data capture.The average estimated root mean squared error (RMSE) for the position in the horizontal direction was 0.22 m and 0.32 m in the vertical direction.The estimated RMSE for the orientation in roll and pitch was 0.06 degrees, and 1.30 degrees in heading.Some studies have shown that the GNSS solution can be improved by using the digital elevation model obtained from high-resolution ALS data.The techniques are often referred to as terrain aided GNSS [40], but was not tested in this study.The estimated positions were transformed to the coordinate system called EUREF89 Norwegian Transversal Mercator projection, which has a mapping scale value close to 1.The same coordinate system was applied for the ALS data.All other processing was performed using Matlab [41].

Processing Method
The data processing was divided into nine different processing steps, which can be seen in Scheme 1.
Remote Sens. 2017, 9, 350 5 of 15 The estimated positions were transformed to the coordinate system called EUREF89 Norwegian Transversal Mercator projection, which has a mapping scale value close to 1.The same coordinate system was applied for the ALS data.All other processing was performed using Matlab [41].

Processing Method
The data processing was divided into nine different processing steps, which can be seen in Scheme 1.The first step in this processing method was to calculate position and orientation of the laser scanner, described in Section 2.3.The second step was to apply the post processed position and orientation to the laser measurement.The laser points were realized in a local coordinate system.The position and orientation together with the installation calibration parameters provided the information to transform the laser points from this local coordinate system into the global coordinate system described in Section 2.3.

Point Group Classification
Step number 3 in Scheme 1 starts a classification process.Figure 3 shows the laser points from one scan rotation before the classification were performed.All laser measurements were classified into "groups of tree stem points," "groups of ground points," and "not classified points" based on a distance increment analysis for each of the 16 laser beams, also called laser channels.This is a sort of spike landmark extraction, where we searched for points that stands out from the rest [42].The technique is often used in leg detecting algorithms in robotics [43], but it was also used for tree classification by Forsman et al. [33].In this method, all laser points were separated into different point groups where each laser channel was treated separately.Two neighboring points were placed into the same point group if the difference between the measured distances was below a given threshold.Based on the evaluation of the field data, the following classification parameters were used: Step 1 Post process position and orientation for the laser scanner Step 2 Apply position and orientation to the laser point cloud Step 3 Classify laser points into: 1. Groups of tree stem points 2. Groups of ground points

Not classified points
Step 4 1. Estimate tree diameter and center point 2. Improve position and orientation using the estimated center points in the ICP process Step 5 Correct for height offset due to poor GNSS condition Step 6 Convert the height system into height above ground level Step 7 Convert estimated tree stem diameters to tree stem diameters at breast height

Step 8
Clustering tree stem center points to center point group based on location and assign a unique tree identification for each center point group Step 9 Calculate average tree stem diameter and position for each center point group Scheme 1.Each step in the proposed method is illustrated in the scheme.IPC in step 4 means iterative closest point algorithm, see text in Section 3.3 for details.
The first step in this processing method was to calculate position and orientation of the laser scanner, described in Section 2.3.The second step was to apply the post processed position and orientation to the laser measurement.The laser points were realized in a local coordinate system.The position and orientation together with the installation calibration parameters provided the information to transform the laser points from this local coordinate system into the global coordinate system described in Section 2.3.

Point Group Classification
Step number 3 in Scheme 1 starts a classification process.Figure 3 shows the laser points from one scan rotation before the classification were performed.All laser measurements were classified into "groups of tree stem points," "groups of ground points," and "not classified points" based on a distance increment analysis for each of the 16 laser beams, also called laser channels.This is a sort of spike landmark extraction, where we searched for points that stands out from the rest [42].The technique is often used in leg detecting algorithms in robotics [43], but it was also used for tree classification by Forsman et al. [33].In this method, all laser points were separated into different point groups where each laser channel was treated separately.Two neighboring points were placed into the same point group if the difference between the measured distances was below a given threshold.Based on the evaluation of the field data, the following classification parameters were used:  A point group was classified as a potential group of tree stem points if the distance increment was less than 0.1 m within one laser channel (see Figure 4a).In addition, the horizontal distance between the first and last point in the group had to be smaller than the maximum tree diameter.The analysis was done individually for each of the 16 laser channels.For each point group, a center point and tree radius were estimated using a circle fit algorithm [44].It is necessary to have three measurements to be able to estimate a center point of the tree stem and a tree diameter.Based on the construction of the laser scanner we can estimate the theoretical maximum distance for a tree detection.Due to the precision and beam divergence, the minimum detectable tree diameter is set to 4 cm.The minimum detectable tree diameter (Minstem) depends on the distance to the scanner, and can be estimated with Min 2 distance sin 360 scan rate puls repetition/16 (1) The tree stems were assumed to be vertical and well defined, i.e., clearly visible from the scanner and not hidden behind obstacles such as branches and other trees.The tree stem center points established from the different laser channels were subsequently merged together.The vertical distance between each laser channel is two degrees and this was used to decide if two center points belonged to the same tree.To be classified as a tree, there must be five consecutive center points with a vertical distance of two degrees (see Figure 4b).The minimum free visible tree stem (Sstem) of a tree stem can be estimated as A point group was classified as a potential group of tree stem points if the distance increment was less than 0.1 m within one laser channel (see Figure 4a).In addition, the horizontal distance between the first and last point in the group had to be smaller than the maximum tree diameter.The analysis was done individually for each of the 16 laser channels.For each point group, a center point and tree radius were estimated using a circle fit algorithm [44].It is necessary to have three measurements to be able to estimate a center point of the tree stem and a tree diameter.Based on the construction of the laser scanner we can estimate the theoretical maximum distance for a tree detection.Due to the precision and beam divergence, the minimum detectable tree diameter is set to 4 cm.The minimum detectable tree diameter (Min stem ) depends on the distance to the scanner, and can be estimated with Min stem = 2 × distance × sin 360 × scan rate puls repetition/16 (1) The tree stems were assumed to be vertical and well defined, i.e., clearly visible from the scanner and not hidden behind obstacles such as branches and other trees.The tree stem center points established from the different laser channels were subsequently merged together.The vertical distance between each laser channel is two degrees and this was used to decide if two center points belonged to the same tree.To be classified as a tree, there must be five consecutive center points with a vertical distance of two degrees (see Figure 4b).The minimum free visible tree stem (S stem ) of a tree stem can be estimated as A group of points was classified as ground if the distance increment was smaller than 0.1 m and the horizontal distance was larger than the maximum tree diameter.Finally, a point group was classified as "not classified" if the distance increment was greater than 0.1 m or if the distance increment was below 0.1 m for less than three consecutive measurements.
Remote Sens. 2017, 9, 350 7 of 15 A group of points was classified as ground if the distance increment was smaller than 0.1 m and the horizontal distance was larger than the maximum tree diameter.Finally, a point group was classified as "not classified" if the distance increment was greater than 0.1 m or if the distance increment was below 0.1 m for less than three consecutive measurements.If the number of consecutive tree stem center points were five or more the center points were accepted as a tree and used in step 9 illustrated in Scheme 1.

Processing Method, Position, and Orientation Improvement
The post processed position and orientation solution described in Section 2.3 had a degraded result due to unfavorable conditions for reception of satellite signals.The proposed method is not very sensitive to position and orientation accuracy, but if the distances between the trees become small, tree identification might become unsuccessful.To avoid unsuccessful identification, it was decided to test scan matching.A version of the iterative closest point algorithm (ICP) called iteratively re-weighed least squares [45] was implemented by adding a Matlab function created by Bergström [46].In the processing, the robust Welsch criterion function was used.The function estimates a threedimension translation and rotation between two different point clouds.The first point cloud was locked and called model points.The other point cloud was fitted to the model points and called data points.For each scan rotation, a point cloud was created and a translation and rotation was estimated between the model points and the new data points.It was found that the ICP became more robust if just the estimated tree center points calculated for each scan rotation were entered.Since the amount of estimated tree center points was limited, we allowed the model points to be built up sequentially.After the new data points were transformed to fit the model points, the new transformed data points were added to the model points.In the end, the model points contained all estimated center points.The process is referred to as step 4 in Scheme 1.

Processing Method, Height Adjustment
To estimate DBH, the heights obtained from the scanner were adjusted.Due to the properties of the scanner and the acquisition, ground points were not present in the entire area for each scan rotation.The ground points from the ALS data were therefore utilized when calculating each center points height above ground and the conversion of the estimated diameter to DBH.The method involves three different height adjustments and was performed separately for each scan rotation.The estimated position accuracy was assumed to be sufficient in the horizontal direction to correspond to the ALS data.The method is further described below.

Height Adjustment due to Poor 268 GNSS Conditions
Due to the poor GNSS conditions under tree canopies, the TLS data were not accurately aligned in height with the ALS data.The points classified as "ground" from the ALS acquisition were Illustrates points from the same laser channel reflecting on a tree stem.The red points have all shorter observation distance increment than the given threshold and forming a group of tree stem points.These points are used to estimate a tree stem center point and diameter marked with black color.(b) Illustrates five different groups of tree stem points with the center point marked as a black point and the red points are the laser points.If the number of consecutive tree stem center points were five or more the center points were accepted as a tree and used in step 9 illustrated in Scheme 1.

Processing Method, Position, and Orientation Improvement
The post processed position and orientation solution described in Section 2.3 had a degraded result due to unfavorable conditions for reception of satellite signals.The proposed method is not very sensitive to position and orientation accuracy, but if the distances between the trees become small, tree identification might become unsuccessful.To avoid unsuccessful identification, it was decided to test scan matching.A version of the iterative closest point algorithm (ICP) called iteratively re-weighed least squares [45] was implemented by adding a Matlab function created by Bergström [46].In the processing, the robust Welsch criterion function was used.The function estimates a three-dimension translation and rotation between two different point clouds.The first point cloud was locked and called model points.The other point cloud was fitted to the model points and called data points.For each scan rotation, a point cloud was created and a translation and rotation was estimated between the model points and the new data points.It was found that the ICP became more robust if just the estimated tree center points calculated for each scan rotation were entered.Since the amount of estimated tree center points was limited, we allowed the model points to be built up sequentially.After the new data points were transformed to fit the model points, the new transformed data points were added to the model points.In the end, the model points contained all estimated center points.The process is referred to as step 4 in Scheme 1.

Processing Method, Height Adjustment
To estimate DBH, the heights obtained from the scanner were adjusted.Due to the properties of the scanner and the acquisition, ground points were not present in the entire area for each scan rotation.The ground points from the ALS data were therefore utilized when calculating each center points height above ground and the conversion of the estimated diameter to DBH.The method involves three different height adjustments and was performed separately for each scan rotation.The estimated position accuracy was assumed to be sufficient in the horizontal direction to correspond to the ALS data.The method is further described below.

Height Adjustment Due to Poor 268 GNSS Conditions
Due to the poor GNSS conditions under tree canopies, the TLS data were not accurately aligned in height with the ALS data.The points classified as "ground" from the ALS acquisition were therefore used to estimate a height offset between the ALS and TLS data.This operation is referred to as step 5 in Scheme 1.Since the Velodyne laser scanner has a vertical field of view of 30 degrees, some laser channels were pointing down and others were pointing more upwards.In Figure 5, it is shown that the laser channels pointing down had a higher reliability when estimating the height difference between the ALS ground and the terrestrial ground points.Based on the findings in Figure 5, the laser channels 0, 2, and 4 having vertical observation angles corresponding to −15, −13, and −11 degrees, respectively, in the laser scanner frame system were used in the height difference estimation.
All estimated tree stem center points were adjusted for the estimated height difference between the TLS and ALS data.
Remote Sens. 2017, 9, 350 8 of 15 therefore used to estimate a height offset between the ALS and TLS data.This operation is referred to as step 5 in Scheme 1.Since the Velodyne laser scanner has a vertical field of view of 30 degrees, some laser channels were pointing down and others were pointing more upwards.In Figure 5, it is shown that the laser channels pointing down had a higher reliability when estimating the height difference between the ALS ground and the terrestrial ground points.Based on the findings in Figure 5, the laser channels 0, 2, and 4 having vertical observation angles corresponding to −15, −13, and −11 degrees, respectively, in the laser scanner frame system were used in the height difference estimation.All estimated tree stem center points were adjusted for the estimated height difference between the TLS and ALS data.

Conversion into Height above Ground
After the height adjustment due to poor GNSS observations, the TLS data refer to the same coordinate and height system as the ALS point cloud.Normalized heights (height above ground) for the TLS data were then computed for all TLS points by subtracted their respective ALS TIN heights at the corresponding x and y coordinates, referred to as step 6 in Scheme 1.

Conversion of Estimated Tree Stem Diameters to Tree Stem Diameters at Breast Height
Most of the TLS points were obtained in the region from ground to 3 m above ground.To extract the diameter at breast height, a simple model that assumed that the stem diameter was reduced by 1 cm per m along the stem was used.This is referred to as step 7 in Scheme 1. Points obtained between 0 to 0.5 m were excluded from the processing due to the irregular taper of the lower stem section.

Clustering of Center Points
In step 8 in Scheme 1, the tree stem center points were clustered based on location and distance analysis.All points closer than 0.3 m to each other were clustered together into one center point

Conversion into Height above Ground
After the height adjustment due to poor GNSS observations, the TLS data refer to the same coordinate and height system as the ALS point cloud.Normalized heights (height above ground) for the TLS data were then computed for all TLS points by subtracted their respective ALS TIN heights at the corresponding x and y coordinates, referred to as step 6 in Scheme 1.

Conversion of Estimated Tree Stem Diameters to Tree Stem Diameters at Breast Height
Most of the TLS points were obtained in the region from ground to 3 m above ground.To extract the diameter at breast height, a simple model that assumed that the stem diameter was reduced by 1 cm per m along the stem was used.This is referred to as step 7 in Scheme 1. Points obtained between 0 to 0.5 m were excluded from the processing due to the irregular taper of the lower stem section.

Clustering of Center Points
In step 8 in Scheme 1, the tree stem center points were clustered based on location and distance analysis.All points closer than 0.3 m to each other were clustered together into one center point group.Each center point group represents a tree with a unique identification.In step 9, the average tree stem diameter and position for each center point group was calculated.

Evaluation
Finally, the estimated tree coordinates and DBH were compared to the corresponding field reference measurements.Based on these results, mean differences, RMSE, and RMSE% were calculated as and RMSE% = RMSE y r 100 where y i is the estimate, y ir the reference, y r is the mean reference value, and n the number of observations.The tree positon accuracy was evaluated using Equations ( 3) and ( 4).

Results
The proposed method was able to detect 14 of the 18 trees on the 250 m 2 plot.The actual data capture time consumption was 30 s, which did not include the startup time, shutdown time, and the walk in and out of the forest to reach the sample plot.Startup and shutdown can be done while walking to the area of interest.The actual walking path during data capture is shown in Figure 6a.A total of 9.5 million points were measured and 0.2 million of these were automatically classified as tree stems.Based on these points, 3249 center points with a corresponding diameter were estimated.The result is shown in Figure 6b.The center points were clustered to center point group based on location and given a unique tree identification.
Remote Sens. 2017, 9, 350 9 of 15 group.Each center point group represents a tree with a unique identification.In step 9, the average tree stem diameter and position for each center point group was calculated.

Evaluation
Finally, the estimated tree coordinates and DBH were compared to the corresponding field reference measurements.Based on these results, mean differences, RMSE, and RMSE% were calculated as and where is the estimate, the reference, is the mean reference value, and the number of observations.The tree positon accuracy was evaluated using Equations ( 3) and (4).

Results
The proposed method was able to detect 14 of the 18 trees on the 250 m 2 plot.The actual data capture time consumption was 30 s, which did not include the startup time, shutdown time, and the walk in and out of the forest to reach the sample plot.Startup and shutdown can be done while walking to the area of interest.The actual walking path during data capture is shown in Figure 6a.A total of 9.5 million points were measured and 0.2 million of these were automatically classified as tree stems.Based on these points, 3249 center points with a corresponding diameter were estimated.The result is shown in Figure 6b.The center points were clustered to center point group based on location and given a unique tree identification.Table 2 summarizes the result of the tree stem diameter estimation.The difference between the field measured DBH and the estimated DBH varied between −1 cm and 3 cm for the detected trees.The undetected trees were all small, with a field measured DBH <10 cm, which was smaller than the Table 2 summarizes the result of the tree stem diameter estimation.The difference between the field measured DBH and the estimated DBH varied between −1 cm and 3 cm for the detected trees.The undetected trees were all small, with a field measured DBH <10 cm, which was smaller than the minimum detectable tree diameter except for one of the four undetected trees.The minimum detectable tree diameter depends on the distance from the scanner to the tree, and is given for each tree in the rightmost column in Table 2. Figure 7 summarizes this minimum detectable tree stem diameter versus distance from the scanner to the tree.For example, if the scan rate is 20, the maximum distance to allow tree detection of a tree with 10 cm in diameter is 7 m. a Calculated based on the field measured tree positions and the scanner location.See Figure 7 for further details.
In Scheme 1, step 9, an average DBH was calculated for each center point group.This was used as the estimated DBH for each detected tree and was compared to the field reference DBH.The result is presented in Table 2 and was used to estimate the mean difference, RMSE, and RMSE%.The estimation was based on Equations (3) to (5).Table 3 summarize these findings.The mean difference between the field reference DBH and the estimated DBH was 0.9 cm, with an RMSE of 1.5 cm.

Number of Trees
Mean Difference (cm) RMSE (cm) RMSE% 14 of 18 0.9 1.5 7.5 An average position was calculated for each center point group.Equations ( 3) and (4) were then used to calculate the mean difference and RMSE.The result is summarized in Table 4.

Discussion
The results showed that small trees were difficult to extract with the proposed method.The main reason was that three measurements were required to be able to calculate the tree stem center point and diameter.Based on Equation (1), it is possible to estimate the minimum tree stem diameter observed from different distances.Figure 7 summarizes this minimum detectable tree stem diameter versus distance from the laser scanner to the tree.
Remote Sens. 2017, 9, 350 11 of 15 observed from different distances.Figure 7 summarizes this minimum detectable tree stem diameter versus distance from the laser scanner to the tree.In the practical test, there were four trees that were not detected (Table 2).For each of the nondetected trees, the maximum detectable distance was calculated and for three of the four trees, the observed distance was too large for the tree to be detected.The smaller trees would be detected if observed distance or the scan rate had been reduced (Figure 7).If the scan rate is reduced from 20 Hz to 10 Hz, the possible detectable tree diameter is reduced with the same ratio, which means that 17 of 18 trees should theoretically be detectable by reducing the scan rate to 10 Hz.For the tree with ID number 6, the possible detectable tree diameter was smaller than the actual diameter of the tree.This means that one should have expected a positive detection.The actual reason for this was not found, but branches or other vegetation blocking the visibility of the stem might be one explanation.
The classification algorithm requires five continuous center points along the stem to be classified as a potential tree.The main reason for this is to filter out points that were falsely classified as a tree stem.Figure 8 shows the relationship between observed distance and the required minimum free visible tree stem.The calculation is based on Equation (2).This requirement was not a big problem in the test area, but different forest conditions with respect to tree species and ages might give other results.Minimum free visible tree stem (m) Distance from scanner to tree (m) In the practical test, there were four trees that were not detected (Table 2).For each of the non-detected trees, the maximum detectable distance was calculated and for three of the four trees, the observed distance was too large for the tree to be detected.The smaller trees would be detected if observed distance or the scan rate had been reduced (Figure 7).If the scan rate is reduced from 20 Hz to 10 Hz, the possible detectable tree diameter is reduced with the same ratio, which means that 17 of 18 trees should theoretically be detectable by reducing the scan rate to 10 Hz.For the tree with ID number 6, the possible detectable tree diameter was smaller than the actual diameter of the tree.This means that one should have expected a positive detection.The actual reason for this was not found, but branches or other vegetation blocking the visibility of the stem might be one explanation.
The classification algorithm requires five continuous center points along the stem to be classified as a potential tree.The main reason for this is to filter out points that were falsely classified as a tree stem.Figure 8 shows the relationship between observed distance and the required minimum free visible tree stem.The calculation is based on Equation (2).This requirement was not a big problem in the test area, but different forest conditions with respect to tree species and ages might give other results.
The classification algorithm requires five continuous center points along the stem to be classified as a potential tree.The main reason for this is to filter out points that were falsely classified as a tree stem.Figure 8 shows the relationship between observed distance and the required minimum free visible tree stem.The calculation is based on Equation (2).This requirement was not a big problem in the test area, but different forest conditions with respect to tree species and ages might give other results.The position and orientation solution in this study were relative good, but not good enough to merge the different scan rotations into one common point cloud.The main reason for this was the poor conditions for GNSS signal reception under the tree canopy.
For many forest inventory applications, the tree position accuracy achieved in this study is sufficient.It is, however, possible to improve the accuracy of the positions.The positions can be improved by different methods:

•
Take actions to improve the GNSS accuracy, such as reduce the distance to the GNSS reference station;

•
Use the ALS point cloud to extract the tree positions and perform matching per scan rotation to estimate a 2-dimensional translation and rotation for each scan rotation (c.f.Hauglin et al. [47]); • Develop the proposed method into a three-dimensional SLAM system; • The trees are assumed to grow vertically.In cases with low orientation accuracy, the trees can be used to estimate an average initial roll and pitch angle.
In the study presented by Bauwens et al. [32], three different tests were performed.The study area was located in Belgium and consisted of 10 different locations.The locations were chosen to maximize the variation in forest types, tree density, and terrain slope.One test was called FARO1 where the reference field was measured from one location with a Faro focus 3D 120.Another test called FARO5 was performed by measuring with the same instrument from five different locations.Finally, a test called ZEB1 was performed by collecting data with a handheld scanner (ZEB1).In terms of accuracy, the result achieved in the current study was between the FARO1 and the FARO5 tests (Table 5).
Table 5. Result for estimated DBH in the current study compared to Bauwens et al. [32].

Study Mean Difference (cm) RMSE (cm) RMSE%
Current study 0.9 1.5 7.5 FARO1 [32] −1.2 3.7 13.4 FARO5 [32] −0.2 1.3 4.7 ZEB1 [32] −0.1 1.1 4.1 An important factor when collecting data in the field is time consumption.The method with the smallest time consumption achieved by Bauwens et al. [32] was the FARO1 method, which required 10 min of field work.The required time only included the acquisition of data to get the measurements into a local coordinate system.In the current study, the TLS acquisition took 30 s.The latter also included the acquisition of GNSS and IMU data, which enabled us to produce tree positions in a global reference system.To be able to use the data in operational forest inventory, it is necessary to have the measurements in the same global coordinate system as the ALS data.It is, however, difficult to compare the methods directly because the forest conditions were different with respect to forest type, tree density, terrain conditions, and field plot size.In the current study, we used a plot with a radius of 8.9 m while Bauwens et al. [32] used a radius of 15 m.

Conclusions
GNSS signal reception conditions are typically degraded under forest canopies, and the proposed method worked with a relatively coarse positioning and orientation solution.This was possible because the laser measurements were divided into limited time frames determined by the scan rotation.For each scan rotation, tree center points and diameters were calculated.This provided a means to circumvent the need for a high accuracy position and orientation the moving scanner platform.To avoid the different trees from being mixed, the estimated center points were treated by an ICP algorithm to improve the homogeneity in the merged dataset.This process improved the tree identification process and made the entire process more robust.Future work should be devoted to improving the position and orientation solution with different GNSS techniques, different scan matching technics, or SLAM.Although the method proposed in the current study is in an initial phase of development, it has the potential to become a robust method, rapid and cost-effective acquisition of forest field plot data.

Figure 1 .
Figure 1.The complete measurement unit with global navigation satellite system (GNSS) antenna, Velodyne VLP16 laser scanner, and the Applanix APX-15 UAV sensor with an inertial measurement unit (IMU) at the bottom.The instrumentation was mounted and carried on a backpack during data capture in the field.

Figure 2 .Figure 1 .
Figure 2. The position dilution of precision (Pdop) during data capture is shown in (a) and (b) shows the number of satellites in the same time period.
of 30 × 360 degrees.In this project, the pulse repetition rate was 300 kHz and the scan rate was 20 Hz.The data were stored and exported with VeloView v3.1.1.

Figure 1 .
Figure 1.The complete measurement unit with global navigation satellite system (GNSS) antenna, Velodyne VLP16 laser scanner, and the Applanix APX-15 UAV sensor with an inertial measurement unit (IMU) at the bottom.The instrumentation was mounted and carried on a backpack during data capture in the field.

Figure 2 .Figure 2 .
Figure 2. The position dilution of precision (Pdop) during data capture is shown in (a) and (b) shows the number of satellites in the same time period.

Scheme 1 .
Scheme 1.Each step in the proposed method is illustrated in the scheme.IPC in step 4 means iterative closest point algorithm, see text in Section 3.3 for details.

• 15 •
Distance measurement increment threshold was set to 0.1 m; • Maximum tree diameter was set to 0.5 m;•Trees were assumed to be vertical.Based on these three classification parameters, all measurements were automatically classified into the three mentioned classes.No manual editing of the point classification was performed.Remote Sens. 2017, 9, 350 6 of Distance measurement increment threshold was set to 0.1 m; • Maximum tree diameter was set to 0.5 m; • Trees were assumed to be vertical.Based on these three classification parameters, all measurements were automatically classified into the three mentioned classes.No manual editing of the point classification was performed.

Figure 3 .
Figure 3. Approximately 10,000 laser points are collected during one scan rotation.Figure 3 shows an example of one scan rotation from an oblique perspective where the points are colored by elevation.The symmetric vertical structures are tree stems.

Figure 3 .
Figure 3. Approximately 10,000 laser points are collected during one scan rotation.Figure 3 shows an example of one scan rotation from an oblique perspective where the points are colored by elevation.The symmetric vertical structures are tree stems.

Figure 4 .
Figure 4.Each of the 16 laser channels was treated separately.(a) Illustrates points from the same laser channel reflecting on a tree stem.The red points have all shorter observation distance increment than the given threshold and forming a group of tree stem points.These points are used to estimate a tree stem center point and diameter marked with black color.(b) Illustrates five different groups of tree stem points with the center point marked as a black point and the red points are the laser points.If the number of consecutive tree stem center points were five or more the center points were accepted as a tree and used in step 9 illustrated in Scheme 1.

Figure 4 .
Figure 4.Each of the 16 laser channels was treated separately.(a) Illustrates points from the same laser channel reflecting on a tree stem.The red points have all shorter observation distance increment than the given threshold and forming a group of tree stem points.These points are used to estimate a tree stem center point and diameter marked with black color.(b) Illustrates five different groups of tree stem points with the center point marked as a black point and the red points are the laser points.If the number of consecutive tree stem center points were five or more the center points were accepted as a tree and used in step 9 illustrated in Scheme 1.

Figure 5 .
Figure5.The height difference between the airborne laser scanning (ALS) ground points and terrestrial laser scanning (TLS) ground points was compared and evaluated based on the 16 different laser channels of the Velodyne laser scanner.The blue bars give the height difference between the ground classified points from TLS data and ground classified points from the ALS data.The orange line tells which laser channels are associated with the corresponding blue bar.The figure shows that laser channels pointing down are more likely to hit the true ground level.

Figure 5 .
Figure5.The height difference between the airborne laser scanning (ALS) ground points and terrestrial laser scanning (TLS) ground points was compared and evaluated based on the 16 different laser channels of the Velodyne laser scanner.The blue bars give the height difference between the ground classified points from TLS data and ground classified points from the ALS data.The orange line tells which laser channels are associated with the corresponding blue bar.The figure shows that laser channels pointing down are more likely to hit the true ground level.

Figure 6 .
Figure 6.The walking path (purple line) during data capture and field measured position of the trees are shown in (a).In (b), all estimated center points are shown.Each tree stem center point is colored by a unique color for each tree.

Figure 6 .
Figure 6.The walking path (purple line) during data capture and field measured position of the trees are shown in (a).In (b), all estimated center points are shown.Each tree stem center point is colored by a unique color for each tree.

Figure 7 .
Figure7.Minimum diameter at breast height as a function of observation distance and scan rate calculated from Equation(1).Each diameter estimation requires minimum three laser points hitting the tree stem.Then minimum detectable diameter is set to 4 cm.

Figure 8 .
Figure 8. Calculated minimum unobstructed view to the tree stem based on observed distance.The

Figure 7 .
Figure7.Minimum diameter at breast height as a function of observation distance and scan rate calculated from Equation(1).Each diameter estimation requires minimum three laser points hitting the tree stem.Then minimum detectable diameter is set to 4 cm.

Figure 8 .
Figure 8. Calculated minimum unobstructed view to the tree stem based on observed distance.The calculation is based on Equation (2).

Figure 8 .
Figure 8. Calculated minimum unobstructed view to the tree stem based on observed distance.The calculation is based on Equation (2).

Table 1 .
Summary of tree diameter at breast height (DBH).

Table 2 .
Summary at the individual tree level: Estimated DBH, the number of laser points used for detection.The minimum (Min.dist.) and maximum (Max dist.)observed distance to the tree, and field measured DBH.Minimum detectable tree diameter calculated for each tree.All distances are given in cm.

Table 3 .
Comparison of estimated DBH and field reference DBH.Mean difference and root mean squared error (RMSE).

Table 4 .
Comparison of estimated tree positions and field reference positions.