LIDAR strip adjustment: Application to volcanic areas

Geomorphology 111 (2009) 123–135

Contents lists available at ScienceDirect

Geomorphology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / g e o m o r p h

LIDAR strip adjustment: Application to volcanic areas Massimiliano Favalli ⁎, Alessandro Fornaciai, Maria Teresa Pareschi Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via della Faggiola 32, 56126 Pisa, Italy

a r t i c l e

i n f o

Article history: Received 22 July 2008 Received in revised form 7 April 2009 Accepted 8 April 2009 Available online 16 April 2009 Keywords: LIDAR Calibration DEM Etna

a b s t r a c t DEMs derived from LIDAR data are nowadays largely used for quantitative analyses and modelling in geology and geomorphology. High-quality DEMs are required for the accurate morphometric and volumetric measurement of land features. We propose a rigorous automatic algorithm for correcting systematic errors in LIDAR data in order to assess sub-metric variations in surface morphology over wide areas, such as those associated with landslide, slump, and volcanic deposits. Our procedure does not require a priori knowledge of the surface, such as the presence of known ground control points. Systematic errors are detected on the basis of distortions in the areas of overlap among different strips. Discrepancies between overlapping strips are assessed at a number of chosen computational tie points. At each tie point a local surface is constructed for each strip containing the point. Displacements between different strips are then calculated at each tie point, and minimization of these discrepancies allows the identification of major systematic errors. These errors are identified as a function of the variables that describe the data acquisition system. Significant errors mainly caused by a non-constant misestimation of the roll angle are highlighted and corrected. Comparison of DEMs constructed using first uncorrected and then corrected LIDAR data from different Mt. Etna surveys shows a meaningful improvement in quality: most of the systematic errors are removed and the accuracy of morphometric and volumetric measurements of volcanic features increases. These corrections are particularly important for the following studies of Mt. Etna: calculation of lava flow volume; calculation of erosion and deposition volume of pyroclastic cones; mapping of areas newly covered by volcanic ash; and morphological evolution of a portion of an active lava field over a short time span. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Airborne Light Intensity Detection and Ranging (LIDAR) is an active remote sensing system consisting of a laser scanner, a Global Positioning System (GPS) and an Inertial Navigation System (INS) generally mounted on a small aircraft. The laser scanner transmits brief laser pulses to the ground surface which are reflected or scattered back to the scanner. The instrument detects the returning pulses and records the travel time of the light from the laser scanner to the ground and back. The position of the aircraft when firing the pulse is determined using a differential GPS. Laser pulse directions are combined with aircraft roll, pitch and heading values determined by the INS to obtain range vectors from the aircraft to ground points. When these vectors are summed with the aircraft positions, they yield accurate coordinates of points on the surface of the terrain (e.g. Axelsson, 1999; Baltsavias, 1999; Wehr and Lohr, 1999; Wagner et al., 2006). There has been a significant increase in the use of LIDAR data to produce accurate DEMs needed for quantitative analyses and modelling in geology and geomorphology (Gritzner et al., 2001; Bishop et al.,

⁎ Corresponding author. Tel.: +39 0508311948. E-mail address: [email protected] (M. Favalli). 0169-555X/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.geomorph.2009.04.010

2003; White and Wang, 2003; Bolongaro-Crevenna et al., 2004; Ganas et al., 2005). In recent years, LIDAR technology has also been extensively applied in volcanology (Mouginis-Mark and Garbeil, 2005; Mazzarini et al., 2007; Davila et al., 2007; Harris et al., 2007; Neri et al., 2008; Fornaciai et al., in press). Accurate DEMs are required for the accurate morphometric and volumetric measurement of volcanic features (Mazzarini et al., 2005; Ventura and Vilardo, 2007). Moreover, updated detailed DEMs are of crucial importance in the study of hazardous surface phenomena such as landslides, debris flows, and floods; and on volcano flanks, pyroclastic density currents and lava flows (Favalli et al., 2009). LIDAR data are greatly affected by systematic errors introduced by the instrumental chain. These errors are sometimes wrongly considered to be negligible, and the total error in a DEM derived from LIDAR data is deduced from instrument specifications alone. Even in cases where systematic errors are considered, only basic corrections such as boresight misalignment are usually calculated. Most errors in LIDAR-derived DEMs can be attributed to systematic and random errors introduced by both the laser scanner system and the GPS/ INS unit during data acquisition procedures (Vosselman and Maas, 2001). These errors produce discrepancies in the overlap regions between neighbouring strips. The primary error sources are related to sensor position errors (mainly depending on inherent GPS errors), orientation errors (errors in the determination of yaw, pitch and roll of the airplane,

124

M. Favalli et al. / Geomorphology 111 (2009) 123–135

Fig. 1. Strips of the 2005 Mt. Etna LIDAR survey. Black points represent tie points used in the strip adjustment. Colours represent the number of overlapping strips.

those in the description of the sensor-airplane geometry, and those in the scanning mirror angle) and errors in range determination. Different approaches for the calibration of LIDAR data have been presented in recent years. They vary with respect to both the parameterisation of the errors in the laser data and the kind of observations that are used to eliminate these systematic errors (Pfeifer et al., 2005; Pfeifer and Briese, 2007). Laser strip adjustment parameterizations can be divided into two groups according to the space in which the error parameterization is performed: the object space or the measurement space. The first group of methods uses only the observed discrepancies in laser scanner data points from two or more overlapping strips; correction functions are determined for each strip, and the parameters of these functions are chosen so as to minimize discrepancies. Strip adjustments such as offsets and rotations or higher order transformations are modelled directly in the object domain (e.g. Kilian et al., 1996; Crombaghs et al., 2000; Kraus and Pfeifer, 2001). More rigorous calibration procedures are used in the second group of methods: modelling of systematic errors is performed directly in the measurement domain and is based on the modelling of the sensor system geometry, whereby each point is related to the original observation (Burman, 2002; Filin, 2003; Kager, 2004; Filin and Vosselman, 2004). The latter group of methods is superior, since it evaluates the discrepancies directly in the space of the measured quantities, allowing great improvement in strip matching using a limited number of parameters (e.g. bore-sight biases and range biases).

