B -spline deconvolution for differential target cross-section determination in full-waveform laser scanning data

ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

B-spline deconvolution for differential target cross-section determination in full-waveform laser scanning data Andreas Roncat a,∗ , Gunther Bergauer b , Norbert Pfeifer c a

Christian Doppler Laboratory ‘‘Spatial Data from Laser Scanning and Remote Sensing’’, Vienna University of Technology, Gußhausstraße 27–29, AT-1040 Vienna, Austria

b

Institute of Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, AT-1040 Vienna, Austria

c

Institute of Photogrammetry and Remote Sensing (IPF), Vienna University of Technology, Gußhausstraße 27–29, AT-1040 Vienna, Austria

article

info

Article history: Received 24 June 2010 Received in revised form 31 January 2011 Accepted 1 February 2011 Available online 1 March 2011 Keywords: Laser scanning Full-waveform Deconvolution Linear estimation

abstract In full-waveform laser scanning, short laser pulses are emitted and travel towards Earth and object surfaces. The sensor samples the waveform of the emitted pulse and its complete backscattered echo as a function of time. This technique allows for the three-dimensional reconstruction of the terrain, natural and man-made objects, and for the derivation of (geo-)physical quantities such as the differential target cross-section. The retrieval of the differential target cross-section requires deconvolution which is an illposed problem. In this study, we present a novel technique for the computation of the differential target cross-section using B-splines. This class of mathematical functions enables a well-posed linear approach for deconvolution. Furthermore, it is not dependent on the symmetry of the temporal profiles of the emitted laser waveform and the received echoes, as approaches previously suggested. In this paper, the algorithm for deconvolution is presented in detail and validated for both synthetic and real-world fullwaveform laser scanner data. © 2011 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved.

1. Introduction The acquisition of geometric information in the fields of photogrammetry, remote sensing, and engineering geodesy has been revolutionized by the advent of laser scanning. One of the main advantages is the direct, area wide observation of all three coordinates of object points, in contrast to the observation of two coordinates with passive imaging, and in contrast to the acquisition of single points by total stations. Additional advantages, e.g., seeing through the gaps of the foliage (‘‘canopy penetration’’) and measurement of the ground below a forest canopy, or the reduced dependency on clear-sky conditions for data acquisition, have promoted the proliferation of laser scanning. This paper is especially concerned with determining the range between the laser scanner and solid targets and other properties of these targets. For details about platforms, scanning mechanisms, and phase-based laser scanners we refer to the literature, e.g. Baltsavias (1999). Time-of-flight laser scanners emit a short pulse of laser energy with a duration in the range of a few nanoseconds. The time lag between emission of the pulse and detection of its echo is measured,

∗

Corresponding author. E-mail addresses: [email protected] (A. Roncat), [email protected] (G. Bergauer), [email protected] (N. Pfeifer).

which leads via the multiplication with the speed of light to the two-way travel distance (Rüeger, 1990). One nanosecond of time thus corresponds to a distance of approx. 15 cm. Due to the diffraction of the laser beam when passing through the optical system it has a certain divergence in the order of 0.1 mrad and higher (Jelalian, 1992). Additional divergence can be desired, e.g., for the classification of the device into a certain laser class and thus ensuring a certain level of eye-safety. As a consequence of the divergence, only a portion of the laser beam may illuminate a target. The remainder of the emitted pulse may proceed and illuminate further targets along the ray. The reflection of the laser beam at different targets causes a portion of the energy to travel back to the detector which records the corresponding echo. Absorption, specular reflection, and the direction characteristics of the scattering of the target account for differences in the strength of the detected echoes (Rees, 2001). In discrete-return laser ranging, hardware components are used for detecting echoes and for measuring their round-trip travel time (e.g. by methods such as constant fraction, threshold detection, etc., Kamerman, 1993). In full-waveform laser scanning, the entire waveform of the backscattered echo is recorded, sampled typically by 1 ns intervals (Mallet and Bretar, 2009). This digitized waveform may be used for computing the return time or related to object properties, e.g. tree height (Duong et al., 2008), land cover, icesheet roughness or biomass (see references in Section 2.2).

0924-2716/$ – see front matter © 2011 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). Published by Elsevier B.V. All rights reserved. doi:10.1016/j.isprsjprs.2011.02.002

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

The target cross-section, which is the product of reflectivity and illuminated target area, characterizes a target illuminated by the laser beam. Other target properties, e.g. depolarization behaviour and scattering at different wavelengths, are not observed with full-waveform laser scanning systems (yet). As it will be further detailed, the received echo is directly related to the convolution of the emitted waveform with the derivative of the target crosssection. This derivative with respect to range is the differential target cross-section. Inferring target properties, e.g., its range to the sensor, or its longitudinal spread within the transmitted laser beam, is therefore possible by analysis of the received waveform. Target cross-sections are (geo-)physical properties of scattering surfaces which do not depend on the characteristics of the employed instrument and can therefore be used for comparing data between two campaigns. It was also shown that cross-sections were linked to surface roughness (Hollaus and Höfle, 2010). As cross-sections are surface parameters, thresholds for classification can be based on target properties, rather than on the appearance in raw sensor data. In this article we present a new method for the estimation of the differential cross-section of targets acquired from echoes recorded by full-waveform laser scanning. This method does not make strict assumptions on the temporal profile of the emitted waveform or the reflecting surface, e.g. symmetry or requirements for signal periodicity as found in current approaches. Examples of asymmetric emitted waveforms can be found in Jutzi (2007). The presented approach specifically provides the product of reflectivity and illuminated area as a function of height (range), without a need to model individual scatterers along the beam. Improved modeling of the measurement process is seen as a prerequisite for improving the quality of derived products. B-splines allow the approximation of arbitrary curves and surfaces in a fast and efficient way. Thus, B-splines have gained huge attention in the fields of computer graphics and computer-aided geometric design over the past decades. In this approach, we use the flexibility of B-splines for curve approximation as well as the convolution properties of (uniform) B-splines for the estimation of the differential target cross-section. The observations are the recorded amplitude values of the emitted and the received waveform, the estimated unknowns are B-spline coefficients describing the differential target cross-section. After outlining the basic physics underlying this technique and presenting existing methods (Section 2), we provide, in Section 3, an in-depth description of our new proposed deconvolution technique. We follow this with worked examples for the deconvolution of synthetic data and for the deconvolution of real-world data and then discuss the results in Section 5. Our conclusions are presented in the last section. 2. Physical principles and related work 2.1. Range equation and target cross-section The microwave radar range equation, in short the range equation which applies to optical and microwaves equally, relates the emitted power PE to the received power at the detector PD (Jelalian, 1992): PD =

σ π D2 ηATM ηSYS . βE2 R2 4π R2 4 PE

