Lidar Attenuation Through a Physical Model of Grass-Like Vegetation

Autonomous vehicles have become increasingly viable and capable in recent years. Several semi-autonomous vehicle functionalities are already available to the public, such as lane departure warnings, collision avoidance, and adaptive cruise control [1]. Tremendous progress has been made in both hardware and software to enable these developments [2].

However, much of the progress in autonomous vehicles has been for applications on paved roads, where there are markers such as road signs, stop lights, and painted lines, as well as environment controls such as gradient limits and lane width. Off-road vehicles have none of these markers, with the added complexity of mud, sand, steep slopes, rocks, deep water, and so forth. To enable autonomous vehicles to drive off-road, it is of critical importance to understand how the available sensory systems interact with these factors. For an overview of the history and challenges of off-road autonomous vehicles, see Ref. [3].

Lidar is a light-based sensor commonly used to inform autonomy algorithms [4]. In this paper, we focus on lidar effectiveness in vegetated environments, and our results will be relevant to any applications concerned with such interactions.

Grass and grass-like vegetation pose a challenge to lidar systems [3]. First, it can be difficult to discern via lidar between potentially drivable grass-like vegetation and more impassable materials, such as rough terrain. Second, even if accurately identified, it can be difficult to determine if grass-like vegetation obscures obstacles [5]. These difficulties marked performances in the DARPA PerceptOR program. Therein, some tested robots completely avoided passable grass-like patches, which led to suboptimal paths, and none of the teams identified logs hidden in tall grass, which led to mission failures [6].

There are lidar-dependent methods to help mitigate these challenges. Seminal lidar-only techniques include, but are not limited to, the identification of vegetation in 2D point cloud data via single-pixel statistical analysis [5,7,8], which will be further discussed in the next subsection, and the identification of light and heavy vegetation through a geometric analysis of 3D point cloud data [9]. More complex approaches classify vehicle surroundings via lidars fuzed with a suite of other sensors, such as monocular and stereo RGB cameras, thermal cameras, and hyperspectral cameras [10,11]. The approaches noted rely on machine learning to accurately classify sensor measurements (Ref. [9] uses a Gaussian mixed models classifier, Ref. [10] uses a neural network, and Ref. [11] uses deep neural networks). With the maturation of deep learning algorithms, we anticipate that machine learning will become increasingly prevalent in this field. This trend may be heralded by the success of convolutional neural networks (CNNs) for vegetation identification using remote sensors in forestry, agriculture, and environmental conservation [12]. CNNs have also shown potential in helicopter landing zone detection with lidar [15]. Consequently, large volumes of training data will be increasingly demanded [14].

Although powerful for field use, the machine learning methods mentioned, especially the neural network methods, do not advance understanding of the sensor’s interaction with grass itself, particularly for the purpose of in silico experimentation. A detailed understanding of lidar beams’ interaction with grass blades is important for accurate simulation modeling.

One detail that can arise in vegetation settings is the mixed pixel effect. Because lidar beams diverge, single beams will hit multiple blades of grass at different depths, creating a perception of objects that are not there. These are essentially unavoidable when firing lidar through grass-like vegetation at close distances, but when obstacles are at least 200 cm behind the leaves, mixed pixels tend to disappear [5], which is relevant to the experimental setup presented here. The same reference suggests that lidar works to detect objects through vegetation at depths of 10 cm to 7 m, depending on the plant type that the beam is traveling through, and our experiment falls within this range.

Another project specific to advancing lidar-grass interaction understanding is the creation of histograms of distance measurements from within and in front of grass, as well as identifying solid surfaces (e.g., rocks) in a patch of grass [7]. Part of the approach was to compare the histograms to a theoretical distribution and detecting deviations.

Neither [5] nor [7], nor any other source to the authors’ knowledge, attempted to design a full factorial experiment to quantify the effects and interaction of distance, density, and target reflectivity. An understanding of these interactions will help increase simulation accuracy and, in turn, increase the accuracy of the machine learning methods.

