A SIGNAL DENOISING METHOD FOR FULL-WAVEFORM LIDAR DATA

The lack of noise reduction methods resistant to waveform distortion can hamper correct and accurate decomposition in the processing of full-waveform LiDAR data. This paper evaluates a time-domain method for smoothing and reducing the noise level in such data. The Savitzky-Golay (S-G) approach approximates and smooths data by taking advantage of fitting a polynomial of degree d, using local least-squares. As a consequence of the integration of this method with the Singular Value Decomposition (SVD) approach, and applying this filter on the singular vectors of the SVD, satisfactory denoising results can be obtained. The results of this SVD-based S-G approach have been evaluated using two different LiDAR datasets and also compared with those of other popular methods in terms of the degree of preservation of the moments of the signal and closeness to the noisy signal. The results indicate that the SVD-based S-G approach has superior performance in denoising full-waveform LiDAR data. * Corresponding author.


INTRODUCTION
A LiDAR system is an active remote sensing instrument that directly measures the distance between the target's surface and the sensor based on the round trip time of a laser pulse (Baltsavias 1999;Wehr and Lohr 1999;Jutzi and Stilla 2003;Wagner et al. 2006).Full-waveform Airborne Laser Scanning or LiDAR systems possess the advantages of continuous sampling to directly retrieve vertical and horizontal structural characteristics of both small and large objects, in addition to the 3D position of ground points.A variety of applications of this ranging technology, e.g. in urban planning, vegetation mapping and forestry, have already been explored with the aim of quantitative estimation of different parameters used in modelling the physical environment.
The range (distance between sensor and target), amplitude and pulse width are the most common attributes which can be derived from LiDAR waveforms by means of deconvolution/decomposition techniques (Wagner et al. 2006;Mallet and Bretar 2009;Heinzel and Koch 2011;Höfle et al. 2012).With full-waveform data, it is feasible 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 considerable interest in active remote sensing (Wagner et al. 2006;Wang et al. 2009).The accuracy of determining these parameters depends upon the shape of the waveform, and any change or shift in the waveform can cause uncertainty in the results.
On the other hand, an unavoidable phenomenon of images and signals, especially in remote sensing data, is noise that influences the resolution of outputs and restricts their capability by imposing distortion, misalignment and stretch.Preprocessing of signals in order to decrease the effect of noise is recommended before deconvolution of waveforms and extraction of further information from the signal (Pirotti 2010).Many methods have been proposed for both noise removal and improvement of the signal-to-noise ratio (SNR) in fullwaveform LiDAR data.Some methods are not appropriate in all cases and need specific assumptions.Typically, in Low-Pass filters, the high frequency segments are assumed to contain noise.However, some filters cause a shift in the signal and a change in shape of the signal (Hassanpour et al. 2012).
The FFT-based and Moving Average (MA) noise reduction methods, despite their ability to remove noise, can distort and smooth the waveform and this can affect the results of subsequent processing, such as the range values in the case of full-waveform LiDAR data.In addition, the MA filter can over smooth peaks and broaden the signals when peaks are narrower than the length of the filter, therefore reducing the amplitude of the signals.The signal accumulation approach (Kraftmakher 2006) assumes both noise with zero mean random characteristics and the signal phase to be known, although this is not the case in most situations and estimating the noise properties is not straightforward and feasible, as shown by Wu et al. (2012) and Okumura et al. (1993).The lasing and nonlasing echoes are necessary in the canonical correlation analysis proposed by Okumura et al. (1993), the nonlasing part, however, is not available for most end users (Wu et al. 2012).Some popular methods in the case of full-waveform data, such as the Wiener filter, cause oscillation and ringing effects near edges (Parrish 2007) as well as being unable to work efficiently in situations of low SNR (Hassanpour et al. 2012).
Existing approaches show varying success.The main purpose of the smoothing filters in full-waveform LiDAR is to reduce the noise level while maintaining the useful signal contents (Orfanidis 2010).This paper describes a method based on 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 fitting a polynomial to signal values such that raw (noisy) samples will be substituted by values obtained from fitted smooth polynomial curves.This method, based on the Savitzky-Golay (S-G) filter (Savitzky and Golay 1964), has the potential of preserving the shape and amplitude of signals as well as reducing the noise level (Hassanpour 2007;Orfanidis 2010;Schafer 2011;Hassanpour et al. 2012).However, this is achieved at the cost of potentially retaining much more noise than the MA smoothing method (Orfanidis 2010).
In Hassanpour (2007;2008), the S-G filter is applied in conjunction with singular vectors derived from a singular value decomposition (SVD) of a noisy Electroencephalography (EEG) signal in the time-frequency domain.The SVD-based noise reduction method is operated on a Hankel Matrix form of the noisy signal in the time domain by Hassanpour et al. (2012), followed by the S-G filter.A higher level of noise reduction can be achieved through eliminating the noise components during the SVD-based enhancement.
In this paper, we report our experience of denoising fullwaveform LiDAR data.The SVD-based S-G method has been investigated and evaluated against the Moving Average approach and the S-G filter.The S-G, and SVD-based S-G filters are first described and following this the results of these methods applied to full-waveform LiDAR data are presented and compared.Experimental results are then discussed and concluding remarks are offered.