The method developed in this work follows the second approach. Errors affecting LIDAR data are identified by comparing two or more overlapping strips without requiring a priori knowledge of the surface. To quantify the discrepancies between strips, we select a number of tie points regularly distributed in the areas of overlap between strips and quantify strip mismatches using a TIN (Triangular Irregular Network)based method. Surface mismatches are calculated not only in the vertical direction, as is commonly done, but in all three directions. This is done because common sources of systematic errors introduce very minor vertical misalignments but much more significant horizontal discrepancies. In Section 4.2 we examine in detail the effects of various corrections on the discrepancies between strips and identify the corrections that minimize the discrepancies. The developed procedure was used to remove systematic errors in raw airborne LIDAR data from Mt. Etna surveys. In Section 5, improvements are assessed by comparing DEMs constructed from corrected LIDAR data with those built from uncorrected LIDAR data. Lastly, in Section 6, we apply the method to studies made at Mt. Etna which illustrate the importance of the corrections. We evaluate metric and submetric vertical changes in time by comparing either data from different surveys (time intervals of some years) or data from the same survey (time intervals of a few hours). We calculate the volume of the 2006 southern lava flow and the volumes of erosion and deposition of the 2002–2003 pyroclastic cones over the 2005–2007 period. We map the changes in the surface roughness of the summit area of the

M. Favalli et al. / Geomorphology 111 (2009) 123–135

125

Fig. 2. Vectors representing the displacement between overlapping strips (green and yellow areas) for the Mt. Etna 2005 survey (a) without corrections and (b) with corrections. The colours of the background image represent the number of strips (see Fig. 1 for legend).

volcano for the 2005–2007 time interval and, lastly, give an example of how it is possible to study and quantify the morphological evolution of an active lava flow using a temporal sequence of DEMs. 2. Error recovery model To remove the effect of systematic errors, the errors that affect LIDAR data must be identified and then the laser point locations must be parameterizes as a function of the variables which are the source of errors (α1, α2, α3…). Then the corrections in these variables that minimize the distance between the measured laser point locations and the real positions are given by: X i

2

jXi;laser ðα 1 ; α 2 ; α 3 N Þ − Xi;real j = min

ð1Þ

where the sum is over all the control points i, Xi,real is the location of the i-th control point and Xi, laser(α1, α2, α3…) is the point location of the i-th control point as measured by the LIDAR. Knowledge of the geometry of the data acquisition system is crucial for obtaining the best parameterization for removal of systematic error. A general and complete form for 3D coordinate transformation from the laser frame to the mapping frame is (Filin, 2003):   → → → → X l = X 0 + RW RG RINS δ + Rm Rs D

ð2Þ

where X 1 is the footprint position in the mapping frame (e.g. local coordinate system: WGS-84, UTM 33 North); X 0 is the location of the phase centre of the GPS receiver; R-terms are matrix rotations: RW is the rotation matrix from the local ellipsoidal system to the geocentric reference frame; RG is the rotation matrix from the reference system defined by the local vertical to the ellipsoidal reference frame; RINS is the rotation matrix from the body reference frame to the reference frame defined by the local vertical; δ is the offset vector between the phase centre of the GPS and the laser firing point defined by the body frame (boresight offset component); Rm is the mounting bias given by a rotation matrix between the altimeter and the body frame (boresight matrix); Rs is the rotation matrix between the laser beam and the laser system; and D = [0,0,ρ] is the range vector measured by the laser system, i.e. the 3D object coordinates in the laser frame. Eq. (2) is indeed a system of three equations, one for each of the x, y and z component of the vectors in it. Skaloud and Lichti (2006) used Eq. (2) to assess systematic errors due to the three bore-sight angles representing the orientation offsets that define the rotation matrix Rm and the range-finder offset Δρ. Filin and Vosselman (2004) estimated GPS biases to determine the magnitude of possible positional and height GPS offsets in data: X 0 in Eq. (2) is split into X 0 + δX 0i, where the vector δX 0 representing the three components of GPS bias is not assumed to be constant throughout the mission but is allowed to vary in time and, as a first approximation, independent offsets are assigned as a function of the strip number i. In the above examples, corrections were modelled by means of translations and rotations. In other cases the distortions may be nonlinear and different parameterizations might be required.

126

M. Favalli et al. / Geomorphology 111 (2009) 123–135

3. Tie point evaluation procedure The removal of systematic errors is generally complex. In practice, when dealing with laser mapping, it is impossible to know the exact location of the laser point, i.e. the footprint position in the real world. As point correspondences are practically impossible under such conditions, shape-based rather then traditional point-based algorithms are needed. Several authors recover the errors by collecting data along different flight paths over flat, locally horizontal surfaces and ‘flattening’ the surface as a function of the systematic errors. Some have generalized these procedures to make use of control surfaces of any slope. Others base their procedure on control height and tie points, in a fashion similar to a photogrammetric block adjustment, or reconstruct the elevation model around distinct landmarks to tie the overlapping strips. In our case we recover the errors directly from the selfconsistency of data collected for the Mt. Etna and Stromboli volcanoes without using any real control points. Systematic errors are detected and removed using only the distortions in the overlapping areas between different strips (Fig. 1). LIDAR data are usually collected along parallel overlapping strips, and the presence of perpendicular strips facilitates data correction. The 2005 Mt. Etna survey data analyzed in this paper consists of 34 NNE–SSW-trending strips covering a ground area of 622 km2, whereas the total area imaged by LIDAR, including overlapping areas, amounts to 1078 km2 (Fig. 1). Thus 58% of the surveyed ground area is covered by two or more strips. A total of 257,470,597 points were collected, for a total average point density of 0.41 pts m− 2. However, the density of ground points varies greatly for different strips and also along the same strip, depending on both the acquisition geometry (mainly the varying distance between the sensor and the ground) and the intensity of the backscattered signal (backscattering properties of the ground). For example, along a strip taken from 3700 m a.s.l. and passing over the summit of the volcano (∼ 3350 m), the density varies from an average of 0.21 pts m− 2 in areas far from the summit to an average of 1.05 pts m− 2 on the summit. To quantify the discrepancies between overlapping strips due to both systematic and random errors, we automatically generate a number of tie points on a regular mesh (Fig. 1). For each strip containing a tie point we locally reconstruct the surface, within a circle of given radius, using only data belonging to the strip. In this way if a tie point belongs to an area of overlap among n different strips, we have n independent surfaces representing the neighbourhood of the point. Surfaces are reconstructed in the form of TINs. Tie points are automatically discarded if the point density of the strip in that location is below a given threshold density. The radius and the density thresholds are selected according to survey characteristics so as to ensure good shape reconstruction. At each tie point, each strip is tested against overlapping strips. For each tie point, displacements Δx, Δy and Δz are calculated between each pair of strips containing the point using a surface to surface matching (Fig. 2). This method is not completely effective in flat horizontal areas, where only the z component of Eq. (2) is informative. The presence of any roughness at scales between the average inter-point distance and the radius scale is enough for the point to be a fully acceptable tie point. Areas characterized by small terrain undulations are ideal, since

