www.mdpi.com/journal/remotesensing Advances in Forest Inventory Using Airborne Laser Scanning

We present two improvements for laser-based forest inventory. The first improvement is based on using last pulse data for tree detection. When trees overlap, the surface model between the trees corresponding to the first pulse stays high, whereas the corresponding model from the last pulse results in a drop in elevation, due to its better penetration between the trees. This drop in elevation can be used for separating trees. In a test carried out in Evo, Southern Finland, we used 292 forests plots consisting of more than 5,500 trees and airborne laser scanning (ALS) data comprised of 12.7 emitted laser pulses per m2. With last pulse data, an improvement of 6% for individual tree detection was obtained when compared to using first pulse data. The improvement increased with an increasing number of stems per plot and with decreasing diameter breast height (DBH). The results confirm that there is also substantial information for tree detection in last pulse data. The second improvement is based on the use of individual tree-based features in addition to the statistical point height metrics in area-based prediction of forest variables. The commonly-used ALS point height metrics and individual tree-based features were fused into the non-parametric estimation of forest variables. By using only four individual tree-based features, stem volume estimation improved when compared to the use of statistical point height metrics. For DBH estimation, the point height metrics and individual tree-based features complemented each other. Predictions were validated at plot level.


Introduction
Trees are important for the carbon balance of the Earth.Forests have great economic and ecological importance.In Finland, about 77% of the country's land area is forested, which is the highest percentage in Europe.International interest in biomass detection is strongly linked to forest health, photosynthetic activity and other processes related to the carbon cycle [1] and the variability of the climate.There is a growing need and constant shortage of data for improved forest monitoring (e.g., [2]).
Approaches aimed at obtaining forest and forestry data from airborne laser scanning (ALS) data have been divided into two groups [3]: (1) area-based approaches (ABAs) and ( 2) individual/ single-tree detection approaches (ITDs).ABA prediction of forest variables is based on the statistical dependency between the variables measured in the field and the predictor features derived from ALS data [4][5][6][7][8][9][10][11][12].The sample unit in the ABA is most often a grid cell, the size of which depends on the size of the field-measured training plot.Stand-level forest inventory results are aggregated by summing and weighting the grid-level predictions inside the stand.When using ITD techniques, individual trees are detected and tree-level variables, such as height and volume, are measured or predicted from the ALS data, i.e., the basic unit is an individual tree.Then the stand-level forest inventory results are aggregated by summing up the tree data.The ABA does not make use of the neighborhood data of laser returns, whereas in ITD [13][14][15][16][17][18][19] the neighborhood data are applied using pattern recognition methods to locate individual trees and to derive the physical features of individual trees, such as height, species, crown diameter and crown volume.On the other hand, ABAs are based on the height and density data acquired by ALS, which are highly correlated with the forest variables.Currently, the ABA is operationally applied in the Nordic countries when carrying out standwise forest inventories.Some 5.5 Mha of Finland's forests have already been inventoried by applying ALS.
Initially, ITD began with the manual interpretation of analogue aerial images [20], followed by attempts to automate this task [21].In recent years, there have been several attempts to improve both image-based and laser-based ITD .An international comparison of the use of ALS for ITD was reported by Kaartinen and Hyyppä [34] and more recently by Vauhkonen et al. [63].In [34], it was concluded that (1) the ITD method is the main factor impacting on the achieved location, tree height and inventory accuracy, and that (2) the impact of an increase in point density is marginal when compared to the effect of the ITD method.More detailed comparisons of the extraction techniques were carried out in [68], and it was found that methods relying on pixel-based Canopy Height Models (CHMs) are not suitable for finding suppressed trees and that full-waveform or laser-point based techniques should be further investigated (as is also proposed in [42,69,70]).
The ABA has been performed with statistical metrics, calculated from a point cloud, and ITD has been done using the dimensions of the detected trees.However, the methods have been converged for a long time (e.g., [24,26,35,41,56,71]). Already in [16,17], individual tree methodology was proposed and tested at stand level.Hyyppä et al. [24,26] proposed an individual tree-based technique for focusing directly on tree clusters.A similar concept was adapted later on by Breidenbach et al. [44].Since the work done by Villikka et al. [35], density-related and height-related features, derived from first and last pulse ALS data, have also been used in the estimation of individual tree variables.In the past, the ABA has relied on large, meticulously collected field data from experimental plots, whereas ITD has provided results either without any field date, or with some reference to it [3].
Conventionally, trees are detected from CHMs, which are interpolated from the point data [3,63].In the present study, we demonstrate that the use of last pulse data is beneficial for tree detection.Firstly, we demonstrate by means of a large, real-life data set that trees can be discriminated better by using last pulse data than by using first pulse data.Secondly, we show by means of examples that the use of individual tree features as predictors increases the accuracy of stem volume in area-based predictions, which suggests that individual tree-based features should be used in addition to point height metrics in the ABA.These two improvements are separately depicted.The use of last pulse is described for tree detection in Sections 2.4 (the method of verification) and 3.1 (results).The improvements in area-based predictions are described in Sections 2.5 (methods) and 3.2 (results).These two improvements are separate cases.

