Decomposition of LiDAR waveforms by B-spline-based modeling

ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

Contents lists available at ScienceDirect

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

Decomposition of LiDAR waveforms by B-spline-based modeling Xiang Shen a,b,c, Qing-Quan Li a,⇑, Guofeng Wu a, Jiasong Zhu a a Key Laboratory for Geo-Environmental Monitoring of Coastal Zone of the National Administration of Surveying, Mapping and Geo-Information & Shenzhen Key Laboratory of Spatial-temporal Smart Sensing and Services, Shenzhen University, Nanhai Road 3688, Shenzhen 518060, PR China b Beijing Key Laboratory of Urban Spatial Information Engineering, Beijing 100038, PR China c College of Information Engineering, Shenzhen University, Nanhai Road 3688, Shenzhen 518060, PR China

a r t i c l e

i n f o

Article history: Received 14 September 2016 Received in revised form 13 March 2017 Accepted 15 March 2017

Keywords: Full waveform Laser scanning Light Detection And Ranging (LiDAR) Signal processing Spline

a b s t r a c t Waveform decomposition is a widely used technique for extracting echoes from full-waveform LiDAR data. Most previous studies recommended the Gaussian decomposition approach, which employs the Gaussian function in laser pulse modeling. As the Gaussian-shape assumption is not always satisfied for real LiDAR waveforms, some other probability distributions (e.g., the lognormal distribution, the generalized normal distribution, and the Burr distribution) have also been introduced by researchers to fit sharply-peaked and/or heavy-tailed pulses. However, these models cannot be universally used, because they are only suitable for processing the LiDAR waveforms in particular shapes. In this paper, we present a new waveform decomposition algorithm based on the B-spline modeling technique. LiDAR waveforms are not assumed to have a priori shapes but rather are modeled by B-splines, and the shape of a received waveform is treated as the mixture of finite transmitted pulses after translation and scaling transformation. The performance of the new model was tested using two full-waveform data sets acquired by a Riegl LMS-Q680i laser scanner and an Optech Aquarius laser bathymeter, comparing with three classical waveform decomposition approaches: the Gaussian, generalized normal, and lognormal distribution-based models. The experimental results show that the B-spline model performed the best in terms of waveform fitting accuracy, while the generalized normal model yielded the worst performance in the two test data sets. Riegl waveforms have nearly Gaussian pulse shapes and were well fitted by the Gaussian mixture model, while the B-spline-based modeling algorithm produced a slightly better result by further reducing 6.4% of fitting residuals, largely benefiting from alleviating the adverse impact of the ringing effect. The pulse shapes of Optech waveforms, on the other hand, are noticeably right-skewed. The Gaussian modeling results deviated significantly from original signals, and the extracted echo parameters were clearly inaccurate and unreliable. The B-spline-based method performed significantly better than the Gaussian and lognormal models by reducing 45.5% and 11.5% of their fitting errors, respectively. Much more precise echo properties can accordingly be retrieved with a high probability. Benefiting from the flexibility of B-splines on fitting arbitrary curves, the new method has the potentiality for accurately modeling various full-waveform LiDAR data, whether they are nearly Gaussian or non-Gaussian in shape. Ó 2017 Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS).

1. Introduction LiDAR (Light Detection And Ranging) has been known as one of most promising remote sensing techniques because of its excellent altimetric accuracy and canopy penetration capability (Baltsavias, 1999; Glennie et al., 2013; Nelson, 2013). Different from the lastgeneration discrete-return LiDAR systems that can only record a very limited number of echoes in a laser shot, full-waveform laser

⇑ Corresponding author.

