Modeling the Use of LiDAR through Adverse Weather

Abstract

Due to the outstanding characteristics of LiDAR imaging systems, they seem essential for the consolidation of novel applications related to computer vision, in fields such as autonomous vehicles, outdoor recognition, and surveillance. However, the final technology implementation still has some uncertainties and needs in-depth work for its use in these real-world applications. Under the presence of adverse weather conditions, for example in fog, LiDAR performance is heavily influenced and the quality of the detection becomes severely degraded. The range is reduced due to the dispersion of the media and the sensor could be saturated due to backscattering or deliver a very limited range. Light propagation modeling through turbid media is used as a tool to understand and study these phenomena. Mie Theory allows the characterization of the optical media and light-particle interactions. Monte-Carlo methods are used to solve the radiative transfer problem related to these situations. When working with those models, the results obtained are in accordance with the ones shown in experimental tests, and it is possible to predict the necessities and problems of the designed systems.

1. Introduction

One of the up-to-date new applications of LiDAR technology is its use in transport, surveillance, and security [1]. A pulsed laser is used to measure ranges (variable distances). Thanks to distance calculations, it is possible to perform 3D mapping using different approaches [2] and to recover the geometry of the scene and not just a projection as in a conventional camera.

LiDAR technology has been presented as a disruptive technology regarding computer vision, as it gives precise and real-time visualization of the surrounding area and its distribution of objects. It offers an outstanding performance toward the required specifications [3]. However, the feasibility of the instrument for day-to-day outdoor uses is still facing numerous questions, one of them related to its detection breakdown when working in adverse weather conditions [4, 5, 6].

The propagation of light through scattering media such as fog, rain, smoke, dust, or others shows two main problems: the attenuation of the pulse and the saturation of the sensor [7, 8, 9]. The first one is a result of the dispersive effect and absorption characteristics of the media, which leads to the loss of energy while pulse light is propagating. The second is related to the backscattering that light undergoes when just entering the media, which may blind the sensor [10]. Due to their nature, both phenomena can be easily studied with models if the optical properties of the media are known. When working with models, the results obtained are in accordance with the ones shown in experimental tests.

At this point, it is worth showing in Figure 1 a point cloud obtained under fog conditions. The corresponding RGB image of the scene without fog is shown along three different views of the point cloud (with fog). Objects found in the scene are indicated on the RGB image with red letters. All point clouds are labeled as in the RGB image to facilitate its interpretation, i.e., the same letters are used inside black tags on the point clouds to point to objects. Z is the direction of propagation, normal to the RGB image and along the tunnel forming the fog chamber, Y refers to the height of the chamber and X to the width, distances are shown in meters [m]. On the left, we present a 3D view selected with an adequate orientation to highlight the aspect of the point cloud. On top, there is a YZ view, i.e., a side view of the scene; on the right there is a ZX view, which is a top view of the scene. These views are useful to notice the spread of the point cloud around a determined object/distance, and the points appearing due to the backscattering of fog, especially just in front of the sensor. In conclusion, the point cloud is rather noisy and the range is limited.

Figure 1.

Different views of a point cloud under the presence of artificial fog. The scene used is shown in the RGB image without fog.

Nowadays, there is no current solution to overcome the problem presented. The described effects are still present and seriously damage the performance of LiDAR systems in adverse weather conditions. However, there are several lines of research trying to find ways to improve this situation [11, 12, 13, 14, 15, 16]. The challenge of this topic makes researchers work with the novelties in optical engineering (optical design, materials of the components, new sensors…) along with the basis of physics of light (propagation of light, light-matter interactions…) [1].

In this chapter, we want to review some basics of the physics of light to properly become aware of the problem. How does light propagate through any turbid media? How is the media characterized? Which is its effect when working with pulsed light? If one wants to face this problem with a plausible solution, it is necessary to know and understand the involved physical phenomena in depth. Modeling allows us to go deeper into what is happening. Thus, it will be also reviewed how models are conventionally approached and which is the state of the art on the topic.

2. Optical properties

2.1 Description of the optical properties

For the description of the physical phenomenon of propagation of a pulse of light through a turbid media (energy transfer), the balance of incoming, outgoing, absorbed, and emitted photons is used. This is a wide, established, and well-known field, with a dense literature corpus [17, 18, 19, 20, 21, 22, 23].

Adverse weather conditions can be thought of as turbid media, as particles of different types and shapes (water, smoke, dust…) are found suspended within the main media, which is air [24]. By definition, a turbid medium is characterized by having localized non-uniformities randomly distributed within it. These optical non-uniformities are usually inclusions of one substance within another with a different index of refraction $n$. These inclusions cause the medium to be optically inhomogeneous and cause it to behave as a scattering media.

