Real and simulated waveform-recording LiDAR data in juvenile boreal forest vegetation

Remote Sensing of Environment 140 (2014) 665–678

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Real and simulated waveform-recording LiDAR data in juvenile boreal forest vegetation A. Hovi ⁎, I. Korpela Department of Forest Sciences, University of Helsinki, P.O. Box 27, Helsinki FI-00014 Finland

a r t i c l e

i n f o

Article history: Received 16 April 2013 Received in revised form 2 October 2013 Accepted 6 October 2013 Available online xxxx Keywords: Monte Carlo ray tracing System waveform Calibration Radiative transfer Forest inventory Footprint

a b s t r a c t Airborne small-footprint LiDAR is replacing field measurements in regional-level forest inventories, but auxiliary fieldwork is still required for the optimal management of young stands. Waveform (WF) -recording sensors can provide a more detailed description of the vegetation than discrete return (DR) systems through accurate characterization of the backscattered signal. Furthermore, knowing the signal shape facilitates comparisons between real data and those obtained with simulation tools. We performed calibration and quantitative validation of a Monte Carlo ray tracing (MCRT) -based LiDAR simulator against real data, and used simulations and real data to study small-footprint WF-recording LiDAR for the classification of juvenile boreal forest vegetation. The simulations were based on geometric-optical models of three species: birch (Betula pendula Roth), raspberry (Rubus idaeus L.), and fireweed (Chamerion angustifolium (L.) Holub). Simulated WF features were in good agreement with the real data (differences of −19% to 11% in radiometric features, −0.23 m to 0.45 m in mean height), and relative interspecies differences were preserved. We used simulated data to study the effects of sensor parameters on the classification among the three species. An increase in footprint size improved the classification accuracy up to 0.30–0.36 m in diameter, while the emitted pulse width and the WF sampling rate showed minor effects. Finally, we used real data to classify four silviculturally important vegetation functional groups (conifers, broad-leaved trees, low vegetation (green), low vegetation (barren) + abiotic material). Classification accuracies of 68.1–86.7% (kappa 0.50–0.80) showed slight improvement compared with existing studies on DR LiDAR and passive optical data. The results of simulator validation serve as a basis for the future use of simulation models, e.g. in LiDAR survey planning or in the simulation of synthetic training data, while the other findings clarify the potential of small-footprint WF data for characterizing vegetation in intensively managed forest stands at seedling and sapling stages in the boreal region. © 2013 Elsevier Inc. All rights reserved.

1. Introduction The utility of small-footprint airborne LiDAR in producing precise quantitative estimates of forest variables (e.g. tree height, stem volume, basal area) is well demonstrated (Naesset, 2007; Naesset et al., 2004). As a result, LiDAR-based methods have begun to replace traditional field-based inventory methods (Maltamo et al., 2011; Naesset, 2004). The requirements for a particular forest inventory system are dependent on the spatial scale and the forest management scheme applied. In Finland, compartmentwise forestry is practiced, in which the basic unit is a forest stand (0.5–10 ha). Clear-cutting is normally performed and a new tree generation established at the end of the rotation period (70–150 yr). In addition, 1–3 intermediate thinnings are applied. Airborne LiDAR, combined with aerial images, is used in regional-level inventories for operational forest management planning to provide estimates of the growing stock per tree species (Packalén, ⁎ Corresponding author. Tel.: +358 9 191 58214. E-mail addresses: aarne.hovi@helsinki.fi (A. Hovi), ilkka.korpela@helsinki.fi (I. Korpela). 0034-4257/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.rse.2013.10.003

Suvanto, & Maltamo, 2009). Currently, however, the use of LiDAR and other remote sensing (RS) data in the inventory of seedling and sapling stands is limited. In particular, information relating to the need and timing of silvicultural treatments (weed-, insect-, and herbivore control, clearing to favor the future trees) must be acquired in the field. This is a subjective process based on visual estimates and/or measurements of tree height, stand density, species composition, and the spatial pattern and condition of the trees (Korhonen et al., 2013). The few studies conducted on the usage of LiDAR in the inventory of young forest stands report underestimation of tree heights and low predictive power for stand density and species composition (Korpela, Tuomola, Tokola, & Dahlin, 2008; Naesset & Bjerknes, 2001). Physical limitations, such as the small size of the trees, and the nonzero spatial extent of the LiDAR pulse (footprint, pulse width) can result in large relative errors in the LiDAR-based tree height or ground elevation measurements compared with tree size, and mixing of the signal from an individual tree with ground and other vegetation. Due to the subjectivity in determining the treatment need based on the forest variables, a direct prediction from the RS data is also an option as demonstrated for LiDAR data in mature forests (Vastaranta et al.,

666

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

2011). In addition, RS data could potentially be used to reduce the number of required field visits. For example, Korhonen et al. (2013) predicted the treatment need, using discrete return (DR) LiDAR and aerial images and could detect the cases in which the RS-based predictions were the most reliable. Auxiliary information on regeneration activities and site quality can also be used to predict, using growth models, the expected state of the stand at some point in time. Regardless of the inventory method, the dynamic nature of the vegetation in young forest stands sets requirements for the overall inventory system. Pulsed LiDAR systems send 5–10-ns-long, collimated beams of laser light towards the target at high frequency (Baltsavias, 1999). The returning signal is detected by fast electronics. The shape of the signal is a convolution of the time-dependent function of the transmitted pulse power with the target cross-section profile and the system response function (Wagner, Ulrich, Ducic, Melzer, & Studnicka, 2006). In addition, both external and internal noise contribute to the signal detected. The shape of the returning signal is equal to that of the emitted pulse if perfectly planar targets are scanned collinear with their normal. In vegetation, the signals are always extended, due to volumetric scattering. The coordinates of the illuminated target are calculated from the two-way traveltime of the pulse, accounting for the sensor position and attitude and the internal geometry of the mechanism that deflects the pulses. DR sensors detect echoes from the returning pulse (range detection) on-the-fly. Coordinates and an intensity measurement are recorded for each echo. Waveform (WF) -recording systems store a digitized amplitude sequence, the waveform, for postprocessing (Blair, Coyle, Bufton, & Harding, 1994; Schutz, Zwally, Shuman, Hancock, & DiMarzio, 2005; Ullrich & Pfenningbauer, 2011). It is evident that WF data provide a more accurate characterization of vegetation than do DR data. Research in forest applications has focused on deriving more dense point clouds by finding additional echoes (Chauve et al., 2009; Reitberger, Schnoerr, Krzystek, & Stilla, 2009) or on using radiometric WF features for classification tasks (Heinzel & Koch, 2011; Reitberger, Krzystek, & Stilla, 2008). The siteand campaign-specific conditions and the lack of DR data from the same study area complicate the assessment of additional benefits compared with the DR data. Wagner (2010) emphasized the feasibility of deriving campaign-independent, physicallybased radiometric quantities (backscattering cross-section and coefficient). This requires a solution for deconvolution, for which Gaussian decomposition (Wagner et al., 2006), use of β-splines (Roncat, Bergauer, & Pfeifer, 2011), and Fourier transforms (Jutzi & Stilla, 2006) have been proposed. Operational area-based forest inventories rely on empirical dependencies established between the forest variables and LiDAR features (Naesset, 2002; Packalén & Maltamo, 2007). Case-specific LiDAR acquisition settings and forest characteristics make it difficult to generalize the models. Less indirect estimation is feasible by means of individual tree detection, in which crowns are delineated and the prediction of stem attributes is done with allometric models that link the measured crown size, tree height, and species with the stem (Hyyppä & Inkinen, 1999; Persson, Holmgren, & Soderman, 2002). Theoretical simulations, although not completely replacing the need for empirical training and validation data, can be useful, e.g. in studying the effects of acquisition settings on the LiDAR signal observed (Holmgren, Nilsson, & Olsson, 2003) or in the planning of sampling strategies (Ene et al., 2012). In addition, simulations allow for testing interpretation algorithms and extending the empirical results to a variety of sensor and acquisition settings. In theory, it is also possible to predict forest variables from the LiDAR data through inversion of the simulation model, i.e. to produce synthetic training data. However, this requires very good absolute accuracy of the simulated quantities. Validation against real data is essential for verifying the correctness of the simulation models, and for testing their robustness to various parameters and assumptions. From the validation standpoint, DR systems are ‘black boxes’, since very little information on the emitted

pulse, echo detection algorithm, or nature of the intensity observations is normally available. In WF systems, the shape of the returning signal is observed. This facilitates comparisons between simulated and real data, although some calibration is still needed to bring the simulated quantities to match the measurement scale of the instrument. DR data acquired with terrestrial laser scanners (TLSs) can be used as references with which airborne WF data are compared or in providing the vegetation structure to be used in simulations (Bittner et al., 2012; Vierling, Xu, Eitel, & Oldow, 2012). The computational burden and the desired level of model detail are the main factors behind the selection of an appropriate simulation model. Recent studies demonstrated the use of Monte Carlo ray tracing (MCRT), combined with explicit leaf- or needle-level vegetation models, in simulation of optical RS signals even from entire forest stands (Disney, Lewis, & Saich, 2006; Disney et al., 2011). MCRT is computationally intensive, and the availability of detailed vegetation models can limit its use. On the other hand, MCRT aims at modeling the real photon paths in the canopy, and simplifying assumptions can be kept to a minimum when modeling the geometric-optical properties of vegetation. Therefore, MCRT is considered as one of the most accurate methods for radiative transfer modeling in vegetation (Widlowski et al., 2008). There have been applications for LiDAR (Disney et al., 2010; Hancock, Disney, Muller, Lewis, & Foster, 2011; Morsdorf, Nichol, Malthus, & Woodhouse, 2009). Comparisons with real data are often lacking, but have been reported for large-footprint spaceborne LiDAR (North, Rosette, Suarez, & Los, 2010). A common assumption behind the MCRT models is that the radiative transfer theory holds, i.e. light propagation in the medium can be modeled with direct rays, and the scattering is determined by the geometric-optical properties of the scene elements. Due to the coherent nature of laser light, wave phenomena (interference, diffraction) can have an effect as well (Wagner, 2010). One example is the speckle effect that causes fluctuations in the signal. In addition, ill-posed model inversion can restrict the practical applicability of MCRT. However, this problem can be tackled through the simulation of signals from different types of vegetation and by the use of a look-up-table approach (e.g. Knyazikhin, Martonchik, Myneni, Diner, & Running, 1998). The rationales of our study were twofold: first, to develop a MCRT simulation model and test its performance in reproducing smallfootprint WF-recording LiDAR signals, and second, to explore the potential of WF data in classifying tree and competing vegetation species in juvenile boreal forest vegetation. The species information is crucial to determining the need for silvicultural treatments in commercially managed seedling and sapling stands. Three specific research aims (RAs) were formulated: RA1 To perform quantitative validation of a MCRT simulation model against real data and to test the model sensitivity to variation in vegetation parameters. RA2 To use the simulation model for testing various WF sensor parameters for the classification of three selected tree and competing vegetation species. RA3 To study real WF signatures for the classification of juvenile forest vegetation representative of a large number of species. 2. Materials and methods 2.1. Outline of the experiments Our experimental framework is shown in Fig. 1. We start by describing the field reference and LiDAR datasets used (Sections 2.2 and 2.3). We then present the model used for simulating WFs from vegetation (Section 2.4). Feature extraction forms the basis for the analyses of both real and simulated data and is described in Section 2.5. Finally, we provide a short description of the analyses performed to answer the RAs. They include a comparison of the

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

