A Dynamic Blind Zone Simulation and Analysis Model for Roadside Light Detection and Ranging Sensor Deployment Considering the Full Roadway Terrain and Vehicle Dynamics

Abstract

Blind zones in light detection and ranging (LiDAR) sensors arise from their limited physical field of view and obstructions caused by static infrastructure or moving objects. Although originally intended for vehicle-based applications, LiDAR sensors are now increasingly deployed in roadside infrastructure for traffic monitoring and connected and automated vehicle safety and mobility applications. However, there is a dearth of robust tools for analyzing their detection range, resolution, and other characteristics in such settings. This study introduces a three-dimensional (3D) blind zone simulation model for analyzing the detection characteristics of roadside LiDAR sensor deployment. The model replicates the impact of static infrastructure conditions and dynamic blind zones during live traffic. Initially, a real-world digital surface model (DSM) captures 3D data of road surfaces and obstructing infrastructure objects. Optical geometry models then assess blind zone severity across various roadway areas. Subsequent 3D vehicle shape and dynamic simulations evaluate blind zone distributions under typical traffic conditions. The model’s effectiveness is validated using field 3D point cloud data and vehicle detection data collected from a roadside LiDAR site on Route 18 in New Brunswick, NJ. Evaluation results demonstrate the model’s capability in analyzing complex static and dynamic blind zone distributions, offering insights for optimizing LiDAR sensor location, height, tilting angle, and manufacturer configuration parameters to minimize sensing blind zones. For code availability, see https://github.com/rutgerstslab/LiDAR-Coverage-Analysis.

Keywords

data and data science information systems and technology sensor data analytics visualization in transportation 3D data visualization intelligent transportation systems connected

Light detection and ranging (LiDAR) sensors utilize laser beams to determine the distance and three-dimensional (3D) location and shape of objects. LiDAR sensors can perform vehicle tracking ( 1 – 4 ), accurate speed estimation ( 1 ), vehicle classification ( 2 , 5 ), and pedestrian detection ( 2 , 6 – 8 ). LiDAR has made significant advancements in functionalities, range, resolution, form factors, and affordability, driven by the growing demand from the self-driving industry. Although LiDAR sensors have been extensively used in on-board advanced driver assistance system and automated driving system platforms, roadside LiDAR technology has emerged as a promising sensor support solution for enhancing connected and automated vehicle (CAV) applications. Several studies have explored the integration of roadside LiDAR into transportation applications, including aiding traffic operations ( 1 ), enhancing transportation safety ( 6 , 8 ), and supporting CAV applications with high-resolution trajectory data ( 6 – 8 ). Overall, roadside LiDAR sensors demonstrate promising potential to enhance real-time situational awareness and decision-making.

High-resolution roadside sensing from LiDAR, computer vision, and other sensors provides excellent opportunities to overcome the biggest weakness of prevailing connected vehicle (CV) system deployment. The current low penetration rate of CV on-board units (OBUs) in general traffic and regulatory actions on shrinking the designated dedicated short-range communication bandwidth by the Federal Communications Commission ( 9 ) means that traditional CAV system designs relying on in-vehicle OBUs may not be feasible in the near term. The high-resolution sensors can track the movements of vehicles, pedestrians, motorcycles, bicycles, and other objects and create basic safety messages and personal safety messages on their behalf without the need for on-board units—that is, full participation of all road users.

Roadside LiDAR sensing operates differently from on-board LiDAR sensing for self-driving. Roadside LiDAR sensors are installed at a fixed roadside location, making their coverage stay at the same location with a fixed infrastructure background. Such coverage ( 10 ) differs from on-board LiDAR sensing, which experiences more dynamic backgrounds because of vehicle movement. This distinction simplifies the reduction and removal of roadside LiDAR background data ( 11 ). Additionally, the fixed installation eliminates the need for extensive real-time simultaneous localization and mapping operations. Compared with on-board LiDAR, roadside LiDAR has a wide variety of vehicles and full 24 × 7 infrastructure coverage as a trade-off for limited geographical coverage. Recent research in transportation technology has significantly improved traffic safety and operational efficiency with roadside LiDAR. Zhao et al. ( 2 ) emphasize the need for real-time cooperative systems with infrastructure-mounted LiDAR sensors to enhance situational awareness for all road users, particularly in scenarios prone to occlusion, where their proposed solution involves deploying multiple LiDAR sensors to extend detection range and mitigate obstructions. Kloeker et al. ( 12 ) introduce a framework for the quality evaluation of smart roadside infrastructure sensors, offering insights into sensor setup configurations for automated driving applications and smart cities. Du et al. ( 13 ) present a life-cycle monitoring framework that employs high-precision positioning and perception data from CAVs to assess the accuracy of roadside sensor data throughout their operational life cycle within cooperative vehicle infrastructure systems.

Despite the promising potential of roadside LiDARs in enabling CAV applications, the design and planning tools or methods for roadside LiDAR installation still need to be improved. Like other line-of-sight (LOS)-based sensors, LiDAR cannot penetrate objects and is therefore affected by the blind zones caused by infrastructure objects and vehicle occlusions. It blocks part of the 3D space whenever something is detected, creating a blind zone that the sensor cannot “see.” A recent study by the research team developed a preliminary two-dimensional (2D) blind zone simulation model for roadside LiDAR design and planning ( 14 ). The simulation model takes into account the 2D geometry of the roadway and uses basic optical geometric models to simulate the blind zone distribution for a roadside installation on a leveled, flat, and clear infrastructure background but with dynamic vehicle trajectory inputs such as Next Generation SIMulation (NGSIM) datasets. The model can assess the impact of basic roadside LiDAR configuration parameters and provide insights into how LiDAR blind zones are distributed within live traffic. However, the model does not consider the 3D terrain and complex infrastructure conditions in the real world and assumes that all vehicles are traveling on a straight road through a rectangular grid system without considering vehicle orientation or attitude because of roadways’ horizontal and vertical curves.

In this paper, we proposed several key features to create a comprehensive LiDAR blind zone simulation model suitable for real-world LiDAR deployment design, planning, and blind zone analysis. First, a digital surface model (DSM) is proposed to reconstruct real-world 3D roadway and roadside infrastructure. Second, the simulated vehicles can now have additional degrees of freedom, including orientation and attitude. Third, flexible 3D shapes of vehicles can now be considered, and their full interaction with 3D road surfaces can be modeled. Last but not least, a sensing overlapping analysis module is introduced to allow multi-sensor configuration. We also designed new computational-efficient occlusion determination methods to enable the above new features. With these new capabilities, the simulated blind zone model is validated with field LiDAR data at the grid level and is used for roadside site planning and design in the DataCity Smart Mobility Testing Ground project. The dynamic DSM-based 3D reconstruction is efficient in both computation and transmission and is ready for large-scale real-time digital twin visualizations. The LOS computing method fully utilizes the efficiency of DSM and only require 2D LOS determination plus one height comparison for z-axis, instead of 3D frustum. Compared with the authors’ prior work, the proposed LiDAR blind zone analysis model mainly has the following improvements. (1) Real-world 3D road surface is sampled in the form of DSM for sensor location evaluation, including multi-sensor deployment. (2) Preprocessing is designed to smooth the DSM for road user placement and movement simulation. (3) Real-world trajectories are mapped and embedded into the 3D road surface with actual position and 3D orientation to reflect the impact of complex behaviors such as lane changing, making turns, and climbing slopes. (4) The redesigned the parallel-computable per grid occlusion determination has stable performance against large-scale trajectories. (5) Different LiDAR sensors are compared and valuable insights are generated for both sensor manufacturing and deployment.

Literature Review

The existing literature has explored the impact of different vehicle-mounted on-board LiDAR sensor configurations on 3D object detection and tracking in self-driving applications. The impact of different configuration parameters and design elements have been considered, including physical design perspectives ( 15 ), information-theoretic surrogate metrics ( 15 ), and Bayesian theory-based evaluations ( 16 ) to mixed-integer linear programming ( 17 ) and bio-inspired cost functions ( 18 ). Hu et al. ( 15 ) introduced an information-theoretic surrogate metric, called Surrogate Maximum Information Gain (S-MIG), to quantitatively and efficiently evaluate different on-board LiDAR configurations and their impact on detection performance, including center (Argo AI), line (Ford), pyramid (Cruise and Pony AI), and trapezoid (Toyota). Their experiment demonstrates a correlation between S-MIG and detection performance. However, S-MIG is more efficient, as it does not require performing the actual detection. Ma et al. ( 16 ) proposed a perception entropy-based metric to model the perception potential of candidate on-board sensor configurations based on Bayesian possibility theory. Dybedal and Hovland ( 17 ) proposed a method and accelerated it with a graphic processing unit (GPU) ( 19 ) to handle cases where several sensors must cover certain sub-volumes (redundancy). The objective function includes maximizing the number of cubes covered by the sensors and minimizing the cost of the sensors needed to cover the entire area. Liu et al. ( 18 ) used a bio-inspired measure called the volume-to-surface-area ratio as a cost function representing the size of non-detectable subspaces of a given LiDAR configuration to address the optimal LiDAR configuration problem for autonomous vehicles. Kim et al. present multi-sensor placement methods for multiple LiDAR sensors on an autonomous vehicle ( 20 ) and roadside infrastructure in urban environments ( 21 ). For autonomous vehicle LiDAR configuration ( 20 ), the objective is to find optimal LiDAR positions that reduce the dead zone and improve the point cloud resolution. The problem is formulated as an optimization problem using the LiDAR occupancy grid and solved by a genetic algorithm. Experimental results using commercial LiDAR devices demonstrate the effectiveness of the proposed method in improving LiDAR placement. For roadside infrastructure ( 21 ), they analyze the dead zone caused by buildings near the intersections with a voxel-based LiDAR occupancy space and LiDAR occupancy voxel grid. The evaluation is done with different simulation scenarios involving both pedestrians and vehicles. The authors also mention future extensions of the method to incorporate sensor fusion with multiple LiDARs, radars, and cameras into the autonomous vehicle system.

Point cloud simulation has recently been introduced to evaluate and optimize different roadside LiDAR sensor placement and configuration scenarios. Ge et al. ( 14 ) proposed a 3D simulation model for roadside LiDAR sensor blind zones. The model takes different location, height, and tilting scenarios of the LiDAR sensors and field vehicle trajectory and size data, such as NGSIM data, and can run LOS simulation for all LiDAR beams to evaluate the spatial and temporal distribution of blind zone areas where LiDAR cannot penetrate in live traffic. However, the model cannot take complex road geometries with curves, 3D terrain changes, and complex roadside infrastructure. The LOS assessment can be quite computationally extensive for large-scale simulation. Cai et al. ( 22 ) proposed an optimization framework for infrastructure LiDAR sensor placement using LiDAR point cloud simulation within the CARLA (CAR Learn to Act) simulator environment. Their approach incorporated various simulation effects, such as motion distortion and ghosting object effect. Although the simulation did not account for different weather conditions, it included trajectories from OPV2V ( 23 ). The optimization is based on a performance metric, the laser density on the infrastructure through predefined LiDAR occupancy boards ( 24 ). The existing simulation-based roadside LiDAR sensor performance evaluation models are limited by their reliance on preexisting models or their inability to accommodate complex road geometries and infrastructure. Gouda et al. ( 25 – 27 ) proposed the use of octree-based point cloud simulation for infrastructure LOS analysis and compared it with traditional voxel-based simulations ( 28 ). The octree-based method surpasses traditional voxel-raycast approaches in accuracy and practicality for evaluating traffic sign visibility in autonomous vehicle scenarios. Octree-based methods are commonly employed in general 3D data representations, particularly for sparse point cloud data. The advantage of using octrees in LOS simulation lies in their ability to determine LOS more efficiently without requiring detailed examination of every smallest voxel. However, these methods do not fully leverage the specific characteristics of road traffic, leaving room for further efficiency improvements when embedded with dynamic trajectories.

Building on the authors’ earlier work ( 14 ), this paper presents a computationally efficient 3D simulation of LiDAR sensor blind zones that can formulate complex roadway infrastructure, 3D terrain conditions, and detailed 3D dynamic traffic; proposes a more efficient parallel LOS computing method; and adds support for multi-sensor situations. Furthermore, the model is independent of traffic simulation or simulator models, so field vehicle trajectory data can be used to perform sensor placement and configuration optimization with real-world data. Meanwhile, the model still preserves the ability to interact with traffic simulation and simulator models through vehicle trajectory data.

Methodology

Notations and Definitions

$M, N$ : the numbers of horizontal and vertical grids.

$Ω_{M \times N}$ : DSM.

$H_{M \times N}$ : an $M \times N$ array that stores the height information of each grid in the DSM.

$h (m, n)$ : an element in $H_{M \times N}$ .

$H_{V, M \times N} (t)$ : an $M \times N$ array-like $H_{M \times N}$ that stores the height information of road user grids at time $t$ .

$h_{V} (m, n | t)$ : an element in $H_{V, M \times N} (t)$ .

$W_{M \times N \times 3}$ : an $M \times N \times 3$ array that stores the color information of the top point in each grid in red-green-blue (RGB).

$P_{c}$ : coordinate system $c$ .

$α$ , $β$ , $γ$ : rotation angle of yaw (z-axis), pitch (y-axis), and roll (x-axis).

$Ψ_{d, M_{d} \times N_{d}}$ : sensor coverage model (SCM).

$H_{d, M_{d} \times N_{d} \times L_{d}}$ : an $M_{d} \times N_{d} \times L_{d}$ array that stores the height information of each laser in each grid.

$h_{d} (m, n)$ : an element in $H_{d, M_{d} \times N_{d} \times L_{d}}$ .

$Λ_{d, M_{d} \times N_{d} \times L_{d}}$ : an $M_{d} \times N_{d} \times L_{d}$ array that stores each laser’s occupancy and detectability status in each grid.

$h^{*}$ : height of interest (HOI).

$Γ (m, n | t, h^{*})$ : an $M \times N$ Boolean array that stores the blind zone status, determined by HOI $h^{*}$ .

