Scan Pattern Characterization of Velodyne VLP-16 Lidar Sensor for UAS Laser Scanning

Many lightweight lidar sensors employed for UAS lidar mapping feature a fan-style laser emitter-detector configuration which results in a non-uniform pattern of laser pulse returns. As the role of UAS lidar mapping grows in both research and industry, it is imperative to understand the behavior of the fan-style lidar sensor to ensure proper mission planning. This study introduces sensor modeling software for scanning simulation and analytical equations developed in-house to characterize the non-uniform return density (i.e., scan pattern) of the fan-style sensor, with special focus given to a popular fan-style sensor, the Velodyne VLP-16 laser scanner. The results indicate that, despite the high pulse frequency of modern scanners, areas of poor laser pulse coverage are often present along the scanning path under typical mission parameters. These areas of poor coverage appear in a variety of shapes and sizes which do not necessarily correspond to the forward speed of the scanner or the height of the scanner above the ground, highlighting the importance of scan simulation for proper mission planning when using a fan-style sensor.


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. 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. Sensors 2020, 20, x 3 of 21  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  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 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.

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.

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 x, can be closely approximated by the point density function where h is the height above ground, l f is the pulse frequency of the scanner (approximately 300,000 pulses per second), v z 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 µs; after all sixteen lasers have fired, a recharge period of 18.43 µs 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 x i , and can possibly occur at the distances found using the equation where r 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 i is an integer, indicating the ith gap outward from the flight line at x = 0. Note that the term for which the arccosine is taken in Equation (2) can be ≥ 1, which yields a complex value of y i . 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 w can be expressed as where the minimum desired return density is p d , expressed as points/m 2 . This equation is a further derivation of the point density function. Inspection of the plot of p(x) (Equation (1)) should be used to inform a sensible choice for a value of p d . The derivation of Equations (1)-(3) are given in the Appendix A.

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 s t 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 R (as shown in Figure 1), is found by Sensors 2020, 20, 7351 For each laser pulse emitted at epoch t, the simulation passes the line parameters → s t (the position of which is calculated according to the forward speed and time t) and direction of simulated pulse r 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 R 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, (X, Y, Z); 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, t ; 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).

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 s i , its nearest neighbor and the distance between the two (d min ) is found. This is accomplished by calculating a distance matrix (of dimensions n × n) for n points in the bin, and for each column in the matrix, the lowest off-diagonal value is retained. The mean of those n nearest neighbor distances is the observed nearest neighbor distance, or d obs : Next, the expected nearest neighbor distance d exp is calculated, which is a function of the area A of the bin: Finally, the z-score (i.e., nearest neighbor index) for each bin is found: where SE = 0.26136/ n 2 /A.
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 ̅ :  The z-score provides a normalized measure of the point dispersion or clustering in the bin. For example, a z-score of ≤ −1.96 (with its corresponding p-value of ≤ 0.05), indicates with 95% or greater likelihood that the point pattern is clustered, while a z-score of ≥ +1.96 (p ≥ 0.95) 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.

Single Swath
As expected and shown in the point density function equation, the point density of the scan pattern is inversely proportional both to the flying height and forward velocity of the scanner. These relationships are shown in Figures 5 and 6, respectively. The simulations are in agreeance with  (2)) made evident through simulation:

1.
Any gaps present at or near nadir to the flight line (i.e., gaps at some low value of x i ) are not always reported by the gap equation.

2.
The severity of the gaps at some lateral distance x i is not reflected in the results of the gap equation. The equation only reports the possibility of gaps occurring at some lateral distance.

