Automatic Estimation of Excavation Volume from Laser Mobile Mapping Data for Mountain Road Widening

Artículo patrocinado por Extraco, Misturas, Lógica, Enmacosa e Ingeniería InSitu, dentro del proyecto SITEGI, cofinanciado por el CDTI. (2012). 

Article sponsored by Extraco, Misturas, Lógica, Enmacosa and Ingeniería Insitu inside the SITEGI project, cofinanced by the CDTI. (2012)

Autores / Authors : Jinhu Wang 1,*, Higinio González-Jorge 2, Roderik Lindenbergh 1, Pedro Arias-Sánchez 2 and Massimo Menenti

1. Department of Geoscience and Remote Sensing, Delft University of Technology, Building 23, Stevinweg 1, Post Box 5048, 2628 CN Delft, The Netherlands; E-Mails: r.c.lindenbergh@tudelft.nl (R.L.); m.menenti@tudelft.nl (M.M.)
2. 2 Department of Natural Resources and Environmental Engineering, School of Mining Engineering, University of Vigo, E-36310 Vigo, Spain; E-Mails: higiniog@uvigo.es (H.G.-J.); parias@uvigo.es (P.A.-S.)

* Author to whom correspondence should be addressed; E-Mail: jinhu.wang@tudelft.nl;

Abstract:

Roads play an indispensable role as part of the infrastructure of society. In recent years, society has witnessed the rapid development of laser mobile mapping systems (LMMS) which, at high measurement rates, acquire dense and accurate point cloud data. 

This paper presents a way to automatically estimate the required excavation volume when widening a road from point cloud data acquired by an LMMS. Firstly, the input point cloud is down-sampled to a uniform grid and outliers are removed. For each of the resulting grid points, both on and off the road, the local surface normal and 2D slope are estimated.

Normals and slopes are consecutively used to separate road from off-road points which enables the estimation of the road centerline and road boundaries. In the final step, the left and right side of the road points are sliced in 1-m slices up to a distance of 4 m, perpendicular to the roadside. Determining and summing each sliced volume enables the estimation of the required excavation for a widening of the road on the left or on the right side. The procedure, including a quality analysis, is demonstrated on a stretch of a mountain road that is approximately 132 m long as sampled by a Lynx LMMS. The results in this particular case show that the required excavation volume on the left side is 8% more than that on the right side. In addition, the error in the results is assessed in two ways. First, by adding up estimated local errors, and second, by comparing results from two different

datasets sampling the same piece of road both acquired by the Lynx LMMS. Results of both approaches indicate that the error in the estimated volume is below 4%. The proposed method is relatively easy to implement and runs smoothly on a desktop PC. The whole workflow of the LMMS data acquisition and subsequent volume computation can be completed in one or two days and provides road engineers with much more detail than traditional single-point surveying methods such as Total Station or GPS profiling. A drawback is that an LMMS system can only sample what is within the view of the system from the road.

1. Introduction