scanners have the ability to record the complete waveform of the backscatter response, which offers the opportunity for researchers and end users to adopt advanced signal processing algorithms that can detect echoes with higher accuracy and reliability (Mallet and Bretar, 2009). More echo features (e.g., the pulse width and return energy) can also be extracted from waveform data (Hancock et al., 2015; Pirotti, 2011), and they have been shown to be greatly helpful in a variety of geoscience applications such as land cover classification (Fieber et al., 2013; Liu et al., 2015; Mallet et al., 2011; Tseng et al., 2015), building extraction (Michelin et al., 2012; Słota, 2015), canopy height retrieval (Gwenzi and Lefsky, 2014;

E-mail address: [email protected] (Q.-Q. Li). http://dx.doi.org/10.1016/j.isprsjprs.2017.03.006 0924-2716/Ó 2017 Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS).

183

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

technique is highly recommended for acquiring reliable solutions (Wang et al., 2009). As a special form of waveform decomposition, the Gaussian decomposition approach assumes that a LiDAR pulse can be modeled by the Gaussian function (Hofton et al., 2000). Owing to its conceptual simplicity and fairly good performance, the Gaussian decomposition has been very popular both in the academia and industry (Słota, 2014). A large number of studies have reported the successful use of the Gaussian decomposition algorithm for processing various LiDAR waveform data captured by satellite laser altimeters (Khalefa et al., 2013), airborne laser scanners (Yi et al., 2015), and terrestrial laser scanners (Hakala et al., 2012). However, some researchers have also pointed out that the Gaussian-shape assumption is not always satisfied for small-footprint LiDAR waveform data (Chauve et al., 2009; Hartzell et al., 2015; Mallet and Bretar, 2009), and they have suggested the use of the lognormal distribution (Chauve et al., 2007), generalized normal distribution (Chauve et al., 2009), Nakagami distribution (Mallet et al., 2010), Burr distribution (Mallet et al., 2010), or a piecewise exponential function (Hernandez-Marin et al., 2007) to model non-Gaussian waveforms in peaked and/or tailed shapes. Previous non-Gaussian waveform modeling methods hold the assumption that all LiDAR pulses fit an a priori probability distribution, and they therefore can only be suitable for processing the waveforms in particular pulse shapes. In this paper, we propose a more universal waveform decomposition approach. Given the high shape similarity between the returned echoes and the corresponding transmitted laser pulse, a received LiDAR waveform is treated as the mixture of finite transmitted waveforms after linear transla-

Hayashi et al., 2013; Nie et al., 2015), and biomass estimation (Allouis et al., 2013; Boudreau et al., 2008; Zhuang et al., 2015). Accurately extracting echoes from LiDAR waveforms is a challenging task. Currently, most of the sophisticated algorithms in the literature fall into two categories: waveform decomposition and deconvolution (Mallet and Bretar, 2009; Wang et al., 2015). In a waveform decomposition method, the backscattered signals are seen as a mixture of finite echo components in similar shapes. The number of the echoes and the initial shape parameters are first estimated (Hofton et al., 2000; Qin et al., 2012), and the echo parameters are then optimized by a non-linear least-square method such as the Levenberg-Marquardt algorithm (Tian et al., 2015; Wagner et al., 2006) or a statistical learning method, e.g., the Expectation-Maximization algorithm (Parrish et al., 2011; Parrish and Nowak, 2009) and the Reversible-Jump Monte-CarloMarkov-Chain algorithm (Hernandez-Marin et al., 2007; Mallet et al., 2010). The waveform deconvolution method, on the other hand, follows a quite different philosophy. A reflected waveform is seen as the convolution product between the transmitted laser pulse and the backscatter cross section (Mallet and Bretar, 2009; Wagner et al., 2006). By reversing the effect of transmitted waveform, a waveform deconvolution algorithm retrieves the backscatter cross section and then detects echoes from it. Typical deconvolution methods that have been successfully introduced into the full-waveform LiDAR field include the Wiener deconvolution (Jutzi and Stilla, 2006), Richardson-Lucy deconvolution (Wu et al., 2011), B-spline deconvolution (Roncat et al., 2011), and Gold deconvolution (Gao et al., 2015). Being an ill-posed problem, the signal deconvolution is inherently unstable, and the regularization

a)

250

b) 100

Transmitted waveform Gaussian fit

= 0.998

R

200

80

150

60

Amplitude

Amplitude

R

2

100

50

10

0

20

30

Received waveform Gaussian fit

= 0.994

40

Ringing effect

20

Ringing effect

0

2

0

40

0

20

40

c)

2500

R

2

R

2000

100

2

Received waveform Gaussian fit

= 0.983

1500

Amplitude

Amplitude

80

d)

2000 Transmitted waveform Gaussian fit

= 0.967

60

Digital bins

Digital bins

1500

1000

1000

500 500 0 0 0

20

40

Digital bins

60

80

0

50

100

150

Digital bins

Fig. 1. Typical airborne LiDAR waveforms and their Gaussian fit results. (a) A transmitted waveform and (b) a received waveform recorded by a Riegl LMS-Q680i laser scanner; (c) a transmitted waveform and (d) a returned topographic waveform captured by an Optech Aquarius LiDAR bathymeter.

184

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

Fig. 3. The workflow of the B-spline-based waveform decomposition algorithm. The marked procedures are the three differences from the Gaussian decomposition method.

Fig. 2. Spline modeling of an Optech LiDAR waveform pair. (a) The B-spline constructed from a transmitted waveform and its scaling result along the range and amplitude axes and (b) the spline curve is then translated horizontally, and the product closely fits to the first echo of the received waveform.

tion and scaling transformation. The B-spline curve-fitting technique is employed to model arbitrary shaped pulses, and the new waveform decomposition method can therefore be applied for a variety of LiDAR waveforms whether they are in regular or irregular shapes. The remainder of this paper is structured as follows. The next section briefly reviews the core idea of the waveform decomposition and provides a few examples of real LiDAR waveforms that cannot be adequately modeled by the classical Gaussian decomposition method. Then, in Section 3, we introduce the new B-splinebased modeling algorithm in detail. Finally, we provide comparative experimental results and conclude our work in the last two sections. 2. Waveform decomposition From the viewpoint of the waveform decomposition technique, a received LiDAR waveform is the superposition of finite echo components and the background noise (Mallet and Bretar, 2009). M X f ðx; hÞ ¼ f ðx; hm Þ þ Abn m¼1

