SPARSITY BASED REGULARIZATION APPROACHES IN RECONSTRUCTING THE RANGE AND CROSS SECTION IN FULL-WAVEFORM LIDAR

Accurate range determination and retrieval of the cross section are two important issues in the processing of full-waveform LiDAR data, especially between closely located targets. The dependency of the received waveform on the emitted pulse can be removed through deconvolution and consequently comparisons between waveforms recorded by different sensors become feasible and meaningful, since the cross-section is independent of system specifications. Common methods, such as the Wiener filter, are reported for producing oscillation of the results around main peaks, along with negative values of the amplitude. Regularization is necessary to approximate a stable solution of deconvolution resistant to noise or error. A sparse solution of these linear inverse problems can be attained by minimizing the one-norm of the solution. Satisfactory deconvolution results can then be achieved by utilizing sparsity constraints. The results of regularization methods with sparse solutions have been evaluated using two synthetic datasets and distinctions are highlighted for comparison with those from two widely used deconvolution methods, in terms of the level of surface response retrieval. The results presented indicate that the one-norm regularization approach can outperform the other methods considered. * Corresponding author.


INTRODUCTION
The round-trip time of a laser pulse between the sensor and the target surface is the basis for direct range measurement in pulsed airborne laser scanning systems, which are commonly referred to as LiDAR (Light Detection and Ranging) systems (Baltsavias 1999;Wehr and Lohr 1999;Jutzi and Stilla 2003;Wagner et al. 2006).Single echoes are recorded in real-time in the case of discrete return systems, whereas they can also be detected in a post-processing step in full-waveform systems by end users.Therefore, potentially higher accuracies of range extraction and even detection of more echoes (weak and/or overlapping echoes) are possible from full-waveform data, depending upon the adopted processing methods (Höfle et al. 2012).In addition, full-waveform data allow to derive the target cross-section, which is not possible in the case of discrete return LiDAR (Höfle and Pfeifer 2007) and which is an issue of great interest in active remote sensing (Wagner et al. 2006;Wang et al. 2009).
Decomposition methods for extracting a parametric description of the pulse properties; e.g. the range, pulse width and amplitude, have been proposed by several authors for different sensors (Hofton et al. 2000;Persson et al. 2005;Wagner et al. 2006).Generally, in the case of Gaussian-based decomposition methods, it is assumed that a set of Gaussian functions can characterize the probability density function of the target response.
The Gaussian decomposition method has been used in many cases to interpret targets related to the backscattered waveform.
Results depict the challenging nature of this method in the case of echoes of low strength level (low SNR).In addition, Gaussian decomposition is not capable of calculating the crosssection in some situations (complex waveforms) due to the complexity of the land surface, deviation of intensity within the footprint of the laser beam and also the necessity of determining the number of targets beforehand (Wang et al. 2009;Lin et al. 2010;Mallet et al. 2010;Zhu et al. 2010).In some cases, negative amplitude values might be obtained or a solution might not even be achieved in the fitting process (Lin et al. 2010;Roncat et al. 2011).Furthermore, The hypothesis of vertical Gaussian distribution of targets within the laser beam is not valid in the case of small-footprint LiDAR data and might result in inaccurate results (Chauve et al. 2007).These problems can influence not only the accuracy of the range, but also the point cloud coordinates and further processes.
More complicated distribution functions (e.g. the Generalized Gaussian and Lognormal) have been utilized by Chauve et al. (2007) and Chauve et al. (2009) in order to improve the density and geometric accuracy of the extracted points.These functions also have the capability of extracting more attributes from the shape of the waveforms.The generalized Gaussian model has been shown to yield a better than expected quality of fitting globally, yet the Lognormal model yields a better result only in the case of asymmetric pulses.
Weak and overlapping signals were studied by Lin et al. (2010) in order to detect targets using a rigorous Gaussian decomposition.Their proposed method had the capability to ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey both extract points in overlapping signals with higher accuracy than the methods in commercial software (e.g. the centre of gravity and Gaussian pulse fitting), and to detect more weak amplitude pulses.However, the high computational cost potentially restricts applicability of the method.Mallet et al. (2010) improved the curve fitting by considering a set of parametric pulse shapes.Results revealed shortcomings of the approach in the selection of an appropriate model in the case of small-footprint data.Utilization of derived shape features individually did not result in superior discrimination, but along with geometric features it tended to achieve better results.
Similar received waveforms can be recorded after the convolution of different target cross-sections with the system waveform, and as a consequence, indiscernible results may be achieved (Neilsen 2011).Therefore, performing deconvolution to remove the sensor effects can help to retrieve the target characteristics.In addition, the target cross-section can be useful in differentiation between targets in different sensor operations since it is independent of the instrument and is related to roughness and geophysical properties of the surface (Roncat et al. 2011).
Allocation of a constant and limited time for detection of the received signal causes a reduction in resolution and leads to a mixed and overlapped waveform for surface objects with a short distance between them (Wu et al. 2011), in addition to the fact that the received waveform is a convolution of the target crosssection and the emitted waveform (Wagner et al. 2006;Zhu et al. 2010).Hence, deconvolution is essential to both removal of the dependency of the received signal from other parameters and to the reconstruction of the target characteristics (Jutzi and Stilla 2006;Shan and Toth 2009;Zhu et al. 2010).Deconvolution approach has the advantage that knowing the number of peaks (as in Gaussian decomposition) is not necessary (Zhu et al. 2010).Also, assumptions on pulse shape are not essential (Mallet and Bretar 2009).
Derivation of the target cross-section from the received waveform is an ill-posed problem and cannot be solved uniquely without enforcing some constraints and restrictions (Wang et al. 2009).Different deconvolution methods exist.Among them, the Wiener filter has been extensively used to restore the target cross-section and enhance the range estimation from full-waveform LiDAR (Jutzi and Stilla 2006;Parrish 2007;Wu et al. 2011).Jutzi and Stilla (2006) differentiated between compound pulses from different target parts by modelling essential components of the sensor and surface to improve estimation of the range.It was recommended not to use the maximum amplitude of the pulse directly in order to extract the range, since this value is altered after superimposing different target responses (Jutzi and Stilla 2006).The approach was utilized by Parrish (2007), with the aim of detecting vertical barriers in an airport area, and an increase in vertical target detection of over 15% was achieved.
Although the Wiener filter is computationally efficient, the required degradation level of the signal cannot be achieved due to the problems of (1) properly characterizing the signals in the frequency domain, (2) emerging negative values, (3) producing ringing effects near edges (Parrish 2007;Wu et al. 2011) and (4) the possibility that noise might be converted to another form of noise in the case of low SNR (Hassanpour et al. 2012).
A regularizing approach to restrain the noise propagation and also improve the computational cost in deconvolution of fullwaveform data was considered by Wang et al. (2009).Some small fluctuations were seen close to the main peaks and the authors believed that noise was the main cause.The approximation of the noise rate in their reported experiment was carried out based on the assumption that the data had was of good quality.However, new methods must be taken into account when the level of noise in unknown, which is always the case regarding real data.
The estimated quantities from existing methods, only based on the received waveform, cannot correspond to the surface properties; hence they are not trustworthy from sensor to sensor.B-Splines, as a linear deconvolution approach, were implemented in Roncat et al. (2011) to compute the cross section, range and other target characteristics.In comparison with Gaussian decomposition, symmetric emitted and received waveforms were not considered as the only basis in deconvolution based on B-splines.However, Roncat et al. (2011) indicated small negative values in the calculation of the cross-section together with echoes with narrower width than the emitted waveform.Zhu et al. (2010) explored the efficacy of the Gold's deconvolution method on retrieving surface characteristics.This method tends to achieve different results by setting a different number of iterations, while the solution is always non-negative.
The optimal number of iterations must be found in this case.However, some detected peaks were in between the true data peaks, which can cause significant errors in range detection and false target detection.Features extracted from waveforms via an approach based on PCA were shown to be related to the amount of biomass in the study area.
In Wu et al. (2012), four steps were introduced as a preprocessing sequence to calibrate Full-Waveform LiDAR data by using simulated and real data, mainly for biomass modelling.Among these steps, deconvolution emerged as the most useful step in nadir scan angle, whereas in the off-nadir situation, the geometric calibration of waveforms showed the most enhanced results.
By using robust signal processing approaches, it is possible both to decrease the side effects of the parameters that affect the received waveforms, and to reconstruct calibrated signals for further processing (Shan and Toth 2009;Wu et al. 2012).Furthermore, attributes based on the shape of the waveform which are ignored by discrete returns could be useful in the modelling of objects characteristics, such as biomass (Muss et al. 2013).
A comparison of the main methods for range detection from full-waveform LiDAR data (namely Gaussian decomposition, deconvolution and a hybrid method of deconvolution and decomposition) using laboratory collected data has been ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey reported by Parrish et al. (2011).The deconvolution and hybrid (the Wiener filter, followed by the Gaussian decomposition) approaches proved to be capable of separating close targets, while the computational cost is the smallest for Gaussian decomposition.It was thus recommended that waveform processing with a fusion of several approaches may yield the optimum possible results.
Regularization methods with the advantage of sparse solutions have been emerged as better choices to work on usual images.
In this paper, we represent comparisons between the Wiener and R-L deconvolution methods on the one hand, and the sparsity based regularization approaches on the other hand, in order to optimally reconstruct the surface response (cross section).The simplified form of full-waveform signal convolution is first presented, following by 2-norm and onenorm Tikhonov-type (variational) regularizations.Experimental results of these methods applied to full-waveform LiDAR data are then presented and discussed, followed by concluding remarks.