3.
At higher forward velocities, some significant bands of gaps may not be reported by the gap equation.
An example of these limitations is shown in Figure 7 (near nadir gaps sometimes not predicted) and 8 (high velocity gaps sometimes not predicted). In conjunction with the gap equation results, inspecting the point density histograms for each bin's z-score (Equation (6)) reveals flaws in using the nearest neighbor index as a proxy for detecting gaps (see Discussion). As observed through simulation (e.g., Figures 7 and 8), in the areas of the greatest degree of point clustering, the gaps in the scan pattern tend to present as either near-rhombuses, narrow and elongated perpendicular to the direction of flight. As the lateral distance from the flight line x i increases, the gaps will sometimes present as somewhat linear clusters of points, which have a less deleterious effect on the coverage of the area. Gaps near nadir to the flight line are elongated to the point of presenting similarly to a widely spaced linear scan pattern, especially as the flying height increases (and the eccentricity of the hyperbolas decrease). An example of these limitations is shown in Figure 7 (near nadir gaps sometimes not predicted) and 8 (high velocity gaps sometimes not predicted). In conjunction with the gap equation results, inspecting the point density histograms for each bin's z-score (Equation (6)) reveals flaws in using the nearest neighbor index as a proxy for detecting gaps (see Discussion). As observed through simulation (e.g., Figures 7 and 8), in the areas of the greatest degree of point clustering, the gaps in the scan pattern tend to present as either near-rhombuses, narrow and elongated perpendicular to the direction of flight. As the lateral distance from the flight line increases, the gaps will sometimes present as somewhat linear clusters of points, which have a less deleterious effect on the coverage of the area. Gaps near nadir to the flight line are elongated to the point of presenting similarly to a widely spaced linear scan pattern, especially as the flying height increases (and the eccentricity of the hyperbolas decrease).  An example of these limitations is shown in Figure 7 (near nadir gaps sometimes not predicted) and 8 (high velocity gaps sometimes not predicted). In conjunction with the gap equation results, inspecting the point density histograms for each bin's z-score (Equation (6)) reveals flaws in using the nearest neighbor index as a proxy for detecting gaps (see Discussion). As observed through simulation (e.g., Figures 7 and 8), in the areas of the greatest degree of point clustering, the gaps in the scan pattern tend to present as either near-rhombuses, narrow and elongated perpendicular to the direction of flight. As the lateral distance from the flight line increases, the gaps will sometimes present as somewhat linear clusters of points, which have a less deleterious effect on the coverage of the area. Gaps near nadir to the flight line are elongated to the point of presenting similarly to a widely spaced linear scan pattern, especially as the flying height increases (and the eccentricity of the hyperbolas decrease). An example of these limitations is shown in Figure 7 (near nadir gaps sometimes not predicted) and 8 (high velocity gaps sometimes not predicted). In conjunction with the gap equation results, inspecting the point density histograms for each bin's z-score (Equation (6)) reveals flaws in using the nearest neighbor index as a proxy for detecting gaps (see Discussion). As observed through simulation (e.g., Figures 7 and 8), in the areas of the greatest degree of point clustering, the gaps in the scan pattern tend to present as either near-rhombuses, narrow and elongated perpendicular to the direction of flight. As the lateral distance from the flight line increases, the gaps will sometimes present as somewhat linear clusters of points, which have a less deleterious effect on the coverage of the area. Gaps near nadir to the flight line are elongated to the point of presenting similarly to a widely spaced linear scan pattern, especially as the flying height increases (and the eccentricity of the hyperbolas decrease).

Rotation Rate and Gaps
The results of the gap equation (Equation (2)) can be initially misleading. As the rotation rate increases, the number of predicted possible locations of gaps increases proportionally. In fact, as the rotation rate doubles, so too does the number of predicted gaps. The predicted gap locations are somewhat linear, occurring every meters along the profile at some rotation rate ; at rotation rate

Rotation Rate and Gaps
The results of the gap equation (Equation (2)) can be initially misleading. As the rotation rate increases, the number of predicted possible locations of gaps increases proportionally. In fact, as the rotation rate doubles, so too does the number of predicted gaps. The predicted gap locations are somewhat linear, occurring every x meters along the profile at some rotation rate r ; at rotation rate 2r, the distance between gaps becomes roughly x/2. But the severity of those gaps is not reflected in the equation. In fact, as the rotation rate increases, the severity of the gaps decreases notably, as shown in Figure 9. This study does not comment on the potential benefits of lowering the rotation rate of the scanner head, but the results do show that a lower rotation rate, especially when coupled with a higher forward velocity, can have detrimental effects on the quality of the resulting coverage.