(1)

Here βE is the divergence of the emitted beam, R is the range from the sensor to the target, σ is the effective target cross-section (in m2 ), D is the aperture diameter of the receiver optic, ηATM is the atmospheric transmission factor, and ηSYS the system transmission factor. The target cross-section is a product of the target area

419

(dA (m2 )), the target reflectivity (ε[]), and the factor 4π /Ω describing the scattering angle of the target (Ω (sr)) in relation to an isotropic scatterer (Jelalian, 1992):

σ =

4π

Ω

ε dA.

(2)

If the area upon which the laser beam is reflected is larger than the footprint of this beam, the target area becomes proportional to the square of the range, reducing the 1/R4 relation in Eq. (1) to 1/R2 . In the above it is assumed that the target is a flat surface orthogonal to the beam. A more detailed analysis, presented by Wagner et al. (2006) includes the temporal shape of the emitted signal PE (t ) and the spatial distribution of the targets along the laser beam σi (R), i.e. the differential target cross-sections. The target crosssection is the integral of the targets’ differential cross-sections. With some simplifications and omission of the transmission factors of Eq. (1) the i-th detected echo can be described as: PD,i (t ) ≈

D2 4π R4i βE2

∫

Ri +δ

 PE

Ri −δ

t−

2R

vg



σi (R)dR.

(3)

Here, vg is the group velocity of the laser ray, and each differential target cross-section σi is within the interval [Ri − δ, Ri + δ]. Eq. (3) includes the convolution of the emitted signal with the differential target cross-section, leading to the received waveform. The emitted waveform itself is often unknown. The output of the detector when receiving the emitted waveform is called the system waveform and can be recorded. It may be used instead of PE (t ) in Eq. (3) since the same detector records both the emitted and the received waveform. In this study, we will use the term emitted waveform because it makes the distinction between the signal internal to the laser scanner and the signal traveling through the atmosphere clearer. 2.2. Existing methods for target cross-section determination A number of methods has been published for analysing recorded waveforms, including methods for localization of an echo (Roncat et al., 2008) or for decomposing the waveform into Gaussian (Duong et al., 2008; Wagner et al., 2006; Hofton et al., 2000; Harding and Carabajal, 2005) or more generalized (Mallet et al., 2009) components. The parameters of these components can then be compared from one epoch to another or related to other observations, e.g. ice sheet roughness (Zwally et al., 2002), land cover (Duong et al., 2009), or forest parameters (Lefsky et al., 2002, 2007). These methods utilize parameters extracted from the received backscattered waveform, but the derived quantities are still significantly influenced by the used instrument and the general setup of the flight campaign. Thus, the estimated values do not represent parameters of the scattering surface and are therefore not comparable for laser scanners operating with, e.g., different lengths of the emitted pulses. In (Wagner et al., 2006) the target cross-section is ‘‘reconstructed’’.1 As mentioned before, the target cross-section is a (geo-) physical parameter of the scatterer itself. The assumptions are that the emitted waveform is a Gaussian function in time and that the scatterers have a Gaussian target cross-section, i.e. that their differential target cross-section is also a Gaussian function in time. In this case, the recorded echo is a Gaussian function in time, too. The variance of the echo in received waveform is the sum of the variances of the emitted waveform and the differential target cross-section.

1 We use the term ‘‘estimated’’ in order to emphasize that an estimation of unknown parameters based on observations is involved.

420

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

Jutzi and Stilla (2006) exploit that convolution in time domain corresponds to multiplication in frequency domain. The emitted waveform and the received waveform are therefore transformed to the frequency domain. The spectrum of the differential target cross-section is then found via the division of the spectra, i.e. received waveform divided by emitted waveform. In this approach, a Wiener Filter is applied for noise reduction in the frequency domain. The time-domain differential target cross-section is eventually obtained by the transformation from frequency to time. Wang et al. (2009) propose to make the deconvolution in the time domain. The emitted waveform, the received waveform, and the differential target cross-sections are treated as continuous, piecewise linear functions, which allows rewriting the deconvolution as a system of linear equations. To gain numerical stability in the solution and to handle noise the equation system is regularized. The approach suggested in Wagner et al. (2006) is routinely applied, but the symmetry requirements are violated typically by the emitted waveform as well as by the targets, leading therefore to biases in the estimated properties. For complicated received waveforms, the choice of the number of individual scatterers and their approximate positions remains difficult. Gaussian components (see e.g. example in Duong et al., 2009) may not necessarily relate to individual scatterers. The approach of Jutzi and Stilla (2006) directly reconstructs the differential target cross-section but does not provide quality measures of the estimated results of the deconvolution.2 Furthermore, since the actual deconvolution is performed in the frequency domain, after the inverse Fourier transform the result might be complex-valued. However, the underlying physical model requires real-valued results. 3. B-spline deconvolution

3.1. B-splines and B-spline curves A B-spline Bnl (u) is a piecewise continuous polynomial function of degree n with C n−1 continuity. The subscript l denotes the shift along the abscissa (u axis). B-splines can be defined in a recursive form known as Cox recursion (Farin, 2002) u − ul−1 ul+n−1 − ul−1

Bnl −1 (u) +

Bn1

ul+n − u n−1 B (u), ul+n − ul l+1

(4)

B0l (u) =

1 0

· · · ul−1 ≤ u < ul · · · elsewhere.

(5)

u − (l − 1)1u n1u

Bnl −1 (u) +

(l + n)1u − u n−1 Bl+1 (u) n1u

and B0l

(u) =



1 0

( u) ⊗

B01

(u) =

∫

∞

B1n−1 (τ )B01 (u − τ )dτ

(6)

with B01 (u) as above. As a consequence, this gives n +n +1

n

n

2 1 2 Bl11+l2 − 1 (u) = Bl1 (u) ⊗ Bl2 (u).

(7)

Since this paper focuses on (de-)convolution, we will mainly refer to Eq. (7) in the subsequent text. The degree of the resulting Bspline curve is the sum of the degrees of the convolution operands plus 1. Fig. 1 illustrates the construction of B-splines by repeated convolution (cf. Zorin and Schröder, 2000). A B-spline curve γ (u) = (γ1 (u), . . . , γm (u))⊤ ⊂ Rm of degree n is defined as a linear combination of B-splines

γ ( u) =

imax −

bi Bni (u).

(8)

i=1

The bi ∈ Rm are called control points. Fitting a B-spline curve to observations is a linear least-squares problem jmax −

 y(uj ) −

imax −

2 bi Bni

( uj )

→ min,

(9)

i =1

