Surface Fitting Filtering of LiDAR Point Cloud with Waveform Information

Full-waveform LiDAR is an active technology of photogrammetry and remote sensing. It provides more detailed information about objects along the path of a laser pulse than discrete-return topographic LiDAR. The point cloud and waveform information with high quality can be obtained by waveform decomposition, which could make contributions to accurate filtering. The surface fitting filtering method with waveform information is proposed to present such advantage. Firstly, discrete point cloud and waveform parameters are resolved by global convergent Levenberg Marquardt decomposition. Secondly, the ground seed points are selected, of which the abnormal ones are detected by waveform parameters and robust estimation. Thirdly, the terrain surface is fitted and the height difference threshold is determined in consideration of window size and mean square error. Finally, the points are classified gradually with the rising of window size. The filtering process is finished until window size is larger than threshold. The waveform data in urban, farmland and mountain areas from “WATER (Watershed Allied Telemetry Experimental Research)” are selected for experiments. Results prove that compared with traditional method, the accuracy of point cloud filtering is further improved and the proposed method has highly practical value. * Corresponding author Funded by State Key Laboratory of Geo-information Engineering,NO. SKLGIE2016-M-3-1


INTRODUCTION
In recent years, with the continuously emerge of commercial small-footprint airborne ground measuring device, a new airborne full-waveform LiDAR measurement system receives extensive attentions, which is capable of recording laser pulse (Lai Xudong, 2013;Qin Yuchu, 2011).Compared with the traditional airborne LiDAR system, it mainly provides the users no longer discrete point cloud, but many sets of onedimensional sampling pulse waveform data.Processing the waveform data not only can obtain the discrete point cloud, but also return information of targets (also known as waveform information) (Mallet and Bretar, 2009).Full-waveform LiDAR technology not only has a significant role on vegetation and biomass information, but greatly helps laser ground measurement, of which the most typical application is DEM generation.Using full-waveform LiDAR data to generate high quality DEM includes two key steps: waveform data decomposition and point cloud filtering.Current waveform data decomposition algorithms mainly include nonlinear least square method (Hofton M, 2000;Chanve A, 2007), expectation maximization (Persson A, 2005;Li Qi, 2008) and markov chain monte carlo algorithm (Hernandez-Martin S, 2007), in which the nonlinear least square method based on Levenberg Marquardt (LM) is widely applied.This algorithm both has the advantages of gradient method and newton method, but has high requirements for initial value, which can easily run into local optimum; Filtering algorithms based on discrete point cloud geometric features can be divided into three main categories (Huang Xianfeng, 2009): morphological method (Qi C, 2007;Sui Lichun, 2010), method based on interpolation (Axelsson P, 2000;Kraus K, 1998) and method based on surface constraint (Su Wei, 2009).In surface constraint method, because of the limitation of geometry information, there may be some error seed points, which bring adverse impacts on surface fitting filtering; the feasibility of differentiating ground and nonground by using waveform parameter is verified by some scholars (Wagner et al., 2008;Doneus M, 2006;Mandlburger G, 2007;Mücke et al, 2010;Jalobeanu and Goncalves, 2014), however, such algorithm remains to be further researched.This paper presents a surface fitting filtering method with waveform information.The basic idea of the method is to use decomposed waveform information and robust estimation theory to help the selection of seed points; and the surface fitting weights are determined by waveform parameters; the window size and mean square error of surface fitting are comprehensively considered to make height difference threshold determined adaptively.The workflow is shown in Figure 1.

WAVEFORM DECOMPOSITION
Waveform data is usually composed of two parts, one is waveform attribute information, including transmitted pulse number, laser wavelength, transmitted pulse width, and threedimensional coordinates of sampling points at the initial and end; another is all the pulse sampling data.Waveform is decomposed into multiple echoes to describe characters of different objects under pulse, in which waveform simulation and waveform modelling are two key technologies.Generally, waveform data is regarded as a combination of some Gaussian component (Wagner W, 2006), so it is modelled by Gaussian function model; waveform simulation is to solve optimal parameters.As mentioned above, the LM algorithm is usually adopted.However, the solutions of LM algorithm are usually local optimal, which make waveform decomposition results inaccurate.In this paper, global convergent LM is introduced to solve this problem and its decomposing steps are as follows: Step 1: pre-processing.In the process of laser pulse scanning, it may generate some background noise on account of the influence of weather or other factors, which appear small amplitude shaking form in the waveform.Waveform may be decomposed incorrectly if such noise is not removed.
Step 2: peak detection.Calculate second derivative of each sampling point.Regard zero position as the initial peak position.The number of detected peaks determines the number of waveform parameters.
Step 3: parameter initialization.As for Gaussian function model, parameters to be initialized include pulse amplitude k A , half width k  and pulse distance k  .
Step 4: parameter optimization.Using global convergent LM can obtain global optimal waveform parameter solutions.The principle of this algorithm can be seen in reference.
Step 5: calculate residual.The residual  is calculated, if residual is less than former, repeat Step2 to Step5.The method continues until residual does not reduce anymore.The formulation of residual is shown in formula(1), where n is the number of sampling points, m is the number of parameters, y is the amplitude of each sampling point in original waveform, () fx is the amplitude of each sampling point in the fitted waveform.
Step 6: generate three-dimensional point cloud.After the iteration, the waveform parameters are obtained.Each waveform component corresponds to an object reflection echo, so the coordinates of the corresponding object are calculated based on waveform parameters.