they require very low spatial sampling. In our case, lava flows, scoria cones and other morphological features present on volcano slopes are good locations for tie points. However, forested or densely vegetated areas should be avoided. In our case, since the Mt. Etna flanks below 2100 m are densely vegetated and we want to automate the correction procedure, we add a weight for each tie point according to the residual of the matching procedure. Very densely vegetated areas for which the surface matching procedure gives high residuals are almost neglected due to their very low weight values. The use of a TIN to represent the surface makes the matching procedure sensitive to ranging noise, especially in regions with high data density. These points also have high residuals in the matching procedure and are therefore almost neglected due to their very low weight values. For example, for the 2005 Mt. Etna survey, our automatic tie point extraction procedure started from a 300-m step grid of tentative tie points; only 2547 of these points were located in regions of overlap between strips and with a sufficiently high point density on each strip (Fig. 1). 4. Calibration method 4.1. Error parameterization We rewrite Eq. (2) as → → → → X l = X 0 + R1 δ + R2 D

ð3Þ

by introducing new rotation matrices R1 = RW RG RINS and R2 = RW RG RINS Rm Rs. In this work we do not attempt to correct errors in X 0 or in the δ offset, even if our procedure can be used for their evaluation. Since the effects of angular inaccuracies are amplified by the object distance, the errors due to angular inaccuracy in the second term on the right side of Eq. (3) are negligible compared to those in the third term, which depend on the flying height (and so are at least three orders of magnitude greater). For the sake of simplicity, we can therefore rewrite Eq. (3) as follows: → → → X l = X f + Rðκ; /; ωÞ D ðρÞ

ð4Þ

where the X 1 coordinates of each point are given as the vector sum of the firing position X f of the given point (with X f = X 0 + R1δ ) and the range vector D (ρ) = [0,0,ρ] rotated by R(κ,ϕ,ω). Here κ and ϕ are the azimuth and pitch angles for the firing direction of the considered point. These angles are equal to those of the airplane if there is no mounting bias; however, there generally is an angular difference, which is taken into account by the matrix Rm in Eq. (2). The angle ω is the “roll” angle of the firing direction; it varies for each fired point and is equal to the sum of the roll angle of the airplane plus ωm (the ωcomponent of the mounting bias), plus the rotation angle of the laser beam (Rs in Eq. (2)). We consider corrections of the type: → → → X l = X f + Rðκ + δκ; / + δ/; ω + δωÞ D ðρ + δρÞ

ð5Þ

Fig. 3. Planimetric and altimetric analyses of displacements for uncorrected data (blue dots) and best-fit corrected data (red dots) for different error sources. a) planimetric displacements between strips based on uncorrected data. The green dashed line indicates the direction of flight. General corrections by Eq. (6): b) Comparison between planimetric displacements for uncorrected and corrected data; c) Altimetric displacements vs. angles of horizontal displacement for uncorrected and corrected data. Constant correction of the swapping angle (δω = αω): d) Schematic illustration of the effects of a constant error in the swapping angle; e) Comparison between planimetric displacements for uncorrected and corrected data; f) Altimetric displacements vs. angles of horizontal displacement for uncorrected and corrected data. Error in the swapping angle of the form δω = βωω: g) Schematic effects of a constant percentage error in the swapping angle; h) Comparison between planimetric displacements for uncorrected and corrected data; i) Altimetric displacements vs. angles of horizontal displacement for uncorrected and corrected. Error in the range parameter of the form δρ = βρρ: j) Schematic effects of a constant percentage error in the range; k) Comparison between planimetric displacements for uncorrected and corrected data; l) Altimetric displacements vs. angles of horizontal displacement for uncorrected and corrected data. This figure is discussed in Section 4.2.

M. Favalli et al. / Geomorphology 111 (2009) 123–135

127

128

M. Favalli et al. / Geomorphology 111 (2009) 123–135

where δκ,δϕ,δω and δρ represent the correction of κ,ϕ,ω and ρ, respectively. The most general corrections we consider are as follows: ρ

ρ

ρ

ρ

ð6aÞ

κ

κ

κ

κ

ð6bÞ

/

/

δρ = α + β ρ + βω ω + β/ / δκ = α + β κ + βω ω + β/ / /

δ/ = α + β / + βω ω ω

ω

ð6cÞ

ω

δω = α + β ω + β/ /

ð6dÞ

The α terms represent 0th order corrections. Considering only the 0th order corrections, if the mounting bias has already been properly accounted for, the angular 0th order corrections should be null. The angular 0th order corrections of Eq. (6) are not the mounting biases which are given, by definition, by expressing the corrections to the matrix in Eq. (5) in the form R(κ,ϕ,ω1)·R(δκ,δϕ,δω)·Rm(ωm) (see the matrix Rm in Eq. (2)) and since the rotation group is noncommutative, these values for δκ,δϕ,δω would be different from the ones we obtain. The β terms represent 1st order corrections: only first order corrections in the internal variables and in the ϕ and ω angles are taken into account. The α and β parameters in Eq. (6) are calculated from the minimization of the following expression: X i

2

wik wil jXi;k ðα; βÞ − Xi;l ðα; βÞj = min

ð7Þ

where the sum is over all the tie points i, k and l refer to two different strips, wik and wil are weights, and Xi,k and Xi,l are the point locations calculated using Eq. (5). Xi,k(α, β)−Xi,l(α, β) is the displacement between the strips k and l at the location of the tie point i, calculated using a surface to surface matching as explained in Section 3. The terms used in Eq. (6) represent the greatest contribution in reducing the errors in Eq. (7): 0th terms (constants) and 1st order terms in the angles ϕ and ω.

by the horizontal black line is incorrectly reconstructed as the slanted red line. As shown in Fig. 3e,f, a correction of the form δω=αω =const in Eq. (6) produces an insignificant reduction in errors. Fig. 3g shows the schematic effect of an error in the swapping angle of the form δω =βωω, i.e. underestimation or overestimation of ω by a constant percentage. The possible causes of such errors are discussed by Katzenbeisser (2003). In this case the displacements are again perpendicular to the flight direction, but the horizontal ground line is erroneously reconstructed as a curve (upward or downward concave, depending on whether ω is overestimated or underestimated). Fig. 3h shows how the form δω=βωω corrects the original data for this kind of systematic error: the initial bilobate distribution collapses to a cloud of points distributed around the origin. A correction in range of the form δρ=βρρ, shown schematically in Fig. 3j, produces a very similar reduction in planimetric errors (Fig. 3k). Although the two different error corrections shown in Fig. 3g,j produce very similar planimetric adjustments, they yield opposite vertical displacements. So while the plots in Fig. 3h,k cannot discriminate which is the best correction (δω, δρ or a linear combination of the two), the vertical adjustments (Fig. 3i,l) clearly reveal that the systematic vertical error is greatly reduced by the δω=βωω alone, whereas the δρ=βρρ correction increases the vertical errors by a factor three. The misestimation of ω in the form δω =βωω is the most evident systematic error and determines the bilobate shape in Fig. 3a. In Fig. 3c, the angular distribution of displacements in uncorrected data shows that they aggregate around the two angles perpendicular to the main flight direction. Each of these angular aggregations has two main peaks due to a constant error in the pitch angle of the form δϕ =αϕ. In contrast, the data corrected using Eq. (6), i.e. the red dots in Fig. 3c, have a more uniform angular distribution (i.e. only the random error component is present) and the vertical displacement values are visibly reduced. In this particular case and in general, the main contributions to systematic errors are due to errors in ϕ and ω. This is the reason why they are given particular emphasis in Eq. (6). 4.3. Error correction