METHODOLOGY
Generally, the received waveform, () gt , can be explained with some simplifications as a convolution of the emitted signal, () ht , and the target response (cross section), () ft, in the form of where () nt is the additive noise element.The matrix form of (1) can be expressed as where H is a two-dimensional matrix of the emitted waveform.
The ill-posed nature of the target cross section restoration from the degraded waveform is due to a breakdown occurring in the existence, singularity or immutability of the result (Gunturk and Li 2012).Regularization methods are a well recognized way to solve this problem so as to yield a solution (Wang et al. 2009;Gunturk and Li 2012).Regularization is necessary to solve the deconvolution since by applying typical approaches, substantial fluctuations can appear as a result of minor noise and error (Bandkh et al. 1997;Morhac 2006;Wang et al. 2009;Morháč and Matoušek 2011;Neilsen 2011).Direct and iterative deconvolution solutions have been defined.Tikhonov regularization (Tikhonov and Arsenin 1977;Tikhonov et al. 1995) is the most well-known method of the direct type, where the following cost function is minimized where 0   is a regularization parameter, and L is a regularization operator (an identity matrix or a high-pass filter).
The results of applying Tikhonov regularization in cases where the assumption of a Gaussian noise distribution is not valid may not be acceptable (Pan 2010).This regularization technique may produce negative values together with oscillation near edges (Morháč and Matoušek 2011) and it does not result in preferred results in some situations (Pan 2010).
In the case of usual images, it has been shown that using the one-norm instead of the two-norm in regularization functions is more appropriate, with the potential of achieving a sparse solution (Gunturk and Li 2012).In this case, the problem can be expressed as || .|| indicates the Euclidean norm, and Elements of f with small values are enforced to zero as a consequence of the one-norm (Figueiredo et al. 2007).
Two equivalent convex constrained forms of the unconstrained convex problem in (4) are as following: and where ,0  are real parameters.

