ANALYSING THE SUITABILITY OF RADIOMETRICALLY CALIBRATED FULL-WAVEFORM LIDAR DATA FOR DELINEATING ALPINE ROCK GLACIERS

With full-waveform (FWF) lidar systems becoming increasingly available from different commercial manufacturers, the possibility for extracting physical parameters of the scanned surfaces in an area-wide sense, as addendum to their geometric representation, has risen as well. The mentioned FWF systems digitize the temporal profiles of the transmitted laser pulse and of its backscattered echoes, allowing for a reliable determination of the target distance to the instrument and of physical target quantities by means of radiometric calibration, one of such quantities being the diffuse Lambertian reflectance. The delineation of glaciers is a time-consuming task, commonly performed manually by experts and involving field trips as well as image interpretation of orthophotos, digital terrain models and shaded reliefs. In this study, the diffuse Lambertian reflectance was compared to the glacier outlines mapped by experts. We start the presentation with the workflow for analysis of FWF data, their direct georeferencing and the calculation of the diffuse Lambertian reflectance by radiometric calibration; this workflow is illustrated for a large FWF lidar campaign in the Ötztal Alps (Tyrol, Austria), operated with an Optech ALTM 3100 system. The geometric performance of the presented procedure was evaluated by means of a relative and an absolute accuracy assessment using strip differences and orthophotos, resp. The diffuse Lambertian reflectance was evaluated at two rock glaciers within the mentioned lidar campaign. This feature showed good performance for the delineation of the rock glacier boundaries, especially at their lower parts.


INTRODUCTION
The field of airborne lidar, being a core source for high-resolution 3D topographic data, has been enriched in the past years by fullwaveform (FWF) digitizing systems in both geometric and physical aspects.Compared to "classical" discrete-return lidar data, the use of FWF data allows for a better target separation in case of multiple echoes per shot (Parrish et al., 2011) as well as for the derivation of additional echo features like amplitude, echo width etc. (Wagner et al., 2006, Mallet et al., 2011).By means of radiometric calibration of such echo features, physically meaningful target features may be derived (Wagner, 2010).
The research questions of this paper are the investigation of the relative quality of the radiometric calibration on the object.Furthermore, the suitability of the output of the radiometric calibration, i.e. the diffuse Lambertian reflectance (Briese et al., 2012), for the delineation of alpine rock glaciers is analysed.
The main contribution of this paper is the answer to the aforementioned research questions.Moreover, we demonstrate the direct georeferencing of the used FWF data and analyse its quality.The data used for this study was captured in the Ötztal Alps (Tyrol, Austria) on several subsequent days in October 2010 by an Optech ALTM 3100 instrument.While glaciers in this area have meanwhile been regularly surveyed by lidar for more than a decade now (Fritzmann et al., 2011), the dataset analysed in this paper was the first FWF lidar dataset in this area.
The relative geometric accuracy was assessed by means of strip differences, an absolute quality check was performed by comparison to orthophotos provided by the local authorities.
The paper is organized as follows: Section 2 contains the theoretic background and the formulas that were used for direct georeferencing and radiometric calibration in this study.The used dataset is described in the subsequent section.Results are given in Section 4 and the paper concludes with discussion and outlook in Section 5.

THEORY AND METHOD
In this section, the theoretical background and calculation formulas are given w.r.t. the extraction of echoes in FWF data (Section 2.1), direct georeferencing of FWF data captured by the sensor used in this study (Section 2.2) and radiometric calibration of such data (Section 2.3).

