Remote Sens. 2011, 3, 1472-1491; doi:10.3390/rs3071472
OPEN ACCESS
Remote Sensing
ISSN 2072-4292
www.mdpi.com/journal/remotesensing
Article
Evaluation of a LIDAR Land-Based Mobile Mapping System
for Monitoring Sandy Coasts
Maja Bitenc 1,*, Roderik Lindenbergh 2, Kourosh Khoshelham 2 and A. Pieter van Waarden 3
1
2
3
Geodetic Institute of Slovenia, Jamova 2, 1000 Ljubljana, Slovenia
Department of Remote Sensing, Delft University of Technology, Kluyverweg 1, 2629 HS Delft,
The Netherlands; E-Mails: r.c.lindenbergh@tudelft.nl (R.L.); K.KhoshElham@tudelft.nl (K.K.)
Dutch Ministry of Public Works and Water Management, Data ICT Dienst, Derde Werelddreef 1,
2622 HA Delft, The Netherlands; E-Mail: pieter.van.waarden@rws.nl
* Author to whom correspondence should be addressed; E-Mail: bitenc.m@gmail.com;
Tel.: +386-1-200-29-00; Fax: +386-1-425-06-77.
Received: 1 April 2011; in revised form: 2 June 2011 / Accepted: 22 June 2011 /
Published: 8 July 2011
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.
Remote Sens. 2011, 3
1473
Keywords: coast; hazards; laser scanning; mobile; quality; DEM/DTM; mapping
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].
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
Remote Sens. 2011, 3
1474
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-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.
Remote Sens. 2011, 3
1475
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.
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
Remote Sens. 2011, 3
1476
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, , and the incoming laser ray, , 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, Dfp, can be computed for each
measured laser point using the point and trajectory position as follows.
Figure 3. Scanning geometry attributes.
LS
R
α
β
Θ
Dfp
–
–
–
–
–
–
–
–
Laser scanner
Laser ray vector from scanner to object
Range, length of vector
Local surface normal of the object
Incidence angle
Laser beam divergence
Scan angle
Diameter laser footprint
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 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 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,
Dfp, is computed from the laser beam divergence, β, the incidence angle, α, and range, R, as written in
Equation (1):
∙
=
(1)
Remote Sens. 2011, 3
1477
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].
�
= ( )�� +
∙(
�� /
−
�� /��
)+
( )
∙
∙
�
(2)
Here � is the vector denoting the position of point P in the mapping frame, m; ( )�� is the vector
denoting the position of the GPS antenna in the mapping frame, m.
is the rotation matrix
relating the IMU frame, imu to the mapping frame, m; �� / denotes the offset between the IMU
sensor, IMU, and the laser scanner, S, expressed in the IMU frame, while similarly, �� /�� indicates
the offset between IMU and GPS antenna. The position of the point P in the laser scanner frame, s, is
given by � , while
denotes the rotation matrix relating the laser scanner frame, s, to the IMU
frame, imu. Finally, ( )
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.
Remote Sens. 2011, 3
1478
High incidence angles result in poor intersection geometry and affect the range measurement
precision [10,17-19]. For pulse laser scanners, as used in this research, the approach in [20] is used to
calculate the range error, δR, as follows:
∙ ∙tan
� =
(3)
2
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):
� =
� 2 , + � 2 ,�
(4)
�Δ =
�2 + �2
(5)
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:
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:
i.
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°.
ii.
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 np,z
computed at each laser point, as explained in Section 3.1, is close to 1; that is: np,z ≈ 1.
Remote Sens. 2011, 3
1479
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 di,j between laser point Pi and its
nearest neighbor Pj is smaller than the minimal size of their footprint radii ri and rj 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:
,
≤�
,
,5
(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:
i. Identical points (IP) from the complete data set.
ii. Identical points (IP) belonging to different scanners (scanner overlap).
iii. 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
Remote Sens. 2011, 3
1480
(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):
=
0
+
1
+
(7)
2
Here X, Y, Z are the coordinates of the terrain laser points (observations) that are included into the
plane computation and a0, a1 and a2 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].
To make the least squares computation more efficient, a new coordinate system is used with the
interpolation grid point (XG,YG) as the origin; therefore, the method is called Moving Least Squares
(MLS) adjustment [26]. This simplifies the plane equation as the plane coefficient a0 becomes the
elevation of the grid point itself:
�
=
(8)
0
The mathematical model of Moving Least Squares for linear surface fitting is then given in matrix
vector notation as in Equation (9):
1
2
⋮
1
1
=
⋮
1
−
2−
⋮
−
1
�
�
�
−
2−
⋮
−
1
�
�
�
0
1
2
(9)
Remote Sens. 2011, 3
1481
or in short form: y = A.x. In Equation (9), Xi, Yi, Zi for i = 1 … n are the coordinates of the n original
laser terrain points included in the plane computation. Then the unknowns in vector (in the short
form) and their variance-covariance matrix ∑ are computed in a least squares adjustment as written
in Equation (10) and Equation (11) respectively:
�
= � Σ −1 �
−1
= � � −1 �
� Σ −1
−1
(10)
(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:
�
=
(12)
To finally predict the DTM quality σDTM, a mathematical model after [25] as written in
Equation (13) is used.
� 2 � = � 20 + � 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 Dfp, 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.
Remote Sens. 2011, 3
1482
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.
Number of points
Number of points
Figure 7. Theoretical precision of the laser point heights, computed by the first random
error model. (a) Measuring precision, (b) Geometrical precision.
(a) Precision σZi,m [m]
(b) Precision σZi,,δR [m]
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-
Remote Sens. 2011, 3
1483
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.
No. of identical point pairs
Min
Height difference
Max
Avg
ΔZ [mm]
Std
ALL
17,754
−47
46
0.1
3.1
Scanner overlap
608
−20
36
0.2
2.5
Drive-line overlap
5,473
−47
46
0.0
3.5
Figure 8. The relation of the absolute height differences abs(ΔZ) and scanning geometry
attributes. (a) Range R (b) Incidence angle α.
(a)
(b)
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
Remote Sens. 2011, 3
1484
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
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.
Empirical ΔZALL [m]
Theoretical
σΔz [m]
min
max
mean
std
RMSE
−0.0470
0.0376
0.0460
4.8515
0.0001
0.0573
0.0031
0.0658
0.0031
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
( —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 ZMLS. 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).
Remote Sens. 2011, 3
1485
Table 3. Statistics of the parameters involved in the Moving Least Squares interpolation.
n (FD1)
min
max
med
rstd
4
333
69
76
ZMLS
(FD2)
[m]
0.026
0.85
0.033
0.008
[m]
−0.19
21.88
1.25
1.28
σa0
[m]
0.0015
2.9
0.0042
0.0030
σe
[m]
0.0000
0.1001
0.0008
0.0006
σDTM
[m]
0.0018
2.9
0.0047
0.0030
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.
Remote Sens. 2011, 3
1486
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.
Remote Sens. 2011, 3
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.
(a)
(b)
(c)
1487
Remote Sens. 2011, 3
1488
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
Remote Sens. 2011, 3
1489
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.
Acknowledgements
The authors would like to thank the Dutch Ministry of Transport, Public Works and Water
Management, for kindly giving us the StreetMapper data set. Besides, authors would like to thank the
Geodetic Institute of Slovenia for supporting this research.
Remote Sens. 2011, 3
1490
References and Notes
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
17.
18.
de Ruig, J.H.M.; Hillen, R. Developments in Dutch coastline management: Conclusions from the
second governmental coastal report. J. Coast. Conserv. 1997, 3, 203-210.
van Son, S.T.J.; Lindenbergh, R.C.; de Schipper, M.A.; de Vries, S.; Duijnmayer, K. Using a
Personal Watercraft for monitoring Bathymetyric Changes at Storm Scale. In Proceedings of
Hydro ’09, Cape Town, South Africa, 10–12 November 2009.
Stockdon, H.F.; Doran, K.S.; Sallenger, A.H. Extraction of lidar-based Dune-Crest Elevations for
use in examining vulnerability of beaches to inundation during hurricanes. J. Coast.Res. 2007, 53,
59-65.
Ellum, C.; El-Sheimy, N. Land-based mobile mapping systems. Photogramm. Eng. Remote
Sensing 2002, 68, 13-17.
Vosselman, G.; Maas, H.G. Airborne and Terrestrial Laser Scanning; Whittles Publishing:
Boca Raton, FL, USA, 2010.
Shan, J.; Toth, C.K. Topographic Laser Ranging and Scanning: Principles and Processing;
Taylor & Francis Group: Boca Raton, FL, USA, 2008.
Petrie, G. An introduction to the technology mobile mapping systems. GEOinformatics Magazine
2010, 13, 32-43.
Barber, D.M.; Mills, J.P. Vehicle Based Waveform Laser Scanning in a Coastal Environment. In
The 5th International Symposium on Mobile Mapping Technology (MMT’07), Padua, Italy,
28–31 May 2007.
Pietro, L.S.; O’Neal, M.A.; Puleo, J.A. Developing terrestrial-lidar-based Digital Elevation
Models for monitoring beach nourishment performance. J. Coast. Res. 2008, 24, 1555-1564.
Lichti, D.D.; Gordon, S.J. Error Propagation in Directly Georeferenced Terrestrial Laser Scanner
Point Clouds for Cultural Heritage Recording. In FIG Working Week 2004, Athens, Greece,
23–30 May 2004.
Habib, A.F.; Al-Durgham, M.; Kersting, A.P.; Quackenbush, P. Error Budget of Lidar Systems
and Quality Control of the Derived Point Cloud. In ISPRS Congress, Beijing, China, 3–11 July
2008; Vol. 37, Part B1, pp 203-209.
Sande, C.v.d.; Soudarissanane, S.; Khoshelham, K. Assessment of relative accuracy of AHN-2
laser scanning data using planar features. Sensors 2010, 10, 8198-8214.
Kraus, K.; Karel, W.; Briese, C.; Mandlburger, G. Local accuracy measures for digital terrain
models. Photogramm. Record 2006, 21, 342-354.
Giaccari, L. Fast K-Nearest Neighbors Search; 25 August 2009. Available online:
http://www.advancedmcode.org/gltree.html (accessed date: 6 July 2011).
Glennie, C.L. Rigorous 3D error analysis of kinematic scanning lidar systems. J. Appl. Geodesy
2007, 1, 147-151.
Alharthy, A.; Bethel, J.; Mikhail, E.M. Analysis and Accuracy Assessment of Airborne Laser
Scanning Systems. In XXth ISPRS Congress, Istanbul, Turkey, 12–23 July 2004; pp. 144-149.
Soudarissanane, S.; Lindenbergh, R.; Menenti, M.; Teunissen, P. Incidence Angle Influence on
the Quality of Terrestrial Laser Scanning Points. In Laser Scanning ’09, Paris, France, 1–2
September 2009.
Remote Sens. 2011, 3
1491
19. Schaer, P.; Skaloud, J.; Landtwing, S.; Legat, K. Accuracy Estimation for Laser Point-cloud
including Scanning Geometry 2007. In 5th International Symposium on Mobile Mapping
Technology (MMT2007), Padua, Italy, 28–31 May 2007.
20. Lichti, D.D.; Gordon, S.J.; Tipdecho, T. Error models and propagation in directly georeferenced
terrestrial laser scanner networks. J. Survey. Eng. 2005, 131, 135-142.
21. Kremer, J.; Hunter, G. Performance of the StreetMapper Mobile LIDAR Mapping System in
―Real World‖ Projects. In Photogrammetric Week’07; Fritsch, D., Ed.; Wichmann Verlag:
Heodelberd, Germany, 2007; pp 215-225.
22. Barber, D.M.; Holland, D.; Mills, J.P. Change detection for topographic mapping using threedimensional data structures. Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci. 2008, 37,
1177-1182.
23. Lu, H. Modelling Terrain Complexity. In Advances in Digital Terrain Analysis; Lecture Notes in
Geoinformation and Cartography; Zhou, Q., Lees, B., Tang, G.-a., Eds.; Springer: Berlin,
Germany, 2008; pp. 159-176.
24. Wackernagel, H. Multivariate Geostatistics: An Introduction with Applications, 3rd ed.; Springer:
Berlin, Germany, 2003.
25. Li, Z.; Zhu, Q.; Gold, C. Digital Terrain Modeling: Principles and Methodology; CRC Press:
New York, NY, USA, 2005.
26. Karel, W.; Kraus, K. Quality parameters of digital terrain models. In Checking and Improving of
Digital Terrain Models/Reliability of Direct Georeferencing; European Spatial Data Research
(EuroSDR): Utrecht, The Netherlands, 2006.
27. StreetMapper.
StreetMapper
Mobile
Mapping
System.
Available
online:
http://www.streetmapper.net (accessed on 6 July 2011).
28. Geomaat. Available online: http://cms.geomaat.pageflow.nl (accessed on 6 July 2011).
29. Kraus, K.; Briese, C.; Attwenger, M.; Pfeifer, N. Quality Measures for Digital Terrain Models. In
XXth ISPRS Congress, Istanbul, Turkey, 12–23 July 2004.
© 2011 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article
distributed under the terms and conditions of the Creative Commons Attribution license
(http://creativecommons.org/licenses/by/3.0/).