**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, *μsL*≥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 (*μa*, *μs*, *p*ð Þ θ, ϕ , 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].

*Modeling the Use of LiDAR through Adverse Weather DOI: http://dx.doi.org/10.5772/intechopen.109079*

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–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*=* ffiffiffiffi *N* <sup>p</sup> where *<sup>N</sup>* 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<sup>6</sup> and 10<sup>9</sup> 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–35]. The scheme shown in **Figure 3** summarizes the whole code structure: 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*�*μt*Δ*<sup>s</sup>* , being *μ<sup>t</sup>* ¼ *μ<sup>a</sup>* þ *μ<sup>s</sup>* - 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 *ξ*, uniformly distributed between 0 and 1:

$$
\Delta s = \frac{-\ln \xi}{\mu\_a + \mu\_s} \tag{21}
$$

**Figure 3.** *Flux diagram of the MC method used.*

#### *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 (*ux*,*uy*,*uz*), corresponding to the cosine of the angle that the photon's direction makes with each axis: X-, Y-, and Zaxis 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 (*ux*,*uy*,*uz*), the new coordinates *x*' ; *y*' ; *z*' ) are given by:

$$\begin{aligned} \mathbf{x'} &= \mathbf{x} + \mathbf{u}\_{\mathbf{x}} \Delta \mathbf{s} \\ \mathbf{y'} &= \mathbf{y} + \mathbf{u}\_{\mathbf{y}} \Delta \mathbf{s} \\ \mathbf{z'} &= \mathbf{z} + \mathbf{u}\_{\mathbf{z}} \Delta \mathbf{s} \end{aligned} \tag{22}$$

#### *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 *μa=μt*—that corresponds to the probability of absorption—if the photon's weight is absorbed and the remaining *μs=μ<sup>t</sup>* 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 (*W*Þ, which encloses the current photon position:

$$\mathbf{W}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \frac{\mu\_a}{\mu\_t} w \tag{23}$$

and the new photon's weight (*w*<sup>0</sup> Þ is updated:

$$w' = \frac{\mu\_s}{\mu\_t} w\tag{24}$$

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 (*θ* and *ϕ*) 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 *θ*. Thus, *ϕ* is uniformly distributed between 0 and 2*π* and may be generated by multiplying a random number uniformly distributed over the interval 0 to 1 by 2*π*. Using Mie Theory, it is possible to compute the phase function as shown in Eq. (15), which then has to be sampled to compute *θ*.

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^{HG}(\theta) = \frac{1}{4\pi} \frac{1 - \text{g}^2}{\left(1 + \text{g}^2 - 2\text{g }\cos\theta\right)^{3/2}} \tag{25}$$

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 θ:

$$\cos \Theta = \frac{1}{2g} \left[ 1 + g^2 - \left( \frac{1 - g^2}{1 + g + 2g\xi} \right)^2 \right] \tag{26}$$

where *ξ* is a random number uniformly distributed between 0 and 1.

If a photon is scattered at an angle (θ,ϕ) from the initial direction *ux*,*uy*,*uz* � ) in which it is traveling, then the new direction (*ux* 0 ,*uy* 0 ,*uz* 0 ) is specified by:

$$u'\_{x} = \frac{\sin\Theta}{\sqrt{1 - u\_x^2}} \left( u\_x u\_x \cos\Phi - u\_y \sin\Phi \right) + u\_x \cos\Theta \tag{27}$$

$$u\_{\gamma}^{\prime} = \frac{\sin \Theta}{\sqrt{1 - u\_x^2}} \left( u\_{\gamma} u\_x \cos \Phi + u\_x \sin \Phi \right) + u\_{\gamma} \cos \Theta \tag{28}$$

$$u\_x' = -\sin\theta\cos\phi\sqrt{1-u\_x^2} + u\_x\cos\Theta \tag{29}$$

#### *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 *ξ* between 0 and 1 is generated. The photon is removed from the system only if *ξ*>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 *W***,** as *W x*ð Þ , *y*, *z* , 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 Φ [*Watts=cm*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 *cm*�3. Finally, the fluence rate Φ [Watts/cm<sup>2</sup> ] is obtained by dividing the power density absorbed by the absorption coefficient and multiplying it by the incident power *P* [*Watts*�:

$$\Phi(\mathbf{x}, \mathbf{y}, z) = \frac{W(\mathbf{x}, \mathbf{y}, z) \cdot N\_p \cdot P}{V \cdot N\_t \cdot \mu\_d} \tag{30}$$

In which *Np* is the number of photons per packet and *Nt* 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.