Extraction of echoes and echo attributes
The system waveform S(t), representing the outgoing laser pulse, and the recorded echo waveform Pr(t) were both modeled as a superposition of Gaussian functions.This technique is known as Gaussian Decomposition (Hofton et al., 2000, Wagner et al., 2006): ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey The system waveform is supposed to be uni-modal with the mode at position ts (in [ns]), a peak amplitude Ŝ (in [DN]) and a standard deviation of ss (in [ns]).The recorded echo waveform may consist of several Gaussians; the initial estimates for the mode positions ti and peak amplitudes Pi are found using classical detectors (Jutzi and Stilla, 2003) whereas sp,i is initially set to ss since it represents the lower bound of the meaningful values of sp,i.The fit is performed using non-linear Levenberg-Marquardt optimization (cf.(Wagner et al., 2006)).The final ( Ŝ, ss) and ( Pi, sp,i) amplitude/width pairs serve as pulse and echo features, resp.They will be input for radiometric calibration in a later stage enabling for the assignment of physically meaningful target features (see Section 2.3).
The distance ri [m] of a target is retrieved by ri = 0.299 792 458 • (ti − ts) 2 1 + 78.7 P 273.15 + T • 10 −6 + ∆r. ( The values P and T refer to the air pressure (in mbar) and temperature (in • C), resp.and are given in the header of the corresponding CSD file; these files contain the relevant information for direct georeferencing whereas the raw waveform data are stored in NDF files (Optech Inc., 2007).The range offset ∆r was empirically calculated from examples of single-echo cases with high signal-to-noise ratio.Comparing the ranges directly determined by the sensor and the ones retrieved by Gaussian decomposition resulted in ∆r = 0.65 m.

Direct georeferencing of Optech FWF data
Given the extracted ri range of an echo, we need further elements of direct georeference to calculate the echo's position vector in a superior coordinate system: • the scan angle α given in the corresponding CSD record, • the misalignment angles ψx, ψy and ψz given in the CSD header, • the IMU offset angles ωr,0, ωp,0 and ω h,0 given in the CSD header, • the raw roll, pitch and heading angle ωr,r, ωp,r and ω h,r given in the corresponding CSD record and • the ellipsoidal coordinates λ (longitude), ϕ (latitude) and h (elevation) of the laser origin, given in the corresponding CSD record.
With the Optech ALTM 3100 being an oscillating-mirror scanner, the range ri and scan angle α determine the echo's position vector in the sensor frame as follows (Parrish, 2007): rs = ri (0, − sin α, − cos α) . (3) Its position vector re in the Earth-centered Earth-fixed frame (ECEF) is retrieved by three subsequent rotations and the position vector re,0 of the laser origin in ECEF as translation: The definitions of the rotation matrices R • • and of the laser origin's position vector re,0 are given in Section A in the appendix.

Radiometric Calibration
The echo amplitude Pi allows for the following physical interpretation (Wagner et al., 2006): where Dr is the diameter of the receiver aperture, βt is the transmitter beamwidth, ηSYS and ηATM denote the system and atmospheric transmission factor, resp.The quantity σi (in m 2 ) is the target's backscatter cross-section, a physical target property in the dimension of an area.It actually corresponds to the target area effectively illuminated by the laser (Woodhouse, 2006, Wagner et al., 2006).
To solve the above equation for σi, it is re-grouped in the form that the unknown but constant quantities are separated from the others and summarized to the calibration constant CCAL (Wagner, 2010, Briese et al., 2012): Note that in some references (e.g.(Wagner et al., 2006)), the system waveform amplitude Ŝ is included in CCAL.However, empirical studies have shown that its value might vary significantly during a lidar campaign (Bretar et al., 2009, Roncat et al., 2011, Lehner et al., 2011), a fact which was also noticeable in the dataset investigated here (see Section 3).Thus, it is favourable to work with the normalized amplitude Pi/ Ŝ instead of Pi during radiometric calibration.
The backscatter cross-section is related to the backscatter coefficient γi and the diffuse Lambertian reflectance ρ d,i (both dimensionless) as follows (Woodhouse, 2006, Wagner et al., 2006): with A lf = (πr 2 i β 2 t )/4 denoting the laser footprint area and ϑ being the incidence angle between the laser ray and the local surface normal.The diffuse Lambertian reflectance is independent of the range ri, in contrast to σi, and of the angle ϑ, in contrast to γi.However, ρ d,i can only be reconstructed reliably in the case of extended targets, i.e. targets intercepting the whole laser beam, resulting in single echoes per laser shot.For perfectly diffuse targets, it is in the range [0, 1].Targets of specular reflection behaviour such as glass and water surfaces may exceed the upper limit.
The calibration constant is computed by means of homogeneous surfaces with known reflectivity (Wagner et al., 2006), artificial reference targets (Kaasalainen et al., 2009) or natural reference targets whose diffuse Lambertian reflectance is determined insitu using a reflectometer in the wavelength of the laser scanner (Briese et al., 2008).The latter approach was followed in this study.

DATA
The datasets presented in this study were acquired from October 7 to October 12, 2010 in the Ötztal Alps (Tyrol, Austria).The used system was an Optech ALTM 3100 instrument (Optech Inc., 2013), operated at a wavelength of 1064 nm and a pulse repetition rate of 70 kHz.FWF data were acquired in 68 flight strips.For quality control, additional flight strips with discrete-return scanning were acquired.For the overall area, a digital surface ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey model with 1 m resolution was derived, depicted by the shaded relief shown in Figure 1.For the investigation of the research questions mentioned in the Introduction, we concentrated on two rock glaciers within the project area: Ölgrube (approx.at 2 500 -2 700 m above the WGS 84 ellipsoid) and Äußeres Hochebenkar (2 300 -2 800 m elevation).The location of these areas is given in Figure 1 and a detailed view with the boundaries of the respective rock glaciers is shown in Figure 2.These boundaries are based on the study of (Krainer and Ribis, 2012).A first digitization of these outlines was performed on behalf of orthophotos from the year 2004 and lidar data from the year 2006, both provided by the Tyrolean government.Extended field surveys and expert knowledge allowed for verification and adjusting of the digitally mapped outlines.Additionally, the mapped outlines have been updated using an orthophoto, a shaded relief (acquisition time September 2010) and a difference layer of the digital terrain models of 2006 and 2010.Processing of the FWF data was done based on the method presented in Section 2. The width of the system waveform S(t) was around 10 ns (full width at half-maximum).As already mentioned, significant variations in the peak amplitudes Ŝ were noticed: the peak amplitude has a range of more than 10% of its average (see Figure 3).Therefore, each echo amplitude Pi was normalized w.r.t. the corresponding Ŝ before further processing.System Waveforms S(t) Figure 3: Overlay of 2 000 randomly chosen system waveforms S(t) of the flight strips acquired on October 10, 2010, centered at their peak position ts.The image illustrates the necessity to take the variation of the amplitude Ŝ into account for radiometric calibration.

RESULTS
FWF processing resulted in one 3D point cloud per flight strip with the additional radiometric attributes Pi/ Ŝ (normalized amplitude) and sp,i (echo width).As a first step, the ranging results from FWF analysis were compared against the ones delivered by the sensor, stored in the CSD files.A flight strip around the village Obergurgl was chosen for this purpose (UTM-32N coordinates around 654 000 E, 5 192 500 N.).Differences occured with a mean of 0.2 cm and a σMAD value (standard deviation of the mean absolute differences w.r.t the median) of 1.8 cm.The latter is a robust estimator for the standard deviation (Ressl et al., 2011).
Large differences were only found in areas of tall vegetation; however, in these areas ranging is much less reliable to compare than in open terrain.See Figure 4 for a visual impression of these results.The retrieved values were found acceptable for further processing.To reduce the discrepancies in relative alignment of the flight strips, a strip adjustment based on least-squares matching was performed using the software package OPALS (OPALS, 2013).
The σMAD value could be reduced from 11.5 cm to 5.5 cm with this procedure, as the histograms in Figure 5 show.However, the histogram of the differences after strip adjustment still show a significant skewness and a noticeable amount of differences being larger than of 15 cm absolute value.Rigorous strip adjustment such as the one presented in (Kager, 2004) may reduce these effects, as additional parameters such as range offsets and scale errors can be integrated in the adjustment model.
The absolute georeference was checked against orthophotos provided by the Tyrolean government.No significant shifts or distortions demanding for terrestrial measurements of control patches were noticed.From the such corrected overall point cloud, a digital surface model in the resolution of 1 m was computed.Moreover, local normal vectors could be reliably derived for each point, enabling for the calculation of the local incidence angle ϑ.The high variations in elevation and the steep slopes in the project area (see Figures 6 and 7) led to high variations in range ri and incidence angle ϑ.As outlined in Section 2.3, only the diffuse Lambertian reflectance ρ d is independent of both ri and ϑ; we therefore focused on ρ d in the investigations of this study.Local in-situ reference measurements of ρ d allowed for the calculation of this ρ d for each lidar point.These calculations were also done using OPALS (OPALS, 2013).Figures 6 and 7 show the digital surface models and maps of ρ d in the areas of the two selected rock glaciers.
The gridded data illustrated in these figures are also available online on the PANGAEA data platform as GeoTIFF files containing the digital surface model (1 m resolution), shaded relief (1 m resolution) and diffuse Lambertian reflectance (2 m resolution); see (Roncat et al., 2013).As the visualization of the diffuse Lambertian reflectances shows, in both cases especially the toe of the rock glaciers can be well delineated using this target feature.In a further analysis step of the performed radiometric calibration, the differences in diffuse Lambertian reflectance between overlapping flight strips were calculated in the same manner as the height differences for strip adjustment.The result is shown in Figure 8.
Concerning the two areas of interest, noticeable differences in the range of ±0.10 or of greater absolut value are visible at strip boundaries and at the outlines of the rock glaciers.While the first might be a result of incomplete strip adjustment, the latter may be referred to different reflection behaviour of the surfaces caused by the difference in acquisition dates and daytimes for the single flight strips.(As already mentioned, the whole campaign took 6 days.)

DISCUSSION AND OUTLOOK
In this study, the workflow for direct georeferencing and radiometric calibration of FWF lidar data was presented and evaluated on behalf of a large lidar campaign conducted in a high-alpine environment.We note that this is, according to our knowledge, the first complete description of how to perform direct georeferencing on the basis of Optech CSD and NDF files including the exploitation of FWF data for the determination of physically meaningful target properties.
Direct georeferencing was in the expected accuracy range and improvement by strip adjustment (using an approximative method) was possible.This resulted in strip differences with a σMAD of 5.5 cm.
Next to geometric calibration, also radiometric calibration was performed.It was shown that the emitted pulse has variations in amplitude in the order of 10%.Thus, the variation was considered in the computation of the radiometric object quantities.However, differences in the computed diffuse reflectance of overlapping strips was still in the order of 10% at some places.It was not further investigated if these differences originate in the radiometric measurement itself, i.e. if this is the measurement accuracy, or if these differences originate in the deviation from the Lambertian refection model, which is the basic assumption in the computation of the diffuse reflectance.
Two rock glaciers were chosen as areas of interest to test the suitability of the diffuse Lambertian reflectance, calculated by means of the FWF data, for delineating the outlines of these glaciers.In those areas the agreement was in the order of a few pixel.The grid width was 1 m.
However, in a both geometric and radiometric sense, remaining errors are noticeable; the correction of such errors on behalf of a rigorous strip adjustment and possibly terrestrial geometric and radiometric control measurements are considered as tasks of future research.
where e 2 is the squared excentricity of the WGS 84 ellipsoid and N = a/(1 − e 2 sin 2 ϕ) is the local radius of curvature in prime vertical; a denotes the major semi-axis.

Figure 1 :
Figure 1: Shaded relief of the project area, derived by means of the FWF data.The highlighted areas denote the rock glaciers of Äußeres Hochebenkar and Ölgrube which we will focus on in Section 4. Coordinates are given in UTM, Zone 32N.

Figure 2 :
Figure 2: Orthophoto of the areas of interest Ölgrube (top) and Äußeres Hochebenkar (bottom).The outlines of the rock glaciers are denoted by the black lines.Orthophoto: www.geoimage.atc , accessed on September 13, 2013.

Figure 4 :
Figure4: Differences in ranging between on-board range estimation of the sensor and FWF analysis using Gaussian Decomposition.

Figure 5 :
Figure 5: Histogram of strip differences before (top) and after strip adjustment (bottom) using least-squares matching.The amount of extreme differences and the σMAD value were greatly reduced by this procedure.

Figure 6 :
Figure 6: Digital surface model and shaded relief of the Ölgrube area (left) and diffuse Lambertian reflectance (right), derived by the results of FWF analysis.The toe of the rock glacier can be well distinguished from the surrounding area in the diffuse reflectance model.

Figure 7 :
Figure 7: Digital surface model and shaded relief of the Äußeres Hochebenkar area (top) and diffuse Lambertian reflectance (bottom), derived by the results of FWF analysis.Similar to Ölgrube (Figure 6), the outline in the lower part and the Western outline of the rock glacier is well represented in the reflectance image.

Figure 8 :
Figure 8: Absolute differences in diffuse reflectance between overlapping flight strips at Ölgrube (top) and Äußeres Hochebenkar (bottom).The relatively high differences at the outlines of the rock glaciers might stem from different acquisition dates and daytimes of the respective flight strips.