In the air, which is a non-absorbing media, the inclusions, widely represented as “particles,” are in charge of the actual absorption of part of the propagating energy and the actual change of direction of light. As a result, the dominant effects in the medium are absorption and scattering. According to these two effects, the medium is represented by several key optical parameters: the absorption coefficient ${\mu }_a$, the scattering coefficient ${\mu }_s$, the scattering phase function $p\left(\theta ,\varphi \right)$, and the asymmetry factor $g$; which respectively describe: the absorbing and scattering power, the probability of scattering in a particular direction $\left(\theta ,\varphi \right)$ of the media, and the degree of scattering in the forward direction.

An absorbing medium is composed of particles (or other structures) that can absorb light and transform it into its internal energy as the beam is propagated along the medium, which results in a gradual reduction of the light intensity. To characterize this phenomenon, one uses a parameter that is called the absorption coefficient ${\mu }_a$ (in units of ${\mathit{\text{length}}}^{-1}$, usually in ${\mathit{cm}}^{-1})$, which quantifies the absorbing effect of the medium.

The deviation of light from its straight trajectory due to localized non-uniformities in the medium is known as scattering. The particles become scattering centers, which, when exposed to light, modify the electromagnetic field and re-emit it in a different direction. Analogously to the absorption case, a parameter called the scattering coefficient ${\mu }_s$ (units of ${\mathit{\text{length}}}^{-1}$, typically $c{m}^{-1})$ is defined, which quantifies the scattering effect.

When a collimated beam of light passes through a volume of the medium, it will lose intensity due to both processes: absorption and scattering. In general, both processes cannot be distinguished. This effect is characterized by what is known as the extinction or attenuation coefficient ${\mu }_t={\mu }_a+{\mu }_s$ (${\mathit{cm}}^{-1})$. The extinction coefficient measures the total loss of a narrow-beam intensity, i.e. the loss due to absorption and the loss corresponding to the part of photons that have not been scattered in the forward direction.

Along with the scattering coefficient, a scattering event needs other parameters to be completely defined. In particular, photons may not be isotropically scattered and may need to have this dependence characterized. To completely define the deflection of the trajectory in space after a scattering event, two angles are used. The deflection scattering angle $\theta$ (which ranges from 0 to $\pi$) defines the deflection of the trajectory in the scattering plane—the plane formed by the direction of the incident light and the direction of the outgoing scattered light, i.e., the cone angle; and the azimuthal angle $\mathrm{\varphi }$(which is defined from 0 to $2\pi$) defines the change in the plane perpendicular to the scattering plane. These angles are shown in Figure 2, where the geometry of a scattering event is schematically depicted.

Figure 2.

Schematics of the scattering geometry.

The directionality of the scattering effect is quantified using a phase function $p\left(\mathrm{\theta },\mathrm{\varphi }\right)$. The phase function corresponds to the angular distribution of the light scattered by a scattering center at a given wavelength. It can be thought of as a probability density function, showing the chances of a photon being scattered in a particular direction.

When the suspended particles in media have no preferred scattering orientation (spherical symmetry), it is known as an isotropic medium. Then, light is scattered equally in all directions. However, it is more frequent that natural materials scatter light preferentially in the backward or forward direction. For those non-isotropic cases, it is interesting to know the amount of energy retained in the forward/backward direction after a single scattering event. If a photon is scattered so that its trajectory is deflected by a certain deflection angle $\theta$, then the component of the new trajectory, which is aligned in the former forward direction, is presented as$cos\left(\mathrm{\theta }\right)$. The mean value of this cosine is known as the anisotropy factor $g$, and it is defined as:

$$g\equiv ⟨\mathit{cos}\mathrm{\theta }⟩={\int }_0^{\pi }p\left(\mathrm{\Theta }\right)\mathit{cos}\mathrm{\Theta }·2\mathit{\text{πsin}}\mathrm{\Theta }d\mathrm{\Theta }$$E1

Its value varies in the range from $-1$ (total backward scattering) to $1$ (total forward scattering), being $g=0$ the value corresponding to isotropic scattering.

In addition, scattering may be elastic and inelastic. Without entering into many details, elastic scattering is associated with an interaction with no energy losses (and, thus, no wavelength change), while inelastic scattering corresponds to a process with energy transfer and thus to a wavelength shift. Generally speaking, elastic scattering is the predominant effect when propagating through turbid media, as approximately only one in each ${10}^7$ photons is inelastically scattered. Elastic interactions between photons and scattering particles are mainly described using two physical models: the Rayleigh and the Mie theories. The description of the process using one or another model is linked to the particle size and the wavelength of the incident light. However, Mie Theory is a general model developed using Maxwell equations and gives exact solutions in all cases, which means that it could be valid for any particle size.

Generally speaking, a turbid medium is described as a system of discrete spherical particles suspended within a base medium, which is exact for the case of fog (water suspended in air) and an approximation for the rest of the cases mentioned for the atmosphere (smoke…). Such spherical particles enable Mie Theory to be used for its characterization. This theory gives quantitative results of the interaction of an electromagnetic plane wave with a single homogeneous sphere, being likely the most important exactly soluble problem in the theory of absorption and scattering by small particles. Mie Theory allows the calculation of the absorption and scattering coefficients and the phase function of a spherical particle of radius $a$ as a function of the incident radiation wavelength $\lambda$, the size parameter $x=2\mathit{\pi a}/\lambda$, and the complex refractive indexes of the particles and the host material. The derivation of the complete theory may be long and tedious, and detailed information can be found in [19, 25].