where in our case the y(uj ) denote the observed amplitude values of the recorded waveform at time uj and the bi denote the unknown control points. In the case of full-waveform laser scanning, the values of u are given in a ‘‘natural’’ way by the recorded time stamps, i.e. uj = 1, 2, . . . , jmax . The dimension of the bi is 1, i.e. bi = bi . With the knot vector spacing 1u ≥ 1, and uj as above, the minimization of Eq. (9) has a unique solution, if, additionally, imax 1u ≥ jmax . With 1u = 1, the system is not overdetermined. It is explicitly noted that no multiplicity at the end points of the knot vector is introduced, which would lead to the interpolation of the first and the last control point. 3.2. Convolution and deconvolution of B-splines n

The ul are the elements of the knot vector u = (. . . , u0 , u1 , u2 , . . .). In the case of a uniform knot vector, i.e. ul − ul−1 = 1u for all l and setting u0 = 0, Eqs. (4) and (5) are simplified to Bnl (u) =

B1n−1

−∞

starting with



(u) =

j =1

In this section, a B-spline-based deconvolution algorithm is presented in detail. We start with the basic definitions, present curve fitting of both emitted and received waveform, and introduce a least-squares deconvolution approach.

Bnl (u) =

respectively. The Bnl (u) have positive values in the interval ((l − 1)1u; (l + n)1u) and values of zero outside of this interval. The advantage of using uniform knot vectors is that the convolution of two B-splines with uniform knot vectors and the same equidistance 1u again yields a B-spline with a uniform knot vector of equidistance 1u. The uniformity of the knot vector enables the calculation of B-splines of such a kind with a different, yet equivalent recursive algorithm using convolution (Zorin and Schröder, 2000)

· · · (l − 1)1u ≤ u < l1u · · · elsewhere.

2 However, they introduce a Gaussian-mixture model for the reconstructed differential target cross-section after the deconvolution. For this calculation, quality measures are given.

n

As stated above, the convolution of two B-splines Bl11 (u) ⊗ n +n +1

2 Bl22 (u) with uniform knot vector gives the B-spline Bl11+l2 − 1 (u). Aiming at full-waveform laser scanning the emitted waveform (ew) and the differential target cross-section (cs) are convolved to result in the received waveform (rw). Terms will be in subscript with the corresponding abbreviation. The degrees new and ncs of the B-spline curves approximating the emitted waveform and the (yet unknown) differential target cross-section, respectively, can be chosen arbitrarily.3 The degree nrw of the approximating Bspline curve representing the received waveform is thus nrw = new + ncs + 1 (see Eq. (7)). The length of the knot vector of this curve is kmax = imax + jmax − 1.

3 However, literature on curve fitting focuses mainly on cubic B-splines. Degrees below 3 lead to curves which are too inflexible for fitting, whereas higher degrees do not provide additional advantages, unless C 3 and higher-order continuity is required.

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

B 10(u)

1

421

B 11(u)

1

u

u -2

-1

0

1

2

3

4

5

6

7

8

-2

9

(a) B01 (u).

-1

0

1

2

3

4

5

6

7

8

9

6

7

8

9

(b) B11 (u).

B 11(0.8 −τ) 1

2

B 10( τ )

1

0

B 1 (u) = B 1 (u) ⊗ B 1 (u)

1

u

τ

-2

-1

0

1

2

3

4

5

6

7

8

-2

9

(c) B01 (τ ) and B11 (0.8 − τ ).

-1

0

1

2

3

4

5

(d) B21 (u). .. . B 17(u)

1

u -2

-1

0

1

2

3

4

5

6

7

8

9

(e) B71 (u). Fig. 1. Construction of uniform B-splines by repeated convolution: Images (a) and (b) show uniform B-splines of degree 0 (B01 (u)) and 1 (B11 (u)). The convolution of these  two B-splines yields the B-spline of degree 2, B21 (u) (Image(c)): the gray-hatched area is equal to the integral B01 (τ )B11 (0.8 − τ )dτ and therefore to B21 (0.8). B21 (u) itself is shown in (d), B21 (0.8) is indicated by the gray circle. The convolution procedure can be repeated until the desired degree is reached, e.g. B-splines of degree 7 will be used in Section 4 to approximate received waveforms. Image (e) shows the B-spline of such a kind, B71 (u).

Defining the B-spline approximation of the emitted waveform as imax −

ε(u) :=

bi,ew Bi ew (u), n

i=1

the one representing the (still unknown) differential target crosssection as jmax −

κ(u) :=

bj,cs Bj cs (u), n

j=1

and the one representing the received waveform as kmax =imax +jmax −1

−

ρ(u) :=

bk,rw Bk rw (u) n

k=1

=

jmax imax − − 

bi,ew Bi ew (u) ⊗ bj,cs Bj cs (u) n



n





Thus, in theory the redundancy only depends on the number of knots in the emitted waveform, but not on the length of the received waveform. However, the number of knots of the received waveform, kmax , has to be ≥ imax to solve the equation system. The normal equation matrix of this system has a band form with a width of new + 1. The disadvantage of the adjustment technique described above is that the parameters of the curve representing the emitted waveform, bi,ew are treated as constants in Eq. (10), neglecting the fact that they were estimated in a least-squares adjustment. This can be overcome by an overall adjustment (i.e. curve fitting and deconvolution) following the general case for adjustment (Mikhail, 1976) after the single steps. In this case, the results of curve fitting and deconvolution serve as initial values for the parameters in the overall adjustment. Eq. (9) gives the setup for the observation equations whereas Eq. (10) gives the one for the constraints: A1 x = y + e

i=1 j=1

 A1 =

and, therefore (see Eq. (7))

ρ(u) =

jmax imax − − 

bi,ew bj,cs

i =1 j =1



n

n



new +ncs +1 Bi+ =Bni+rwj−1 j−1

bk,rw =



A2 = (Gcs , Gew , −I)



so that the whole equation system reads as follows:





bi,ew bj,cs .

(10)

i,j:i+j−1=k

The values of the control points of the differential target crosssection, bj,cs , can be determined using the scheme of Eq. (10). Interpreting this equation as an observation equation for least squares adjustment, the control points of the recorded waveform are the observations. The control points of the emitted waveform are known quantities, and the control points of the differential target cross-section are the unknowns. The observations are acquired in a regular sequence with the same device, thus equal weighting is suggested. The least-squares deconvolution has a redundancy of #observations − #unknowns = (imax + jmax − 1) − jmax

= imax − 1.

o Mrw

⊤ ⊤ x⊤ = b⊤ ew , bcs , brw

this yields for the values of the received waveform’s control points

−

A2 x = o o o



Bi ew (u) ⊗ Bj cs (u) ,



Mew o

A⊤ 1 A1 A⊤ 2

A2 o

  x

λ

⊤ y⊤ = y⊤ ew , yrw





 =





A⊤ 1 y o

(11)