Rotation Rate and Gaps
The results of the gap equation (Equation (2)) can be initially misleading. As the rotation rate increases, the number of predicted possible locations of gaps increases proportionally. In fact, as the rotation rate doubles, so too does the number of predicted gaps. The predicted gap locations are somewhat linear, occurring every meters along the profile at some rotation rate ; at rotation rate 2 , the distance between gaps becomes roughly /2. But the severity of those gaps is not reflected in the equation. In fact, as the rotation rate increases, the severity of the gaps decreases notably, as shown in Figure 9. This study does not comment on the potential benefits of lowering the rotation rate of the scanner head, but the results do show that a lower rotation rate, especially when coupled with a higher forward velocity, can have detrimental effects on the quality of the resulting coverage.

Effect of Yaw Angle
Adding yaw, i.e., crabbing the scanner with respect to the direction of flight, leads to the coverage gaps becoming narrower along track and elongated across track, with respect to the gaps that would occur under similar conditions without crabbing. This is shown in Figures 10 and 11. The primary deleterious effect of the coverage gaps-increasing the odds of missing small areas of interest, e.g., linear features such as sidewalks or power lines, or small areas of interest such as targets [18]-is negated with the addition of only a small amount of yaw. In other words, the yaw of the scanner leads to a more homogeneous pattern of returns. This has a minimal effect on the width of a single swath, as its width is a function of the cosine of the yaw angle. For example, even a yaw angle of 30° would only yield a 1 − cos 30°= 13% reduction in the width of that swath. This finding aligns with the results of [7], which achieved an analogous result by developing an automatic tilting platform for the VLP-16 in a static, terrestrial pose.

Effect of Yaw Angle
Adding yaw, i.e., crabbing the scanner with respect to the direction of flight, leads to the coverage gaps becoming narrower along track and elongated across track, with respect to the gaps that would occur under similar conditions without crabbing. This is shown in Figures 10 and 11. The primary deleterious effect of the coverage gaps-increasing the odds of missing small areas of interest, e.g., linear features such as sidewalks or power lines, or small areas of interest such as targets [18]-is negated with the addition of only a small amount of yaw. In other words, the yaw of the scanner leads to a more homogeneous pattern of returns. This has a minimal effect on the width of a single swath, as its width is a function of the cosine of the yaw angle. For example, even a yaw angle of 30 • would only yield a 1 − cos 30 • = 13% reduction in the width of that swath. This finding aligns with the results of [7], which achieved an analogous result by developing an automatic tilting platform for the VLP-16 in a static, terrestrial pose.

Effect of Yaw Angle
Adding yaw, i.e., crabbing the scanner with respect to the direction of flight, leads to the coverage gaps becoming narrower along track and elongated across track, with respect to the gaps that would occur under similar conditions without crabbing. This is shown in Figures 10 and 11. The primary deleterious effect of the coverage gaps-increasing the odds of missing small areas of interest, e.g., linear features such as sidewalks or power lines, or small areas of interest such as targets [18]-is negated with the addition of only a small amount of yaw. In other words, the yaw of the scanner leads to a more homogeneous pattern of returns. This has a minimal effect on the width of a single swath, as its width is a function of the cosine of the yaw angle. For example, even a yaw angle of 30° would only yield a 1 − cos 30°= 13% reduction in the width of that swath. This finding aligns with the results of [7], which achieved an analogous result by developing an automatic tilting platform for the VLP-16 in a static, terrestrial pose.