FILTERING USING WAVEFORM INFORMATION
Point cloud filtering refers to processing and recognizing the irregular distributed discrete points on different objects, namely distinguishing ground point cloud (bare area and road) and nonground point cloud (buildings, cars and vegetation etc.) (Sui Lichun, 2010).The traditional filtering methods generally use geometric features of discrete point cloud, which cannot well recognize low vegetation.Waveform parameters obtained by waveform decomposition will be in favour of improving filtering effect.As mentioned, scholars conducted some attempts by utilizing waveform information.In mountainous area, they verified the feasibility of using the pulse half width to distinguish ground and vegetation.However, the filtering method fully using waveform information remains to be further explored, the experiments in cities and other areas remain to be further conducted.Traditional surface fitting filtering method (Su Wei, 2009) regarded the lowest point in window as ground seed point.Block terrain is fitted by seed points in several windows.However, some points on buildings will be selected improperly when it contains large size building in window, so that the fitted terrain should be inaccurate.In order to solve this problem, waveform parameters and robust estimation theory are used to select reliable ground seed points.Weighted surface fitting based on waveform parameters and adaptive determination of height difference threshold should make up the deficiencies.The workflow is shown in Figure 2. It is also an iterative process.With increasing of the window size, the non-ground parts will be filtered gradually.The key technologies contain detection of abnormal seed points, weighted surface fitting and adaptive determination of height difference threshold.

Detection of abnormal seed points
Pulse half width is helpful to distinguish vegetation and building edge (Doneus M, 2006).Take Riegl LMS-Q560 for example, it is probable that the point with pulse half width lager than 1.7ns belongs to vegetation or building edges; and belongs to flat ground with pulse half width smaller than that.The pulse half width threshold  is determined by the pulse width  , which is given in system specification parameters.And the formula is 2 2ln 2   . Therefore, the seed points are first screened using pulse half width.A further detection is conducted based on robust estimation theory.The initial value of robust iteration should make sure high breakdown pollution rate.So, the trend surface model is calculated using robust estimation with high breakdown pollution rate (Yang Yuanxi, 1996).And the equivalent weight of each seed point after robust iteration is used to detect abnormal.The specific steps are as follows (Chang Yifeng, 2011): Step 1: calculate the strong eliminating function as follows: Where, 0 i p is the function value of the th i seed point; the residual i v is height difference between the th i seed point and the median of all seed points, which is ; 0 c is critical value of weight function, and it is generally set as 1.0 or 1.5.
Step 2: calculate initial value of trend surface model using 0 i p .Given observation equation (3) Where, iterative equivalent weight matrix is is equivalent weight value of seed point, which is determined by IGGIII weight function: Where, ( ) ( ) ( ) ; the range of 0 k and 1 k are 1.0~2.0 and 2.5~4.0 separately.Step 3: ( 1) ( ) max kk xx    is the iteration stopping criterion.If the equivalent weight value of seed point is 0, it is labeled as abnormal one after iteration stopped.
During iteration, the weight values of untrusted seed points are gradually decreased by judging standard residual, so that, the accuracy of detecting abnormal seed points is increased.

Weighted surface fitting
After detection of abnormal seed points, the remained seed points are used to fit the terrain.The weight value of seed point is calculated according to pulse half width to obtain a better fitting result.The pulse half width get from waveform decomposition reflect the flatness of terrain.The greater  is, the smaller the possibility that this point is on a flat surface, otherwise the possibility is greater.Therefore, the weights of seed points are set according to formula (6).
Where, j p is the weight of the th j seed point in the block, its pulse half width is j  , n is the number of seed points in the block.The significance of setting weight in this way is that magnifying the contribution of points which have small pulse half width and decreasing the influence of the points which have large pulse half width.