ð1Þ

where M is the number of the components in the mixture, hm represents the descriptive parameters of the m-th model, and Abn refers to the mean noise level, which is usually estimated by averaging the amplitudes of the first or the last few waveform samples and is then removed from raw waveform signals (Duong et al., 2008; Hofton et al., 2000). The Gaussian decomposition, which is currently the most widely used algorithm both in literature and practice, is a special case of the waveform decomposition method by assuming that the laser pulses are all Gaussian in shape (Hofton et al., 2000; Wagner et al., 2006).

f ðx; hm Þ ¼ f ðx; lm ; rm ; am Þ G

G



¼ am e

ðxlm Þ2 2r2 m

ð2Þ

where the superscript G refers to the Gaussian model, and lm , rm , and am represent the mean, standard deviation, and amplitude of the m-th Gaussian component, respectively. The Gaussian decomposition method cannot always perform well on real LiDAR waveforms, because the Gaussian-shape assumption sometimes is overly simplistic. Fig. 1 shows two examples. The first two subgraphs illustrate a waveform pair (including a transmitted waveform and a received waveform) acquired by a Riegl LMSQ680i laser scanner. The waveforms were very accurately fitted by the Gaussian decomposition algorithm, except the tails of the echoes. The marked signals in the figure are the ringing artifacts (Lin et al., 2010; Roncat et al., 2008), and they cannot be modeled by the Gaussian function. The last two subgraphs represent typical examples of Optech Aquarius LiDAR waveforms. It can be found that both the transmitted and received waveforms are noticeably rightskewed and were poorly fitted by the Gaussian mixture model. 3. Method The signals received by a LiDAR system are the result of complex interaction between the outgoing laser properties and the ter-

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

starting with

Table 1 Experimental data.



Data set (Sensor)

Riegl LMS-Q680i

Optech Aquarius

Wavelength (nm) Emitted pulse width (ns) Time sampling (ns) Flying height above ground (m) Number of waveforms

1550 4.7 1 1000 69,766

532 8.3 1 300 83,119

I X ci bi;P ðxÞ

ð3Þ

i¼0

where ci is the control coefficient (often called the control point in the literature), P is the degree of the spline, and bi;P ðxÞ is the basis function of the B-spline and is given by a recursion formula (Gallier, 1999).

x  xi xiþpþ1  x bi;p ðxÞ ¼ bi;p1 ðxÞ þ biþ1;p1 ðxÞ xiþpþ1  xiþ1 xiþp  xi

bi;0 ðxÞ ¼

1 if xi 6 x < xiþ1 0

otherwise

ð5Þ

where xi is called the knot. In the proposed method, a cubic B-spline (i.e., the degree of the spline equals 3) is first created from the transmitted waveform (corresponding to the received waveform to be processed).

rain scattering (Hancock et al., 2015; Hartzell et al., 2015). The later one usually cannot be directly measured, and its effect in the waveform decomposition is implicitly modeled by alterable echo widths and amplitudes. In other words, a received waveform can be approximately seen to be the mixture of one or more standard probability distribution functions after linear stretches both in the range and amplitude directions. As shown in Fig. 1, the pulse shape of a received LiDAR waveform is highly similar to that of the corresponding transmitted waveform. It is therefore safely assumed that a received waveform can be accurately modeled by the superposition of many transmitted pulses after linear stretches. The proposed waveform decomposition algorithm is built on this assumption, and the B-spline technique is employed for fitting arbitrary shaped pulses and for facilitating the calculation of echo parameters (including the echo width, the amplitude, etc.). The B-spline method has been widely used in curve-fitting applications (Gruen and Akca, 2005; Guan et al., 2014). Its formula is given by

BðxÞ ¼

185

ð4Þ

f tw ðxÞ ¼ Bðx; xtw ; c tw Þ

ð6Þ

where xtw and ctw are the knots and the control coefficients, respectively. As shown in Fig. 2a, the B-spline knots xtw are first selected to uniformly cover the region of the transmitted pulse, and the control coefficients ctw are then calculated by fitting LiDAR waveform samples (the 2-D points constituted by sampling positions and signal amplitudes). Next, each detected echo in the received waveform is also modeled by a cubic B-spline. B

f ðxÞ ¼ Bðx; xm ; c m Þ

ð7Þ

where the superscript B refers to the B-spline-based model, xm and cm are the knot vectors and the control coefficient vectors (the values are unknown yet) of the m-th returned echo, respectively. As shown in Fig. 2, given the assumption that the echo shape of the received waveform can be approximated by the corresponding transmitted pulse after a linear translation and scaling, the knots and the control coefficients of the two B-splines defined by Eqs. (6) and (7) should practically meet

xm ¼ km xtw þ lm

ð8Þ

and

c m ¼ jm c tw

ð9Þ

respectively, where km and lm are the scale and translation coefficients in the range direction, respectively, and jm refers to the scale factor in the amplitudes. From Eq. (8), it can be seen that the echo shape in the B-spline-based modeling method is controlled by three