4.2. Strip adjustment The arrows in Fig. 2 show the planimetric displacements between different strips for a small region of the 2005 Mt. Etna survey. Fig. 2a shows the displacements between the uncorrected strips, and Fig. 2b shows the displacements after correction using Eq. (6). The planimetric displacements between the strips for the uncorrected 2005 Mt. Etna survey data are plotted in Fig. 3. Fig. 3b and c compare planimetric and altimetric displacements before (blue dots) and after (red dots) correction. In Fig. 3c the vertical displacement is plotted against the counter clockwise angle that the planimetric displacement forms with the x-axis. Fig. 3a clearly shows that the planimetric displacements are mainly perpendicular to the flight direction (green line): they are distributed in two symmetric lobes perpendicular to the main flight direction. In Fig. 3c these displacements show up as an aggregation of the vertical displacement around the same angles. Displacements in this direction are mainly due to either errors in the swapping angle or errors in the range. The swapping angle is equal to ω (as defined in the previous section) minus the roll angle of the airplane. For the sake of simplicity, in the following discussion we assume that the roll angle of the airplane is negligible compared to the swapping angle, so that ω is equal to the swapping angle. Fig. 3d schematically shows the effect of a constant error in the swapping angle. This is, for example, the error introduced by roll mounting bias (e.g. Skaloud and Lichti, 2006). If we add a constant error δω=αω to the real ω, the position of the black tree in Fig. 3d is erroneously reconstructed as that of the red tree: an altimetric displacement is introduced, as well as a planimetric displacement perpendicular to the flight direction. Indeed, the true ground represented

The systematic correction described by Eq. (6) was applied to LIDAR data collected during three surveys performed with an Optech ALTM 3033 and during another performed with an Optech ALTM Gemini laser scanner. Sensor characteristics are summarized in Table 1. Table 2 lists the LIDAR surveys (the lines with bold characters in the table); for some of these surveys, separate systematic corrections were performed for each flight. For each survey and flight, the table reports the total number of collected data points, the number of strips, the number of evaluated tie points, RMS planimetric displacement (Δx) and RMS altimetric displacement (Δz) for the uncorrected data and for the 0th and 1st order corrections described by Eq. (6). All the surveys/flights, except those in 2007, have an uncorrected RMS planimetric (altimetric) error of around 2 m (0.4 m) and a fully corrected RMS planimetric (altimetric) error of around 0.5 m (0.2 m). In contrast, the 2007 survey has rather low uncorrected RMS planimetric and altimetric errors (Δxuncorr = 0.83 m and Δzuncorr = 0.11 m), and the full correction does not significantly reduce these errors (Δxcorr = 0.60 m Δzcorr = 0.09 m). This difference is due to the different Optech sensors used in the surveys: the 2007 survey employed an ALTM Gemini sensor, whereas all other surveys were performed using an ALTM 3033. For the all surveys with the ALTM 3033 sensor, the main source of error is underestimation of the ω angle, i.e. δω =βωω with βω equal to −0.0025 in all cases. This means that the ALTM 3033 sensor used in all three surveys seems to overestimate the laser firing angle by a fixed percentage (0.25%). The errors in Table 2 are calculated in the overlap region between strips: they are therefore an overestimation of the actual error. In fact, errors are mainly dependent on the swapping angle, which is higher in the overlapping areas between adjacent strips. Nevertheless, the mean error in the overlapping areas is a reliable indicator of DEM quality.

M. Favalli et al. / Geomorphology 111 (2009) 123–135 Table 1 Characteristics of the ALTM 3033 and Gemini systems. System

ALTM 3033

Operating altitude Elevation accuracy

265–3000 m 150–4000 m 15 cmb 1200 m 10 cm b 1000 m (1 sigma) (1 sigma) 25 cm b 2000 m 15 cm b 2000 m (1 sigma) (1 sigma) 35 cm b 3000 m 20 cm b 3000 m (1 sigma) (1 sigma) 1/2000×flying altitude 1/5.500 × flying altitude (1 sigma) (1 sigma) 1 cm 1 cm Variable from 0° to Variable from 0° to ±25° ± 20° Variable from 0 m to Variable from 0 m to 0.93 × altitude 0.68 × altitude Varies with scan angle Variable, maximum 70 Hz 1064 nm 1064 nm 33 kHz 33 kHz (max. altitude AGL 4000 m) 50 kHz (max. altitude AGL 3000 m) 70 kHz (max. altitude AGL 2500 m) 100 kHz (max. altitude AGL 2000 m)

Horizontal accuracy Range resolution Scan angle Swath width Scan frequency Laser wavelength Laser repetition rate

ALTM Gemini

Errors reported in Table 2 are not the commonly reported vertical root-mean-square errors (RMSE) (USGS,1986; Hodgson et al., 2003) but the root-mean-squared horizontal (Δx) and vertical (Δz) offsets calculated from the displacement vectors existing between surfaces reconstructed from different overlapping strips. Such offsets do not express point matching errors but rather surface matching errors. The reason for this choice is that the adopted correction procedure is based on the minimization of both positional and height offsets. The case of the Etna NW 2004 flight, the standard vertical RMSE between overlapping strips decreases from 0.55 m for the uncorrected data to 0.24 m for the corrected data.