Test Site and Field Data
The boreal managed forest study area, 5 × 5 km in size, is situated in Evo, Southern Finland.A total of 5,532 trees from 292 plots, each with a radius of 10 m, were used.The plots were measured in 2009.The average stand size in the study area is slightly less than 1 ha.The terrain elevation varies between 125 m and 185 m above sea level.Scots Pine (Pinus sylvestris) and Norway Spruce (Picea abies) are the dominant tree species, contributing respectively 40% and 35% of the total volume.The percentage of deciduous trees is 25% of the total stem volume.The sampling of the field plots in our study was based on pre-stratification of the existing stand inventory data.All the trees in each plot having a DBH greater than 5 cm were recorded and tree height, DBH, the lower limit of the living crown, crown width, and tree species were determined.The plot-level data, summarized in Table 1, were obtained by averaging or summing the tree data.The tree locations were calculated using the coordinates of the plot centers and using the direction and distance of the trees relative to these centers.The plot centers were measured using a Trimble GEOXM 2005 Global Positioning System (GPS) device (Trimble Navigation Ltd., Sunnyvale, CA, USA), and the locations were post-processed with local base-station data.This resulted in an average error of approximately 0.6 m.The tree heights were measured using Vertex clinometers.The DBH was measured using steel calipers.

Airborne Laser Data Collection
The ALS data were collected in the summer of 2009 using a Leica ALS50-II system operating at a pulse rate of 150 kHz.The data were acquired using a flight altitude of 400 m with FOV of 30 degrees.ALS data comprised of 12.7 emitted laser pulses per m 2 , including overlaps of strips.A footprint of 8.8 cm in diameter, collected by a beam having a divergence of 0.22 mrad, was achieved.The system was configured to record multiple returns per pulse, i.e., the first of many returns, the last of many returns, an only return and intermediate returns.

Laser Data Pre-Processing
The ground points were classified using TerraScan and based on the method explained in [72].Strip adjustment was performed using TerraMatch.A DTM was then calculated using the ground points.Canopy heights were calculated by subtracting the corresponding ground elevation from the elevations of the laser returns.