Fig. 4. Overview of the experimental areas. (a) Riegl and (b) Optech data sets. The point clouds are colored by elevation.

186

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

Fig. 5. The experimental results of the B-spline model. Above: Riegl data set; below: Optech data set. The point clouds are colored by the echo width (a and c) or signal amplitude (b and d).

variables, which is the same as the classical Gaussian decomposition (cf. Eq. (2)).

f ðx; hm Þ ¼ f ðx; lm ; km ; jm Þ B

B

ð10Þ

Fig. 3 illustrates the overall workflow of the proposed method. Beyond the aforementioned B-spline modeling procedure, two small differences between the B-spline-based method and the Gaussian decomposition should be pointed out. The first issue is how to estimate the initial values of ðlm ; km ; jm Þ in Eq. (10). Similar to the Gaussian decomposition method, the number of the echoes and their initial parameters

can be estimated through a peak detection algorithm (or a more sophisticated approach) in the preprocessing stage (Wang et al., 2015). The translation factor in the range direction can be initialized as

l0m ¼ P0m  Ptw

ð11Þ

where the superscript 0 represents that the parameters are approximate values, P0m is the rough position of the m-th echo, and Ptw is the peak position of the B-spline curve constructed by the transmitted waveform.The initial scale factor in the range is given by

187

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

a)

Gaussian

100

Generalized normal

Received waveform

Lognormal

140%

Gaussian

80

Generalized normal

120%

Lognormal

60

100%

Amplitude

Relative fit error

B-spline

80%

B-spline

40

60%

Ringing effect

20 40% 20%

0

0

20

40

60

80

100

Digital bins

0% Optech data set

Riegl data set

b)

Fig. 6. Statistical results of the least-square-fit residuals. The root mean square error of Gaussian decomposition results was served as a reference and was set to 100%.

80 Received waveform Gaussian

a)

Lognormal

Amplitude

Average width of transmitted pulses

0.12

Gaussian Generalized normal

0.1

Lognormal

Probability

Generalized normal

60

0.08

B-spline

40

20

B-spline

Ringing effect

0.06

0 0.04

0

20

40

60

80

100

Digital bins 0.02 0

Fig. 8. Waveform decomposition examples from the Riegl data set. (a) Two separate and (b) overlapping echoes exist in the waveforms.

5

4

3

7

6

8

9

10

Echo width (ns)

The scale factor of the amplitudes is initialized as

b) 0.025

j0m ¼

Average width of transmitted pulses Gaussian

Probability

Lognormal B-spline

0.015

0.01

PBm ¼ km Ptw þ lm W Bm ¼ km W tw

0.005

0

ð13Þ

where A0m and Atw are the maximal amplitudes of the m-th echo and the B-spline curve formed by the transmitted waveform, respectively. The second issue is to retrieve accurate echo properties from the model parameters ðlm ; km ; jm Þ. From Eq. (8), the following formulas can be easily derived.

Generalized normal

0.02

A0m Atw

ð14Þ

ABm ¼ jm Atw 5

6

7

8

9

10 11 12 13 14 15 16 17 18 19 20

Echo width (ns) Fig. 7. Histograms of the echo widths retrieved by different decomposition methods. (a) Riegl and (b) Optech data sets.

where W tw is the pulse width of the transmitted waveform and is defined by the FWHM (Full width at half maximum). 4. Experiments 4.1. Data

k0m

¼1

ð12Þ

which means that the pulse width of the m-th echo is initialized by that of the transmitted waveform.

Two airborne small-footprint LiDAR data sets were used in the experiments. The first one was acquired by a Riegl LMS-Q680i system, which is a near-infrared laser scanner with a high Pulse Repetition Frequency (PRF) up to 400 kHz. A small amount of data

188

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

Table 2 The detailed decomposition results of the two Riegl LiDAR waveforms in Fig. 8. The subscripts 1 and 2 represent two different echoes, P, W, and A refer to the position, width, and amplitude of an echo, respectively, R2 is the coefficient of determination in the least-square fit, and RMS is the root mean square difference between the LiDAR waveform samples and their corresponding fitted results. Data

Model

P1

W1

A1

P2

W2

A2

R2

RMS

Fig. 8a

Gaussian Generalized normal Lognormal B-spline

19.18 19.18 19.04 19.10

4.46 4.49 4.46 4.41

51.01 50.74 50.98 51.06

30.34 30.34 30.26 30.31

4.41 4.34 4.41 4.35

80.23 81.14 80.25 79.89

0.994 0.994 0.994 0.996

1.31 1.30 1.24 1.06

Fig. 8b

Gaussian Generalized normal Lognormal B-spline

17.23 17.03 19.59 17.75

5.68 5.38 11.62 6.65

19.68 18.42 22.43 19.95

23.62 23.59 23.65 23.69

5.08 4.79 4.06 4.85

64.97 68.38 50.56 63.20

0.986 0.988 0.991 0.990