$l (l = 1, \dots, L_{d})$ : laser of LiDAR $d$ .

$t_{0}$ : the time when there is no traffic.

$H_{GT}$ : the lowest point cloud in each grid from actual LiDAR, used as ground truth.

Basic Concept

The proposed method is designed for complex field deployment scenarios in roadside instrumentation. The DSM is used to formulate the 3D terrain and roadway infrastructure ( 30 ). Correspondingly, to simulate vehicles traveling over roads with horizontal and vertical curves, the vehicle attitude (pitch) is added as an additional dimension in vehicle movement simulation. Combined with the vehicle shape and trajectory data, the proposed method will be able to model the 3D dynamics of vehicles in the complex infrastructure conditions in the field. However, the increased complexity caused by the high-resolution DSM formulation poses computational challenges to the proposed LiDAR LOS simulation. We propose new model formulations and parallel computing algorithms for efficient processing of the DSM and LiDAR data for occlusion detection.

DSM

As Figure 1 shows, a DSM is a 3D grid-based surface representation of a real-world terrain and infrastructure condition generated from LiDAR or photogrammetry point cloud. In the field, a road infrastructure’s point cloud data can be collected using vehicle-based LiDAR scanning or drone-based aerial scanning. LiDAR scanning obtains the precise 3D locations of infrastructure points scanned by the laser beams of LiDAR sensors. Drone-based aerial scanning generates a point cloud by matching and triangulating features between overlapping aerial images or video frames.

Figure 1.

Top view of digital surface model (DSM) and sensor coverage model (SCM) grid scheme (Vehicle Icons: 29).

From the availability perspective, the digital elevation model (DEM) is widely available across the country, and the proposed DSM-based method is compatible with DEMs. DEMs can be separated into digital terrain models and DSMs based on the elevation they are collected from. From an efficiency perspective, the proposed method utilized the DSM structure that accurately describes the road surface with dynamic traffic information. Compared with voxel-based methods, the proposed DSM-based method is more efficient because it stores the height information of a 2D map instead of all the 3D information. The occlusion determination is also simplified to 2D plus one height comparison instead of considering the entire frustum. Compared with octree-based methods, generating and updating the DSM with dynamic traffic is less complex, and the DSM-based information can be summarized in quadtree. As the trade-off, the DSM-based method has trouble handling complicated 3D information, which is acceptable for road traffic simulation. We define DSM as $Ω_{M \times N} = {H_{M \times N}, W_{M \times N \times 3}}$ , where $M$ and $N$ are the numbers of horizontal and vertical grids defined for the entire study area. Each grid is a square grid with a side length of $Δ U$ . Uniform square grids can simplify the subsequent calculations compared with rectangular grids. $H_{M \times N}$ stores the height information of each grid, and $W_{M \times N \times 3}$ stores the RGB information of the top color of each grid. The input data for constructing a DSM is a 3D point cloud ${(x_{i}, y_{i}, z_{i} | r_{i}, g_{i}, b_{i}) | i = 1, \dots, P}$ , where $P$ is the total number of 3D points in the point cloud. Each point $i$ contains its 3D coordinates $(x_{i}^{P}, y_{i}^{P}, z_{i}^{P})$ , where $P$ represents the projected or world coordinate system in which the point is collected and its color information $(r_{i}, g_{i}, b_{i})$ .

The DSM is constructed in three steps. First, we project the point cloud ${(x_{i}^{P}, y_{i}^{P}, z_{i}^{P})}$ from its original coordinate system into the world longitude and latitude coordinate system defined for DSM. Second, the point cloud is flattened by grouping the points into the DSM gridding scheme and identifying the height and color of the top points to be inserted into $H_{M \times N}$ and $W_{M \times N \times 3}$ . Third, the hypothetical LiDAR sensors are placed into the DSM for LiDAR blind zone evaluation. The following sections describe the details of these processes.

Multi-Source Point Cloud and Map Coordinate System Matching and Unification

The input data to construct the proposed DSM model includes three different data sources: (1) the satellite imagery from the World Map of ArcGIS in NAD83 (North American Datum of 1983) world coordinates, (2) colored photogrammetry point cloud collected from a drone with centimeter-level accuracy with projected coordinates in meters, and (3) high-resolution millimeter-level uncolored LiDAR point cloud collected from vehicle-based LiDAR scanning. The first data source uses world coordinates in longitudes and latitudes, whereas the other two use their own projected coordinate system, with their origin set at different locations within the study area. Furthermore, the validation and evaluation study for the proposed model also requires the output to be matched with the relative coordinate systems used by LiDAR sensor data centered at the locations of the LiDAR sensors. Since such multi-source scenarios are quite common when working with field geospatial data, a generic method is proposed below to align different coordinate systems with reference points across all models.

The proposed alignment method is an optimization formulation with respect to a series of coordinate system transformation parameters, including rotation angle yaw $α$ , pitch $β$ , roll $γ$ , translation value $Δ x$ , $Δ y$ , $Δ z$ , and squeeze ratio $s$ to be performed with the coordinate system of each point cloud or map data source for alignment. The objective function is to minimize the deviations of the coordinates for the same reference points in different coordinate systems. The optimization will rely on the coordinates of the same ground reference points in all the coordinate systems of the data sources. Through empirical experiments, the minimal number of ground reference points is three. This study uses several reference points evenly distributed within the study area. Assuming the list of points and their coordinates in different coordinate systems are the following:

\begin{matrix} {P_{J}^{P_{c}}} = {(x_{j}^{P_{c}}, y_{j}^{P_{c}}, z_{j}^{P_{c}}) | j = 1, \dots, J} \end{matrix}

(1)

where $J$ is the number of reference points selected and $c$ is the index of each coordinate system where the coordinates are collected for each reference point. Let $c = 0$ represent the base coordinate system that all coordinate systems will project into. Our study uses the NAD83 world coordinate system in ArcGIS World Map.

The problem becomes finding the correct 3D transformation matrix to transform the coordinates of the reference points from a coordinate system other than the target coordinate system ( $c \neq 0$ )—that is, ${x_{j}^{P_{c}}, y_{j}^{P_{c}}, z_{j}^{P_{c}} | j = 1, \dots, J, c \neq 0}$ —to the target coordinate system ( $c = 0$ ), so that transformed coordinates are as close as possible to their coordinates in the target coordinate system. The transformation considering different transformation scenarios can be formulated as follows:

\begin{matrix} [\begin{matrix} {\tilde{x}}_{j}^{P_{c \to 0}} \\ {\tilde{y}}_{j}^{P_{c \to 0}} \\ {\tilde{z}}_{j}^{P_{c \to 0}} \\ 1 \end{matrix}] = S_{4 \times 4} T_{4 \times 4} R_{4 \times 4} [\begin{matrix} x_{j}^{P_{c}} \\ y_{j}^{P_{c}} \\ z_{j}^{P_{c}} \\ 1 \end{matrix}] \end{matrix}

(2)

where $S_{4 \times 4}$ , $T_{4 \times 4}$ , and $R_{4 \times 4}$ are three matrices to account for general scenarios of coordinate system transformations, including scaling, translation, and rotation. The scaling matrix $S$ applies a scaling factor $s$ to all three coordinates:

\begin{matrix} S = [\begin{matrix} s & 0 & 0 & 0 \\ 0 & s & 0 & 0 \\ 0 & 0 & s & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] . \end{matrix}

(3)

The translation matrix accounts for the combined movement along the three axes as follows:

\begin{matrix} T = [\begin{matrix} 1 & 0 & 0 & Δ x \\ 0 & 1 & 0 & Δ y \\ 0 & 0 & 1 & Δ z \\ 0 & 0 & 0 & 1 \end{matrix}] . \end{matrix}

(4)

Finally, the rotation matrix accounts for the combined rotations around all three axes, including yaw (along the z-axis), pitch (along the y-axis), and roll (along the x-axis).

\begin{array}{l} R = R_{z} \cdot R_{y} \cdot R_{x} = [\begin{matrix} \cos α & - \sin α & 0 & 0 \\ \sin α & \cos α & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] \cdot \\ [\begin{matrix} \cos β & 0 & \sin β & 0 \\ 0 & 1 & 0 & 0 \\ - \sin β & 0 & \cos β & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] \cdot [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \cos γ & - \sin γ & 0 \\ 0 & \sin γ & \cos γ & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] . \end{array}

(5)

With the above matrices defined, the generic optimization problem can be formulated as follows:

min ‖ ({\tilde{x}}_{j}^{P_{c \to 0}}, {\tilde{y}}_{j}^{P_{c \to 0}}, {\tilde{z}}_{j}^{P_{c \to 0}}) - (x_{j}^{P_{0}}, y_{j}^{P_{0}}, z_{j}^{P_{0}}) ‖

where

\begin{matrix} [\begin{matrix} {\tilde{x}}_{j}^{P_{c \to 0}} \\ {\tilde{y}}_{j}^{P_{c \to 0}} \\ {\tilde{z}}_{j}^{P_{c \to 0}} \\ 1 \end{matrix}] = S_{4 \times 4} T_{4 \times 4} R_{4 \times 4} [\begin{matrix} x_{j}^{P_{c}} \\ y_{j}^{P_{c}} \\ z_{j}^{P_{c}} \\ 1 \end{matrix}] . \end{matrix}

(6)

Then, with the calculated matrix $S_{4 \times 4} T_{4 \times 4} R_{4 \times 4}$ , we can transform any points and results from their coordinate system $P_{c}$ back to the real-world coordinate system $P_{0}$ . We can also place any real-world coordinates—for example, the longitude and latitude of the potential sensor location—into the DSM $Ω$ with ${(S \cdot T \cdot R)}^{- 1}$ . However, the direct calibration of the nine parameters, including $s, Δ x, Δ y, Δ z, α, β,$ and $γ$ , can be quite cumbersome and requires a significant number of reference points to be selected and matched across different datasets. Several key simplifications are applied based on the characteristics of the datasets in this study:

Objective function simplification. Since the purpose is to group 3D points into 2D grids, we can further simplify the objective function to focus on correcting the longitude and latitude. Since the trajectory placement and blind zone analysis are based on relative elevation, this approach does not affect accuracy for a single DSM analysis. However, when expanding to a larger area, additional effort is required to align different DSMs to the same elevation because of the absence of real-world elevation data:

min ‖ ({\tilde{x}}_{j}^{P_{c \to 0}}, {\tilde{y}}_{j}^{P_{c \to 0}}) - (x_{j}^{P_{0}}, y_{j}^{P_{0}}) ‖

3D rotation simplification. $R_{y} \cdot R_{x}$ projects the 3D points to the top view of the real-world $P_{0}$ , so among any three matched points, { $P_{j}^{P_{0}}, P_{j + 1}^{P_{0}}, P_{j + 2}^{P_{0}}$ } are similar triangles and satisfy:

\begin{matrix} s = \frac{| \vec{P_{i}^{P_{0}} P_{j}^{P_{0}}} |}{| \vec{P_{i}^{P_{c}} P_{j}^{P_{c}}} |}, i \neq j, i, j \in {1, \dots, J} . \end{matrix}

(7)

The Nelder–Mead simplex method ( 31 ) can be used to minimize the variance of $s$ . Then we have the approximated values of $s$ and $R_{y} \cdot R_{x}$ . In most cases, $R_{y} \cdot R_{x}$ can be ignored because the direction of the $z$ -axis is usually correct. Second, applying $R_{z}$ will make the points match the north direction. Therefore, it satisfies:

\begin{matrix} s = \frac{\vec{P_{i}^{P_{0}} P_{j}^{P_{0}}}}{\vec{P_{i}^{P_{c}} P_{j}^{P_{c}}}}, i \neq j, i, j \in {1, \dots, J} . \end{matrix}

(8)

The initial value for $s$ can be set at 1:

\begin{matrix} s_{0} \approx 1 . \end{matrix}

(9)

Translation offsets calculation. With the rotation factors determined through Equations 7 and 8, $Δ x$ , $Δ y$ , $Δ z$ can be calculated by minimizing the $RMS E_{P_{0}, P_{c}}$ . Assuming the 3D model is acceptably accurate in the distance, we can provide the initial values of the translation parameters using the center point coordinate ( $x_{O^{P_{c}}}^{P_{0}}, y_{O^{P_{c}}}^{P_{0}})$ of model $P_{c}$ in the coordinate system of $P_{0}$ :