Simulation of autonomous vehicles is often used to increase the speed of development and has the added advantage of reducing cost, and it is frequently useful to have simulations be faithful to the physical world. Furthermore, the massive amount of training data for neural networks can be supplemented by synthetically generated data, a practice used in various disciplines, including lidar and autonomous driving [13,15–19]. In order to ensure the equations and algorithms driving the simulation are accurate, it is necessary to strategically test different aspects of the simulation with physical experiments. Once the behavior of the physical sensor is well understood, the simulated sensor can be refined to be more accurate and applicable in a variety of simulated environments. This process falls under the umbrella of verification and validation (commonly referred to as V&V), which may be qualitative or quantitative.

An additional factor at play is that much of the available lidar hardware and software implementations are proprietary, creating opaqueness about their exact performance specifications. Inherent differences in design and manufacturing, combined with rapid technological development, can create discrepancies in unit performance, meriting the need for a suite of standardized experiments to encapsulate important aspects of lidar devices. Nascent industry benchmarking efforts are underway, but much more is needed in this area [20].

A probability density function modeling lidar penetration through grass was presented in Eq. (5) of Ref. [21], and we have integrated this density function to compute probabilities. We have also performed experiments in silico with an in-house simulator [22,23]. The analytical equations predict significantly lower penetration measurements at 1 m through the grass than the simulator (A selection of these comparative measurements are reprinted in Sec. 3.5). Therefore, a physical model was employed to understand discrepancies between the equation and the simulator. There are several potential reasons why the equation and the simulator do not match, and these will be discussed in Sec. 2.1.

In order to generate physical models of grass, we needed to accurately represent its distribution. After a literature search, we were unable to find a compelling model to use for the spatial distribution of grass stems that was based on real data and endorsed by ecologists. One source used Poisson processes as null distributions [24]. Other sources used Poisson processes for a null model for trees, which we considered relevant since we want this work to apply to young saplings as well [25,26]. Furthermore, a significant motivation for this work is the validation of in-house sensor simulators [22]. Exact virtual replication is facilitated with the use of a basic geometric shape distributed in a well-defined pattern. The following research uses a similar stem shape and distribution: the work in Ref. [21] accounts for beam divergence, while Ref. [5] uses a simpler model without beam divergence and includes experimental data, and Ref. [27] uses a uniform distribution even in the context of animal grazing and savanna dynamics, suggesting that the Poisson process—which distributes points uniformly—is applicable in a diversity of real contexts.

We setup an artificial model of grass, placed the lidar on one side of the model, and placed a calibrated, diffuse-reflecting lidar target on the other side. Our response of interest for the statistical model was whether or not a beam fired through the grass hit the target.

The factors we investigate in the model are the stem density (stem count per unit area), the distance from the lidar sensor to the grass, and the reflectivity of the target. We also investigate differences between two lidar devices.

We investigated VLP-16 and HDL-32E lidar units, both manufactured by Velodyne. The information in this section can be found in greater detail in Ref. [28] for the VLP-16 or Ref. [29] for the HDL-32E, but we include a brief summary here for convenience. Lidar works by sending out pulses of laser light and measuring the return time of each beam to calculate the distance traveled. There are various ways to accomplish this, but a common method is taking measurements from a set of lasers that spins while the beams pulse. Each revolution measures several thousand points, with each representing a point at which a beam reflected off an object in the space around the lidar device.

An added complexity is that the beam diverges as it travels. Consequently, the diameter of a single blade of grass may be a fraction as wide as the beamspot, even at lidar-to-grass distances of a few meters, and any given beam may partially reflect off multiple blades at different depths, or off of objects around the grass. Furthermore, small variations in the rate-of-revolution cause different readings as the beams hit various places in the vegetation. As the lidar sensor moves further from the vegetation, the beam diverges more and also returns with less intensity, eventually weakening the reading.

Our primary data come from a HDL-32E lidar unit, model year 2020. We have secondary data from a Velodyne VLP-16 lidar unit, model year 2017. Both units were operating with the most recent firmware available as of July 2021.

Both lidar sensors have strong, last, and dual return modes. If the entire beam spot hits a single target, they will all coincide. The discrepancy arises when a single beam is split across multiple targets at differing distances.