with Mew · · · Jacobian Matrix of Eq. (9) w.r.t. unknown control points bew Mrw · · · Jacobian Matrix of Eq. (9) w.r.t. unknown control points brw Gew · · · Jacobian Matrix of Eq. (10) w.r.t. unknown control points bew Gcs · · · Jacobian Matrix of Eq. (10) w.r.t. unknown control points bcs Since observations and constraints are linear and bi-linear equations, resp., the algorithm is expected to converge fast and almost surely. The investigations on empirical data already carried out did not show evidence for divergence.

422

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428 0.9

0.8

noise-free σ = 0.01 σ = 0.02 σ = 0.05

0.7

0.8 0.7

0.6

0.6 amplitude [DN]

amplitude [DN]

0.5 0.4 0.3

0.5 0.4 0.3

0.2 0.1

0.2

0

0.1 0

-0.1 0

1

2

3

4

5

6

7

0

8

1

2

3

timestamp [ns]

4

5

6

7

8

9

timestamp [ns]

(a) Synthetic emitted waveform ε(u).

(b) Differential target cross-section representing asymmetric scatterer

κ(u). 0.7

0.6

amplitude [DN]

0.5

0.4

0.3

0.2

0.1

0 0

5

10

15

timestamp [ns]

(c) Differential target cross-section representing three symmetric scatterers κ(u). Fig. 2. Construction of synthetic examples for B-spline convolution and deconvolution. B-splines are shown as solid lines, while their individual components are depicted as dashed curves. Image (a) shows the emitted waveform, (b) shows the differential target cross-section of a synthetic asymmetric scatterer, and (c) the differential target cross-section of three synthetic symmetric scatterers. Fig. 3 shows the received waveforms corresponding to these data. ‘‘×’’, ‘‘+’’, ‘‘·’’ and ‘‘◦’’ markers indicate sampled values without and with added Gaussian noise of different intensity.

4. B-spline deconvolution for differential target cross-section determination In this section the deconvolution of received waveforms, i.e. the computation of differential target cross-sections, will be demonstrated. The radiometric calibration of the laser scanner (i.e. the determination of the constant factors of Eq. (1)), and the consideration of atmospheric effects (Höfle and Pfeifer, 2007) will not be treated here. First, two synthetic examples will be presented in Section 4.1. In both cases, the emitted waveform and the differential target cross-section are given. All of them are of degree 3, leading to the received waveform by convolution with degree 7. Second, we show a real-world example for B-spline deconvolution in Section 4.2 and analyse the same example by Gaussian Decomposition in Section 4.3. For both synthetic examples, the amplitudes were normalized in the form that the value of the respectively greatest bi,ew and bj,cs was set to 1. Both examples use the same emitted waveform which is slightly asymmetrical (see Fig. 2(a), solid line). Among its control points there are three with non-zero values (1u = 1 ns) so that the whole emitted waveform has a duration of 7 ns. The full width at

half maximum (FWHM) which is commonly used to characterize the pulse length is 1.8 ns. For the reconstruction of the differential target cross-section, both the synthetic emitted waveform and the synthetic waveform were sampled in 1 ns intervals. To present information about the stability of the proposed algorithm in the presence of noise, four sampling variants were made:

• noise-free (‘‘×’’ markers in Figs. 2 and 3) • added Gaussian noise with σ = 0.01 (‘‘+’’ markers in Figs. 2 and 3)

• added Gaussian noise with σ = 0.02 (‘‘·’’ markers in Figs. 2 and 3)

• added Gaussian noise with σ = 0.05 (‘‘◦’’ markers in Figs. 2 and 3). The Gaussian noise was added additively to the sampled values. For each variant, the investigation was based on statistic quantities, namely the square root of the variance a posteriori of the adjustment, s0 , the root mean square error, r.m.s., and the normalized root mean square error, r.m.s.norm . The latter are defined as follows:

 r.m.s.(f (u), g (u)) :=

1 umax − umin

∫

umax umin

(f (u) − g (u))2 du

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428 1.2

0.7 noise-free σ = 0.01 σ = 0.02 σ = 0.05

1

noise-free σ = 0.01 σ = 0.02 σ = 0.05

0.6 0.5 amplitude [DN]

0.8 amplitude [DN]

423

0.6

0.4

0.4 0.3 0.2

0.2 0.1 0

-0.2

0

0

2

4

6

8

10

12

14

16

18

-0.1 0

5

10

timestamp [ns]

15

20

25

timestamp [ns]

(a) Synthetic received waveform: asymmetric scatterer.

(b) Synthetic received waveform: three symmetric scatterers.

Fig. 3. Synthetic received waveforms ρ(u) resulting from the convolution of the emitted waveform and the differential target cross-sections of Fig. 2. The solid line in image (a) is the convolution of the solid lines in Fig. 2(a) and (b), whereas the solid line in image (b) is the convolution of the solid lines in Fig. 2(a) and (c). ‘‘×’’, ‘‘+’’, ‘‘·’’ and ‘‘◦’’ markers indicate sampled values without and with added Gaussian noise of different intensity. Table 1 Error analysis: Curve fitting of emitted waveform (a) and received waveforms (b), (c). (a) Emitted waveform

Table 2 Error analysis: Deconvolution (a) and forward modeling (b) of the received waveform containing one asymmetric echo. (a) Deconvolution

Noise level (σ ), Symbol

s0

r.m.s. (ˆε(u), ε(u))

r.m.s.norm (ˆε(u), ε(u))

Noise level (σ ), Symbol

s0

r.m.s. (κ( ˆ u), κ(u))

r.m.s.norm (κ( ˆ u), κ(u))

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0148 0.0167 0.0439

0.0000 0.0048 0.0126 0.0307

0.0000 0.0155 0.0411 0.0997

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0413 0.1643 0.0942

0.0000 0.0184 0.0639 0.0709

0.0000 0.0473 0.1646 0.1825

(b) Received waveform asymmetric echo

(b) Forward modeling

Noise level (σ ), Symbol

s0

r.m.s. (ρ( ˆ u), ρ(u))

r.m.s.norm (ρ( ˆ u), ρ(u))

Noise level (σ ), symb.

r.m.s. (ρ( ˆ u), ρ(u))

r.m.s.norm (ρ( ˆ u), ρ(u))

r.m.s. (ρ( ¯ u), ρ( ˆ u))

r.m.s.norm (ρ( ¯ u), ρ( ˆ u))

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0109 0.0203 0.0535

0.0000 0.0030 0.0161 0.0416

0.0000 0.0082 0.0437 0.1132

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0058 0.0296 0.0348

0.0000 0.0158 0.0806 0.0947

0.0000 0.0067 0.0304 0.0156

0.0000 0.0183 0.0820 0.0432

(ρ( ˆ u), ρ(u))

r.m.s.