1.53 1.44 1.23 1.29

were extracted from the whole survey area for the test, which contain 69,766 received waveforms and the corresponding emitted waveforms. The second data set was captured by an Optech Aquarius laser bathymeter. The instrument was operated at a PRF of 33 kHz. The survey area is over the coastal regions of an inhabited island located in the East China Sea. Only topographic waveforms were used in the experiments, because the processing of the water-column scatter contribution is beyond the scope of this work. Table 1 lists more technical information of the LiDAR systems and their configurations in the flight, and Fig. 4 illustrates the study areas by discrete LiDAR point clouds. It can be seen that many typical terrain features, such as deciduous trees and built structures, exist in both the two data sets.

4.2. Results and discussion Fig. 5 illustrates the point clouds produced by the B-spline based decomposition approach, colored by the echo width or amplitude. It can be seen that the LiDAR waveforms from bare earth and roads have relatively consistent echo widths, which are close to those of the corresponding transmitted waveforms (4.7 ns and 8.3 ns for the Riegl and Optech data sets, respectively). Widened pulses mainly exist in vegetation areas and building edges, and they usually have low intensity values. To analyze the performance of the B-spline-based modeling method, the following three waveform decomposition approaches were also employed in the comparative experiments.

b)

a) 80

80 Received waveform Generalized-normal fit

Received waveform Gaussian fit Individual components

40

20

0

Individual components

60

Amplitude

Amplitude

60

40

20

0

20

40

60

80

0

100

0

20

40

100

80 Received waveform

Received waveform Lognormal fit Individual components

B-spline fit Individual components

60

Amplitude

60

Amplitude

80

d)

c) 80

40

20

0

60

Digital bins

Digital bins

40

20

0

20

40

60

Digital bins

80

100

0

0

20

40

60

80

100

Digital bins

Fig. 9. Illustrations of the detected echoes of the LiDAR waveform in Fig. 8b. The decomposition results of (a) Gaussian, (b) generalized-normal, (c) lognormal, and (d) B-spline models.

189

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

a) 2000 Received waveform Gaussian Generalized normal

1500

Lognormal

Amplitude

B-spline 1000

500

0 0

20

40

60

80

100

Digital bins

b) 1200

Received waveform Gaussian

1000

Generalized normal Lognormal

Amplitude

800

B-spline

600 400 200 0 0

20

40

60

80

100

Digital bins Fig. 10. Waveform decomposition examples from the Optech data set. (a) Separate and (b) overlapping echoes exist in the waveforms.

(1) The Gaussian model. It is the most widely used method in the field (Hofton et al., 2000; Wagner et al., 2006). (2) The generalized normal model. It is especially suited to fit symmetric pulses in a peaked or flattened shape (Chauve et al., 2007). (3) The lognormal model. It is more suited to fit asymmetric pulses in a right-skewed shape (Chauve et al., 2007). Figs. 6 and 7 show the statistical results of least-square-fit residuals and retrieved echo widths, respectively. All waveform decomposition approaches used in the test, except the generalized normal model, yielded similar results in the Riegl data set. The Bspline based method achieved slightly better performance than

the Gaussian decomposition by further reducing 6.4% of residuals. As Optech LiDAR waveforms are remarkably right skewed (cf. Fig. 1), the accuracy of the Gaussian decomposition is not satisfied and much worse than that of the lognormal model. The mean fit residual of the B-spline method is 45.5% and 11.5% smaller than those of the Gaussian and lognormal decomposition models, respectively, which indicates that more accurate echo parameters can be retrieved. As shown in Fig. 7b, the pulse widths of most echoes retrieved by the B-spline model are close to those of transmitted waveforms, while the other three approaches produced overestimated values. The following parts of this section provide some waveform examples processed by different decomposition models. Fig. 8 illustrates two typical LiDAR waveforms from the Riegl data set and their waveform modeling products. The curves fitted by the four methods are almost the same in shape, except for the regions with the ringing artifacts. Because the ringing effect in the transmitted waveform can be effectively modeled by flexible spline curves, its negative impact on the received waveform was considerably alleviated in the B-spline-based decomposition method. However, the fit residuals are still noticeable, which shows that the relationship between the ringing artifacts in the transmitted and received waveforms is too complicated to be completely explained by a linear scaling. Table 2 shows the retrieved echo parameters and the leastsquare fit results of the waveforms in Fig. 8. By comparing the root mean square (RMS) difference (or the coefficients of determination R2 ) of the Gaussian and B-spline-based models, it can be concluded that the new method performs marginally better in processing Riegl LMS-Q680i LiDAR waveforms. The echo parameters retrieved by the two methods are highly consistent, with one remarkable exception. For the first echo in the waveform of Fig. 8b, the position difference is about 0.5 digital bins, which roughly corresponds to 0.15 m in the range. Given the more accurate least-square fit results, the echo properties retrieved by the B-spline based method are more reliable than those of the Gaussian decomposition method. It should be pointed out that for the waveform in Fig. 8b, the lognormal model yielded a slight better least-square-fit result when comparing to the B-spline model, but the retrieved echo parameters were highly unreliable. As shown in Fig. 9c and Table 2, the first echo of the lognormal model has an abnormally stretched width to fit the ringing artifacts, and it completely overlaps with the second echo. Two Optech waveform examples and their decomposition results are given in Fig. 10. The laser pulses are noticeably nonsymmetric (right-skewed), and they were poorly fitted by the Gaussian mixture model. The B-spline-based modeling method yielded much better fit results, especially in the starting and ending regions of laser echoes. It should be mentioned that the 60th–80th waveform samples in Fig. 10a were not well fitted by the four test algorithms. Even by manual interpretation, it is difficult to determine whether these signals come from a very weak