Overlapping Swaths
With the exception of one-way mission plans along some corridor, such as a utility easement or transportation corridor, most lidar data collection missions using UAS will likely feature parallel, overlapping swaths of data. These overlapping swaths not only provide common targets in each swath which can be used for strip adjustment and accuracy assessment, but also can be used to "fill in" the gaps present in the VLP-16 scan pattern. Using the optimal flight line separation equation, it is possible to plan a mission in which a desired minimum point density is achieved across the mission area most efficiently. Figure 12

Overlapping Swaths
With the exception of one-way mission plans along some corridor, such as a utility easement or transportation corridor, most lidar data collection missions using UAS will likely feature parallel, overlapping swaths of data. These overlapping swaths not only provide common targets in each swath which can be used for strip adjustment and accuracy assessment, but also can be used to "fill in" the gaps present in the VLP-16 scan pattern. Using the optimal flight line separation equation, it is possible to plan a mission in which a desired minimum point density is achieved across the mission area most efficiently. Figure 12

Overlapping Swaths
With the exception of one-way mission plans along some corridor, such as a utility easement or transportation corridor, most lidar data collection missions using UAS will likely feature parallel, overlapping swaths of data. These overlapping swaths not only provide common targets in each swath which can be used for strip adjustment and accuracy assessment, but also can be used to "fill in" the gaps present in the VLP-16 scan pattern. Using the optimal flight line separation equation, it is possible to plan a mission in which a desired minimum point density is achieved across the mission area most efficiently. Figure 12

Spatial Analysis
The question remains of which (if any) z-score threshold is indicative of point clustering resulting in undesirable gaps in coverage. The z-score interpretation presented above is based on the distribution of expected nearest neighbor distances in a theoretically random distribution within each bin; the nearest neighbor index is used to estimate the probability that a distribution of points in a given area is clustered, random, or dispersed. The points in the VLP-16′s scan pattern are not random events, which undermines the probabilistic interpretation of the z-score. In fact, the z-scores appear to be fully relative, only useful for comparison within a given flight configuration. Inspection of Figure 9 shows that both the 20 Hz and 5 Hz flights have comparable peak minimum z-scores of −18 and −21, respectively. This belies the severity of the gap problem in the 5 Hz flight while overstating the gap issue in the 20 Hz flight, which is all but negligible.
Final considerations for mission planning include optimizing the expected incidence angles and ranges to the objects in the scene. Figure 13 shows the average absolute values of scan angles of each return within 1-m bins across the profile. Wider spacing of flight lines may lead to fewer flight lines needed to cover an area, but as the scan angle (and thus range) of the returns increases, so too will their error. Another thought is to simply exclude returns beyond a certain range during postprocessing of the point cloud for the sake of a higher accuracy across the resultant point cloud.

Spatial Analysis
The question remains of which (if any) z-score threshold is indicative of point clustering resulting in undesirable gaps in coverage. The z-score interpretation presented above is based on the distribution of expected nearest neighbor distances in a theoretically random distribution within each bin; the nearest neighbor index is used to estimate the probability that a distribution of points in a given area is clustered, random, or dispersed. The points in the VLP-16 s scan pattern are not random events, which undermines the probabilistic interpretation of the z-score. In fact, the z-scores appear to be fully relative, only useful for comparison within a given flight configuration. Inspection of Figure 9 shows that both the 20 Hz and 5 Hz flights have comparable peak minimum z-scores of −18 and −21, respectively. This belies the severity of the gap problem in the 5 Hz flight while overstating the gap issue in the 20 Hz flight, which is all but negligible.
Final considerations for mission planning include optimizing the expected incidence angles and ranges to the objects in the scene. Figure 13 shows the average absolute values of scan angles of each return within 1-m bins across the profile. Wider spacing of flight lines may lead to fewer flight lines needed to cover an area, but as the scan angle (and thus range) of the returns increases, so too will their error. Another thought is to simply exclude returns beyond a certain range during postprocessing of the point cloud for the sake of a higher accuracy across the resultant point cloud.  Figure 12). Note the y-axis difference in the topmost histogram and minor differences in the legend. Figure 14 shows an example scan profile where the maximum range has been limited to 60 m, down from the manufacturer's reported maximum range of 100 m. The profile width decreases drastically, but depending on the desired accuracy of the data, this may be advantageous. Note that the optimal flight line separation equation cannot be used if the maximum range is limited in this manner.  Figure 12). Note the y-axis difference in the topmost histogram and minor differences in the legend. Figure 14 shows an example scan profile where the maximum range has been limited to 60 m, down from the manufacturer's reported maximum range of 100 m. The profile width decreases drastically, but depending on the desired accuracy of the data, this may be advantageous. Note that the optimal flight line separation equation cannot be used if the maximum range is limited in this manner.  Figure 12). Note the y-axis difference in the topmost histogram and minor differences in the legend. Figure 14 shows an example scan profile where the maximum range has been limited to 60 m, down from the manufacturer's reported maximum range of 100 m. The profile width decreases drastically, but depending on the desired accuracy of the data, this may be advantageous. Note that the optimal flight line separation equation cannot be used if the maximum range is limited in this manner.