r.m.s.norm (ρ( ˆ u), ρ(u))

0.0000 0.0079 0.0157 0.0522

0.0000 0.0339 0.0674 0.2237

(c) Received waveform three echoes Noise level (σ ), Symbol

s0

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0095 0.0103 0.0644

and r.m.s.norm (f (u), g (u)) :=

r.m.s.(f (u), g (u)) r.m.s.(g (u), 0)

• r.m.s. errors of the forward-modeled received waveform, i.e. convolution of fitted emitted waveform and reconstructed differential target cross-section ρ( ¯ u) = εˆ (u)⊗κ( ˆ u), in comparison to the synthetic (r.m.s.(ρ( ¯ u), ρ(u)), r.m.s.norm (ρ( ¯ u), ρ(u))) and fitted received waveform (r.m.s.(ρ( ¯ u), ρ( ˆ u)), r.m.s.norm (ρ( ¯ u), ρ( ˆ u))). See Tables 2(b) and 3(b), resp. 4.1. Deconvolution of synthetic data

,

resp. In the examples treated here, umin is always equal to 0 and umax is equal to 1uimax , 1ujmax or 1ukmax , resp. In detail, the error analysis consists of:

• s0 of the estimated emitted waveform (ˆε (u)) and received waveform (ρ( ˆ u)), root mean square errors r.m.s.(ˆε (u), ε(u)) and r.m.s.(ρ( ˆ u), ρ(u)), and normalized r.m.s. errors r.m.s.norm (ˆε (u), ε(u)) and r.m.s.norm (ρ( ˆ u), ρ(u)) of these fitted curves (see Table 1), where ε(u) and ρ(u) denote the synthetically generated emitted and received waveform

• s0 of the deconvolution (κ( ˆ u)) and root mean square errors r.m.s.(κ( ˆ u), κ(u)) and r.m.s.norm (κ( ˆ u), κ(u)) of the differential target cross-section reconstructed by deconvolution in comparison to the synthetic differential target cross-section (see Tables 2(a) and 3(a), resp.)

The differential target cross-section of one synthetic asymmetric scatterer is shown in Fig. 2(b) (full duration: 7 ns). It is highly asymmetric with a steeper ascending slope. The convolution of emitted waveform and differential target cross-section leads to an asymmetric synthetic received waveform (full duration: 14 ns; see Fig. 3(a)). The synthetic differential target cross-section representing three symmetrical scatterers is shown in Fig. 2(c) (full duration: 14 ns). The convolution of emitted waveform and differential target cross-section leads to a synthetic received waveform containing three echoes of different echo width (full duration: 21 ns; see Fig. 3(b)). The whole differential target cross-section is symmetric, but the asymmetric emitted waveform leads to an asymmetric received waveform. Table 1 gives an overview of the results of the curve fittings. In both cases, an error-free sampling in 1 ns intervals yields

424

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428 0.9

0.7

κ κˆ

0.8 0.7

0.5 amplitude [DN]

amplitude [DN]

0.6 0.5 0.4 0.3 0.2

0.4 0.3 0.2 0.1

0.1 0

0

-0.1

-0.1

0

1

2

3

4

5

6

7

8

9

0.04

0

5

0

5

10

15

10

15

0.1 κ− κˆ

0.02 κ− κˆ

κ κ ˆ

0.6

0

0.05 0

-0.05

-0.02

-0.1 0

1

2

3

4

5

6

7

8

9

timestamp [ns]

timestamp [ns]

(a) Asymmetric scatterer.

(b) Three symmetric scatterers.

Fig. 4. Top: Synthetic (κ(u), solid line) vs. reconstructed differential target cross-sections (κ( ˆ u), dashed line), noise level σ = 0.01. Bottom: Difference κ(u) − κ( ˆ u) of the two curves. 0.6

0.8 0.6 0.4 0.2

0.4 0.3 0.2 0.1 0

-0.2 0 0.05

-0.1 0.050

4

6

8

10

12

14

16

18 ρ− ρˆ

2

0 2

4

6

8

10

12

14

16

0 -0.05 0.05 0

2

4

6

8

10

12

14

16

18

0 -0.05

5

10

15

20

25

5

10

15

20

25

5

10

15

20

25

5

10

15

20

25

0 -0.05 0 0.05

18 ρˆ − ρ¯

ρˆ − ρ¯

ρ ρˆ ρ¯

0.5

0

-0.05 0 0.05

ρ− ρ¯

amplitude [DN]

ρ ρˆ ρ¯

1

ρ− ρ¯

ρ− ρˆ

amplitude [DN]

1.2

0 -0.05 0 0.05 0 -0.05

0

2

4

6

8 10 12 timestamp [ns]

14

16

18

(a) Asymmetric scatterer.

0

timestamp [ns]

(b) Three symmetric scatterers.

Fig. 5. Top: Synthetic (ρ(u), solid line), approximated (ρ( ˆ u), dashed line) and forward-modeled received waveforms (ρ( ¯ u), dotted line), noise level σ = 0.01. Bottom: Differences of the respective curves: ρ(u) − ρ( ˆ u), ρ( ˆ u) − ρ( ¯ u), and ρ(u) − ρ( ¯ u) (from top to bottom).

an approximation of the synthetic curves without noticeable deviation from the original curves, as was expected. Thus, in this case also the differential target cross-section can be reconstructed by deconvolution without significant deviation. In the case of curve fitting under noise, both the s0 and r.m.s. error increase according to the growing noise level. This is also valid for the deconvolution. For the asymmetric scatterer, Fig. 4(a) shows that the general shape, including the asymmetry, of the differential target crosssection is estimated well. The noise levels of 1% of the emitted waveform and the received waveform result in 5% of the deconvolved differential target cross-section (value of 0.0473 in Table 2(a)). In the forward modeling, however, the error level of 0.01 is maintained, indicating that the deconvolution remains a task with poorly defined optimal values. With increasing noise level, the general shape of the differential target cross-section is maintained (figures are not shown). However, the deviations grow according to the noise level. For the case of three symmetric scatterers, the estimated differential target cross-section contains three modes (Fig. 4). It is noticeable that the maxima of the synthetic, error-free and the estimated differential target cross-section deviate in the order of the noise level while the lateral deformation is more pronounced. The modes are shifted by less than 1 ns.

Table 3 Error analysis: Deconvolution (a) and forward modeling (b) of the received waveform containing three echoes. (a) Deconvolution Noise level (σ ), Symbol

s0

r.m.s. (κ( ˆ u), κ(u))

r.m.s.norm (κ( ˆ u), κ(u))

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0543 0.0978 0.1035

0.0000 0.0306 0.0463 0.0980

0.0000 0.1270 0.1919 0.4066

(b) Forward modeling Noise level (σ ), symb.