The strong return mode returns the strongest portion of the returned beam. Because the beam gets weaker with distance, this often means it returns the closest object. However, hitting a more reflective or retro-reflective object will often override the distance factor.

The last return mode returns the portion of the returned beam that took the longest to return.

The dual return mode returns both last and strongest, as long as the separation between the two is at least 1 m.

Return mode is relevant because of the beam divergence. The horizontal divergence angles for the HDL-32E and the VLP-16 are reported as 2.79 mrad and 3.0 mrad, respectively. This equates to beamspots of 20 mm and 47 mm at lidar to target distances of 3.5 m and 8.5 m, respectively, for the HDL-32E, and 21 mm and 51 mm for the VLP-16. Meanwhile, the model grass blades have a diameter of 6.35 mm. This means that every beam that hits a blade of grass will split, potentially differentiating strong, and last return modes. Last return would prioritize portions of the beamspot that hit the target (provided they make it sufficiently intact through the patch), and strongest return would be more variable, based on the reflectivity of the grass, the width of the grass, how many blades were hit, and the reflectivity of the target.

Preliminary testing showed that last return mode generally gave a higher percentage of beams hitting the target beyond the grass than strongest return mode. In the setting of measuring a target through vegetation, we used the more relevant last return mode.

Although Eq. (5) of Ref. [21] assumes first return mode, the Velodyne lidar units available to the experimenters did not have a first return mode. Although first and last returns frequently coincide, some approximation is introduced in this comparison. Equation (5) of Ref. [21] also assumes circular beam spots and Velodyne lidar has rectangular beam spots, which could introduce a small amount of noise. Reflectivity is an additional factor the model does not account for, and the experimental setup contains the two main reflectivities of the wooden dowels and the calibrated target (a diffuse reflector). Additionally, Velodyne’s signal processing algorithms of strongest and last return are proprietary and unknown, complicating comparison to a theoretical model and further suggesting the need for benchmarking experiments. Nonetheless, the authors consider comparison to [21] to be insightful, considering that there are so few theoretical models available with such detail as Goodin’s and given the limitations of the Velodyne lidar units.

The beams of the HDL-32E fire between −30.67 deg and 10.67 deg away from horizontal, with step sizes of 1.67 deg From these, we selected four that were on-target at 10 m away. The beams of the VLP-16 fire between −15 deg and 15 deg away from horizontal, with step sizes of 2 deg. From these, we selected two that were on-target at 10 m away.

We built 12 patches of model grass measuring 1 m by 0.5 m using wooden dowels that were 6.35 mm (0.25 in.) in diameter and 1.22 m (48 in.) tall. We placed the dowels into extruded acrylic sheets that were drilled by a Thunder Laser Nova 51. The sheets measured 6.35 mm by 59.69 cm by 120.65 cm (0.25 in. by 23.5 in. by 47.5 in.).

The coordinates of the centers of stems were distributed in R by a Poisson process [30], conditional upon sufficient separation distance (details below). Although the dowels were labeled as 6.35 mm (0.25 in.), we drilled 6.096 mm (0.24 in.) holes because the cutting laser added additional width to the hole. We also enforced a 2.54 mm (0.1 in.) separation between the perimeters of the circles for structural integrity, so in total we required the stem centers to be at least 8.636 mm (0.34 in.) apart.

The stem centers were generated at mean densities of 40 and 120 stems per square meter. Specifically, for each patch i of density λ, we generated the number of stems n_i ∼ Poisson (λ/2). The λ is halved because we are working with a patch of 0.5 m². We then picked the (x_j, y_j) coordinates of the stem centers with x_j ∼ uniform (0, 1) and y_j ∼ uniform (0, 0.5), since the patches are 0.5 m by 1 m. Each time we generated a point (x_j, y_j), we rejected it if it was within a Euclidean distance of 8.636 mm from any previously added point. The edges of the dowels were allowed to cross over the edge of the patch, but the centers of the dowels had to stay in the boundary of the patch. We continued rejection sampling until we had n_i points for each patch i.