Table 3 Decomposition results of the waveforms in Fig. 10. All symbols have the same meanings as Table 2. Data

Model

P1

W1

A1

P2

W2

A2

R2

RMS

Fig. 10a

Gaussian Generalized normal Lognormal B-spline

27.79 27.79 27.40 27.03

9.22 8.75 9.18 8.28

458.37 473.28 461.92 480.56

52.58 52.56 52.34 51.78

10.12 9.93 10.18 9.23

1826.61 1851.48 1826.84 1899.15

0.983 0.983 0.991 0.997

37.90 37.78 28.04 16.99

Fig. 10b

Gaussian Generalized normal Lognormal B-spline

27.88 28.32 27.54 28.82

5.79 7.32 6.87 9.36

645.68 772.63 889.20 1197.83

36.02 37.00 36.44 38.25

16.26 13.23 14.30 10.71

1086.55 1130.92 1056.37 779.62

0.996 0.998 0.999 0.999

15.27 10.55 9.82 5.95

190

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191

b)

a) 1200

Received waveform Gaussian fit

1200

Received waveform Generalized-normal fit

1000

Individual components

1000

Individual components

Amplitude

Amplitude

800 600 400 200

800 600 400 200 0

0 0

20

40

60

80

0

100

20

40

60

Digital bins

Digital bins

c)

d)

80

Received waveform

1200

Received waveform

1000

Lognormal fit Individual components

1000

B-spline fit Individual components

800

Amplitude

Amplitude

1200

600 400 200

100

800 600 400 200 0

0 0

20

40

60

80

100

0

20

40

60

80

100

Digital bins

Digital bins

Fig. 11. Illustrations of the detected echoes of the LiDAR waveform in Fig. 10b. The decomposition results of (a) Gaussian, (b) generalized-normal, (c) lognormal, and (d) Bspline models.

echo or are the ringing artifacts accompanied with the strong pulse ahead. In the experiment, we simply assumed that it is an ineffective echo. Table 3 shows the detailed decomposition and least-square fit results of the Optech waveforms in Fig. 10. For the waveform shown in Fig. 10b, the echo properties extracted by the four test methods are highly inconsistent, which are further illustrated by Fig. 11. In the results of the Gaussian, generalized-normal, and lognormal models, the pulse widths of the two detected echoes are significantly different, and they are obviously unreasonable because the first and the second ones are much smaller and larger than that of the transmitted waveform (8.3 ns), respectively. The echo-position inconsistencies between the four test methods are noticeable. The differences between the lognormal and B-spline modeling results are as large as 1.28 ns and 1.81 ns for the first and the second echoes (which equal to 0.38 m and 0.54 m in the range), respectively. Given the smaller RMS values of fit residuals and more consistent echo widths, it can be safely concluded that the echo parameters retrieved by the B-spline method are more accurate and reliable than those of the other three models. 5. Conclusions This paper presented a new waveform decomposition method based on the B-spline modeling technique. The shape similarity between a returned echo and the corresponding transmitted pulse was observed, and a received LiDAR waveform can therefore be approximated by the mixture of finite transmitted waveforms after

linear translation and scaling transformations. B-splines are used for modeling the LiDAR waveforms in arbitrary shapes, and they can simplify the solution of retrieving echo parameters from irregular waveforms. The experimental results show that the B-splinebased modeling method can more accurately fit LiDAR waveforms and extract more reliable echo parameters than the traditional Gaussian decomposition method as well as the generalized normal and lognormal distribution-based methods, whether the laser pulses are in nearly Gaussian or noticeable skewed shapes. Benefiting from the high flexibility of B-splines on curves fitting, our proposed method can be universally used for processing the LiDAR waveforms in unknown distributions. Acknowledgments This work was supported in part by the Beijing Key Laboratory of Urban Spatial Information Engineering (No. 2014206), the National Natural Science Foundation of China (No. 41601493), the China Postdoctoral Science Foundation (No. 2015M570724), and the Shenzhen future industry development funding program (No. 201507211219247860). References Allouis, T., Durrieu, S., Vega, C., Couteron, P., 2013. Stem volume and above-ground biomass estimation of individual pine trees from LiDAR data: contribution of full-waveform signals. IEEE J. Sel. Topics Appl. Earth Observations Remote Sens. 6 (2), 924–934.

