Evaluation of a LIDAR Land-Based Mobile Mapping System for Monitoring Sandy Coasts

Abstract

The Dutch coast is characterized by sandy beaches flanked by dunes. Understanding the morphology of the coast is essential for defense against flooding of the hinterland. Because most dramatic changes of the beach and the first dune row happen during storms, it is important to assess the state of the coast immediately after a storm. This is expensive and difficult to organize with Airborne Laser Scanning (ALS). Therefore, the performance of a Land-based Mobile Mapping System (LMMS) in mapping a stretch of sandy Dutch coast of 6 km near the municipality of Egmond is evaluated in this research. A test data set was obtained by provider Geomaat using the StreetMapper LMMS system. Both the relative quality of laser point heights and of a derived Digital Terrain model (DTM) are assessed. First, the height precision of laser points is assessed a priori by random error propagation, and a posteriori by calculating the height differences between close-by points. In the a priori case, the result is a theoretical laser point precision of around 5 cm. In the a posteriori approach it is shown that on a flat beach a relative precision of 3 mm is achieved, and that almost no internal biases exist. In the second analysis, a DTM with a grid size of 1 m is obtained using moving least squares. Each grid point height includes a quality description, which incorporates both measurement precision and terrain roughness. Although some problems remain with the scanning height of 2 m, which causes shadow-effect behind low dunes, it is concluded that a laser LMMS enables the acquisition of a high quality DTM product, which is available within two days.

1. Introduction

The Dutch coast typically consists of a relatively flat sandy beach lined on the land side by dunes, which are partly covered by marram grass. This coastal area is important for the Netherlands for many reasons, e.g., as recreational and nature area, and as protection against sea floods and storms. The last usage is especially crucial, because the most densely populated areas in the Netherlands are located just behind the coastal defense and are partly below the mean sea level. Therefore, it is essential to continuously monitor and maintain the coast in order to protect the Dutch hinterland from the sea. In 1990 a national coastal policy was adopted, with the aim of maintaining the seaward position of the coastline, as it was on 1 January 1990 [1]. To successfully maintain this so-called Basal Coast Line a suitable acquisition technique to measure beach morphology and its changes needs to be employed. Because high energy events like storms may cause large changes, as for example shown in Figure 1, the main interest is to monitor coastal topography on the temporal and spatial scale of storm impacts [2]. Therefore, a flexible system is needed that can access a damaged area immediately after the storm and provide the results of morphologic changes as quickly as possible (in one day). Besides, to assess in detail the beach erosion caused by heavy storm events, high spatial resolution measurements are needed. In this research we focus on assessing possible damage at the intertidal and supra tidal area of the beach, i.e., those parts that are accessible from land, at least during part of the day. A jetski supplied with echo sounder, Global Positioning System (GPS) and Inertial Measurement Unit (IMU) is a good example of a flexible system to assess morphological changes in the near shore zone, [2].

Figure 1. A real example of dune erosion on the Dutch coast on Ameland and the possible consequence; photo taken by Johan Krol [Archive Natuurcentrum Ameland].

Figure 1. A real example of dune erosion on the Dutch coast on Ameland and the possible consequence; photo taken by Johan Krol [Archive Natuurcentrum Ameland].

Since 1996 the Dutch Ministry of Transport, Public Works and Water Management (RWS, Rijkswaterstaat) annually measures the beach topography by means of Airborne Laser Scanning (ALS). The ALS technique has limitations in projects where cost-effective capturing of 3D data and dense point coverage of vertical features are required (e.g., steep dune slopes). Besides, the ALS data in general cannot be provided on demand; firstly because flying permissions are needed, and secondly, after-storm weather conditions may hinder or prevent the acquisition. To summarize, the ALS method offers good results in terms of quality and reliability [3], but is not flexible.

One of the potential alternative techniques is a Land-based Mobile Mapping System (LMMS). LMMS is a complex real-time, multi-tasking and multi-sensor system, which integrates: (i) a number of line scanners and/or digital cameras for surface mapping; (ii) GNSS for positioning; and (iii) additional sensors like for example Inertial Navigation System (INS) to measure the attitude of the vehicle. Those sensors are usually mounted on a rigid platform, placed on the roof of a vehicle. The LMMS mapping sensors can be of different type and orientation, which makes every LMMS system unique in terms of performance and thus quality of acquired data. For an overview of the early LMMS see [4]. More recent LMMS and system providers are described in [5,6,7]. The LMMS discussed in this research comprises three laser scanners as mapping sensors and an integrated GPS/INS as the main navigation sensor.

