Non Destructive Method for Biomass Prediction Combining TLS Derived Tree Volume and Wood Density

Abstract

This paper presents a method for predicting the above ground leafless biomass of trees in a non destructive way. We utilize terrestrial laserscan data to predict the volume of the trees. Combining volume estimates with density measurements leads to biomass predictions. Thirty-six trees of three different species are analyzed: evergreen coniferous Pinus massoniana, evergreen broadleaved Erythrophleum fordii and leafless deciduous Quercus petraea. All scans include a large number of noise points; denoising procedures are presented in detail. Density values are considered to be a minor source of error in the method if applied to stem segments, as comparison to ground truth data reveals that prediction errors for the tree volumes are in accordance with biomass prediction errors. While tree compartments with a diameter larger than 10 cm can be modeled accurately, smaller ones, especially twigs with a diameter smaller than 4 cm, are often largely overestimated. Better prediction results could be achieved by applying a biomass expansion factor to the biomass of compartments with a diameter larger than 10 cm. With this second method the average prediction error for Q. petraea could be reduced from 33.84% overestimation to 3.56%. E. fordii results could also be improved reducing the average prediction error from −17.24% to −7.30%. Only P. massoniana results had a low prediction error of 2.75% utilizing the total TLS-estimated volume, which was not improved by the biomass expansion method (3.82%).

1. Introduction

Trees take up the greenhouse gas carbon dioxide (CO₂) from the atmosphere. Therefore forests are a sink of CO₂. A measure of this sink can be derived from inventories in both managed and unmanaged forests [1]. The carbon storage of a tree is directly correlated to its biomass. In inventories, the above ground biomass is often predicted from volume tables utilizing the diameter at 1.3 m height (DBH) and height as reference variable. The latter two forestry parameters are easy to obtain. Additionally, power functions, known as allometric models, use these parameters as predictor variables to estimate the biomass of trees [2]. Those functions have been published [3] or can be accessed from online databases [4].

Models like these are not always available and sometimes the prediction quality needs to be improved to take account of regional and site dependent growth patterns [5]. New models have to be collated in this case, relying traditionally on destructive weight measurements.

In the last two decades Light Detection and Ranging (LiDAR) has become an efficient tool for the retrieval of forestry parameters [6]. Terrestrial laserscanning (TLS) technologies have been used to estimate the above ground volume of trees. Combined usage of volume and density values, derived from analysis of drilling cores, enables the analyst nowadays to calculate the biomass of trees in a non destructive way.

1.1. Related Work

1.1.1. TLS Derived Volume Estimations

All methods presented in this subsection estimating the above ground volume of trees directly are considered to have the potential to predict biomass by including density analyses. Works focusing on the extraction of parameters like the DBH and height can also be used if accurate allometric models are available. In addition non forestry related tree modeling approaches are included to a minor extent, as these works can still have an impact in computational forestry if their accuracy is validated in future research.

Early approaches for the estimation of forestry parameters from terrestrial point clouds have already been carried out in 2003. Simonse et al. [7] first determined a digital terrain model (DTM) by selecting the lowest points from the input cloud and extracted according to the DTM a layer corresponding to the hight of 1.3 m. In this layer the authors fitted circles in stem sections with the application of a Hough transform. Ground truth comparison was applied to estimated tree positions. For most trees a difference of less than 20 cm to the reference data could been carried out. Comparison regarding the DBH resulted in a standard deviation of 2.8 cm.

In 2004, Thies et al. [8] used registered TLS data derived from three scan positions per tree to perform a more accurate stem modeling. Two stem of different tree species, namely Fagus sylvatica (beech) and Prunus avium (cherry), were modeled by fitting cylinders into the point cloud. The non-linear fitting method for estimating the cylinder parameters minimized the square sum of residuals to the cylinder surface. The beech stem was modeled with 75 overlapping cylinders with a length 0.5 m each. The cherry stem was modeled with 59 cylinders in total. Both stems could be modeled to the height of the first fork. A root mean square error (RMSE) of 1.7 cm of the fitted cylinders was calculated and also a ground truth comparison to DBH measurements was applied, showing a deviation of −1.3 cm and 0.6 cm for both trees respectively.

These publications represent two different trends in retrieving forestry parameters from TLS data applied during the last decade. Various works focus on single trees manually extracted to gain a high level of information. Others determine less accurate parameters, often just DBH and height, but cover a larger area of interest, e.g., hectare sized plots [5].

Plotwise detection methods for forestry parameters are well reported within the literature:

Four plots were scanned by Bienert et al. [9] in both single and multiple scan mode. A DTM was extracted by using density variation along the z-axis. By a clustering operation followed by a circle fitting routine trees were detected in a 10 cm thick slice 1.3 m above the DTM. Output data partially consisted of the relative coordinates of each tree. After isolating the detected trees, DBH and height were also estimated for each individual tree. A DBH comparison to manual collected reference data resulted in a standard deviation between 1.21 cm and 2.47 cm.

Moscal et al. [10] developed an algorithm called point cloud slicing (PCS) algorithm where DBH, height, basal area and volume of trees could be calculated. The data was collected in single scan mode and then transferred to the voxel domain. After generating a DTM, a slice at 1.3 m height was produced to estimate DBH by cylinder fitting applying the least squares method. Height was predicted as the difference between minimum and maximum z-coordinates and volume was obtained by summing up the volumes of all voxels. Validation was accomplished by comparison to single tree measurement of DBH, stem location and tree height. TLS methods captured 91.17% of the variation of DBH, but on average only 57.27% of tree height was predicted.

Concerning modeling individual trees at a higher level of detail, various works have been reported in literature:

Pfeifer et al. [11] were the first authors found in literature to fit cylinders representing both stem and branches into multiple scan mode point clouds by non-linear least squares fitting. Single trees were extracted and cylinders fitted by five parameter determination. As proposed by Thies et al. [8], cylinders have been shifted backward and forward to get an approximate position of the next cylinder. NLS methods need accurate starting parameters, while nowadays being reported to be more accurate than Hough transforms [12]. Stem and branches were only partially detected. Quality assessment was performed through visual inspection and also by the calculation of a RMSE of 1.8 cm between the cylinders and their allocated points.

Eysn et al. [13] modeled 120 trees in a 0.65 ha sized area using 34 scan positions. The individual scan positions were registered into a common coordinate system. Yet each scan was processed independently. Through the manual digitization of stem and branch axes in AUTOCAD software all scans were handled as a two dimensional range or intensity map. After this manual step an automatic computation of cylinder radii was applied. Validation was performed by investigating deviations between the model and their corresponding TLS point clouds applied on five randomly selected trees. The mean deviation was between 2 mm and 6 mm for these trees, with a standard deviation of 1 cm for one tree and 2 cm for the others.

Bayer et al. [14] used a manual skeletonization method on extracted point clouds of single trees in leaves. The tree skeleton served as a basis for further automatic computation of branch angles, branch length and branch bending. Furthermore, by using α-shapes [15] on points allocated to a branch section the authors could calculate the space requirement of the branches by summing up the volumes of all tetrahedrons being part of the α-shape. Comparison with ground truth data was not applied. This work is mentioned, as further research on α shapes of branches might result in deeper knowledge about the photosynthesis potential of trees which is related to the growth of biomass.

Dassot et al. [16] determined the skeletons of 42 trees manually with Polyworks software. After an automatic cylinder fitting routine utilizing 25 cm long skeleton segments as input data, wrongly fitted cylinders were identified and removed manually. Gaps in the resulting cylinder model were linearly interpolated. Comparison of the modeled volume to ground truth measurements revealed a relative error in the range of 10% for the stems and in the range of 30% for branches with a diameter larger than 7 cm. Buksch et al. [17,18,19] developed an automatic method to determine a tree skeleton using point cloud data. The authors generated an octree search structure [20] of the input points. By using the neighborhood information of the octree cells a graph was extracted and contained cycles in the graph were removed. The resulting graph represented the tree skeleton. Goodness of fit evaluation was performed by calculating and depicting distances between the points and the predicted skeleton.

Xu et al. [21] found a way to calculate tree models by applying the Dijkstra algorithm [22], which is normally used for the solution of shortest paths determination in way finding routines. Every point was linked with neighboring points resulting in a connected graph. Along this graph for every node the shortest path from a preselected root point was calculated with the Dijkstra algorithm. The lengths of these paths were quantized and clustered into bins. Then a skeleton was formed by connecting centroids of adjacent bins. This method was further applied by Livny et al. [23]. Both publications aim for computer vision modeled trees rather than on the estimation of forestry parameters. All three presented tree skeletonization methods work automatically and the skeletons can be used as input data in tree modeling algorithms relying on manual derived tree structures, i.e., [13,14,16].

Belton et al. [24] extracted evergreen trees manually from TLS raw data. By performing a principal component analysis (PCA) they were able to estimate parameters like the curvature of the surface represented by each points neighborhood. After clustering operations through a Gaussian mixture model [25] utilizing these parameters leaf points were removed from the input point cloud. By connecting ellipses fitted in horizontal slices the tree skeleton was extracted. This skeleton was used as input data for a cylinder fitting routine. The volume of the tree was calculated by summing up the cylinders’ volumes. The authors compared the results (74 m³) to the estimated volume of an allometric function (34 m³).

Raumonen et al. [26] performed a PCA on small subsets of the point cloud. The eigenvectors were used to automatically extract skeleton lines of the underlying branch segments. Neighboring subsets belonging to the same skeleton line were merged to subpoint clouds representing branch sections. Into these sections cylinders were modeled. Validation with an artificial point cloud with known volume revealed a complete modeling of the stem and up to 90% of branches. In nowadays literature the method is commonly referred to as quantitative structure modeling (QSM).

Kaasalainen S. et al. [27] applied the QSM method to both real and artificial time-series. Five scanning campaigns covering two complete growth periods were carried out to analyze the biomass change of one Acer platanoides. In a laboratory experiment a Populus tremula branch was scanned four times additionally. Each time sub-branches had been cut off before new time stamp data was collected. The difference in volume change between the model and reference data was 20–40 mL (12% of the estimated total branch volume). For the A. platanoides scans validation could not be performed at this high accuracy level, as no destructive reference data could be collected. An overall trend in both total branch length and total branch volume was detected.

Hackenberg et al. [5] used high quality TLS data derived from trees growing under controlled conditions. Their method utilized spheres to follow the complete branching structure of a tree. The spheres’ surfaces cut the tree point cloud, resulting in spatially unconnected sub-point-clouds representing the cross sectional areas of the branches. Into these sub point clouds circles where fitted and their radii were used as radii for preliminary fitted cylinders. These cylinders were enhanced with a non-linear least squares fitting routine. Comparison between the resulting models and the input point cloud revealed that the models did cover up to 99% of the tree with a fit quality in sub-millimeter accuracy.

Latest research reports that nowadays both presented approaches can be combined to build tree models including the thin branching structure on a larger scale:

Calders et al. [28] further developed the QSM method to work on plot level. A 40 m radius plot containing 75 Eucalyptus trees was covered with 5 scans. A semi automatic extraction of those trees allowed the modeling proposed by Raumonen et al. [26]. 65 tree models could be linked to ground truth data gained by destructive harvesting. Biomass of those trees was inferred from total volume estimates (QSM) and basic density. Comparison between ground truthed biomass and TLS derived biomass revealed a root mean squared error (RMSE) of 16.1%.

The following publications [6,29] give a comprehensive overview of articles concerning LiDAR-based computational forestry.

1.1.2. Non Remote Sensing Based Biomass Estimations

Fayolle et al. [30] obtained biomass from the destructive sampling of 138 trees in Cameroon. The trees belonged to 47 identified taxa. The authors fitted single-species as well as multi-species allometric models to the data. The models utilized DBH and density as independent variables. Results revealed that single-species models were not necessarily better than the multi-species models when those included wood density as a predictor variable.