\begin{matrix} {\begin{matrix} Δ x_{0} = x_{O^{P_{c}}}^{P_{0}} \\ Δ y_{0} = y_{O^{P_{c}}}^{P_{0}} \\ Δ z_{0} = 0 \end{matrix} . \end{matrix}

(10)

The final expanded optimization formulation is as follows:

minimize

\begin{matrix} f (α, β, γ, Δ x, Δ y, Δ z, s) = ‖ ({\tilde{x}}_{j}^{P_{c \to 0}}, {\tilde{y}}_{j}^{P_{c \to 0}}) - (x_{j}^{P_{0}}, y_{j}^{P_{0}}) ‖ \end{matrix}

(11)

subject to

\begin{matrix} {\begin{matrix} 0 \leq α < 2 π \\ - \frac{π}{2} \leq β \leq \frac{π}{2} \\ 0 \leq γ < 2 π \\ s > 0 \end{matrix} \end{matrix}

(12)

where $x_{j}^{P_{c \to 0}}, y_{j}^{P_{c \to 0}}$ are calculated by $S \cdot T \cdot R$ :

\begin{matrix} [\begin{matrix} x_{1}^{P_{s}'} & \dots & x_{j}^{P_{s}'} & \dots & x_{J}^{P_{s}'} \\ y_{1}^{P_{s}'} & \dots & y_{j}^{P_{s}'} & \dots & y_{J}^{P_{s}'} \\ z_{1}^{P_{s}'} & \dots & z_{j}^{P_{s}'} & \dots & z_{J}^{P_{s}'} \\ 1 & 1 & 1 & 1 & 1 \end{matrix}] = S \cdot T \cdot R \cdot \\ [\begin{matrix} x_{1}^{P_{s}} & \dots & x_{j}^{P_{s}} & \dots & x_{J}^{P_{s}} \\ y_{1}^{P_{s}} & \dots & y_{j}^{P_{s}} & \dots & y_{J}^{P_{s}} \\ z_{1}^{P_{s}} & \dots & z_{j}^{P_{s}} & \dots & z_{J}^{P_{s}} \\ 1 & 1 & 1 & 1 & 1 \end{matrix}] . \end{matrix}

(13)

Digital Surface Calculation

With the matched and unified 3D model and point clouds, we can flatten the point clouds into 2D representations along the z-axis with respect to a DSM gridding system at the targeted resolution. As Figure 2 shows, for each grid, we need to determine the RGB color and height of the top points in DSM. The flattening process starts with defining the 2D $M \times N$ grids for DSM $Ω_{M \times N}$ at the same resolution as the LiDAR coverage area. Then $M$ and $N$ are calculated as follows:

\begin{matrix} {\begin{matrix} M = ⌈ \frac{X_{\max} - X_{\min}}{Δ U} ⌉ \\ N = ⌈ \frac{Y_{\max} - Y_{\min}}{Δ U} ⌉ \end{matrix} \end{matrix}

(14)

where $X_{\min}$ , $Y_{\min}$ , $X_{\max}$ , and $Y_{\max}$ are the minimum and maximum values of $x^{P_{c}}$ and $y^{P_{c}}$ , $Δ U$ is the grid size, and $(x^{P_{c}}, y^{P_{c}}, z^{P_{c}})$ is any point in the 3D point cloud model $P_{c}$ .

Figure 2.

Real-world DSM sampling. (a) Grid scheme of experiment site. (b) 3D laser status.

Then, we group all the 3D points in the point cloud by the corresponding grids with parallel computing and then calculate the largest z-value to be set as the height for a specific grid. For each point $(x_{p}^{P_{c}}, y_{p}^{P_{c}}, z_{p}^{P_{c}})$ in the 3D point cloud model $P_{c}$ , their grid position can be determined as follows:

\begin{matrix} {\begin{matrix} m_{p}^{P_{c}} = g_{x} (x_{p}^{P_{c}}) = ⌈ \frac{x_{p}^{P_{c}} - X_{\min}}{Δ U} ⌉ \\ n_{p}^{P_{c}} = g_{y} (y_{p}^{P_{c}}) = ⌈ \frac{y_{p}^{P_{c}} - Y_{\min}}{Δ U} ⌉ \\ k_{p}^{P_{c}} = \frac{z_{p}^{P_{c}}}{Δ Z} \end{matrix} \end{matrix}

(15)

where $m_{p}^{P_{c}}$ and $n_{p}^{P_{c}}$ are the calculated grid indices of point $p$ , $k_{p}^{P_{c}}$ is its scaled z-level at a resolution of $Δ Z$ , and $g_{x} (x)$ and $g_{y} (y)$ are gridding functions as defined above. Then, for every grid $(m, n)$ in DSM $Ω$ , the elements of $H (m, n)$ can be determined by first identifying the point $P_{(m, n)}^{*}$ with the maximal height within the grid $(m, n)$ as follows:

\begin{matrix} h (m, n) = max (k_{i} | \forall i, g_{x} (x_{i}) = m, g_{y} (y_{i}) = n) \end{matrix}

(16)

Then let $i^{*}$ be the index of the point $P_{(m, n)}^{*}$ with the maximal height in the grid $(m, n)$ . The elements of $W (m, n)$ can be determined:

\begin{matrix} W (m, n) = {\begin{matrix} r (m, n) = r_{i^{*}} \\ g (m, n) = g_{i^{*}} \\ b (m, n) = b_{i^{*}} \end{matrix} . \end{matrix}

(17)

Point Cloud Z-Value Smoothing

Real-world point cloud datasets may also contain holes and spikes because of the limitation of 3D reconstruction technologies or the appearance of vehicles and other moving objects shown during scanning. Such issues can worsen with drone-based videogrammetry technologies. Significant holes on the ground can appear because of the lack of matched features for the 3D triangulation of pixels in the image. To solve this issue, we apply the bicubic interpolation ( 32 ) to fill the empty holes on the road surface and surrounding infrastructure (e.g., medians and barriers):

\begin{matrix} \tilde{h} (m, n) = [\begin{matrix} 1 & m & m^{2} & m^{3} \end{matrix}] A_{4 \times 4} {[\begin{matrix} 1 & n & n^{2} & n^{3} \end{matrix}]}^{T} . \end{matrix}

(18)

The coefficient matrix $A_{4 \times 4}$ is calculated as follows:

\begin{matrix} A = [\begin{matrix} a_{00} & a_{01} & a_{02} & a_{03} \\ a_{10} & a_{11} & a_{12} & a_{13} \\ a_{20} & a_{21} & a_{22} & a_{23} \\ a_{30} & a_{31} & a_{32} & a_{33} \end{matrix}] \\ = A_{0} [\begin{matrix} h (m, n) & h (m, n + 1) & h_{n} (m, n) & h_{n} (m, n + 1) \\ h (m + 1, n) & h (m + 1, n + 1) & h_{n} (m + 1, n) & h_{n} (m + 1, n + 1) \\ h_{m} (m, n) & h_{m} (m, n + 1) & h_{mn} (m, n) & h_{mn} (m, n + 1) \\ h_{m} (m + 1, n) & h_{m} (m + 1, n + 1) & h_{mn} (m + 1, n) & h_{mn} (m + 1, n + 1) \end{matrix}] A_{0}^{T} \end{matrix}

(19)

where

\begin{matrix} A_{0} = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ - 3 & 3 & - 2 & - 1 \\ 2 & - 2 & 1 & 1 \end{matrix}] . \end{matrix}

(20)

The partial derivatives $h_{m} (m, n)$ , $h_{n} (m, n)$ , and $h_{mn} (m, n)$ for the square gridding are estimated by the following differencing equations:

\begin{matrix} h_{m} (m, n) = \sum_{i = 1}^{3} \sum_{j = 0}^{3} \frac{a_{ij} i m^{i - 1} n^{j}}{Δ m} = \sum_{i = 1}^{3} \sum_{j = 0}^{3} a_{ij} i m^{i - 1} n^{j} \end{matrix}

(21)

\begin{matrix} h_{n} (m, n) = \sum_{i = 1}^{3} \sum_{j = 0}^{3} \frac{a_{ij} m^{i} j n^{j - 1}}{Δ n} = \sum_{i = 1}^{3} \sum_{j = 0}^{3} a_{ij} m^{i} j n^{j - 1} \end{matrix}

(22)

\begin{matrix} h_{mn} (m, n) = \sum_{i = 1}^{3} \sum_{j = 0}^{3} \frac{a_{ij} i m^{i - 1} j n^{j - 1}}{Δ m Δ n} = \sum_{i = 1}^{3} \sum_{j = 0}^{3} a_{ij} i m^{i - 1} j n^{j - 1} . \end{matrix}

(23)

In the proposed simulation model, this is implemented with the Python SciPy package.

For spike removal, the goal is to smooth out abnormal heights that exist in $H_{M \times N}$ . Spikes will cause some grids to be noticeably higher than the neighboring grids. To smooth the irregular “hills” in the DSM $Ω$ caused by shadows of moving vehicles or pedestrians, we can apply the Gaussian filter (33) $G_{2 D} (i, j; σ)$ with respect to an origin grid $(m, n)$ accordingly to smooth the surface and reduce the noise:

\begin{matrix} G_{2 D} (i, j; σ) = \frac{1}{2 π σ^{2}} e^{- \frac{i^{2} + j^{2}}{2 σ^{2}}} \end{matrix}

(24)

where $Δ m$ and $Δ n$ are the offsets in the numbers of grids from the smoothing center $(m, n)$ . Then, the smoothed height in $\hat{h} (m, n)$ can be calculated as follows:

\begin{matrix} \hat{h} (m, n) = \sum_{i = 0}^{1} \sum_{j = 0}^{1} G_{2 D} (i, j; σ) \cdot \tilde{h} (m + i, v + j) \end{matrix}

(25)

where $σ$ is the standard deviation that determines the width of the Gaussian kernel and $k$ is the kernel radius equal to $4 \times σ$ .

Road Sign and Overhead Power Line Removal

During the generation of DSM, the grids with road signs and overhead power lines will use the top of the signs/wires for the grid heights, which cannot correctly reflect the occlusion. To mitigate those incorrect blind zones, we can filter out the narrow signs and wires from the DSM with morphological operation, opening ( 34 ), which is usually used to first erode an image and then dilate the eroded image to remove small objects and thin lines from an image while preserving the shape and size of larger objects in the image. We can directly apply it to the DSM and update the grid heights like modifying pixel colors of an image.

Dynamic Vehicle Placement

We can further expand the DSM with the time dimension $t$ to update the dynamic vehicle information onto $H (t)$ . The trajectory data we used contain vehicle ID $v$ , class $clas s_{v}$ , timestamp $t$ , and position in 2D $(x_{v (t)}, y_{v (t)})$ . We calculate the 2D orientation $α_{v (t)}$ of vehicle $v$ at time $t$ based on the timestamps before and after:

\begin{matrix} α_{v (t)} = {\begin{matrix} \frac{π}{2}, if x_{v (t + 1)} = x_{v (t - 1)} Λ y > 0 \\ - \frac{π}{2}, if x_{v (t + 1)} = x_{v (t - 1)} Λ y < 0 \\ \arctan (\frac{y_{v (t - 1)} - y_{v (t - 1)}}{x_{v (t + 1)} - x_{v (t - 1)}}), if x_{v (t + 1)} \neq x_{v (t - 1)} \\ α_{lane}, if ∄ v (t - 1) \end{matrix} \end{matrix}

(26)

where $α_{lane}$ is the direction of the current lane.

To simplify the calculations, here we assume that the slope change is not significant so that we can simulate the vehicle pitch altitude by aggregating the vehicle heights over the DSM.

Because of the lack of 3D vehicle shape, we generate random heights ${h_{front}}_{v}$ , ${h_{mid}}_{v}$ , and ${h_{rear}}_{v}$ based on the vehicle class within reasonable range combined with length $l_{v}$ , width $w_{v}$ detected from the LiDAR data, where ${h_{front}}_{v}$ , ${h_{mid}}_{v}$ , and ${h_{rear}}_{v}$ are the three values for different parts of the vehicle. This design enables further detailing with negligible computing resources once a more detailed vehicle model is extracted.

We can create a DSM-like array $H_{V, M \times N} (t)$ with all zeros to store the vehicle heights in each grid $(m, n)$ at time $t$ . For vehicle $v$ at time $t$ , we can calculate the grids it runs on and put the vehicle height onto $H_{V} (t)$ .

First, we can reuse Equation 15 to calculate the vehicle center in the grid coordinate system $(m_{v (t)}, n_{v (t)})$ based on its coordinates $(x_{v (t)}, y_{v (t)})$ :

\begin{matrix} {\begin{matrix} m_{v (t)} = ⌊ \frac{x_{v (t)}}{Δ X} ⌋ \\ n_{v (t)} = ⌊ \frac{y_{v (t)}}{Δ Y} ⌋ \end{matrix} . \end{matrix}

(27)

With width $w_{v}$ , length $l_{v}$ , and vehicle center $(m_{v (t)}, n_{v (t)})$ , we can calculate the four corner coordinates and locate the grids that vehicle $v$ occupies on the road surface before orientation:

\begin{matrix} {\begin{matrix} W_{\frac{1}{2}} = (- \frac{w_{v}}{2}, \frac{w_{v}}{2}, \frac{w_{v}}{2}, - \frac{w_{v}}{2}) \\ L_{\frac{1}{2}} = (- \frac{l_{v}}{2}, - \frac{l_{v}}{2}, \frac{l_{v}}{2}, \frac{l_{v}}{2}) \end{matrix} \end{matrix}

(28)

\begin{matrix} {\begin{matrix} m_{v (t), corners} = m_{v (t)} + W_{\frac{1}{2}} \\ n_{v (t), corners} = n_{v (t)} + L_{\frac{1}{2}} \end{matrix} \end{matrix}

(29)

where $m_{v (t), corners}$ , $n_{v (t), corners}$ are the four sets of corner coordinates calculated from the center ( $m_{v (t)}, n_{v (t)})$ with four sets of half-widths $W_{\frac{1}{2}}$ and lengths $L_{\frac{1}{2}}$ .

Then, for any grid $(m_{v (t), body}, n_{v (t), body})$ within the four vehicle corners of vehicle $v$ at time $t$ , we can apply 2D rotation matrix $R$ with orientation $α_{v (t)}$ to get the oriented grid coordinate $(m_{v (t), body}^{'}, n_{v (t), body}^{'})$ :

\begin{matrix} R = \\ [\begin{matrix} \cos α_{v (t)} & \sin α_{v (t)} & (1 - \cos α_{v (t)}) * m_{v (t)} - \sin α_{v (t)} * n_{v (t)}) \\ - \sin α_{v (t)} & \cos α_{v (t)} & \sin α_{v (t)} * m_{v (t)} + (1 - \cos α_{v (t)}) * n_{v (t)} \end{matrix}] \end{matrix}

(30)

\begin{matrix} [\begin{matrix} m_{v (t), body}^{'} \\ n_{v (t), body}^{'} \end{matrix}] = ⌊ R \cdot [\begin{matrix} m_{v (t), body} \\ n_{v (t), body} \\ 1 \end{matrix}] ⌋ \end{matrix}

(31)

where $R$ is a rotation matrix from OpenCV (35). In our application, we updated $R$ with 1 as the scale and the vehicle center grid as the rotation center.

After getting the vehicle grids, we can place the heights (e.g., ${h_{front}}_{v}$ , ${h_{mid}}_{v}$ , and ${h_{rear}}_{v}$ ) into the corresponding grids in $h_{V} (t)$ to simulate the occlusion with the detailed shape of the vehicle. For every grid $m_{v (t), body}^{'}, n_{v (t), body}^{'}$ :

\begin{matrix} h_{V} (m_{v (t), body}^{'}, n_{v (t), body}^{'} | t) = h_{v} \end{matrix}

(32)

where $h_{v}$ is the height of vehicle $v$ and can be detailed into different parts and $h_{V} (m_{v (t), body}^{'}, n_{v (t), body}^{'}; t)$ is a grid in $H_{V} (t)$ .

After updating the height of every vehicle in time $t$ onto $H_{V} (t)$ , we can add $H_{V} (t)$ onto DSM $H (t_{0})$ and generate a dynamic DSM $H (t)$ for blind zone determination:

\begin{matrix} H (t) = H (t_{0}) + H_{V} (t) \end{matrix}

(33)

where $H (t_{0})$ is the raw DSM with only terrain and infrastructure and $H (t)$ is the DSM fused by both raw DSM and the vehicle heights at time $t$ .

LiDAR SCM

In practice, the roadside LiDAR instrumentation can involve multiple sensors to create coverage overlaps to minimize the impact of occlusions and ensure continuous coverage along a corridor. Similar to the definition in Ge et al. ( 14 ), the “blind zone” of roadside LiDAR sensors refers to (1) the built-in occlusions caused by preconfigured vertical field of view (FOV) designed for mounting on top of vehicles and (2) occluded areas caused by the infrastructure or the moving vehicles blocking the LOSs. To analyze the blind zone in the proposed model, we introduce the LiDAR SCM, whose spatial area is a subset of DSM, to record the coverage details of the LiDAR at different grids within DSM.

We denote the LiDAR SCM of sensor $d$ as $Ψ_{d, M_{d} \times N_{d}} = {H_{d, M_{d} \times N_{d} \times L_{d}}, Λ_{d, M_{d} \times N_{d} \times L_{d}}}$ , where $M_{d}$ and $N_{d}$ are the total number of vertical and horizontal grids in the DSM within the LiDAR sensor coverage. $H_{d, M_{d} \times N_{d} \times L_{d}}$ stores the height location of each laser $l (l = 1, \dots, L_{d})$ , where $L_{d}$ is the total number of laser beams for LiDAR sensor $d$ . $Λ_{d, M_{d} \times N_{d} \times L_{d}}$ stores the occlusion and detection status of segments on a laser beam above all grids within its scanning range. Such data structure in $Λ_{d, M_{d} \times N_{d} \times L_{d}}$ is crucial to improve the computational efficiency without creating a matrix that contains the full 3D space and laser beam index ( $M_{d} \times N_{d} \times Z_{d} \times L_{d}$ ) for blind zone analysis and simulation.

Single-Sensor Single-Beam Grid Height Calculation

The simulation of the coverage of a LiDAR sensor starts with the computation of the height of each beam $l$ of a LiDAR sensor $d$ across each grid $(m, n)$ in the SCM of LiDAR sensor $d$ : $Ψ_{d, M_{d} \times N_{d}} = {H_{d, M_{d} \times N_{d} \times L_{d}}, Λ_{d, M_{d} \times N_{d} \times L_{d}}}$ . We use six parameters to capture the configuration: the 3D coordinates $(x_{d}^{P_{d}}, y_{d}^{P_{d}}, z_{d}^{P_{d}})$ of the LiDAR sensor $d$ in a coordinate system $P_{d}$ and $α_{d}$ , $β_{d}$ , and $γ_{d}$ for its yaw, pitch, and roll angles. It should be noted that the roll of LiDAR sensors should be zero to maximize its detection range, even with tilting. Meanwhile, the coordinate system of $P_{d}$ is the one used by LiDAR to generate the raw point cloud. It is often set by default by manufacturers to use the LiDAR sensor location $(x_{d}^{P_{d}}, y_{d}^{P_{d}}, z_{d}^{P_{d}})$ as the origin of the sensor coverage area. Therefore, we use it as the origin of the DSM. Then, the SCM $h_{d}$ of each sensor $d$ is calculated based on the 3D triangular relationship as derived in the authors’ prior work ( 14 ). If sensor $d$ is placed on the DSM without tilting,

\begin{matrix} h_{d} (m, n, l) = \frac{z_{d}^{P_{c}}}{Δ U} + r \tan ϕ_{l} \end{matrix}

(34)

where $ϕ_{l}$ is laser $l$ ’s altitude angle of LiDAR device model $d$ and can be achieved from either the user manual or configuration file, $\frac{z_{d}^{P_{c}}}{Δ U}$ is from Equation 15, and $r$ is the horizontal grid distance between $(m, n)$ and sensor $d$ :

\begin{matrix} r = \sqrt{m^{2} + n^{2}} . \end{matrix}

(35)

If sensor $d$ is placed with tilting angle $β$ , then the following equations can be used to calculate the height:

\begin{matrix} h_{d} (m, n, l) = {\begin{matrix} \frac{z_{d}^{P_{c}}}{Δ U} + min (h_{+}^{'}, h_{-}^{'}) n < 0 \\ \frac{z_{d}^{P_{c}}}{Δ U} + max (h_{+}^{'}, h_{-}^{'}), n \geq 0 \end{matrix} \end{matrix}

(36)

where $h_{+}^{'}$ and $h_{-}^{'}$ are the two roots of the following quadratic equation, and:

\begin{matrix} h_{+, -}^{'} = r \sin θ_{+, -} \sin β - r \tan ϕ_{l} \cos β . \end{matrix}

(37)

The two roots are reflected in the formula of $\sin θ$ as follows:

\begin{matrix} \sin θ_{+, -} = \\ \frac{- 2 \frac{m^{2}}{n^{2}} \tan ϕ_{l} \sin β \cos β \pm \sqrt{{(2 \frac{m^{2}}{n^{2}} \tan ϕ_{l} \sin β \cos β)}^{2} - 4 (\frac{m^{2}}{n^{2}} \cos^{2} β + 1) (\frac{m^{2}}{n^{2}} \tan^{2} ϕ_{l} \sin^{2} β - 1)}}{2 \frac{m^{2}}{n^{2}} \cos^{2} β + 2} \end{matrix}

(38)

Integration of Multi-Sensor SCMs to DSM

It should be noted that the above two steps are most computationally efficient when the calculation of LiDAR laser beam heights in SCM is executed in the original coordinate system of the raw LiDAR point cloud. Then, the grids of SCMs containing the calculated results are converted into the DSM gridding system without the need to convert the massive raw or simulated LiDAR point cloud data into the DSM coordinate system. Given the SCMs of multiple sensors in their original coordinate systems, $Ψ_{d, M_{d} \times N_{d} \times L_{d}}^{P_{d}} = {H_{d, M_{d} \times N_{d} \times L_{d}}^{P_{d}}, Λ_{d, M_{d} \times N_{d} \times L_{d}}^{P_{d}}}, d = 1, \dots, D$ , where $D$ is the total number of LiDAR sensors considered. For each LiDAR sensor, the $H_{d}^{P_{d}}$ was first calculated within its coordinate system $P_{d}$ , then the resulting $M_{d} \times N_{d}$ grid will be transformed into the DSM coordinate system $P_{0}$ following the $S_{d} \cdot T_{d} \cdot R_{d}$ transformation. Given the coordinate systems are both distance based, there is no scaling needed. Therefore, $S_{d}$ is an identity matrix. The resulting transformation matrix is as follows:

\begin{matrix} T_{d} \cdot R_{d} = [\begin{matrix} \cos α_{d} & - \sin α_{d} & 0 & x_{d}^{P_{0}} \\ \sin α_{d} & \cos α_{d} & 0 & y_{d}^{P_{0}} \\ 0 & 0 & 1 & z_{d}^{P_{0}} \\ 0 & 0 & 0 & 1 \end{matrix}] . \end{matrix}

(39)

Then, the coordinate transformation equation from $P_{d}$ to $P_{0}$ is as follows:

\begin{matrix} [\begin{matrix} \begin{matrix} x^{P_{0}} \\ y^{P_{0}} \\ z^{P_{0}} \end{matrix} \\ 1 \end{matrix}] = T_{d} \cdot R_{d} [\begin{matrix} \begin{matrix} x^{P_{d}} \\ y^{P_{d}} \\ z^{P_{d}} \end{matrix} \\ 1 \end{matrix}] . \end{matrix}

(40)

Given $m = | \frac{x}{Δ U} |$ , $n = | \frac{y}{Δ U} |$ , and $k = | \frac{z}{Δ U} |$ , with a small enough $Δ U$ , the above relationship between coordinates can also be approximately applied to each grid:

\begin{matrix} [\begin{matrix} \begin{matrix} m^{P_{0}} \\ n^{P_{0}} \\ k^{P_{0}} \end{matrix} \\ 1 \end{matrix}] = ⌊ T_{d} \cdot R_{d} [\begin{matrix} \begin{matrix} m^{P_{d}} \\ n^{P_{d}} \\ k^{P_{d}} \end{matrix} \\ 1 \end{matrix}] ⌋ . \end{matrix}

(41)

SCM Laser Beam Segment-Based Blind Zone Status Determination

The proposed blind zone status determination model can be applied to static infrastructure and roadway infrastructure with dynamically moving objects such as vehicles, pedestrians, motorcycles, bicycles, and so forth. It should be noted that the blind zone status can vary in 3D space above a grid $(m, n)$ . A LiDAR laser beam above the grid can be visible above a certain height threshold, $τ$ but can be occluded underneath. Given that most roadside LiDAR sensors have multiple beams from 16 to 32 to 64 to 128, we need to understand the occlusion status of all beams above the grid to determine the blind zone height of the entire grid.

As shown in Figure 3, to analyze each LiDAR $d$ , we first extract the $H_{d, M_{d} \times N_{d} \times L_{d}}$ for all the matched $Ψ_{d, M_{d} \times N_{d}}$ grids. Figure 3a shows the bird’s-eye view of the gridding schemes between DSM and SCM. Because of limited detection range, SCM focuses more on the roadway and pedestrian pathway surfaces and is often smaller than the full DSM covering the area of interests that may contain the full roadway with curvatures.

Figure 3.

Real-world laser beam sampling and DSM comparison. (a) Grid scheme of experiment site. (b) 3D laser status.

Within the SCM, $O (x_{d}^{P_{0}}, y_{d}^{P_{0}}, z_{d}^{P_{0}})$ is the position of LiDAR, and $(m^{d}, n^{d})$ are the indexes of a grid in which the blind zone status concerning LiDAR sensor $d$ is to be modeled. The blue line in Figure 3b indicates the theoretical LOSs by $L_{d}$ laser beams from $O$ if no occlusion. The yellow grids are identified as LOS grids along the theoretical LOS with Bresenham’s line algorithm ( 36 ). As illustrated in Figure 3b, when each laser beam $l_{d}$ of LiDAR sensor $d$ traverses the space above each grid $(m^{d}, n^{d})$ along the Bresenham line, each LOS grid can divide the $L$ lasers into $L$ LOS segments. Each LOS segment in the same LOS grid may have different occlusion or detection status. We classify the status of each LOS segment above an individual LOS grid with two criteria based on its occupancy and detectability:

Occupancy. Any elevation of a grid can be occupied by a moving object (O) (e.g., a vehicle, a bicycle, a motorcycle, a pedestrian, etc.), by an infrastructure object (I), and be vacant for road users to traverse (V).

Detectability. Any grid elevation is detectable by a laser beam when the laser beam can penetrate the entire path from the LiDAR sensor origin to the grid without being occluded (D); otherwise, the grid is not detectable (N).

We denote each segment $ε_{m_{d}, n_{d}, l_{d}}$ of theoretical LiDAR beam path $l_{d}$ over the space of $(m^{d}, n^{d})$ as a LOS segment, as illustrated in Figure 3b. Then, we introduce a two-letter system to represent the occupancy and detectability status sequentially $[V | O | I] [D | N]$ of each LOS segment. The following list provides all feasible combinations, as illustrated in Figure 3b.

OD: An occupied (O) detectable (D) LOS segment.

ON: An occupied (O) but undetectable (N) LOS segment.

VD: A vacant (V) detectable (D) LOS segment.

VN: A vacant (V) but undetectable (N) LOS segment.

ID: A LOS segment occupied by infrastructure objects (I) that is detectable (D).

IN: A LOS segment occupied by infrastructure objects (I) that is undetectable (N).

With the above analysis, the lowest LOS segment among all the laser beams that can reach the current grid determines the blind zone height in the current grid. VN and ON LOS segments are the targeted occlusion and detectability status to define blind zones considered in this paper. It should be noted that V/O is determined by ground truth rather than the LiDAR detection results. VN and ON status can be derived from grid occupancy statuses in ground truth infrastructure and moving object conditions from a simulation-based design and planning study or alternative sensors in field evaluation studies (e.g., drone-based aerial sensors).

Figure 4 shows the detailed process of the proposed quantized laser status determination methods. The quantization means that even though the LOS segments may have different distances in different grids because of the 3D traversing angles of laser beams, we consider only their grid-by-grid detectability and occlusion status and thus irrelevant to the segment lengths.

Figure 4.

Quantized laser status. (a) 3D laser status. (b) Extracted laser status. (c) Laser status by laser ID. (d) Quantized two-dimensional (2D) laser status. (e) Quantized 3D laser status (Vehicle Icons: 39).

Figure 4a shows sample laser beam paths from Figure 3b. The laser path’s coloring shows the different occlusion statuses of the LOS segments. Figure 4, b–d, shows how the blind zone status matrix $Λ_{d, M_{d} \times N_{d} \times L_{d}}$ is constructed to better fit the characteristics of LiDAR data and spatial analysis needs.

We only need to record the laser ID and its occlusion and detection status in each grid, since we can always retrieve the corresponding height information from $H_{d}$ by laser ID and grid for further analysis. Such data structure is consistent with the raw LiDAR data structure and spatiotemporal resolution of SCM and DSM, as Figure 4, c and d , shows. The data structure also allows us to reduce the need to conduct computationally intensive 3D analytics over sparse matrices for occlusion determination. It only takes two rounds to determine the blind zone status without calculating the lasers’ intersections with the complex 3D structures of vehicles, motorcycles, bicycles, and pedestrians. For ease of explanation, we use $(m, n)$ to refer to the matched $(m^{d}, n^{d})$ grid in SCM $Ψ_{d}$ and the $(m^{P_{0}}, n^{P_{0}})$ grid in DSM for single and multiple LiDAR sensors.

The objective of blind zone determination is to segment the 3D space into six categories: ID, VD, OD, IN, VN, and ON. To achieve this objective, we only need to determine the occupancy and detectability status by comparing the SCM height matrix $H_{d}$ with the DSM height matrix $H$ . The static occlusion determination only considers the static DSM height matrix $H (t_{0})$ , whereas the dynamic occlusion determination also considers the dynamic heights $H_{V} (t)$ of road user $V$ added onto the static DSM $H (t)$ , where $t_{0}$ is the time with no traffic, $t$ is the time that there is traffic on the DSM $H$ , as described in dynamic vehicle placement. $H_{V} (t)$ is a DSM-like array that stores the heights of road users at time $t$ based on their occupied grids. When simulating the blind zones caused by moving traffic, only the statuses of those SCM grids altered by vehicles with respect to DSM need to be updated.

The occupancy status in ground truth or reference data for each LOS segment $ε_{m, n, l}$ at grid $(m, n)$ of LiDAR beam $l$ can be determined as follows. Let the two individual digits of $ε_{m, n, l} = OD$ , where $λ (m, n, l | t) = O (m, n, l | t) * 10 + D (m, n, l | t)$ :

\begin{matrix} O (m, n, l | t) = \\ {\begin{matrix} 0, & Vacant (V), if h_{d} (m, n, l) > h (m, n, t) \\ 1, & Occupied by Objects (O), if h (m, n, t_{0}) < h_{d} (m, n, l) \leq h (m, n, t) \\ 2, & Occupied by Iinfrastructure (I), if h_{d} (m, n, l) \leq h (m, n, t_{0}) \end{matrix} . \end{matrix}

(42)

Meanwhile, the detectability status for each LOS segment $ε_{m, n, l}$ can be determined as follows:

\begin{matrix} D (m, n, l | t) = \\ {\begin{matrix} 0, & Not detectable (N), if \exists a LOS grid g : h_{d} (m, n, l | g) \leq h (m, n, t | g) \\ 1, & Detectable (D), if \forall LOS grids : h_{d} (m, n, l | g) > h (m, n, t | g) \end{matrix} \end{matrix}

(43)

where $g$ is the index of the LOS grids between LiDAR $d$ and grid $(m, n)$ , excluding grid $(m, n)$ .

With the occupancy and detectability status, we can infer the category based on Table 1 and fill up $Λ_{d}$ .

Table 1.

Dynamic SCM Zone Type Determination

	Detection
Occupancy	Undetectable (0)	Detectable (1)
I (0)	IN (00)	ID (01)
O (1)	ON (10)	OD (11)
V (2)	VN (20)	VD (21)

Note: SCM = sensor coverage model; I = infrastructure objects; O = occupied; V = vacant; N = undetectable; D = detectable.

Grid Blind Zone Status Determination

The final blind zone status in a grid is determined as follows:

\begin{matrix} Γ (m, n | t, h^{*}) = \underset{l}{V} D (m, n, l | t, h^{*}) \end{matrix}

(44)

where $V$ is the logical OR operator for the status of all laser beams $l$ and $h^{*}$ is the HOI threshold above ground. For every grid $(m, n)$ , if any $l$ within $h^{*}$ above the ground is detectable at time $t$ , then the blind zone detectability status in the grid $(m, n)$ at $t$ is $D$ for $h^{*}$ . If all the $l$ between the ground and $h^{*}$ above the ground are undetectable at time $t$ , then the blind zone status in the grid $(m, n)$ at $t$ is $N$ for $h^{*}$ .

Combined with the trajectory, the proposed model can determine the grid status based on real-world detection demands. For roadway traffic-sensing applications, it is not necessary for the LiDAR devices to see the ground. Instead, as long as the road user is higher than the lowest laser, the LiDAR can sense the user. Furthermore, for most scenarios, the road users will have their lanes to travel, and there can be a lenient requirement for the grids where the traffic will not travel through. In other words, we can get the statistically lowest road user height $h_{\min}$ for each grid based on $H_{V} (t)$ :

\begin{matrix} h_{\min} (m, n) = {\begin{matrix} 0, \\ \min_{t : H_{V} (t) > 0} H_{V} (m, n | t), \end{matrix} \begin{matrix} if {max}_{t} h_{V} (m, n | t) = 0 \\ otherwise \end{matrix} . \end{matrix}

(45)

Then, we can use $h_{\min}$ to determine whether each grid is an object-oriented blind zone. This status holds regardless of whether the blind zone is caused by other traffic (dynamic) or the background (static). Applying Equation 44 with the grid-level $h_{\min}$ , for every grid $(m, n)$ ,

\begin{matrix} Γ (m, n | t, h_{\min} (m, n)) = \underset{l}{V} D (m, n, l | t, h_{\min} (m, n)) . \end{matrix}

(46)

The result of Equation 46 ensures the lowest road user can always be sensed by at least one laser and is fundamental for sensor(s) placement.

Experiment Design

Data Source

In the experiment design, the proposed dynamic blind zone detection model is implemented on a DSM generated from a real-world 3D model with vehicle trajectories retrieved from roadside LiDAR data. The proposed SCM is deployed over grids of 0.3 m × 0.3 m, consistent with the DSM. The LiDAR sensor Velodyne Alpha Prime is installed on the roof of a roadside infrastructure near Boyd Park, New Brunswick, NJ, as depicted in Figure 5, a–c. This location is part of the New Jersey Department of Transportation DataCity Smart Mobility Testing Ground project and is equipped with one 128-beam LiDAR and two 2K cameras.

Figure 5.

Data source. (a) Three-dimensional (3D) model. (b) Surface elevation. (c) Light detection and ranging (LiDAR) point cloud data.

As Figure 5a shows, the 3D point cloud data used in the DSM are from a 3D model generated by videogrammetry in the projected coordinate system of NAD83 (EPSG [European Petroleum Survey Group]: 3615) ( 37 ) from a drone flight executed on June 11, 2023.

The dynamic LiDAR data were collected on August 17, 2023, between 5:50 and 6:00 p.m. The trajectories retrieved by Ouster’s Bluecity Solution ( 38 ) comprise 1,843 unique road users within 5,901 frames (99 frames without detections). The data fields used as the simulation inputs include location information (local_x, local_y) and shape information (length, width, and type). Vehicle types were used to infer the vehicle heights based on typical design vehicle values ( 39 ).

The real-world LiDAR background data used as ground truth for validation were collected from the same LiDAR sensor. The background determination is based on the LiDAR data processing results when no vehicles are detected. The specific time is 1:37 a.m. on July 26, 2023.

As depicted in Figure 5, b and c , the surface elevation shows Route 18 splitting into two levels, and the occlusion caused by the higher-level infrastructure and traffic poses a significant challenge in sensing lower-level traffic. The proposed paper generates the DSM from the 3D model and validates it with LiDAR data. The static blind zone is validated and evaluated based on the 3D laser distribution of field LiDAR data and simulated LiDAR data. Meanwhile, the dynamic blind zone is evaluated based on actual trajectories collected from the same location.

DSM Validation Metrics

Height Difference between DSM and LiDAR

As the base of the simulation, it is essential to evaluate the accuracy of the DSM before performing any evaluation on the simulation model. If the DSM coordinate systems perfectly match that of the LiDAR sensor, all the LiDAR data should stay right on the DSM. For each grid $(m, n)$ , the height difference between DSM and LiDAR in the grid $Δ H_{H - H_{GT}}$ can be calculated as follows:

\begin{matrix} Δ h_{H - H_{GT}} (m, n) = \\ {\begin{matrix} h (m, n) - h_{GT} (m, n), if h (m, n) and h_{GT} (m, n) are not NaN \\ NaN, if H (m, n) or H_{GT} (m, n) is NaN \end{matrix} \end{matrix}

(47)

\begin{matrix} {RMSE}_{Δ h_{H - H_{GT}}} = \\ \sqrt{\frac{\sum_{m = 1}^{M} \sum_{n = 1}^{N} ((Δ h_{H - H_{GT}} (m, n) is not NaN)) {(Δ h_{H - H_{GT}} (m, n))}^{2}}{\sum_{m = 1}^{M} \sum_{n = 1}^{N} (Δ h_{H - H_{GT}} (m, n) is not NaN) (1)}} \end{matrix}

(48)

where $h (m, n)$ is the height of the DSM in the grid $(m, n)$ and $h_{GT} (m, n)$ is the lowest height of the LiDAR point cloud in the grid $(m, n)$ that can be calculated similarly to Equation 16 but with min instead of max. $Δ h_{H - H_{GT}} (m, n)$ shows how much $h (m, n)$ is higher than $h_{GT} (m, n)$ or lower if negative. Because of the sparsity of LiDAR data, $H_{GT}$ are empty in most of the grids, and if grid $(m, n)$ is empty, $h_{GT} (m, n) = NaN$ . NaN is “not a number” and is not a finite value. To calculate the ${RMSE}_{Δ h_{H - H_{GT}}}$ , we can only verify the grids in which $h_{GT} (m, n)$ are not NaN. Therefore, $\sum_{m = 1}^{M} \sum_{n = 1}^{N} Δ h_{H - H_{GT}} {(m, n)}^{2}$ calculates the sum of squares of all the $Δ h_{H - H_{GT}} (m, n)$ that is not NaN, $\sum_{m = 1}^{M} \sum_{n = 1}^{N} (Δ h_{H - H_{GT}} (m, n) is not NaN) (1)$ gets the count of all the $Δ h_{H - H_{GT}} (m, n)$ that is not NaN. With the 2D results of $Δ h_{H - H_{GT}} (m, n)$ presented in a map-like figure, we can calibrate the yaw and pitch angles of the virtual sensor location simulated in the DSM. ${RMSE}_{Δ h_{H - H_{GT}}}$ can be used to calibrate the virtual LiDAR sensor location, especially for the height of the sensor. The smaller the ${RMSE}_{Δ h_{H - H_{GT}}}$ , the closer the virtual LiDAR sensor is to the real LiDAR device and the closer the DSM is from the real-world background.

However, in practice, we found that there were some outliers caused by distant high buildings, which significantly increased the ${RMSE}_{Δ h_{H - H_{GT}}}$ and misled the calibration. We found that using the difference between the 15th and 85th percentiles of $Δ h_{H - H_{GT}}$ can effectively mitigate the outliers and focus on the main part of the road. The silver bullet for the calibration of LiDAR angle $S_{β}$ and height $S_{h}$ can be represented as follows:

\begin{matrix} S_{β} = | Δ H_{{85}} - Δ H_{{15}} | \end{matrix}

(49)

\begin{matrix} S_{h} = \frac{Δ H_{{85}} + Δ H_{{15}}}{2} \end{matrix}

(50)

where $β$ is the tilting (pitch) angle, and $h$ is the height of the LiDAR. The smaller the $S_{β}$ , the closer the $β$ is to the correct value. The closer the $S_{h}$ is to 0, the closer the $h$ is to the correct value. When calculating the 15th and 85th percentiles, grid $(m, n)$ is ignored if $Δ h_{H - H_{GT}} (m, n) is NaN$ .

Grid-Based Percentage Error of DSM

After getting the height offset and the tilting (pitch) offset of the LiDAR sensor, we can further look into the detailed grid matches to approximate the 2D location $(x, y)$ and the rotation (yaw) angle $α$ of the LiDAR. Instead of comparing the overall height difference, we check the percentage of the grids within a difference threshold $δ$ :

\begin{matrix} B_{H} (m, n) = {\begin{matrix} True, \\ False, \\ NaN, \end{matrix} \begin{matrix} if | H_{GT} (m, n) - H (m, n) | \leq δ \\ if | H_{GT} (m, n) - H (m, n) | > δ \\ if H_{GT} (m, n) or H (m, n) is NaN \end{matrix} \end{matrix}

(51)

ϵ_{H} = \frac{\sum_{m = 1}^{M} \sum_{n = 1}^{N} (ε_{H} (m, n) is false) (1)}{\sum_{m = 1}^{M} \sum_{n = 1}^{N} (ε_{H} (m, n) is not NaN) (1)} \times 100 %

(52)

where $B_{H} (m, n)$ shows if the imported DSM matches the LiDAR data in the current grid $(m, n)$ . The $δ$ value we use in this paper is 0.3 m (approximately 1 ft). The percentage error $ϵ_{H}$ shows a mismatch ratio at the current $x$ , $y$ , and $α$ . The smaller the $ϵ_{H}$ , the closer the $(x, y, α)$ are to the correct values. With the percentage error $ϵ_{H}$ , we can further calibrate the virtual sensor location, especially for the horizontal position and the yaw angle.

SCM Validation

With the empty frame of LiDAR point cloud data without any vehicle, we can compare it with the simulation model by sampling the LiDAR point cloud data into the same resolution as the simulation model using Equation 15 and getting the lowest height of the point cloud in each grid.

Grid-Based SCM Verification

The SCM provides the whole 3D spatial distribution $H_{d}$ of the lasers, whereas the LiDAR ground truth data are sparse and can only show up in the grids where they reach and get occluded by something. Therefore, we can only verify the simulation model for true positives and false negatives. If we calculate and mark the overlapping grids of the SCM and the DSM, we can get a similar ground pattern like the curves LiDAR projects onto the ground. However, because of the imperfection of the DSM and real-world interferences, it is barely possible to get point-to-point matches. For every grid $(m, n)$ , the simulated ground pattern grids $K (m, n; τ; t_{0})$ with thickness $τ$ can be determined as:

\begin{matrix} K (m, n; τ; t_{0}) = {\begin{matrix} True \\ False \end{matrix} \begin{matrix} if \min (H_{d} (m, n; t_{0})) < H (m, n; t_{0}) + τ \\ if \forall l : H_{d} (m, n; t_{0}) \geq H (m, n; t_{0}) + τ \end{matrix} \end{matrix}

(53)

where $H_{d} (m, n; t_{0})$ are the heights of all the unoccluded lasers in the grid $(m, n)$ based on the detectability status of $Λ_{d} (m, n, t_{0})$ . If the lowest SCM $H_{d} (m, n; t_{0})$ that is detectable in $Λ_{d}$ is within $τ$ above static DSM $H (m, n; t_{0})$ , then the current grid $(m, n)$ is a ground pattern grid, and $K (m, n; τ; t_{0})$ is true. If for all the laser $l$ , $H_{d} (m, n, l)$ are either undetectable in $Λ_{d}$ or beyond $τ$ above DSM $H (m, n; t_{0})$ , then the current grid $(m, n)$ is considered not a ground pattern grid with a blind zone of at least $τ$ above the DSM.

Then, we can compare the lowest LiDAR point cloud data $h_{GT} (m, n)$ in each grid $(m, n)$ with $K (m, n; τ; t_{0})$ to verify whether the static blind zone determined by the simulation model is correct. For each grid $(m, n)$ :

\begin{matrix} B (m, n; τ; t_{0}) = \\ {\begin{matrix} False, \\ NaN \\ True \end{matrix} \begin{matrix} if K (m, n; τ; t_{0}) = False and H_{GT} (m, n) < H (m, n) + τ \\ if H_{GT} (m, n) is NaN \\ Otherwise \end{matrix} . \end{matrix}

(54)

$B (m, n; τ; t_{0})$ is relatively lenient and only considers the simulation is wrong when the SCM blind zone is confirmed not a blind zone by ground truth LiDAR data. Then, we can calculate the percentage error $ϵ_{H_{d}} (τ)$ of $B (m, n; τ)$ as follows:

ϵ_{H_{d}} (τ) = \frac{\sum_{m = 1}^{M} \sum_{n = 1}^{N} (ε (m, n; τ) is false) (1)}{\sum_{m = 1}^{M} \sum_{n = 1}^{N} (ε (m, n; τ) is not NaN) (1)} \times 100 % .

(55)

The lower the $ϵ_{H_{d}}$ , the closer the SCM is to ground truth.

To be noted, $K$ , $B,$ and $ϵ$ can be used for both static SCM verification and dynamic SCM verification. This paper only uses it for static SCM verification because of the lack of actual vehicle shape data.

Static Blind Zone Evaluation

The background causes a static blind zone, and the static blind zone pattern is relatively permanent once the roadside sensor is installed.

Blind Zone Height

The blind zone height (14) $h_{BZ} (m, n; t)$ determines the lowest objects detected in the grid $(m, n)$ at time $t$ . As Figure 6 shows, the blind zone height is determined by the lowest unoccluded laser. Combined with the background DSM, the blind zone height above ground can better reflect the blind zone impact on road user detection:

\begin{matrix} h_{BZ} (m, n; t) = min (h_{d} (m, n; t)) - h (m, n; t_{0}) \end{matrix}

(56)

where $min (h_{d} (m, n; t))$ is the height of the lowest unoccluded laser in the grid $(m, n)$ at time $t$ and $h (m, n; t_{0})$ is the static DSM height in the grid $(m, n)$ . For static blind zone evaluation, $t = t_{0}$ . The lower the blind zone height $h_{BZ} (m, n; t)$ is, the more the sensor can sense in gird $(m, n)$ . If $h_{BZ} (m, n; t)$ is 0, there is no blind zone in gird $(m, n)$ .

Figure 6.

Blind zone height and effective laser count (Vehicle Icons: 39).

Laser Count

Given HOI $h^{*}$ above DSM $H$ , we can get the laser count $C (m, n; t; h^{*})$ within $h^{*}$ above the ground (14). Based on different detection targets, we can set $h^{*}$ accordingly:

\begin{array}{l} C (m, n; t; h^{*}) = \sum_{l = 1}^{L} (h_{d} (m, n, l; t) \geq h (m, n; t_{0}) \\ and h_{d} (m, n, l; t) < h (m, n; t_{0}) + h^{*}) (1) \end{array}

(57)

where $h_{d} (m, n, l; t)$ is the SCM after static occlusion determination at time $t$ . $h^{*}$ is the HOI based on the target vehicle heights.

Laser count is meaningful for evaluating the detection and reconstruction potential of the sensor data. Having more laser beams within the HOI $h^{*}$ increases the effectiveness of traffic sensing and enhances the details available for 3D reconstruction.

Percentage of Unoccluded Lasers

In addition to the laser count, we can emphasize the impact of occlusion by calculating the percentage of unoccluded lasers $ς_{H_{d}}$ (14). $ς_{H_{d}}$ is defined by the remaining laser count divided by the expected laser count:

ς_{H_{d}} (m, n; t; h^{*}) = \frac{ℂ_{rem} (m, n; t; h^{*})}{ℂ_{\exp} (m, n; h^{*})} \times 100 %

(58)

where $C_{rem} (m, n; t; h^{*})$ is the result from the laser count calculated with the occlusion at time $t$ and $C_{\exp} (m, n; h^{*})$ is the count of the lasers that are $τ$ above DSM $h (m, n)$ without considering the occlusion. The higher the $ς_{H_{d}}$ is, the less the sensing is affected by the occlusion.

Dynamic Blind Zone Evaluation

Unlike the limitation of the static blind zone evaluation, in dynamic blind zone evaluation, we can get ground truth points in the blind zone covered lanes. As Figure 5 shows, the right-most lane on the lower level is in the blind zone. The importation of the vehicles brings the hope of getting some points on the vehicles when they pass the blind zone.

Temporal Mean and Percentile Blind Zone Height

We can use the temporal mean and percentile of blind zone height ( 14 ) over time to evaluate the severity of the blind zone affected by the dynamic traffic with regard to detecting the traffic. Take mean blind zone height, for example:

\begin{matrix} \bar{h_{BZ} (m, n)} = \frac{\sum_{t = t_{0}}^{T} h_{BZ} (m, n; t)}{T} \end{matrix}

(59)

where $T$ is the total time and $h_{BZ} (m, n; t)$ is the blind zone height in the grid $(m, n)$ at time $t$ .

Temporal Mean Laser Count

We can use the temporal mean of laser count $\bar{C (m, n; h^{*})}$ over time (14) to evaluate the severity of the blind zone affected by the dynamic traffic with regard to reconstructing and tracking different types of road users:

\begin{matrix} \bar{C (m, n; h^{*})} = \frac{\sum_{t = t_{0}}^{T} C (m, n; t; h^{*})}{T} \end{matrix}

(60)

where $C (m, n; t; h^{*})$ is the laser count in the grid $(m, n)$ at time $t$ .

Aggregated Time in Blind Zone and Maximum Consecutive Time in Blind Zone

Aggregated time in the blind zone $T_{agg}$ and the maximum consecutive time in the blind zone $T_{\max}$ from the authors’ prior paper (14) can be used to evaluate the temporal dynamic blind zone impact over road user tracking with the new DSM-based blind zone determination. $T_{agg}$ is quantized as a percentage of the total time, whereas $T_{\max}$ is calculated as an absolute value. A grid $(m, n)$ is considered in the blind zone at HOI $h^{*}$ at time $t$ when $h_{BZ} (m, n; t)$ is more than $h^{*}$ .

Analysis of Results

DSM

Coordinate System Calibration

As Figure 7 shows, with three pairs of reference points, the 3D photogrammetry model is matched with the real-world coordinate system. Figure 7a shows the selected three distant points. Figure 7b shows the flattened model placed over Google Maps. By doing this, the calibrated 3D model can be verified to have the correct north and up directions. With the three pairs of matched coordinates from Google Maps, we can calculate the geo-information for the 3D model. Figure 7c shows the GeoTIFF file generated from the 3D model that can be directly imported into ArcGIS/ArcMap. GeoTIFF is based on the TIFF format and is used as an interchange format for georeferenced raster imagery. From Figure 7c, we can see that NJ-18 and its ramps are properly aligned and the lanes are well matched at the boundaries.

Figure 7.

Coordinate system calibration result. (a) Manually picked reference points. (b) Flattened three-dimensional (3D) model matched with Google Maps. (c) GeoTIFF with coordinates visualized in ArcMap.

DSM Generation

Figure 8, a–d, shows the 3D views of the DSM filtered by different Gaussian kernels. The 3D model was generated from high-angle drone video and had a lot of noise. The noise affects the DSM quality, as Figure 8a shows. By applying the Gaussian filter, we can see the results with different $σ$ values, as Figure 8, b–d, shows. The purpose of applying the Gaussian filter is to mitigate the glitches in the DSM while keeping the correct details as many as possible. From the 3D views of the DSM, the $σ$ of 2 outperforms other $σ$ values by smoothing the road while maintaining the overall geometry.

Figure 8.

Gaussian smoothed digital surface model (DSM) with different $σ$ . (a) Raw DSM. (b) DSM (Gaussian $σ$ = 1). (c) DSM (Gaussian $σ$ = 2). (d) DSM (Gaussian $σ$ = 3).

Figure 9 shows the road sign and overhead power line removal. Figure 9b shows the original DSM that contains the road signs. It is correct, but it will not cause blind zones that affect the real-world traffic sensing. Figure 9, b and c , shows the DSM processed by morphological operations that no longer reflects road signs or overhead power lines while deliberately retaining the thick vertical pillars. By removing those higher objects from the DSM, a more realistic blind zone can be simulated.

Figure 9.

Road sign and overhead power line removal with image opening. (a) Street view. (b) Digital surface model (DSM). (c) DSM after opening. (d) Combined DSM.

DSM Validation

Figure 10a shows the top view of the 3D model for DSM, focusing on the sight of interest (SOI) of the simulation. The white blocks mean that there are no available heights in the grids. Figure 10b shows the top view of the LiDAR data over the SOI. Figure 10d shows the height difference between DSM and LiDAR $Δ H_{H - H_{GT}}$ . The color bar describes how dynamic the heights in Figure 10, a–c, are. Figure 10c shows some irregularly high height differences at the bottom right. Comparing it with Figure 10, a and b , we found a large area of dark red (60 ft+) in the DSM, but the LiDAR data are mixed with green (20 ft) and red (60 ft). ${RMSE}_{Δ H_{H - H_{GT}}}$ for the above calibration is 8.704 ft, but during the calibration, we found that the ${RMSE}_{Δ H_{H - H_{GT}}}$ was not perfect as the evaluation metric because of the outliers. To better reflect the majority of the matched points, we used the 15th and 85th percentiles of $Δ H_{H - H_{GT}}$ as the evaluation matrix. The 15th percentile is −0.862 ft, and the 85th percentile of HDDL is 0.919 ft. It should be mentioned, the max $Δ H_{H - H_{GT}} = 81.566$ ft and the min $Δ H_{H - H_{GT}} = - 40.795$ ft under this configuration, which explains why the ${RMSE}_{Δ H_{H - H_{GT}}}$ is not reliable for the calibration. The $ϵ_{H}$ for $δ = 1$ is 28.9%, which matches with the 15th–85th percentile, meaning that 71.1% of the grids have a height difference smaller than 1 ft. Figure 10d provides a more straightforward view, showing whether the grid is within $δ$ or beyond $δ$ . From Figure 10d, we can see that the majority of the true grids are distributed on the road, which is quite satisfying, whereas the false grids are on the infrastructure, trees, or the edges, which may be a side effect of the Gaussian filter. Figure 10e shows the result of the calibrated LiDAR installation in the DSM coordinate system. Figure 10f shows the histogram of $Δ H_{H - H_{GT}}$ , where most of the $Δ H_{H - H_{GT}}$ are close to 0.

Figure 10.

Digital surface model (DSM) calibration. (a) DSM over the sight of interest (SOI). (b) Light detection and ranging (LiDAR) data over the SOI. (c) $Δ H_{H - H_{GT}}$ . (d) B_H(m, n). (e) Calibrated LiDAR data over DSM. (f) $Δ H_{H - H_{GT}}$ histogram.

SCM

SCM Validation

Figure 11a shows the ground pattern of the simulation model with $τ = 0.5$ , which means that the simulation model can sense any object in the yellow grids that is higher than $τ = 0.5$ ft. Figure 11b verifies whether there are any LiDAR data in the simulated blind area with $τ$ , where the original LiDAR coverage can be found in Figure 10b. The yellow grids are the ones that can find matched LiDAR data. Purple means that the simulation model cannot sense the grids, whereas in the LiDAR data, there are points in the grids within the height $τ$ . The $ϵ_{H_{d}}$ for $τ = 0.5$ ft is 24.3%, meaning 75.7% of the LiDAR points are reflected in the simulation model. If we look into the details of Figure 11b, we can see that most of the roads are validated, and the major differences are located in the bottom-right corner, the bottom-left corner, and the top. The bottom right is as expected because, in Figure 10, a and b , the DSM has shown some differences with the LiDAR data. The top left is not reflected in the simulation model because it is supposed to be occluded. The 3D model was collected after the actual LiDAR was installed, and therefore the virtual LiDAR is occluded by the real LiDAR in the DSM.

Figure 11.

Sensor coverage model validation. (a) Ground pattern with $τ$ = 0.5, before morphological operations. (b) B(m, n; $τ$ ; t₀) with $τ$ = 0.5, before morphological operations. (c) Ground pattern with $τ$ = 0.5, after morphological operations. (d) B(m, n; $τ$ ; t₀) with $τ$ = 0.5, after morphological operations.

Similarly, because we have slightly changed the position of the virtual LiDAR to avoid it being completely occluded by the real LiDAR and other roadside instruments, the blind zones caused by nearby poles and trees are at different angles in the bottom left. However, it does show a problem, as has been red-marked, that a large sign installed on a thin pole will raise the DSM to the top of the sign, which is not correct in the real world. Similar problems will also happen to the gantries, which are not shown in the figures.

Figure 11, c and d , shows the simulated ground pattern and the validation after removing the road signs with morphological operations. The red-marked blind zones in Figure 11b are correctly mitigated. However, compared with Figure 11b, the actual LiDAR ground pattern, some thin poles were also removed from the DSM, and some actual blind zones were incorrectly mitigated. Therefore, the morphological operations should be used cautiously when necessary. For example, removing the overhead wires in the city street is necessary. Otherwise, the wires in the DSM will occlude the lasers like high walls.

Static Blind Zone Evaluation

With the advantage of simulations, static blind zones can be calculated by placing no traffic in the background. With regard to their contribution to traffic safety and management, both static and dynamic blind zones can affect sensor coverage and hinder the roadside system’s ability to predict near-miss events. However, static blind zones can be easily calculated and subsequently avoided or mitigated during the design phase. In contrast, dynamic blind zones require actual or simulated trajectories that accurately reflect traffic patterns, such as traffic flow, large vehicle ratios, and vehicle distances. To analyze dynamic blind zones during the design phase, accurate microscopic traffic simulations are essential.

Figure 12 shows the laser count $ℂ$ and percentage of unoccluded lasers $ς_{H_{d}}$ at three different HOI for sedans/coupes (4 ft), SUVs/pickups (6 ft), and trucks (14 ft). From Figure 12, a, c, and e , we can see that with the increase of $h^{*}$ , we get more laser beams within the HOI. For most of the lanes, we have at least a small segment with more than 10 laser beams in the HOI, which helps capture the whole skeleton of the passing traffic. If we downgrade the expectation to detection only, the majority of the road is covered by at least one laser beam, except for the bottom-left “stubborn” area on the lower level that stays completely blank (0 laser beam within HOI) for $h^{*} \leq 6$ . Even if we increase the HOI, there is only limited improvement at $h^{*} = 14$ . The traffic in the blank area will not be detected because of the lack of laser beams. We can further investigate it through $ς_{H_{d}}$ . Vertically comparing Figure 12, b, d, and f , we see that $ς_{H_{d}}$ at different $h^{*}$ share a similar pattern, which means that almost all the missing traffic in the blank area in $C$ , no matter the heights, will have occlusion issues. Horizontally comparing Figure 12, a, c, and e , with Figure 12, b, d, and f , the blank area in $C$ matches with the low percentage area in $ς_{H_{d}}$ , which means that the LiDAR could have sensed the traffic in the blank area but cannot because of the occlusion.

Figure 12.

Laser count and percentage of unoccluded lasers within $τ$ (ft). (a) $C$ (m, n; t; h*), h* = 4. (b) $ς_{H_{d}}$ (m, n; t; h*), h* = 4. (c) $C$ (m, n; t; h*), h* = 6. (d) $ς_{H_{d}}$ (m, n; t; h*), h* = 6. (e) $C$ (m, n; t; h*), h* = 14. (f) $ς_{H_{d}}$ (m, n; t; h*), h* = 14.

Dynamic Blind Zone Simulation

Figure 13 shows dynamic LiDAR simulation results for four frames at 2-s intervals and a special case frame. The left side of Figure 13 illustrates the ground pattern of the LiDAR laser beams affected by the traffic, whereas the right side shows the corresponding blind zone height, ranging from 0 to 14+ ft above the ground.

Figure 13.

Dynamic light detection and ranging (LiDAR) simulation.

For the upper-level traffic, because of both the shorter distance and the higher vehicle positions, the dynamic blind zone caused by large vehicles on the upper level may affect both the upper-level sensing and the lower-level sensing, as the special case at time $t$ shows in Figure 13. However, smaller vehicles will not occlude other vehicles with the current LiDAR position and oriented angle. Concerning lower-level traffic, large vehicles tend to cause a significant blind zone, but fortunately, they are on the lower level, far away from the LiDAR. Thanks to the additional elevation, the LiDAR has a better angle for sensing lower-level traffic, and the smaller vehicles don’t cause larger blind zones, even though the distances are much longer.

Dynamic Blind Zone Evaluation

We can further look into actual traffic to analyze the dynamic blind zone. Figure 14 shows the blind zone heights aggregated by frames over $t = 10$ min. The color scheme is set to show from 0 to 15+ ft because almost all the road users become invisible when the blind zone is higher than 15 ft. The blank area in Figure 12e results in a blind zone exceeding 14 ft in Figure 14 because the lowest laser beam determines the blind zone height. The static blind zone remains the same height, and the dynamic blind zone updates the height over time. Road users become invisible when entering a blind zone higher than themselves. If the blind zone is lower than 4 ft, most vehicles will likely be detected from the LiDAR data. When the blind zone is higher than 4 ft but lower than 6 ft, sedans and coupes may remain invisible in the LiDAR data. When the blind zone exceeds 14 ft, only trucks have the chance of being revealed; sedans, coupes, SUVs, and pickups may become undetectable in the LiDAR data. Figure 14a shows the mean $H_{BZ}$ over $t$ , from which we can find that, in addition to the static blind zone, there is a new sector caused by the traffic instead of the static background. After investigation, we found it was caused by a stopped truck waiting for the green light on the upper level, as the special case in Figure 13 shows. Figure 14, b and c , shows the 15th and 50th percentiles of the $H_{BZ}$ . Comparing Figure 14, b and c , we see that the medium traffic doesn’t create too many new blind zones over the light traffic blind zone, but there is a trend that the left part of the upper level and the right part of the lower level will have higher dynamic blind zone. If we look at Figure 14d, the 85th percentile of $H_{BZ}$ , the blind zone issue becomes more challenging, and approximately half of the road has a $H_{BZ}$ higher than 6 ft. This is just the 85th percentile of peak-hour traffic.

Figure 14.

Mean and 15th, 85th, and 100th percentiles of temporal aggregated blind zone height (above ground). (a) $\bar{H_{BZ}}$ . (b) 15th percentile $H_{BZ}$ . (c) 50th percentile $H_{BZ}$ . (d) 85th percentile $H_{BZ}$ .

Although the blind zone height in Figure 14d is not optimistic, we can still count on the tracking algorithm to predict if there is any road user on the occluded grid. Figure 15 shows the grid-based blind zone duration and the maximum consecutive time in the blind zone at HOI levels of 4, 6, and 14 ft, representing the heights of sedans/coupes, SUVs/pickups, and trucks. Ignoring the elephant in the room, since the static blind zone will last “permanently,” we can explore how the traffic affects the dynamic blind zone. Figure 15, a, c, and e , shows the percentage of total time in the blind zone. As Figure 15a shows, even for $h^{*} = 4$ ft, the aggregated time in the blind zone $T_{agg}$ remains below 20% for most of the new dynamic blind zone road, with a few areas staying in the blind zone for 40%–60% of the total time. Figure 15, c and e , shows that $T_{agg} (h^{*} = 6 ft)$ and $T_{agg} (h^{*} = 14 ft)$ mitigate the time in the dynamic blind zone. Figure 15, b, d, and f , shows the maximum consecutive time in the blind zone $T_{\max}$ . The capped $T_{\max}$ values are mainly on the upper level because of the red light. The vehicles waiting for the green light caused temporary static blind zones. The higher vehicles such as trucks and buses cause a larger range of blind zones that affect the lower levels too. Other than the static blind zone, most of the lower level has a $T_{\max}$ shorter than 5 s, and the traffic in the blind zone is not too hard to predict.

Figure 15.

Time in the blind zone and maximum consecutive time in the blind zone. (a) Time in the blind zone—4 ft above ground. (b) Maximum consecutive time in the blind zone—4 ft. (c) Time in the blind zone—6 ft above ground. (d) Maximum consecutive time in the blind zone—6 ft. (e) Time in the blind zone—14 ft above ground. (f) Maximum consecutive time in the blind zone—14 ft.

The major issues with this LiDAR setup are the infrastructure-caused static blind zone on the lower level and the large vehicle-caused blind zone. Limited by regulations, power availability, and geometry, the positions and heights of the sensors are restricted. However, we can explore potential solutions, such as trailer sites with solar power.

Multi-Sensor Blind Zone Evaluation

We can place more sensors and evaluate the coverage to enhance the sensing capability and gather more information from the site.

By adding one LiDAR looking from the lower level, the blind zone caused by the infrastructure can be significantly reduced. Figure 16 shows the LiDAR ground pattern generated using Equation 53. Figure 16, a and b , shows the ground pattern using only one LiDAR. In Figure 16a, the infrastructure causes significant static blind zones because of elevation differences. In Figure 16b, to make the LiDAR able to cover the upper level, it has to be installed much higher, making the traffic below the LiDAR invisible. If we combine both LiDAR devices, as Figure 16c shows, the coverage is significantly improved, and dense LiDAR lasers cover the whole road surface.

Figure 16.

Multi light detection and ranging (LiDAR) ground pattern with one trailer-based location. (a) LiDAR ground pattern—infrastructure location. (b) LiDAR ground pattern—trailer location. (c) Fused LiDAR ground pattern.

Figure 17 shows the lowest blind zone height map for two sensors. The improvement achieved by using two sensors is significant compared with using any single LiDAR sensor in this complex scenario. The static blind zone caused by terrain and infrastructure is well covered. The mesh-like pattern in Figure 17c results from resampling during mapping the two results onto the same DSM and should be disregarded; the actual blind zone does not exhibit this pattern. However, this introduces new challenges in sensor placement optimization, as multi-sensor placement problems can become exponentially more complex.

Figure 17.

Blind zone height map of multi-sensor setup. (a) Light detection and ranging (LiDAR) blind zone height map—infrastructure location. (b) LiDAR blind zone height map—trailer location. (c) Fused LiDAR blind zone height pap.

Dynamic Blind Zone Height Map of Different LiDAR Devices

As Figure 18 shows, model VLP-16 has 16 laser beams distributed between −15° and 15°, HDL-32E has 32 beams distributed between −30° and 10°, and Alpha Prime has 128 beams distributed between −25° and 15°. With all of them installed at 23 ft with 10° tilting toward the roadway, it is obvious that HDL-32E benefits from its −30° vertical FOV range and has a smaller blind spot beneath the LiDAR than the other two models. However, with significant more lasers, Alpha Prime has better robustness against occlusion-caused blind zones thanks to the smoother blind zone boundaries. For roadside LiDAR deployment, it is important for the manufacturers to not only increase the laser beams but also expand the vertical FOV further downwards.

Figure 18.

Dynamic blind zone height map of different light detection and ranging (LiDAR) devices (VLP-16, HDL-32E, and Alpha Prime).

Statistical Blind Zone Analysis of Different LiDAR Devices

Figure 19 shows some interesting patterns for comparison of different LiDAR devices. Unlike Figure 18, Alpha Prime leads the performance because of its 128 lasers, especially for the laser count within 4 ft above ground. It also relieves some areas from the blind zone duration. From the average blind zone height and laser count within 4 ft above ground, neither the 16-beam LiDAR VLP-16 nor the 32-beam LiDAR HDL-32E can satisfy the demand for road traffic reconstruction or even just tracking under the multi-level roadway condition.

Figure 19.

Average blind zone height, laser count, time in the blind zone, and maximum consecutive time in the blind zone of different light detection and ranging (LiDAR) devices.

Conclusions

This paper has presented an improved LiDAR blind zone evaluation model that simulates laser beams over a background generated using videogrammetry to evaluate and determine the optimal LiDAR configuration including six parameters: installation $x$ , $y$ , $z$ , rotation angle $α$ , tilting angle $γ$ , and laser beam angles $ϕ$ s. We have validated and evaluated the model based on actual LiDAR data and road-user trajectories. The model accepts road users’ shapes and 3D orientation for a more detailed blind zone analysis and can calculate and evaluate the mitigated blind zone after fusing data from different LiDAR sensors.

With regard to applications, the inclusion of real-world background data allows for a meaningful evaluation of the blind zone with different LiDAR positions and angles. Trajectories help determine whether the blind zone affects traffic sensing. All simulation and evaluation results can be directly projected onto real-world coordinate systems such as NAD83 and WGS84, validated by generating 2D GeoTIFF images and importing them into ArcMaps. In addition to virtual LiDAR simulation for design, the real-world coverage map with the placement of actual LiDAR devices also allows geo-referenced LiDAR sensor search, which can speed up data query for specific incidents or accidents.

With regard to performance, to enhance efficiency in dynamic LiDAR blind zone simulation, we sampled the 3D background model into a DSM-like structure. We also proposed a new data structure that stores the occlusion status of each laser beam in each grid rather than in 3D space with accurate heights. This approach can save significant memory compared with traditional voxel-based structures, and even in comparison with the grid-laser-based structure from our prior paper, it saves seven times the memory for the structure when integrating LiDAR data with 3D models, such as in Digital Twin applications. At a resolution if 0.3 m, the sensor coverage analysis operates at approximately 53.35–58.97 Frames Per Second (FPS) for each 16-beam LiDAR device, 30.05–32.25 FPS for each 32-beam LiDAR device, and 8.25–9.36 FPS for each 128-beam LiDAR device, using a desktop with an Intel 13900KF CPU, 16 GB of 5,200 MHz DDR5 RAM, and an Nvidia RTX 3090 GPU. The processing involves updating the DSM to reflect dynamic traffic over a static background, determining the blind zone height, and generating the dynamic blind zone height results. Combining this with the 10 FPS frame rate of the Velodyne Alpha Prime and real-time edge processing at 10 FPS, the maximum system latency is expected to be less than 321 ms. This latency is calculated from the moment a road user appears to when the blind zone is modeled and analyzed, assuming that the processing units are on the same local area network that internal transmission time can be ignored. With 321 ms latency, a vehicle traveling at 40 mph can move 5.74 m between two frames. This can increase to 7.175 m at 50 mph, 8.61 m at 60 mph, and 10.045 m at 70 mph. In our experiment, the processed trajectory has an additional 0.11- to 0.15-s delay because of Internet transmission. The primary source of latency is the frame rate of the LiDAR devices. Unlike cameras, which can easily achieve 60 FPS or even 120–240 FPS, LiDAR devices typically operate at 5–20 rotations per second. Restricted by the implementation, the processing time with Python may not be competitive with more efficient code implemented in C++. The methodology proposed in this paper that uses a DSM-like structure for 3D analysis is theoretically computing more efficiently than voxel-based methods and has better robustness against an increasing number of road users compared with mesh-based methods.

With regard to validation, we propose two matrices to validate the generated DSM and one matrix to validate the simulated static blind zone with actual LiDAR data. With regard to evaluation, we updated two matrices from the prior paper and proposed one new matrix to evaluate the static blind zone affected by the background (terrain and infrastructures). We reapplied three matrices with new “in blind zone” determination criteria to evaluate the dynamic blind zone affected by traffic.

The DSM-based LiDAR 3D blind zone analysis model sacrifices actual 3D information in favor of computational and memory efficiency. As a trade-off, the model struggles to handle multilayered 3D environments, such as tunnels or bridges with under-bridge traffic. Although the road users placed on the DSM reflect 3D orientation based on each vehicle’s direction and the terrain, the model does not accurately account for the precise position of the four wheels. Consequently, it fails to capture vertical movements caused by small obstacles with accuracy. Another limitation is that the blind zone determination operates at the DSM resolution. For instance, if the DSM is sampled on a 0.3 m × 0.3 m grid, the vehicle will also be sampled at a resolution of 0.3 m to avoid unnecessary computational or storage costs caused by mismatched accuracies. For a site-specific blind zone analysis, a DSM-based background and actual trajectory data are generated and imported to assess the detection performance of using LiDAR devices on Route 18, NJ. The proposed model calculates the 3D blind zones based on the background (if provided) for static blind zones and traffic trajectories for dynamic blind zones. The output visualizes the 3D blind zones on 2D maps, color coded by different evaluation matrices, and projects the maps onto a real-world coordinate system (if applicable). In addition to the LiDAR point cloud data simulation, it not only shows the seen objects; the results also provide insights into what is really covered and what is not, even if no objects are seen. The optical geometry-based 3D SCM can be easily sampled at any resolution to match any background model. Thanks to the ease of creating a DSM background, the model can be transferred to evaluate the sensor blind zones on a virtual background generated for different scenarios.

The major blind zone findings include the following:

LiDAR Installation

○ Height. The proper installation height for LiDAR sensors can help reduce occlusions and allow laser beams to penetrate gaps between vehicles. However, increasing the installation height comes at the cost of a larger blind spot directly below the sensor.

○ Angle. Tilting the LiDAR sensor effectively reduces the blind spot directly underneath the sensor and increases beam density near the LiDAR. However, excessive tilting can create new blind zones on the far side.

○ Position. Ideally, the sensor should be positioned as close to the road as possible. If large vehicles travel in the nearest lane, it may be beneficial to position the sensor slightly farther from the road.

LiDAR Design

○ Vertical FOV. The vertical FOV is important for roadside deployment. With a larger FOV downwards, the blind spot beneath the LiDAR will become smaller. Meanwhile, the LiDAR can be placed closer to the road.

○ Laser count. The lasers become sparse on distant objects. For roadside deployment, existing 16- or 32-beam LiDAR devices have no robustness against complex terrain, and the coverage is significantly impaired. But if the manufacturers can focus all the lasers from a 64-beam LiDAR downwards, it may maintain a good balance between cost and effective coverage.

External Factors

○ Large vehicle. Higher and longer vehicles tend to create larger and longer-lasting blind zones.

○ Congestion. In light traffic conditions, blind zones are less frequent and severe. Congestion worsens dynamic blind zones because of increased traffic density and prolonged periods that grids that remain in blind zones. Sensor configurations should consider more severe blind zone situations during peak traffic hours rather than during free flow.

○ Terrain. Complex road terrain and infrastructure can cause significant static blind zones, making a multi-sensor setup necessary.

○ Surrounding. Small objects close to the sensor can create larger blind zones than large vehicles farther away. It is important to ensure that nearby devices do not obstruct the laser path between the LiDAR and the road/traffic.

Overall, it is generally advisable to increase the height to twice that of the tallest vehicle that may appear in the ROI. The resulting increased blind spot can be mitigated by adjusting the tilting angle, until the sparse upper lasers hit the road. When feasible, the sensor pole should be placed as close to the road as possible.

The exploration has been both interesting and challenging. Future work includes further refining the background model with trajectory data, evaluating the impact of additional weather and environmental factors, assessing the effect of system latency on safety-related applications, designing an efficient optimization method to fine-tune the six parameters for various objectives—such as maximizing detection coverage, maximizing 3D reconstruction coverage, and minimum costs with acceptable coverage trade-offs—and exploring the migration to on-board sensor coverage and camera coverage analysis. From a phased perspective, background model refinement and system latency assessment are short-term goals that depend on interdisciplinary expertise. In contrast, placement optimization for different tasks can be viewed as a long-term goal, requiring more efficient methods to analyze the overall impact of blind zones on traffic safety and to determine optimal configurations.

Footnotes

Acknowledgements

The corresponding author, Yi Ge, is deeply grateful to his wife, Yuting Wang, for her patience and support, especially during the revision process, which coincided with their honeymoon.

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: Peter J. Jin, Yi Ge; data collection: Yi Ge, Anjiang Chen; analysis and interpretation of results: Yi Ge, Peter J. Jin; draft manuscript preparation: Yi Ge, Peter J. Jin. All authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper is based on work partially supported by Federal Highway Administration Research Project 21-60168, Middlesex County Resolution 21-821-R and the National Science Foundation under Grant No. 2133516. This paper is also partially based on work supported by the U.S. Department of Homeland Security under Grant Award 22STESE00001-03-02.

ORCID iDs

Yi Ge

Peter J. Jin

Anjiang Chen

The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Department of Homeland Security.

References

Zhang

Xiao

Coifman

Mills

J. P.

Vehicle Tracking and Speed Estimation from Roadside Lidar. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, Vol. 13, 2020, pp. 5597–5608.

Zhao

Liu

Zheng

Detection and Tracking of Pedestrians and Vehicles Using Roadside LiDAR Sensors. Transportation Research Part C: Emerging Technologies, Vol. 100, 2019, pp. 68–87.

Cui

Sun

Zhao

Automatic Vehicle Tracking with Roadside LiDAR Data for the Connected-Vehicles System. IEEE Intelligent Systems, Vol. 34, No. 3, 2019, pp. 44–51.

Zhang

Jin

P. J.

Weighted Bayesian Gaussian Mixture Model for Roadside Lidar Object Detection. arXiv Preprint arXiv:2204.09804, 2022.

Zheng

Zhang

Tian

Automatic Vehicle Classification Using Roadside LiDAR Data. Transportation Research Record: Journal of the Transportation Research Board, 2019. 2673: 153–164.

Zheng

Tian

A Novel Method of Vehicle–Pedestrian Near-Crash Identification with Roadside LiDAR Data. Accident Analysis & Prevention, Vol. 121, 2018, pp. 238–249.

Sun

Zheng

Yue

Automatic Background Filtering Method for Roadside LiDAR Data. Transportation Research Record: Journal of the Transportation Research Board, 2018. 2672: 106–114.

Zhang

Zheng

Wang

Fan

Chen

Automatic Background Construction and Object Detection Based on Roadside LiDAR. IEEE Transactions on Intelligent Transportation Systems, Vol. 21, No. 10, 2019, pp. 4086–4097.

Federal Communications Commission. Dedicated Short Range Communications (DSRC) Service. https://www.fcc.gov/wireless/bureau-divisions/mobility-division/dedicated-short-range-communications-dsrc-service. Accessed January 22, 2024.

10.

3D Dynamic Blind Zone Modeling and Analysis for High-Resolution Roadside LiDAR Sensors. Doctoral dissertation. Rutgers University-School of Graduate Studies, New Brunswick, NJ, 2024.

11.

Zhang

T. T.

Chen

Sartipi

Jin

P. J.

Hash-Based Gaussian Mixture Model (HGMM) for Roadside LiDAR Smart Infrastructure Applications. IEEE Transactions on Intelligent Transportation Systems, Vol. 25, No. 10, 2024, pp. 12968–12970.

12.

Kloeker

Liu

Wei

Eckstein

Framework for Quality Evaluation of Smart Roadside Infrastructure Sensors for Automated Driving Applications. Proc., 2023 IEEE Intelligent Vehicles Symposium (IV), IEEE, New York, June 4, 2023, pp. 1–7.

13.

Shi

Zhao

A Lifelong Framework for Data Quality Monitoring of Roadside Sensors in Cooperative Vehicle-Infrastructure Systems. Computers and Electrical Engineering, Vol. 100, 2022, p. 108030.

14.

Jin

P. J.

Zhang

T. T.

Chen

Roadside LiDAR Sensor Configuration Assessment and Optimization Methods for Vehicle Detection and Tracking in Connected and Automated Vehicle Applications. Transportation Research Record: Journal of the Transportation Research Board, 2024. 2678: 184–203.

15.

Liu

Chitlangia

Agnihotri

Zhao

Investigating the Impact of Multi-Lidar Placement on Object Detection for Autonomous Driving. Proc., IEEE/CVF Conference on Computer Vision and Pattern Recognition, IEEE Xplore, 2022, pp. 2550–2559.

16.

Liu

Perception Entropy: A Metric for Multiple Sensors Configuration Evaluation and Design. arXiv Preprint arXiv:2104.06615, 2021.

17.

Dybedal

Hovland

Optimal Placement of 3D Sensors Considering Range and Field of View. Proc., 2017 IEEE International Conference on Advanced Intelligent Mechatronics (AIM), IEEE, New York, July 3, 2017, pp. 1588–1593.

18.

Liu

Arief

Zhao

Where Should We Place Lidars on the Autonomous Vehicle? - An Optimal Design Approach. Proc., 2019 International Conference on Robotics and Automation (ICRA), IEEE, New York, 2019, pp. 2793–2799.

19.

Dybedal

Hovland

GPU-Based Optimisation of 3D Sensor Placement Considering Redundancy, Range and Field of View. Proc., 2020 15th IEEE Conference on Industrial Electronics and Applications (ICIEA), IEEE, New York, 2020, pp. 1484–1489.

20.

Kim

T. H.

Park

T. H.

Placement Optimization of Multiple Lidar Sensors for Autonomous Vehicles. IEEE Transactions on Intelligent Transportation Systems, Vol. 21, No. 5, 2019, pp. 2139–2145.

21.

Kim

T. H.

G. H.

Yun

H. S.

Yun

K. S.

Park

T. H.

Placement Method of Multiple Lidars for Roadside Infrastructure in Urban Environments. Sensors, Vol. 23, No. 21, 2023, p. 8808.

22.

Cai

Jiang

Zhao

Liu

Analyzing Infrastructure Lidar Placement with Realistic Lidar Simulation Library. Proc., 2023 IEEE International Conference on Robotics and Automation (ICRA), IEEE, New York, 2023, pp. 5581–5587.

23.

Xiang

Xia

Han

OPV2V: An Open Benchmark Dataset and Fusion Pipeline for Perception with Vehicle-to-Vehicle Communication. Proc., 2022 International Conference on Robotics and Automation (ICRA), IEEE, New York, 2022, pp. 2583–2589.

24.

Kini

R. R.

Sensor Position Optimization for Multiple LiDARs in Autonomous Vehicles. https://www.diva-portal.org/smash/record.jsf?pid=diva2%3A1525703&dswid=7752.

25.

Gouda

Mirza

Weiß

Ribeiro Castro

El-Basyouny

Octree-Based Point Cloud Simulation to Assess the Readiness of Highway Infrastructure for Autonomous Vehicles. Computer-Aided Civil and Infrastructure Engineering, Vol. 36, No. 7, 2021, pp. 922–940.

26.

Gouda

Pawliuk

Mirza

El-Basyouny

Using Convex Hulls with Octree/Voxel Representations of Point Clouds to Assess Road and Roadside Geometric Design for Automated Vehicles. Automation in Construction, Vol. 154, 2023, p. 104967.

27.

Gouda

Chowdhury

Weiß

Epp

El-Basyouny

Automated Assessment of Infrastructure Preparedness for Autonomous Vehicles. Automation in Construction, Vol. 129, 2021, p. 103820.

28.

Gouda

El-Basyouny

Evaluation of Traffic Sign Occlusion for Autonomous Vehicles Using Point Cloud Data. Transportation Research Record, 2024. 0(0). https://doi.org/10.1177/03611981241255359.

29.

Vehicle Icons in Figure 1–4 and 6: flaticon.com.

30.

AEVEX. DEM, DSM & DTM: Digital Elevation Model—Why It’s Important. https://geodetics.com/dem-dsm-dtm-digital-elevation-models/. Accessed June 11, 2023.

31.

Gao

Han

Implementing the Nelder–Mead Simplex Algorithm with Adaptive Parameters. Computational Optimization and Applications, Vol. 51, No. 1, 2012, pp. 259–277.

32.

Keys

Cubic Convolution Interpolation for Digital Image Processing. IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. 29, No. 6, 1981, pp. 1153–1160.

33.

Babaud

Witkin

A. P.

Baudin

Duda

R. O.

Uniqueness of the Gaussian Kernel for Scale-Space Filtering. IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. No. 1, 1986, pp. 26–33.

34.

Comer

M. L.

Delp

E. J.

III . Morphological Operations for Color Image Processing. Journal of Electronic Imaging, Vol. 8, No. 3, 1999, pp. 279–289.

35.

LearnOpenCV. Image Rotation and Translation Using OpenCV. https://learnopencv.com/image-rotation-and-translation-using-opencv/.

36.

Bresenham

J. E.

Algorithm for Computer Control of a Digital Plotter. Proc., Seminal Graphics: Pioneering Efforts That Shaped the Field, Association for Computing Machinery, New York, July 1, 1998, pp. 1–6.

37.

SkyeBrowse. Public Safety Drone Mapping Software. https://www.skyebrowse.com.

38.

BlueCity. https://ouster.com/products/software/bluecity.

39.

American Association of State Highway and Transportation Officials (AASHTO). A Policy on Geometric Design of Highways and Streets, 7th ed. American Association of State Highway and Transportation Officials (AASHTO), Washington, D.C., 2018.