Using laser LMMS it is in principle possible to quickly obtain 3D geo-referenced data of a large extended area, such as a beach [8]. The main requirement is that the area to be mapped is accessible by the moving platform. On the beach, a 4WD car is most suited. Dunes however have often strong local relief and are more vulnerable. Therefore it is more difficult to access dunes, even by a quad vehicle. Above static terrestrial laser scanning [9], LMMS has the advantage that tenths of kilometers of coast can be sampled in one day. High frequency laser pulse measurements enable high spatial resolution. In comparison with ALS, higher point density is expected from LMMS, because the measured ranges are smaller. On the other hand, more data voids might occur behind elevated features when measuring from the ground. Besides, attention must be paid to the intersection geometry of the laser beam with the relatively horizontal beach. When scanning a horizontal surface, the geometry gets poorer further away from the trajectory. This decreases the laser point position precision. In order to test the laser LMMS performance on the Dutch coast, RWS initiated a pilot-project. Particular interest of the RWS is the level of obtainable accuracy and processing time of a final topographic product, which is a Digital Terrain Model (DTM). The RWS requirements are twofold. First a vertical DTM accuracy of at least 10 cm at a grid spacing of 1 × 1 m is required, and second, it is required that the results are available close to real-time.

In this research the quality of the derived LMMS laser point cloud and DTM is analyzed. The primary objective is to investigate the contribution of various error components to the performance of the laser LMMS and to estimate the quality of the individual laser points. The secondary objective is, to evaluate the quality of the derived DTM.

2. An Overview of Processing Steps

The methodology chosen for this research consists of several processing steps as outlined in Figure 2. In general it is important to know the quality of individual laser points prior to using points in further processing, like computing a DTM. The quality of a laser point position depends on the quality of parameters included in the mathematical georeferencing model, which is used to calculate the point position. Thus, for a priori estimation of the quality of individual laser points, the random errors of LMMS measurements and calibration parameters are propagated through the georeferencing model.

The quality of the calibration parameters is usually known on forehand and depends on the calibration procedure. The actual quality of the individual laser range measurements depends on the instantaneous scanning geometry at the time of acquisition Here the instantaneous scanning geometry is reconstructed by combining the location of the vehicle with the laser measurements itself. By considering for each laser point a small neighborhood of nearby points, it is possible to estimate the angle between terrain and incoming laser ray, which is basically the so-called incidence angle. If, for a given scan point, an incidence angle is not perpendicular, the scanning geometry was not optimal and therefore a higher value of the precision if incorporated. How much higher this additional value should be, is computed using a theoretical model given in [10]. This procedure results for each scan point in a theoretical height precision, depending on: (i) the system specifications after calibration; and (ii) the scan geometry for that point.

Figure 2. Overview of the processing steps to evaluate the quality of LMMS data and products.

Figure 2. Overview of the processing steps to evaluate the quality of LMMS data and products.

To evaluate those theoretical models, a proper Quality Control (QC) using the real data is needed. The result is an empirical (a-posteriori) quality of a laser point position. In [11,12], the existing QC procedures are explained in detail. However, standard and efficient procedures for validating the quality of derived laser points and further on the DTM are still missing. Because the external reference measurements are not available in this research, a QC is developed that estimates the relative precision of the laser points. Namely, the height differences between the so-called identical points, points with overlapping footprints, are taken as a measure for the relative height precision.

Finally, the grid point heights and their precision are obtained by the Moving Least Squares method developed in [13]. The method allows incorporating the theoretical precision of laser point heights in the adjustment model by using a weight matrix W. In the following these processing steps are described in detail.

3. Quality of Laser Point Heights

In Section 3.1, first the scanning geometry at the time of each laser point acquisition is reconstructed by applying simple geometrical rules. The scanning geometry is used firstly in Section 3.2 together with the random errors of LMMS measurements and calibration parameters to estimate the theoretical absolute precision of LMMS laser point heights, and secondly in Section 3.3 to estimate the empirical relative precision of LMMS laser point heights.

3.1. Reconstructing the Scanning Geometry

The instantaneous scanning geometry of a laser point is determined by the range, R, the distance between scanner and object, and by the incidence angle, α, the angle between the objects local surface normal, $\stackrel{\rightharpoonup }{n}$, and the incoming laser ray, $\stackrel{\rightharpoonup }{p}$, see Figure 3. When a beam hits a surface perpendicularly, the incidence angle is 0° while when a beam is almost parallel to a surface the incidence angle approximates 90°. Range and incidence angle together determine the footprint, the area on the object illuminated by the laser pulse. Both incidence angle and footprint size, D_fp, can be computed for each measured laser point using the point and trajectory position as follows.

Figure 3. Scanning geometry attributes.

Figure 3. Scanning geometry attributes.

The range, R, can be computed for each laser point once the sensor position at the time of acquisition is known. The laser scanner position is linearly interpolated using the consecutive trajectory positions, determined by the GPS/INS unit on the scanning vehicle. Here it is assumed that the trajectory position directly represents the laser scanner position; i.e., offsets between the laser scanners and the GPS/INS are not taken into account.

To estimate the incidence angle, α, the direction of the local surface normal has to determined first. The normal vector $\stackrel{\rightharpoonup }{n}$ is computed as follows. For each laser point the four closest points are determined using a k-Nearest Neighbor algorithm [14]. A plane is fitted to all five points using Least Squares resulting in the normal$\stackrel{\rightharpoonup }{n}$ of the plane at the laser point. The number k = 4 of neighboring laser points participating in plane fitting is chosen such that the computed normal vector represents just a local surface.

The area on the surface illuminated by the laser beam is approximated by a circle, whose diameter, D_fp, is computed from the laser beam divergence, β, the incidence angle, α, and range, R, as written in Equation (1):