667

Fig. 1. Schematic representation of the materials and methodological steps. The boxes with rounded corners illustrate processes, and the rectangles are the input and output data for these. The related section number is given in the brackets.

simulated and real data and a sensitivity analysis of the simulator (RA1, Sections 2.6.1 and 2.6.2), a study of the sensor parameter effects on classification of the three species simulated (RA2, Section 2.6.3), and a study of the real WF data for vegetation classification (RA3, Section 2.6.4).

sample. A subjective estimate of the canopy cover (0–100%) was recorded for a subset of the S12 samples. Table 1 lists the number and h values of the samples used in the analyses.

2.3. LiDAR data 2.2. Ground reference data The study area is in Hyytiälä, Finland (61°50′N, 24°20′E) and covers 2 × 6 km. Boreal coniferous forests prevail, the dominant species being Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.), with a varying degree of deciduous (mainly birches, Betula L. spp.) mixture (see also Korpela, 2006). The forests are commercially managed and mostly even-aged. Field data were collected in the summers of 2010 and 2012, shortly after the LiDAR campaigns (Section 2.3). Seedling and sapling stands were chosen to represent variation in stand age and site type (Fig. 2), and vegetation samples were collected from these. The number of stands was 14 in 2010 and 13 in 2012. The average tree height (h) was usually 1.5–5 m, but some recently established stands were also selected for sampling of low vegetation. The samples (tree base or the ground at the center of the other vegetation sample) were positioned in 3D with a Network GNSS at an accuracy of better than 3 cm. Two types of sampling schemes were used. In 2010, rectangular plots (100–200 m2) were established at subjectively selected locations. All trees were positioned and recorded for h and species. These samples are referred to as P10. Isolated vegetation samples were collected separately in 2010 and 2012. These are referred to as S10 and S12 and include trees and circular samples of low vegetation (radius (r) = 0.5 m in S10, r = 1.0 m in S12). Maximum h was measured for each

Data from two Leica ALS60 campaigns (15:00–17:00 GMT, July 19, 2010 and 17:00–19:00 GMT, July 5, 2012) were used. The data include nine datasets differing by acquisition height and WF sampling rate (Table 2). The ALS60 sensor operates at a wavelength (λ) of 1064 nm. The beam divergence is 0.15 mrad or 0.22 mrad, measured at 1/e or 1/e2, i.e. the point at which the irradiance drops to 36.8% or 13.5% from its maximum, respectively. Up to four DR echoes can be recorded per pulse. An echo consists of target coordinates and eight-bit values for the intensity and the automatic gain control (AGC). The algorithm of the intensity recording is unknown. In our data, the intensity was highly correlated with the peak amplitude of the WF from well-defined surfaces (R2 = 0.99). The returning signal is directed to a separate WF sampling unit, and the first DR echo detected triggers the WF recording. The continuous 256-ns-long (38 m) WF record contained a sample of 3.3–4.8 m that preceded the first echo. The maximum scan zenith angle was always 15°, and the flight lines were oriented in azimuths of 160° and 340° with some perpendicular strips. The WF sampling rate was 1 ns, except for two flying altitudes in the 2010 campaign, which were repeated at both 1- and 2-ns sampling. The WF data were recorded for every pulse, except in the datasets of 2010_12_1 and 2012_05_1, in which every second pulse pair was recorded, due to data transfer constraints. The groundsampling density was 0.9–5.2 WFs per m2. In the 2012 campaign, the peak power of the pulses was adjusted for constant irradiance at the

Fig. 2. Examples of the selected stands. A 12-yr-old sowed pine stand on a barren site (left), an 11-yr-old planted spruce with a naturally regenerated broad-leaved mixture on a fertile site (middle), and a recently clear-felled stand with the vegetation cover consisting of 1-m-tall fireweed, raspberry and grasses on semifertile soil (right).

668

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

Table 1 Overview of the field samples collected in 2010 and 2012 for the species represented by at least 30 pulse observations in any of the LiDAR datasets. P10 includes samples from fixed-area plots measured in 2010, while S10 and S12 depict separate, isolated vegetation samples from 2010 and 2012, respectively. The number of samples is given, as well as the minimum and maximum heights (h) in meters. Trees with h b 1.5 m were not used in the analyses. Species Trees Scots pine Norway spruce Birches European aspen European mountain-ash Gray alder Goat willow Other willows Common juniper

Scientific name

P10

S10

S12

min h

max h

Pinus sylvestris Picea abies Betula pendula, B. pubescens Populus tremula Sorbus aucuparia Alnus incana Salix caprea Salix sp. Juniperus communis

374 406 1517

153 254 282

201 118 248

1.5 1.5 1.5

7.4 6.0 7.5

49 950 47 84 106 1

20 194 79 59 48 24

56 83 75 56 86 –

1.5 1.5 1.5 1.5 1.5 1.5

4.6 4.6 7.0 4.5 4.8 2.8

–

75

102

0.6

2.5

– – – – – –

122 18 113 51 – –

120 24 – – 46 23

0.5 0.6 0.4 0.2 0.2 0.6

1.8 1.3 1.6 0.8 0.6 1.8

– – – – – –

64 20 30 94 15 1

50 28 – – – 17

0.2 0.1 0.3 0.1 0.1 0.5

0.5 0.2 0.5 0.4 0.4 0.8

Low vegetation associated with fertile sites Fireweed Chamerion angustifolium American red raspberry Rubus idaeus Western brackenfern Pteridium aquilinum Reedgrass Calamagrostis sp. Wavy hairgrass Deschampsia flexuosa Fresh grass (low) – Fresh grass (high) – Low vegetation associated with barren sites Heather Calluna vulgaris Reindeer lichens Cladina/Cladonia sp. Whortleberry Vaccinium myrtillus Lingonberry Vaccinium vitis-idaea Black crowberry Empetrum nigrum Marsh Labrador tea Rhododendron tomentosum Schreber's big red Pleurozium schreberi stem moss Mosses (misc) – Ripening grass –

–

38

–

0

0.2

– –

11 –

46 49

0 0.6

0.1 1.0

Abiotic material Stone surface Mineral soil Logging residue Sand roada

– – – –

34 21 64 67

– 19 30 152

0 0 0.1 0

0 0 0.8 0

a

– – – –

Used as reference surfaces when comparing WF features in different datasets.

ground from all four acquisition heights. Since the beam divergence in the ALS60 is constant, the footprint diameter was directly proportional to the flying altitude. Campaign-level boresight and range (R) calibrations were carried out, followed by relative strip matching for correcting the roll, pitch, and Z discrepancies. The XY accuracy of the DR data was evaluated as

better than 0.25 m, using artificial and natural features. The Z accuracy was evaluated, using GNSS samples from level surfaces. Small, cmlevel systematic shifts in Z were corrected for. The trajectory data, and thus the 3D pulse paths, were also available. The ALS60 applies the AGC that regulates the receiver gain (Korpela, 2008). It ensures that the recorded intensity values are within an eight-

Table 2 Overview of the LiDAR datasets. The name reflects the year, acquisition height in hectometers, and the WF sampling rate in nanoseconds. Dataset

Mean and (SD) of R in the field samples, m

CV of intensity due to R variation, %

Pulse repetition frequency, kHz

Laser peak power, kW

Emitted pulse width (FWHM), ns

Average footprint diameter at ground (1/e), m

WF sampling rate, ns

Density of pulses with a WF record, m−2

2010_12_1 2010_19_1 2010_19_2 2010_30_1 2010_30_2 2012_05_1 2012_10_1 2012_20_1 2012_27_1

1200 (18) 1957 (23) 1957 (21) 3022 (39) 3021 (38) 522 (17) 1029 (19) 2030 (26) 2787 (33)

3.3 2.6 2.4 2.9 2.8 7.3 4.1 2.8 2.6

173.6 119.0 119.0 83.3 83.3 152.6 99.2 58.9 44.9

3.76 6.93 6.93 7.08 7.08 0.53 1.42 3.60 5.30

7.2 7.2 7.2 10.1 10.1 7.7 10.1 10.3 10.5

0.18 0.29 0.29 0.45 0.45 0.08 0.15 0.30 0.42

1 1 2 1 2 1 1 1 1

4.7 2.7 2.8 0.9 0.9 5.2 5.2 3.1 2.1

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

bit range across the scanned area. The AGC affects the recorded signal levels. To account for this, we used a normalization model modified from that in Korpela (2008): Ical ¼

I obs 1 þ ðAGC−aÞ  b

ð1Þ