We set the lidar 3 m and 8 m from the front of the grass patch, with the light then passing through 1 m of simulated grass and then traveling another meter before hitting the calibrated target, as shown in Fig. 3. Although the functional range of lidar systems is significantly higher than 10 m, granular resolution with multiple beams was desired. The data at further distances would be less rich, since the geometry of the beams would preclude multiple beams sweeping through the patch simultaneously, and at a further distance the firing rate would give only a few data points per sweep through the patch. At great enough distances, the beams would regularly miss the patch on a revolution altogether. Lastly, the desired functional use case of a vehicle determining whether it can drive through vegetation preferred a closer, data-rich distance, as opposed to being able to make the determination with less accuracy from further away. At the time of the experiment little was known about the details of lidar-grass interactions, so close distances were used to establish a baseline.

A photograph of the HDL-32E experimental setup is shown in Fig. 1. A photograph of the secondary experiment with the VLP-16 is shown in Fig. 2. The patches, distances, and targets were identical in both setups, but the HDL-32E measurements were taken with metal framing in place in addition to marked measurements on the floor to make the entire structure more exact. Precise markings were made on the floor for both the VLP-16 and HDL-32E data. Measurements were taken indoors with the windows covered and lights off.

We used a 2³ factorial designs with mixed models to account for repeated measures on patches. We investigate the probability of hitting a target through grass as it varies in response to three factors: distance from lidar to target, grass density, and target reflectivity. There are 12 patches, with six patches of each density. Each patch was tested across a 2 × 2 combination of distances and targets, for a total of 48 experimental groups. The VLP-16 data have 24 groups, since only the 5 m measurements were available. At each combination of factors, we take 100 measurements. Each measurement corresponds to a frame, and we measure for 10 s at ten frames per second.

Within a single frame of the lidar, we count the number of beams fired at the target, n_f, and the number of beams that reached the target, n_t. To model this logistically we require a binary response. For a frame with n_f and n_t counts with a treatment tuple of density λ₁, distance d₁, and target t₁, we create a data frame with n_f − n_t rows having a response of 0 (representing failure to reach the target) and n_t rows having a response of 1 (representing success), so each row of the data represents one pulse of one beam.

Note that n_f can vary slightly based on where the beam spot first hits the edge of the target during each revolution. For example, when the target is 5 m from the lidar, the back of the 0.5 m patch (4 m from the lidar) spans about 7.15 deg. For the HDL-32E at 600 rpm the beams are spaced at 0.166 deg (see Ref. [29], p. 55), so the beam could hit about 43 times. For the VLP-16, at 600 rpm the beams are spaced at 0.199 deg (see Ref. [28], p. 51), so the beam could hit about 36 times. If the beam starts slightly past the edge, then it is possible to lose some hits, but such variation is low at the distances in this experiment [31].

We outline the hierarchical randomization as follows:

1. Generate six patches of low density (labeled as 1–6) and six patches of high density (labeled as 7–12).

2. Randomly select one of the 12 patches without replacement.

3. Randomly select a target (low or high reflectivity).

4. Randomly select a distance (near or far).

5. Measure for 10 s (100 frames) at the given patch-target-distance combination.

6. Repeat step 5 from the distance not yet measured.

7. Repeat steps 4–6 with the target not yet measured.

8. Perform steps 2–7 repeatedly until all 12 patches have been measured.

Complete randomization was not feasible due to the intensive nature of constructing the patches and moving the equipment.

By measuring at multiple distances and targets through each patch we maximize our use of the patches, which are difficult to align to the lidar, time-consuming to place dowels in, and require the most expensive consumable materials (particularly acrylic, which we bought to gain accuracy with the laser cutter). Additionally, patches of the same density may allow light through at different rates, even under the same experimental conditions. Furthermore, we break the standard regression assumptions by losing statistical independence if we repeatedly measure the same patch. This suggests the use of a model that rigorously accounts for the dependence introduced by repeated use of the patches.