$$D_{fp}=\frac{R·\beta }{cos\alpha }$$

(1)

3.2. Theoretical Quality of Laser Points

In order to estimate the theoretical or expected precision of the final 3D laser point coordinates the observation equations are required to propagate the random errors. The LMMS mathematical model, the so-called geo-referencing model, relates the system measurements to the ground coordinates of the laser points. This relationship is embodied in the so-called LIDAR equation, as written in Equation (2) [4,8,15].

$${\overrightarrow{r}}_P^m = \overrightarrow{r}{\left(t\right)}_{GPS}^m+R{\left(t\right)}_{imu}^m·({\overrightarrow{r}}_{IMU/S}^{imu} - {\overrightarrow{r}}_{IMU/GPS}^{imu})+ R{\left(t\right)}_{imu}^m·R_S^{imu}·{\overrightarrow{r}}_P^S$$

(2)

Here ${\overrightarrow{r}}_P^m$ is the vector denoting the position of point P in the mapping frame, m; $\overrightarrow{r}{\left(t\right)}_{GPS}^m$ is the vector denoting the position of the GPS antenna in the mapping frame, m. $R{\left(t\right)}_{imu}^m$ is the rotation matrix relating the IMU frame, imu to the mapping frame, m; ${\overrightarrow{r}}_{IMU/S}^{imu}$denotes the offset between the IMU sensor, IMU, and the laser scanner, S, expressed in the IMU frame, while similarly, ${\overrightarrow{r}}_{IMU/GPS}^{imu}$ indicates the offset between IMU and GPS antenna. The position of the point P in the laser scanner frame, s, is given by ${\overrightarrow{r}}_P^S$, while $R_S^{imu}$ denotes the rotation matrix relating the laser scanner frame, s, to the IMU frame, imu. Finally, $R{\left(t\right)}_{imu}^m$ denotes the rotation matrix relating the IMU frame, imu, to the mapping frame, m. For derivations of the LIDAR equation, see e.g., [4,8,15].

By linearizing the geo-referencing Equation (2), the equations of the first order error model are obtained. Now the random errors of the LMMS measurements (i.e., range, scan angle, IMU angles and GPS position) and calibration parameters (i.e., lever-arms and boresight misalignment angles) are propagated through the model. Values of those random errors are taken from the sensors’ specification or are approximated by a value found in the literature [4,8,15]. In addition, the real measurements of range, scan angle and the IMU angles are included into the first order error model. The result is the theoretical or overall expected (a priori) precision of the derived 3D laser point coordinates. In this paper only the Z-component is further analyzed. The height precision of a laser point i due to LMMS measurement errors is called measuring precision and is denoted as σ_Zi,m.

The previously mentioned random range error taken from specifications is valid when the laser beam falls perpendicular to the target. In practice, the incidence angle is changing over the acquisition area and is usually non-perpendicular as shown in Figure 4.

Figure 4. The range error δR due to the non-perpendicular scanning geometry and the influence of δR on vertical and horizontal laser point positioning error.

Figure 4. The range error δR due to the non-perpendicular scanning geometry and the influence of δR on vertical and horizontal laser point positioning error.

High incidence angles result in poor intersection geometry and affect the range measurement precision [10,17,18,19]. For pulse laser scanners, as used in this research, the approach in [20] is used to calculate the range error, δR, as follows:

$$\delta R=\frac{R·\beta ·\text{tan}\alpha }{2}$$

(3)

Here R denotes the range, β the beam divergence and α the incidence angle.

This range error δR is then propagated through the geo-referencing equation. In this way, the height precision of a laser point i with respect to the (non-perpendicular) scanning geometry σ_Zi,δR (geometrical precision) is computed.

Finally, the overall theoretical height precision of a single point i is written as in Equation (4):

$${\sigma }_{Zi}=\sqrt{{\sigma }_{Zi,m}^2+{\sigma }_{Zi,\delta R}^2}$$

(4)

Using the law of error propagation, the theoretical precision, σ_ΔZ, of the height difference, ΔZ, between two laser points q and r is computed as:

$${\sigma }_{\Delta Z}= \sqrt{{\sigma }_{Zq}^2+{\sigma }_{Zr}^2}$$

(5)

3.3. Empirical Quality of Laser Points

The relative quality describes the relation between two points acquired in the same region in a short time period (point-to-point quality) [21]. Here a real laser LMMS data set is used and the empirical quality of point heights is estimated by employing a QC procedure.

In general, the idea of validating the relative quality of laser data is based on checking the compatibility of laser points in areas, where data overlap [22]. In [11,12], some QC procedures are explained. However, the acquisition area discussed in this research, i.e., the sandy beach, does not include (many) steady points or lines, that are sufficiently well defined in the laser LMMS point cloud. In other words, the beach area lacks artificial sharp edges or planes, which could be extracted from the laser points and used in a relative QC procedure. Besides, the terrain on the beach is changing smoothly. Thus, finding and aligning the breaklines of beach morphology is not a promising method either. Instead, the advantage of high LMMS laser point density is used and a point-to-point comparison is made. Namely, the height differences between laser points that lie so close together that their footprints partly overlap, i.e., the height differences between so-called identical points, are analyzed. Not all measured laser points are considered in the process of finding those identical points. The following two conditions are set for the selection of laser points:

• As the footprint diameter goes to infinity when the incidence angle is 90° only laser points that have an incidence angle less than 89.9° are considered: α_P < 89.9°.
• Because just the vertical component of two points is compared, points should lie on an almost horizontal plane in order to avoid the influence of surface slope on the height difference. This requirement is considered to be fulfilled if the vertical component of the unit normal vector n_p,z computed at each laser point, as explained in Section 3.1, is close to 1; that is: n_p,z ≈ 1.

Now pairs of closest points in 3D are found using the kNN algorithm [14] where k = 1. A closest point pair enters the set of identical point pairs, if the 3D distance d_i,j between laser point P_i and its nearest neighbor P_j is smaller than the minimal size of their footprint radii r_i and r_j respectively. Additionally a second threshold is needed, because the footprint size is here approximated with a circle. In general, it is expected that for bigger incidence angles in combination with the increasing range the footprint takes an elliptical shape. For more rough surfaces, the footprint can have any shape in 3D. Therefore, the circle approximation results in an overestimated footprint size when the incidence angle increases. To reduce the influence of this effect, a threshold of 5 cm for the radii of the footprints was experimentally found suitable. Thus:

$$d_{i,j}\le Min\left(min\left(r_i,r_j\right), 5 cm\right)$$

(6)

where i,j = 1 … n, for i ≠ j with n the number of laser points. The height differences ΔZ between the identical points are considered as an empirical relative quality measure. It is expected that the mean of signed height differences ΔZ should approximately equal zero.

LMMS is characterized by a high laser point density as compared to ALS. This high point density has several reasons. First, from an operational viewpoint, the drive paths can be arbitrary close together, resulting in overlapping drive-lines, while the vehicle can also scan at low driving speeds. Besides, usually multiple laser scanners are mounted on a vehicle and are measuring at the same time. It is not clear a priori if points from different drive-lines have the same quality. That is because the acquisition time is different and the configuration of GPS satellites may have changed. Also different scanners may result in points of different quality. Therefore, the height differences ΔZ of identical points are investigated for three different cases:

• Identical points (IP) from the complete data set.
• Identical points (IP) belonging to different scanners (scanner overlap).
• Identical points (IP) belonging to overlapping drive-lines (drive-line overlap).

For each case the height differences of identical points are analyzed in order to estimate noise levels and possibly identify systematic errors. Besides, the correlation with geometric attributes, i.e. the range and incidence angle, of laser points is investigated.

4. DTM Interpolation and Quality

Terrain laser points, which were extracted from the raw data by the data provider Geomaat, are used to interpolate the DTM. The importance of DTM applications makes it inevitable to provide DTMs with adequate quality measures at a high level of detail, as is for example described in [13]. The idea is to inform the user about the DTM quality and warn them of weakly determined areas. Thus, in the following an approach to evaluate the quality of each grid point height is described [6].

In general, the quality of a DTM depends on a number of individual influencing factors, see [23]. The ones investigated here are: the number of terrain points (FD1), height precision of individual terrain point (FD2), terrain point distribution (FD3), terrain roughness (FR) and interpolation method (FI). When the DTM is constructed from the existing laser data, the first three influencing factors (FD1, FD2 and FD3) are usually known or can be estimated. The fourth influencing factor, the terrain roughness (FR), is related to the interpolation method (FI).

There are many different algorithms to interpolate a DTM. The more common are Nearest Neighbor, Inverse Distance Weighting, Moving Least Squares and Kriging [6]. Following the research in [13], it is chosen here to use linear interpolation. Both linear interpolation and Kriging allow proper error propagation. The reason not to use Kriging is because of the characteristics of the input data. In most cases, the amount of observations contributing to one interpolated grid cell elevation is large, that is, in the order of 50 points per cell. Using Ordinary Kriging only a few points close by to the grid point would significantly contribute to the interpolation. More far points would get screened off, compare the discussion on the screen effect in [24]. A natural and easy way to incorporate all observations into a common grid point elevation value is by linear interpolation.

That is, each grid point elevation and its precision are estimated by linear interpolation (FI). Rules of error propagation based on variances and co-variances of the original terrain laser points are applied, to estimate the quality of the grid points. The output is then, strictly speaking, the precision of a grid point, which is denoted by a standard deviation σ_DTM. In other words, the systematic errors are assumed to be zero [13]. First a grid of 1 × 1 m size is laid over the terrain laser points. For grid cells, which include four or more terrain laser points, a tilted plane is fitted in a Least Squares fashion by a first order polynomial as given in Equation (7):

$$Z= a_0+a_{1}X+a_{2}Y$$

(7)

Here X, Y, Z are the coordinates of the terrain laser points (observations) that are included into the plane computation and a₀, a₁ and a₂ are the unknown plane coefficients. The graphical representation of each term in Equation (7) is shown in Figure 5.

Figure 5. A graphic representation of terms given in Equation (7); after [25].