where Iobs is the raw and Ical the normalized DR intensity, parameter a is a reference AGC value to which raw values are normalized, and parameter b gives the increase in signal level per unit increase in AGC. The model was applied to the WF data by replacing Ical and Iobs with the amplitude values, and the entire WF sequence was then normalized, using Eq. (1). Since the AGC had no effect on the WF noise level, the minimum noise was first subtracted from the observed WF and added back to the normalized WF. Thus, only the noise-exceeding part of the WF was normalized. The value of parameter a was set at the average AGC value in each dataset, and the optimal values for b were searched for by minimizing the nested coefficient of variation (CV) of the WF peak amplitude data, using homogeneous targets that covered a range of surface reflectivity. The targets in both 2010 and 2012 were bitumen roof, old asphalt, fine sand, and grass. In addition, reflectance tarps (see Section 2.4.4) were used in 2010. The standard deviations (SDs) of R were small (the CV of R was within 3.2%), limiting the effect of R on intensity to below 7.3% in CV (Table 2). Here we assumed that the effect of R is relative to the power of 2.5 (Korpela, Orka, Hyyppä, Heikkinen, & Tokola, 2010). Since the AGC effect was much stronger, we omitted the intensity variation due to R in the normalization. 2.4. Simulation of WF LiDAR data from vegetation 2.4.1. Overview of the simulation approach Simulations were performed in two steps. First, the time (or R) -dependent differential scattering cross-section of the target was modeled, using MCRT (Section 2.4.2). Secondly, the cross-section was convolved with the LiDAR system WF to obtain an estimate of the recorded signal (Section 2.4.3). Sections 2.4.4–2.4.6 describe the various components needed in the simulations. The system WF and noise level

Fig. 3. Illustration of the MCRT sampling scheme. The dashed small arrows represent the populations of all possible photon paths. If a ray i intersects a scene element, a new ray is cast towards the receiver and the contribution to radiant intensity (Ii,n) is calculated. In addition, a ray with power Pi,n is generated to sample the next scattering order. Here, the ray tree escapes the scene after two interactions. Note that in reality the receiver and transmitter are at the same location in LiDAR sensors.

669

for each LiDAR campaign were estimated, using natural and artificial calibration targets (Section 2.4.4). To compare the echo height of the simulated and real data, a range detection algorithm was applied to the simulated WFs (Section 2.4.5). Representations of the geometricoptical properties of the vegetation were used as input for MCRT (Section 2.4.6). 2.4.2. Modeling the target scattering cross-section with MCRT We developed a MCRT radiative transfer model for estimating the ratio of scattered radiant intensity towards the receiver [W sr−1] to the incident irradiance at the target surface [Wm−2], i.e. the differential scattering cross-section of the target (σ, [m2 sr−1]). The idea in MCRT is to estimate the energy distribution of light by randomly sampling it with infinitely narrow rays, i.e. to provide an approximate solution for the rendering equation (Kajiya, 1986). Monte Carlo techniques are used for random sampling of all possible photon paths. For a more detailed discussion and fundamentals of MCRT, see Disney, Lewis, and North (2000). Our model works in the forward mode, i.e. the photon paths are traced from the light source to the receiver. The spatial density of the rays cast towards the scene follows the Gaussian radiation field of the laser beam. For each ray, the first intersection with a scene element is searched, and when found, a new ray is cast towards the receiver (Fig. 3). If the receiver is visible, the ray contributes to the radiation observed. The scene elements are assumed to be bi-Lambertian with known reflectance (ρ) and transmittance (τ). The contribution of ray i of scattering order n to the total radiant intensity towards the receiver is denoted by Ii,n. It is calculated from the ρ, view zenith angle (θ) in relation to the surface normal, and the incident power of the ray (Pi,n − 1) (cf. Smolander & Stenberg, 2003):

Ii;n ¼

ρ cosθ P i;n−1 π

ð2Þ

The subsequent scattering order is sampled as follows. First, a decision is made between a reflection or a transmission. Probabilities for reflection and transmission are calculated as the ratios of ρ and τ to the single scattering albedo (ρ + τ) of the scene element in question. Secondly, the vector of the scattered ray is randomly generated from the angular distribution of the photons scattered by a Lambertian surface. The power of the ray (Pi,n) is calculated by multiplying Pi,n − 1 with the single scattering albedo. Again, intersections with scene elements are tested and the contribution to the radiant intensity is calculated with

Fig. 4. Effect of receiver iFOV at ground (width of the circular area viewed by the sensor) on the simulated σ in 8 × 8-m raspberry canopies (height 1.3 m) with varying leaf area index (LAI). Each line represents the average σ simulated from 10 instances. The scan angle was 0°, scanning R 1000 m, and the footprint diameter at ground was 0.15 m.

670

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

Eq. (2), if necessary. The sampling is stopped if the ray escapes the scene, or if the predefined number of scattering orders is reached. By setting the total power of the emitted rays at 1, the output from the MCRT is directly the target σ. To account for the time dependency of the LiDAR signal, σ was binned into 0.05-ns-long intervals (corresponding to 0.75 cm in R), based on the lengths of the individual ray paths. We call this the target cross-section profile (σ(t)). The model was implemented in Java programming language. It was first validated with simple plane surfaces in which an analytical solution for σ could be obtained. To test the approach of modeling the radiant intensity with Eq. (2), more complex test scenes with circular ‘leaves’ in different orientations were generated, and the solution was compared with that of a more intuitive approach, in which no rays are sent to the sensor directly, but instead the rays Pi,n that hit a circular aperture with a certain instantaneous field of view (iFOV) contribute to the σ. Our model showed no bias, and was computationally much more efficient. The number of rays and scattering orders needed was deduced from tests, using repeated simulations with modeled birch and raspberry canopies (see Section 2.4.6). The variation in the simulated solutions decreased with the number of rays emitted per pulse, but began to saturate at 30 000. The first 10 scattering orders contributed to over 97% of the total σ, the first-order contribution being 57–88%. The branching ratio (number of subrays sent from each intersection point to sample higher order scattering) did not considerably affect the solution. Thus, we used 30 000 rays per pulse, 10 scattering orders, and a branching ratio of 1 in all simulations. The simulation times with these parameters in the vegetation modeled were 10–40 s per pulse on a standard desktop computer (Win7, 16 GB RAM, 3.33 GHz CPU). Instead of using a fixed angular iFOV, we assumed that the sensor always ‘sees’ a circular area 5 m in diameter (Fig. 3). In a raspberry canopy, the targets outside a 5-m-wide circle contribute b 3% to the total σ (Fig. 4), which is the potential bias of our approach. The iFOV of commercial LiDAR systems is often not known to users, but it can neither be very small due to mirror movement between pulse emission and transmission, nor very large due to the solar illumination that increases with the iFOV. 2.4.3. Derivation of the LiDAR signal from the scattering cross-section and the system WF The physical basis for the formation of the LiDAR signal is presented here. The final result is an equation that converts the σ(t) from the MCRT simulation into the voltage pattern recorded by the sensor (Eq. 7). The instrument emits light at power Pe(t) [W]. The instantaneous irradiance at the target [W m−2] is obtained by multiplying Pe(t) with the atmospheric attenuation factor (ηatm, [unitless]) and by dividing with the target surface area (A, [m2]). The A is assumed constant, because the receiver iFOV is much larger than the LiDAR footprint at the ground. The radiant intensity towards the receiver (I(t), [W sr−1]) is then the convolution of the irradiance with σ(t) [m2 sr−1]: Iðt Þ ¼

ηatm P e ðt Þ ⊗σ ðt Þ: A

ð3Þ

The power observed by the receiver (Po(t), [W]) is obtained from I(t) by multiplying it with ηatm and the solid angle subtended by the receiver (πD2/4R2, [sr], where D is the receiver aperture [m] and R the scanning range [m]): πD2 P o ðt Þ ¼ Iðt Þηatm 2 : 4R

ð4Þ

The voltage recorded by the instrument (Vo(t), [digital number (DN), unitless]) is obtained by convolving Po(t) with the system

response function (Γ(t), [unitless]), multiplying with the system gain (g, [W−1]), and adding a noise term (c, [DN, unitless]): V o ðt Þ ¼ g  P o ðt Þ⊗Γ ðt Þ þ c:

ð5Þ

We call the combined effect of g, Pe(t), A, ηatm, D, and Γ(t) the system waveform (S(t), [unitless]): 2

Sðt Þ ¼ g 

P e ðt Þ 2 πD η ⊗Γ ðt Þ: A atm 4

ð6Þ

By substituting S(t) into Eq. (5) and rearranging the terms, the final equation for Vo(t) becomes: V o ðt Þ ¼ Sðt Þ⊗σ ðt Þ

1 þ c: R2

ð7Þ

In the simulations, the convolution of σ(t) with S(t) (both discretized at 0.05-ns intervals) and the division with R2 were performed first. The noise added to each time discrete sample was randomly generated from a Gaussian distribution. No autocorrelation between WF samples was assumed. Finally, Vo(t) was resampled at the given WF sampling rate (0.25–4 ns), and each sample was rounded to the nearest integer value. The S(t) was assumed constant over an acquisition. We refer to the process of determining S(t) and the mean and SD of c for each acquisition as simulator calibration. It is described in Section 2.4.4. 2.4.4. Calibration using reference targets The ALS60 does not record the emitted WF, and thus, for determining S(t), we had to rely on the WFs detected from welldefined planar surfaces of known reflectance (calibration targets). If the pulse path is perpendicular to a perfectly planar surface, σ(t) simplifies to a unique value σ at some t that corresponds to the R, and σ can be derived from the target reflectance. Real surfaces are rough and the pulse is not always perpendicular to the surface. Since the maximum scan zenith angle in our data was 15°, this assumption was not strongly violated for near-level surfaces. For these, the S(t) was solved from Eq. (7) by replacing the convolution with multiplication: Sðt Þ ¼

V o ðt Þ−c 2 R : σ

ð8Þ