5. Mt. Etna case study The model developed to remove systematic errors from raw airborne LIDAR data was tested on Mt. Etna (Sicily, Italy), the highest volcano in Europe. Although Mt. Etna is not characterised by highly explosive eruptions, the frequency of eruptions is very high. Besides an almost continuous activity at the summit craters, flank eruptions occur every few years, and potentially extensive lava flows can invade the densely inhabited lower slopes (Behncke et al., 2005). During the 20th century, the village of Mascali was reached in 1928; in 1971, a lava flow narrowly missed the village of Fornazzo; in 1981, a lava flow reached the outskirts of the village of Randazzo; in 1991–1993 a lava flow threatened the village of Zafferana (Barberi et al., 1993), and during the 2001–2003 flank eruptions of Mt. Etna, several tourist facilities were damaged and some villages on both the southern and northern flanks of the volcano were threatened (Behnke and Neri, 2003; Andronico et al., 2005). An updated and accurate DEM is of crucial importance for assessing lava flow hazard at Mt. Etna (Favalli et al., 2009). We here focus on three LIDAR data surveys completed at Mt. Etna: the first was performed on 16th and 17th September 2004, the second on 29th and 30th September 2005, both using an Optech ALTM 3033 laser altimeter, and the third on 20th and 21st June 2007 using an Optech ALTM Gemini laser altimeter (http://optech.on.ca). The LIDAR point density is highly heterogeneous in the Mt. Etna test area because the flanks of the volcano are covered by surfaces with different characteristics: vegetated and inhabited areas, very rough lava fields and areas covered with ash deposits. The intensity of backscatter from these surfaces varies significantly: vegetated and inhabited areas usually show very high backscatter intensity, whereas ash deposits are characterized by a very weak backscatter signal (Mazzarini et al., 2007; Fornaciai et al., in press).

129

The input raw data were corrected for systematic errors following the procedure described above. Delaunay triangulation of input points was used to create two sets of 2004, 2005 and 2007 DEMs, one from corrected data and the other from uncorrected data. Fig. 4 shows three examples of strip displacement effects and the resulting hillshaded relief of the 2005 DEM after systematic error correction. The line appearing in frames a and c occurs at the separation between two different strips, where there is a vertical misalignment of the overlapping zones. This offset causes the ‘comb effect' in frame e. The corrections applied to the data removed these artefacts (see frames b, d and f). Residual errors after correction were assessed by comparing the corrected and uncorrected DEMs for areas i) where no significant natural modifications occurred between 2004 and 2005: i.e. excluding the summit crater zone and the 2004 active lava flow (Mazzarini et al., 2005); ii) lacking a vegetation cover, so that the level of noise in the DEM is low, and iii) with sufficiently high point densities. Difference images were produced for both the corrected and uncorrected DEMs in order to compare them pixel by pixel (Fig. 5). Such images are widely used to visualize the spatial distribution of error between a DEM and a particular reference: the difference image removes the topography from the DEM and only shows the true deviation in elevation from the reference. In this case, the 2004 DEM was chosen as reference and differences between 2005 and 2004 DEMs were derived both for the uncorrected (Fig. 5a) and corrected (Fig. 5b) data. Fig. 5a shows NNE–SSW oriented error bands corresponding to mismatches in adjacent strips in either the 2004 or 2005 DEM. These error bands are absent in Fig. 5b, where they have been removed by the correction procedure, and persisting errors are mainly due to areas of low point density in one of the two DEMs. 6. Geomorphological applications The quality of DEMs strongly affects quantitative and qualitative geomorphological analysis. To highlight the importance of an accurate DEM for the measurement of volcanological features, we use both raw and corrected LIDAR data to calculate volume changes caused by lava flow emplacement during the 2006 Mt. Etna eruption (Fig. 6) and the amount of erosion and ash-scoria deposition on the 2002–2003 pyroclastic cones (Fig. 7) over the 2005–2007 period. The volume is estimated from the DEMs by summing elevation changes and multiplying this value by the DEM cell area (1 m2), whereas errors were assessed on the basis of volume differences in regions not affected by lava emplacement. The amount of lava measured corresponds to that of the lava field on the southern flank of the volcano, emitted from 27th October to 27th November 2006 by the vents at 3050 and 3180 m a.s.l. The volume change in the corrected DEMs was 4.56 ± 0.72 × 106 m3, Table 2 Summary of LIDAR surveys (lines with bold characters) and relative flights. Data

Δz

N

strips tie points

104,037,769 44,436,535 59,601,234

20 11 9

1012 319 693

ETNA 2004 2.09 0.39 0.73 0.27 0.55 0.16 ETNA NW 2004 1.63 0.32 0.67 0.26 0.39 0.11 ETNA SE 2004 2.27 0.42 0.64 0.27 0.47 0.15

257,470,597 34 150,286,028 22 107,184,569 12

2547 1371 1176

ETNA 2005 2.38 0.39 0.67 0.26 0.48 0.16 ETNA NW 2005 2.06 0.32 0.61 0.22 0.45 0.12 ETNA SE 2005 2.71 0.46 0.67 0.27 0.40 0.17

231,153,619

38

6635

6 3 3

904 129 112

6,201,783 3,063,288 3,138,495

N

Δx

N Total points

0th order

1st order

Δx

Δx

Δz

0.68 0.11

Δz

ETNA 2007

0.83 0.11

STR 2005 STR 2005 STR 2005

2.21 0.36 0.52 0.25 0.36 0.17 2.93 0.48 0.47 0.27 0.33 0.17 2.94 0.55 0.57 0.37 0.42 0.23

0.60 0.09

The number of points, strips and tie points are reported in the left columns, whereas horizontal and vertical RMS errors for uncorrected data and for different order corrections are reported on the right side. ETNA 2004 data were collected on 16th and 17th September 2004, ETNA 2005 data on 29th and 30th September 2005, ETNA 2007 data on 20th and 21st June 2007, and STR 2005 data refers to a survey of the Stromboli volcano (Sicily, Italy) completed on 30th September 2005.

130

M. Favalli et al. / Geomorphology 111 (2009) 123–135

Fig. 4. Three examples of DEM improvement in the case of the 2005 Mt. Etna survey. The effects of discrepancies between overlapping uncorrected strips are evident in the upper frames (a, c and e), whereas these errors are absent in the lower frames (b, d and f). NEC is the North East Crater.

whereas that in the uncorrected DEMs was 4.63 ± 0.85 × 106 m3. This relatively small difference is due to the large extent of the study area, where significant errors occur only in limited parts.

The importance of accurate DEMs becomes evident when we measure the topographic and volumetric change of small features such as the upper edges of pyroclastic cones built on the southern

Fig. 5. DEM difference between 2004 and 2005 data. a) difference between uncorrected 2004 and 2005 DEMs. The clearly visible NNE–SSW bands are due to errors between adjacent strips. b) difference between corrected 2004 and 2005 DEMs. The error bands disappear, and the residual errors are small except in the areas of very low data density in one of the DEMs.

M. Favalli et al. / Geomorphology 111 (2009) 123–135

131

Fig. 6. Lava field emplaced on the southern flank of the Mt. Etna volcano during the 2006 eruption.