Figure 5. A graphic representation of terms given in Equation (7); after [25].

To make the least squares computation more efficient, a new coordinate system is used with the interpolation grid point (X_G,Y_G) as the origin; therefore, the method is called Moving Least Squares (MLS) adjustment [26]. This simplifies the plane equation as the plane coefficient a₀ becomes the elevation of the grid point itself:

$$Z_G= a_0$$

(8)

The mathematical model of Moving Least Squares for linear surface fitting is then given in matrix vector notation as in Equation (9):

$$\left[\begin{array}{c}{Z}_1\\ Z_2\\ ⋮\\ Z_n\end{array}\right]=\left[\begin{array}{ccc}1& X_1-X_G& Y_1-Y_G\\ 1& X_2-X_G& Y_2-Y_G\\ ⋮& ⋮& ⋮\\ 1& X_n-X_G& Y_n-Y_G\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]$$

(9)

or in short form: y = A.x. In Equation (9), X_i, Y_i, Z_i for i = 1 … n are the coordinates of the n original laser terrain points included in the plane computation. Then the unknowns in vector $\widehat{x}$ (in the short form) and their variance-covariance matrix ${\sum }_{\widehat{x}\widehat{x}}$ are computed in a least squares adjustment as written in Equation (10) and Equation (11) respectively:

$$\widehat{x}={\left(A^T {\text{Σ}}_{yy}^{-1}A\right)}^{-1}{A}^T {\text{Σ}}_{yy}^{-1}y$$

(10)

$${\Sigma }_{\widehat{x}\widehat{x}}={\left(A^T {\Sigma }_{yy}^{-1}A\right)}^{-1}$$

(11)

where ∑_yy is the covariance matrix of observations, in which the theoretical height precision of the laser points σ_Zi computed in Section 3.2 is used.

Besides, the vertical distances between the original terrain points and the modeled plain are computed. The Root Mean Square Error (RMSE) of these residuals e is calculated for each plane, as written in Equation (12), and is an indication of the goodness of fit:

$$RMSE= \sqrt{\frac{e^{T}e}{n}}$$

(12)

To finally predict the DTM quality σ_DTM, a mathematical model after [25] as written in Equation (13) is used.

$${\sigma }_{DTM}^2= {\sigma }_{a0}^2+{\sigma }_e^2$$

(13)

Here the standard deviation of the constant plane coefficient σ_a0 represents the quality of the original data and accounts for the precision of the original laser points (FD2), their density (FD1) and distribution (FD3). The second term σ_e represents the quality loss due to the representation of the terrain surface by a plane. In this research, the RMSE is considered as a measure of the terrain surface roughness (FR) with respect to the plane modeled by the chosen random-to-grid moving least squares interpolation (FI). Therefore, σ_e simply equals the RMSE as computed in Equation (12).

5. Results and Discussion

5.1. Data Description

The LMMS data set was acquired on the Dutch coast near Egmond aan Zee using the StreetMapper system owned by provider Geomaat [27,28]. The acquisition took place on 27 November 2008 at the time of low tide. Within two hours a stretch of beach of 6 km long and 180 m wide was covered. The point cloud consists of about 56 million laser points. As experienced by Geomaat, the 3D laser point coordinates and the classification into terrain and non-terrain points could be done within two days. In this research, a smaller representative test area of 213 × 101 m was chosen, as presented in Figure 6. The data set consists of 1,220,825 laser points. Each record of a laser point has 15 attributes: Nine original attributes written in the *.las file of laser points and 6 additional attributes computed as explained in Section 3. The attributes are: 3D laser point position X,Y,Z, intensity I, class number C, scan angle Θ, time of point acquisition T, drive-line number DL and scanner number SC, range R, incidence angle α, footprint diameter D_fp, range error due to the scanning geometry δR, measuring precision σ_Zm and geometrical precision σ_Z,δR. Besides the 3D laser point data, the data of trajectories given in *.trj file is used. The positions of eight trajectories within the test area are shown in Figure 6.

Figure 6. Left, location of the test area at the west coast of The Netherlands. Right, aerial view of the test area [GoogleMaps]. Coordinates are in RDNAP. The black dashed lines mark the trajectories driven downward, i.e., from the north to the south, and the solid lines mark the trajectories driven in the opposite direction.

Figure 6. Left, location of the test area at the west coast of The Netherlands. Right, aerial view of the test area [GoogleMaps]. Coordinates are in RDNAP. The black dashed lines mark the trajectories driven downward, i.e., from the north to the south, and the solid lines mark the trajectories driven in the opposite direction.

5.2. Results of Theoretical Precision

In Figure 7 the distribution of the measuring precision σ_Zi,m and geometrical precision σ_Zi,δR are shown. Laser points within the whole test area are considered. Because the geometrical precision σ_Zi,δR increases to infinity when the incidence angle approaches 90°, the computation of statistical measures in Figure 7(b) considers just points that have incidence angle α smaller than 89.9°. The median measuring precision σ_Zi,m is 2.63 cm, while the median geometrical precision σ_Zi,δR is much smaller and is 0.86 cm. On the other hand the dispersion of σ_Zi,δR is as expected much bigger than for the σ_Zi,m. The minimum measuring precision σ_Zi,m is 2.54 cm, which is due to the main error contributor that is the GPS error.