Each physical grass model has inherent features specific to the particular arrangement of its stems that can increase or decrease light penetration. These features are not connected to the experimental conditions of interest, and can be accounted for with a random intercept. In the fitted regression model, all output values at a particular patch of grass are adjusted by the intercept value corresponding to that particular patch. These random intercepts come from viewing the model patches of grass as members of a larger population of all possible patches of grass generated via the process described in Sec. 2.2.

We took a preliminary set of measurements with the VLP-16 before refining the process and taking a second set of measurements with the HDL-32E. The measurements for the second set were more precise due to the use of better equipment for the alignment of the lidar device with the patch and the target. The targets, distances, and patches of grass were identical to the first set of measurements taken with the VLP-16. There are two main reasons we cannot fold a comparison of the two lidar units into a single regression model. The VLP-16 measurements could only be processed properly at the closer distance of 5 m, and the equipment additions for the HDL-32E were substantial. However, the VLP-16 measurements are complete enough to provide valuable insights and direction for further research. We will analyze the HDL-32E and the VLP-16 in their own respective models and informally compare coefficients.

Reading in the raw data was done in Veloview 4.1.3 for the HDL-32E and Veloview 3.5.0 for the VLP-16 [34,35]. Defining the field of view was done in python via Spyder using pandas and numpy [36–40].

Retroreflectors, which are easy to detect on lidar due to high intensity readings, were placed at the edge of the acrylic piece, flush with the front of the patch, and 4.85 cm (1.91 in.) from the front corners of the patch, to define a known marking point to be able to define the relevant field of view of the lidar. Using the retroreflectors and the guidance rails we had set up, for each set of 100 frames we restricted to all beams between the corners of the patch closest to the target, to ensure every beam had traveled the full depth of the grass. This gave us Θ_i, the field of view for experimental block i. The shared field of view, θ_FOV, across all trials with a shared lidar-to-grass distance is defined as θ_FOV = min_i(Θ_i). The physical targets we used were wider than 0.5 m, so every beam within the described field of view was capable of hitting the target (see Fig. 3).

The final data we used for the model contained the treatment information (density, patch, target, distance), as well as the frame index (from 1 to 100) within that treatment set, and the outcome variable hit encoded as 1 for success (the beam hit the target) and 0 for failure (the beam was stopped in the grass model). Density, target, and distance were coded as factors with 0 and 1 representing low and high, respectively.

In this section, we show example grass patch simulations and exploratory data analysis plots. We also investigate how the number of beams fired and hit-rate success of those beams varied over factor combinations and individual frames. The work throughout the remainder of the paper was done in R via RStudio using tidyverse and lme4 [41–43].

The random orientation of the grass stems had a significant impact on the outcome. In the low-density patches (labeled as 1–6) at 40 stems per square meter, patch 1 had the fewest stems, the most skewed stem arrangement, and unsurprisingly the highest proportion of beams hitting the target, while patch 4 had the most stems and the lowest proportion. However, in the high-density patches (labeled as 7–12) at 120 stems per square meter, patch 7 allowed fewer beams through than any other patch (nearly 0% on average), while patch 11 allowed more than any other high-density patch, and almost as many as patch 4, even though the numbers of stems between patches 7 and 11 were scarcely different (see Fig. 5). This supported our assumption that patch-specific random effects are significant factors in attenuation rate. The four mentioned patches are shown in Fig. 4.

There was high variation in hit rates across low-density patches of grass, but low variation in hit rates across high-density patches. In Fig. 5, we have averaged hit rates over all 100 frames of each patch-target-distance treatment set for the HDL-32E.

Several things become clear from Fig. 5. Every low-density patch had a higher hit rate than every high-density patch. The wide variation is evident across patches, particularly low-density patches. All low-density patches except for 4 showed hit rate decreasing as distance increased, but the high-density patches were less consistent. The effect of target reflectivity is unclear. Table 1 summarizes the hit rates averaged across all patches.