Adaptive height difference threshold
The height difference between discrete point and the fitting surface in the block can finally determine its attribute.When the height difference is greater than threshold T , this point will be marked as non-ground point, or it will be marked as ground point.However, different blocks have different shape of terrain.
It is obviously improper that making the height differences threshold constant, so after overall considering the influence of window size and surface fitting error, we propose a method to adaptively determine height difference threshold as shown in formula (7).Its basic idea is that if window size and mean square error are both small, the fitting result will present the real terrain, namely the fitting precision is high.Therefore, the height difference threshold can be set strictly to ensure the low object points can be filtered effectively.On the contrary, as the window size and mean square error increase gradually, the fitted surface may be untrusted.Therefore, height difference threshold should be set loosely to avoid filtering error.P are both set as 0.5; if height difference threshold is too small, some ground points will be hard to remain; if the threshold is too large, several non-ground points could not be filtered.Therefore, a range from 0.3 to 2 is given to the height difference threshold, namely, min T is 0.3, and max T is 2.

EXPERIMENTS AND ANALYSIS
Small footprint LiDAR waveform data from "WATER (Watershed Allied Telemetry Experimental Research)" are selected for experiments.These data are collected by LMS-Q560 airborne full-waveform LiDAR system, which is developed by RIEGL.The beam divergence is 0.5mrad, wavelength is 1550nm and pulse width is 4ns.We select data of urban area and crop land from the surveying area in Zhangye city, and data of mountain land from surveying area in Dayekou forest.The information of data is shown in Table 1.1~Figure 3 shows the filtering results of urban area, crop land and mountain land respectively, in which, (a) is the original point cloud, (b) is the filtered point cloud.Furthermore, we filtered these data sets by traditional method, adaptive TIN method and our method, and I class error rate, II class error rate and total error rate are counted separately.I class error rate is the rate that ground points are classified as non-ground points; II class error rate is the rate that non-ground points are classified as ground points; total error rate is the weighted sum of I and II class error rate.The adaptation is represented by I and II class error rate, and the feasibility is reflected by total error rate (Huang Xianfeng, 2009).Reference filtering results are edited by manual.The maximum window size should be determined according to the objects in test area no matter using surface fitting method or adaptive TIN method.In our experiments, the maximum window size is set as 90m in urban area, 50m in crop land and 30m in mountain land.For different methods, we choose the same window size in one test data to ensure fair.The error rates of experimental results are shown in Table 2~Table 4     Those 2 test data are both characterized by flat terrain.The wrong seed points on buildings and dense vegetation make the fitted surface unmatched with real terrain in traditional method, so that large amount of ground points are filtered and some nonground points are remained, therefore, I class error rate and II class error rate are both relatively high.Our method makes up the deficiencies by detection of abnormal seed points.In mountain land data, the height difference threshold is relative large as the terrain relief, therefore, it should remain many nonground points, and its small I class error rate is just a false appearance.Our method improves this problem by integrating waveform information and self-adaptive height difference threshold.And comparison with adaptive TIN method shows that, the result in urban area by our method is better, and results in crop land and mountain land are nearly the same with adaptive TIN method.

CONCLUSIONS
Full-waveform LiDAR is an active technology in the field of photogrammetry and remote sensing.It is characterized by full record of return pulse and the information obtained by it could be in favor of DEM generation.Nowadays, filtering using waveform information is not researched in depth, and there are shortcomings in selection of seed points in surface fitting method.We propose a weighted surface fitting filtering method using waveform information.Compared with traditional method, it has more accurate filtering result, and the problem of improper selection of seed points and deviation of surface fitting could be resolved.Waveform information could be not only used to DEM generation, but also be used for point cloud classification.Deep mining of waveform characteristics could improve the classification accuracy further, which should be our next research emphasis.

Figure 2 .
Figure 2. Flow chart of filtering . (a) original point cloud (b) filtered point cloud Figure 3. Filtering result of urban area (a) original point cloud (b) filtered point cloud Figure 4. Filtering result of crop land (a) original point cloud (b) filtered point cloud Figure 5. Filtering result of mountain land

Table 1 .
Information of test data Figure

Table 2 .
Filtering results of urban area

Table 3 .
Filtering results of crop land

Table 4 .
Filtering results of mountain landThe experimental results show that, our method realizes point cloud filtering effectively as the non-ground parts are filtered in all the above three test areas.I class error rate, II class error rate and total error rate in urban area and crop land are obviously improved by our method compared with traditional method.