The calibration targets were 5 × 5-m reflectance tarps, and asphalt and fine sand surfaces. Hemispherical-directional reflectance factors (HDRFs) relative to a Spectralon® reference panel, measured in the nadir view, were used for the determination of σ for each target. The measurements were made in 2008 for the tarps (solar elevation 30°), and in 2009 for the other targets (38–43°). Three tarps (HDRF 0.14, 0.21, 0.41), asphalt (0.21), and sand (0.34) were used for calibration of the 2010 LiDAR data. In 2012, only the asphalt and sand were used. The WFs from the calibration targets were normalized for the AGC (Eq. 1). The normalization notably decreased the intratarget variation in peak amplitude. The normalized WFs were resampled at 0.05 ns by linear interpolation. Eq. (8) was then applied separately to each resampled WF, and the final estimate of S(t) was obtained by averaging over all pulses and targets. The noise (c) in each dataset was estimated, using the same pulses. The last 50 WF amplitude values from each pulse were selected and the mean and SD of these were calculated. We observed no variation in the mean c between different surfaces. 2.4.5. Echo detection from the simulated WFs Since there can be considerable differences in the R values produced by different echo detection algorithms (Wagner, Ulrich, Melzer, Briese, & Kraus, 2004), it was essential to correctly reproduce the constant fraction discriminator (CFD) algorithm of the ALS60. We applied the CFD to the simulated data, using the WFs directly. In real ALS60 data,

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

671

Fig. 5. Geometric-optical vegetation models. A 3-m-high birch (A), individual raspberry and fireweed shoots (B), and a 6 × 6-m area of fireweed (C) and raspberry (D) canopies (heights of 1.5 and 1.2 m). The scale differs in each subfigure.

the discrete R data are resolved by separate circuits that are probably faster than 1 GHz. We assumed, however, that the signal processed by the DR circuits is similar to that fed to the WF digitizer and that the synchronization is perfect. Since the exact values for the CFD parameters were unknown, we experimented with different values and evaluated the performance by feeding real WFs from different targets (planar surfaces, mature forest, birch saplings) to our algorithm. The R values obtained by our CFD were compared with those of the DR echoes in the real data. The absolute differences in R were below 0.10 m in over 95% of the pulses, and there was no bias. Finally, a peak amplitude-based R correction was applied as in the real Leica ALS50/60/70 data. 2.4.6. Geometric-optical vegetation models Three common tree and competing vegetation species were selected (Fig. 5). These were European white birch (Betula pendula Roth), American red raspberry (Rubus idaeus L.), and fireweed (Chamerion angustifolium (L.) Halub). Geometric-optical representations of these were created, using regression models of the shoot or branch structure, digitized leaf shapes, and literature values for leaf optical properties and leaf area index (LAI). The birches were created, using an empirical regression model (Lintunen, Sievänen, Kaitaniemi, & Perttunen, 2011) that contains information on the branching and shoot structure, as well as individual leaf area and the number of leaves per shoot. For raspberry and fireweed, we performed measurements of shoots and fitted regression models explaining the number, size, and positions of the leaves along the shoot. For deriving the leaf shape, sample leaves were digitized and scaled to match the leaf area given by the regression models. Leaf zenith angles were generated, using theoretical leaf angle distributions (LADs) (Weiss, Baret, Smith, Jonckheere, & Coppin, 2004). A spherical LAD was used for birch and planophile for raspberry and fireweed. The leaf azimuth angles were random. Directional-hemispherical reflectance and transmittance factors (0.47 and 0.50), measured for birch leaves at λ of 1064 nm (Lukeš, Stenberg, Rautiainen, Mõttus, & Vanhatalo, 2013), were used for all species as an approximation of leaf ρ and τ. The ρ values of the shoot elements in raspberry and fireweed were assumed equal to leaf ρ. For the birch branches and the trunk, a ρ of 0.50 was used. Zero τ was assumed for all shoot, branch, and trunk elements. The ground was modeled as a Lambertian planar surface with a ρ of 0.30, which corresponds to an average observed for mosses at λ of 1064 nm (Lang et al., 2002).

A 6 × 6-m area was simulated at a time. An individual birch was placed in the middle (Fig. 5A). In the case of raspberry and fireweed, the number of shoots (Fig. 5B) was adjusted to obtain a given LAI, and the area was then filled with shoots at random locations (Fig. 5C,D). We used the LAI values by Kuusk, Lang, and Nilson (2004), estimated for clear-cut areas in Sweden and Estonia (3.35 for raspberry, 2.0 for fireweed). Small random variations in LAI were allowed between each vegetation instance. 2.5. Waveform feature extraction Radiometric and geometric WF features were calculated for the first DR echo detected. The extraction was initiated by delineating a noiseexceeding segment around the echo. The radiometric features were calculated for that segment. These were the number of peaks (Npeaks), peak amplitude (Apeak), full width at half maximum (FWHM), and the total energy (E). In addition, the echo height from ground (H) was calculated, based on the ground elevation measured with GNSS. The process was similar for real and simulated data, except that the real data were normalized for the AGC (Eq. 1), because the simulated sensor did not have the AGC. 2.6. Short description of the data analyses 2.6.1. Comparison of features in real and simulated data Five LiDAR datasets (2010_12_1, 2012_05_1, 2012_10_1, 2012_20_1, 2012_27_1) were used to validate the simulation model. Data representing the three species modeled were simulated. First, real pulses that had intersected a field sample were selected. A new instance of vegetation was created for each pulse, using the field-measured h. For birch, trees with a h from 2.5 m to 5 m were used. The simulation was performed, using the same scan/pulse geometry and tree/sample locations as in the real data. The WF sampling rate in the simulations was 1 ns. The H histograms, and the mean and SD of each of the WF features were compared between the simulated and real data. When the WF features were examined, echoes with H N 1.5 m were selected in birch samples to avoid confusion with the ground vegetation. See the results in Section 3.1. 2.6.2. Simulator sensitivity to vegetation parameters To test the effects of the various assumptions made in the vegetation models, the model parameters were altered and the simulated Apeak

672

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

Fig. 6. Real and simulated echo height histograms for 2.5–5-m-high birch, 0.6–1.8-m-high raspberry, and 0.8–2.5-m-high fireweed in the 2010_12_1 dataset.

data were examined. The simulations were performed, using the S(t) of the 2010_12_1 dataset, a sampling rate of 1 ns, and the birch model without the ground surface. Nadir scanning was used. The R was 1200 m and the position of the sensor varied randomly ± 0.30 m in XY. The mean leaf inclination angle, leaf ρ and τ, and the total leaf area were the parameters investigated. Variation in total leaf area was obtained by scaling individual leaves with factors of 0.5–1.5. See Section 3.2.

each parameter configuration. When the emitted pulse width was altered, the amplitude values were scaled correspondingly to keep the pulse energy constant. The simulated first-echo data were classified, using linear discriminant analysis (LDA) and features Apeak, FWHM, and E. The training and validation datasets were the same. The classification was performed at the level of individual pulses and fivepulse samples (mean features used). These are referred to as ‘1-pulse’ and ‘5-pulse’ data, respectively. See Section 3.3.

2.6.3. Effects of sensor parameters on species classification accuracy in the simulated data The effects of footprint size, WF sampling rate, and the emitted pulse width on species classification accuracy were studied for the three species modeled. The flying altitude was constant at 1200 m, and the scan zenith angle varied randomly from 0° to 15°. The S(t) of the 2010_12_1 dataset, 1-ns sampling rate, and beam divergence of 0.15 mrad (approx. 0.18-m footprint) were used as the default. One parameter was altered at a time, and 500 pulses were simulated with

2.6.4. Analyses of real data for vegetation classification Pulses that had produced a single DR echo within a predefined radius from the center of a field sample (Section 2.2) were selected. For trees, a radius of 0.5 m was used, and echoes below 1 m were excluded. The species or species groups with less than 30 pulses were omitted. In total, 28 species or species groups fulfilling the criterion were found (Table 1). The WF features for the selected species were analyzed. A classification into four vegetation functional groups, important from a silvicultural standpoint, was carried out, using LDA.

Table 3 Comparison of echo height (H, m), number of peaks (Npeaks), peak amplitude (Apeak, digital number), echo width (FWHM, ns), and total energy (E, digital number) between the simulated and real data. Values given are the mean and (SD). Relative mean difference is denoted by Δ. Birch

Raspberry Reference

Δ, %

Simulated

Reference

Δ, %

2010_12_1 (N = 1733, 452, 316 (birch, raspberry, fireweed)) H 2.74 (0.70) 2.59 (0.66) 6 Npeaks 1.47 (0.52) 1.12 (0.33) 31 Apeak 59.3 (12.6) 72.9 (13.7) −19 FWHM 12.1 (3.4) 12.2 (3.3) −1 E 1027 (227) 1184 (121) −13

0.60 (0.19) 1 (0) 110.9 (13.6) 9.8 (1.2) 1327 (94)

0.68 (0.24) 1.00 (0.05) 108.3 (9.0) 8.9 (0.7) 1300 (104)

−11 0 2 9 2

0.53 (0.40) 1.02 (0.12) 87.8 (13.4) 10.7 (2.0) 1182 (184)

0.70 (0.31) 1 (0) 96.0 (10.7) 9.9 (1.3) 1264 (136)

−25 2 −9 8 −7

2012_05_1 (N = 480, 1301, 1173) H 2.63 (0.61) Npeaks 1.43 (0.51) Apeak 76.2 (17.6) FWHM 11.0 (3.1) E 1202 (223)

2.38 (0.59) 1.15 (0.36) 76.3 (16.6) 11.7 (3.2) 1153 (137)

10 24 0 −6 4

0.61 (0.24) 1.00 (0.03) 137.0 (23.3) 10.1 (1.3) 1665 (219)

0.73 (0.27) 1 (0) 142.2 (15.1) 9.1 (0.5) 1682 (164)

−17 0 −4 11 −1

0.38 (0.29) 1 (0) 118.0 (20.4) 10.2 (1.2) 1480 (267)

0.61 (0.31) 1.00 (0.03) 111.9 (13.7) 9.6 (1.0) 1386 (194)

−37 0 6 6 7

2012_10_1 (N = 681, 2210, 1758) H 2.58 (0.61) Npeaks 1.32 (0.47) Apeak 68.9 (14.5) FWHM 15.1 (4.1) E 1377 (248)