Despite sharing a common angular field of view, the quantity n_f of beams returning to the lidar (the number of hits on-target and misses hitting the grass combined) in a frame varies even under perfect conditions because the pulses do not hit the precise leading edge of the target exactly, and the firing rate is not in perfect sync with the RPMs. The most predictable factor in a varying beam count is distance. Being twice as far away at the same RPM and firing rate will cut the number of beams fired in a fixed physical window in half. Therefore, in this section, we refer to “adjusted beams fired” as the exact number of beams for 10 m measurements and the number of beams divided by two for 5 m measurements. A beam only counts as “fired” if the lidar unit registers the beam on its return. In Fig. 6, we can see variation in the adjusted number of beams detected across experimental groups. The more stems are in the patch, the fewer beams we measure returning to the lidar at all (whether on-target or not).

Importantly, we did not see significant trends in frames 1 through 100 in either proportion of beams hitting the target or in total number of beams detected. Figure 7 shows this for four out of of the 48 experimental groups, with unadjusted beam counts per laser. This observation supported our assumption that we did not need to add levels of hierarchical structure in the mixed model for the frame or the laser.

The data discussed in the previous sections were from the HDL-32E, with a more stable setup alignment and consequently twice the number of experimental groups (both distance 0 and distance 1). However, a preliminary experiment with the VLP-16 still gave us useful insights.

The VLP-16 analog of the plots shown in Fig. 7 show similar patterns. There was variation across experimental blocks, but the across-frame, across-laser counts of beams fired, and the proportion of beams on-target were consistent with our assumption to exclude frame and laser as additional layers of random effect.

The analog of the plots shown in Figs. 5 and 6 are presented in Figs. 8 and 9. The case for more stems resulting in a lower total count of lidar beams recorded was slightly more convincing with the VLP-16 data (see Fig. 9), but because of higher variation in beam counts across experimental groups, we were hesitant to make any firm conclusions.

The VLP-16 hit rates showed a trend with target reflectivity, while the HDL-32E showed more mixed results. Compare Figs. 5 and 8. The hit rates were overall higher with the VLP-16 than with the HDL-32E, which in turn affects the coefficient estimates, since the logistic function becomes less stable as event probability approaches zero. The summary of hit rates averaged across all patches for the VLP-16 can be seen in Table 2.

Coefficient estimates rounded to three decimal places are in Table 4 (with estimates and standard errors computed from the glmer command). These coefficients are on the scale of log-odds. The bottom row is the estimate of the standard deviation of the random effects, i.e., the patches. The last two columns are computed from the 200 bootstrap replicates.

The fixed effect intercept represents the hit rate for the setup with the low-density patch on a 10% reflective target from 5 m away. Its coefficient estimate is −0.213 and its 95% compatibility interval ranges from −0.864 to 0.509, log odds. These values equate to pass-through ratios of 0.447 and 0.297 to 0.625, respectively.

The other probabilities are computed by taking the inverse logit of the sum of the relevant coefficients. Table 3 shows the resulting values.

The patch intercept shown in the last row of Table 4 represents the estimate of deviation of passes per patch density. The estimated value of 0.851 suggests that specific arrangements of grass are significantly more transparent to lidar beams than others, even when controlling other variables. A similar conclusion holds in Table 5.

The VLP-16 bootstrap analysis over 200 replicates is shown in Table 5. Note the similarities between the models in density coefficients, but the differences in target coefficients. Compare Figs. 5 and 8 to visualize the difference in target effects. Furthermore, Fig. 8 shows the difference in slopes between high and low density patches as target varies. On average, higher density patches are more greatly affected by increasing target reflectivity.

The model was analyzed using standard diagnostic tools such as residual analysis to make sure that it was reasonable.

The existence of a two-way interaction shows that the effect of a factor differs across the level of another factor—e.g., in Fig. 8, the slope for target reflectivity in high-density patch runs is higher than the slope in low-density patch runs. This means that, for the VLP-16, we saw target reflectivity affecting hit rate more in high-density patches than low-density patches.

A three-way interaction signifies that a two-way interaction differs across the level of a third factor. This is shown in Fig. 10. In low-density patches on the left panel, increasing target reflectivity decreased hit rate at both distances. However, in high-density targets, increasing target reflectivity increased hit rate for both distances, and moreover these slopes are not equal in magnitude.

The cause of this complex interplay between factors is unknown, but it is hypothesized that beam refraction may be involved.