Improved Tree Detection
Most of the current approaches for tree detection are based on finding trees from the CHM, which is calculated as a maximum of canopy height values within each raster cell.Thus, the CHM corresponds to the maximum canopy height of the first pulse data.Our approach is based on using the canopy penetration capability of the last pulse returns with overlapping trees.When trees overlap, the surface model corresponding to the first pulse stays high, whereas with last pulse, even a small gap results in a drop in elevation, i.e., the trees can be more readily discriminated.Use of first pulse works, when the whole laser beam penetrates the gap between the crowns, so that the drop is detectable after filtering.With last pulse, the drop in elevation is substantially larger and the drop can be detected even with overlapping trees since the last pulse is more sensitive to lower canopy levels.Since the objective was to demonstrate the usefulness of last pulse data for tree detection using the above-mentioned phenomenon, we used the same tree detection method and verification for both the first and last pulse data.From the last pulse data, we calculate three different products corresponding to the minimum of last returns, the mean of last returns and the maximum of last returns.The tree-detection algorithm in this test was based on finding the local maximum.Finding the local maximum is computationally cost-effective, and this is needed in operational inventories where data processing currently costs about €0.20/ha in Finland.Based on [68], finding the local maximum was one of the most useful methods based on individual tree extraction in the EuroSDR/ISPRS Tree Extraction test.
The comparison consisted of the following steps: 1. Different raster models with 0.5 m × 0.5 m pixel size were created for tree location.The created models were as follows: minimum of last returns (Lmin), maximum of last returns (Lmax), mean of last returns (Lmean), and first returns maximum (Fmax).2. The raster models were smoothed by means of a Gaussian filter.A 3 m × 3 m window size was selected in order to eliminate minor tree level fluctuations and to avoid the merging of overlapping trees.Given Finland's forest conditions, larger window sizes lead to the merging of overlapping trees, especially when trying to locate trees from within the suppressed tree storey.3. Local maxima were sought from the smoothed surface model in a 3 m × 3 m window, and trees were considered to have been detected if the local maxima were greater than 2 m above the ground.
The extracted tree locations were compared with the tree locations measured in the field.The laser-detected trees were automatically matched with trees measured in the field, based on the Hausdorff distance on XY plane.The matching technique is described in detail in [73].Tree pairs with a distance from each other of less than 2.5 m were considered to be a correct match.The value of 2.5 m was set by taking the distance between the treetop and the roots into consideration.
The methodology and the applied automatic accuracy assessment are further demonstrated in Figure 1, which shows four raster models for one sample plot with the detected trees marked as '+' and field-measured trees as 'o'.The Fmax surface model corresponds to the commonly accepted way of finding trees.When comparing Lmin and Fmax models, it seems to be easier to discriminate trees from Lmin rather than from Fmax, even though visual processing of laser data in the forest is inferior to the best automated techniques [68].

Using Point Height Metrics and Individual Tree-Based Features in Area-Based Predictions
In order to demonstrate the usefulness of individual tree features in area-based predictions, which was the second objective of this paper, individual trees were located from the CHM (Fmax model), based on a method employing the minimum curvature object detection [66].This method was among the best in an extended analysis of the EuroSDR/ISPRS Tree Extraction test [68].Since the two objectives of this paper were processed in parallel, it was not possible to improve the minimum curvature object detection method with the last pulse data.The DBH was modelled using the tree height and crown area of each tree as the predictors [75].Tree volume was calculated based on Laasasenaho's volume equations [74] with laser-derived tree heights (the values of an unfiltered CHM at the tree locations) and DBH as the inputs.Thus, no local data were used to calibrate this conversion.
Plotwise mean height and the mean DBH value were obtained from the arithmetic mean values of the extracted individual tree heights and DBH.The volumes were obtained by summing up the individual tree volumes calculated from Laasasenaho's equations [74].Both individual tree-based and point height metrics (Table 2) were used as the inputs for the Random Forest (RF) classifier [66].RF is a non-parametric regression method in which the prediction is obtained by aggregating regression trees, each constructed using a different random sample of the training data, and choosing splits of the trees from among the subsets of the available features, randomly chosen at each node.The samples that are not used in training are called "out-of-bag" observations.They can be used to estimate the feature's importance (as applied in Section 3.2) by randomly permutating out-of-bag data across one feature at a time and then estimating the increment in error due to this permutation.The greater the increment, the more important the feature ( [66,76]).