Recommendations for Mission Planning
The study reveals trends in which mission parameters may lead to coverage gaps. The gaps are more likely to occur as overlap decreases (or overlap is not present), forward velocity increases, and/or the rotation rate of the scanner head decreases. The first two relations make sense: overlap obviously helps fill in any potential gaps; as forward velocity decreases, point density increases. The change in the scan pattern with respect to rotation rate of the scanner head, is not so straightforward.
The gaps were shown to have more potential to provide degrading quality of coverage as the rotation rate of the scanner head decreases. Velodyne allows the user to adjust the rotation rate between 5-20 Hz. A lower rotation rate may lead to increased accuracy of the azimuth reading in the scanner frame and/or decreased power consumption, though this was not examined in this study. Concerns of gaps in coverage, especially with respect to flying height, forward velocity, and the size of objects of interest in the scene, should be considered alongside the potential benefits of decreasing

Recommendations for Mission Planning
The study reveals trends in which mission parameters may lead to coverage gaps. The gaps are more likely to occur as overlap decreases (or overlap is not present), forward velocity increases, and/or the rotation rate of the scanner head decreases. The first two relations make sense: overlap obviously helps fill in any potential gaps; as forward velocity decreases, point density increases. The change in the scan pattern with respect to rotation rate of the scanner head, is not so straightforward.
The gaps were shown to have more potential to provide degrading quality of coverage as the rotation rate of the scanner head decreases. Velodyne allows the user to adjust the rotation rate between 5-20 Hz. A lower rotation rate may lead to increased accuracy of the azimuth reading in the scanner frame and/or decreased power consumption, though this was not examined in this study. Concerns of gaps in coverage, especially with respect to flying height, forward velocity, and the size of objects of interest in the scene, should be considered alongside the potential benefits of decreasing the rotation rate of the scanner head. but random, back-and-forth yaw/crab of the airframe helps to fill in gaps.
The study also suggests that the maximum range of returns be set at some limit below the VLP-16 s empirical maximum range of 120 m. In the authors' experience, this has been found to produce vastly more accurate results. Bear in mind that this finding has been confined to primarily research missions where maximizing area of coverage was not a concern. This concern is a product of range-dependent error and generalizes to all lidar sensors. If time and resources permit, narrowing the field of view should be considered.
The factor that ensures the most evenly dispersed scan pattern when using the fan-style lidar sensor is yaw, i.e., crabbing the scanner with respect to the direction of flight (see Section 3.3). Intentionally misaligning the lidar sensor's rotational axis and the airframe's roll axis when mounting the payload will yield this result at minimal cost to the width of the lidar swath. This is not an option to the end user with a turnkey UAS lidar system in hand; to the authors' knowledge, most UAS lidar payloads are designed such that the rotation axis of the scanner and the roll axis of the airframe are parallel. One possible solution is to exploit wind; planning a mission's flight lines to be perpendicular (or at least not nearly parallel) to the direction of wind will force the UAS autopilot to crab the airframe in order to stay on course.