Following standard practice, lower-order interaction terms are kept in the model if a higher-order interaction term is significant, regardless of whether that lower-order interaction term is significant or not. This implies the important fact that for prediction purposes, the three-factor logistic model with all eight parameters $β0,…,β7$ is no more parsimonious than using the eight cells means directly (Table 3). A similar result holds for the four-parameter VLP-16 model and Table 6. This is significant for any future attempts to model lidar measurements of grass: testing for complex interactions should be included.

We simulated the experiment with in-house lidar simulation software [23]. This study used a similar point process to that in Sec. 2.2 but omitted a rejection scheme. The results are also listed under the simulation predictions column in Table 7. Neither the analytical model nor the simulation applied target reflectivity in computation of their predictions.

We found that the precise layout of stems is an important contributing factor, as evidenced by the substantial standard deviation of the random patch effect estimate (see the last row of Tables 4 and 5).

Findings regarding distance are limited to the HDL-32E; their distance seemed to have a small but measurable effect.

The compatibility intervals for the interaction coefficients are wide, particularly for density with distance, and further investigation is needed here. The significant second order interaction term is notable. Looking at Table 1, we see that within low-density patches, hit rates changed similarly, by approximately 15% points, with distance across targets. However, within high-density patches, low-reflectivity targets saw a slight decrease while high-reflectivity targets increased by a smaller amount. This distance-target interaction changing across levels of density is a sufficient, but not necessary, condition for the statistical significance of the three-way interaction.

Both target reflectivity and overall pass-through rate varied significantly by lidar unit. While the VLP-16 showed reflectivity to have a measurable impact on hit probability, the HDL-32E showed a much smaller effect. Additionally, most patches had a lower proportion of beams passing through with the HDL-32E than the VLP-16. There are a number of reasons this could be the case. Perhaps the HDL-32E is better at seeing small objects and is more robust to reflectivity. However, other factors could have been at play, from lens condition to motor age to firmware.

An important point is how different individual lidar units can be. Consistent measurements with one unit might not transfer to another unit. This experiment showed notable differences between two Velodyne-manufactured lidar units, a VLP-16 and an HDL-32E. Our preliminary experiments suggested that two different VLP-16 units can show markedly different data as well. Importantly, the distance measurements were much more consistent across different lidar units than some of the other phenomena we measured. We were not specifically investigating distance measurements themselves, but from what we saw in data processing, both devices were accurate in that category, both within each device and across devices.

The analytical and simulation predictions qualitatively agree with the basic relationship between stem density and pass-through rate. The predictions also demonstrate a qualitative interaction between patch density and distance. However, the predicted pass-through rates are consistently and considerably higher than those observed in our experiment. Moreover, they fall outside of the compatibility intervals for nearly every test configuration. Three exceptions are in the last two rows of Table 3 and in the last row of Table 6, but these rows also have exceptionally wide compatibility interval ranges. These contrasts suggest that there are deeper concepts at play in our experimental process, the equation, the simulation, or all three. Of note is that both the simulation and the analytical equation predicted an increased hit rate with greater distance (Table 7), while the experimental data show a marked decrease in hit rate from a greater distance (Table 1), with only a slight exception in one case. However, the simulator was only available for the VLP-16, and the distance measurements in this experiment were only available for the HDL-32E, so there could be some discrepancy there based on divergence, beam angle, and other factors. The authors believe reflectivity is an important contributor as well in the real data.

Weather conditions are of great interest to anyone using lidar in outdoor environments, and the additional complication of fog, rain, and snow with the lidar-grass interaction would be an important direction of future investigation.

Further work is needed to understand the discrepancy mentioned previously between the discussed analytical equation, ERDC’s in-house simulator, and our results here. The magnitude between the results from the closed-form expression and the experimental data suggest missing information, and the authors believe reflectance could be important. Reflectance of machined wood versus various types of vegetation was not able to be quantified here. Furthermore, as described in Sec. 2.1, it is known that the analytical equation bears hypotheses that could not be met in the experimental setup, which does explain some of the discrepancies.