2.32 (0.55) 1.10 (0.30) 73.1 (13.7) 15.6 (4.0) 1438 (141)

11 20 −6 −3 −4

0.70 (0.19) 1 (0) 124.2 (14.3) 12.3 (0.8) 1797 (147)

0.64 (0.26) 1 (0) 131.9 (11.2) 11.7 (0.5) 1894 (150)

9 0 −6 5 −5

0.46 (0.18) 1 (0) 111.3 (15.5) 12.4 (0.6) 1652 (207)

0.42 (0.25) 1 (0) 106.1 (12.5) 11.9 (0.7) 1572 (209)

10 0 5 4 5

2012_20_1 (N = 312, 1104, 952) H 2.77 (0.64) Npeaks 1.44 (0.50) Apeak 64.2 (12.1) FWHM 17.2 (4.8) E 1432 (273)

2.40 (0.58) 1.15 (0.36) 70.6 (13.7) 17.3 (4.5) 1610 (181)

16 26 −9 −1 −11

0.69 (0.15) 1 (0) 126.1 (10.0) 12.5 (0.7) 1845 (106)

0.61 (0.25) 1 (0) 132.1 (10.6) 12.0 (0.5) 1954 (152)

13 0 −4 4 −6

0.49 (0.13) 1 (0) 110.8 (11.6) 12.8 (0.5) 1704 (138)

0.48 (0.25) 1 (0) 107.1 (10.8) 12.3 (0.8) 1661 (215)

2 0 3 4 3

2012_27_1 (N = 200, 836, 618) H 2.83 (0.63) Npeaks 1.54 (0.50) Apeak 62.2 (9.7) FWHM 18.0 (4.8) E 1502 (229)

2.38 (0.54) 1.14 (0.35) 72.7 (12.6) 18.5 (4.8) 1752 (165)

19 35 −14 −3 −14

0.71 (0.14) 1 (0) 129.9 (8.7) 12.7 (0.7) 1937 (89)

0.70 (0.24) 1 (0) 132.9 (10.9) 12.2 (0.5) 2010 (162)

2 0 −2 4 −4

0.50 (0.10) 1 (0) 114.8 (10.7) 13.0 (0.5) 1793 (114)

0.54 (0.24) 1 (0) 109.0 (10.3) 12.5 (0.8) 1717 (215)

−7 0 5 4 4

Reference

Δ, %

Fireweed

Simulated

Simulated

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

673

Fig. 7. Effect of footprint size on the simulated (top) and real (bottom) WFs from birch crowns in the 2012 LiDAR data. The first 100 pulses that produced a DR echo more than 1.5 m above the ground were chosen. Footprint diameter corresponds to the 1/e definition. The x axis shows the time (ns), and the y axis the amplitude (DN).

The classification was performed at the level of individual pulses (‘1-pulse’ data) and vegetation samples, using mean features per sample (‘mean-of-sample’ data). See Section 3.4. 3. Results 3.1. Comparison of features in real and simulated data The H histograms of the simulated echoes were generally in good agreement with the real data (Fig. 6, Table 3), especially for birch, although there was an overestimation of the mean H by 0.15–0.45 m. The overestimation was also observed in vertical energy distributions calculated from the WF data directly. Since these were not affected by the echo detection method applied, the H discrepancies were most likely caused by differences in the real vs. modeled vegetation. For the other species, the mean H was simulated with reasonable accuracy,

but the SD was underestimated towards the highest flying altitudes (Table 3). In the two high-altitude datasets, 2012_20_1 and 2012_27_1, this was also clear in the H histograms (not shown). The larger SD in the real data could have been due to the spatial distribution of the individual shoots, which in reality probably is more clumped at the scale of the largest footprints. The characteristics of the WF features were well reproduced in the simulations (Table 3). The relative interspecies order of the mean Apeak, FWHM, and E was preserved, except for E in 2012_27_1. The mean FWHM tended to be underestimated for birch and overestimated for the other species. However, the differences between the simulated and real data were within 11%. The SD of FWHM was similar to that in the real data. The mean Apeak and E were also reproduced with reasonable accuracy. The largest differences were for birch in datasets 2010_12_1 and 2012_27_1, in which an underestimation of 13–19% was observed. The values for raspberry were similar to those in the

Fig. 8. A–C: Effect of vegetation model parameters on the simulated peak amplitude from the crown of a 3-m-high birch. D–F: Dependency between the peak amplitude and the visually estimated canopy cover (~leaf area) in birch (D), raspberry (E), and fireweed (F) samples in the 2012_10_1 dataset.

674

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

Fig. 9. Effect of footprint size, emitted pulse width, and the WF sampling rate on the overall classification accuracy between simulated birch, raspberry, and fireweed vegetation. The x axis is linear and the range from zero to the maximal values is shown.

real data, whereas for fireweed a slight underestimation was observed in the 2010 data and overestimation in the 2012 data. The SD of E was always overestimated for birch. Further examination revealed that there were more multireturn pulses in the simulated than in the real data (Fig. 7). The birch model had a vertical gap between the crown base and the ground, which in reality is usually filled with grasses and small seedlings. Different near-ground behavior was also observed in the H histograms (Fig. 6). Probably for the same reason, the number of peaks within an echo (Npeaks) was also overestimated in the simulated birch data. For the other two species, the SD of Apeak and E tended to be overestimated in the low-altitude datasets and underestimated towards the higher altitudes. This could have been due to the spatial arrangement of the shoots in the geometric-optical models, being probably too heterogeneous at the scale of the smallest footprints, but too homogeneous at the larger scale. The findings considering the echo H histograms also support this. The performance of the simulator in reproducing signals acquired with different footprint sizes was important in considering RA2. This could be evaluated in the 2012 data, in which the footprint size on the ground varied. The variation in Apeak decreased and the mean FWHM increased with increasing footprint size in the real data (Fig. 7). The decrease in Apeak variation was most probably caused by an averaging effect, and the increase in FWHM by the greater depth of volumetric scatterers in the large-footprint data. Similar behavior was also seen in the simulated data. The emitted pulse width was slightly smaller in the 2012_05_1 dataset (Table 2), which affected the echo width as well.

3.2. Simulator sensitivity to vegetation parameters The leaf inclination angle of the simulated birch strongly impacted the Apeak (Fig. 8A). The mean leaf angle in the spherical LAD is 57.3° (Weiss et al., 2004). For that angle, the simulated Apeak in Fig. 8A was approximately at a DN of 55. Fig. 8A shows that even a 10° change from this angle causes the simulated Apeak to vary 20–30%. The leaf ρ did not have as profound an effect (Fig. 8B). When ρ was altered from 0.4 to 0.5, the average Apeak increased from approx. 55 to 70. The response to τ was also linear, but smaller (not illustrated), because τ affects only the second and subsequent scattering orders. When the individual leaf area was scaled from 0.5 to 1.5, Apeak increased from approx. 50 to 80 (Fig. 8C). The result can be compared with the real data in which the canopy cover percentage (~leaf area) was estimated in the field. The response of Apeak to change in leaf area is similar in the simulated (Fig. 8C) and real data (Fig. 8D–F). Due to the subjective nature of the field canopy cover estimates, the dependencies observed are only tentative. However, this finding provides evidence that the simulator is functioning logically.

Fig. 10. Mean WF features per species or species group in the LiDAR datasets. To facilitate comparisons, the values are normalized to those obtained from sand road samples (value = 1) in each dataset.

3.3. Effects of sensor parameters on the species classification accuracy in simulated data Fig. 9 shows the effect of the three parameters investigated on the overall classification accuracy among the species modeled. In general, and as expected, higher performance was obtained with the less noisy 5-pulse data than the 1-pulse data. The classification accuracy even reached 100% with the largest simulated footprints. The three species differ quite clearly in structure, which partly explains the good classification performance. It is also possible that the diversity in structure and reflectance properties of the real data was not fully captured by the models, reducing the intraclass variation. However, the classification accuracies were also high in the real data, being 78.6–88.3% and 89.6–95.8% for the 1-pulse and the 5-pulse mean features. The simulated beam divergence was 0.05–0.7 mrad, corresponding to footprint diameters of 0.06–0.84 m at 1200 m. Low classification

Table 4 Overall classification accuracy (%) and kappa in LDA classification of four vegetation functional groups: 1) conifers, 2) broad-leaved trees, 3) low vegetation (green), and 4) low vegetation (barren) + abiotic material. Features used were the peak amplitude (Apeak) and total echo energy (E). Dataset

2010_12_1 2010_19_1 2010_19_2 2010_30_1 2010_30_2 2012_05_1 2012_10_1 2012_20_1 2012_27_1

1-Pulse

Mean-of-sample

Correct-%

Kappa

Correct-%

Kappa

Pulses per sample

68.1 75.7 71.3 75.3 72.5 76.1 83.1 85.4 86.7

0.51 0.62 0.54 0.57 0.50 0.66 0.76 0.79 0.80

74.9 77.2 77.0 77.1 73.8 75.6 83.3 84.2 84.8

0.59 0.62 0.62 0.60 0.52 0.67 0.77 0.78 0.79

3.1 2.2 2.5 1.4 1.5 8.7 9.7 5.9 4.6

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