There are various modern surveying and mapping techniques that make it possible to quickly obtain 3D surface geometry data. Notably, in recent years, there has been rapid development in Light Detection and Ranging (LiDAR) systems, for both airborne and mobile applications. In this study, a car-based system is used, which is generally referred to as a Laser Mobile Mapping Systems (LMMS). Such a system typically integrates one or more LiDAR profilers with a Global Navigation Satellite System (GNSS) and an odometer for positioning, and an Inertial Measurement Unit (IMU) for attitude control [1]. While driving, the LiDAR profilers determine distances between the car and its surroundings. By combining these distances with the position and attitude of the vehicle, a 3D point cloud is obtained that represents the area surrounding the road as visible from the car. Traditionally, the road and roadside geometry are measured by surveyors on a point-by-point basis using Total Station or GNSS. Unlike these methods, LMMS does not measure the 3D positions of well-defined points in the terrain. Rather, it acquires many random 3D points by sampling the whole road surface and the surface of the roadside, visible by the LMMS system. The quality of the individual 3D LMMS points is worse than the quality provided by Total Station or GNSS, but an LMMS is able to provide full coverage instead of acquiring only single points. In addition, the acquisition is fast and automatic, and there is no need for surveyors to leave the car. These latter advantages also distinguish LMMS from static Terrestrial Laser Scanning, in which a point cloud is obtained from a panoramic scanner mounted on a tripod [2–4]. Such panoramic scans need to be combined with GNSS positioning data in post-processing to acquire a georeferenced point cloud. To conclude, the LMMS is currently the fastest ground based method for acquiring 3D surface information in large areas. It has already been applied for high way surveying [5,6], sandy coast morphology and riverine erosion measurements [7,8], railway monitoring and rail geometry extraction [7,9], road environment management [8,10–12], and highway documentation [13–17]. Road management starts with the planning phase and ends with the rehabilitation or maintenance phase [7]. When the road has been constructed, documentation for road information is needed for an increasing number of other applications, such as noise modeling, road safety, road maintenance, location-based services and navigation [13]. Road documentation and management mainly consist of recording the road geometry and monitoring the road environment [14,18]. Road geometry refers to parameters used for the design of a road, such as design speed, stop sign distance, line of sight, number of lines, line width, longitudinal and transverse slope, road pavement materials and so on. Road environment refers to the surroundings of the road on both sides, including buildings, trees, vegetation, traffic signs, traffic light poles and other objects. The information discussed above plays an important role in the ongoing maintenance of the road, especially for mountain roads where there is a risk for rock and stone fall. Additionally, the geometric features observed on a mountainous road could be helpful for monitoring the flow of rainwater in the case of heavy rain and may be used to assist natural disaster prevention [19]. Moreover, steep and unstable road sides may cause landslides, resulting in further road damage [20]. This paper presents an approach for automatically estimating the roadside material volume to be excavated for road widening. Firstly, LMMS point cloud data in a mountainous area is down sampled and the outliers and noisy points are removed. Based on the data, normal vectors, and 2D slopes are estimated at every point. Then, an automatic iterative floating window approach, taking advantage of point height, normal vector and slope, is used to filter and segment the road points. After that, a local neighbourhood feature is defined based on the vectors between a query point and its neighbouring points to obtain the outline and skeleton of the road. These steps finally allow us to compute the volume of the roadside material that would have to be moved to widen the road by 4 m. In addition, an analysis of the quality of the results is presented notably by comparing results from different datasets both sampling the region of interest. Finally, conclusions and future work are discussed.

2. Methodology

This section explains the procedure for estimating the amount of roadside material that has to be removed for road widening. It consists of five steps. (1) Point cloud data preprocessing; the original point cloud data have a very high point density and need to be downsampled before processing. Additionally, outliers are removed, and a Digital Surface Model (DSM) is generated; (2) Surface normal estimation; (3) Slope and aspect estimation; (4) Road detection and segmentation; taking the normal and slopes as estimated in steps (2) and (3) as input, an automatic iterative filtering approach is used to segment the road points from the downsampled point cloud data; (5) Volume computation; the material volume that needs to be moved to widen the road is computed based on step (4). A flowchart of the procedure is shown in Figure 1.

2.1. Point Cloud Data Preprocessing

The raw point cloud data have a very high point density, so for quick and efficient processing, downsampling is a necessity. The approach followed here is to represent the point cloud data by voxels [21–23]. In this paper, a uniform-sized voxel filtering approach was used, as implemented in the Point Cloud Library (PCL) [22]. The main concept of the approach is to create a 3D voxel using a given voxel size. All points in a voxel are represented by their centroid [13]. The outliers of the downsampled point cloud data have to be removed before estimating geometric features from the point cloud data. In this procedure, the query point neighbourhood concept is introduced, as shown in Figure 2. Here, P represents the set of points within radius Rquery of query point pquery, so P is defined as:

where p(i) represents the individual points from the point cloud.

Figure 1. Overall methodology of data processing.

Figure 2. Neighbourhood of a query point within a certain radius.

Outliers are the measurements located at jump edges or discontinuous boundaries where there should be no points [24]. There is abundant literature on the approaches used to remove outliers [25–27]. The approach used in this paper is based on a statistical analysis of each point’s neighbourhood [28,29]. For each point p_query Î P calculated. After that, for each point in the point cloud, the mean distance and standard deviation of the distances to their k nearest neighbours are determined. The main objective is to retain only those points whose mean distance to the closest k neighbors is similar to the average mean distance. As this describes a measure of the underlying point cloud density surrounding a point, the criterion for keeping a point is simply formulated as follows:

where, α is a desired density restrictiveness factor, while 􀟤􀯞 and 􀟪􀯞 are the mean and standard deviation of the distance from a query point to its neighbors respectively, and 􀜲􀗛 is the set of remaining points.

2.2. Local Surface Normal Estimation

The surface normal at a discrete point is a vector perpendicular to the tangential plane of the local surface at that point. Various methods exist to estimate a normal at a certain point in 3D point cloud data [29–33]. The simplest is based on 3D plane fitting. With this method, the problem of determining the normal of a point becomes a least-square 3D plane estimation problem for a suitable spatial neighbourhood. Figure 3 depicts the concept of normal estimation in discrete 3D point cloud data. 

When estimating a local normal using least-squares, the goal is to approximate the tangent plane at a point of interest and take that plane’s normal. The best fitting plane is obtained by determining those planar parameters that minimize the squared distances between suitable neighbouring points and the plane to be estimated. Suppose we have a point of interest with Euclidean coordinates (x, y, z)^T set of k neighbouring points. The least-squares method yields a normal vector (n_x, n_y, n_z)^T which the error is minimized.

Figure 3. Normal estimation at a query point using a local plane fitting approach

2.3. Local Slope Computation

The 2D slope, also known as the 2D gradient, is the vector field of a surface. The vector direction points in the direction of the greatest change in height, and the vector’s magnitude is equal to the rate of change. Based on the previously generated gridded DSM, the height at every grid point is interpolated from neighbouring points by inverse distance interpolation with power 2. If the grid size is d(grid), then the slope Si in the direction of each of the eight neighboring grid cells is given by: 

Here, H(i) is the elevation of the i-th neighbour of the query point, and hquery is the elevation of the query point itself. The variable di is the distance between the i-th neighbour and the query point. Note that di is the square root of the grid size in diagonal direction.

2.4. Road Detection

Based on the normal vectors and slopes that were computed at each point in the previous steps, an

automatic iterative point cloud data filtering approach is used to detect road surface points from the

point cloud data. The main steps are as follows:

• Input an initial grid and window size;
• Generate a virtual reference 3D gridded layer. The elevation of each grid point is interpolated from its neighbouring points, and also the grid point’s normal as unit vector and its direction to zenith. Based on this layer, a window of predefined size is created to move over the grid and point cloud;
• In the current moving window, compare for each grid point the 3D layer grid elevation and the normal vector direction with that of the point cloud; Calculate the height and angle differences between the 3D layer and the point cloud to verify if the differences are beyond the threshold;
• If the difference is less than the distance and direction threshold, then the point is accepted as a road point, else the point is regarded as off-road point;
• Go to step (2) and generate a new virtual layer using a smaller grid size, then iteratively process the point cloud again;
• The loop ends when the grid size reaches a pre-defined smallest size.

2.5. Calculation of Excavation Volume

In this step, based on the segmented road points from step 4, the road’s outlines are determined.

First, a local feature descriptor is defined based on the neighbourhood concept.

Figure 4. Illustration of local neighbourhood feature descriptor.

After the road outline is determined, the road central line is estimated based on the location of the

road edges. Now suppose the road needs to be widened to four lanes. As a consequence, a certain

volume of the road has to be removed or added to extend the flat road surface.

As shown in Figure 5, roadsides were divided into slices. For each slice, the volume is computed.

Summing all of the sliced volumes together gives the total volume that needs to be excavated or filled.

Note that the sign of a sliced volume indicates whether material needs to be removed (positive sign) or

added (negative sign). Based on the road central line and edge lines, as well as the expanded width of

the road, the expanded edge points are found.

Figure 5. Computation of excavation volume.

Figure 6. Slice volume computational geometry.

3. Implementation and Testing of the Method
3.1. Software Implementation

Continue on: http://carreteras-laser-escaner.blogspot.com/2014/11/automatic-estimation-of-excavation_33.html