Tree Detection Accuracy
Figure 2 depicts the percentage of correctly-matched trees using four different raster models for tree location.As mentioned earlier, the conventional way of detecting trees is to use the Fmax surface model (i.e., the CHM).The use of last pulse data gave a higher degree of discrimination between the trees than the use of first pulse data.The use of the raster corresponding to the minimum of last returns resulted in the highest discrimination between trees.An improvement of 6% in individual tree detection is better than that obtained by increasing the pulse density from 2 to 8 pulses per m 2 reported in [34].Figure 3 shows that the improvement in tree detection increases when the density of the forest stand increases.Obviously, the number of commission errors (the percentage of the number of falsely detected trees that would not be matched with field-measured trees) increases, as shown in Figure 4.The method also reveals many small saplings, which are not even recorded in field surveys because of their small diameter (>5 cm DBH is required to qualify for being recorded).For this reason, in future commission errors should be analyzed in detail, taking this into account.The proposed method is prone to creating real commission errors when there are gaps within individual tree crowns.Thus, the filtering has to be specified for the forest type.Figure 5 shows the percentage of correctly matched trees as a function of DBH.The smaller the DBH class, the better the performance of the proposed method (Lmin raster) compared to the use of the CHM (Fmax).With the DBH class 5-10 cm, the last pulse resulted in 10% better detection of trees.The results confirm that there is also substantial information for tree detection in last pulse data.Currently, in raster-based processing, this information has been largely neglected.The obtained results would even suggest the use of last pulse data for detection, but we assume that a hybrid model utilizing both the first and last pulse data should be developed, even when processing is done at raster level.The advantages of first pulse data obviously include the lower number of commission errors and the high quality of tree separation when the crowns are not overlapping, whereas the advantage of last pulse is in the separation of trees whose crowns overlap.A hybrid model, utilizing the advantages of both pulse types should be developed.Previously, last pulse data has been demonstrated to be usable for various applications.Liang et al. [77] used a simple technique (the difference of the first and last pulse returns under leaf-off conditions) to discriminate between deciduous and coniferous trees at individual tree level.In Matikainen et al. [78], a comparison between first pulse and last pulse laser scanner data in building detection was carried out.Both visual and numerical quality evaluations showed that the correctness of the results improved when last pulse data were used instead of first pulse data.In building detection using leaf-off data, the clearest difference between the first pulse and last pulse results was the proportion of the area classified as trees, which was considerably larger when the first pulse DSM was used [78].With leaf-off data, more complex decision rules have to be developed for tree detection.With leaf-off data, the coniferous trees behave similarly to the leaf-on data, but the response of the deciduous trees can vary and needs to be studied in detail.A hybrid technique, utilizing both the first and last pulse data, may provide a working solution for deciduous canopies.
The applied 2.5 m maximum distance for tree matching is a possible error source when matching small trees.A variable distance, based on tree height, could be used in future.Matching using variable distance based on the DBH of the tree is reported in [79].

The Accuracy of Area-Based Prediction of Forest Variables
We wanted to show that the use of individual tree features as predictors increases the accuracy of stem volume in area-based predictions.Table 3 summarizes the results for mean height, mean DBH and stem volume when using (a) all features (both individual tree-based and point height metrics); (b) point height metrics and (c) individual tree-based features.With regard to volume estimation, there is a substantial improvement when using individual tree-based features as input for area-based prediction when compared to point height metrics.The results were repeated using two point densities and two DBH versus height models, and the combined data set, consisting of individual tree-based and point height metrics, provided the most accurate forest variable predictions in all cases.Since individual tree-based features can easily be added to the ABA method, the results propose a potential improvement in the accuracy of the estimation of forest variables in the ABA.Individual tree-based features improved the ABA's accuracy since they had very high correlation, e.g., with the reference stem volume.When calculating the importance of the features, most of the individual tree-based features were among the best features.Figure 6 shows that, when estimating mean height, the best laser-derived feature was the mean height derived by using the individual tree technique.When estimating DBH, the best laser-derived features were (1) mean canopy height and (2) penetration to the ground (as used in [4][5][6]80,81] and which was confirmed to correlate powerfully with standing volume and biomass in [82]); (3) mean tree height (derived from the extracted individual trees) and (4) mid percentiles.For the estimation of stem volume, the best laser-derived feature was the stem volume derived from extracted individual trees, followed by the basal area derived from extracted individual trees.It is possible to easily derive further laser point height metrics and individual tree-based features.The selected features are those reported quite early on, in [11][12][13].Especially for the estimation of DBH, more features should be extracted to improve accuracy.In the field of statistics it is well-known that the use of predictors with high powers of explanation yields better estimates.
Another matter of concern is pulse density.It is probable that individual tree-based features lose some of their explanatory powers when applied at lower pulse densities, and with point densities of about 1 point/m 2 , the improvement is assumed to be more modest.However, in [34], pulse density did not have a significant impact on the accuracy of individual tree-based results Several ALS inventory studies have been carried out in the same Evo area as the present study.In [83], the ABA stem volume estimation RMSE (relative stem volume estimation error, Root Mean Squared Error) was 27.1% using 282 field plots for training the k-NN method.In [67], the obtained RMSE for stem volume ranged from 24.8% to 27.2%, which is comparable to the 25.2% obtained in the present study with point height metrics as predictors.An improvement of RMSE value from about 25% to 20.3-20.4% can be considered to be substantial.Future studies should test the sensitivity of individual tree algorithms in tree finding and the sensitivity of point density for the estimation of forest variables.Since the overall estimation methods based on RF are robust, our preliminary understanding is that estimation is not as sensitive to the tree finding algorithm as has previously been reported in ITD literature.It should be borne in mind that the applied tree extraction method, namely the minimum curvature objector detection method [66], was among the best in the extented EuroSDR/ISPRS comparison [68].Yet, even that method can be assumed to be improved upon by incorporating last returns, as explained in this paper.We believe that the use of physically-strong features [84], together with statistically proven non-parametric estimation techniques, can also result in savings in practical forest inventory, even when using lower point density data (<1 emitted laser pulses per m 2 ).