r.m.s. (ρ( ˆ u), ρ(u))

r.m.s.norm (ρ( ˆ u), ρ(u))

r.m.s. (ρ( ¯ u), ρ( ˆ u))

r.m.s.norm (ρ( ¯ u), ρ( ˆ u))

0 (×) 0.01 (+) 0.02 (·) 0.05 (◦)

0.0000 0.0104 0.0182 0.0558

0.0000 0.0446 0.0778 0.2390

0.0000 0.0067 0.0168 0.0172

0.0000 0.0294 0.0707 0.0713

The difference from the forward models ρ¯ to the fitted received waveforms ρˆ , Fig. 5, demonstrate that the beginning and the end of the signals are the weakest part in the determination. Compared to the example of one single scatterer, it is noted that the deviations increase for the three scatterer case.

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428 250

ε (u)

amplitude [DN]

εˆ (u)

425

Table 4 r.m.s. values of B-spline curve fitting and forward modeling for the example of Section 4.2.

200

(f (u), g (u))

r.m.s.(f (u), g (u))

r.m.s.norm (f (u), g (u))

150

(ˆε(u), ε(u)) (ρ( ˆ u), ρ(u)) (ρ( ¯ u), ρ(u)) (ρ( ¯ u), ρ( ˆ u))

18.96 0.483 0.479 0.084

0.169 0.039 0.039 0.007

100

Table 5 Results of Gaussian Decomposition for the example of Section 4.2. Echo ID

50

0

310

315

320

325

330

335

340

345

timestamp [ns]

30

ρ

(u) ρˆ (u)

Emitted waveform 0 Received waveform 1 2 3

Position (µ) (ns)

Amplitude (A) (DN)

σ (ns)

327.90

223.67

1.78

3701.70 3709.48 3720.44

9.32 22.09 7.30

1.39 4.00 2.81

amplitude [DN]

25

20

15

10

5

0 3690

3695

3700

3705

3710

3715

3720

3725

3730

3735

timestamp [ns]

Fig. 6. Emitted waveform (top) and received waveform (bottom) recorded with a Riegl LMS-Q560 instrument. The sampled values are depicted by the ‘‘×’’ markers, the B-spline approximations εˆ (u) and ρ( ˆ u) by solid curves.

4.2. Deconvolution of real-world data The next example (Fig. 6) is based on real-world data, acquired during a flight mission in early 2007 in Burgenland, Eastern Austria (Doneus et al., 2008; Roncat et al., 2010). The emitted waveform and the received waveform are recorded by a Riegl LMS-Q560 instrument (Riegl, 2010). The emitted waveform is unimodal and is modeled with a control polygon of length 2. For the mathematical reasons given before, we choose a cubic B-spline curve to model the emitted waveform and the differential target cross-section. Furthermore, empirical tests carried out with data of the same flight campaign showed that curve fitting performed best (in the sense of lowest s0 values) with cubic B-spline curves in the case of emitted waveforms. To suppress noise, the spacing of the B-spline control points is chosen as twice the sampling interval, here 2 ns. This leads to a full duration of the emitted waveform of 10 ns, i.e. the power is strictly zero outside this interval. The FWHM is 3.8 ns. The s0 of B-spline curve fitting for the emitted waveform is 17.5.4 Besides determining the vertices of the B-spline control polygon, an additional parameter has to be determined in the curve fitting due to hardware reasons: During A/D conversion in the waveform recording process of the LMS-Q560, a constant offset (DC offset) is added to the raw amplitude values since background radiation passing through the filter in front of the detector might

4 For comparing this result with the ones of Section 4.1, please note that the amplitudes here are higher by a factor more than one hundred than in the case of the synthetic examples.

cause negative amplitude values. However, these are outside of the allowed range of values of the recorder. Adding the DC offset keeps the amplitudes in the intended range of values which is typically [0 . . . 255] (1 byte). As a consequence of the choice of the degrees new = 3 and ncs = 3, the received waveform is modeled as B-spline curve of degree nrw = new + ncs + 1 = 7. According to the emitted waveform, control points are spaced in 2 ns intervals as well. With 12 non-zero control points in total, the received waveform has a duration of 38 ns, for which the amplitude values are larger than zero. The temporal profile of the received waveform is multimodal. Looking closer at the shape of the descending slope of the received waveform at approx. 3708–3715 ns suggests to assume two overlapping scatterers at this point. The distance of these scatterers is lower than the respective echo width so that no extra mode is visible. The same applies to the section from 3720–3730 ns (cf. Fig. 7). Again, a DC offset is determined within the least-squares adjustment for curve fitting. This curve fitting of the received waveform results in s0 = 0.54. The deconvolution yields a cubic B-spline curve representing the differential target cross-section with a length of 26 ns (see Fig. 7, left). According to Section 3.2, the results of the separated curve fittings and deconvolution serve as initial values for an overall adjustment following the principle of the general case for adjustment. This procedure gives the control polygons for the curves εˆ (u), κ( ˆ u) and ρ( ˆ u) (and the DC offsets) as results. Despite of the huge overlap of the scatterers in the received waveform ρ( ˆ u), the deconvolved B-spline curve κ( ˆ u) contains five clearly separated maxima, giving empirical evidence for the practicability of our inversion approach. The r.m.s. and r.m.s.norm values of curve fitting and forward modeling are summarized in Table 4. Both the curve fitting of the received waveform and its discrete forward model agree well with the originally recorded received waveform. There is practically no difference between the forward-modeled and the fitted received waveform (r.m.s.norm (ρ( ¯ u), ρ( ˆ u)) < 1%, see Table 4). 4.3. Comparison to Gaussian Decomposition The analysis is continued with the example of Section 4.2, but concentrate now on the deconvolution results of Gaussian Decomposition. A standard implementation of Gaussian Decomposition, suggested by Wagner et al. (2006), models the received waveform with three Gaussian scatterers. The detailed results of parameter estimation, i.e. curve fitting, are given below (see Fig. 8 and Table 5).

426

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428 0.08

250

κˆ (u) κˆ corr. (u) constr. adjustment

0.07

200 amplitude [DN]

0.06 amplitude [DN]

0.05 0.04 0.03

150

100

0.02 0.01

50

0 -0.01 3360

0 3365

3370

3375

3380

3385

3390

3395

3405

310

320

325

330

335

340

345

time stamp [ns] ρ(u ) ˆ (u ) ρ ρ(u )

25

30

20 25

15 amplitude [DN]

10 5

ρ− ρ

0 3690

3695

3700

3705

3710

3715

3720

3725

3730

3735

2

20

15

10

0 -2 3690

ˆ − ρ ρ

315

timestamp [ns]

30 amplitude [DN]

3400

3695

3700

3705

3710

3715

3720

3725

3730