Instead of using a binary indicator of whether a beam passed through a patch, one could use the entire distribution of distances as the response variable. In theory, with an infinitely thin beam, measured distances in a patch of grass should follow an exponential distribution [8], but the beams in lidar diverge so quickly on most relevant scales that that theory falls apart. Another option would be to devise a sensible metric for comparing two-point clouds, using two- or three-dimensional space instead of collapsing to a single distance distribution. This could help with grass identification and give some idea about its density, even from a distance. A binary response or a distributional response each might make sense depending on context. One approach would be to use a binary response when looking through grass, and a distributional one when deciding whether to drive through it.

The target 10 m away receives half the number of beams fired at the target as the target 5 m away. Beta regression using a response of per-frame hit rates in (0, 1) would balance the size of the experimental groups, if such balance is desired.

Although we do not look at it here, stem diameter is an interesting additional factor because of how it interacts with beam diameter. A single stem of near-zero width (a strand of silk, for example) would have virtually no effect, but a threshold point must exist where the beam would begin detecting a sufficiently wide stem. What this threshold point is and how it interacts with the reflectivity of the stem is unknown. None of this takes into account leaves and variable blade thickness. Distance interacts as well, since the beam gets wider as it travels, thus making a given stem thinner relative to the beam, especially from distances still feasible in practice and measurable with a Velodyne unit. Light intensity also drops off with distance, however, further affecting how the lidar reads each stem.

A notable source of variation is the lidar devices themselves. Further research is needed to discover the level of variation across lidar units within a single company, model, or model year, while controlling external factors.

We sought to understand how grass density, target reflectivity, and lidar distance affected whether lidar beams penetrated a grass-like vegetation model. We found that the particular layout of stems can have a significant effect on the outcome, as does the individual lidar unit used to make the measurements. It is important to note that distance and density interact in a complicated way and in the case of the HDL-32E there was significant three-way interaction. Beam divergence, stem diameter, light refraction, and other properties could be combining to exhibit this interaction. Based on this analysis, models, simulations, and algorithms related to lidar and vegetation should include higher-order interaction terms or they could be suffering in predictive power.

We were surprised to see that target reflectivity alone affected the HDL-32E so inconsistently, with many patches of grass experiencing fewer registered beams on the higher-reflectivity target, while the VLP-16 behaved in a very straightforward and expected manner, with every patch of grass showing an increase in on-target hit rate as target reflectivity increased. The authors believe this discrepancy between units motivates the need for standardized characterization experiments with different lidar devices, model years, and firmware editions, since finely-tuned algorithms could give significantly different results based on what could be perceived as a minor software or hardware upgrade to an autonomous vehicle. These experiments should seek to quantify how individual lidar units react to changes in target reflectivity, grass density, distance, return mode, stem diameter, or any other parsimonious set of variables that are suitable to encapsulate the differences between lidar units that are most meaningful for the application at hand.

The mixed models approach gave us a more in-depth view than if we had simply taken averages over multiple frames at a time, all while avoiding pseudoreplication, and this approach is generalizable to other situations where lidar units are measuring the same objects multiple times under different scenarios. With a sufficient number of grass patches, generalized estimating equations could be used in a similar way, depending on project goals. If within-frame or within-lidar measurements show higher correlation than across-frame or across-lidar measurements, respectively, then a hierarchical mixed model could be used with nested random effects.

The authors would like to thank the NSF and ORISE for funding, Ben Webb and Siege Robotics for their time, machinery, and expertise, and Ann Marie Weideman for her insights on statistical modeling.

Mr. Petty was supported by an appointment with the National Science Foundation (NSF) Mathematical Sciences Graduate Internship (MSGI) Program sponsored by the NSF Division of Mathematical Sciences. This program is administered by the Oak Ridge Institute for Science and Education (ORISE) through an interagency agreement between the U.S. Department of Energy (DOE) and NSF. ORISE is managed for DOE by ORAU. All opinions expressed in this paper are the author’s and do not necessarily reflect the policies and views of NSF, ORAU/ORISE, or DOE. This research was also supported in part by the Engineer Research and Development Center (ERDC) within the U.S. Army Corps of Engineers.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.