Figure 7. Theoretical precision of the laser point heights, computed by the first random error model. (a) Measuring precision, (b) Geometrical precision.

Figure 7. Theoretical precision of the laser point heights, computed by the first random error model. (a) Measuring precision, (b) Geometrical precision.

5.3. Results of Height Differences of Identical Points

By analyzing the attributes of the identical point pairs, i.e., the scanner and drive-line number, it was found that the majority of the identical point pairs belong to the same scanner and the same drive-line. Fewer identical points are found in the scanner overlap and drive-line overlap. In Table 1 the results of height differences for the three cases are presented. The mean (Avg) of height differences ΔZ is very close to zero for all three cases. This indicates that there is no bias in the data, i.e., neither scanner calibration errors nor offset between drive-lines is present. The standard deviations (Std) are equal to or smaller than 3.5 mm, which denotes the relative precision of LMMS laser points.

Table 1. Statistics of the height differences of identical points.

		ALL	Scanner overlap	Drive-line overlap
No. of identical point pairs		17,754	608	5,473
Height difference ΔZ [mm]	Min	−47	−20	−47
	Max	46	36	46
	Avg	0.1	0.2	0.0
	Std	3.1	2.5	3.5

Table 1. Statistics of the height differences of identical points.

		ALL	Scanner overlap	Drive-line overlap
No. of identical point pairs		17,754	608	5,473
Height difference ΔZ [mm]	Min	−47	−20	−47
	Max	46	36	46
	Avg	0.1	0.2	0.0
	Std	3.1	2.5	3.5

Figure 8. The relation of the absolute height differences abs(ΔZ) and scanning geometry attributes. (a) Range R (b) Incidence angle α.

Figure 8. The relation of the absolute height differences abs(ΔZ) and scanning geometry attributes. (a) Range R (b) Incidence angle α.

For the correlation analysis between the height differences ΔZ and the scanning geometry, one value of range and incidence angle is taken per IP pair. Namely, the maximum value of the two points in the IP pair is taken, because it is expected that the larger those values are, the bigger the ΔZ. In Figure 8 the absolute height differences abs(ΔZ) are shown in relation with the IP attributes range max(R) and incidence angle max(α); gray dots are identical point pairs belonging to the same class and the blue dots represent identical point pairs belonging to different classes. To show the changes of abs(ΔZ) more clearly, the mean (black squares) is computed in certain intervals, i.e., for each bin of IP attribute values. The standard deviation values represented by the red error bars show the spread of the absolute height differences abs(ΔZ) within the corresponding bins. In Figure 8(a), a large mean of absolute height differences can be observed at the range of about 5 m. Most of the identical points with the IP range around 5 m and the absolute height difference bigger than 1 cm belong to different classes (see blue dots). When the IP range increases from 7 m to 45 m the mean values slowly increase from 0.8 mm to 2.5 mm. In Figure 8b the average absolute height differences are almost constant for incidence angles changing from 72° to 87°. The average values are around 1.1 mm. Thus, the analyses of the correlation between height differences of identical points lying on a horizontal surface and the geometric attributes, range and incidence angle, do not show a clear trend.

5.4. Comparison of Empirical and Theoretical Height Precision

The theoretical precision of height differences σ_Δz is computed for each height difference between identical points, according to Equation (5). In Table 2 the statistics of empirical and theoretical precision measures are given. The comparison of theoretical and empirical precision of height differences between laser points shows large differences. The theoretical RMSE is approximately 28 times larger. This could partly be expected, because the estimation of the theoretical height precision relies on many assumptions (e.g., about calibration parameters, scan angle error, error due to non-perpendicular scanning geometry).

Table 2. Statistics of the height differences of identical points.

		min	max	mean	std	RMSE
Empirical	ΔZALL [m]	−0.0470	0.0460	0.0001	0.0031	0.0031
Theoretical	σΔz [m]	0.0376	4.8515	0.0573	0.0658	0.0872

Table 2. Statistics of the height differences of identical points.

		min	max	mean	std	RMSE
Empirical	ΔZALL [m]	−0.0470	0.0460	0.0001	0.0031	0.0031
Theoretical	σΔz [m]	0.0376	4.8515	0.0573	0.0658	0.0872

The theoretical error is an indication of the absolute height precision of the points based on nominal precision of the individual sensors and the scanning geometry. To approve the theoretical precision measures, an external source of reference measurements should be employed.

The empirical error is an indication of relative precision and internal consistency of the points in the point cloud. Smaller empirical errors may be due to the constraints applied in the selection of identical point pairs; e.g., points lying close to each other. Those points were on average acquired in a short time interval and share almost the same scanning geometry.

5.5. Results of DTM Interpolation and Precision Estimation

Within the test area the terrain laser points, as classified by Geomaat, were used in the following DTM analysis. In Table 3 the statistics of the factors FD1 (n—number of points per grid cell) and FD2 ($\overline{{\sigma }_{Zi}}$—the average laser point height precision per grid cell), which directly influence the DTM quality, are presented. Besides, an overview of the results from the Moving Least Squares interpolation is given. Robust statistics (median med and robust standard deviation rstd) are used to compare the values of different variables.