EXPERIMENTAL RESULTS AND DISCUSSIONS
Sparse solutions of the target cross section retrieval for two synthetic datasets have been obtained based on BPDN and LASSO, and these are compared with solutions from two extensively used methods, namely the Wiener filter and R-L deconvolution.

Datasets
Two synthetic LiDAR pulses have been generated with different characteristics, for both the emitted pulse and the cross section.
The emitted pulse and these two cross sections are depicted in Figures 1-3.As can be seen in Figure 2, the synthetic cross ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey section contains three separate segments (returns) with enough range difference between them.However, four segments of the synthetic cross section in Figure 3 are closer together, which makes the retrieval task more problematic since the convolution result of this signal with the emitted pulse contains overlapping pulses.The received waveforms, generated via (1), are obtained from convolution of the transmitted and cross section samples.
Additive Gaussian noise has also been added to the synthetic received waveforms, in this case white Gaussian noise with a SNR of 2 dB.The received waveforms and the noisy versions are illustrated in Figures 4 and 5.The third pulse in Figure 5 is generated as a contribution of two last pulses from the cross section sample 2 in convolution with the transmitted pulse, while other two pulses can still be separated more easily because of the longer range between them.

Experimental results
For two different datasets, the results two sparse solutions, together with the Wiener filtering and R-L deconvolution outputs are represented in Figures 6-9.The BPDN and LASSO results with their proper parameters for dataset 1, which are depicted in Figure 6, are almost the same without any oscillations or under/ overestimation of the pulse amplitude.On the other hand, in comparison to R-L deconvolution, the Wiener filtering produces a result with more of the oscillations and fake peaks that can lead to misinterpretation in future steps.Specification of an appropriate number of iterations in R-L deconvolution can be an issue, and some false peaks might be created as a consequence of determining a wrong number of iterations.For the closer peaks of dataset 2, LASSO resulted in a slightly better result in term of the under/overestimation of the signal amplitude, while both LASSO and BPDN could successfully retrieve the cross section for a case in which the received waveform constrained overlapping peaks.The results of the Wiener and R-L methods, however, illustrate an imperfection in retrieval of the cross section, especially in the case of the Wiener filter, both in terms of oscillation and under/ overestimation of the signal amplitude.
Wu et al. (2011) compared three deconvolution approaches, the Wiener filter, Richardson-Lucy (R-L) and nonnegative least squares (NNLS) to recover the spatial resolution of fullwaveform signals with regard to the target cross-section retrieval and biomass level classification accuracy.Results indicated that the R-L deconvolution was superior in both quantitative and qualitative terms, whereas previous studies(Harsdorf and Reuter 2000;Jutzi and Stilla 2006;Nordin 2006;Roncat et al. 2011) did not assess results quantitatively.

Figure 6 .
Figure 6.Reconstructed cross sections from BPDN and LASSO for data sample 1

Figure 8 .
Figure 8. Reconstructed cross sections from BPDN and LASSO for data sample 2