Conclusions
The fan-style lidar sensor ( Figure 1) has become widely adopted for UAS lidar mapping. The fan-style geometry of lasers in this class of sensors creates a unique pattern of returns which warrants careful study for the purpose of planning successful UAS lidar data collection missions. For this study, an often-used fan-style lidar sensor, the Velodyne VLP-16, was selected. Through the development and use of a novel VLP-16 simulation software, the study examines the sensor's scan pattern of the VLP-16 both qualitatively-i.e., manually examining simulated point clouds-and quantitatively, via spatial statistics. This study also presents novel characterization equations for the VLP-16 s point density, optimal separation of flight lines, and possible gaps in coverage, all of which can be generalized to other fan-style lidar sensors. Flight planning recommendations for minimizing coverage gaps are offered based on examination of the results. The factor found to have the greatest effect in minimizing the deleterious effects of the coverage gaps is the addition of yaw, i.e., crabbing the sensor with respect to the direction of flight.

Acknowledgments:
The authors thank the Andrew Ortega and Chad Tripp of the University of Florida UAS Research Program for their work in conceptualizing, integrating, and constructing the UAS and lidar mapping payload used to collect field data used to verify the software's and equations' results. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Characterization Equations
Appendix A.1. Point Density Probability Density Function The first assumption that is presented is the angle of the laser firing from the lidar unit can be modeled as a uniform random variable θ, i.e., θ ∼ U(−π, π). This is valid since the lasers are being fired at a constant rate from the scanner head which is rotating at a constant angular velocity. The angle θ is defined as rotation about the scanner's z-axis, Z N , as depicted in Figure A1. It is also assumed that the sensor is travelling parallel to the ground along its axis of rotation Z N -i.e., at a constant height above ground-and at a constant velocity v Z . The frequency of the laser firings is designated as the variable l f .
The cumulative distribution function (CDF) of θ is: Note that this function limits θ to the range [−π/2, π/2] because no pulses emitted outside of this range will intersect with the ground below (see Figure A1). The next step is to find the probability density function (PDF) of another random variable , defined as the lateral distance of each laser return (in the direction) from the centerline of the laser return pattern (i.e., the line parallel and directly nadir to the flight line of the scanner). is a function of , and can be defined as where ℎ is the height of the sensor above the ground. The cumulative density function (CDF) of can then be defined as Substituting using Equation (A2) and rearranging (noting that h is always positive) yields ( ) = tan ≤ ℎ .
The next step is to find the probability density function (PDF) of another random variable X, defined as the lateral distance of each laser return (in the X N direction) from the centerline of the laser return pattern (i.e., the line parallel and directly nadir to the flight line of the scanner). X is a function of θ, and can be defined as where h is the height of the sensor above the ground. The cumulative density function (CDF) of X can then be defined as F X (x) = P(X ≤ x).
Substituting using Equation (A2) and rearranging (noting that h is always positive) yields Inverse tangent is an increasing one-to-one function, which allows for further rearranging: Substituting from the CDF of θ expressed in Equation (A2) yields (A6) Finally, the derivative of Equation (A6) with respect to x yields the PDF To convert Equation (A7) to a point density function as a function of x, or lateral distance from the centerline of the laser return pattern, it is multiplied by the laser pulse frequency l f and divided by the forward velocity v z . (Note that l f is divided by 2 because half of the laser firings will not intersect the ground.) Thus the CDF of X can then be defined as F X (x) = P(h tan θ cos α ≤ x). (A20) Substituting and rearranging further yields tan −1 (x/(h cos α)) π + 1 2 . (A21) The derivative yields the PDF d dx F X (x) = f X (x) = h cos α π(h 2 cos 2 α + x 2 ) . (A22)