flank of the volcano during the 2002–2003 Etna eruption. The 2007– 2005 DEM difference provided a grid dataset in which cells report the difference in altitude: a negative value corresponds to erosion, a positive value to deposition, and a very low or zero value to no change. Fig. 7 shows how the use of corrected DEMs dramatically changes the extent of both erosion and deposition areas. In Fig. 7c, the cones are split into two different deposition–erosion zones by a horizontal line coinciding with the boundary between different strips. This clearly artificial subdivision is reduced in Fig. 7d. Moreover, Fig. 7c shows an unnatural accumulation zone on the steeper eastern flank of the north cone; in Fig. 7d such a zone is reduced and largely converted into a more likely zone of erosion. Numerically, the volume of eroded material changes from 116 ± 32 × 103 m3 using the uncorrected DEMs to 137 ± 30 × 103 m3 using the corrected DEMs, whereas the volume of deposited material changes from 59 ± 14 × 103 m3 to 43 ± 8 × 103 m3, respectively. Raw and corrected LIDAR data were also used to generate maps of local topographic roughness, which is a measure of the texture of a surface. Several different parameters can be used to quantitatively model surface topographic roughness (Shepard et al., 2001). One of the most commonly used parameters is the root mean square (RMS) height around the mean or the standard deviation of height. We here calculate roughness using a 10 × 10 m moving window, which allows us to effectively calculate roughness up to a scale of 10 m. At each pixel, we

detrend the data within the moving window by subtracting the interpolating plane and then computing the RMS height inside the window. Subtraction of the interpolating plane allows us to measure roughness, irrespective of the average local slope. In this way, for example, the roughness of a plane is always 0, regardless of its inclination. Fig. 8b shows the roughness map for the summit area of Mt. Etna using the corrected 2005 LIDAR data. Smooth surfaces with a roughness of less than 12 cm consist of fine, loose material: pyroclastic or epiclastic ash. Lava flows have moderate to high roughness (12– 50 cm), whereas vegetated areas have very high roughness (N50 cm). Roughness, a measure of surface texture, is highly dependent on variations in elevation over short horizontal distances. In our case, outside the vegetated areas, surface roughness is of the order of a few decimetres; error bands due to the misalignment of neighbouring strips in the case of the uncorrected 2005 LIDAR data are thus easily detected as bands with higher roughness values (thin bands highlighted by the white arrows in Fig. 8a). Fig. 8c,d show the roughness maps created from uncorrected and corrected 2007 LIDAR data respectively. In the case of the 2007 data, vertical errors are very small (~11 cm) even for the uncorrected data, and the roughness map presents no artefacts. Comparison between Fig. 8b and d reveals that the roughness distribution in the summit area of the volcano changed significantly over the 2005–2007 period, mainly due to ash deposition in new areas. These regions are mapped in Fig. 8f (as red areas) where the local roughness

132

M. Favalli et al. / Geomorphology 111 (2009) 123–135

Fig. 7. Effects of LIDAR data correction. a) and b) DEMs of 2002–2003 Mt. Etna pyroclastic cones derived from 2004 and 2005 data, respectively, corrected for the systematic error. c) 2005–2004 difference between uncorrected DEMs. There is an evident E–W line corresponding to an artefact in the erosion/accumulation between the northern and southern portions of the cones. d) 2005–2004 difference between corrected DEMs. The E–W line is almost absent and the erosion/accumulation distribution appears natural.

decreases by more than 40 cm between 2005 and 2007. The smoothing of these surfaces is mainly due to ash infilling of pre-existing incisions. This phenomenon affects a significant portion of the area around the summit craters, especially on the northern flank. In Fig. 8e the same change in roughness is mapped using uncorrected data: some artefacts, such as the NNE–SSW bands, appear. This method for mapping new loose deposits in volcanic areas seems promising, but comparison between Fig. 8e and f highlights the need for accurate source data. The 2004 LIDAR survey was performed during an ongoing effusive eruption at Mt. Etna (Fig. 9a). The 2004 eruption started on 7th September from a WNW–ESE trending fracture near the summit area, at about 3000 m a.s.l. (Neri and Acocella, 2006). On 10th September a new vent erupted along the same fracture, giving rise to a channelled lava

flow (flow 1, Fig. 9b). Two lateral channels developed from the main flow (flows 1.1 and 1.2, Fig. 9b). A second flow (flow 2, Fig. 9b) formed south of flow 1 (Mazzarini et al., 2005). This active lava flow recorded during the LIDAR survey on the morning of 17th September extends across three of the nine NNE–SSW strips acquired on that day. Strips were collected at different times, and the overlap areas were thus recorded at two different times. It is therefore possible to study the temporal evolution of the lava field in the areas of overlap. Unfortunately, in our case, only strips 15 and 16 overlap significantly in the area of interest, and only the central portion of flow 1 and the terminal portion of flow 2 are imaged in both strips (Fig. 9a,b). Fig. 9c,d shows the difference between the DEMs created from the two strips in the case of uncorrected and corrected data respectively. Strip 15 was acquired at

M. Favalli et al. / Geomorphology 111 (2009) 123–135

133

Fig. 8. Roughness maps for: a) the uncorrected 2005 DEM, b) the corrected 2005 DEM, c) the uncorrected 2007 DEM and d) the corrected 2007 DEM. White arrows in a) highlight the NNE–SSW error bands due to the misalignment of neighbouring strips. Regions in which the local roughness decreases by more than 40 cm between 2005 and 2007 are mapped in red using uncorrected data in e) and corrected data in f).

6:13 a.m. GMT and strip 16 at 7:29 a.m. GMT. Elevation changes due to lava flow are of the order of a few meters, i.e. they have the same order of magnitude as the errors due to misplacement of uncorrected strips, which vary from close to zero along the NW margin of the overlap area to 1–2 m along the SE margin (Fig. 9c). After the correction of data these errors are almost completely removed (Fig. 9d), allowing us to correctly describe the morphological evolution of this portion of lava field in the 76-minute interval between the acquisition of strips. Fig. 9d shows two open-channel lava flows (flows 1.2 and 2) in which the bulk volume flux of the flowing lava varies greatly over the down-flow distance and thus in time. Similar behavior is commonly observed for open-channel flows at Mt. Etna (e.g. Bailey et al., 2006).

The characteristic distance between successive flux pulses is ~ 350 m for lava flow 1.2 and ~ 100 m for lava flow 2. All the lava flows are active to some extent. The flow in channel 1.1 has stopped almost completely: the lava front does not advance significantly, but near the bifurcation point there is still a small flux pulse moving along the channel. Although flow 1 also shows low activity along the main channel, the emplacement of a small overflow is clearly visible (white arrow in Fig. 9d). Flow 1.2 is much more active, it receives a large part of the flux pulses arriving along the main channel. All five lobes of flow 1.2 appear to be advancing; the main thrust was observed along the main axis of the channel, which advanced 23 m with an average speed of 18 m h− 1 and emplaced a