Savitzky-Golay filter
As the first stage in describing the SVD-based S-G filter, the mathematical background of the S-G filter is briefly presented.The S-G filter output of length N and order d, applied to () xn takes the form of (Orfanidis 2010 The fitted polynomial ( ˆm x ) of order d can be written as where i c are the polynomial coefficients.
The polynomial basis vectors i s are defined as The filtered signal, x is given as According to Hassanpour et al. (2012), the following cost function should be minimized in order to adjust the polynomial degree and window length where The similarity between the enhanced and noisy signal and avoidance of large and abrupt differences between consecutive samples are taken into account in this cost function.

SVD-based S-G filter
As previously mentioned, the SVD-based S-G filter represents one of the extended approaches based on the S-G smoothing that can improve the level of signal enhancement.The method commences with a singular value decomposition of a matrix H, which is here the Hankel matrix of the noisy signal: where the variables with subscripts s and n represent the signal and noise subspaces.
The S-G filter will be applied on the U and V matrices of the signal subspace and the enhanced Hankel matrix e H , which represents the best estimate of the signal subspace, can be calculated from the enhanced U and V matrices: As can be seen in Figure 2, higher values of the cost function result from setting lower thresholds on singular values, while better signal smoothing, with lower oscillation and ringing effects, can be recovered as a consequence of an appropriate threshold value.However, increasing the threshold value to more than the optimal value leads to an unveiling of more superfluous detail and to similar results as those obtained via the S-G method.

EXPERIMENTAL RESULTS AND DISCUSSIONS
The S-G based filters have been implemented and applied to LiDAR datasets acquired over different terrain.Comparisons of their performance with the MA approach were conducted and results will be reported in this section, following a brief description of the dataset.

Datasets
Two full-waveform LiDAR data sets were employed in the experimental evaluation of the SVD-based S-G filtering approach.The first data set, recorded with an Optech ALTM 3100 incorporating a laser of 1065 nm wavelength, covered an urban area near Verona, Italy.The second, acquired with a Riegl LMS-Q560 sensor with a 1550 nm laser, was over a forested area near Chowilla, South Australia.Shown in Figure 3 are two representative samples of the received waveforms within the two data sets.

Comparison of different methods
The performance of the SVD-based S-G method applied to the recorded noisy waveform data has been compared with the MA smoothing approach and the S-G method.In the SVD-based S-G filter, small oscillations in the waveform can be moderated and diminished after utilizing the SVD, and noise will be removed using the S-G low-pass filter.Here, the smoothness factor is set to 0.7 to yield results with satisfactory degrees of smoothness, and the maximum polynomial degree (d) and window length (N) are set at 5 and 20, respectively.High values of the polynomial degree can result in retention of more of the higher frequencies, increasing cut-off frequency and hence allowing more noise to be preserved (Orfanidis 2010).The performance of the S-G low-pass filter is improved by exploiting the SVD of the signal and separating the noise and signal subspaces.As can be seen from Figure 4, although the value of the cost function for the S-G filter is clearly smaller than for the SVD-based S-G filter, the resulting signal from the former still remains noisy, with the filtered signal being only slightly smoother than the original waveform.It is noticeable that although these two methods can preserve the main moments of the waveform (the amplitude and the pulse width), the SVD-based S-G method is more robust in alleviating noise, and small abrupt fluctuations exhibited in the result of the S-G method can be smoothed more efficiently.
However, it is noteworthy that the smoothed waveform, resulted from the SVD-based S-G filter, is affected by fluctuation near to the end of the signal.This problem is related to the form of the Hankel matrix and can be reduced largely by substituting this matrix with a Toeplitz form of the noisy signal.The enhanced version of the SVD-based S-G filter, after replacing the Hankel matrix by a more generalized Toeplitz form, is depicted in Figure 5.
Figure 6 illustrates the difference between the MA and SVDbased S-G filtering results for a complex waveform with three distinct peaks and small distances between peaks.The SVDbased S-G filters demonstrate the ability to largely preserve the pulse width and signal amplitude.On the other hand, the amplitude of the waveform cannot be recovered properly using the MA approach, particularly in instances of sharper peaks (mostly from the ground returns).Thus, estimates of the signal amplitude can be significantly underestimated by the MA approach.In the case of the second peak, the amplitude is markedly influenced by the higher amplitude of the neighbouring peak, resulting in a significantly distorted form of the signal, which will cause incorrect interpretation of the data.Another approach to moderate oscillations yielded by the SVDbased SG filter on the Hankel form of the noisy signal is to weight the singular values.For this purpose, a normalized weight between 0 and 1 can be considered with evenly distributed values based on the total number of signal bins.This approach is yet to be implemented.
Finally, by considering a non-square Hankel matrix and reducing the number of its columns (and consequently its rank), it is possible to obtain more consistent results.Indeed, a smaller number of columns of the Hankel matrix can help to lessen the oscillation and more effectively remove noise by producing a smaller number of singular values.The influence of the number of columns is insignificant on the obtained results according to their cost function values for the modified thresholds on singular values (Table 1).As seen in Table 1, the threshold value is raised by increasing the number of columns of the Hankel matrix, allowing more singular values to contribute to denoising the signal and, accordingly, to increasing the cost function value.Although differences between these values are not considerable, it is recommended that the smallest number of columns be used such that the number is not smaller than the highest filter window length.Figure 7 illustrates the result of the related signal with 20 selected columns of the Hankel matrix.There is no evidence of oscillation and/or fake peaks at the tail-end of the signal, compared to the results from full-rank Hankel matrix in Figures 4-6.
No From the comparison between the SVD-based SG filter with a square and non-square Hankel matrix, it emerges that the result is a more robust preservation of the signal shape, with the cost function value being slightly smaller, as per Figure 8.It is noticeable that while the MA result leads to a smaller discrepancy value with the noisy waveform, it also markedly underestimates the amplitude of the first pulse in addition to causing a near distortion of the second pulse.lower rank matrix approximation of the Hankel form of the noisy signal, achieves better noise reduction and improvement in the results.This has the importance of reducing the computational cost by initially eliminating small singular values.The number of matrix columns then depends upon the maximum window length of the filter, while the threshold value does not play a significant role in signal distortion.

EnhFWF
corresponds to the enhanced waveform Noisy FWF corresponds to the noisy waveform [0,1]   denotes the smoothness factor

V
represent the enhanced left hand and right hand orthogonal matrices after applying the S-G filter s  = matrix of remained singular values.The same cost function applied in the S-G filtering is employed for the recovery of optimal parameters values.The key to the success of the SVD-based S-G filter is the separation of the singular values of matrix Σ into signal and noise subspaces since this can affect the resulting amplitude and the range of the filtered waveform.Different values of the threshold on singular values result in dissimilar shapes of the signal and misleading results can be obtained from inappropriate parameter values.The singular values shown in Figure. 1, obtained from the full-waveform signal in Figure.2, are normalized in the range of [0, 1] in order to easily choose the threshold value.A proper threshold value can be set on a point where the slope of the curve becomes more constant.

Figure 1 .Figure 2 .
Figure 1.Normalized singular values of the matrix Σ of a given noisy signal Figure.3 (a) is a waveform over flat, open terrain, whereas that in Figure.3(b) was collected over an area yielding more complex, multi-layer returns from both the vegetation canopy and understory, and the ground.

Figure 4 .Figure 5 .
Figure 4. Comparison of the SVD-based S-G and S-G filters

Figure 7 .
Figure 7. SVD-based S-G smoothing result with its specified number of columns of the Hankel matrix ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, VolumeII-5/W2, 2013  ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey NoisyH is divided into signal and noise subspaces by assigning a threshold to  and eliminating singular values of smaller than a specified magnitude by setting them to zero: