- Email: [email protected]
Deriving airborne laser scanning based computational canopy volume for forest biomass and allometry studies
ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Contents lists available at ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
Deriving airborne laser scanning based computational canopy volume for forest biomass and allometry studies Jari Vauhkonen a,⇑, Erik Næsset b, Terje Gobakken b a b
University of Helsinki, Department of Forest Sciences, Latokartanonkaari 7, P.O. Box 27, FI-00014 University of Helsinki, Finland Norwegian University of Life Sciences, Department of Ecology and Natural Resource Management, Høgskoleveien 12, P.O. Box 5003, NO-1432 Ås, Norway
a r t i c l e
i n f o
Article history: Received 10 January 2014 Received in revised form 19 June 2014 Accepted 2 July 2014
Keywords: Light Detection and Ranging (LiDAR) Forest inventory Tree allometry Delaunay triangulation Alpha shape Simplicial homomorphism Persistent homology
a b s t r a c t A computational canopy volume (CCV) based on airborne laser scanning (ALS) data is proposed to improve predictions of forest biomass and other related attributes like stem volume and basal area. An approach to derive the CCV based on computational geometry, topological connectivity and numerical optimization was tested with sparse-density, plot-level ALS data acquired from 40 field sample plots of 500–1000 m2 located in a boreal forest in Norway. The CCV had a high correspondence with the biomass attributes considered when derived from optimized filtrations, i.e. ordered sets of simplices belonging to the triangulations based on the point data. Coefficients of determination (R2) between the CCV and total above-ground biomass, canopy biomass, stem volume, and basal area were 0.88–0.89, 0.89, 0.83– 0.97, and 0.88–0.92, respectively, depending on the applied filtration. The magnitude of the required filtration was found to increase according to an increasing basal area, which indicated a possibility to predict this magnitude by means of ALS-based height and density metrics. A simple prediction model provided CCVs which had R2 of 0.77–0.90 with the aforementioned forest attributes. The derived CCVs always produced complementary information and were mainly able to improve the predictions of forest biomass relative to models based on the height and density metrics, yet only by 0–1.9 percentage points in terms of relative root mean squared error. Possibilities to improve the CCVs by a further analysis of topological persistence are discussed. Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
1. Introduction Carbon accounting systems and policies promoting the use of renewable energy have drawn an increasing attention toward extensive yet detailed mapping of forest biomass. Wall-to-wall estimation and mapping of forest biomass calls for various airborne and satellite remote sensing (RS) based solutions operated at local to national scales. The most recent overviews of forest biomass assessments based on RS data are provided by Koch (2010), Zolkos et al. (2013) and Popescu and Hauglin (2014). RS of forest biomass is based on indirect relationships between remotely sensed and field measured attributes. An increasingly popular approach to provide auxiliary RS data is airborne laser scanning (ALS) that produces three-dimensional (3D) vegetation height profiles from which various forest biophysical properties such as mean height (Nilsson, 1996; Næsset, 1997a; Magnussen and Boudewyn, 1998; Magnussen et al., 1999), basal area (Means ⇑ Corresponding author. Present address: University of Eastern Finland, Yliopistokatu 7, P.O. Box 111, FI-80101 Joensuu, Finland. E-mail address: jari.vauhkonen@uef.fi (J. Vauhkonen).
et al., 2000), stand volume (Næsset, 1997b; Means et al., 2000), and stem size distribution (Gobakken and Næsset, 2004) can be derived. Since the tree stems constitute the single most important forest above-ground biomass (AGB) component, these findings even apply to AGB predictions. However, ALS also allows assessing individual biomass components (Popescu, 2007; Næsset and Gobakken, 2008; Kotamaa et al., 2010; Hauglin et al., 2012), and may in fact produce more detailed information with respect to canopy biomass (i.e. branches and foliage) than conventional methods based on field measurement and modeling, as verified by destructive sampling (Hauglin et al., 2013; Kankare et al., 2013). ALS-assisted large-area assessments most often rely on areabased analyses (e.g. Næsset, 2002) of the distributions of height values extracted from sparse-density data (<1 pulses per m2; e.g. Næsset, 2007; Maltamo et al., 2011; Woods et al., 2011). Notably, such analyses employ only one of the three dimensions available in the data. Information on the horizontal distribution of forest patches may be provided by analyzing the neighborhoods of the computation units such as grid cells or segments, or their internal variation as depicted by features derived from interpolated canopy height models (Hyyppä et al., 2012; Pippuri et al., 2012; Packalén
http://dx.doi.org/10.1016/j.isprsjprs.2014.07.001 0924-2716/Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
58
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
et al., 2013). Yet, the 3D information in the point data is usually ignored. In light of the theories on tree allometry and stem development (Shinozaki et al., 1964a,b), the method development should ideally aim at quantifying foliage mass due to its strong allometric links with other tree attributes. Optimizing the use of the existing predictor space in various machine learning frameworks has however been a more popular approach than expanding it with such an allometric reasoning. Earlier studies carried out at the individual tree level (Vauhkonen, 2010a) propose triangulations of the 3D point data to preserve the three-dimensional properties and produce foliage mass equivalent information for various analyses. A particularly interesting approach is the concept of 3D alpha shapes (Edelsbrunner and Mücke, 1994), or a-shapes, in which a predefined parameter a is used as a size-criterion to determine the level of detail in the obtained triangulation. Volume and complexity metrics based on a-shapes were used as proxy variables for predicting the species and stem attributes of individual trees by Vauhkonen et al. (2008, 2009, 2010a,b), Reitberger et al. (2009), Rentsch et al. (2011), and Yao et al. (2012), while Vauhkonen et al. (2012) provided further insight on using these metrics to improve area-based estimations of species-specific plot volumes. An important step in the analyses based on a-shapes is the selection of a, for which the selected value should enable separation of void spaces from those populated by canopy biomass. Vauhkonen et al. (2008, 2009, 2010a,b, 2012) evaluated a range of a0 s, but rather than selecting one particular value, they used metrics quantifying the difference of the shape obtained with the range of a0 s to the convex hull of the point data, which corresponds to a-shapes with a ? 1. The other studies mentioned in the previous paragraph do not explain the selection of a. Korhonen et al. (2013), attempting to model individual tree crown volume, selected a quasi-optimal a so that the resulting a-shape enclosed the point data within a single connected component. However, they found the volume based on the convex hull to have a closer correspondence with the field-measured reference crown volume, which also included void space due to the applied measurement principle. Most interestingly, Vauhkonen (2010b) proposed iterating over a sequence of a’s and by that vertically delineating individual tree crowns, which was further tested in an area-based application by Maltamo et al. (2010). This approach has a formal basis in topological connectivity and is particularly related to filtration of simplicial complexes (e.g. Delfinado and Edelsbrunner, 1995). Triangulating ALS point data corresponds to subdividing the underlying space of the points into simplices, which results in weighted simplicial complexes (e.g. Edelsbrunner, 2011), with weights quantifying the (empty) space delimited by the points. We hypothesize that filtering these complexes according to the weights could provide a means to reconstruct the volume populated by biomass and separate that volume from the canopy voids (Fig. 1). Thus, the purpose of this study was to develop a filtering and optimization procedure for deriving a ‘‘biomass-optimal’’ computational canopy volume (CCV) from the triangulations of plot-level ALS point clouds. Two criteria were employed to filter the triangulations of the point data, and the resulting filtration parameters and CCVs were evaluated on field data with respect to the total and canopy biomass, stem volume and basal area as target biophysical properties. 2. Methods 2.1. Study area and data The study area is located in the municipality of Aurskog-Høland, in southeast Norway (59°500 N, 11°340 E, 120–390 m a.s.l.). The dominant tree species are Norway spruce (Picea abies (L.) Karst.)
Fig. 1. Derivation of the CCV from the ALS point data of an example plot (schematic). Each step constitutes the input of the next consecutive step.
and Scots pine (Pinus sylvestris L.). Younger stands have a larger portion of deciduous species than mature stands. Birch (Betula pubescens Ehrh.) is the dominant deciduous species. For further details, see Magnussen et al. (2010), who studied the same area.
2.1.1. Field plots Altogether 40 circular field plots (36 of 1000 m2 and 4 of 500 m2) were measured during the fall of 2007 and winter of 2008. These plots were located along five systematically spaced lines, but the exact plot locations within the lines were determined subjectively to obtain a sample of 20 spruce dominated and 20 pine dominated stands. Furthermore, 5 and 15 plots covered young and old growth for each of the two species, respectively. Each plot center was geographically referenced with a survey grade Global Positioning System and Global Navigation Satellite System receiver (Topcon LegacyE) applying differential post-processing by the Pinnacle software package version 1.00 (Anon, 1999). Based on the positional standard errors reported by Pinnacle, the estimated accuracy of the planimetric plot coordinates ranged from 0.1 to 0.35 m, with an average of 0.12 m. All trees with a diameter at breast height (DBH) >5 cm were measured for the DBH and species. For the dominant tree species, the height of two trees in each of five diameter classes – covering the entire range of diameters – was measured using a Vertex hypsometer. Up to five trees of each minor species, selected with a probability proportional to stem basal area, were measured for
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
height. The total number of height measurements per plot ranged from 11 to 20 with an average of 17. Six nonlinear height–diameter regression models were developed and the missing tree heights were predicted with these models. The models were species (Norway spruce, Scots pine, and deciduous spp.) site productivity (high, low) specific. The plot-level forest attributes considered in this study were AGB and foliage biomass (FB), the latter comprising branches and foliage, basal area (G), and stem volume (V). Single-tree biomass fractions were estimated for all trees using Marklund’s (1988) biomass models based on the DBH, height, and tree species. The biomass models for birch were used for all deciduous trees. The tree-level biomass fractions were summed to AGB and FB. Plotlevel G was computed based on summing from the diameter measurements. Individual stem volumes were estimated using speciesspecific models based on the DBH and height (Braastad, 1966; Brantseg, 1967; Vestjordet, 1968), and plot-level V was computed as the sum of the individual stem volumes. See Magnussen et al. (2010) for further details regarding the measurements and computations. 2.1.2. ALS data The ALS data were acquired under leaf-on conditions 8–10 June 2005 with a supplementary flight conducted on 6 September 2005 to fill in a minor gap in the data. The acquisition was carried out using a fixed-wing Piper PA31-310 aircraft flying at an altitude of approximately 1850 m above ground level at a speed of 75 ms1. The Optech ALTM 3100 scanner was operated with a pulse repetition frequency of 50 kHz. The scan frequency was 71 Hz, resulting in a point density of approximately 0.7 m2 on the ground. The maximum scan angle was 15° but pulses emitted at an angle >13° were discarded during the subsequent data processing. Standard pre-processing including strip adjustment and ground classification was performed by the contractor (see Maltamo et al. (2010) for details). As a result of the pre-processing, the echo heights were normalized with respect to the ground level. The main results of this study were based on all echoes reflected from the areas of the field plots above a ground threshold of 0.5 m. In addition, only the first echoes (i.e. ‘only’ and ‘first-of-many’ of up to four echoes recorded per pulse) extracted employing a ground threshold of 2 m were included in a benchmark calculation. 2.2. An overview of the analysis The areas of interest (AOIs) of our approach were the individual plots (or grid cells in the case of prediction) corresponding to typical area-based analyses of ALS data (Næsset, 2002). We defined the CCV of an individual AOI as the sum of the volumes of individual simplices (tetrahedra) belonging to a triangulated ALS point cloud of that particular AOI (Fig. 1). The triangulation and the number of potential filtrations (Steps 2–3 of Fig. 1) were determined solely based on the topological properties of the point data (see Section 2.3). The fundamental problem was thus to determine a proper filtration (Step 4 of Fig. 1) such that the simplices considered in the computation of the CCV represented the areas populated by the forest biomass. As hypothesized in Section 1, we assumed that an AOI-specific filtration was required to focus the computation on the simplices representing forest biomass and to exclude those representing empty space. According to the reasoning of the previous paragraph and in order to apply the proposed approach, one should be able to filter the triangulations of all the AOIs considered such that the CCVs produced were related to the field-observed biomasses as closely as possible. To define such filtrations, we applied the practical two-stage procedure proposed by Næsset (2002), in which (i) models to predict the forest attributes of interest for the individual AOIs
59
are fit based on a set of training field plots; and (ii) the resulting models are applied to all the AOIs of the entire inventory area to produce wall-to-wall predictions. In our case, the corresponding stages were: (i) to optimize the filtrations with respect to the selected forest biophysical attributes observed from the training plots and to model the parameters producing these optimal filtrations; and (ii) to predict corresponding filtrations to all the AOIs based on information extracted from the point data. Lacking independent validation data, we followed the two-stage approach with cross-validation, i.e. predicted step (ii) to the training field plots and evaluated this prediction as if these plots represented totally independent data with unknown field observations. 1. Define filtering rules to be able to adjust the level of detail in the triangulations (Section 2.3). 2. Optimize the degree of filtering such that the CCVs extracted from the filtrations were in a linear relationship with the reference biomass values observed on the corresponding training field plots (Section 2.4). 3. Analyze the filtration parameters obtained in Step 2 with respect to forest biomass and structural attributes (Section 2.5). The main aim of this analysis was to assess the conditions under which the filtration parameters could be estimated in a wall-towall ALS inventory based on a set of training field plots. 4. Predict the filtration parameters and extract the CCV for the AOIs of the entire inventory area (Section 2.6). In the absence of a prediction grid or independent validation data, we used the set of training plots to mimic the AOIs of the inventory area, and evaluated the predictions obtained in the same data used in the earlier steps following the principle of cross-validation (Section 2.7). 2.3. Deriving triangulations and their filtrations The terminology and notation in the following are adapted from Delfinado and Edelsbrunner (1995), to which the reader is referred for additional details. A triangulation T of a finite set of points S is a subdivision of the convex hull of S into simplices according to the rules of the particular triangulation (see, e.g., Preparata and Shamos, 1985). T is a simplicial complex that is subsettable to a finite number of subcomplexes Ci, i = 1, 2, . . ., s, where i indexes the position and s is the number of different subcomplexes in a list of potential filter values F defining which simplices belong to the subcomplex (see the next paragraph). When F is ordered (here ascendingly), the sequence Ci constitutes a filtration, in which each complex is a sub-complex of its successor, and VCi is the volume of the underlying space of subcomplex Ci according to this filtration. The point data of each individual plot were triangulated separately and filtered by two alternative criteria. First, the triangulations were filtered with respect to the tetrahedral volume (v), i.e. F was an ascending sequence of the individual v’s and Ci included the tetrahedra with v 6 F[i]. Second, the filtration corresponded to the family of 3D a-shapes (Edelsbrunner and Mücke, 1994), each of which consist of those simplices of T, which have an empty circumscribing sphere with a squared radius smaller than a. Although an a-shape is defined for every non-negative real a, there are only finitely many different a-complexes (see Delfinado and Edelsbrunner, 1995, p. 782) and F was thus an ascending sequence of those a-values which altered the a-shape. All computations were based on Delaunay triangulations and implemented with C++ and an open-source Computational Geometry Algorithms Library (Pion and Teillaud, 2013; Da and Yvinec, 2013). Fig. 2 gives an example and illustrates the differences between the two types of filtrations. The filtered v-complexes constitute of individual tetrahedra simply sorted according to their volume, whereas the a-complexes are constructed via more delicate
60
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
analysis on simplices exposed to the a parameter. The CCV extracted from both the filtration types seem to be related to the canopy volume originating from the actual locations of the trees, but also on the degree of the filtering applied (Fig. 2).
AGB-v, FB-v, AGB-a, and FB-a are used to denote these op – filter pairs.
2.4. Optimizing the filtrations with respect to forest biomass
The analyses of Section 2.3 provided a range of filtration parameters and thus CCVs that were in a quasi-optimal linear relationship with forest biomass in at least one out of the 100 solutions. The correspondence of the CCVs and the forest attributes of interest was further quantified by fitting an allometric power equation (White and Gould, 1965):
To optimize the filtration parameters, all the studied plots were considered jointly. For convenience, let V C i ;p in the following text denote the CCV based on a filtration of the ALS point data of plot p. The optimization was based on minimizing the residual errors of a linear regression model:
op ¼ b1 þ b2 V C i ;p þ ep ;
2.5. Statistical analyses
k
yp ¼ bxp þ ep ;
ð2Þ
ð1Þ
where op was the reference biomass, b1 and b2 the model parameters, and ep the residual error of plot p. Both AGB and FB were tested as op, since the CCV based on ALS point data was expected to represent canopy biomass (branches and foliage) rather than that originating from the stems. The initial values of b1 and b1 were obtained by fitting the model to V C i ;p with i selected randomly for each plot. In the optimization, V C i ;q of plot q which produced the largest residual error was altered to V C i 10;q in successive rounds of 500 iterations, until either the error, determined as 1 – the obtained coefficient of determination (R2), decreased below 0.0009 or did not change during the last round of iterations. The whole optimization procedure was repeated 100 times with op as either AGB or FB, and Ci filtered by v or a. In the following text,
where yp was the considered forest attribute (either AGB, FB, G, or V), xp the optimized CCV, b and k the model parameters, and ep the residual error of plot p. The model fitting was carried out with the nls function of R (R Core Team, 2012). The ranges of the filtration parameters belonging to the optimal solutions were analyzed. The analysis was carried out with respect to both the absolute values and percentiles of the distributions of these values. The filtration parameters were examined with respect to forest structural attributes, of which the plot basal area had the highest correlation with them. The sensitivity of the CCV to the filtration parameters was analyzed by randomly selecting 10,000 realizations of the filtration parameters from the optimal range and recording the obtained R2 between the CCV and the predicted biomass attributes.
Fig. 2. Spatial distribution of a filtered 3D triangulation of an example plot. Circles depict the relative sizes and locations of the centroids of the tetrahedra belonging either to the whole triangulation (subplot A; V C i = 11251.6 m3) or filtrations a < 21.0 m3 (B; V Ci = 5388.6 m3), v < 3.7 m3 (C; V C i = 1538.9 m3), or a < 11.3 m3 (D; V C i = 2620.8 m3). The filter values in B–D correspond to the maximum a, median v, and median a, respectively, of this plot. The red crosses depict the relative sizes and positions of the field-measured trees and the outer circle indicates the plot border. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
2.6. Predicting the filtration parameters and biomass attributes The area-based two-stage procedure proposed by Næsset (2002) was tested to evaluate the feasibility of predicting both the filtration parameters and biomass attributes in a practical ALS-based forest inventory. The filtration parameters were correlated against the most common ALS-based predictor variables (cf. Næsset, 2002), i.e. the mean and standard deviation of the height values, the proportion of echoes above 2 m, and the 5th, 10th, 20th, . . ., 90th, and 95th percentiles and the corresponding proportional densities of the ALS-based canopy height distribution, calculated according to Korhonen et al. (2008, pp. 502–503). The CCV was then extracted according to the filtration parameters modeled by the ALS variables. Further, combinations of CCV with other ALS variables were inserted in Eq. (2) to find out the added value of the CCV over other potential predictors. The x variables tested in Eq. (2) were (1) CCV; (2) a height (H) and a density (D) variable derived from the canopy height distribution (see above); and (3) H + D + CCV. The plot-level mean height and the proportion of echoes above 2 m, both based on the first echoes, were found yield the highest correlations with the predicted biomass attributes and were thus used in all related experiments as H and D, respectively. 2.7. Evaluation and performance measures The correspondence of the estimates with the reference data was evaluated by graphical assessments and in terms of R2, root mean squared error (RMSE) and mean difference (bias), the latter two calculated according to:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 ^ i¼1 ðyi yi Þ RMSE ¼ n Pn bias ¼
^ yi Þ n
i¼1 ðyi
ð3Þ
ð4Þ
^i and y are the prewhere n is the number of observations, and y dicted and reference values, respectively. The relative RMSE and bias were calculated by dividing the absolute RMSE and bias values by the mean value of the reference attribute. A paired t-test was used to test the significance of the mean difference between the observation vectors. In the following text, the term ‘significant’ refers to the statistical significance of this difference at the 95% confidence level. 3. Results The following results are based on using all echoes and a ground threshold of 0.5 m. A brief comparison with the first return data and a threshold of 2 m is given as the last paragraph of this section. The optimization (Section 2.3) required 6100 or 5100 iterations, on average, when the filtrations were based on v or a, respectively. The average R2 of Eq. (1) was 0.99 or 0.92 when op was either AGB or FB, respectively. Altogether 88, 85, 92, and 90 of the 100 iterations based on AGB-v, FB-v, AGB-a, and FB-a, respectively, had R2 P 0.985. The following results are based only on those iterations. The ranges of the filtration parameters included in the optimal solutions are shown in Fig. 3. The ranges varied between the plots and their mean values increased significantly according to an increasing plot basal area (R2 = 0.42–0.61 depending on the filtration). Further, the length of the range was significantly correlated with plot basal area (R2 = 0.12–0.30) in the filtrations based on v. The relationship with the basal area was more pronounced when examined in terms of percentiles (Fig. 3). The mean values (R2 = 0.74–0.88) and ranges (R2 = 0.20–0.44) of all the filtrations
61
were significantly correlated with the basal areas. The betweenplot variation in the data density (0.94–3.66 echoes/m2) correlated insignificantly with the absolute values, and the correlation with the percentile values was only slightly significant (R2 = 0.09– 0.13). The mean values of the filtration parameters were not significantly different whether the optimization was based on either AGB or FB (see also Fig. 3). The CCVs extracted with the mean values of the filtration parameters were in 1:1 relationship (R2 = 1.0 based on all the filtrations) with those biomass attributes that were considered in the optimization. The relationships between these CCVs and the other forest attributes considered are shown in Fig. 4. The CCV derived from the optimized filtrations had R2 of 0.88–0.89, 0.89, 0.83– 0.97, and 0.88–0.92 with AGB, FB, V, and G, respectively, depending on the applied filtration. The CCVs derived from the filtrations based on v were generally larger than those based on a, but there were only minor differences between the relationships with the CCVs and the forest attributes (Fig. 4). The sensitivity analysis of deriving the CCVs with filtration parameters sampled from various ranges around the means indicated that these had to be close to the quasi-optimal values to obtain high accuracies. A worst-case R2 > 0.76 could be obtained for all the considered forest attributes, when the filtration parameters were at most within one standard deviation around the mean (Fig. 5). The correspondence was lower when sampling around two standard deviations (R2 = 0.67–0.96) or the whole range (R2 = 0.60– 0.96). However, the average result obtained from the 10,000 realizations always yielded an R2 > 0.73. Of the considered forest attributes, G had the lowest R2 values. The filtrations based on a generally provided higher R2 values (Fig. 5). The filtration parameters were found to be significantly correlated with several ALS-based height and density metrics. As could be expected based on Fig. 3, the correlations between these metrics were higher with the percentiles of the filtration parameter distributions than the absolute values. Several individual height and density variables had R2 > 0.6 with the percentiles of the filtration parameter distributions (detailed results not shown). The a-values were generally found to have higher correlations and better predictability than v, and thus the following analyses are presented for a-based filtrations only. A model 33.9 + 27.2 ln (H10) + 47.4 D2, where ln(H10) is the natural logarithm of the 10th height percentile and D as defined earlier (Section 2.4), had R2 = 0.83 for predicting the percentile of the a-based filtration parameter. With such a model, altogether 20 or 32 of the 40 predictions laid within one or two standard deviations of the reference values (Fig. 6), respectively. The CCV derived with the predicted a had R2 of 0.90, 0.77, 0.90, and 0.84 with AGB, FB, V, and G, respectively (Table 1). The obtained CCVs always had a higher R2 than any individual height and density metric, but H + D always produced a better fit than CCV with respect to both R2 and RMSE (Table 1). However, the performance of H + D could be improved by including the CCV as an additional variable in the model in all other cases except FB. The improvement of the other predictions was in order of 0.6–1.9 percentage points, least for AGB and most for G. The residuals of all models were distributed normally and did not show trends toward forest structural attributes or the between-plot variation in the pulse density. None of the mean differences were statistically significant. When repeating the experiment in the data set including only the first echoes extracted above a higher ground threshold, the only significant differences were observed with respect to the running times of the optimization algorithm and the accuracy of the results related to FB. The optimization required less iterations (2800–3500), which were also faster due to the lower number of simplices considered. The proportion of FB was presented more
62
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Fig. 3. The ranges (lines) and means (dots) of the absolute values (left) and percentiles (right) of the optimized filtration parameters of AGB-v (top), FB-v, AGB-a, and FB-a (bottom). The x-axes are ordered according to the plot basal area, which increases from left (13.6 m2/ha) to right (42.0 m2/ha). The maximum v beyond the scale of the two topmost graphs was 102.3 m3 produced by a single outlier.
accurately in all performance measures – most remarkably, the RMSE obtained by the H + D + CCV model was approximately one percentage point lower than that reported in Table 1. Otherwise, the results based on the two point data sets showed decimal-level differences in terms of the performance measures presented above.
4. Discussion The results of this study indicate that the developed new feature type, computational canopy volume, has a high potential to improve ALS-based estimations of forest biomass attributes. This is a logical outcome, since the canopy volume is directly related to the amount of foliage and canopy bulk density (cf. Riaño et al., 2003, 2004). The approach to derive the CCV adopted here was based on triangulating the 3D point data and filtering the triangulations for empty space. When the degree of the filtration was determined in an optimization, the derived CCVs were more closely related to the biomass attributes that were considered. When the filtration parameters were predicted by means of a simple
linear regression model, the accuracies were lower and approximately at the same level than those obtained with simple height and density metrics. However, the accuracy levels are overall higher than those typically reported in similar boreal forest conditions, likely owing to the larger plot size used here, which may partially explain the small improvement obtained by our approach. Yet the determination of the filtration parameters may require further considerations when applied in practical inventories (see further text), the present results are considered as an important step in the development of ALS-based auxiliary variables for forest biomass inventories and should warrant further testing of the presented approach. The approach to describe individually delineated tree crowns with volumetric properties extracted from the ALS data was first tested by Vauhkonen et al. (2008), who proposed a density of 3 pulses m2 as a minimum data requirement for such an analysis. However, the results of this study together with those presented by Vauhkonen et al. (2012) show that the operational densities of <1 m2 are feasible when triangulations of areas larger than individual tree segments are considered. Earlier, determining forest
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
63
Fig. 4. The relationships between AGB and FB (top), V (middle), and G (bottom), and CCV optimized with AGB (left) or FB (right). CCVs based on v and a are depicted by circles and crosses, respectively. The coefficients shown are based on fitting Eq. (2) to the data.
canopy volume by simply multiplying the laser-based canopy height and area, potentially applying a correction for mean canopy cover, has been proposed by Nelson et al. (1984), Næsset (1997b), Lefsky et al. (1999), Riaño et al. (2003, 2004), Chen et al. (2007), and Hollaus et al. (2009). Further, an iterative surface wrapping technique based on delineated tree crowns (Kato et al., 2009) or stacking horizontal slices of the height values (Tang et al., 2013), both operating with the point data but in 2.5D space, have been proposed. However, deriving the canopy volume from the 3D triangulations of the point data differs from the previous approaches and has certain benefits over them. A main benefit is the ability to adjust the level of detail in the triangulated point cloud by means of the filtrations, thus accounting for canopy gaps and detailed properties existing in the data.
The filtration parameters were determined in an optimization, in which the aim was to remove the empty space within the point cloud and focus on the remaining canopy structures populated by biomass. The applied optimization forced a linear relationship between the CCV and biomass. Whether the true form of the allometric relationship between these attributes was known, a corresponding model could have been fit in the optimization. However, since the correspondence of the CCV with the actual canopy volume would in any case depend on the definition and measurement of these attributes (see Korhonen et al., 2013), the developed CCV is seen as a computational feature complementing the ALS-based height and density information rather than an estimate of canopy volume with a strictly allometric definition. For example, our experiment with the first echoes and a higher ground
64
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Fig. 5. Sensitivity of the filtration parameters (left: v, right: a) to R2 between the extracted CCV and the forest attributes considered. In the analysis, 10,000 realizations of the filtrations parameters were sampled from within the mean of the optimal solutions ± one (1sd) or two (2sd) standard deviations or the full range.
Fig. 6. Predicted vs. reference a (black dots). One and two standard deviations of the reference values are depicted by the solid and broken lines, respectively.
threshold showed that the CCV can only be determined for a certain part of the canopy and still be linked accurately with the biomass attributes. Two types of filtrations and optimization schemes were tested for deriving the CCV. It should be noted that the CCVs and the filtrations overall could likely be related to several other forest biophysical properties than the biomass-based variables considered. For example, the spatial distribution of the tetrahedra could potentially be related to the spatial pattern of the trees (see Fig. 2), the prediction of which has been found difficult based on both the area- and tree-level analyses of ALS data (cf. Packalén et al., 2013). The results promote the usability of the a-shapes (Edelsbrunner and Mücke, 1994), which constitute an interesting geometric modeling tool also based on earlier related studies (e.g. Vauhkonen et al., 2012). The approach has a profound mathematical formalism on topological connectivity rather than the alternative approach based on the tetrahedral volume, which was
considered as a straightforward and potentially more easily implementable approach. Except for the optimization, the implementations were based on efficient data structures (Pion and Teillaud, 2013; Da and Yvinec, 2013) and are computationally feasible for practical applications. Compared to FB, the use of AGB in the optimization provided better results when predicting the other forest attributes (G, V). The use of the first echo data improved the predictions related to FB, yet AGB was preferred also in those situations. Although it was initially assumed that the ALS data would rather depict the canopy biomass, this observation may be related to deriving the reference values by models, which are typically found more inaccurate with respect to the individual biomass components than the total AGB. Further, it is acknowledged that species-specific differences should be accounted for to improve the estimations, since the amount and allocation of foliage and branches differ between the species considered (see also Næsset and Gobakken, 2008). Since the largest biomasses in our data originated from spruce stands and smallest from pine stands, a plot-level analysis of species would not have been independent of the tree size. Besides species, a future topic to be studied is the plot size, which was purposely large to allow a comprehensive analysis of the studied phenomenon. It is reasonable to assume that different species, plot sizes and ALS sensors used would have a similar effect on the developed features as earlier described in the literature for other predictor variables (e.g. Næsset and Gobakken, 2008; Gobakken and Næsset, 2008). A proper selection of the filtration parameters turned out to be highly sensitive toward the obtained results (cf. Figs. 5 and 6). A simple approach based on relating the a-parameter to the height and density metrics of ALS data was tested as a means to determine the degree of filtration for a practical prediction of forest biomass. The approach was found to work, yet producing a suboptimal result. Improvements could potentially be obtained by further examining the inherent algebraic properties of the triangulations such as simplicial homomorphism and topological persistence (Zomorodian, 2005). For example, Martynov (2008), computing persistent homology groups at different spatial resolutions and assuming that important features and structures persisted over a wide range of scales, de-noised geometric models of power line towers reconstructed from the ALS point data. Yet such an analysis is acknowledged to be considerably simpler with manmade objects than vegetation, which has natural structural variations, computing persistence homology from vegetation point clouds is identified as an interesting topic for future studies. Besides ALS-based large-area inventories, the proposed approach could have a high practical importance toward applications requiring geometrically explicit models of forest structure, which are currently constructed by field-intensive terrestrial laser scanning (TLS). Such applications include mapping and modeling of vegetation clumping and effective leaf area index (LAI, Hopkinson et al., 2013), vertical distribution of photosynthetically active radiation (PAR) and the fraction of PAR (fPAR, van Leeuwen et al., 2013), light transmittance and the spatial distribution of foliage (Béland et al., 2014), and forest primary production (Seidel et al. 2013) and growth (Metz et al., 2013). Particularly the approach proposed here could provide a considerably higher spatial coverage when applied with the ALS data, while the input on acquiring the TLS data could be focused on collecting the optimization samples for this approach.
5. Conclusion A novel ALS-based feature type, computational canopy volume (CCV), showed potential to improve forest biomass estimations
65
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66 Table 1 Accuracies of predicting the forest attributes considered with combinations of CCV, H, and D. Attribute
Predictors
R2
RMSE
RMSE%
Bias
Bias%
AGB (t/ha)
CCV H+D H + D + CCV
0.90 0.91 0.92
18.8 17.0 16.4
15.9 14.4 13.8
0.3 0.1 0.0
0.3 0.1 0.0
FB (t/ha)
CCV H+D H + D + CCV
0.77 0.85 0.86
6.8 5.4 5.4
23.4 18.6 18.6
0.05 0.00 0.00
0.2 0.0 0.0
V (m3/ha)
CCV H+D H + D + CCV
0.90 0.91 0.92
32.7 31.1 29.4
15.8 15.1 14.2
0.4 0.5 0.03
0.2 0.2 0.0
G (m2/ha)
CCV H+D H + D + CCV
0.84 0.90 0.93
3.5 2.7 2.2
13.9 10.8 8.9
0.01 0.03 0.01
0.1 0.1 0.0
from sparse-density, area-based ALS data. The best-case R2 obtained between the CCV and total biomass, canopy biomass, stem volume, and basal area were 0.88–0.89, 0.89, 0.83–0.97, and 0.88–0.92, respectively. The R2 depended on the applied filtration and was in all cases based on optimization with forest biomass attributes. The magnitude of the required filtration increased according to an increasing basal area, and a simple prediction model with related ALS height and density metrics provided CCVs which had R2 of 0.77–0.90 with the aforementioned forest attributes. The derived CCVs mainly improved the predictions of forest biomass, yet only by 0–1.9 percentage points in terms of RMSE%, compared to models based on the conventional height and density metrics. A benchmark calculation showed that similar results could be obtained by restricting the analysis to specific echo types (here first echoes), thus potentially reducing the computational burden involved. The main bottleneck is the prediction of the filtration parameters, which could potentially be improved by further analyses of topological persistence of the triangulations and their filtrations.
Acknowledgements This study was enabled by the Visiting Researcher Grant (2012/ 2099) of the Norwegian University of Life Sciences and additionally supported by University of Helsinki Funds. The authors wish to thank Blom Geomatics (Norway) for providing and processing the ALS data. The authors also wish to thank K. M. Hauglin and V. Lien, both at the Norwegian University of Life Sciences, for the help in the fieldwork.
References Anon, 1999. Pinnacle User’s Manual. Javad Positioning Systems, San Jose, CA, p. 123. Béland, M., Baldocchi, D.D., Widlowski, J.-L., Fournier, R.A., Verstraete, M.M., 2014. On seeing the wood from the leaves and the role of voxel size in determining leaf area distribution of forests with terrestrial LiDAR. Agric. Forest Meteorol. 184, 82–97. Braastad, H., 1966. Volume tables for birch. Meddelelser fra Det norske Skogforsøksvesen 21, 265–365 (In Norwegian with an English summary). Brantseg, A., 1967. Volume functions and tables for scots pine. South Norway. Meddelelser fra Det norske Skogforsøksvesen 22, 689–739 (In Norwegian with an English summary). Chen, Q., Gong, P., Baldocchi, D., Tian, Y.Q., 2007. Estimating basal area and stem volume for individual trees from Lidar data. Photogramm. Eng. Remote Sens. 73, 1355–1365. Da, T.K.F., Yvinec, M., 2013. 3D alpha shapes. CGAL User and Reference Manual, CGAL Editorial Board, 4.3 ed.. Delfinado, C.J.A., Edelsbrunner, H., 1995. An incremental algorithm for Betti numbers of simplicial complexes on the 3-sphere. Comput. Aided Geom. Des. 12, 771–784.
Edelsbrunner, H., 2011. Alpha shapes – A Survey. Book Chapter Manuscript, Tessellations in the Sciences. Springer-Verlag,. Edelsbrunner, H., Mücke, E.P., 1994. Three dimensional alpha-shapes. ACM Trans. Graphics 13, 43–72. Gobakken, T., Næsset, E., 2004. Estimation of diameter and basal area distributions in coniferous forest by means of airborne laser scanner data. Scand. J. Forest Res. 19, 529–542. Gobakken, T., Næsset, E., 2008. Assessing effects of laser point density, ground sampling intensity, and field sample plot size on biophysical stand properties derived from airborne laser scanner data. Can. J. Forest Res. 38, 1095–1109. Hauglin, M., Dibdiakova, J., Gobakken, T., Næsset, E., 2013. Estimating single-tree branch biomass of Norway spruce by airborne laser scanning. ISPRS J. Photogramm. Remote Sens. 79, 147–156. Hauglin, M., Gobakken, T., Lien, V., Bollandsås, O.M., Næsset, E., 2012. Estimating potential logging residues in a boreal forest by airborne laser scanning. Biomass Bioenergy 36, 356–365. Hollaus, M., Wagner, W., Schadauer, K., Maier, B., Gabler, K., 2009. Growing stock estimation for alpine forests in Austria: a robust Lidar-based approach. Can. J. Forest Res. 39, 1387–1400. Hopkinson, C., Lovell, J., Chasmer, L., Jupp, D., Kljun, N., van Gorsel, E., 2013. Integrating terrestrial and airborne lidar to calibrate a 3D canopy model of effective leaf area index. Remote Sens. Environ. 136, 301–314. Hyyppä, J., Yu, X., Hyyppä, H., Vastaranta, M., Holopainen, M., Kukko, A., Kaartinen, H., Jaakkola, A., Vaaja, M., Koskinen, J., Alho, P., 2012. Advances in forest inventory using airborne laser scanning. Remote Sens. 4, 1190–1207. Kankare, V., Räty, M., Yu, X., Holopainen, M., Vastaranta, M., Kantola, T., Hyyppä, J., Hyyppä, H., Alho, P., Viitala, R., 2013. Single tree biomass modelling using airborne laser scanning. ISPRS J. Photogramm. Remote Sens. 85, 66–73. Kato, A., Moskal, L.M., Schiess, P., Swanson, M.E., Calhoun, D., Stuetzle, W., 2009. Capturing tree crown formation through implicit surface reconstruction using airborne Lidar data. Remote Sens. Environ. 113, 1148–1162. Koch, B., 2010. Status and future of laser scanning, synthetic aperture radar and hyperspectral remote sensing data for forest biomass assessment. ISPRS J. Photogramm. Remote Sens. 65, 581–590. Korhonen, L., Peuhkurinen, J., Malinen, J., Suvanto, A., Maltamo, M., Packalén, P., Kangas, J., 2008. The use of airborne laser scanning to estimate sawlog volumes. Forestry 81, 499–510. Korhonen, L., Vauhkonen, J., Virolainen, A., Hovi, A., Korpela, I., 2013. Estimation of tree crown volume from airborne lidar data using computational geometry. Int. J. Remote Sens. 34, 7236–7248. Kotamaa, E., Tokola, T., Maltamo, M., Packalén, P., Kurttila, M., Mäkinen, A., 2010. Integration of remote sensing-based bioenergy inventory data and optimal bucking for stand-level decision making. Eur. J. Forest Res. 129, 875–886. van Leeuwen, M., Coops, N.C., Hilker, T., Wulder, M.A., Newnham, G.J., Culvenor, D.S., 2013. Automated reconstruction of tree and canopy structure for modeling the internal canopy radiation regime. Remote Sens. Environ. 136, 286–300. Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A., Harding, D., 1999. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-Fir Western Hemlock forests. Remote Sens. Environ. 70, 339–361. Magnussen, S., Boudewyn, P., 1998. Derivations of stand heights from airborne laser scanner data with canopy-based quantile estimators. Can. J. Forest Res. 28, 1016–1031. Magnussen, S., Eggermont, P., LaRiccia, V.N., 1999. Recovering tree heights from airborne laser scanner data. Forest Sci. 45, 407–422. Magnussen, S., Næsset, E., Gobakken, T., 2010. Reliability of LiDAR derived predictors of forest inventory attributes: a case study with Norway spruce. Remote Sens. Environ. 114, 700–712. Maltamo, M., Bollandsås, O.M., Vauhkonen, J., Breidenbach, J., Gobakken, T., Næsset, E., 2010. Comparing different methods for prediction of mean crown height in Norway spruce stands using airborne laser scanner data. Forestry 83, 257–268. Maltamo, M., Packalén, P., Kallio, E., Kangas, J., Uuttera, J., Heikkilä, J., 2011. Airborne laser scanning based stand level management inventory in Finland. In:
66
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Proceedings of SilviLaser 2011 – 11th International Conference on LiDAR Applications for Assessing Forest Ecosystems, 16–20 October, 2011, Hobart, Australia.. Means, J.E., Acker, S.A., Fitt, B.J., Renslow, M., Emerson, L., Hendrix, C.J., 2000. Predicting forest stand characteristics with airborne scanning lidar. Photogramm. Eng. Remote Sens. 66, 1367–1371. Metz, J., Seidel, D., Schall, P., Scheffer, D., Schulze, E.-D., Ammer, C., 2013. Crown modeling by terrestrial laser scanning as an approach to assess the effect of aboveground intra- and interspecific competition on tree growth. Forest Ecol. Manage. 310, 275–288. Næsset, E., 1997a. Determination of mean tree height of forest stands using airborne laser scanner data. ISPRS J. Photogramm. Remote Sens. 52, 49–56. Næsset, E., 1997b. Estimating timber volume of forest stands using airborne laser scanner data. Remote Sens. Environ. 51, 246–253. Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sens. Environ. 80, 88–99. Næsset, E., 2007. Airborne laser scanning as a method in operational forest inventory: status of accuracy assessments accomplished in Scandinavia. Scand. J. Forest Res. 22, 433–442. Næsset, E., Gobakken, T., 2008. Estimation of above- and below-ground biomass across regions of the boreal forest zone using airborne laser. Remote Sens. Environ. 112, 3079–3090. Nelson, R., Krabill, W., MacLean, G., 1984. Determining forest canopy characteristics using airborne laser data. Remote Sens. Environ. 15, 201–212. Nilsson, M., 1996. Estimation of tree heights and stand volume using an airborne lidar system. Remote Sens. Environ. 56, 1–7. Packalén, P., Vauhkonen, J., Kallio, E., Peuhkurinen, J., Pitkänen, J., Pippuri, I., Strunk, J., Maltamo, M., 2013. Predicting the spatial pattern of trees with airborne laser scanning. Int. J. Remote Sens. 34, 5154–5165. Pion, S., Teillaud, M., 2013. 3D triangulations. CGAL User and Reference Manual, CGAL Editorial Board, 4.3 ed. . Pippuri, I., Kallio, E., Maltamo, M., Peltola, H., Packalén, P., 2012. Exploring horizontal area-based metrics to discriminate the spatial pattern of trees and need for first thinning using airborne laser scanning. Forestry 85, 305–314. Popescu, S.C., 2007. Estimating biomass of individual pine trees using airborne lidar. Biomass Bioenergy 31, 646–655. Popescu, S.C., Hauglin, M., 2014. Estimation of biomass components by airborne laser scanning. In: Maltamo, M., Næsset, E., Vauhkonen, J., (Eds.). Forestry Applications of Airborne Laser Scanning – Concepts and Case Studies. Managing Forest Ecosystems. vol. 27, Springer, The Netherlands, pp. 157-175 http:// dx.doi.org/10.1007/978-94-017-8663-8_8 Preparata, F.P., Shamos, M.I., 1985. Computational Geometry: An Introduction. Springer-Verlag, New York, USA, p. 390. R Core Team, 2012. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0, . Reitberger, J., Schnörr, C., Krzystek, P., Stilla, U., 2009. 3D segmentation of single trees exploiting full waveform lidar data. ISPRS J. Photogramm. Remote Sens. 64, 561–574. Rentsch, M., Krismann, A., Krzystek, P., 2011. Extraction of non-forest trees for biomass assessment based on airborne and terrestrial LiDAR data. In: Stilla, U., Rottensteiner, F., Mayer, H., Jutzi, B., Butenuth, M., (Eds.). Photogrammetric
Image Analysis 2011. Lecture Notes in Computer Science. vol. 6952. pp. 121– 132. Riaño, D., Meier, E., Allgower, B., Chuvieco, E., Ustin, S.L., 2003. Modeling airborne laser scanning data for the spatial generation of critical forest parameters in fire behavior modeling. Remote Sens. Environ. 86, 177–186. Riaño, D., Chuvieco, E., Condés, S., González-Matesanz, J., Ustin, S.L., 2004. Generation of crown bulk density for Pinus sylvestris L. from lidar. Remote Sens. Environ. 92, 345–352. Seidel, D., Leuschner, C., Scherber, C., Beyer, F., Wommelsdorf, T., Cashman, M.J., Fehrmann, L., 2013. The relationship between tree species richness, canopy space exploration and productivity in a temperate broad-leaf mixed forest. Forest Ecol. Manage. 310, 366–374. Shinozaki, K., Yoda, K., Hozumi, K., Kira, T., 1964a. A quantitative analysis of plant form: the pipe model theory. I. Basic analyses. Jpn. J. Ecol. 14, 97–105. Shinozaki, K., Yoda, K., Hozumi, K., Kira, T., 1964b. A quantitative analysis of plant form: the pipe model theory. II. Further evidence of the theory and its application in forest ecology. Jpn. J. Ecol. 14, 133–139. Tang, S., Dong, P., Buckles, B.P., 2013. Three-dimensional surface reconstruction of tree canopy from lidar point clouds using a region-based level set method. Int. J. Remote Sens. 34, 1373–1385. Vauhkonen, J., 2010a. Estimating single-tree attributes by airborne laser scanning: methods based on computational geometry of the 3-D point data. Doctoral Thesis Summary. University of Eastern Finland, Faculty of Science and Forestry. Dissertationes Forestales 104. Vauhkonen, J., 2010b. Estimating crown base height for scots pine by means of the 3-D geometry of airborne laser scanning data. Int. J. Remote Sens. 31, 1213– 1226. Vauhkonen, J., Korpela, I., Maltamo, M., Tokola, T., 2010a. Imputation of single-tree attributes using airborne laser scanning-based height, intensity, and alpha shape metrics. Remote Sens. Environ. 114, 1263–1276. Vauhkonen, J., Seppänen, A., Packalén, P., Tokola, T., 2012. Improving speciesspecific plot volume estimates based on airborne laser scanning and image data using alpha shape metrics and balanced field data. Remote Sens. Environ. 124, 534–541. Vauhkonen, J., Tokola, T., Maltamo, M., Packalén, P., 2010b. Applied 3D texture features in ALS-based forest inventory. Eur. J. Forest Res. 129, 803–811. Vauhkonen, J., Tokola, T., Maltamo, M., Packalén, P., 2008. Effects of pulse density on predicting characteristics of individual trees of Scandinavian commercial species using alpha shape metrics based on airborne laser scanning data. Can. J. Remote Sens. 34, S441–S459. Vauhkonen, J., Tokola, T., Packalén, P., Maltamo, M., 2009. Identification of Scandinavian commercial species of individual trees from airborne laser scanning data using alpha shape metrics. Forest Sci. 55, 37–47. Vestjordet, E., 1968. Volum av nyttbart virke hos gran og furu basert på relativ høyde og diameter i brysthøyde eller ved 2,5 m fra stubbeavskjær. Meddelelser fra Det norske Skogforsøksvesen 25, 411–549 (In Norwegian with an English summary). White, J.F., Gould, S.J., 1965. Interpretation of the coefficient in the allometric equation. Am. Naturalist 99 (904), 5–18. Woods, M., Pitt, D., Penner, M., Lim, K., Nesbitt, D., Etheridge, D., Treitz, P., 2011. Operational implementation of a LiDAR inventory in Boreal Ontario. Forestry Chronicle 87, 512–528. Yao, W., Krzystek, P., Heurich, M., 2012. Tree species classification and estimation of stem volume and DBH based on single tree extraction by exploiting airborne full-waveform LiDAR data. Remote Sens. Environ. 123, 368–380. Zolkos, S.G., Goetz, S.J., Dubayah, R., 2013. A meta-analysis of terrestrial aboveground biomass estimation using lidar remote sensing. Remote Sens. Environ. 128, 289–298. Zomorodian, A.J., 2005. Topology for computing. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, U.K, p. 258.
Contents lists available at ScienceDirect
ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs
Deriving airborne laser scanning based computational canopy volume for forest biomass and allometry studies Jari Vauhkonen a,⇑, Erik Næsset b, Terje Gobakken b a b
University of Helsinki, Department of Forest Sciences, Latokartanonkaari 7, P.O. Box 27, FI-00014 University of Helsinki, Finland Norwegian University of Life Sciences, Department of Ecology and Natural Resource Management, Høgskoleveien 12, P.O. Box 5003, NO-1432 Ås, Norway
a r t i c l e
i n f o
Article history: Received 10 January 2014 Received in revised form 19 June 2014 Accepted 2 July 2014
Keywords: Light Detection and Ranging (LiDAR) Forest inventory Tree allometry Delaunay triangulation Alpha shape Simplicial homomorphism Persistent homology
a b s t r a c t A computational canopy volume (CCV) based on airborne laser scanning (ALS) data is proposed to improve predictions of forest biomass and other related attributes like stem volume and basal area. An approach to derive the CCV based on computational geometry, topological connectivity and numerical optimization was tested with sparse-density, plot-level ALS data acquired from 40 field sample plots of 500–1000 m2 located in a boreal forest in Norway. The CCV had a high correspondence with the biomass attributes considered when derived from optimized filtrations, i.e. ordered sets of simplices belonging to the triangulations based on the point data. Coefficients of determination (R2) between the CCV and total above-ground biomass, canopy biomass, stem volume, and basal area were 0.88–0.89, 0.89, 0.83– 0.97, and 0.88–0.92, respectively, depending on the applied filtration. The magnitude of the required filtration was found to increase according to an increasing basal area, which indicated a possibility to predict this magnitude by means of ALS-based height and density metrics. A simple prediction model provided CCVs which had R2 of 0.77–0.90 with the aforementioned forest attributes. The derived CCVs always produced complementary information and were mainly able to improve the predictions of forest biomass relative to models based on the height and density metrics, yet only by 0–1.9 percentage points in terms of relative root mean squared error. Possibilities to improve the CCVs by a further analysis of topological persistence are discussed. Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
1. Introduction Carbon accounting systems and policies promoting the use of renewable energy have drawn an increasing attention toward extensive yet detailed mapping of forest biomass. Wall-to-wall estimation and mapping of forest biomass calls for various airborne and satellite remote sensing (RS) based solutions operated at local to national scales. The most recent overviews of forest biomass assessments based on RS data are provided by Koch (2010), Zolkos et al. (2013) and Popescu and Hauglin (2014). RS of forest biomass is based on indirect relationships between remotely sensed and field measured attributes. An increasingly popular approach to provide auxiliary RS data is airborne laser scanning (ALS) that produces three-dimensional (3D) vegetation height profiles from which various forest biophysical properties such as mean height (Nilsson, 1996; Næsset, 1997a; Magnussen and Boudewyn, 1998; Magnussen et al., 1999), basal area (Means ⇑ Corresponding author. Present address: University of Eastern Finland, Yliopistokatu 7, P.O. Box 111, FI-80101 Joensuu, Finland. E-mail address: jari.vauhkonen@uef.fi (J. Vauhkonen).
et al., 2000), stand volume (Næsset, 1997b; Means et al., 2000), and stem size distribution (Gobakken and Næsset, 2004) can be derived. Since the tree stems constitute the single most important forest above-ground biomass (AGB) component, these findings even apply to AGB predictions. However, ALS also allows assessing individual biomass components (Popescu, 2007; Næsset and Gobakken, 2008; Kotamaa et al., 2010; Hauglin et al., 2012), and may in fact produce more detailed information with respect to canopy biomass (i.e. branches and foliage) than conventional methods based on field measurement and modeling, as verified by destructive sampling (Hauglin et al., 2013; Kankare et al., 2013). ALS-assisted large-area assessments most often rely on areabased analyses (e.g. Næsset, 2002) of the distributions of height values extracted from sparse-density data (<1 pulses per m2; e.g. Næsset, 2007; Maltamo et al., 2011; Woods et al., 2011). Notably, such analyses employ only one of the three dimensions available in the data. Information on the horizontal distribution of forest patches may be provided by analyzing the neighborhoods of the computation units such as grid cells or segments, or their internal variation as depicted by features derived from interpolated canopy height models (Hyyppä et al., 2012; Pippuri et al., 2012; Packalén
http://dx.doi.org/10.1016/j.isprsjprs.2014.07.001 0924-2716/Ó 2014 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.
58
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
et al., 2013). Yet, the 3D information in the point data is usually ignored. In light of the theories on tree allometry and stem development (Shinozaki et al., 1964a,b), the method development should ideally aim at quantifying foliage mass due to its strong allometric links with other tree attributes. Optimizing the use of the existing predictor space in various machine learning frameworks has however been a more popular approach than expanding it with such an allometric reasoning. Earlier studies carried out at the individual tree level (Vauhkonen, 2010a) propose triangulations of the 3D point data to preserve the three-dimensional properties and produce foliage mass equivalent information for various analyses. A particularly interesting approach is the concept of 3D alpha shapes (Edelsbrunner and Mücke, 1994), or a-shapes, in which a predefined parameter a is used as a size-criterion to determine the level of detail in the obtained triangulation. Volume and complexity metrics based on a-shapes were used as proxy variables for predicting the species and stem attributes of individual trees by Vauhkonen et al. (2008, 2009, 2010a,b), Reitberger et al. (2009), Rentsch et al. (2011), and Yao et al. (2012), while Vauhkonen et al. (2012) provided further insight on using these metrics to improve area-based estimations of species-specific plot volumes. An important step in the analyses based on a-shapes is the selection of a, for which the selected value should enable separation of void spaces from those populated by canopy biomass. Vauhkonen et al. (2008, 2009, 2010a,b, 2012) evaluated a range of a0 s, but rather than selecting one particular value, they used metrics quantifying the difference of the shape obtained with the range of a0 s to the convex hull of the point data, which corresponds to a-shapes with a ? 1. The other studies mentioned in the previous paragraph do not explain the selection of a. Korhonen et al. (2013), attempting to model individual tree crown volume, selected a quasi-optimal a so that the resulting a-shape enclosed the point data within a single connected component. However, they found the volume based on the convex hull to have a closer correspondence with the field-measured reference crown volume, which also included void space due to the applied measurement principle. Most interestingly, Vauhkonen (2010b) proposed iterating over a sequence of a’s and by that vertically delineating individual tree crowns, which was further tested in an area-based application by Maltamo et al. (2010). This approach has a formal basis in topological connectivity and is particularly related to filtration of simplicial complexes (e.g. Delfinado and Edelsbrunner, 1995). Triangulating ALS point data corresponds to subdividing the underlying space of the points into simplices, which results in weighted simplicial complexes (e.g. Edelsbrunner, 2011), with weights quantifying the (empty) space delimited by the points. We hypothesize that filtering these complexes according to the weights could provide a means to reconstruct the volume populated by biomass and separate that volume from the canopy voids (Fig. 1). Thus, the purpose of this study was to develop a filtering and optimization procedure for deriving a ‘‘biomass-optimal’’ computational canopy volume (CCV) from the triangulations of plot-level ALS point clouds. Two criteria were employed to filter the triangulations of the point data, and the resulting filtration parameters and CCVs were evaluated on field data with respect to the total and canopy biomass, stem volume and basal area as target biophysical properties. 2. Methods 2.1. Study area and data The study area is located in the municipality of Aurskog-Høland, in southeast Norway (59°500 N, 11°340 E, 120–390 m a.s.l.). The dominant tree species are Norway spruce (Picea abies (L.) Karst.)
Fig. 1. Derivation of the CCV from the ALS point data of an example plot (schematic). Each step constitutes the input of the next consecutive step.
and Scots pine (Pinus sylvestris L.). Younger stands have a larger portion of deciduous species than mature stands. Birch (Betula pubescens Ehrh.) is the dominant deciduous species. For further details, see Magnussen et al. (2010), who studied the same area.
2.1.1. Field plots Altogether 40 circular field plots (36 of 1000 m2 and 4 of 500 m2) were measured during the fall of 2007 and winter of 2008. These plots were located along five systematically spaced lines, but the exact plot locations within the lines were determined subjectively to obtain a sample of 20 spruce dominated and 20 pine dominated stands. Furthermore, 5 and 15 plots covered young and old growth for each of the two species, respectively. Each plot center was geographically referenced with a survey grade Global Positioning System and Global Navigation Satellite System receiver (Topcon LegacyE) applying differential post-processing by the Pinnacle software package version 1.00 (Anon, 1999). Based on the positional standard errors reported by Pinnacle, the estimated accuracy of the planimetric plot coordinates ranged from 0.1 to 0.35 m, with an average of 0.12 m. All trees with a diameter at breast height (DBH) >5 cm were measured for the DBH and species. For the dominant tree species, the height of two trees in each of five diameter classes – covering the entire range of diameters – was measured using a Vertex hypsometer. Up to five trees of each minor species, selected with a probability proportional to stem basal area, were measured for
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
height. The total number of height measurements per plot ranged from 11 to 20 with an average of 17. Six nonlinear height–diameter regression models were developed and the missing tree heights were predicted with these models. The models were species (Norway spruce, Scots pine, and deciduous spp.) site productivity (high, low) specific. The plot-level forest attributes considered in this study were AGB and foliage biomass (FB), the latter comprising branches and foliage, basal area (G), and stem volume (V). Single-tree biomass fractions were estimated for all trees using Marklund’s (1988) biomass models based on the DBH, height, and tree species. The biomass models for birch were used for all deciduous trees. The tree-level biomass fractions were summed to AGB and FB. Plotlevel G was computed based on summing from the diameter measurements. Individual stem volumes were estimated using speciesspecific models based on the DBH and height (Braastad, 1966; Brantseg, 1967; Vestjordet, 1968), and plot-level V was computed as the sum of the individual stem volumes. See Magnussen et al. (2010) for further details regarding the measurements and computations. 2.1.2. ALS data The ALS data were acquired under leaf-on conditions 8–10 June 2005 with a supplementary flight conducted on 6 September 2005 to fill in a minor gap in the data. The acquisition was carried out using a fixed-wing Piper PA31-310 aircraft flying at an altitude of approximately 1850 m above ground level at a speed of 75 ms1. The Optech ALTM 3100 scanner was operated with a pulse repetition frequency of 50 kHz. The scan frequency was 71 Hz, resulting in a point density of approximately 0.7 m2 on the ground. The maximum scan angle was 15° but pulses emitted at an angle >13° were discarded during the subsequent data processing. Standard pre-processing including strip adjustment and ground classification was performed by the contractor (see Maltamo et al. (2010) for details). As a result of the pre-processing, the echo heights were normalized with respect to the ground level. The main results of this study were based on all echoes reflected from the areas of the field plots above a ground threshold of 0.5 m. In addition, only the first echoes (i.e. ‘only’ and ‘first-of-many’ of up to four echoes recorded per pulse) extracted employing a ground threshold of 2 m were included in a benchmark calculation. 2.2. An overview of the analysis The areas of interest (AOIs) of our approach were the individual plots (or grid cells in the case of prediction) corresponding to typical area-based analyses of ALS data (Næsset, 2002). We defined the CCV of an individual AOI as the sum of the volumes of individual simplices (tetrahedra) belonging to a triangulated ALS point cloud of that particular AOI (Fig. 1). The triangulation and the number of potential filtrations (Steps 2–3 of Fig. 1) were determined solely based on the topological properties of the point data (see Section 2.3). The fundamental problem was thus to determine a proper filtration (Step 4 of Fig. 1) such that the simplices considered in the computation of the CCV represented the areas populated by the forest biomass. As hypothesized in Section 1, we assumed that an AOI-specific filtration was required to focus the computation on the simplices representing forest biomass and to exclude those representing empty space. According to the reasoning of the previous paragraph and in order to apply the proposed approach, one should be able to filter the triangulations of all the AOIs considered such that the CCVs produced were related to the field-observed biomasses as closely as possible. To define such filtrations, we applied the practical two-stage procedure proposed by Næsset (2002), in which (i) models to predict the forest attributes of interest for the individual AOIs
59
are fit based on a set of training field plots; and (ii) the resulting models are applied to all the AOIs of the entire inventory area to produce wall-to-wall predictions. In our case, the corresponding stages were: (i) to optimize the filtrations with respect to the selected forest biophysical attributes observed from the training plots and to model the parameters producing these optimal filtrations; and (ii) to predict corresponding filtrations to all the AOIs based on information extracted from the point data. Lacking independent validation data, we followed the two-stage approach with cross-validation, i.e. predicted step (ii) to the training field plots and evaluated this prediction as if these plots represented totally independent data with unknown field observations. 1. Define filtering rules to be able to adjust the level of detail in the triangulations (Section 2.3). 2. Optimize the degree of filtering such that the CCVs extracted from the filtrations were in a linear relationship with the reference biomass values observed on the corresponding training field plots (Section 2.4). 3. Analyze the filtration parameters obtained in Step 2 with respect to forest biomass and structural attributes (Section 2.5). The main aim of this analysis was to assess the conditions under which the filtration parameters could be estimated in a wall-towall ALS inventory based on a set of training field plots. 4. Predict the filtration parameters and extract the CCV for the AOIs of the entire inventory area (Section 2.6). In the absence of a prediction grid or independent validation data, we used the set of training plots to mimic the AOIs of the inventory area, and evaluated the predictions obtained in the same data used in the earlier steps following the principle of cross-validation (Section 2.7). 2.3. Deriving triangulations and their filtrations The terminology and notation in the following are adapted from Delfinado and Edelsbrunner (1995), to which the reader is referred for additional details. A triangulation T of a finite set of points S is a subdivision of the convex hull of S into simplices according to the rules of the particular triangulation (see, e.g., Preparata and Shamos, 1985). T is a simplicial complex that is subsettable to a finite number of subcomplexes Ci, i = 1, 2, . . ., s, where i indexes the position and s is the number of different subcomplexes in a list of potential filter values F defining which simplices belong to the subcomplex (see the next paragraph). When F is ordered (here ascendingly), the sequence Ci constitutes a filtration, in which each complex is a sub-complex of its successor, and VCi is the volume of the underlying space of subcomplex Ci according to this filtration. The point data of each individual plot were triangulated separately and filtered by two alternative criteria. First, the triangulations were filtered with respect to the tetrahedral volume (v), i.e. F was an ascending sequence of the individual v’s and Ci included the tetrahedra with v 6 F[i]. Second, the filtration corresponded to the family of 3D a-shapes (Edelsbrunner and Mücke, 1994), each of which consist of those simplices of T, which have an empty circumscribing sphere with a squared radius smaller than a. Although an a-shape is defined for every non-negative real a, there are only finitely many different a-complexes (see Delfinado and Edelsbrunner, 1995, p. 782) and F was thus an ascending sequence of those a-values which altered the a-shape. All computations were based on Delaunay triangulations and implemented with C++ and an open-source Computational Geometry Algorithms Library (Pion and Teillaud, 2013; Da and Yvinec, 2013). Fig. 2 gives an example and illustrates the differences between the two types of filtrations. The filtered v-complexes constitute of individual tetrahedra simply sorted according to their volume, whereas the a-complexes are constructed via more delicate
60
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
analysis on simplices exposed to the a parameter. The CCV extracted from both the filtration types seem to be related to the canopy volume originating from the actual locations of the trees, but also on the degree of the filtering applied (Fig. 2).
AGB-v, FB-v, AGB-a, and FB-a are used to denote these op – filter pairs.
2.4. Optimizing the filtrations with respect to forest biomass
The analyses of Section 2.3 provided a range of filtration parameters and thus CCVs that were in a quasi-optimal linear relationship with forest biomass in at least one out of the 100 solutions. The correspondence of the CCVs and the forest attributes of interest was further quantified by fitting an allometric power equation (White and Gould, 1965):
To optimize the filtration parameters, all the studied plots were considered jointly. For convenience, let V C i ;p in the following text denote the CCV based on a filtration of the ALS point data of plot p. The optimization was based on minimizing the residual errors of a linear regression model:
op ¼ b1 þ b2 V C i ;p þ ep ;
2.5. Statistical analyses
k
yp ¼ bxp þ ep ;
ð2Þ
ð1Þ
where op was the reference biomass, b1 and b2 the model parameters, and ep the residual error of plot p. Both AGB and FB were tested as op, since the CCV based on ALS point data was expected to represent canopy biomass (branches and foliage) rather than that originating from the stems. The initial values of b1 and b1 were obtained by fitting the model to V C i ;p with i selected randomly for each plot. In the optimization, V C i ;q of plot q which produced the largest residual error was altered to V C i 10;q in successive rounds of 500 iterations, until either the error, determined as 1 – the obtained coefficient of determination (R2), decreased below 0.0009 or did not change during the last round of iterations. The whole optimization procedure was repeated 100 times with op as either AGB or FB, and Ci filtered by v or a. In the following text,
where yp was the considered forest attribute (either AGB, FB, G, or V), xp the optimized CCV, b and k the model parameters, and ep the residual error of plot p. The model fitting was carried out with the nls function of R (R Core Team, 2012). The ranges of the filtration parameters belonging to the optimal solutions were analyzed. The analysis was carried out with respect to both the absolute values and percentiles of the distributions of these values. The filtration parameters were examined with respect to forest structural attributes, of which the plot basal area had the highest correlation with them. The sensitivity of the CCV to the filtration parameters was analyzed by randomly selecting 10,000 realizations of the filtration parameters from the optimal range and recording the obtained R2 between the CCV and the predicted biomass attributes.
Fig. 2. Spatial distribution of a filtered 3D triangulation of an example plot. Circles depict the relative sizes and locations of the centroids of the tetrahedra belonging either to the whole triangulation (subplot A; V C i = 11251.6 m3) or filtrations a < 21.0 m3 (B; V Ci = 5388.6 m3), v < 3.7 m3 (C; V C i = 1538.9 m3), or a < 11.3 m3 (D; V C i = 2620.8 m3). The filter values in B–D correspond to the maximum a, median v, and median a, respectively, of this plot. The red crosses depict the relative sizes and positions of the field-measured trees and the outer circle indicates the plot border. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
2.6. Predicting the filtration parameters and biomass attributes The area-based two-stage procedure proposed by Næsset (2002) was tested to evaluate the feasibility of predicting both the filtration parameters and biomass attributes in a practical ALS-based forest inventory. The filtration parameters were correlated against the most common ALS-based predictor variables (cf. Næsset, 2002), i.e. the mean and standard deviation of the height values, the proportion of echoes above 2 m, and the 5th, 10th, 20th, . . ., 90th, and 95th percentiles and the corresponding proportional densities of the ALS-based canopy height distribution, calculated according to Korhonen et al. (2008, pp. 502–503). The CCV was then extracted according to the filtration parameters modeled by the ALS variables. Further, combinations of CCV with other ALS variables were inserted in Eq. (2) to find out the added value of the CCV over other potential predictors. The x variables tested in Eq. (2) were (1) CCV; (2) a height (H) and a density (D) variable derived from the canopy height distribution (see above); and (3) H + D + CCV. The plot-level mean height and the proportion of echoes above 2 m, both based on the first echoes, were found yield the highest correlations with the predicted biomass attributes and were thus used in all related experiments as H and D, respectively. 2.7. Evaluation and performance measures The correspondence of the estimates with the reference data was evaluated by graphical assessments and in terms of R2, root mean squared error (RMSE) and mean difference (bias), the latter two calculated according to:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 ^ i¼1 ðyi yi Þ RMSE ¼ n Pn bias ¼
^ yi Þ n
i¼1 ðyi
ð3Þ
ð4Þ
^i and y are the prewhere n is the number of observations, and y dicted and reference values, respectively. The relative RMSE and bias were calculated by dividing the absolute RMSE and bias values by the mean value of the reference attribute. A paired t-test was used to test the significance of the mean difference between the observation vectors. In the following text, the term ‘significant’ refers to the statistical significance of this difference at the 95% confidence level. 3. Results The following results are based on using all echoes and a ground threshold of 0.5 m. A brief comparison with the first return data and a threshold of 2 m is given as the last paragraph of this section. The optimization (Section 2.3) required 6100 or 5100 iterations, on average, when the filtrations were based on v or a, respectively. The average R2 of Eq. (1) was 0.99 or 0.92 when op was either AGB or FB, respectively. Altogether 88, 85, 92, and 90 of the 100 iterations based on AGB-v, FB-v, AGB-a, and FB-a, respectively, had R2 P 0.985. The following results are based only on those iterations. The ranges of the filtration parameters included in the optimal solutions are shown in Fig. 3. The ranges varied between the plots and their mean values increased significantly according to an increasing plot basal area (R2 = 0.42–0.61 depending on the filtration). Further, the length of the range was significantly correlated with plot basal area (R2 = 0.12–0.30) in the filtrations based on v. The relationship with the basal area was more pronounced when examined in terms of percentiles (Fig. 3). The mean values (R2 = 0.74–0.88) and ranges (R2 = 0.20–0.44) of all the filtrations
61
were significantly correlated with the basal areas. The betweenplot variation in the data density (0.94–3.66 echoes/m2) correlated insignificantly with the absolute values, and the correlation with the percentile values was only slightly significant (R2 = 0.09– 0.13). The mean values of the filtration parameters were not significantly different whether the optimization was based on either AGB or FB (see also Fig. 3). The CCVs extracted with the mean values of the filtration parameters were in 1:1 relationship (R2 = 1.0 based on all the filtrations) with those biomass attributes that were considered in the optimization. The relationships between these CCVs and the other forest attributes considered are shown in Fig. 4. The CCV derived from the optimized filtrations had R2 of 0.88–0.89, 0.89, 0.83– 0.97, and 0.88–0.92 with AGB, FB, V, and G, respectively, depending on the applied filtration. The CCVs derived from the filtrations based on v were generally larger than those based on a, but there were only minor differences between the relationships with the CCVs and the forest attributes (Fig. 4). The sensitivity analysis of deriving the CCVs with filtration parameters sampled from various ranges around the means indicated that these had to be close to the quasi-optimal values to obtain high accuracies. A worst-case R2 > 0.76 could be obtained for all the considered forest attributes, when the filtration parameters were at most within one standard deviation around the mean (Fig. 5). The correspondence was lower when sampling around two standard deviations (R2 = 0.67–0.96) or the whole range (R2 = 0.60– 0.96). However, the average result obtained from the 10,000 realizations always yielded an R2 > 0.73. Of the considered forest attributes, G had the lowest R2 values. The filtrations based on a generally provided higher R2 values (Fig. 5). The filtration parameters were found to be significantly correlated with several ALS-based height and density metrics. As could be expected based on Fig. 3, the correlations between these metrics were higher with the percentiles of the filtration parameter distributions than the absolute values. Several individual height and density variables had R2 > 0.6 with the percentiles of the filtration parameter distributions (detailed results not shown). The a-values were generally found to have higher correlations and better predictability than v, and thus the following analyses are presented for a-based filtrations only. A model 33.9 + 27.2 ln (H10) + 47.4 D2, where ln(H10) is the natural logarithm of the 10th height percentile and D as defined earlier (Section 2.4), had R2 = 0.83 for predicting the percentile of the a-based filtration parameter. With such a model, altogether 20 or 32 of the 40 predictions laid within one or two standard deviations of the reference values (Fig. 6), respectively. The CCV derived with the predicted a had R2 of 0.90, 0.77, 0.90, and 0.84 with AGB, FB, V, and G, respectively (Table 1). The obtained CCVs always had a higher R2 than any individual height and density metric, but H + D always produced a better fit than CCV with respect to both R2 and RMSE (Table 1). However, the performance of H + D could be improved by including the CCV as an additional variable in the model in all other cases except FB. The improvement of the other predictions was in order of 0.6–1.9 percentage points, least for AGB and most for G. The residuals of all models were distributed normally and did not show trends toward forest structural attributes or the between-plot variation in the pulse density. None of the mean differences were statistically significant. When repeating the experiment in the data set including only the first echoes extracted above a higher ground threshold, the only significant differences were observed with respect to the running times of the optimization algorithm and the accuracy of the results related to FB. The optimization required less iterations (2800–3500), which were also faster due to the lower number of simplices considered. The proportion of FB was presented more
62
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Fig. 3. The ranges (lines) and means (dots) of the absolute values (left) and percentiles (right) of the optimized filtration parameters of AGB-v (top), FB-v, AGB-a, and FB-a (bottom). The x-axes are ordered according to the plot basal area, which increases from left (13.6 m2/ha) to right (42.0 m2/ha). The maximum v beyond the scale of the two topmost graphs was 102.3 m3 produced by a single outlier.
accurately in all performance measures – most remarkably, the RMSE obtained by the H + D + CCV model was approximately one percentage point lower than that reported in Table 1. Otherwise, the results based on the two point data sets showed decimal-level differences in terms of the performance measures presented above.
4. Discussion The results of this study indicate that the developed new feature type, computational canopy volume, has a high potential to improve ALS-based estimations of forest biomass attributes. This is a logical outcome, since the canopy volume is directly related to the amount of foliage and canopy bulk density (cf. Riaño et al., 2003, 2004). The approach to derive the CCV adopted here was based on triangulating the 3D point data and filtering the triangulations for empty space. When the degree of the filtration was determined in an optimization, the derived CCVs were more closely related to the biomass attributes that were considered. When the filtration parameters were predicted by means of a simple
linear regression model, the accuracies were lower and approximately at the same level than those obtained with simple height and density metrics. However, the accuracy levels are overall higher than those typically reported in similar boreal forest conditions, likely owing to the larger plot size used here, which may partially explain the small improvement obtained by our approach. Yet the determination of the filtration parameters may require further considerations when applied in practical inventories (see further text), the present results are considered as an important step in the development of ALS-based auxiliary variables for forest biomass inventories and should warrant further testing of the presented approach. The approach to describe individually delineated tree crowns with volumetric properties extracted from the ALS data was first tested by Vauhkonen et al. (2008), who proposed a density of 3 pulses m2 as a minimum data requirement for such an analysis. However, the results of this study together with those presented by Vauhkonen et al. (2012) show that the operational densities of <1 m2 are feasible when triangulations of areas larger than individual tree segments are considered. Earlier, determining forest
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
63
Fig. 4. The relationships between AGB and FB (top), V (middle), and G (bottom), and CCV optimized with AGB (left) or FB (right). CCVs based on v and a are depicted by circles and crosses, respectively. The coefficients shown are based on fitting Eq. (2) to the data.
canopy volume by simply multiplying the laser-based canopy height and area, potentially applying a correction for mean canopy cover, has been proposed by Nelson et al. (1984), Næsset (1997b), Lefsky et al. (1999), Riaño et al. (2003, 2004), Chen et al. (2007), and Hollaus et al. (2009). Further, an iterative surface wrapping technique based on delineated tree crowns (Kato et al., 2009) or stacking horizontal slices of the height values (Tang et al., 2013), both operating with the point data but in 2.5D space, have been proposed. However, deriving the canopy volume from the 3D triangulations of the point data differs from the previous approaches and has certain benefits over them. A main benefit is the ability to adjust the level of detail in the triangulated point cloud by means of the filtrations, thus accounting for canopy gaps and detailed properties existing in the data.
The filtration parameters were determined in an optimization, in which the aim was to remove the empty space within the point cloud and focus on the remaining canopy structures populated by biomass. The applied optimization forced a linear relationship between the CCV and biomass. Whether the true form of the allometric relationship between these attributes was known, a corresponding model could have been fit in the optimization. However, since the correspondence of the CCV with the actual canopy volume would in any case depend on the definition and measurement of these attributes (see Korhonen et al., 2013), the developed CCV is seen as a computational feature complementing the ALS-based height and density information rather than an estimate of canopy volume with a strictly allometric definition. For example, our experiment with the first echoes and a higher ground
64
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Fig. 5. Sensitivity of the filtration parameters (left: v, right: a) to R2 between the extracted CCV and the forest attributes considered. In the analysis, 10,000 realizations of the filtrations parameters were sampled from within the mean of the optimal solutions ± one (1sd) or two (2sd) standard deviations or the full range.
Fig. 6. Predicted vs. reference a (black dots). One and two standard deviations of the reference values are depicted by the solid and broken lines, respectively.
threshold showed that the CCV can only be determined for a certain part of the canopy and still be linked accurately with the biomass attributes. Two types of filtrations and optimization schemes were tested for deriving the CCV. It should be noted that the CCVs and the filtrations overall could likely be related to several other forest biophysical properties than the biomass-based variables considered. For example, the spatial distribution of the tetrahedra could potentially be related to the spatial pattern of the trees (see Fig. 2), the prediction of which has been found difficult based on both the area- and tree-level analyses of ALS data (cf. Packalén et al., 2013). The results promote the usability of the a-shapes (Edelsbrunner and Mücke, 1994), which constitute an interesting geometric modeling tool also based on earlier related studies (e.g. Vauhkonen et al., 2012). The approach has a profound mathematical formalism on topological connectivity rather than the alternative approach based on the tetrahedral volume, which was
considered as a straightforward and potentially more easily implementable approach. Except for the optimization, the implementations were based on efficient data structures (Pion and Teillaud, 2013; Da and Yvinec, 2013) and are computationally feasible for practical applications. Compared to FB, the use of AGB in the optimization provided better results when predicting the other forest attributes (G, V). The use of the first echo data improved the predictions related to FB, yet AGB was preferred also in those situations. Although it was initially assumed that the ALS data would rather depict the canopy biomass, this observation may be related to deriving the reference values by models, which are typically found more inaccurate with respect to the individual biomass components than the total AGB. Further, it is acknowledged that species-specific differences should be accounted for to improve the estimations, since the amount and allocation of foliage and branches differ between the species considered (see also Næsset and Gobakken, 2008). Since the largest biomasses in our data originated from spruce stands and smallest from pine stands, a plot-level analysis of species would not have been independent of the tree size. Besides species, a future topic to be studied is the plot size, which was purposely large to allow a comprehensive analysis of the studied phenomenon. It is reasonable to assume that different species, plot sizes and ALS sensors used would have a similar effect on the developed features as earlier described in the literature for other predictor variables (e.g. Næsset and Gobakken, 2008; Gobakken and Næsset, 2008). A proper selection of the filtration parameters turned out to be highly sensitive toward the obtained results (cf. Figs. 5 and 6). A simple approach based on relating the a-parameter to the height and density metrics of ALS data was tested as a means to determine the degree of filtration for a practical prediction of forest biomass. The approach was found to work, yet producing a suboptimal result. Improvements could potentially be obtained by further examining the inherent algebraic properties of the triangulations such as simplicial homomorphism and topological persistence (Zomorodian, 2005). For example, Martynov (2008), computing persistent homology groups at different spatial resolutions and assuming that important features and structures persisted over a wide range of scales, de-noised geometric models of power line towers reconstructed from the ALS point data. Yet such an analysis is acknowledged to be considerably simpler with manmade objects than vegetation, which has natural structural variations, computing persistence homology from vegetation point clouds is identified as an interesting topic for future studies. Besides ALS-based large-area inventories, the proposed approach could have a high practical importance toward applications requiring geometrically explicit models of forest structure, which are currently constructed by field-intensive terrestrial laser scanning (TLS). Such applications include mapping and modeling of vegetation clumping and effective leaf area index (LAI, Hopkinson et al., 2013), vertical distribution of photosynthetically active radiation (PAR) and the fraction of PAR (fPAR, van Leeuwen et al., 2013), light transmittance and the spatial distribution of foliage (Béland et al., 2014), and forest primary production (Seidel et al. 2013) and growth (Metz et al., 2013). Particularly the approach proposed here could provide a considerably higher spatial coverage when applied with the ALS data, while the input on acquiring the TLS data could be focused on collecting the optimization samples for this approach.
5. Conclusion A novel ALS-based feature type, computational canopy volume (CCV), showed potential to improve forest biomass estimations
65
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66 Table 1 Accuracies of predicting the forest attributes considered with combinations of CCV, H, and D. Attribute
Predictors
R2
RMSE
RMSE%
Bias
Bias%
AGB (t/ha)
CCV H+D H + D + CCV
0.90 0.91 0.92
18.8 17.0 16.4
15.9 14.4 13.8
0.3 0.1 0.0
0.3 0.1 0.0
FB (t/ha)
CCV H+D H + D + CCV
0.77 0.85 0.86
6.8 5.4 5.4
23.4 18.6 18.6
0.05 0.00 0.00
0.2 0.0 0.0
V (m3/ha)
CCV H+D H + D + CCV
0.90 0.91 0.92
32.7 31.1 29.4
15.8 15.1 14.2
0.4 0.5 0.03
0.2 0.2 0.0
G (m2/ha)
CCV H+D H + D + CCV
0.84 0.90 0.93
3.5 2.7 2.2
13.9 10.8 8.9
0.01 0.03 0.01
0.1 0.1 0.0
from sparse-density, area-based ALS data. The best-case R2 obtained between the CCV and total biomass, canopy biomass, stem volume, and basal area were 0.88–0.89, 0.89, 0.83–0.97, and 0.88–0.92, respectively. The R2 depended on the applied filtration and was in all cases based on optimization with forest biomass attributes. The magnitude of the required filtration increased according to an increasing basal area, and a simple prediction model with related ALS height and density metrics provided CCVs which had R2 of 0.77–0.90 with the aforementioned forest attributes. The derived CCVs mainly improved the predictions of forest biomass, yet only by 0–1.9 percentage points in terms of RMSE%, compared to models based on the conventional height and density metrics. A benchmark calculation showed that similar results could be obtained by restricting the analysis to specific echo types (here first echoes), thus potentially reducing the computational burden involved. The main bottleneck is the prediction of the filtration parameters, which could potentially be improved by further analyses of topological persistence of the triangulations and their filtrations.
Acknowledgements This study was enabled by the Visiting Researcher Grant (2012/ 2099) of the Norwegian University of Life Sciences and additionally supported by University of Helsinki Funds. The authors wish to thank Blom Geomatics (Norway) for providing and processing the ALS data. The authors also wish to thank K. M. Hauglin and V. Lien, both at the Norwegian University of Life Sciences, for the help in the fieldwork.
References Anon, 1999. Pinnacle User’s Manual. Javad Positioning Systems, San Jose, CA, p. 123. Béland, M., Baldocchi, D.D., Widlowski, J.-L., Fournier, R.A., Verstraete, M.M., 2014. On seeing the wood from the leaves and the role of voxel size in determining leaf area distribution of forests with terrestrial LiDAR. Agric. Forest Meteorol. 184, 82–97. Braastad, H., 1966. Volume tables for birch. Meddelelser fra Det norske Skogforsøksvesen 21, 265–365 (In Norwegian with an English summary). Brantseg, A., 1967. Volume functions and tables for scots pine. South Norway. Meddelelser fra Det norske Skogforsøksvesen 22, 689–739 (In Norwegian with an English summary). Chen, Q., Gong, P., Baldocchi, D., Tian, Y.Q., 2007. Estimating basal area and stem volume for individual trees from Lidar data. Photogramm. Eng. Remote Sens. 73, 1355–1365. Da, T.K.F., Yvinec, M., 2013. 3D alpha shapes. CGAL User and Reference Manual, CGAL Editorial Board, 4.3 ed.
Edelsbrunner, H., 2011. Alpha shapes – A Survey. Book Chapter Manuscript, Tessellations in the Sciences. Springer-Verlag,
66
J. Vauhkonen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 96 (2014) 57–66
Proceedings of SilviLaser 2011 – 11th International Conference on LiDAR Applications for Assessing Forest Ecosystems, 16–20 October, 2011, Hobart, Australia.
Image Analysis 2011. Lecture Notes in Computer Science. vol. 6952. pp. 121– 132. Riaño, D., Meier, E., Allgower, B., Chuvieco, E., Ustin, S.L., 2003. Modeling airborne laser scanning data for the spatial generation of critical forest parameters in fire behavior modeling. Remote Sens. Environ. 86, 177–186. Riaño, D., Chuvieco, E., Condés, S., González-Matesanz, J., Ustin, S.L., 2004. Generation of crown bulk density for Pinus sylvestris L. from lidar. Remote Sens. Environ. 92, 345–352. Seidel, D., Leuschner, C., Scherber, C., Beyer, F., Wommelsdorf, T., Cashman, M.J., Fehrmann, L., 2013. The relationship between tree species richness, canopy space exploration and productivity in a temperate broad-leaf mixed forest. Forest Ecol. Manage. 310, 366–374. Shinozaki, K., Yoda, K., Hozumi, K., Kira, T., 1964a. A quantitative analysis of plant form: the pipe model theory. I. Basic analyses. Jpn. J. Ecol. 14, 97–105. Shinozaki, K., Yoda, K., Hozumi, K., Kira, T., 1964b. A quantitative analysis of plant form: the pipe model theory. II. Further evidence of the theory and its application in forest ecology. Jpn. J. Ecol. 14, 133–139. Tang, S., Dong, P., Buckles, B.P., 2013. Three-dimensional surface reconstruction of tree canopy from lidar point clouds using a region-based level set method. Int. J. Remote Sens. 34, 1373–1385. Vauhkonen, J., 2010a. Estimating single-tree attributes by airborne laser scanning: methods based on computational geometry of the 3-D point data. Doctoral Thesis Summary. University of Eastern Finland, Faculty of Science and Forestry. Dissertationes Forestales 104