3735

5

1 0

0 3690

-1 3690

3695

3700

3705

3710

3715

3720

3725

3730

3695

3700

3705

3710

3715

3720

3725

3730

3735

3735 time stamp [ns]

timestamp [ns]

Fig. 7. Top: Differential target cross-section retrieved by deconvolution, before (κ( ˆ u)) and after constrained adjustment and correction (κˆ corr. (u)). Bottom: Forward-modeled received waveform ρ( ¯ u) (dashed line) in comparison to the sampled values ρ(u) (‘‘×’’ markers) and the B-spline curve ρ( ˆ u) retrieved by constrained adjustment (solid line). Images at the bottom show the differences ρ(u) − ρ( ¯ u) and ρ( ˆ u) − ρ( ¯ u) (from top to bottom).

Fig. 8. Results of Gaussian Decomposition for the emitted (top) and received (bottom) waveform. The sampled values are depicted by the ‘‘×’’ markers, the dashed curves show the Gaussian functions representing three individual scatterers and the solid curve shows the sum of the Gaussian functions. The leftmost dashed curve shows the first Gaussian model which is very narrow, see also Fig. 9.

0.1

Ai,cs =

Ai,rw Aew

,

µi,cs = µi,rw − µew ,

2 σi2,cs = σi2,rw − σew .

(12) In the example presented here, one of the Gaussian functions representing individual scatterers has a σi,rw lower than the one of the emitted waveform, σew . In this case, deconvolution does not yield a ‘‘proper’’ Gaussian function any more, i.e. a function of the −x 2

x2

general form e , but one of the form e (see leftmost dashed curve in Fig. 9). The variance of the differential target cross-section is −0.15 ns. As a consequence, the integral of this function, the target cross-section, becomes incalculable for this scatterer, since the integral is unbound. This example shows one requirement for the Gaussian Decomposition approach: The model complexity, i.e. number of scatterers, has to be determined before fitting, either using heuristics or an optimization (e.g. Mallet et al., 2009, Duong et al., 2008 and Roncat et al., 2008). In this example, three scatterers were chosen because of three discernible peaks. The width of the second Gaussian is overestimated, and therefore the width of the first Gaussian

0.09 0.08 amplitude [DN]

The deconvolution of two Gaussian functions results in a Gaussian function again. The amplitude, position and variance of the Gaussian function representing the emitted waveform are 2 denoted as Aew , µew and σew and the ones representing a scatterer within the received waveform as Ai,rw , µi,rw and σi2,rw , resp. The deconvolution results in:

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0

3360 3365 3370 3375 3380 3385 3390 3395 3400 3405 time stamp [ns]

Fig. 9. Deconvolution results of Gaussian Decomposition of the example given in 2

2

Section 4.2. The leftmost dashed curve has ex form instead of e−x form since the 2 variance σi2,cs = σi2,rw − σew is negative. The solid curve shows the sum of the two other Gaussian functions representing a physically meaningful result.

becomes too small for being physically meaningful. As shown in Fig. 7, top, five scatterers all together are found by B-spline deconvolution. While a slight asymmetry in these scatterers can also be seen, especially their width and location are the basis for the better description after forward modeling (Fig. 7, bottom). The second Gaussian therefore does not describe a scattering surface.

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

427

5. Discussion

6. Conclusion

The results of Section 4 show that deconvolution based on Bsplines is feasible. The synthetic examples gave evidence that, on the one hand, both the curve-fitting and the deconvolution approach work without noticeable deviation in the absence of noise, and on the other hand, that in the presence of noise the algorithm remains stable. Noise was in the order of 5%. The B-splinebased deconvolution is neither restricted to symmetric emitted waveforms nor symmetric scatterers whereas the well-established approach of Gaussian Decomposition is dependent on these symmetry assumptions. The analysis of a real example containing a multi-modal received waveform confirms the feasibility of our deconvolution technique. The adjustment algorithm consists of two steps: 1. Sequential analysis. (a) Curve fitting of emitted waveform (Eq. (9)). (b) Curve fitting of received waveform (Eq. (9)). (c) Deconvolution (Eq. (10)). 2. Integrated analysis (Eq. (11)).

In this article, a novel deconvolution approach based on Bsplines for the retrieval of the differential target cross-section in full-waveform laser scanning data was presented. The algorithm includes least-squares B-spline curve fitting of the emitted waveform and of the received waveform, and convolution conditions including additionally the differential target cross-section. On behalf of simulated and real-world data, it was shown that the deconvolution is feasible, also in the presence of noise. In a real-world example, forward modeling of the received waveform by convolution of the emitted waveform with the estimated differential target crosssection resulted in an error below 1%. Compared to established approaches such as Gaussian Decomposition, our approach does neither need initial values for the adjustment nor depend on symmetric temporal profiles of emitted waveforms and/or scatterers. The method is therefore, from a theoretical point of view, preferable for deriving the differential target cross-section. Being a bi-linear approach, it is easier to implement than previously suggested approaches and provides also a direct solution of the inversion problem.

Whereas in steps 1(a) and (b) the originally recorded sampled values are treated as observations, in step 1(c) the B-spline parameters of the received waveform are treated as observables. This modus operandi ensures convergence of the adjustment in one step since the mathematical model is linear by definition. However, in step 1(c) the parameters of the emitted waveform are treated as constants, neglecting the statistics, i.e. the variance of these parameters, although not all of these parameters may have been determined equally precisely. This weakness can be overcome by step 2, an integrated analysis. It is an overall bi-linear equation system for solving the curve-fitting of the emitted and received waveform and the deconvolution after step 1, following the general case for adjustment. In this system, the already determined parameters are treated as approximation values of the unknowns bi,ew , bj,cs and bk,rw . It has to be mentioned at this point that the mathematical model does not make any restrictions to the sign of the determined parameters. This is no problem as long as the curve resulting from the deconvolution remains positive. In Fig. 4(b), small negative values of the differential target cross-section can be seen. In this case, it is clear that they are a result of measurement noise since the ground truth, i.e. the synthetic differential target cross-section, does not contain any negative parts. However, also ‘‘too narrow’’ echoes may appear, i.e. their echo width is smaller than the width of the emitted waveform. Currently, the only explanation of this phenomenon is that the echoes’ amplitudes are small in comparison to the noise level (as visible in the figures of Mücke et al., 2010). In Gaussian Decomposition, such small echo widths lead to negative variances σ 2 of the Gaussian function representing the differential target cross-section of this scatterer. As a consequence, the integral of this function, the target cross-section, is incalculable. The appearance of ‘‘too narrow’’ echoes also affects B-spline deconvolution and results in a differential target cross-section containing negative parts. Continuing the example of Section 4.2, physically meaningful differential target cross-sections must not be negative. One approach is to overcome this problem by sampling the retrieved differential target cross-section κ( ˆ u) with a short discretization interval (here 0.1 ns) and set the negative values to zero (see Fig. 7 (top), dotted curve). This new differential target cross-section, κˆ corr. (u), is now discretely convolved with an equivalently sampled version of the emitted waveform εˆ (u), yielding the discrete forward model of the received waveform ρ( ¯ u). Fig. 7 (bottom) shows the comparison of the originally recorded received waveform ρ(u), the result of the constrained adjustment ρ( ˆ u) and the discrete forward model ρ( ¯ u) = εˆ (u) ⊗ κˆ corr. (u). As visible in Fig. 7 and Table 4, this forward model is also well-suited to describe the received waveform. The difference in r.m.s. is about 1%.