X. Shen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 182–191 Baltsavias, E.P., 1999. Airborne laser scanning: basic relations and formulas. ISPRS J. Photogramm. Remote Sens. 54 (2–3), 199–214. Boudreau, J., Nelson, R.F., Margolis, H.A., Beaudoin, A., Guindon, L., Kimes, D.S., 2008. Regional aboveground forest biomass using airborne and spaceborne LiDAR in Quebec. Remote Sens. Environ. 112 (10), 3876–3890. Chauve, A., Mallet, C., Bretar, F., Durrieu, S., Deseilligny, M.P., Puech, W., 2007. Processing full-waveform lidar data: modeling raw signals, International Archives of Photogrammetry. Remote Sens. Spatial Sci. 36 (Part 3/W52), 102– 107. Chauve, A., Vega, C., Durrieu, S., Bretar, F., Allouis, T., Deseilligny, M.P., Puech, W., 2009. Advanced full-waveform lidar data echo detection: assessing quality of derived terrain and tree height models in an alpine coniferous forest. Int. J. Remote Sens. 30 (19), 5211–5228. Duong, V.H., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2008. Single and two epoch analysis of ICESat full waveform data over forested areas. Int. J. Remote Sens. 29 (5), 1453–1473. Fieber, K.D., Davenport, I.J., Ferryman, J.M., Gurney, R.J., Walker, J.P., Hacker, J.M., 2013. Analysis of full-waveform LiDAR data for classification of an orange orchard scene. ISPRS J. Photogramm. Remote Sens. 82, 63–82. Gallier, J., 1999. Curves and Surfaces in Geometric Modeling: Theory and Algorithms. Morgan Kaufmann. Gao, S., Niu, Z., Sun, G., Zhao, D., Jia, K., Qin, Y.C., 2015. Height extraction of maize using airborne full-waveform LIDAR data and a deconvolution algorithm. IEEE Geosci. Remote Sens. Lett. 12 (9), 1978–1982. Glennie, C.L., Carter, W.E., Shrestha, R.L., Dietrich, W.E., 2013. Geodetic imaging with airborne LiDAR: the Earth’s surface revealed. Rep. Prog. Phys. 76 (8), 086801. Gruen, A., Akca, D., 2005. Least squares 3D surface and curve matching. ISPRS J. Photogramm. Remote Sens. 59 (3), 151–174. Guan, H.Y., Li, J., Yu, Y.T., Wang, C., Chapman, M., Yang, B.S., 2014. Using mobile laser scanning data for automated extraction of road markings. ISPRS J. Photogramm. Remote Sens. 87, 93–107. Gwenzi, D., Lefsky, M.A., 2014. Modeling canopy height in a savanna ecosystem using spacebome lidar waveforms. Remote Sens. Environ. 154, 338–344. Hakala, T., Suomalainen, J., Kaasalainen, S., Chen, Y.W., 2012. Full waveform hyperspectral LiDAR for terrestrial laser scanning. Opt. Express 20 (7), 7119– 7127. Hancock, S., Armston, J., Li, Z., Gaulton, R., Lewis, P., Disney, M., Danson, F.M., Strahler, A., Schaaf, C., Anderson, K., Gaston, K.J., 2015. Waveform lidar over vegetation: an evaluation of inversion methods for estimating return energy. Remote Sens. Environ. 164, 208–224. Hartzell, P.J., Glennie, C.L., Finnegan, D.C., 2015. Empirical waveform decomposition and radiometric calibration of a terrestrial full-waveform laser scanner. IEEE Trans. Geosci. Remote Sens. 53 (1), 162–172. Hayashi, M., Saigusa, N., Oguma, H., Yamagata, Y., 2013. Forest canopy height estimation using ICESat/GLAS data and error factor analysis in Hokkaido, Japan. ISPRS J. Photogramm. Remote Sens. 81, 12–18. Hernandez-Marin, S., Wallace, A.M., Gibson, G.J., 2007. Bayesian analysis of Lidar signals with multiple returns. IEEE Trans. Pattern Anal. Mach. Intell. 29 (12), 2170–2180. Hofton, M.A., Minster, J.B., Blair, J.B., 2000. Decomposition of laser altimeter waveforms. IEEE Trans. Geosci. Remote Sens. 38 (4), 1989–1996. Jutzi, B., Stilla, U., 2006. Range determination with waveform recording laser systems using a Wiener Filter. ISPRS J. Photogramm. Remote Sens. 61 (2), 95– 107. Khalefa, E., Smit, I.P.J., Nickless, A., Archibald, S., Comber, A., Balzter, H., 2013. Retrieval of savanna vegetation canopy height from ICESat-GLAS spaceborne LiDAR with terrain correction. IEEE Geosci. Remote Sens. Lett. 10 (6), 1439– 1443. Lin, Y.C., Mills, J.P., Smith-Voysey, S., 2010. Rigorous pulse detection from fullwaveform airborne laser scanning data. Int. J. Remote Sens. 31 (5), 1303–1324. Liu, C.X., Huang, H.B., Gong, P., Wang, X.Y., Wang, J., Li, W.Y., Li, C.C., Li, Z., 2015. Joint use of ICESat/GLAS and landsat data in land cover classification: a case study in Henan Province, China. IEEE J. Sel. Topics Appl. Earth Observations Remote Sens. 8 (2), 511–522.