2.2 Mie theory

The interaction of light with a spherical particle can be described and quantified using Mie Theory. Some conditions, however, have to be fulfilled to apply the theory. It has to be supposed that the media is homogeneous and that the particles that are embedded within it are spherical, homogeneous, and act independently—so they are distant enough from each other to consider only far-field scattering effects. Their radius and refractive index need to be known. As with most problems in theoretical optics, the scattering of light by a homogeneous sphere is treated as a formal problem of Maxwell’s theory with the appropriate boundary conditions [19, 25].

Suppose that one or more particles are placed in a beam of electromagnetic (EM) radiation. The rate at which EM energy is received by a detector downstream from the particles is denoted by $U$. If the particles are removed, the power received by the detector is $U_0$, where $U_0>U$ . We say that the presence of the particles has resulted in the extinction of the incident beam. If the medium in which the particles are embedded is non-absorbing (such as air), the difference $U_0-U$ accounts for absorption and scattering by the embedded particles (water droplets). Although the specific details of extinction depend on many parameters, certain general features are shared in common by all particles.

Now, consider extinction by a single arbitrary particle embedded in a non-absorbing medium and illuminated by a plane wave. If an imaginary sphere of radius r is constructed around the particle, the net rate at which EM energy crosses the surface $A$ of the sphere is $W_A$. If $W_A<0$, energy is absorbed within the sphere (being $W_A$ the rate at which energy is absorbed by the particle).

$W_A$ maybe conveniently written as the sum of:

$$W_A=-W_S+W_{\mathit{ext}}$$E2

$W_S$ is the rate at which energy is scattered across the surface $A$, and, therefore, $W_{\mathit{ext}}$ is just the sum of the energy absorption rate and the energy scattering rate: $W_{\mathit{ext}}=W_S+W_A$.

Now it is possible to define $C_x$ as the ratio of $W_x$ (being $x$: $A$, $S$ or $\mathit{ext}$) to $I_i$ (incident irradiance):

$$C_x=\frac{W_x}{I_i}$$E3

The $C_x$ quantities are called cross sections of the particle, and they have area dimensions. Let the total energy scattered in all directions be equal to the energy of the incident wave falling on the area $C_S$; likewise, the energy absorbed inside the particle may be defined as the energy incident in the area $C_A$, and the energy removed from the original beam may be equal to the energy incident in the area $C_{\mathit{ext}}$, which gives an idea of the amount of energy removed from the incident field due to scattering and/or absorption generated by the particle. The law of conservation of energy then requires that:

$$C_{\mathit{ext}}=C_A+C_S$$E4

When solving Maxwell’s equations for the defined problem, the scattered EM field is written as an infinite series in the vector spherical harmonics $M_n$ and $N_n$, which are the EM normal modes of the spherical particle. Thus, the scattered field is expressed as a superposition of these normal modes, each weighted by the appropriate coefficient $a_n$ or $b_n$, known as scattering coefficients.

It is found that:

$$C_S=\frac{W_s}{I_i}=\frac{2\pi }{k^2}\sum _{n=1}^{\infty }\left(2n+1\right)\left({\left|a_n\right|}^2+{\left|b_n\right|}^2\right)$$E5

$$C_{\mathit{ext}}=\frac{W_{\mathit{ext}}}{I_i}=\frac{2\pi }{k^2}\sum _{n=1}^{\infty }\left(2n+1\right)\mathit{Re}\left\{a_n+b_n\right\}$$E6

$$C_A=C_{\mathit{ext}}-C_S$$E7

Assuming that the series expansion of the scattered field is uniformly convergent, it is proved that the series can be terminated after:

$$n=x+4x^{1/3}+2$$E8

where $x=\mathit{ka}$ is the size parameter, $\text{with}k\text{being}$the wavenumber and $a$ the radius of the sphere.

We also need to obtain explicit expressions for the scattering coefficients:

$$a_n=\frac{S_n^{\prime }\left(y\right)S_n\left(x\right)-m{S}_n\left(y\right)S_n^{\prime }\left(x\right)}{S_n^{\prime }\left(y\right){\zeta }_n\left(x\right)-m{S}_n\left(y\right){\zeta }_n^{\prime }\left(x\right)}$$E9

$$b_n=\frac{m{S}_n^{\prime }\left(y\right)S_n\left(x\right)-S_n\left(y\right)S_n^{\prime }\left(x\right)}{{\mathit{mS}}_n^{\prime }\left(y\right){\zeta }_n\left(x\right)-S_n\left(y\right){\zeta }_n^{\prime }\left(x\right)}$$E10

where:

$$S_n\left(z\right)=\sqrt{\frac{\mathit{\pi z}}{2}}{J}_{n+0.5}\left(z\right)$$E11

$${\zeta }_n\left(z\right)=\sqrt{\frac{\mathit{\pi z}}{2}}{H}_{n+0.5}^{\left(2\right)}\left(z\right)$$E12

With $J_{n+0.5}\left(z\right)$ being the half-integral-order spherical Bessel function of first kind and $H_{n+0.5}^{\left(2\right)}\left(z\right)$ the half-integral-order Hänkel function of the second kind. The variables $x$ and $y$, in this case, correspond to $x=\mathit{ka}$ and $y=\mathit{mka}$, and $m$ is the relative refractive index between the sphere and the medium in which it is embedded; and finally, $S_n^{\prime }\left(z\right)$ and ${\zeta }_n^{\prime }\left(z\right)$ denote the derivatives of the corresponding functions.

Once the cross section of a single interaction has been computed, one needs to characterize the media. ${\mu }_a$ can be understood as the sum of contributions of the absorption cross sections of the absorbers per unit volume, i.e., the product of the absorption cross section $C_A$($c{m}^2$) by the density of absorbers ${\rho }_a$ ($\#/c{m}^{-3}$):

$${\mu }_a=C_A{\rho }_a$$E13

Analogously, knowing the number of scattering particles per unit volume ${\rho }_s$ ($\#/c{m}^{-3}$) and their $C_S,$ it is possible to compute the scattering coefficient ${\mu }_s$ of the propagating medium as:

$${\mu }_s=C_S{\rho }_s$$E14

The phase function and the asymmetry factor can also be computed using Mie Theory:

$$p\left(\theta \right)=\frac{2\pi }{k^2{C}_S}\left({\left|S_1\right|}^2+{\left|S_2\right|}^2\right)$$E15

$$g=\frac{2\pi }{k^2{C}_S}\left(\sum _{n=1}^{\infty }\frac{2n+1}{n\left(n+1\right)}\mathit{Re}\left\{a_n{b}_n^{\ast }\right\}+\sum _{n=1}^{\infty }\frac{n\left(n+2\right)}{n+1}\mathit{Re}\left\{a_n{a}_{n+1}^{\ast }+b_n{b}_{n+1}^{\ast }\right\}\right)$$E16

being $S_1$ and $S_2$:

$$S_1=\sum _{n=1}^{\infty }\frac{2n+1}{n\left(n+1\right)}\left(a_n{\pi }_n\mathit{cos}\mathrm{\theta }+b_n{\tau }_n\mathit{cos}\mathrm{\theta }\right)$$E17

$$S_2=\sum _{n=1}^{\infty }\frac{2n+1}{n\left(n+1\right)}\left(b_n{\pi }_n\mathit{cos}\mathrm{\theta }+a_n{\tau }_n\mathit{cos}\mathrm{\theta }\right)$$E18

From which $a_n$ and $b_n$ are the coefficients computed in Eqs. (9) and (10), and ${\pi }_n$ and ${\tau }_n$ two angle-dependent functions known as Mie angular functions. These angular functions are generated with the associated Legendre polynomials, and they can be calculated from the recurrence relations:

$${\pi }_n=\frac{\left(2n-1\right)\mu }{n-1}{\pi }_{n-1}-\frac{n}{n-1}{\pi }_{n-2}$$E19

$${\tau }_n=\mathit{n\mu }{\pi }_n-\left(n+1\right){\pi }_{n-1}$$E20

where$\mathrm{\mu }=\text{cosθ}$, and the first terms of π_n are ${\mathrm{\pi }}_0=0$ and ${\mathrm{\pi }}_1=1$.

3. Propagation model

Propagation through turbid media is characterized by multiple scattering. This is the situation in which along the total traveling distance $L$ of light, each photon undertakes on average many collisions with particles of the medium. Under these conditions, ${\mu }_{s}L\ge 1$ is fulfilled and light properties such as phase, polarization, and ray trace are severely degraded, so only radiometric information may be considered. The basic theory that allows the calculation of light distributions in multiple scattering media with absorption is the radiation transfer theory (RTT). Its core is the radiation transfer equation (RTE)—linear transport or Boltzmann equation, a balance equation characterizing the flow of photons (or any particle) in a given volume element [17].

Some examples of multiple scattering media are biological tissue, nebulous media, colloids, murky water, clouds, and also adverse weather such as fog and rain. For adverse weather conditions, particles suspended in the air act as light scatterers, being responsible for the reduction of the visibility near the ground surface. By using Mie Theory and knowing the approximate distribution of the size particles [26] (for example, modified gamma distributions are the standard choice for fitting the models of fog droplet size distributions [27]), it is possible to give values to the optical parameters (${\mu }_a$, ${\mu }_s$, $p\left(\mathrm{\theta },\mathrm{\varphi }\right)$, and $g$) and characterize the medium for solving the RTE.

Once the media has been characterized, the RTE has to be solved to describe the energetic distribution of light. However, the analytical solution to this equation is complicated and often impossible to solve when boundaries, inhomogeneities, or nonstationary effects are involved. The Monte-Carlo (MC) method is used as the conventional tool to arrive at a statistical solution in these cases. MC is a stochastic method, which offers robustness and versatility for solving this kind of problem. By modeling, we can predict the shape, range, and intensity of a pulsed LiDAR working through turbid media.

The MC method refers to a technique first proposed by Metropholis and Ulam [28] to simulate physical processes using a stochastic model. Regarding the radiative transfer problem, the MC method is based on recording photons’ histories as they are scattered and absorbed, using the global optical properties of the medium [29, 30]. Simulations show the expected movement of individual photons, which are treated as particles of light that move according to certain probability density functions. The photon moves in a straight path and may come across obstacles. At the surface of the obstacles, it may undergo absorption or scattering. Then, the photon continues its flight until it is absorbed or leaves the medium (some examples are shown in Refs. [31, 32, 33]). When this process is repeated for many individual photons, MC simulations provide a flexible approach towards light transport that yields maps of the light distributions in turbid media induced by a light source.

MC solutions can be obtained for any desired accuracy, which is proportional to $1/\sqrt{N}$ where $N$ is the number of photons propagated. Thus, relative errors less than a few tenths of a percent will require the propagation of substantial numbers of photons (between ${10}^6$ and ${10}^9$ photons) and may require large amounts of computer time [33].

3.1 Model structure

The modeling of pulsed light propagation is based on solving the RTE using MC method [29, 30, 31, 32, 33, 34, 35]. The scheme shown in Figure 3 summarizes the whole code structure:

Figure 3.

Flux diagram of the MC method used.

Next, we will briefly describe the main steps of the MC method.

3.1.1 Photon initialization

There are different MC approaches. Variance reduction techniques are used to reduce the number of photons necessary to achieve the desired accuracy for a Monte-Carlo calculation. One of them is implicit capture. In this approach, to improve the efficiency of the MC program, many photons (a packet) are propagated along each pathway.

The MC method thus begins by launching a packet of photons, with a size called weight, into the medium. Usually, one may think that only one photon follows each pathway, and at each step, the photon may be either absorbed or scattered. The packet approach is used to improve the efficiency of the MC program, as many photons (a packet) are propagated along each pathway at the same time. If a packet of photons followed each pathway, then some portion of the packet would be absorbed in each step. So, after each propagation step, w is reduced by the probability of absorption.

When a collimated beam normally incident on a slab vertically is simulated, the packets’ (which from now on will be called photons) initial direction is chosen downward into the medium, orthogonal to the surface. The initial coordinates of the photon are usually identical for all photons and the weight of the photon is initially set to unity.

3.1.2 Propagation distance

Once launched, the photon moves a distance ∆s. The most efficient method is to choose a different step-size $∆s$ for each photon step. The probability density function for the step size follows Beer’s law, so the probability of scattering is proportional to $e^{-{\mu }_t\mathrm{\Delta }s}$, being ${\mu }_t={\mu }_a+{\mu }_s$ - and is chosen in such a way that it is the distance at which the photon interacts with the obstacle. The $∆s$ value may be generated with this probability density function as a function of a random number $\xi$, uniformly distributed between 0 and 1:

$$\mathit{\Delta s}=\frac{-ln\xi }{{\mu }_a+{\mu }_s}$$E21

3.1.3 Moving the photon

A photon is uniquely described by five variables: three spatial coordinates for the position and two directional angles for the direction of travel. However, it is convenient to describe the photon’s spatial position with three Cartesian coordinates ($x;y;z)$ and the direction of travel with three direction cosines ($u_x,u_y,u_z$), corresponding to the cosine of the angle that the photon’s direction makes with each axis: X-, Y-, and Z- axis respectively. In this case, the required formulas for propagation are simpler. For a photon located at ($x;y;z$) traveling a distance ∆s in the direction ($u_x,u_y,u_z$), the new coordinates $(x^{’};y^{’};z^{’}$) are given by:

$${\mathrm{x}}^{\prime }=\mathrm{x}+{\mathrm{u}}_{\mathrm{x}}∆\mathrm{s}$$

$${\mathrm{y}}^{\prime }=\mathrm{y}+{\mathrm{u}}_{\mathrm{y}}∆\mathrm{s}$$E22

$${\mathrm{z}}^{\prime }=\mathrm{z}+{\mathrm{u}}_{\mathrm{z}}∆\mathrm{s}$$

3.1.4 Absorption

There exist different methods to consider absorption during propagation. The most followed approach is known as the technique of implicit capture. It assigns a weight to each photon as it enters the medium. After each propagation step, the photon is split into two parts: a fraction is absorbed, and the rest is scattered. Thus, in every encounter, the photon interacts with the scatters. Upon interaction, a fraction ${\mu }_a/{\mu }_t$—that corresponds to the probability of absorption—if the photon’s weight is absorbed and the remaining ${\mu }_s/{\mu }_t$ fraction of the photon’s weight is scattered and continues to propagate. The absorbed fraction is placed in a bin of the MC absorption matrix ($\mathit{W})$, which encloses the current photon position:

$$\mathit{W}\left(x,\mathrm{y},\mathrm{z}\right)=\frac{{\mu }_a}{{\mu }_t}w$$E23

and the new photon’s weight ($w^{\prime })$ is updated:

$$w^{\prime }=\frac{{\mu }_s}{{\mu }_t}w$$E24

As the weight of the photon falls below a certain threshold, a photon termination strategy (see 3.1.6) needs to be implemented.

3.1.5 Scattering

The angle at which the photon is bent when it strikes an obstacle is defined by the two angles of scattering ($\theta$ and $\varphi$) and the normalized phase function, which describes the probability density function for the angles at which a photon is scattered. As has been studied in previous sections, each type of obstacle is characterized by a different phase function.

For isotropic scatterers, such as spherical particle approximation, the phase function has no azimuthal dependence and only depends on $\theta$. Thus, $\varphi$ is uniformly distributed between 0 and $2\pi$ and may be generated by multiplying a random number uniformly distributed over the interval 0 to 1 by $2\pi$. Using Mie Theory, it is possible to compute the phase function as shown in Eq. (15), which then has to be sampled to compute $\theta$.

Nevertheless, an approximation with a lower computational cost is close enough for most cases. In practice, it is more convenient to use semiempirical approximations of the scattering phase function, with much lower computational cost. This choice is a compromise between realism and mathematical tractability. The most common example is the Henyey-Greenstein (HG) phase function, which includes the average cosine $g$:

$$p^{\mathit{HG}}\left(\mathrm{\theta }\right)=\frac{1}{4\pi }\frac{1-g^2}{{\left(1+g^2-2g\mathit{cos}\mathrm{\theta }\right)}^{3/2}}$$E25

HG is an analytical function originally derived for modeling scattering by interstellar dust and is widely used in biomedical optics and other fields. Then, one only has to invert the probability density function to obtain the generating function of $\mathrm{\theta }$:

$$\mathit{cos}\mathrm{\theta }=\frac{1}{2g}\left[1+g^2-{\left(\frac{1-g^2}{1+g+2\mathit{g\xi }}\right)}^2\right]$$E26

where $\xi$ is a random number uniformly distributed between 0 and 1.

If a photon is scattered at an angle ($\mathrm{\theta },\mathrm{\varphi }$) from the initial direction $(u_x,u_y,u_z$) in which it is traveling, then the new direction (${u_x}^{\prime },{u_y}^{\prime },{u_z}^{\prime }$) is specified by:

$$u_x^{\prime }=\frac{\mathit{sin}\mathrm{\theta }}{\sqrt{1-u_z^2}}\left(u_x{u}_z\mathit{cos}\mathrm{\varphi }-u_y\mathit{sin}\mathrm{\varphi }\right)+u_x\mathit{cos}\mathrm{\theta }$$E27

$$u_y^{\prime }=\frac{\mathit{sin}\mathrm{\theta }}{\sqrt{1-u_z^2}}\left(u_y{u}_z\mathit{cos}\mathrm{\varphi }+u_x\mathit{sin}\mathrm{\varphi }\right)+u_y\mathit{cos}\mathrm{\theta }$$E28

$$u_z^{\prime }=-\mathit{sin}\text{θcosϕ}\sqrt{1-u_z^2}+u_z\mathit{cos}\mathrm{\theta }$$E29

3.1.6 Photon termination

A photon will be terminated either if it exits the considered space or if it is absorbed. However, using the technique of implicit capture, photon weight will never be a mathematical 0. Thus, a photon is terminated when its weight falls below a given threshold, despite the termination of the process using a threshold introduces a systematic negative bias into the system regarding energy conservation. To reduce this bias, the Roulette method is used [36].

The Roulette method works in the following way. A predefined number $N$ between 2 and 10 is usually chosen in practice. Once the weight of the photon reduces below a sufficiently small threshold, a uniform random number $\xi$ between 0 and 1 is generated. The photon is removed from the system only if $\xi >1/N$. The photon that survives is continued with its current weight increased by a factor of $N$. The result is that photons are usually terminated, but energy is conserved by the occasional surviving photon being given extra weight. Since millions of photons are run, the statistically averaged result is correct.

3.1.7 Observable

Once all the photons have been run, the data are stored in the bin of $\mathit{W}$, as $W\left(x,y,z\right)$, in units of weight/bin. In order to obtain physical values, some changes are needed. The results will be delivered in the form of fluence rate $\mathrm{\Phi }$ [$\mathit{\text{Watts}}/c{m}^2$].

Firstly, we need to convert from weight to photons. Then we normalize it by the appropriate voxel volume ($V$) and the total number of photons, which leads to the photon absorption density in units of $c{m}^{-3}$. Finally, the fluence rate $\mathrm{\Phi }$ [Watts/cm²] is obtained by dividing the power density absorbed by the absorption coefficient and multiplying it by the incident power $P$ [$\mathit{\text{Watts}}]$:

$$\mathrm{\Phi }\left(x,y,z\right)=\frac{W\left(x,y,z\right)·N_p·P}{V·N_t·{\mu }_a}$$E30

In which $N_p$ is the number of photons per packet and $N_t$ the total number of photons.

3.1.8 Temporal approach

Given that pulsed LiDAR is time-dependent, we are interested in time-resolved simulations. MC propagating code is easily modified by adding a record of the time history of each photon in order to have the time-tracking of light propagation. By using this history, it is possible to generate the temporal profile of optical power. The time that each package remains in space is obtained by dividing the length of the path traveled by the speed of light in that medium [30, 37].

In summary, modeling the radiative transfer problem in scattering media involves a potentially simple methodology: a photon packet is emitted, it travels a distance, and something happens to it that affects its energetic weight in successive events. MC method gives us a statistical solution, which becomes a powerful tool to design and characterize our systems.

4. Applications

Models provide a tool to predict light behavior in different situations. There exist two basic analytical models that correspond to two particular situations in which the RTE has an analytical solution. The solutions can be used to compare the results obtained with stochastic models, such as MC, and verify that MC can be also applied to solve RTT without solving the RTE.

On the one hand, it is well known that when absorption prevails over scattering (${\mu }_a\gg {\mu }_s$), Beer-Lambert law can be applied [22]. On the other hand, under certain conditions of the diffusion regime ($10{\mu }_a\approx {\mu }_s$), in source-free and homogeneous media conditions, there is an analytical solution based on the theory of photon density wave in steady conditions [38]. Figure 4 shows, in both cases, results obtained from MC simulations (red dots) and the analytical solution (blue line). They are superimposed showing the validity of our implementation of the MC approach.

Figure 4.

Left: Beer-Lambert Law. Propagation of light through a medium with ${\mu }_a=10c{m}^{-1},{\mu }_s=0.05c{m}^{-1},g=0.9$. Right: Diffusion regime. Propagation of light through a medium with ${\mu }_a=0.05c{m}^{-1},{\mu }_s=20c{m}^{-1},g=0.2$.

As discussed, if one wants to study pulsed light propagation of a LiDAR system, time-resolved simulations have to be used. The MC code needs to be modified by adding a record of the time history of each photon to have the time-tracking of light propagation. By using this history, it is possible to generate the temporal profile of optical power. For this case, there also exist situations in which there is an analytical solution to the RTE [39]. For example, in the case of an infinite diffusing medium with a point source, as shown in Figure 5, the setting between the analytical and the simulated solution is almost perfect.

Figure 5.

Time-resolved propagation of light through a medium with ${\mu }_a=0c{m}^{-1},{\mu }_s=2c{m}^{-1},g=0.9.$

Finally, Figure 6 shows how models based on MC methods are able to predict pulsed light interactions, which is, in fact, the topic of interest in this chapter. It shows the temporal profile of a pulse of light sent of Gaussian shape and 1 mJ of energy (in orange) and light received back to the lighting (in blue). The medium is characterized by the following properties: ${\mu }_a=0c{m}^{-1},{\mu }_s=0.01c{m}^{-1},g=0.9$, which would correspond to an almost clear media, and a 100% reflective plane is placed at a distance of 0.5 m. It can be observed from this simple figure the working principle of LiDAR. The position of the reflected light peak, with respect to the initial pulse, allows us to obtain information on the position of the object. The shape of the pulse is maintained and the integral of both pulse profiles meets the expected energy of the propagation calculated from the range equation.

Figure 6.

Light signal sent and received through a medium with ${\mu }_s=0.01c{m}^{-1},g=0.9$.

One of the problems shown by LiDAR technology when working through adverse weather conditions is the saturation of the sensor. This saturation is the result of the backscattering that light undergoes when just entering the media, and it may saturate the sensor. Thus, to start with, it is interesting to analyze simple backscattering interactions in different media in order to have an idea of how the first light that would arrive at our detection in foggy conditions would influence our sensor. As the final objective is to use these MC methods to improve a long-range active sensing technique that interacts with a target, we also focused on studying the effect of reflected light that reaches back the plane of illumination, i.e., the light that comes back from an object hidden behind the turbid environment. Therefore, we want to know how the pulse may be attenuated and what the range of the system will be.

Figure 7a shows different profiles of backscattered light simulated using an MC code. The greater the scattering coefficient, the greater the amount of backscattered energy and the more extended in time. It can be also observed that the backscattering signal does not have the same shape as the initial pulse, which is defined as a perfect Gaussian-shaped pulse [40]. Figure 7b shows experimental data of a pulsed LiDAR working under different artificial fog densities, simulating different meteorological visibilities. One can see that stochastic models (Figure 7a) reproduce the expected behavior of light (Figure 7b).

Figure 7.

a) Detected backscattering in a medium with $g=0.9$ and different ${\mu }_s$, for a light pulse with an initial energy of 1 mJ. b) Experimental results of detected backscattering in a foggy medium for different visibilities.