Acknowledgements We would like to thank Wolfgang Wagner, Josef Jansa, Helmut Kager, Richard Kidd, Thomas Melzer (all IPF) and Andreas Ullrich (Riegl Laser Measurement Systems) for the interesting discussions and their helpful suggestions. Furthermore, we thank the reviewers for the very constructive contributions in improving the quality of this article. The ALS data acquisition in the Leithagebirge was part of the research projects ‘‘Celts in the Hinterland of Carnuntum’’ and ‘‘LiDAR supported archaeological prospection of woodland’’ which were financed by the Austrian Science Fund (P16449-G02; P18674G02). The first author has been supported by a Karl Neumaier Ph.D. scholarship. References Baltsavias, E.P., 1999. Airborne laser scanning: existing systems and firms and other resources. ISPRS Journal of Photogrammetry and Remote Sensing 54 (2–3), 164–198. Doneus, M., Briese, C., Fera, M., Janner, M., 2008. Archaeological prospection of forested areas using full-waveform airborne laser scanning. Journal of Archaeological Science 35 (4), 882–893. Duong, H.V., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2008. Single and two epoch analysis of icesat full waveform data over forested areas. International Journal of Remote Sensing 29 (5), 1453–1473. Duong, H.V., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2009. ICESat full waveform altimetry compared to airborne laser scanning altimetry over the netherlands. IEEE Transactions on Geoscience and Remote Sensing 47 (10), 3365–3378. Farin, G., 2002. Curves and Surfaces for CAGD. A Practical Guide, 5th ed. Morgan Kaufmann Publishers, San Francisco, USA. Harding, D., Carabajal, C., 2005. ICESat waveform measurements of within-footprint topographic relief and vegetation vertical structure. Geophysical Research Letters 32 (21), L21S10. doi:10.1029/2005GL023471. Höfle, B., Pfeifer, N., 2007. Correction of laser scanning intensity data: data and model-driven approaches. ISPRS Journal of Photogrammetry and Remote Sensing 62 (6), 415–433. Hofton, M., Minster, J., Blair, J., 2000. Decomposition of laser altimeter waveforms. IEEE Transactions on Geoscience and Remote Sensing 38 (4), 1989–1996. Hollaus, M., Höfle, B., 2010. Terrain roughness parameters from full-waveform airborne LIDAR data. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 38 (Part 7B), 287–292. Jelalian, A.V., 1992. Laser Radar Systems. Artech House, Boston. Jutzi, B., 2007. Analyse der zeitlichen Signalform von rückgestreuten Laserpulsen. Ph.D. Thesis. Technical University Munich. Jutzi, B., Stilla, U., 2006. Range determination with waveform recording laser systems using a wiener filter. ISPRS Journal of Photogrammetry and Remote Sensing 61 (1), 95–107. Kamerman, G.W., 1993. Laser radar. In: Fox, C.S., Accetta, J.S., Shumaker, D.L. (Eds.), Active Electro-Optical Systems. In: The Infrared and Electro-Optical Systems Handbook, vol. 6. SPIE Optical Engineering Press, (Chapter 1).

428

A. Roncat et al. / ISPRS Journal of Photogrammetry and Remote Sensing 66 (2011) 418–428

Lefsky, M.A., Cohen, W.B., Parker, G.G., Harding, D.J., 2002. LIDAR remote sensing for ecosystem studies. BioScience 52 (1), 19–30. Lefsky, M., Keller, M., Pang, Y., de Camargo, P., Hunter, M., 2007. Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveforms. Journal of Applied Remote Sensing 1 (1), 013537. doi:10.1117/1.2795724. Mallet, C., Bretar, F., 2009. Full-waveform topographic LIDAR: state-of-the-art. ISPRS Journal of Photogrammetry and Remote Sensing 64 (1), 1–16. Mallet, C., Lafarge, F., Bretar, F., Roux, M., Soergel, U., Heipke, C., 2009. A stochastic approach for modelling airborne LIDAR waveforms. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 38 (Part 3/W8), 201–206. Mikhail, E.M., 1976. Observations and Least Squares. IEP-A Dun-Donnelley, New York. Mücke, W., Briese, C., Hollaus, M., 2010. Terrain echo probability assignment based on full-waveform airborne laser scanning observables. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 38 (Part 7A), 157–162. Rees, W.G., 2001. Physical Principles of Remote Sensing, 2nd ed. Cambridge University Press. Riegl, 2010. Homepage of the company RIEGL Laser Measurement Systems GmbH. www.riegl.com (accessed: 31.01.11). Roncat, A., Bergauer, G., Pfeifer, N., 2010. Retrieval of the backscatter cross-section

in full-waveform LIDAR data using B-splines. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 38 (Part 3B), 137–142. Roncat, A., Wagner, W., Melzer, T., Ullrich, A., 2008. Echo detection and localization in full-waveform airborne laser scanner data using the averaged square difference function estimator. The Photogrammetric Journal of Finland 21 (1), 62–75. Rüeger, J.M., 1990. Electronic Distance Measurement: An Introduction, 3rd ed. Springer-Verlag, Heidelberg, Germany. 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 (2), 100–112. Wang, Y., Zhang, J., Roncat, A., Künzer, C., Wagner, W., 2009. Regularizing method for the determination of the backscatter cross section in lidar data. Journal of the Optical Society of America A 26 (5), 1071–1079. Zorin, D., Schröder, P., 2000. Subdivision for modeling and animation. SIGGRAPH 2000 Course Notes. URL: http://www.multires.caltech.edu/pubs/sig00notes. pdf (accessed: 31.01.11). Zwally, H., Schutz, B., Abdalati, W., Abshire, J., Bentley, C., Brenner, A., Bufton, J., Dezio, J., Hancock, D., Harding, D., Herring, T., Minster, B., Quinn, K., Palm, S., Spinhirne, J., Thomas, R., 2002. ICESat’s laser measurements of polar ice, atmosphere, ocean, and land. Journal of Geodynamics 34 (3), 405–445.