134

M. Favalli et al. / Geomorphology 111 (2009) 123–135

Fig. 9. LIDAR survey of an active lava flow: a) Location of the 2004 lava flow (red area) and of strips 15 and 16 of the 2004 LIDAR survey. SC = Summit Craters, VDB = Valle del Bove. b) 2004 active lava flow as recorded by the September 17th, 2004 LIDAR survey. The coloured areas refer to the areas covered by strips 15 and 16; the area of overlap between the two strips is shown in green. The black square is the selected area shown in c) and d). c) and d) Difference between the DEMs obtained from strips 16 and 15 using uncorrected and corrected data respectively.

volume of 1550 ± 250 m3, corresponding to an average volume flux of 0.34 ± 0.06 m3 s− 1. The main pulses along flow 1.2 are moving at relatively low velocities (200 m h− 1) compared to the pulse velocities usually recorded near the vents (up to 1 km per hour; Bailey et al., 2006). Flow 2 shows the same unsteady flow characteristics of flow 1.2, but the flux pulses appear to be smaller and the dispersed flow front is advancing very slowly (~4 m h− 1). Multitemporal LIDAR data on active lava flows acquired over short intervals are ideal for the study and detailed quantification of all morphological changes such as overflows, the advance of flow fronts, and volume flow pulses along channels. To properly assess morphological variations, however, the DEMs must be adequately corrected, as the comparison between Fig. 9c and d clearly shows. 7. Conclusions LIDAR data are frequently affected by systematic errors that are usually overlooked but relatively easy to remove. Despite the widespread use of this technique, a rigorous and well-established treatment of systematic errors is lacking. We have presented a methodology to detect, estimate and correct systematic errors affecting LIDAR datasets. In the described case studies, the correction parameters were defined as misestimation in roll, pitch, yaw and range. They were calculated minimizing the distance between a

series of tie points distributed in the area of overlap between adjacent strips. The advantages of this method are as follows: 1. Straightforward, fully automated software implementation is possible. 2. Errors can be recovered by analyzing the self-consistency of collected data. Systematic errors are detected and removed using only the discrepancies between the areas of overlap between different strips. Since this procedure does not require the presence of ground control points, LIDAR data from remote areas, as in the case of some volcanoes, can be corrected. However, for comparing laser scans acquired in different periods, some ground referencing may be needed. 3. The error recovery procedure is not based on minimizing the standard RMSE between DEMs but takes full advantage of discrepancies in all three spatial directions. 4. The illustrated procedure eliminates the causes, not merely the effects, of errors. In other words, modelling of errors is performed directly in the measurement domain and does not make use of adjustments such as offsets and rotations in the object domain. It is therefore possible to greatly improve data using only a limited number of parameters (e.g. bore-sight biases and range biases). 5. This calibration method can be applied to all common airborne laser scanners. As an example, we applied this method to LIDAR data collected with two different sensors: the ALTM 3033 and ALTM Gemini sensors.

M. Favalli et al. / Geomorphology 111 (2009) 123–135

6. The method greatly improves the quality and accuracy of DEMs derived from LIDAR data. Examples of the improvements in DEM quality were given in Section 5. We presented two volume estimates calculated as the difference between DEMs constructed using uncorrected and corrected data. The first example deals with the volume of lava emitted during the 2006 eruption on the southern flank of Mt. Etna. In this case the correction does not significantly improve volume estimation because the computation area is large and systematic errors within it are statistically negligible. In the second example, we measured the volume of material eroded from and accumulated on the 2002–2003 coalescent pyroclastic cones over the 2005–2007 period. The improvement due to data correction is evident: the artificial horizontal line coinciding with the boundary between different strips is greatly reduced. The third example shows how accurately corrected LIDAR-derived DEMs can be used to derive roughness maps at small scales (up to 10 m) and to study changes in surface roughness through time. In particular, we present a methodology to identify areas newly covered by ash. Lastly, the fourth example shows how adequately corrected LIDAR data can be used to assess the morphological evolution of an active lava flow. This is done by comparing multitemporal LIDAR data acquired over short intervals. We have reported an example relative to the 2004 lava field, for which we had two overlapping DEMs created from LIDAR data acquired 76 min apart on 17th September, 2004. This example shows how a very complex phenomenon such as an active lava flow can be quantitatively investigated using a time sequence of LIDAR data. Acknowledgments This work was partially funded by the Italian Dipartimento della Protezione Civile in the frame of the 2007–2009 Agreement with Istituto Nazionale di Geofisica e Vulcanologia — INGV. One of us (A.F.) benefited from the MIUR-FIRB project “Piattaforma di ricerca multidisciplinare su terremoti e vulcani (AIRPLANE)” n. RBPR05B2ZJ. References Andronico, D., Branca, S., Calvari, S., Burton, M., Caltabiano, T., Corsaro, R.A., Del Carlo, P., Garfì, G., Lodato, L., Miraglia, L., Murè, F., Neri, M., Pecora, E., Pompilio, M., Salerno, G., Spampinato, L., 2005. A multi-disciplinary study of the 2002–03 Etna eruption: insights into a complex plumbing system. Bulletin Volcanologique 67, 314–330. Axelsson, P., 1999. Processing of laser scanning data — algorithms and applications. ISPRS Journal of Photogrammetry and Remote Sensing 54 (2–3), 138–147. Bailey, J.E., Harris, A.J.L., Dehn, J., Calvari, S., Rowland, S.K., 2006. The changing morphology of an open lava channel on Mt. Etna. Bulletin Volcanologique 68, 497–515. Baltsavias, E.P., 1999. Airborne laser scanning: basic relations and formulas. ISPRS Journal of Photogrammetry and Remote Sensing 54, 199–214. Barberi, F., Carapezza, M.L., Valenza, M., Villari, L., 1993. The control of lava flow during the 1991–1992 eruption of Mt. Etna. Journal of Volcanology and Geothermal Research 56, 1–34. Behnke, B., Neri, M., 2003. The July–August, 2001 eruption of Mt. Etna (Sicily). Bulletin Volcanologique 65, 461–476. Behncke, B., Neri, M., Nagay, A., 2005. Lava flow hazard at Mount Etna (Italy): new data from a GIS-based study. In: Manga, M., Ventura, G. (Eds.), Kinematics and Dynamics of Lava Flows, Geol. Soc. Am. Spec. Pap, pp. 189–209. Bishop, M., Shroder, J., Colby, J., 2003. Remote sensing and geomorphometry for studying relief production in high mountains. Geomorphology 55, 345–361. Bolongaro-Crevenna, A., Torres-Rodríguez, V., Sorani, V., Frame, D., Ortiz, M.A., 2004. Geomorphometric analysis for characterizing landforms in Morelos State, Mexico. Geomorphology 67, 407–422. Burman, H., 2002. Laser strip adjustment for data calibration and verification. International Archives of the Photogrammetry and Remote Sensing 34 (3), 67–72. Crombaghs, M.J.E., Brügelmann, R., de Min, E.J., 2000. On the adjustment of overlapping strips of laseraltimeter height data. International Archives of Photogrammetry and Remote Sensing 33 (B3/1), 224–231. Davila, N., Capra, L., Gavilanes-Ruiz, J.C., Varley, N., Norini, G., Vazquez, A.G., 2007. Recent lahars at Volcán de Colima (Mexico): drainage variation and spectral classification. Journal of Volcanology and Geothermal Research 165, 127–141.