191

Mallet, C., Bretar, F., 2009. Full-waveform topographic lidar: state-of-the-art. ISPRS J. Photogramm. Remote Sens. 64 (1), 1–16. Mallet, C., Bretar, F., Roux, M., Soergel, U., Heipke, C., 2011. Relevance assessment of full-waveform lidar data for urban area classification. ISPRS J. Photogramm. Remote Sens. 66 (6), S71–S84. Mallet, C., Lafarge, F., Roux, M., Soergel, U., Bretar, F., Heipke, C., 2010. A marked point process for modeling lidar waveforms. IEEE Trans. Image Process. 19 (12), 3204–3221. Michelin, J.C., Mallet, C.m., David, N., 2012. Building edge detection using smallfootprint airborne full-waveform lidar data. ISPRS Ann. Photogramm. Remote Sens. Spatial Inform. Sci. I-3, 147–152. Nelson, R., 2013. How did we get here? An early history of forestry lidar. Can. J. Remote. Sens. 39, S6–S17. Nie, S., Wang, C., Zeng, H.C., Xi, X.H., Xia, S.B., 2015. A revised terrain correction method for forest canopy height estimation using ICESat/GLAS data. ISPRS J. Photogramm. Remote Sens. 108, 183–190. Parrish, C.E., Jeong, I., Nowak, R.D., Smith, R.B., 2011. Empirical comparison of fullwaveform lidar algorithms: range extraction and discrimination performance. Photogramm. Eng. Remote Sens. 77 (8), 825–838. Parrish, C.E., Nowak, R.D., 2009. Improved approach to LIDAR airport obstruction surveying using full-waveform data. J. Surv. Eng.-ASCE 135 (2), 72–82. Pirotti, F., 2011. Analysis of full-waveform LiDAR data for forestry applications: a review of investigations and methods. Iforest—Biogeosciences Forestry 4, 100– 106. Qin, Y.C., Vu, T.T., Ban, Y.F., 2012. Toward an optimal algorithm for LiDAR waveform decomposition. IEEE Geosci. Remote Sens. Lett. 9 (3), 482–486. Roncat, A., Bergauer, G., Pfeifer, N., 2011. B-spline deconvolution for differential target cross-section determination in full-waveform laser scanning data. ISPRS J. Photogramm. Remote Sens. 66 (4), 418–428. 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. Photogramm. J. Finland 21 (1), 62–75. Słota, M., 2015. Full-waveform data for building roof step edge localization. ISPRS J. Photogramm. Remote Sens. 106, 129–144. Słota, M., 2014. Decomposition techniques for full-waveform airborne laser scanning data. Geomatics Environ. Eng. 8 (1), 61–74. Tian, J.Y., Wang, L., Li, X.J., 2015. Sub-footprint analysis to uncover tree height variation using ICESat/GLAS. Int. J. Appl. Earth Obs. Geoinf. 35, 284–293. Tseng, Y.H., Wang, C.K., Chu, H.J., Hung, Y.C., 2015. Waveform-based point cloud classification in land-cover identification. Int. J. Appl. Earth Obs. Geoinf. 34, 78– 88. 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 J. Photogramm. Remote Sens. 60 (2), 100–112. Wang, C., Li, Q., Liu, Y., Wu, G., Liu, P., Ding, X., 2015. A comparison of waveform processing algorithms for single-wavelength LiDAR bathymetry. ISPRS J. Photogramm. Remote Sens. 101, 22–35. Wang, Y.F., Zhang, J.Z., Roncat, A., Kunzer, C., Wagner, W., 2009. Regularizing method for the determination of the backscatter cross section in lidar data. J. Opt. Soc. Am. A 26 (5), 1071–1079. Wu, J.Y., van Aardt, J.A.N., Asner, G.P., 2011. A comparison of signal deconvolution algorithms based on small-footprint LiDAR waveform simulation. IEEE Trans. Geosci. Remote Sens. 49 (6), 2402–2414. Yi, D.H., Harbeck, J.P., Manizade, S.S., Kurtz, N.T., Studinger, M., Hofton, M., 2015. Arctic sea ice freeboard retrieval with waveform characteristics for NASA’s Airborne Topographic Mapper (ATM) and Land, Vegetation, and Ice Sensor (LVIS). IEEE Trans. Geosci. Remote Sens. 53 (3), 1403–1410. Zhuang, W., Mountrakis, G., Wiley, J.J., Beier, C.M., 2015. Estimation of aboveground forest biomass using metrics based on Gaussian decomposition of waveform lidar data. Int. J. Remote Sens. 36 (7), 1871–1889.