Ketterings et al. [31] predicted parameters a and b of the allometric model biomass = a(DBH)^b from data that does not require destructive measurements. The parameter b was calculated from the site-specific relationship between height and DBH, a was estimated from the average wood density.

Henry et al. [32] determined biomass for trees with a DBH larger than 20 cm in a non destructive way. At one and two meter intervals diameter measurements were carried out. These measurements lead to volume estimations, which further lead to biomass predictions by including wood density in the modeling. According to their conclusions, biomass is highly influenced by wood density, which is supposed to be a site-specific parameter.

2. Rationale

The goal in this paper is to present a non destructive method to estimate the above-ground biomass of trees. The prediction is based on density measurements combined with volume assessments derived from TLS data. Instead of using high quality input data obtained under ideal conditions as in Hackenberg et al. [5] scans of lower quality are utilized here. The scans were affected by wind, precipitation and occlusion effects. We present scans of three different tree species:

• Quercus petraea, a deciduous broadleaved tree species;
• Erythrophleum fordii, an evergreen broadleaved tree species;
• Pinus massoniana, an evergreen coniferous tree species.

The analyses will include a comparison between ground-truthed and predicted biomass results for 36 trees in total. Various proposed methods in literature [24,26] do not include detailed ground-truthing on a large number of trees at his level. Instead of relying mainly on manually derived measurements [16,32] the cylinder fitting in a fully manually isolated and semi-manually denoised point cloud of a single tree is performed fully automatically. Further all used tree-parameters can be extracted from this model.

The tree models built automatically with the suggested method will include the thin branching structure of trees. To our knowledge the only other existing peer-reviewed method automatically achieving this accuracy is the one of Raumonen et al. [26].

3. Materials and Methods

3.1. Definitions.

The presented method for the estimation of above ground biomass of trees combines TLS derived tree volumes and densities values based on measurements with high-frequency densitometry on simulated drilling cores [33]. The estimated volumes include both wood and bark compartments. Therefore the density values are calibrated to basic density based on complete stem discs including wood and bark. Various definitions used in the method description are given in Table 1.

Table 1. Definitions.