135

Favalli, M., Mazzarini, F., Pareschi, M.T., Boschi, E., 2009. Topographic control on lava flow paths at Mt. Etna (Italy): implications for hazard assessment. Journal of Geophysical Research 114, F01019. doi:10.1029/2007JF000918. Filin, S., 2003. Recovery of systematic biases in laser altimetry data using natural surfaces. ISPRS Journal of Photogrammetric Engineering and Remote Sensing 69, 1235–1242. Filin, S., Vosselman, G., 2004. Adjustment of airborne laser altimetry strips. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 34, 285–289 (Part B3). Fornaciai, A., Bisson, M., Landi, P., Mazzarini, F., Pareschi, M.T., in press. A LIDAR survey of Stromboli volcano (Italy): DEM-based geomorphology and intensity analysis. International Journal of Remote Sensing. Ganas, A., Pavlides, S., Karastathis, V., 2005. DEM-based morphometry of range-front escarpments in Attica, central Greece, and its relation to fault slip rates. Geomorphology 65, 301–319. Gritzner, M., Marcus, A., Aspinall, R., Custer, S., 2001. Assessing landslide potential using GIS, soil wetness modeling and topographic attributes, Payette River, Idaho. Geomorphology 37, 149–165. Harris, A., Favalli, M., Mazzarini, F., Pareschi, M.T., 2007. Best-fit results from application of a thermo-rheological model for channelized lava flow to high spatial resolution morphological data,. Geophysical Research Letters 34, L01301. doi:10.1029/ 2006GL028126. Hodgson, M.E., Jensen, J.R., Schmidt, L., Schill, S., Davis, B., 2003. An evaluation of LIDARand IFSAR-derived digital elevation models in leaf-on conditions with USGS Level 1 and Level 2 DEMs. Remote Sensing of Environment 84, 295–308. Kager, H., 2004. Discrepancies between overlapping laser scanning strips — simultaneous fitting of aerial laser scanner strips. International Archives of Photogrammetry and Remote Sensing 35 (B/1), 555–560. Katzenbeisser, R., 2003. About the calibration of LIDAR sensors. ISPRS Workshop “3-D Reconstruction from Airborne Laser-Scanner and InSAR, 8-10 October 2003, Dresden. Kilian, J., Haala, N., Englich, M., 1996. Capture and evaluation of airborne laser scanner data. International Archives of Photogrammetry and Remote Sensing 31 (B3), 383–388. Kraus, K., Pfeifer, N., 2001. Advanced DTM generation from LIDAR data. International Archives of Photogrammetry and Remote Sensing 34 (3/W4), 23–30. Mazzarini, F., Pareschi, M.T., Favalli, M., Isola, I., Tarquini, S., Boschi, E., 2005. Morphology of basaltic lava channels during the Mt. Etna September 2004 eruption from airborne laser altimeter data. Geophysical Research Letters 32, l04305. doi:10.1029/2004Gl021815. Mazzarini, F., Pareschi, M.T., Favalli, M., Isola, I., Tarquini, S., Boschi, E., 2007. Lava flow identification and aging by means of LIDAR intensity: the Mt. Etna case. Journal of Geophysical Research 112, B02201. doi:10.1029/2005JB004166. Mouginis-Mark, P.J., Garbeil, H., 2005. Quality of TOPSAR topographic data for volcanology studies at Kilauea Volcano, Hawaii: an assessment using airborne lidar data. Remote Sensing of Environment 96, 149–164. Neri, M., Acocella, V., 2006. The 2004–2005 Etna eruption: implications for flank deformation and structural behaviour of the volcano. Journal of Volcanology and Geothermal Research 158, 195–206. Neri, M., Mazzarini, F., Tarquini, S., Bisson, M., Isola, I., Behncke, B., Pareschi, M.T., 2008. The changing face of Mount Etna's summit area documented with Lidar technology. Geophysical Research Letters 35, L09305. doi:10.1029/2008GL033740. Pfeifer, N., Briese, C., 2007. Geometrical aspects of airborne laser scanning and terrestrial laser scanning. International Archives of Photogrammetry and Remote Sensing 36 (3/W52), 283–287. Pfeifer, N., Elberink, S., Filin, S., 2005. Automatic tie elements detection for laser scanner strip adjustment. International Archives of Photogrammetry, Remote Sensing 36 (3/W3), 79–84. Shepard, M.K., Campbell, B.A., Bulmer, M.H., Farr, T.G., Gaddis, L.R., Plaut, J.J., 2001. The roughness of natural terrain: a planetary and remote sensing perspective. Journal of Geophysical Research 106, 32,777–32,795. Skaloud, J., Lichti, D., 2006. Rigorous approach to bore-sight self-calibration in airborne laser scanning. ISPRS Journal of Photogrammetric Engineering and Remote Sensing 61, 47–59. USGS (United States Geological Survey), 1986. Standards for Digital Elevation Models, Open File Report 86-004. Ventura, G., Vilardo, G., 2007. Emplacement mechanism of gravity flows inferred from high resolution Lidar data: The 1944 Somma–Vesuvius lava flow (Italy). Geomorphology 95, 223–235. Vosselman, G., Maas, H.G., 2001. Adjustment and filtering of raw laser altimetry data. OEEPE Workshop on Airborne Laser Systems and Interferometric SAR for Detailed Digital Elevation Models, vol. 40, pp. 62–72. Wagner, W., Ullrich, 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 and Remote Sensing 60, 100–112. Wehr, A., Lohr, U., 1999. Airborne laser scanning — an introduction and overview. ISPRS Journal of Photogrammetry and Remote Sensing 54, 68–82. White, S.A., Wang, Y., 2003. Utilizing DEMs derived from LIDAR data to analyze morphologic change in the North Carolina coastline. Remote Sensing of Environment 85, 39–47.