In Figure 9 the 3D surface of the interpolated DTM is shown. In this raster image each pixel represents a 1 × 1 m grid cell and the pixel color shows the corresponding grid point height Z_MLS. The grid point elevation is changing from −0.19 m at the coastline to up to 22 m in the dunes (see Table 3). The white holes in the DTM are results of the shadow-effect (white holes in green area) and most probably of the presence of water bodies on the beach (white holes in the blue area).

Table 3. Statistics of the parameters involved in the Moving Least Squares interpolation.

	n (FD1)	$\overline{{\sigma }_{Zi}}$ (FD2)	ZMLS	σa0	σe	σDTM
		[m]	[m]	[m]	[m]	[m]
min	4	0.026	−0.19	0.0015	0.0000	0.0018
max	333	0.85	21.88	2.9	0.1001	2.9
med	69	0.033	1.25	0.0042	0.0008	0.0047
rstd	76	0.008	1.28	0.0030	0.0006	0.0030

Table 3. Statistics of the parameters involved in the Moving Least Squares interpolation.

	n (FD1)	$\overline{{\sigma }_{Zi}}$ (FD2)	ZMLS	σa0	σe	σDTM
		[m]	[m]	[m]	[m]	[m]
min	4	0.026	−0.19	0.0015	0.0000	0.0018
max	333	0.85	21.88	2.9	0.1001	2.9
med	69	0.033	1.25	0.0042	0.0008	0.0047
rstd	76	0.008	1.28	0.0030	0.0006	0.0030

Figure 9. Raster image of interpolated DTM grid points visualized in 3D.

Figure 9. Raster image of interpolated DTM grid points visualized in 3D.

In Figure 10 the number of laser terrain points, n, per 1 × 1 m grid cell is shown. The positions of the drive-lines are drawn as in Figure 6. Grid cells that contain less than three terrain points are colored pink. As explained in Section 4, the Moving Least Squares adjustment of a plane is only possible, if a grid cell contains at least four terrain points. As a consequence, the pink grid cells in Figure 10, 11% of the total, did not enter the MLS computation. The grid cells on the driving path have the highest point density, with a maximum of 333 points per square meter. The number of points per grid cell on the flat beach drops rapidly with increasing distance to the sensor.

Figure 10. Number of points n per 1 × 1 m grid cell. The black lines mark the trajectories driven by the scanning vehicle, compare Figure 6.

Figure 10. Number of points n per 1 × 1 m grid cell. The black lines mark the trajectories driven by the scanning vehicle, compare Figure 6.

The spatial variation of the DTM data quality component σ_a0 is presented in Figure 11(a). The size of σ_a0 depends mostly on the number of points n (FD1) and the quality of the individual terrain laser point σ_zi (FD2). The standard deviation σ_a0 is the smallest, i.e. below 2 mm, within the driving path of drive-line DL5 (light blue strip on the right side of the figure). The main reason is the high number of terrain laser points, which is more than 250 points along this drive-line (see Figure 10). The green color indicates grid cells that have a standard deviation σ_a0 of about 4.2 mm. The black colored points correspond to 10% of the analyzed grid cells with a standard deviation σ_a0 of larger than 2.32 cm. These lie mostly in the dune area and on the edges of the drive-lines.

The second component in Equation (13) for DTM quality computation is the terrain roughness σ_e. It is represented by the root mean square error (RMSE) of the vertical residuals between terrain laser points and the fitted planes. With a given grid size, which is here 1 × 1 m, and a functional model to represent the surface, which is here the tilted plane, σ_e depends mainly on the complexity of the terrain surface. The quality of the measurements is assumed to be high enough. The spatial variability of the terrain roughness component σ_e is shown in Figure 11(b). This pattern is almost independent of the laser point height precision, although vertical stripes are recognizable in the direction of the trajectories. The pattern shows more distinctly the morphology of the terrain (see also Figure 6). In the pre-dune area higher values of the terrain roughness component σ_e are present. About 10 % of the grid cells have a terrain roughness component σ_e larger than 4.5 mm. These are depicted in black color and are located mostly in the pre-dune and dune area. Therefore, these parts of the test area are considered to have a rougher topography, and are thus more difficult to model with a planar surface of 1 × 1 m size.

Finally, in Figure 11(c) the spatial variation of the grid point height precision σ_DTM over the test area is shown. The height precision σ_DTM of the grid points, as computed by Equation (13), varies between 0.0018 and 2.9 m. The average precision of grid points σ_DTM equals to 4.7 mm. For comparison, the precision of the observations σ_Zi is on average 2.4 cm. The grid cells having a standard deviation σ_DTM higher than 2.56 cm are colored dark red and black. They represent approximately 10 % of all the grid cells. The green color shows grid points having a height precision σ_DTM smaller than 1 cm. Most of the beach area has a high DTM quality, which decreases with increasing distance from the trajectory. For example, the precision at the edges of the drive-line DL11 (the leftmost line) decreases and is in some areas worse than 2.56 cm, mostly due to the number of points n (compare to the dark blue areas in Figure 10). The DTM quality is lower also in the dune area, due to the low quality of the terrain laser points (compare to Figure 11(a)) and the high terrain roughness (compare to Figure 11(b)).