Conclusions
This paper reports two improvements for laser-based forest inventory.The first improvement is based on using the last returns for discriminating between overlapping trees.Using last pulse data and the test, which included 292 forests plots and more than 5,500 trees, an improvement of 6% for individual tree detection was obtained when compared to using first pulse data.In the 5-10 cm diameter breast height class, the use of last pulse resulted in a 10% better detection of trees than when using first pulse data.The results confirm that there is substantial information for tree detection in last pulse data, which should not be neglected even when using raster-based processing.
The second improvement is based on the use of individual tree-based features, in addition to statistical point height metrics, in the area-based prediction of forest variables.By using individual tree-based features as the input in non-parametric estimation, the root mean squared error, when compared to solely point height metrics, were reduced from about 25% to 20% at plot level.Point height metrics and individual tree-based features complemented each other in basal area estimation.The results confirmed the high usability of individual tree level features in the area-based estimation of forest variables.

Figure 1 .
Figure 1.Four raster models for one example plot (radius 10 m) with the detected trees marked as '+' and the field-measured trees marked as 'o'.The trees designated by A, B, C and D were detectable on the Lmin image but not on the Fmax image.Lmin refers to the minimum of last returns, Lmax to the maximum of last returns, Lmean to the mean of last returns and Fmax to the first returns maximum.The return height is color coded.

Figure 2 .
Figure 2. The percentage of correctly matched trees when tested with 5,532 trees surveyed in the field.

Figure 3 .
Figure 3.The percentage of correctly matched trees as a function of plot density when tested with 5,532 trees surveyed in the field.The number of trees refers to a plot with a radius of 10 m.

Figure 4 .
Figure 4.The percentage of commission errors for different surface models when using the local maximum finding as the tree detection algorithm.

Figure 5 .
Figure 5.The percentage of correctly matched trees as a function of diameter breast height (DBH).

Figure 6 .
Figure6.The feature importance (in the left column) and a scatter plot of the predicted versus observed values (in the right column) for mean height, mean DBH and volume estimates based on all features.The order of the feature index is the same as in Table2(24-27 are the individual tree-based features).

Table 1 .
Descriptive statistics of the field plots accessed in the study.

Table 2 .
The features used in predicting the forest attributes.
LBBasal area of the plot, derived from the extracted DBH LV Volume of the plot, derived from the extracted DBH and height

Table 3 .
The bias, RMSE, and correlation coefficient (R) between the predicted and observed values were calculated as measures of the accuracy of area-based inventory based on three different feature sets.