accuracies were associated with small footprints, and the accuracy saturated after approx. 0.3 mrad (0.36 m) (Fig. 9). The emitted pulse width in the simulations was 1.4–14 ns (FWHM). The response to the varying pulse width was less apparent than the footprint size (Fig. 9). Further examination revealed that while an increase in pulse width reduces the intraspecies variation in the WF features, it also diminishes the interspecies differences in the mean features. This was particularly apparent for Apeak and E. Thus, the separability did not improve. The simulated WF sampling rate was 0.25–4 ns. The sampling rate did not considerably affect the classification accuracy, although a small decrease was observed with 3–4-ns sampling, using the 5-pulse data (Fig. 9). 3.4. Analyses of real data for vegetation classification Fig. 10 shows the mean species-specific WF features. The interspecies differences were consistent across the datasets. The common tree species could not be well-separated, except for pine and juniper, which gave rise to lower Apeak and E values than the other species. Interspecies variation in FWHM was small. Spruce did not differ from the broad-leaved trees in any of the features, which is a drawback, because information on broad-leaved trees in conifer plantations is important. Two types of low vegetation could be distinguished. These were ‘green’ (fireweed, raspberry, brackenfern, grasses), and ‘barren and low’ vegetation (shrubs, mosses, ripening grass) (Table 1, Fig. 10). The former are characterized by high Apeak and E and moderate FWHM, the latter by somewhat lower Apeak and E and a narrow FWHM. Abiotic material resembles the latter group. What should be noted is that for tree crowns, Apeak is typically dampened, due to volumetric scattering. Therefore, Apeak values from the ground are equal to or higher than those from trees. Volumetric scattering did not affect E highly, and therefore E is a better indicator of reflectance. The Npeaks was 1 for the low vegetation and only slightly higher for trees (not shown). No differences between tree species were observed. The species were divided into four vegetation functional groups: 1) conifers, 2) broad-leaved trees, 3) low vegetation (green), and 4) low vegetation (barren) + abiotic material. First, a classification was performed for the 1-pulse data, using one WF feature at a time and the different combinations. Using E and Apeak together gave the best performance in almost all datasets (overall accuracy 68.1–86.7%, Table 4). Adding FWHM improved the results marginally (b 1 percentage point) in some of the datasets. The ratio of E to Apeak was highly correlated with FWHM, and thus the FWHM furnished no additional information. Secondly, the classification was repeated for the mean-of-sample data, using E and Apeak. Accuracies were improved in the 2010, but not in the 2012 data. Because of the larger sample size for low vegetation in 2012, the low vegetation classes were overrepresented in the 1-pulse data of 2012. These classes were also the easiest to classify, and thus the overall accuracy did not improve by using the mean-of-sample data, although the accuracies (correct-%) for individual groups were improved. Adding the maximum H within a sample as an explanatory variable increased the overall accuracy by 0.4–7.8 percentage points. We expected better separation between green low vegetation and trees. However, the LiDAR underestimated the h of small trees considerably. Thus, the echoes from tree crowns were still confused with those from low vegetation. 4. Discussion 4.1. Calibration and validation of the simulation model with field reference data (RA1) We presented a model for simulating WF LiDAR data. MCRT and detailed geometric-optical vegetation models were used to estimate the target differential scattering cross-section, which was then convolved with the LiDAR system WF to obtain the signal detected.

675

We used targets of known reflectance to calibrate the model output. Availability of very accurately geo-referenced and up-to-date field samples enabled comparisons of the simulations with the real observations. Qualitative comparisons have been reported for spaceborne large-footprint LiDAR (North et al., 2010), but to our knowledge, quantitative comparisons with small-footprint airborne LiDAR have not been reported before. Absolute validation of MCRT models is difficult, due to the compound effects of uncertainties in the radiative transfer model and in the vegetation models used (Hancock et al., 2011). In addition, the sensor model must be embedded in the chain. However, our results showed that the radiometric and geometric features in the simulated WF data were in good agreement with the field observations, in spite of the uncertainties and assumptions that are discussed below. The results can serve as a basis for attempts to calibrate and validate simulation models so that they could be used, e.g. in the planning of LiDAR surveys (Disney et al., 2010). The 3D structure of the vegetation was based on empirical regression models. The countrywide modeling dataset used for birch was representative of small trees (Lintunen et al., 2011). Thus, the birch model was likely the most accurate available. The raspberry and fireweed models were fitted to our own measurements of the shoot structure, and the amount of shoots was adjusted to reasonable LAI values found in the literature. The behavior of both the geometric and radiometric WF features with different footprints suggests that the spatial structure of the shoots in the real vegetation was not fully captured by our models of fireweed and raspberry. Furthermore, the near-ground vegetation affects the total reflected radiation, and the vegetation can vary widely from site to site. For example, the fireweed samples in 2012 were mainly collected from a newly established burned stand, and there was an absence of low grass under the fireweed canopy. As a consequence, our model produced small overestimations of the signal strength for fireweed in the 2012 data. There were small discrepancies in the real vs. simulated H histograms. Our echo detection algorithm was unbiased, and thus the ranging process was simulated correctly. However, there could have been small coordinate shifts introduced in postprocessing in the real data. In addition, the typical 0.20–0.25-m XY errors in the real data could have affected the results for birch, since the vertical distance from the tree trunk affects the H values (Korpela, Hovi, & Morsdorf, 2012). There were also small geometric differences in the modeled vs. real vegetation, such as the lack of grasses on the ground, and possibly small-scale structural differences (gap-size distribution). Furthermore, errors could have been introduced in the field measurements, e.g. by the antenna penetration into the ground by a few centimeters, or by growth of the vegetation during the time between the LiDAR campaign and the field measurement (max. 3 weeks). The sensitivity analyses showed that the LAD very strongly impacted the simulated signals. We used a spherical LAD for birch and planophile for the other species, based on visual examination of the plants in situ and in photographs. A spherical LAD is often used for tree canopies, but recent studies by Pisek, Ryu, and Alikas (2011) suggested rather a planophile LAD for boreal trees. However, their measurements were performed on large trees. We also tested the planophile LAD with birch, but the spherical case agreed more favorably with the real data. Leaf ρ and τ for birch were obtained from measurements in the study area (Lukeš et al., 2013), but data for raspberry and fireweed were not available. Thus, we used the values of birch for all species. The sensitivity analysis (Section 3.2) suggested that ρ and τ did not affect the results as much as the LAD does. The variation in ρ or τ between fresh plant leaves at 1064 nm is typically between 0.4 and 0.55 (Hosgood et al., 2005). However, this uncertainty constitutes an additional error source. We assumed a bi-Lambertian scattering behavior for the leaves, which is a strong assumption. Leaves often have a specular ρ component, and its magnitude can vary widely, depending on the leaf surface properties (Bousquet, Lacherade, Jacquemoud, & Moya, 2005). In addition, there are other minor error

676

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

sources, such as the use of a constant 5-m area viewed by the receiver and the simulation of only 10 scattering orders. However, our tests showed that the bias by these factors is only a few percent. For simplicity, we omitted the modeling of the wave phenomena. For example, the coherent backscattering can have an effect (Hapke, DiMucci, Nelson, & Smythe, 1996), increasing the intensity of the firstorder scattering in the hotspot view illumination geometry of the LiDAR. We used reference tarps and other artificial surfaces for the radiometric calibration (Section 2.4.4), employing their HDRF values observed at nadir. Clearly the measurement configuration did not correspond to the LiDAR geometry. A rise in reflected intensity near the hotspot geometry was observed for the tarps under laboratory conditions (Kaasalainen, Ahokas, Hyyppä, & Suomalainen, 2005), which certainly affected our results. The absolute values of the target backscattering reflectance factors remained uncertain in our data, but we noted that the WF peak amplitude values were linearly correlated with the HDRF values measured. This suggests that the relative backscattering effects did not vary among the targets used. We used a simplified sensor model, in which a large linear dynamic range of the photon detector and subsequent circuits was assumed. Comparison of the amplitude values with the HDRF data measured suggested no significant deviation from linearity in the range of HDRF values tested. The DR circuit in the ALS60 sensor applies some type of threshold mechanism, so as not to record unwanted DR and WF data. The WF recording is triggered by the first DR echo. The real signals may thus be affected by transmission losses preceding the first echo, e.g. if a branch is intersected by the pulse without triggering an echo. A dataset-specific noise term with random (Gaussian) variation was added to the simulated signal. Thus, the mean noise in the simulations did not vary between targets of different ρ. In reality, some autocorrelation over the vegetation patterns can be expected, due to solar illumination and the AGC, although we did not observe a substantial autocorrelation of noise. In all, we suggest further studies on the calibration and validation of the simulation model. The validation and testing of the various assumptions of the MCRT model can be performed by comparing simulated shoot/branch scattering phase functions (Smolander & Stenberg, 2003) with measurements made in the laboratory (e.g. Mõttus et al., 2012), where the positions and orientations of individual scattering elements can be determined more accurately. The same calibration and validation procedure as used here could also be repeated for a larger variety of species. The uncertainty related to the model parameters could be reduced by more accurate measurements of the vegetation (LAI, LAD, etc.), as well as the scattering properties of the calibration targets. Fixing some of the parameters would allow for an iterative solution of the unknowns, e.g. the individual leaf BRDF. On the other hand, considering the practical applicability, the vegetation models should be readily available. Modeling environments already exist that could provide the input for MCRT models (e.g. Perttunen, Sievänen, & Nikinmaa, 1998). TLS could be utilized as well (Côté et al., 2009; Pfeifer, Gorte, & Winterhalder, 2004; Raumonen et al., 2013). Independent of the method, the vegetation models should be accurate for the study area (vegetation population) in question, or alternatively, the scattering model should be robust to the uncertainty in the vegetation parameters. Improvement in the accuracy of the available vegetation models can be expected, since interest in computationally intensive modeling methods is increasing and measurement data of the required input parameters accumulate. 4.2. Effects of sensor parameters on species classification accuracy (RA2) We used the simulator to examine the effect of footprint size, emitted pulse width, and WF sampling rate on the species classification accuracy, employing the three species modeled. We showed that the accuracy increases with increasing footprint size and that the gain in accuracy was highest with small footprints (b 0.3–0.36 m). A similar

