1. Introduction
Airborne laser scanning (or lidar) from small unoccupied aircraft systems (UAS) has enjoyed wide adoption in the fields of topographic mapping [
1], agriculture [
2], forestry [
3], and urban mapping [
4], to name but a few. However, this rapid adoption has outpaced careful study of the sensor’s accuracy and repeatability, leading to few publications on the topic [
5,
6]. Alongside sensor accuracy and repeatability lies another chief concern for the acquisition of UAS lidar mapping data: the resulting pattern of laser pulse returns, i.e., the scan pattern, of the multi-beam, “fan-style” lidar sensor often used in UAS lidar mapping (
Figure 1). A thorough study of the scan pattern of the multi-beam lidar sensor (and the only study of its kind, to the authors’ knowledge) is presented in Morales et al. [
7]; however, the problem is approached with the lidar sensor in a static, terrestrial pose. The scan pattern of the multi-beam lidar sensor, particularly from the dynamic, aerial pose (
Figure 1), is the principal topic of the following study.
Mission planning for UAS lidar mapping is a new problem that is scarcely addressed in existing literature. The most basic considerations for mapping a target area with a given lidar sensor—e.g., flying height, swath width, and percent sidelap of adjacent swaths—can be deduced from trigonometry and arithmetic. However, the low altitude pose, coupled with the lightweight and relatively low-cost lidar sensors most commonly used for UAS lidar mapping, present new questions that must be addressed to ensure the quality of the lidar data collected.
Perhaps the most crucial aspect of a lidar mapping mission is point density. Methods for calculating point density are sensor-dependent—or, more specifically, are dependent upon the pulse rate of the sensor and the geometry of the sensor’s emitted laser pulses. A small number of solutions exist for calculating the expected point density of a UAS lidar mapping mission that account for specific model of lidar sensor as well as flight parameters such as flying height, forward speed, and percent overlap [
8,
9,
10]. Another key aspect of a UAS lidar mapping mission is the relatively high variation in ranges to target and incidence angles of the laser pulses with objects in the scene. Consider a conventional airborne lidar mapping mission conducted from a piloted aircraft flown many hundreds of meters above ground level. With a typical swath width of 15 to 20° and a flying height on the order of hundreds of meters, the resulting variation of ranges and incidence angles is much lower than those exhibited in a UAS lidar mapping mission, with flying heights on the order of tens of meters and swath widths of 90° or more. Range and incidence angle are key components of the lidar error budget [
11,
12], and careful study must be made of the expected pattern of laser returns—i.e., the scan pattern—of any lidar sensor used for high-accuracy mapping.
The multi-beam “fan-style” lidar sensor presents the unique problem examined in this study. Take for example the Velodyne VLP-16 lidar sensor, currently one of the popular sensors in the UAS lidar mapping field. The VLP-16 scanner head is comprised of sixteen lasers in a “fan” of diverging channels, resulting in a 30° vertical field of view. This “fan” of lasers rotates about the scanner’s vertical axis for a 360° horizontal field of view [
13]. This configuration of lasers is functionally similar to the VLP-16′s predecessors, HDL-64, which was designed for self-driving automobiles [
14]. This fan-style configuration is present in other sensors used in UAS lidar mapping, for example the Velodyne HDL-32 [
5], Quanergy M8 [
15], and the Ouster OS-1-16 [
16].
To ensure maximum coverage along a flight line when using a fan-style sensor in the aerial pose, the most sensible way to orient the scanner is on its side, such that the scanner’s vertical axis is parallel with the direction of travel (
Figure 1). The resulting scan pattern is a series of hyperbolas (
Figure 2). However, the resulting scan pattern does not produce a homogeneous distribution of laser returns (or points), and in some cases, can produce relatively large gaps in coverage (
Figure 3). These gaps in coverage are not a function of point density, but rather point clustering, where many laser pulses are striking the target within very close proximity to each other. This stands in contrast to point dispersion, e.g., a series of laser pulses which produce a more homogeneous pattern across the target.
The VLP-16 emits lasers from sixteen channels oriented between −15° and +15° in 2° intervals from the scanner’s horizontal plane, or xy-plane, resulting in a 30° vertical field of view. The vertical angle ω of each channel is fixed and is defined as the counterclockwise angle of the channel with respect to the scanner xy-plane. The channels fully rotate about the scanner’s vertical axis, or z-axis, for a 360° horizontal field of view. The lasers are fired one at a time according to a precise timing sequence in which one laser is fired every 2.304 μs; after all sixteen lasers have fired, there is a recharge period of 18.43 μs [
17]. Therefore, each laser firing has a unique time, and because the scanner head is in constant rotation, each laser firing has a unique azimuth α, where azimuth is defined as the clockwise angle of rotation about the + z-axis, with the + y-axis set to zero [
17]. Assuming the configuration described above, over flat terrain, the resulting scan lines from the VLP-16 are affine hyperbolas; as the scanner travels forward parallel to the terrain, these overlapping sets of skewed hyperbolas result in areas of varying point dispersion (
Figure 3). Laser pulses are emitted at a constant angular interval (or, put more precisely, emitted at a constant time interval while rotating at a constant angular velocity) along the scanner’s 360° field of view. Under the presented configuration with the VLP-16 turned on its side, point density along a nominally flat surface is a function of linear distance across the scanning profile. The point density, along with the point clustering, are the primary considerations for planning a UAS data collection mission with the VLP-16.
Maximizing not only the accuracy but also the optimal dispersion of the laser returns is crucial for detection of small features, especially for horizontal accuracy assessment of a UAS lidar mapping system. Direct comparison cannot be made between a lidar point cloud and point features in the real world; geometric modeling of surfaces must be utilized to compare lidar data and the scanned area [
18,
19,
20] For horizontal accuracy assessment, virtual points can be constructed via the intersection of surfaces. For practical reasons, especially when not in the presence of common surfaces used for accuracy assessment such as roofs, these surfaces are of finite size. To properly model the surfaces, one must ensure that (1) enough pulses reach the surface to accurately model it [
21] and (2) those pulses are within reasonable thresholds of range and incidence angle. Therefore, the scan pattern of a UAS sensor must be characterized before data collection to ensure the resulting data will meet the standards necessary to perform the accuracy assessment at hand.
This study presents the development and use of VLP-16 simulation software for to examine the scan pattern of the VLP-16 both qualitatively—i.e., manually examining simulated point clouds—and quantitatively, via spatial statistics. Characterization equations for the VLP-16′s point density, optimal separation of flight lines, and possible gaps in coverage are also presented. These parameters can be modified, which facilitate using the software and point cloud simulation for other fan-style scanning systems.
2. Materials and Methods
The primary objective of this study was to characterize and analyze the scan pattern of the VLP-16 as a typical fan-style lidar sensor. To this end, the authors derived characterization equations to predict point density and expected coverage gaps. The authors also created simulation software which models the sensor in the typical UAS pose (
Figure 1) over a target plane, a proxy for flat terrain. Both the characterization equations and simulation software are generalized for those sensor and flight parameters that have the greatest effect on the resulting scan pattern, i.e., flying height, forward speed, rotation rate of the scanner head, and orientation of the scanner with respect to the target plane.
Field data for verification of the equations’ and simulation software’s output were collected on 15 June 2018. The UAS used is a DJI S1000 octocopter with a custom-built lidar mapping payload comprised of the Velodyne VLP-16 lidar sensor (Velodyne Lidar, San Jose, CA, USA) and the original equipment manufacturer (OEM) component equivalent to the NovAtel SPAN-IGM-S1 integrated Global Navigation Satellite System and Inertial Navigation System (GNSS/INS) paired with a NovAtel GPS-702-GG GNSS antenna and a secondary Garmin GPS18x GPS antenna (used for timestamping the lidar pulses) (NovAtel Inc., Calgary, Alberta, Canada). The position data collected by the INS underwent a post-processed kinematic (PPK) solution using a nearby static GNSS reference station to refine its accuracy. All GNSS/INS observations—the raw GNSS observations, the PPK-adjusted GNSS observations, and the raw INS attitude observations—were then processed in a forward/reverse, multi-pass, tightly coupled solution in Waypoint Inertial Explorer software to produce a navigation solution which adjusts and minimizes the errors of the position and attitude observations of the mapping payload onboard the UAS. A custom Python algorithm, using the tightly coupled navigation solution, then (1) transformed the raw Velodyne lidar data form spherical to Cartesian coordinates in the scanner frame as described in [
22], then (2) performed direct georeferencing of the transformed lidar data using the method described in El-Sheimy [
23]. Both the mapping payload and the custom algorithms used to generate the lidar point clouds in this study can be considered typical for this category of lidar data acquisition systems.
In addition to the comparison of the field data with the simulated data and characterization equations, this study also explored a spatial analysis of the simulated scan pattern data to further characterize its non-uniformity. Further detail of these methods is given below.
2.1. Analytical Characterization of VLP-16 Scan Pattern
The point density of the laser return pattern of the VLP-16, as a function of lateral distance from the flight line
, can be closely approximated by the point density function
where
is the height above ground,
is the pulse frequency of the scanner (approximately 300,000 pulses per second),
is the forward velocity of the airframe, and
is the angular difference between the scanner’s z-axis and the direction of travel (herein referred to as the scanner’s “yaw”). This equation is simplified under the assumption that the laser pulses are emitted at a uniform rate, although in practice, this is not the case. Each of the sixteen channels of the VLP-16 emits a pulse once per 2.304
; after all sixteen lasers have fired, a recharge period of 18.43
follows. This assumption, however, has a negligible effect on the results, as demonstrated by comparing the equation’s results to the scan pattern simulation presented in the Results section.
The bands of gaps present in the scan pattern (
Figure 3) occur at certain lateral distances from the flight line
, and can possibly occur at the distances found using the equation
where
is the rotation rate of the scanner head (typically 5–20 Hz),
is the angular separation between adjacent channels (2° for the VLP-16); and
is an integer, indicating the
gap outward from the flight line at
Note that the term for which the arccosine is taken in Equation (2) can be
, which yields a complex value of
.
To assure a minimum density of laser returns along the profile of the laser return pattern of a mission with parallel flight lines, the maximum flight line separation
can be expressed as
where the minimum desired return density is
, expressed as points/m
2. This equation is a further derivation of the point density function. Inspection of the plot of
(Equation (1)) should be used to inform a sensible choice for a value of
.
The derivation of Equations (1)–(3) are given in the
Appendix A.
2.2. Simulation of the VLP-16 Scan Pattern
To test the VLP-16 scan pattern in a controlled environment, software was created in the MATLAB programming language to simulate the laser return pattern of the scanner in the aerial pose. The scanner is modeled as a point
at a user-specified height above a horizontal target plane passing through the origin, traveling at a constant speed (also specified by the user) along a vector parallel to the target plane. The modeled emitted laser pulses are modeled as lines passing through the scanner point toward the target plane. The “vertical” angles
of the lasers (referred to as “channels” by the manufacturer) are fixed, as described above. The “horizontal” angles of the lasers, or azimuths
, are determined by the orientation of the scanner head at each epoch as it rotates about the scanner’s vertical axis. The direction of each modeled pulse in the scanner-oriented coordinate system (SOCS), and the subsequent rotation of SOCS into the mapping frame
(as shown in
Figure 1), is found by
For each laser pulse emitted at epoch
, the simulation passes the line parameters
(the position of which is calculated according to the forward speed and time
) and direction of simulated pulse
to a solver function, which finds the intersection of the line with the target plane, using the common algebraic method. The solutions of these line-plane intersections—a set of points, i.e., a simulated point cloud—are recorded and saved in a text file. This is a simplified sensor model that does not account for the interior geometry of the sensor [
5], the manufacturer reported range accuracy of the sensor [
22], or the incidence-angle-dependent range bias observed in [
11]. For the purposes of this study, however, the simplified sensor model is sufficient to study coverage gaps whose dimensions are two to three orders of magnitude greater than the previously mentioned error sources.
The software takes as input key mission parameters with regards to the orientation and operation of the scanner: height above ground, forward speed, scanner head rotation rate, and yaw (i.e., the difference in the scanner’s +z-axis and the direction of travel). (Yaw is accounted for by generalizing the rotation matrix shown in Equation (4)). The output is a simulated point cloud that is a sample of the laser pulse return that could be expected from a mission flown over flat ground under these mission parameters. The point cloud output is a comma-delimited text file, where each row represents one return, or point. For each point, the following attributes are attached: the mapping frame coordinates of the solved line-plane intersection, ; the azimuth of the scanner head when the pulse was fired, ; the vertical angle, or channel, from which the laser was fired, ; the time at which the laser was fired, ; the range of the laser return (i.e., the distance between the scanner point and the solved line-plane intersection), ; and the direction in the mapping frame that the laser was fired,
The scan pattern analysis and gap equation seek to detect and predict the location of these bands as a function of the mission parameters of flying height, flying speed, and the rotation rate of the scanner head (which can be adjusted anywhere from 5 Hz to 20 Hz). Predicting where these bands of high clustering occur as a function of mission parameters will lead to recommendations for optimal values for the parameters mentioned above, as well as side lap (i.e., the overlap of adjacent swaths).
2.3. Spatial Analysis of Simulated Data
The scan pattern, for all its peculiarities, is in fact a repeating pattern. This allows for the extraction of a narrow across-track profile or “representative profile” (
Figure 4) which can be used for statistical analysis. The apparent clustering and dispersion of the points along the profile is analyzed by binning the representative profile of the simulated point cloud along the across-track axis and calculating the nearest neighbor index (Equations (4)–(6)) for each bin. The nearest neighbor index is a measure applied to point patterns that indicates whether the pattern is clustered, random, or dispersed. First, for each simulated return (point) in the bin
, its nearest neighbor and the distance between the two (
) is found. This is accomplished by calculating a distance matrix (of dimensions
) for
points in the bin, and for each column in the matrix, the lowest off-diagonal value is retained. The mean of those
nearest neighbor distances is the observed nearest neighbor distance, or
:
Next, the expected nearest neighbor distance
is calculated, which is a function of the area
of the bin:
Finally, the z-score (i.e., nearest neighbor index) for each bin is found:
where
.
The z-score provides a normalized measure of the point dispersion or clustering in the bin. For example, a z-score of (with its corresponding p-value of ), indicates with 95% or greater likelihood that the point pattern is clustered, while a z-score of () indicates a 95% or greater likelihood that the pattern is dispersed, i.e., approaching even spacing amongst points (Clark and Evans, 1954; O’Sullivan and Unwin, 2010). It should be mentioned here that this presented interpretation of the z-score is reflective of the source material; however, it is not appropriate for analyzing data of this nature (see the Discussion section).
The analysis of the scan pattern extends beyond its geometric pattern. The point cloud simulation software records other attributes of each simulated return, such as the position of the scanner, range from scanner, time, and direction of the simulated laser pulse. These attributes can provide a simulation of not only of point density and dispersion, but also the expected absolute accuracy of the points. The information that can be gleaned from these attributes—e.g., range and incidence angle—are key components of the lidar error budget, as mentioned in the Introduction. Error propagation within the scanner’s frame (i.e., error related to the scanner itself, relative to the scanner’s coordinate system, not accounting for error from the navigation sensor data) is essential to understanding the reliability of the resulting point cloud. For this study, the simulated returns within each bin are analyzed according to attributes unique to each return. Within each bin, histograms can be generated to show the distribution of the number of returns as a function of average range and average incidence angle for each bin. These two measures are directly proportional over the ideal target plane used in this study’s simulation; as the incidence angle of a laser pulse “increases” (i.e., departs from being normal to the target), so too increases the range that pulse must travel to reach the target. An analog to the accuracy of the points within a given scanned area can be visualized by finding the distribution of ranges and scan angles across the scan profile.