In Figure 12 the correlation between the height precision of grid points σ_DTM (y-axis), the number of points n (x-axis) and the data quality component σ_a0 (colorbar), is presented. A comparison of the colorbar and the y-axis scale shows, that the size of the grid point height precision σ_DTM depends mainly on the data quality component σ_a0. Besides, one can observe that if approximately 50 or more points are included in the grid point computation, the standard deviation of the grid point heights σ_DTM drops below 1 cm.

Figure 11. The two components directly employed in the computation of the grid point height precision and the final grid point height precision. Values are shown per 1 × 1 m grid cell. (a) Data quality component σ_a0, (b) Terrain roughness component σ_e, (c) Grid point height precision σ_DTM.

Figure 11. The two components directly employed in the computation of the grid point height precision and the final grid point height precision. Values are shown per 1 × 1 m grid cell. (a) Data quality component σ_a0, (b) Terrain roughness component σ_e, (c) Grid point height precision σ_DTM.

Figure 12. Correlation between the grid point height precision σ_DTM and the number n of terrain laser points; color-coded by the data quality component.

Figure 12. Correlation between the grid point height precision σ_DTM and the number n of terrain laser points; color-coded by the data quality component.

6. Summary, Conclusion and Recommendations

6.1. Summary

In this article the quality of LMMS laser points and the derived DTM was evaluated and analyzed. The quality of the individual LMMS laser point was estimated as a height precision. Two approaches were used. In the first a priori approach a theoretical height precision of the laser points was defined by error propagation of the random errors through the geo-referencing formula. The random errors considered in this first order error model result arise from:

• The LMMS measurements resulting in the measurement precision of laser points.
• The non-perpendicular scanning geometry resulting in the geometrical precision of laser points.

The error components entering the geo-referencing equation and therefore influencing the measuring precision are: GPS, IMU, range, scan angle, lever arm from laser scanner to IMU and from IMU to GPS antenna, and boresight misalignment angle components. In the second approach, a relative QC procedure was developed, employing height differences between the so-called identical points, to assess the empirical laser point height precision.

6.2. Conclusions

The median a priori measuring precision was 2.88 cm, which was dominated by the GPS error. The QC results showed very small numbers: the average precision of the height differences between identical points lying on the horizontal plane was 3 mm. Because the mean of the signed height differences was almost zero for three different analyzed cases, it was concluded that almost no systematic errors were present in the derived laser point cloud. Within the QC of identical points an attempt was made to show the influence of the scanning geometry on the laser point quality. Results show that height differences between identical points do not depend on the range, nor on the incidence angle.

The comparison of the theoretical precision and empirical precision of height differences shows large differences. The reason is that the theoretical precision includes also the GPS/INS positioning error and is therefore considered as an absolute precision. On the contrary, the empirical precision accounts for the relative precision between two laser points. Thus, choosing the nearest neighbors might already minimize the empirical error.

The average height precision at the grid points of the DTM was equal to 4.7 mm. The computation was performed within the weighted MLS adjustment, using the theoretical height precision of the terrain laser point as weights. It was found that the main influencing factor on the grid point height precision is the density of the terrain laser points, because a higher number of observations (i.e., terrain laser points) enabled partly the elimination of the observation noise. The consequence is that the height precision at the grid points improved with an increasing number of terrain laser points, and exceeded the theoretical height precision of the individual terrain laser points. Rijkswaterstaat required a 1 × 1 m DTM having a precision better than 10 cm. Thus, it was concluded that those requirements can be easily met by employing a laser LMMS.

6.3. Recommendations

To further investigate the reliability of the method of identical points, more research should be done on the identical points assumption. It is also recommended to analyze a larger number of identical points. This would enable additional analysis of many subclasses of identical points. The theoretical precision of the individual laser points should be verified by using reference data. To verify the influence of the scanning geometry on the laser point quality it is recommended to additionally measure control points across the drive-lines and compare them with the point cloud.

The adjustment method for the DTM quality estimation includes just the grid cells with more than three terrain laser points and gives, strictly speaking, the precision of the grid points. Using another method that allows the computation of the precision for all grid cells will result in a slightly higher coverage. To make optimum use of the available data it is recommended that this method is made adaptive to both point density and surface relief. Besides, areas without any terrain laser points, resulting from the shadow-effect or surfaces covered with water, must be separately analyzed, see for example [29]. On the other hand, to assess the absolute positional and height accuracy of the DTM product, external reference data of higher accuracy should be used.

A last recommendation for the projects that are concerned with the assessment of sandy beach morphology is to place laser scanners on a higher platform. The StreetMapper platform of 2 m above the ground resulted in quite some data gaps due to occlusions behind the pre-dunes. Based on the DTM visibility analysis, as given in [25], for a particular area of interest the optimal height of the laser scanner(s) above the ground could be calculated.