effect was observed in the real data. The results suggest that very small footprint is not optimal in the vegetation examined (Table 1). Backscattering of a very small footprint correlates more with the reflectance than with the structure, which is what separates the species. We see potential in the usage of LiDAR data that combines different footprint sizes. Such data potentially characterize the gap-size distribution of the canopy and could be used for enhanced classification. We propose research in this topic, not only in juvenile forests but also for taller vegetation. The emitted pulse width and WF sampling rate had minor effects on the classification performance. It appears that a 2-ns sampling rate suffices and this would reduce the storage and data transfer needs. We used only rather simple WF features. The sampling rate can induce a different effect if more complex feature extraction procedures are implemented, e.g. by Gaussian decomposition (Wagner et al., 2006). The same applies to the effect of the emitted pulse width, which was negligible in our data, but can be different in large trees or tall vegetation, because shorter pulses enhance the R resolution. 4.3. Airborne WF LiDAR in the mapping of juvenile forest vegetation (RA3) We investigated radiometric WF features for the separation of species and species groups. Some interspecies differences were found, and they were consistent across the LiDAR datasets. The LDA classification accuracies were 68.1–86.7% (kappa 0.50–0.80) for the four vegetation functional groups when 1-pulse data were classified and were similar when average features per vegetation sample were used. Korpela et al. (2008) reported classification accuracies of 61–79% (kappa 0.44–0.67) for four operational classes (conifer, broad-leaved, low vegetation, abiotic material), using high-resolution multispectral aerial images and features extracted from DR LiDAR data. The favorable accuracy obtained here could have been due to the capability of the WF data in providing measures of both vegetation structure at the footprint scale (echo width or amplitude to total energy ratio) and the reflectance in the near-infrared (total energy). Due to volumetric scattering, the intensity data in the DR sensors do not necessarily represent either of these properties. Classification results with real data were in line with the simulations. The lowest performance was obtained with the small-footprint datasets 2010_12_1 (footprint 0.18 m) and 2012_05_1 (0.08 m), although the footprint effects were not as distinct as in the simulations. This could have been due to small structural differences in the vegetation of the operational classes, e.g. between the real coniferous and broad-leaved trees, whereas the simulated species differed more in structure. The results can be different with large trees in which there are more pulses per object and different distribution metrics of the LiDAR features are available (Orka, Naesset, & Bollandsas, 2009). The differences between the 1-ns and 2-ns WF sampling rates were small, similarly as with the simulated data. Our focus was not so much in practical inventory procedures that aim at deducing the optimal silvicultural regime (e.g. Korhonen et al., 2013). We anticipate that the estimation of the treatment need, using WF data, can be direct, i.e. in situ training data are collected for the treatment need. The LiDAR features, auxiliary data on past silvicultural regimes, site conditions, and other expert knowledge are then combined in an imputation scheme (Packalén & Maltamo, 2007). Alternatively, it may be possible to predict the stand variables first, e.g. the density and height of the vegetation components, and deduce the treatment need from these data. A prerequisite for this approach is that the stand variables must be reliably measurable in LiDAR data and correlate reasonably with the treatment need variables. In addition, the resulting longer modeling chain may be prone to error propagation. One option is to incorporate the spatial dependencies in the estimation, e.g. using multidimensional rasters of LiDAR features or predicted forest variables to determine the local treatment need (e.g. Korpela et al., 2008; Wulder & Boots, 1998). In addition to the radiometric features,

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

height-based features could potentially be used to discriminate the fastand slow-growing species (Korpela et al., 2008; Naesset & Bjerknes, 2001). Furthermore, improvement can possibly be obtained by combining the usage of WF data in characterizing small-scale vegetation structural differences with spectral information from passive optical sensors. Passive sensors have already been tested in combination with DR LiDAR (Korhonen et al., 2013; Korpela et al., 2008). It is evident that WF LiDAR data extend the list of potential classification features. In juvenile forest stands, accurate elevation models are crucial for the height estimation of vegetation, and a more enhanced digital elevation model (DEM) estimation is possible, using the echo width features (Wagner, Hollaus, Briese, & Ducic, 2008). A positive correlation between the H and the FWHM of the echo was observed in our data as well. Since we restricted our analyses to consider the basic pulse-vegetation phenomena, we propose further studies in both DEM estimation (algorithm development) and in the practical use of the WF data in seedling and sapling stand inventories (design of interpretation techniques at the plot or stand levels). 5. Conclusions Our practical objective was to explore the potential for smallfootprint airborne WF-recording LiDAR in characterizing juvenile boreal forest vegetation. We focused on the species classification task, which is important from the management standpoint. We used simulations along with real datasets, and thus the documentation and validation of our simulator was an important part of the study. The validation and sensitivity analysis of the MCRT-based simulator showed that the simulated WFs were in good correspondence with the real data. This supports the status of MCRT as one of the most accurate methods for simulating RS signals. WF-recording LiDAR is particularly suitable for linking simulations with real observations, because the user has control over the photons transmitted. Although limited to only three species, the findings reported here are to our knowledge among the first to report direct comparisons between real smallfootprint WF data and those obtained with MCRT simulation. The results serve the need in the research community to calibrate and validate simulation models to enhance their operational applicability. We used the simulator to determine the effects of sensor parameters on the classification performance in the three species modeled. The results showed that footprint diameters under 0.3–0.36 m should be avoided, whereas the choice of the pulse width or sampling rate is not crucial to the classification performance in the type of vegetation and WF features used here. Combinations of various footprints could be useful for enhanced description of vegetation structure, and we suggest further research in this topic. Finally, we used real WF data to classify four vegetation functional groups, important from the silvicultural standpoint. We observed small improvements over existing studies performed with DR LiDAR and passive optical sensors. This was most likely attributed to the capability of WF LiDAR to characterize both vegetation structure at the footprint scale and the reflective properties. To develop operationally applicable inventory methods, we suggest research on 1) DEM algorithms that improve the DEM quality, which is important in short vegetation, and 2) interpretation methods that take into account the spatial patterns of the LiDAR features. Acknowledgements The study was funded by TEKES, University of Helsinki, and by the Academy of Finland (Timestamped and free photons — quantitative airborne remote sensing of forests). Arthur Rohrbach, Yuji Kuwano, and Ron Roth from Leica Geosystems and Juha Hyyppä from FGI were vital in bringing the sensor to Hyytiälä in 2010. We had access to FGI's reflectance tarps. People at Finnmap arranged the flights and data postprocessing. Paula Litkey familiarized us with the WF analysis, and

677

the discussions with the LAI Detectives group were helpful for an understanding of radiative transfer modeling. Anna Lintunen helped with the implementation of the 3D birch model. Lauri Korhonen provided early access to the manuscript of Korhonen et al. (2013) that helped in understanding the operational aspects of seedling/sapling stand inventory. Jari Vauhkonen gave valuable comments on the text during the writing process. References Baltsavias, E. (1999). Airborne laser scanning: Basic relations and formulas. ISPRS Journal of Photogrammetry and Remote Sensing, 54, 199–214. Bittner, S., Gayler, S., Biernath, C., Winkler, J. B., Seifert, S., Pretzsch, H., et al. (2012). Evaluation of a ray-tracing canopy light model based on terrestrial laser scans. Canadian Journal of Remote Sensing, 38, 619–628. Blair, J. B., Coyle, D. B., Bufton, J. L., & Harding, D. J. (1994). Optimization of an airborne laser altimeter for remote sensing of vegetation and tree canopies. In T. I. Stein (Ed.), Proceedings of IGARSS'94, Pasadena, CA, 8–12 Aug (pp. 939–941). Bousquet, L., Lacherade, S., Jacquemoud, S., & Moya, I. (2005). Leaf BRDF measurements and model for specular and diffuse components differentiation. Remote Sensing of Environment, 98, 201–211. Chauve, A., Vega, C., Durrieu, S., Bretar, F., Allouis, T., Deseilligny, M. P., et al. (2009). Advanced full-waveform lidar data echo detection: Assessing quality of derived terrain and tree height models in an alpine coniferous forest. International Journal of Remote Sensing, 30, 5211–5228. Côté, J. -F., Widlowski, J. -L., Fournier, R. A., & Verstraete, M. M. (2009). The structural and radiative consistency of three-dimensional tree reconstructions from terrestrial lidar. Remote Sensing of Environment, 113, 1067–1081. Disney, M. I., Kalogirou, V., Lewis, P., Prieto-Blanco, A., Hancock, S., & Pfeifer, M. (2010). Simulating the impact of discrete-return lidar system and survey characteristics over young conifer and broadleaf forests. Remote Sensing of Environment, 114, 1546–1560. Disney, M. I., Lewis, P., Gomez-Dans, J., Roy, D., Wooster, M. J., & Lajas, D. (2011). 3D radiative transfer modelling of fire impacts on a two-layer savanna system. Remote Sensing of Environment, 115, 1866–1881. Disney, M. I., Lewis, P., & North, P. R. J. (2000). Monte Carlo ray tracing in optical canopy reflectance modelling. Remote Sensing Reviews, 18, 163–196. Disney, M. I., Lewis, P., & Saich, P. (2006). 3D modelling of forest canopy structure for remote sensing simulations in the optical and microwave domains. Remote Sensing of Environment, 100, 114–132. Ene, L. T., Naesset, E., Gobakken, T., Gregoire, T. G., Ståhl, G., & Nelson, R. (2012). Assessing the accuracy of regional LiDAR-based biomass estimation using a simulation approach. Remote Sensing of Environment, 123, 579–592. Hancock, S., Disney, M., Muller, J., Lewis, P., & Foster, M. (2011). A threshold insensitive method for locating the forest canopy top with waveform lidar. Remote Sensing of Environment, 115, 3286–3297. Hapke, B., DiMucci, D., Nelson, R., & Smythe, W. (1996). The cause of the hot spot in vegetation canopies and soils: Shadow-hiding versus coherent backscatter. Remote Sensing of Environment, 58, 63–68. Heinzel, J., & Koch, B. (2011). Exploring full-waveform LiDAR parameters for tree species classification. International Journal of Applied Earth Observation and Geoinformation, 13, 152–160. Holmgren, J., Nilsson, M., & Olsson, H. (2003). Simulating the effects of lidar scanning angle for estimation of mean tree height and canopy closure. Canadian Journal of Remote Sensing, 29, 623–632. Hosgood, B., Jacquemoud, S., Andreoli, G., Verdebout, J., Pedrini, A., & Schmuck, G. (2005). Leaf Optical Properties EXperiment 93 (LOPEX93). Report EUR 16095 EN (revised 2005). Ispra, Italy: European Commission, Joint Research Centre, Institute for Remote Sensing Applications. Hyyppä, J., & Inkinen, M. (1999). Detecting and estimating attributes for single trees using laser scanner. Photogrammetric Journal of Finland, 16, 27–42. Jutzi, B., & Stilla, U. (2006). Range determination with waveform recording laser systems using a Wiener Filter. ISPRS Journal of Photogrammetry & Remote Sensing, 61, 95–107. Kaasalainen, S., Ahokas, E., Hyyppä, J., & Suomalainen, J. (2005). Study of surface brightness from backscattered laser intensity: Calibration of laser data. IEEE Geoscience and Remote Sensing Letters, 2, 255–259. Kajiya, J. T. (1986). The rendering equation. Siggraph, 1986, 143–150. Knyazikhin, Y., Martonchik, J. V., Myneni, R. B., Diner, D. J., & Running, S. W. (1998). Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from MODIS and MISR data. Journal of Geophysical Research, 103, 32257–32275. Korhonen, L., Pippuri, I., Packalén, P., Heikkinen, V., Maltamo, M., & Heikkilä, J. (2013). Detection of the need for seedling stand tending using high-resolution remote sensing data. Silva Fennica, 47 (article id 952, 20 p.). Korpela, I. (2006). Geometrically accurate time series of archived aerial images and airborne lidar data in a forest environment. Silva Fennica, 40, 109–126. Korpela, I. (2008). Mapping of understory lichens with airborne discrete-return LiDAR data. Remote Sensing of Environment, 112, 3891–3897. Korpela, I., Hovi, A., & Morsdorf, F. (2012). Understory trees in airborne LiDAR data — Selective mapping due to transmission losses and echo-triggering mechanisms. Remote Sensing of Environment, 119, 92–104. Korpela, I., Orka, H. O., Hyyppä, J., Heikkinen, V., & Tokola, T. (2010). Range and AGC normalization in airborne discrete-return LiDAR intensity data for forest canopies. ISPRS Journal of Photogrammetry and Remote Sensing, 65, 369–379.