Using the simulations and knowing the specifications of our sensor, one can adjust electronic gain and optical power to avoid saturation of it, or even come up with solutions related to the optical design of the system.

Next, in Figure 8 we present the temporal profiles under the same conditions as in Figure 7a, but with the presence of an object placed at 0.5 m with a reflectivity of 100%. If one supposes that the sensor is not saturated under any of the presented conditions, it is observed that around the expected position of the object appears the peak that corresponds to its presence along with the backscattering contribution. Using the simulations, it is possible to evaluate what is the level of scattering at which the target is no longer detected (it is not distinguishable from the backscattering signal), so it would not be possible to obtain its image or compute its distance.

Figure 8.

Detected backscattering in a medium with $g=0.9$ and different ${\mu }_s$, for a light pulse with an initial energy of 1 mJ, with the presence of a perfect reflecting target.

The same type of study can be performed to establish the range limit of the system taken into account between visibility and surface properties of the object. Data from an experimental test are shown in Figure 9, and it can be seen how both parameters may influence the results.

Figure 9.

Signal detected from a returning pulse of light pointing toward (LEFT) a metallic object at 12 m of distance and (RIGHT) a diffusive object at 14 m of distance, for two different visibilities (20 m and 75 m).

5. Conclusions

Over recent years, LIDAR technology has become a panacea in the fields of optomechanical engineering and optoelectronics. It especially seems to hold a relevant role in novel applications related to outdoor environments. One of the keys to its success is the amount of information it can provide despite relying on a very simple method. Being based on a simple working principle (counting elapsed time between events in magnitudes carried out by light, e.g., reflected energy from a pulse of light sent to a target), it allows the performance of complete 3D mapping with outstanding characteristics.