Abbreviation	Name	Description	Calculated by
bc	biomass content	weighted mean value of all ratios between dry weight and fresh weight of all stem discs of a tree	$\frac{\mathit{dryweight}}{\mathit{freshweight}}$
densb	basic density ( $\frac{g}{cm}$	derived from a simulated increment core using high frequency densitometry	HF-densitometry calibrated [ 34]
densf	fresh density ( $\frac{g}{cm}$	weighted mean value of the fresh density values of all stem discs of a tree	$\frac{\mathit{freshweight}}{\mathit{freshvolume}}$
wX	fresh weight (kg)	the fresh weight in kilogram of all compartments of a tree, whose diameter is larger than X cm	directly measured
BEF	biomass expansion factor	value to estimate total biomass from partially measured biomass [ 35,36]	$\frac{w_0}{w_{10}}$
${\mathit{volume}}_{\mathrm{X}}^{\mathit{TLS}}$	TLS derived volume (l)	volume of all compartments of a tree, whose diameter is larger than X cm	directly measured (Simple Tree software)
${\mathit{volume}}_{\mathrm{X}}^{\mathit{weight}}$	volume (l)	weight derived volume of all compartments of a tree, whose diameter is larger than X cm	$\frac{w_X}{den{s}_f}$
biomassTLS	biomass (kg)	total TLS derived tree biomass	${\mathit{volume}}_0^{\mathit{TLS}}\ast {\mathit{dens}}_b$
biomassTLS*	biomass (kg)	total tree biomass estimated with biomass expansion factor derived from TLS	${\mathit{volume}}_0^{\mathit{TLS}}\ast {\mathit{dens}}_b*\mathit{BEF}$
biomassweight	biomass (kg)	total tree biomass derived from biomass content and fresh weight	w0 * bc

Table 1. Definitions.

Abbreviation	Name	Description	Calculated by
bc	biomass content	weighted mean value of all ratios between dry weight and fresh weight of all stem discs of a tree	$\frac{\mathit{dryweight}}{\mathit{freshweight}}$
densb	basic density ( $\frac{g}{cm}$	derived from a simulated increment core using high frequency densitometry	HF-densitometry calibrated [ 34]
densf	fresh density ( $\frac{g}{cm}$	weighted mean value of the fresh density values of all stem discs of a tree	$\frac{\mathit{freshweight}}{\mathit{freshvolume}}$
wX	fresh weight (kg)	the fresh weight in kilogram of all compartments of a tree, whose diameter is larger than X cm	directly measured
BEF	biomass expansion factor	value to estimate total biomass from partially measured biomass [ 35,36]	$\frac{w_0}{w_{10}}$
${\mathit{volume}}_{\mathrm{X}}^{\mathit{TLS}}$	TLS derived volume (l)	volume of all compartments of a tree, whose diameter is larger than X cm	directly measured (Simple Tree software)
${\mathit{volume}}_{\mathrm{X}}^{\mathit{weight}}$	volume (l)	weight derived volume of all compartments of a tree, whose diameter is larger than X cm	$\frac{w_X}{den{s}_f}$
biomassTLS	biomass (kg)	total TLS derived tree biomass	${\mathit{volume}}_0^{\mathit{TLS}}\ast {\mathit{dens}}_b$
biomassTLS*	biomass (kg)	total tree biomass estimated with biomass expansion factor derived from TLS	${\mathit{volume}}_0^{\mathit{TLS}}\ast {\mathit{dens}}_b*\mathit{BEF}$
biomassweight	biomass (kg)	total tree biomass derived from biomass content and fresh weight	w0 * bc

3.2. Data Collection

P.massoniana data were collected in March 2013 in Baiyun forest farm near Pingxiang city, Guangxi province in sub-tropical, southern China (22°6'2'' N; 106°45'01'' E). The research plot is located in a wind exposed area, which did affect the scan quality due to crown movement during scanning operations.

The laserscans were taken with the Z+F IMAGER 5010 [37] in scan mode high (10,000 pixel/360°). The Z+F IMAGER 5010 belongs to scanners using phase-shift measurements as the range measurement principle. Those scanners record only one return per emitted beam. In a crown canopy environment, where multiple returns are expected, this circumstance leads to a large number of wrong measurements, also called ghost points. Scan data taken with scanners based on pulsed time-of-flight measurements tends to contain less ghost points, as those scanners can correctly record multiple returns [29]. Nevertheless phase shift scanners have a higher data acquisition rate. Faster completed scans of the target tree reduce both operation costs and the likelihood of high wind velocities during the scanning process.

Four 400 m² sized subplots were scanned with a total of 19 scans each. In each subplot three trees were selected for destructive sampling to measure the fresh weight and to perform density analysis (Section 3.2 and Table 2).

Table 2. Overview of the applied sample and scan parameters used for data collection.

Species	Q. petraea	E. fordii	P. massoniana
Number of trees	12	12	12
Scan mode	superhigh	superhigh	high
Average scans per tree	8	8	6

Table 2. Overview of the applied sample and scan parameters used for data collection.

Species	Q. petraea	E. fordii	P. massoniana
Number of trees	12	12	12
Scan mode	superhigh	superhigh	high
Average scans per tree	8	8	6

In October 2013, twelve E. fordii were selected to be scanned with the Z+F IMAGER 5010c [38] in scan mode superhigh (20,000 pixel/360°) with eight scan positions each. The trees were located in Baiyun forest farm. Weight and density measurements were performed in the same manner as the P. massoniana measurements.

In southern Germany near the town of Pirmasens (49°17'05'' N; 7°39'15'' E) in February 2014, twelve Q.petraea were processed in the same way. The laserscans were affected by precipitation and wind.

Tree and scan characteristics having a great impact on the modeling quality of trees are summarized in Table 3. Tree parameters, i.e., DBH and height are given in the tables of Appendix A.

Table 3. Overview of species characteristics.

Species	Q. petraea	E. fordii	P. massoniana
leaf/needle condition	off (leaf)	on (leaf)	on (needle)
bark	rough	plain	rough
epicormic shoots	yes	no	no
wind	yes	yes	yes
precipitation	yes	no	no

Table 3. Overview of species characteristics.

Species	Q. petraea	E. fordii	P. massoniana
leaf/needle condition	off (leaf)	on (leaf)	on (needle)
bark	rough	plain	rough
epicormic shoots	yes	no	no
wind	yes	yes	yes
precipitation	yes	no	no

With Z+F LaserControl software [39] multiple scans were co-registered into a common coordinate system. Afterwards trees were extracted manually with the 3D-volume-selection filter from the total scan data and exported in an ASCII file format. Further denoising and building of the tree models for volume prediction was performed with a java software of simpletree project [40].

Ground Truth and Density Measurements

We use height as a proxy for standing length. This is defined as the length of the path along the cylinder axes from the base to a defined point.

After harvesting stem discs were taken at predefined absolute lengths (0.2 m, 1.3 m, 3 m, 5 m, 7.5 m, 10 m, 12.5 m, 15 m, 17.5 m and 20 m). Not all tree stems were analyzed up to a height of 20 m. The analysis stopped at lower heights, in case the trees were shorter (Appendix A), forked or shattered during the felling operations.

For each stem disc fresh and basic density were determined using standard methods. Fresh volumes (including bark) were calculated from thickness (calipers) and surface area (flatbed scanner). Fresh weight was determined in the field, dry weight after drying for 72 h at 104 °C [34]. In order to simulate nondestructive sampling, density values of hypothetical increment cores were simulated (dens_b). Radial density profiles were measured with high frequency densitometry [41] on each of the stem discs from breast height (1.3 m). The resulting relative density values were calibrated to basic density [34] based on complete stem discs including wood and bark. Of those density profiles either an arithmetic mean (E. fordii, P. massoniana) or a 1cm length segment at 2/3 of the radial distance from pith to bark (Q. petraea) was used as representative density for the tree [33,42]. The respective segments had been determined to ideally represent the mean density of tree stems. This is a possible source of error in the calculation of branch biomass.

Additionally, segments (here a segment is the stem part between two discs) were weighed. Weighed branches and segments were classified into four classes containing all compartments with a diameter larger than 10 cm, larger than 7 cm, larger than 4 cm and the total combined above ground wood compartments respectively.

Compartments with a diameter larger than 7 cm will be referred to as solid wood, the ones with a diameter up to 7 cm as twig wood. Fresh density and weight were used to predict the volume of single segments and the volume of complete diameter classes, here-on referred to as weight predicted volume.

3.3. TLS Volume Estimation—Using Software Simpletree

3.3.1. Cylinder Creation

The principal idea of cutting sphere surfaces with tree point clouds [5] is utilized in the program simpletree [40]. The implemented method starts at the stem base of the tree and fits a circle into a small slice with the lowest z-coordinates of the point cloud P. The circle is enlarged and transformed in 3d into a sphere. The sphere surface is cut with P, extracting a sub point cloud S located in the ϵ-neighborhood of the sphere. Spatial clustering of this sub point cloud subdivides S into i clusters S_i. Each S_i represents a cross sectional area of the stem or a branch. Into each S_i a circle is fitted and the procedure is repeated. Whenever a new sphere is detected, its center point and radius are utilized as end point and radius of a cylinder. The start point of the cylinder is the sphere center point of the predecessor sphere of the processing chain. Then a tree model is generated by adding parent and child relations to the cylinders. Fit quality is improved with various methods, e.g., with non linear least squares fit. A detailed method description is available in Hackenberg et al. [5]. The authors processed the largest cross sectional area first, other clusters were processed later with threshold adjustment, as the order of the detected branches was expected to change during this step.

Based on Raumonen et al. [26] we adjusted the method by detecting stem and large branch points before the circle fitting routine. For each point the normal is computed by performing a principal component analysis of a point’s ϵ-neighborhood (ϵ = 2.5 cm), if the neighborhood contains at least 9 points. All points are taken into account which satisfy the following conditions:

• the angle between its normal and the z-axis is larger than 60°; points located on stems not growing in the direction of the z-axis are deleted by this criterium;
• the third principal component accounts to a maximum of 10% of the spatial variation, the second principal component to a minimum of 35%. Only points in a coplanar neighborhood satisfy both eigenvalue conditions [24].

The qualified points are clustered with Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm [43]. The largest cluster contains points located on the stem and the main branches. Not all points located on the stem and the main branches will pass the two thresholds tests utilized before the DBSCAN operation. To not loose this information every point in a 5 cm neighborhood of each pre-qualified point is also considered to be a stem point. Using only two different sets of thresholds during the cylinder detection, one for the detected stem points and one for the remaining points, was found out to be more robust than the threshold adjustment proposed in Hackenberg et al. [5].

3.3.2. Denoising

As the quality of scans is considered worse than those published in Hackenberg et al. [5], some additional denoising procedures were performed.

3.3.2.1. Selection of Scans

Before the extraction with Z+F LaserControl software [39] all scans were visually examined. Up to three wind affected scans were deleted. This minimized a well known source of error in the laser scans. Wind is the most troublesome factor since it decreases the quality of tree description, especially in the upper part of the trees [29]. A minimum of at least five remaining scans was considered to contain enough information to prevent occlusion effects.

3.3.2.2. Filter by Intensity

E. fordii and P. massoniana scans contained many points located on needles and leaves. The albedo of a leaf or needle differs from the albedo of woody biomass. Reflected by differing intensity values, this circumstance allows the removal of noise points by applying a threshold test for the intensity. Kukko et al. [44] showed, that intensity values decrease with an increasing incidence angle of the laser beam. Another experiment by Pesci et al. [45] gives a quantitative evaluation about the magnitude of intensity reduction dependent on the irregularity of the reflecting surface. The more planar the surface is, the higher is the decrease. In addition the distance of the object to the scanner affects the intensity of the reflected beam significantly [46]. Those side effects have not been taken mathematically into account, as a modeling approach would require the knowledge of many unknown factors. Thereby this denoising approach has to be performed manually with visual validation.

Figure 1 depicts three intensity histograms with different thresholds (125, 150 and 175) as well as the resulting point clouds of an E. fordii. The large distribution on the left side of the histograms is accounted to mainly noise/leaf points, the distribution on the right to woody biomass. The noise distribution has patterns of a normal distribution, while the biomass distribution has not. We assume that the previous mentioned influence of incident angle, distance to the scanner and surface structure is too strong and complex for a normal distribution. Also both distributions do overlap, therefore an automatically determined threshold, e.g., the local minimum between the two distributions, is not considered appropriate. The higher the threshold is set, the more noise is deleted, but with it also points located on the stem or branches. We manually selected the threshold with a slider. Operating the slider updates a visualization of the histogram and the resulting point cloud in real time. The used threshold was chosen in a manner that by visual inspection no effect on the stem points could be detected.

Figure 1. Filtering by intensity value threshold for E. fordii (E#10, height ≈ 19.3 m): (a) histogram of intensities, threshold 125; (b) all points with intensity below 125 are deleted; (c) histogram of intensities, threshold 150; (d) all points with intensity below 150 are deleted; (e) histogram of intensities, threshold 175; (f) all points with intensity below 175 are deleted.

Figure 1. Filtering by intensity value threshold for E. fordii (E#10, height ≈ 19.3 m): (a) histogram of intensities, threshold 125; (b) all points with intensity below 125 are deleted; (c) histogram of intensities, threshold 150; (d) all points with intensity below 150 are deleted; (e) histogram of intensities, threshold 175; (f) all points with intensity below 175 are deleted.

3.3.2.3. Filter by Neighborhood Information

Ghost points tend to be spatially isolated. Deleting all points whose nearest neighbor is further away than 5 mm deletes mainly ghost points. We implemented an Euclidean clustering operation [47] and removed the clusters containing only one point. This operation uses as input parameter only a maximum distance. A cluster is initially generated by a seed point. All points within a distance smaller than this input distance are also added. For every new added point also its neighbors are added. The cluster is full when for all added points their neighbors are already included and another cluster is created.

3.3.2.4. Isolating the Target Tree

The used point clouds contained parts of neighboring trees. An Euclidean clustering operation with a threshold of 10 cm often isolates the target tree’s clusters from neighboring trees’ clusters [5].

3.3.2.5. Inclusion of Artificial Gaps

If isolating the target tree did fail, artificial gaps were included manually as proposed in Hackenberg et al. [5]. All points contained in a user defined box are deleted from the point cloud. This procedure was also used to isolate the stem and delete branches if needed.

3.3.2.6. Buffering with Original Point Cloud

After the proposed denoising operations are performed, the density of the resulting point cloud is low. For every point in the low density point cloud also its neighbor points not further away than 5 cm are added to the cloud. This is similar to the euclidean clustering operation, but for newly added points the procedure is not repeated. The enriched point cloud contains only little noise in the spatial neighborhood of the stem and the branches.

3.3.2.7. Combination of Two Point Clouds

In P. massoniana scans stem and branches were treated separately. First a thinned out version of the target tree was produced. The stem was extracted manually from this cloud with the user selection box and buffered with the original point cloud. The stem in high density was merged with the complete tree in low density, duplicated points were deleted.

For the P. massoniana scans, we used the following denoising chain:

First an intensity filter was applied with threshold set in a manner that only a small number of points on woody biomass was deleted. Afterwards noise was removed by utilizing nearest neighbor distance information. By an euclidean clustering operation unconnected artefacts of other trees were removed. To prevent needle noise from being added to the cloud again in the following operation, the stem was isolated manually from the branches and stored in a second point cloud. This was performed by including artificial gaps and clustering afterwards. Only the stem point cloud was spatially enriched by buffering with the original point cloud. Afterwards enriched stem points were merged with the branch points again.

The Q. petraea scans were denoised with this chain:

First a filter was applied by utilizing the neighborhood distance informations. If necessary, artificial gaps to the neighboring trees were included. Then the remaining neighbor trees’ artefacts were removed by euclidean clustering.

For the E. fordii point clouds, the following chain was used:

Leaf points were deleted by applying a threshold for the intensity values. In addition a second filter for the maximum distance to the nearest neighbor was applied. Remaining leaf points were deleted together with neighboring tree artefacts with euclidean clustering. Only the main branching structure was left over. Smaller twigs were also removed, thereby the complete remaining point cloud was buffered with the original one without re-adding leaf noise.

3.3.3. Output

We produced output data consisting of volumes of the four diameter classes. In addition the volume of each stem segment was written in an extra file per tree.

4. Results

We performed consistent analyses on the three different tree species. The axis-scale of corresponding plots is set accordingly to those of Q. petraea. Biomass prediction results and used density values are summarized in appendix A.

4.1. Quercus Petraea

The Q. petraea point clouds contain a large number of ghost points (Figure 2(a)). Visual inspection as depicted in Figure 2(b) indicates a successful removal of the ghost points. The creation of cylinder models covering most parts of the tree (Figure 2(d)) is thereby enabled.

Figure 2. Processing the scan of a Q. petraea (Q#1, height ≈ 25 m): (a) input point cloud containing a large number of ghost points and compartments of other trees; (b) denoised point cloud; (c) detected stem and main branch points are colored brown; (d) the cylinder model of the tree; (e) the automatically extracted stem.

Figure 2. Processing the scan of a Q. petraea (Q#1, height ≈ 25 m): (a) input point cloud containing a large number of ghost points and compartments of other trees; (b) denoised point cloud; (c) detected stem and main branch points are colored brown; (d) the cylinder model of the tree; (e) the automatically extracted stem.

Visualization of the cylinder model revealed two different types of errors:

• points located on the upper part of the stem are not detected as stem points (Figure 2(c)), because the lower point density in combination with a higher curvature of the stem surface prevents a successful pass of the test for stem points; as the thresholds are then adjusted for the detection of smaller twigs, upper stem cylinders might have underestimated diameters (Figure 3);
• points located at the tips of branches often scatter due to windy conditions during the scanning campaign; this circumstance leads to overestimated cylinders, an example can be seen in Figure 4.

Figure 3. Close up view of a Q. petraea Q#8: (a) sub-point cloud of high quality in the mid of the stem; (b) the fitted stem in this area; (c) sub-point cloud of low quality in the crown; (d) erroneous stem cylinders with an underestimated diameter in this segment of the tree.

Figure 3. Close up view of a Q. petraea Q#8: (a) sub-point cloud of high quality in the mid of the stem; (b) the fitted stem in this area; (c) sub-point cloud of low quality in the crown; (d) erroneous stem cylinders with an underestimated diameter in this segment of the tree.

Figure 4. Impact of scan quality on the modeling of twig wood: (a) sub-point cloud of high quality of a P. avium presented in Hackenberg et al. [5]; (b) the modeled twigs; (c) sub-point cloud of low quality of Q. petraea Q#8; (d) erroneous twig cylinders with an overestimated diameter in this section of the tree.

Figure 4. Impact of scan quality on the modeling of twig wood: (a) sub-point cloud of high quality of a P. avium presented in Hackenberg et al. [5]; (b) the modeled twigs; (c) sub-point cloud of low quality of Q. petraea Q#8; (d) erroneous twig cylinders with an overestimated diameter in this section of the tree.

Comparison between weight predicted volume and TLS-derived volume reveals an increasing overestimation when taking into account the smaller diameter classes. ${\mathit{Volume}}_{10}^{\mathit{TLS}}$ is predicted with an average overestimation of 2.46% and a standard deviation of 11.38%. For the total above ground volume (${\mathit{volume}}_0^{\mathit{TLS}}$) the mean value for the error is 32.68% with a standard deviation of 25.53%. The prediction errors for all four diameter classes are visualized in a form of a box plot in Figure 5.

Figure 5. Box plots of the relative volume prediction error for Q. petraea. From left to right increasingly thinner branches are included in the analysis.

Figure 5. Box plots of the relative volume prediction error for Q. petraea. From left to right increasingly thinner branches are included in the analysis.

The results of a detailed analysis for the fit quality of stem segments can be seen in Figure 6. Up to a height of 10 m the TLS derived volume is slightly higher than the weight predicted volume, expressed by the small positive prediction errors (green graph). Above 10 m segments tend to be underestimated, with increasing height the underestimation gets larger. The last depicted error for the segments between 17.5 m and 20 m with a value of −55.93% is only derived from 2 segments, as for the 10 other Q. petraea ground truth analysis of the stem terminated at lower heights due to shattered crowns after the felling operation. The uncertainty of prediction quality increases above a height of 10 m, expressed by larger standard error bars.

Figure 6. Prediction quality analyses of the stem of 12 Q. petraea. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

Figure 6. Prediction quality analyses of the stem of 12 Q. petraea. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

The error regarding the total volume of the stem up to a certain height ranges from −3.87% up to 2.75%. The only two underestimation values are observed above a height of 15 m.

Figure 7 shows two linear models with biomass^weight as independent variable. The slope of the linear model utilizing biomass^TLS as the dependent variable indicates a large overestimation by the method. The computed average error for all 12 trees is 33.84%. The prediction error gets smaller with a value of 3.56% by using biomass^TLS∗ as the dependent variable (BEF = 1.14).

Figure 7. Two linear models utilizing both biomass^weight of 12 Q. petraea as independent and either the biomass^TLS (f(x) = 1.31x + 13.24, adj. r² = 0.46 ) or biomass^TLS∗ (f(x) = 0.96x + 35.03, adj. r² = 0.54) as dependent variable.

Figure 7. Two linear models utilizing both biomass^weight of 12 Q. petraea as independent and either the biomass^TLS (f(x) = 1.31x + 13.24, adj. r² = 0.46 ) or biomass^TLS∗ (f(x) = 0.96x + 35.03, adj. r² = 0.54) as dependent variable.

4.2. Erytrophleum Fordii

The noise removal for E. fordii could be performed efficiently. Leaf hits were removed without much manual interaction. A negative effect of the denoising procedure is the deletion of points representing the small twigs in the crown (Figure 8(b)). Still the main branching structure of solid wood could be preserved and was considered to be of high quality after visual inspection.

Into the remaining point cloud cylinder models without visually detectable over- or underestimation were fitted (Figure 8(c)). Comparing weight predicted volume estimators with TLS derived volumes showed an underestimation of −7.85% for ${\mathit{volume}}_{10}^{\mathit{TLS}}$. The standard deviation of the relative error for this class is 6.08%. Both average error and standard deviation seem to be stable when taking into account the finer branching structure up to all compartments with a diameter larger than 4 cm. Adding the finest branch class including twigs smaller than 4 cm diameter increases the underestimation to −17.47% and the standard deviation to 7.41% (Figure 9).

Figure 8. Input and output data of a E. fordii (E#1, height ≈ 14.6 m): (a) input point cloud containing other trees, ghost points and leaf hits; (b) input point cloud after denoising procedures; (c) the resulting cylinder model.

Figure 8. Input and output data of a E. fordii (E#1, height ≈ 14.6 m): (a) input point cloud containing other trees, ghost points and leaf hits; (b) input point cloud after denoising procedures; (c) the resulting cylinder model.

Figure 9. Box plots of the relative volume prediction error for E. fordii. From left to right increasingly thinner branches are included in the analysis.

Figure 9. Box plots of the relative volume prediction error for E. fordii. From left to right increasingly thinner branches are included in the analysis.

The segment analysis of Q. petraea was carried out in the same manner as that performed for E. fordii (Figure 10). Up to a height of 10 m an underestimation of −3.81% up to −7.20% was computed. The underestimation gets larger for higher segments, reaching an underestimation of −20.31% for the segments between 12.5 m and 15 m.

Figure 10. Prediction quality analysis of the stem of 12 E. fordii. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

Figure 10. Prediction quality analysis of the stem of 12 E. fordii. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

The analysis of the stem from its root up to a certain height revealed a prediction error between −5.22% and −7.20% up to the maximum analyzed height of 15 m.

A calculated prediction error of −17.24% for biomass^TLS is emphasized by the slope value (0.73) of the according linear model depicted in Figure 11. A lower underestimation of −7.30% was achieved by utilizing biomass^TLS∗ as dependent variable (BEF = 1.37) for the linear modeling.

Figure 11. Two linear models utilizing both biomass^weight of 12 E. fordii as independent and either the biomass^TLS (f(x) = 0.73x + 30.29, adj. r² = 0.95 ) or biomass^TLS∗ (f(x) = 0.85x + 22.49, adj. r² = 0.96) as dependent variable.

Figure 11. Two linear models utilizing both biomass^weight of 12 E. fordii as independent and either the biomass^TLS (f(x) = 0.73x + 30.29, adj. r² = 0.95 ) or biomass^TLS∗ (f(x) = 0.85x + 22.49, adj. r² = 0.96) as dependent variable.

4.3. Pinus Massoniana

The P. massoniana scans contained hits on the needles (Figure 12(a)). These noise points could only be partially deleted (Figure 12(b)).

Figure 12. Processing the scan of a P. massoniana (P#10, height ≈ 16.4 m): (a) input point cloud containing a large number of ghost points; (b) the same cloud with ghost points and needle hits partially deleted; (c) the resulting cylinder model; (d) the extracted stem model.

Figure 12. Processing the scan of a P. massoniana (P#10, height ≈ 16.4 m): (a) input point cloud containing a large number of ghost points; (b) the same cloud with ghost points and needle hits partially deleted; (c) the resulting cylinder model; (d) the extracted stem model.

Visually detectable errors in the fitted tree models (Figure 12(c),12(d)) occur mainly during the modeling of twigs (Figure 13):

• some twig parts are not modeled as the cylinder following ends in an occlusion gap; points containing essential information are lost in denoising operations (Figure 13(a),13(b));
• other areas contain still a large number of noise points; therefore cylinders are overestimated (Figure 13(c),13(d)).

Figure 13. Close up view of P. massoniana (P#10, height ≈ 16.4m): (a) a sub point cloud containing low amount of information of twig; (b) the cylinder model of the twig ends in an occlusion gap; (c) a sub point cloud still retaining a large number of noise points; (d) overestimated cylinders representing twigs.

Figure 13. Close up view of P. massoniana (P#10, height ≈ 16.4m): (a) a sub point cloud containing low amount of information of twig; (b) the cylinder model of the twig ends in an occlusion gap; (c) a sub point cloud still retaining a large number of noise points; (d) overestimated cylinders representing twigs.

Analyzing the volume prediction quality of different diameter classes shows an approximately stable overestimation of 3.87% for ${\mathit{volume}}_{10}^{\mathit{TLS}}$ and 3.27% for ${\mathit{volume}}_0^{\mathit{TLS}}$ (Figure 14). The standard devation reaches its minimum for predicting ${\mathit{volume}}_4^{\mathit{TLS}}$ and has the maximum of 11.79% when predicting ${\mathit{volume}}_0^{\mathit{TLS}}$.

Figure 14. Box plots of the relative volume prediction error for P. massoniana. From left to right increasingly thinner branches are included in the analysis.

Figure 14. Box plots of the relative volume prediction error for P. massoniana. From left to right increasingly thinner branches are included in the analysis.

The segment analysis could only be performed for 11 trees, as for one tree during ground truth measurement stem segments could not be related to the absolute lengths accordingly. The average prediction error for segments taken from heights below 10 m ranges from −2.18% up to 4.32% and then drops to −9.04% for the segments taken between 12.5 m and 15 m (refer to Figure 15). The analysis of the stem up to a certain height reveals errors ranging from −2.18% up to 2.35%.

Figure 15. Prediction quality analyses of the stem of 12 P. massoniana. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

Figure 15. Prediction quality analyses of the stem of 12 P. massoniana. The green graph describes the averaged prediction error for all segments taken from the same height. The brown graph represents the prediction error of the cumulative volume of all stem segments up to a certain height. Both graphs include the standard error in form of error bars.

Biomass estimators have an average relative error of 2.75% for the prediction according to biomassTLS . The prediction utilizing the expanded volume (BEF = 1.25) overestimates with an error of 3.82%. The worse prediction of the second method is also denoted by the smaller ${\mathrm{r}}_{\mathit{adj}}^2$ value of the linear model depicted in Figure 16.

Figure 16. Two linear models utilizing both the weight predicted biomass of 12 P. massoniana as independent and either the TLS derived biomass from the total volume (f(x) = 0.98x + 4.06, adj. r² = 0.88 ) or the TLS derived biomass from the volume of compartments with a diameter larger than 10 cm expanded with a biomass expansion factor (f(x) = 1.00x + 6.46, adj. r² = 0.85 ) as dependent variable.

Figure 16. Two linear models utilizing both the weight predicted biomass of 12 P. massoniana as independent and either the TLS derived biomass from the total volume (f(x) = 0.98x + 4.06, adj. r² = 0.88 ) or the TLS derived biomass from the volume of compartments with a diameter larger than 10 cm expanded with a biomass expansion factor (f(x) = 1.00x + 6.46, adj. r² = 0.85 ) as dependent variable.

5. Discussion

For biomass prediction utilizing the total TLS derived volume we consider only the P. massoniana results as acceptable. The average prediction error for this species is below 3%. Although under- and overestimation errors for the modeling of twigs could be detected visually (Figure 13), those errors seem to balance each other out. The modeling of the stems was considered to be accurate.

E. fordii estimates did underestimate the weight predicted biomass with an average error of ∼ −17%. This is partially caused by the deletion of small twigs during the denoising operations. Biomass and volume estimates do coincide quite well regarding the error (−17.24% and −17.47% respectively). The volume estimates for ${\mathit{volume}}_{10}^{\mathit{TLS}}$, ${\mathit{volume}}_7^{\mathit{TLS}}$ and ${\mathit{volume}}_4^{\mathit{TLS}}$ revealed all an error of ∼-8%. This leads to the conclusion, that only twigs smaller 4 cm diameter have been deleted. They account for ∼ 11% of the total biomass of an E. fordii. The tree models of this species were considered to be the best models by visual inspection. Good visual results and the underestimation of ∼-8% for these diameter classes made us assume that there was a systematical error occurrence during ground truth measurements. As the analysis was completed one year after the data collection, the source of this error is difficult to detect. We checked the accuracy of the two used scales and the quality of the volume prediction of the stem discs, but no systematical error could be found. The non detectable error and the good visualization results are contradictory and make a final interpretation of the underestimation impossible.

Q. petraea scans provided the worst results. The overestimation is ∼34% and the standard deviation of the error is the largest within all discussed tree species. We account the error to mainly the rainy conditions during the scanning campaign. Also wind in combination with the height of the trees made the scan quality in the crown low. In Figure 4(c) a trained human eye cannot see the real branching structure, which made it impossible for us to find an automatism to detect such a structure. Our method found largely overestimated and erroneous cylinders for the twigs, which lead to the huge overestimation of ∼34% (Figure 4(d)). Underestimated stem cylinders (Figure 3(d)) in comparison did not affect the total tree volume in a significant manner.

For all three species analyses of stem segments (Figure 6, Figure 10 and Figure 15) revealed that above a height of 10 m the segments’ volumes tend to be underestimated. The angle between the stem surface and the laser beam gets smaller with increasing height. This results in both lower point density in non occluded areas as well as larger occlusion areas (Figure 3(c)). A more detailed analysis of our results revealed that the underestimating error does not occur for the majority of the segments. Eight out of 30 Q. petraea segments above a height of 10 m were hugely underestimated (Figure 3(d)), all other 22 segments revealed a prediction quality comparable to their lower counterparts (Figure 3(b)).

The second approach to estimate biomass^TLS∗ yields much better results of biomass prediction (∼3%, ∼-8% and ∼3% respectively), because the prediction errors for ${\mathit{volume}}_{10}^{\mathit{TLS}}$ are much lower than the errors for ${\mathit{volume}}_0^{\mathit{TLS}}$. This could be partially based on the fact that the TLS modelled biomass of little twigs is based on density values of stem samples. However, the BEFs used in this approach to extrapolate the total above ground volume of the trees from the ${\mathit{volume}}_{10}^{\mathit{TLS}}$ was determined during the destructive sampling of the trees. In order to keep our method non-destructive, an approproate BEF would have to be derived from other sources adding further uncertainty to the estimator.

Conformity between biomass and volume estimator errors lead to the conclusion that the density values [33] are of high accuracy if representing stem density.

6. Conclusions

Our proposed method for the estimation of the above ground leafless biomass of trees yielded good results for the tree compartments with a diameter above 10 cm. The determination of stem volume and tree density seem to be reliable. Large errors occured in the volume estimation of the thinner tree compartments located in the tree crowns. Lower quality of the point clouds in the crown, mainly due to occlusion, wind and precipitation, lead to a huge overestimation (in leaf-less Q. petraea) or underestimation (in evergreen E. fordii) of twig volume.

Laserscanners like the FARO Laser Scanner Focus^3D can be mounted on a tower and be operated by remote control, which makes it possible to position the scanner closer to the tree crown, improving the point cloud quality in this area. Future research will show if point clouds improved in this or in other ways can yield good biomass predictions for all above ground compartments of trees with errors similar to P. massoniana results. Furthermore our results underline that it is necessary to conduct the scanning operations under good weather conditions to achieve maximum data quality. The fact that the scanning operations are costly and time-consuming, however, makes long-term planning of a scan campaign necessary, so good weather conditions cannot always be met.

In cases where our method yields good results it represents an excellent alternative to the destructive sampling of trees and can for example be used to calibrate allometric models in a non-destructive way.

7. Future Work

All used point clouds are available on the project homepage [40]. We want to encourage other analysts in our field to use this data and share the results to gain further knowledge about weaknesses and strength of different methods.