678

A. Hovi, I. Korpela / Remote Sensing of Environment 140 (2014) 665–678

Korpela, I., Tuomola, T., Tokola, T., & Dahlin, B. (2008). Appraisal of seedling stand vegetation with airborne imagery and discrete-return LiDAR — An exploratory analysis. Silva Fennica, 42, 753–772. Kuusk, A., Lang, M., & Nilson, T. (2004). Simulation of the reflectance of ground vegetation in sub-boreal forests. Agricultural and Forest Meteorology, 126, 33–46. Lang, M., Kuusk, A., Nilson, T., Lükk, T., Pehk, M., & Alm, G. (2002). Reflectance spectra of ground vegetation in sub-boreal forests. Web page. Available online. http://www. aai.ee/bgf/ger2600/ (from Tartu Observatory, Estonia. Accessed 6 Feb, 2013) Lintunen, A., Sievänen, R., Kaitaniemi, P., & Perttunen, J. (2011). Models of 3D crown structure for Scots pine (Pinus sylvestris) and silver birch (Betula pendula) grown in mixed forest. Canadian Journal of Forest Research, 41, 1779–1794. Lukeš, P., Stenberg, P., Rautiainen, M., Mõttus, M., & Vanhatalo, K. M. (2013). Optical properties of leaves and needles for boreal tree species in Europe. Remote Sensing Letters, 4, 667–676. Maltamo, M., Packalén, P., Kallio, E., Kangas, J., Uuttera, J., & Heikkilä, J. (2011). Airborne laser scanning based stand level management inventory in Finland. Proc. SilviLaser 2011, 16–19 Oct, Hobart, Australia (10 p.). Morsdorf, F., Nichol, C., Malthus, T., & Woodhouse, I. (2009). Assessing forest structural and physiological information content of multi-spectral LiDAR waveforms by radiative transfer modelling. Remote Sensing of Environment, 113, 2152–2163. Mõttus, M., Rautiainen, M., & Schaepman, M. E. (2012). Shoot scattering phase function for Scots pine and its effect on canopy reflectance. Agricultural and Forest Meteorology, 154, 67–74. Naesset, E. (2002). Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sensing of Environment, 80, 88–99. Naesset, E. (2004). Accuracy of forest inventory using airborne laser scanning: evaluating the first Nordic full-scale operational project. Scandinavian Journal of Forest Research, 19, 554–557. Naesset, E. (2007). Airborne laser scanning as a method in operational forest inventory: Status of accuracy assessments accomplished in Scandinavia. Scandinavian Journal of Forest Research, 22, 433–442. Naesset, E., & Bjerknes, K. (2001). Estimating tree heights and number of stems in young forest stands using airborne laser scanner data. Remote Sensing of Environment, 78, 328–340. Naesset, E., Gobakken, T., Holmgren, J., Hyyppä, H., Hyyppä, J., Maltamo, M., et al. (2004). Laser scanning of forest resources: the Nordic experience. Scandinavian Journal of Forest Research, 19, 482–499. North, P. R. J., Rosette, J. A.B., Suarez, J. C., & Los, S. O. (2010). A Monte Carlo radiative transfer model of satellite waveform LiDAR. International Journal of Remote Sensing, 31, 1343–1358. Orka, H. O., Naesset, E., & Bollandsas, O. M. (2009). Classifying species of individual trees by intensity and structure features derived from airborne laser scanner data. Remote Sensing of Environment, 113, 1163–1174. Packalén, P., & Maltamo, M. (2007). The k-MSN method for the prediction of speciesspecific stand attributes using airborne laser scanning and aerial photographs. Remote Sensing of Environment, 109, 328–341. Packalén, P., Suvanto, A., & Maltamo, M. (2009). A two stage method to estimate species-specific growing stock. Photogrammetric Engineering and Remote Sensing, 75, 1451–1460. Persson, Å., Holmgren, J., & Soderman, U. (2002). Detecting and measuring individual trees using an airborne laser scanner. Photogrammetric Engineering and Remote Sensing, 68, 925–932. Perttunen, J., Sievänen, R., & Nikinmaa, E. (1998). LIGNUM: A model combining the structure and the functioning of trees. Ecological Modelling, 108, 189–198.

Pfeifer, N., Gorte, D., & Winterhalder, D. (2004). Automatic reconstruction of single trees from terrestrial laser scanner data. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXV-B5, 114–119. Pisek, J., Ryu, Y., & Alikas, K. (2011). Estimating leaf inclination and G-function from leveled digital camera photography in broadleaf canopies. Trees — Structure and Function, 25, 919–924. Raumonen, P., Kaasalainen, M., Åkerblom, M., Kaasalainen, S., Kaartinen, H., Vastaranta, M., et al. (2013). Fast automatic precision tree models from terrestrial laser scanner data. Remote Sensing, 5, 491–520. Reitberger, J., Krzystek, P., & Stilla, U. (2008). Analysis of full waveform LIDAR data for the classification of deciduous and coniferous trees. International Journal of Remote Sensing, 29, 1407–1431. Reitberger, J., Schnoerr, C., Krzystek, P., & Stilla, U. (2009). 3D segmentation of single trees exploiting full waveform LIDAR data. ISPRS Journal of Photogrammetry and Remote Sensing, 64, 561–574. Roncat, A., Bergauer, G., & Pfeifer, N. (2011). B-spline deconvolution for differential target cross-section determination in full-waveform laser scanning data. ISPRS Journal of Photogrammetry and Remote Sensing, 66, 418–428. Schutz, B. E., Zwally, H. J., Shuman, C. A., Hancock, D., & DiMarzio, J. P. (2005). Overview of the ICESat Mission. Geophysical Research Letters, 32, L21S01. Smolander, S., & Stenberg, P. (2003). A method to account for shoot scale clumping in coniferous canopy reflectance models. Remote Sensing of Environment, 88, 363–373. Ullrich, A., & Pfenningbauer, M. (2011). Echo digitization and waveform analysis in airborne and terrestrial laser scanning. In Dieter Fritsch (Ed.), Photogrammetric Week 2011, Stuttgart, Germany, 5–9 September (pp. 217–228). Vastaranta, M., Holopainen, M., Xiaowei, Y., Hyyppä, J., Hyyppä, H., & Viitala, R. (2011). Predicting stand thinning maturity from airborne laser scanning data. Scandinavian Journal of Forest Research, 26, 187–196. Vierling, L. A., Xu, Y., Eitel, J. U. H., & Oldow, J. S. (2012). Shrub characterization using terrestrial laser scanning and implications for airborne LiDAR assessment. Canadian Journal of Remote Sensing, 38, 709–722. Wagner, W. (2010). Radiometric calibration of small-footprint full-waveform airborne laser scanner measurements: Basic physical concepts. ISPRS Journal of Photogrammetry and Remote Sensing, 65, 505–513. Wagner, W., Hollaus, M., Briese, C., & Ducic, V. (2008). 3D vegetation mapping using small-footprint full-waveform airborne laser scanners. International Journal of Remote Sensing, 29, 1433–1452. Wagner, W., Ulrich, A., Ducic, V., Melzer, T., & Studnicka, N. (2006). Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. ISPRS Journal of Photogrammetry & Remote Sensing, 60, 100–112. Wagner, W., Ulrich, A., Melzer, T., Briese, C., & Kraus, K. (2004). Form single-pulse to full-waveform airborne laser scanners: potential and practical challenges. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXV-B3, 201–206. Weiss, M., Baret, F., Smith, G. J., Jonckheere, I., & Coppin, P. (2004). Review of methods for in situ leaf area index (LAI) determination Part II. Estimation of LAI, errors and sampling. Agricultural and Forest Meteorology, 121, 37–53. Widlowski, J. -L., Robustelli, M., Disney, M., Gastellu-Etchegorry, J. -P., Lavergne, T., Lewis, P., et al. (2008). The RAMI on-line model checker (ROMC): A web-based benchmarking facility for canopy reflectance models. Remote Sensing of Environment, 112, 1144–1150. Wulder, M., & Boots, B. (1998). Local spatial autocorrelation characteristics of remotely sensed imagery assessed with the Getis statistic. International Journal of Remote Sensing, 19, 2223–2231.