However, for the complete settling of the technology, there are still obstacles pending to be solved. One of the most challenging is related to its outdoor performance. As with any other optical sensor, under the presence of adverse weather conditions, such as fog, the system performance is heavily altered and the quality of the detection becomes severely degraded.

Usually, commercial systems are like black boxes, only returning the point cloud or 3D map. When facing bad weather, data are unreliable and its behavior is unknown. However, the LiDAR system can provide a lot more information, which might enable optimized features.

In order to propose reasonable solutions, the knowledge of optical physics behind the phenomena degrading the performance of the system may be of significant interest. The propagation of light through scattering media, such as adverse weather conditions, shows two main problems: the saturation of the sensor and the attenuation of the pulse. Both are related to the dispersive effect and absorption characteristics of the media. Due to their nature, both phenomena can be easily studied with models if the optical properties of the media are known.

Using Mie Theory, and taking into necessary conditions, we can obtain an approximation of the media properties (absorption coefficient, scattering coefficient, phase function, and anisotropy factor). The derivation of the complete theory may be long and tedious; however, its application is straightforward.

Once all the media properties are derived, they can be used to solve the radiation transfer theory. Essentially, it is the basic theory that allows the calculation of light distributions in multiple scattering media with absorption. Its core is the radiation transfer equation, which computes the balance that characterizes the flow of photons in a given volume element. However, the analytical solution to this equation is complicated and often impossible to solve when boundaries, inhomogeneities, or nonstationary effects are involved. The Monte-Carlo (MC) method is used as the conventional tool to arrive at a statistical solution in these cases. By using the MC method, we can solve the RTT of a pulsed LiDAR when working through a turbid medium. Along this chapter, we have shown how the model succeeds in modeling different kinds of scenarios: a media where absorption prevails over scattering, predictions of pulsed light interactions, dynamics of backscattering, light response to different kinds of objects and media, etc. As a result, we are able to predict the behavior of our system in the different scenarios where it breaks down. For example, it is possible to estimate the range of the system, the response toward different objects or the characteristics of the blinding backscattering.

The next step would be related to studying other features that are considered very relevant when designing LiDAR imaging systems, for example:

• Numerical characterization of the backscattering phenomena of active illuminators due to the first interaction with the environment.

• Scattering effects in propagation and its wavelength dependence.

• Optimal configuration sensing to improve contrast-to-noise ratio for imaging (CNR) and radiometric detection.

• Using simulations to find other properties that could be taken into account to discern between signal photons and backscattering photons (for example polarization).

Moreover, one has to not lose focus on the final implementation of the system. Simulations can be used to guide us, to try exotic ideas without the need to build up the system. However, experimental results are always expected as the final product of this whole process.

Acknowledgments

Authors thank Agència de Gestió d’Ajuts Universitaris i de Recerca for Grants 2021FI_B2 0007. Authors also acknowledge support received from DSTL under contract DSTLX1000145661 and the one from Misnitry Sicience and Innovation through Projects PID2020-119484RB-I00 and PDC2021-121038-I00.

Conflict of interest

The authors declare no conflict of interest.