We are IntechOpen,
the world’s leading publisher of
Open Access books
Built by scientists, for scientists
5,500
136,000
170M
Open access books available
International authors and editors
Downloads
Our authors are among the
154
TOP 1%
12.2%
Countries delivered to
most cited scientists
Contributors from top 500 universities
Selection of our books indexed in the Book Citation Index
in Web of Science™ Core Collection (BKCI)
Interested in publishing with us?
Contact book.department@intechopen.com
Numbers displayed above are based on latest data collected.
For more information visit www.intechopen.com
Chapter
Monte Carlo Radiative Transfer
Modeling of Underwater Channel
Rafael M.G. Kraemer, Luís M. Pessoa
and Henrique M. Salgado
Abstract
The radiative transfer equation (RTE) is a theoretical framework that can be
used for predicting and interpreting underwater light fields in terms of the constituents of natural water bodies. However, the RTE is a complex integrodifferential
equation and deriving exact solutions for it is a difficult task. In this chapter, we aim
to present some details regarding Monte Carlo simulations and how this method
may be applied to solve the RTE numerically. By solving the RTE, one may accurately predict the received power and estimate the channel bandwidth and several
other measurable parameters with regard to multiple water conditions. Simulations
will also be presented.
Keywords: RTE, underwater, optical, wireless, channel, propagation, photons,
UOWC
1. Introduction
When compared with free space optical (FSO) communication channel, the
underwater optical wireless communication (UOWC) shows unique characteristics.
Because of this, the channel models usually used in FSO are not suited for UOWC,
and different channel models must be developed to describe the different photon
interactions under play.
Mobley [1] classified the optical properties of water into two different groups,
inherent optical properties (IOPs) and apparent optical properties (AOPs). The two
major IOPs that will attenuate light propagation in UOWC are absorption and
scattering. Hence, any attempt in modeling the UOWC channel will need to consider absorption and scattering effects within specific link configurations.
Several theoretical models were employed by researchers to model the UOWC
channel, being the Beer-Lambert law the most employed due to its simplicity [2].
Despite its simplicity, the Beer-Lambert law contains two implicit assumptions: that
the transmitter and receiver are perfectly aligned and that all scattered photons are
lost. The assumption that all photons that undergo scattering are lost will severely
underestimate the received optical power, especially in water environments where
scattering is the dominant IOP.
Another theoretical model for modeling the underwater channel is the radiative
transfer equation (RTE) [3]. The RTE simply says that as a beam of radiation travels
through a medium, in this case water, it will lose energy to absorption, gain energy
by emission, and redistribute energy by scattering. Yet, as the RTE is an
1
Wireless Mesh Networks - Security, Architectures and Protocols
integrodifferential equation involving different independent variables, an exact
analytical solution is hard to find. In view of this, solving the RTE numerically is a
preferred approach, being the most popular one the Monte Carlo simulation.
2. The radiative transfer equation
One important law of geometrical radiometry is the n-squared law for radiance.
To derive it, we consider two mediums separated by a transparent surface; in that
case, Snell’s law states that the angles of an incident ray from medium 1 to medium 2
are related by:
n1 sin ðθ1 Þ ¼ n2 sin ðθ2 Þ,
(1)
where n is the index of refraction of the medium. Considering that a ray is
simply a narrow beam of photons traveling in almost the same direction, and if
we consider two narrow rays, incident and refracted, that will have solid angles
given by:
ΔΩ1 ¼ sin θ1 Δθ1 Δϕ,
(2)
ΔΩ2 ¼ sin θ2 Δθ2 Δϕ,
(3)
and as seen in [3], squaring each side of Eq. (1) and multiplying by the Azimuthal spread, Δϕ, we obtain:
n21 cos θ1 Ω1 ¼ n22 cos θ2 Ω2 ,
(4)
which is known as the Straubel’s invariant. Consider now the radiances of the
two rays, incident and refracted, defined as:
L1 ¼
ΔP1
,
ΔA1 ΔΩ2
(5)
L2 ¼
ΔP2
,
ΔA2 ΔΩ2
(6)
where ΔP1, 2 is the spectral radiant power and ΔA1, 2 is the cross section area of
incident and refracted rays. The ratio between the refracted and incident rays is
called the Fresnel transmittance:
ΔP1
¼ T,
ΔP2
(7)
L2
L1
¼T 2,
2
n2
n1
(8)
from which we can obtain:
which is the n-squared law for radiance. For the case of T ¼ 1, which is the case
of normal incidence on an air-water surface, Eq. (8) becomes:
L2 L1
¼ :
n22 n21
2
(9)
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
As seen in [3], this result is known as the fundamental theorem of radiometry
and it states that the radiance divided by the refraction index squared is constant
along any path; however, because all real substances will cause some absorption and
scattering on the incident photons, the validity of theorem is restricted for paths in
vacuum.
If we consider that a ray of radiance L1 starts in a medium with index of
refraction n1 and goes through successive refractions and ends up with radiance LS
in a medium of refractive index ns , we can come up with a formulation for successive crossings of a ray from one medium to another:
s 1
LS Y
n2
T ð j; j þ 1Þ s2 ,
¼
L1
n1
j¼1
(10)
where T ð j; j þ 1Þ is the Fresnel transmittance of the interface between the media
with refractive indexes nj and njþ1 . From this, we can see that the radiance along a
path will change due to variations in the real index of refraction along that same
path. The index of refraction in a water body will change from point to point by
random molecular motions, by organics or inorganic particulate matter and by
turbulent fluctuations in temperature and salinity.
!
!
!
Considering a beam of
traveling from point x to point x þΔ x and
radiance
!
taking the limit as Δr ¼ Δ x ! 0:
ΔðL=n2 Þ DðL=n2 Þ
¼
,
Δr!0
Δr
Dr
(11)
lim
being D=Dr the total rate of change along the path. The total rate of change can
be expressed in terms of the advective derivative, or the substantive derivative, as
referred by Mobley:
D
1D
¼
,
Dr v Dt
(12)
!
where v is the speed of light in the medium at position x at time t for
wavelength λ:
!
!
v x ; t; λ ¼ c=n x ; t; λ ,
(13)
D
∂ !
¼ þ v ∇,
Dt ∂t
(14)
The operator D=Dt is:
and the gradient operator ∇ is defined in the usual way.
From these developments, we can rewrite Eq. (11) as:
D
1D 1∂ ^
¼
¼
þ ξ ∇,
Dr v Dt v ∂t
(15)
!
valid for photons traveling with speed v in the direction ^ξ ¼ v =v:
Writing the expression for the change in L=n2 along a path that combines all the
physical phenomena causing that change, we obtain the most general form of the
RTE for unpolarized radiance:
3
Wireless Mesh Networks - Security, Architectures and Protocols
1∂ L
L
^
þξ∇ 2 ¼
2
v ∂t n
n
L
c 2 þ LE þ LI þ LS :
n
(16)
In Eq. (17), c is the beam attenuation coefficient, LE , LI , and LS are the path
functions for elastic scattering, inelastic scattering, and spontaneous emission,
respectively. The path function for spontaneous emission is also referred in [3] as
the source path function and can be considered as the contribution to total radiance
from a light source and in the case of UOWC systems a laser diode or a light
emitting diode (LED).
For UOWC systems, we are interested in the time-independent RTE in horizontally homogenous water bodies with a constant index of refraction. Because of this,
the factor n 2 will divide both sides of Eq. (15). Hence, with ∂L=∂t ¼ 0 and using the
notation adopted by Mobley where x3 ¼ z and ξ3 ¼ cosθ, we obtain:
^ξ ∇L ¼ ξ3 ∂L ¼ cosθ dLðθ; ϕÞ :
∂x3
dz
(17)
It is also usual [3] to combine LI and LS as an effective source function
S ¼ LI þ LS . The RTE then becomes:
cos θ
dLðθ; ϕÞ
¼
dz
cLðθ; ϕÞ þ LE ðθ; ϕÞ þ Sðθ; ϕÞ,
and the path function for elastic scattering, LE :
ð
E
Lðθ0 ; ϕ0 Þβðθ0 ; ϕ0 ! θ; ϕÞdΩðθ0 ; ϕ0 Þ,
L ðθ; ϕÞ ¼
(18)
(19)
4π
Combining all these terms in Eq. (18), we can come up with the standard form
of the RTE. As a simple statement, the IOPs and boundary conditions will go into
the RTE and the final radiance will come out as a result. The standard form of the
RTE, as shown in [4], is given by:
dLðθ; ϕÞ
¼
cosθ
cdz
Lðθ; ϕÞ þ ω0
ð
4π
e
β ðθ0 ; ϕ0 ! θ; ϕÞLðθ0 ; ϕ0 ÞdΩðθ0 ; ; ϕ0 Þ þ Sðθ; ϕÞ,
(20)
where Lðθ; ϕÞ represents the radiance in the direction ðθ; ϕÞ, ω0 is the single
scattering albedo, c is the beam attenuation coefficient, and e
β ðθ; ϕÞ is the scattering
phase function. Each one of these terms will be explained in more detail in the
following sections.
2.1 Inherent optical properties
During the interaction of a photon with a water molecule, one of two things may
happen: the photon may be absorbed, leaving the water molecule in a state with
higher internal energy, or the photon may undergo scattering. Such scattering
occurs when there is a change in the direction of the propagation of the photon or in
the photon energy—or in both.
When scattering happens, if the molecule that interacted with the photon
returns to its original internal energy state by emitting a photon of the same energy
as the absorbed one, the scattering process is called elastic scattering; on the other
4
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
hand, if the excited molecule emits a photon of lower energy, i.e., a photon with a
longer wavelength, the incident photon goes through a process called inelastic
scattering. As seen in [3], the absorption coefficient, that will be introduced later in
this chapter, accounts for both the conversion of radiant energy into heat and the
loss of power at wavelength λ by inelastic scattering to other wavelengths. Mobley
[3] refers to this as true absorption. In view of this, in the remaining of this chapter,
we will treat the scattering as elastic, the loss of energy due inelastic scattering being
included in the absorption coefficient.
To derive a mathematical coefficient for scattering and absorption, we assume
the model shown in Figure 1.
A volume of water ΔV with thickness Δr is illuminated by a light beam with
wavelength λ: Some part of the incident power, PI ðλÞ, will be absorbed within the
volume of water, PA ðλÞ; some part, PS ðλ; θÞ, is scattered out of the beam at an angle θ
and the remaining light power, PT ðλÞ, is transmitted. Then by conservation of energy:
PI ðλÞ ¼ PA ðλÞ þ PS ðλ; θÞ þ PT ðλÞ:
(21)
Using Eq. (22), we can define the ratio between the absorbed and incident
power as the absorptance and the ratio between scattered power and incident power
as scatterance. As already stated, when solving the RTE the IOPs usually employed
are the absorption and the scattering coefficients which are the absorptance and
scatterance per unit distance of the medium. Taking the limit of these terms as the
water thickness Δr becomes infinitesimally small, we have:
P A ðλ Þ
:
∆r!0 PI ðλÞ∆r
aðλÞ ¼ lim
P S ðλÞ
:
∆r!0 PI ðλÞ∆r
bðλÞ ¼ lim
(22)
(23)
The beam attenuation coefficient present in the RTE (Eqs. (16) and (20)) is
written as:
cðλÞ ¼ aðλÞ þ bðλÞ:
(24)
Another useful IOP showing in the RTE is the single scattering albedo, ω0 , that is
defined as the ratio between the scattering and the beam attenuation coefficient:
ω0 ¼
bðλÞ
:
cðλÞ
(25)
In waters where the beam attenuation is due primarily to scattering, ω0 is near
one and when the beam attenuation is due primarily to absorption, ω0 is near zero.
Figure 1.
Geometry of IOPs for a volume of water ΔV.
5
Wireless Mesh Networks - Security, Architectures and Protocols
In [3], it is possible to find the absorption, scattering, and the attenuation
coefficients for several water types and wavelengths. The usual parameters considered for UOWC systems are the ones presented in Table 1.
As mentioned before, the ratio between the scattered power and incident power
is referred as scatterance, that is the fraction of incident power scattered out of the
beam through an angle θ into a solid angle ΔΩ. Defining the angular scatterance per
unit distance and unit solid angle, βðθ; λÞ, as:
PS ðθ; λÞ
:
∆r!0 ∆Ω!0 PI ðλÞ∆r∆Ω
βðθ; λÞ ¼ lim lim
(26)
In [3], we can see that the spectral power scattered into the given solid angle ΔΩ
is just as the spectral radiant intensity, IS ðθ; λÞ, scattered into the angle θ times the
solid angle:
PS ðθ; λÞ ¼ IS ðθ; λÞΔΩ,
(27)
if the incident power falls on an area ΔA, then the corresponding irradiance
Ei ðλÞ ¼ PI ðλÞ=ΔA. Recalling that ΔV is the volume of water that is illuminated by
the incident beam:
IS ðθ; λÞ
:
ΔV!0 Ei ðλÞΔV
βðθ; λÞ ¼ lim
(28)
As said in [5], the form of Eq. (28) suggests the name volume scattering function
(VSF), that is, the scattered intensity per unit incident irradiance per unit volume of
water. Integrating βðθ; λÞ over all directions gives the total scattered power per unit
of incident irradiance and unit volume of water, or as defined in Eq. (23), the
scattering coefficient:
bðλÞ ¼ 2π
ðπ
0
βðθ; λÞsinθdθ:
(29)
Normalizing Eq. (26) with the scattering coefficient, we obtain the scattering
phase function (SPF):
βðλ; θÞ
e
:
β ðθ; λÞ ¼
b ð λÞ
(30)
The scattering phase function can be interpreted as a probability density function (PDF) for scattering from an incident direction ðθ0 ; ϕ0 Þ to a final direction
ðθ; ϕÞ [5]. In general, the SPF must be solved numerically, or using tabulated data
like the ones given by Petzold [6]. The scattering measurements made by Petzold
cðλÞ m
Clear waters
0:15
0.25
Coastal waters
0:4
0.55
Harbor I waters
1:1
0.83
Harbor II waters
2:19
0.83
Table 1.
Measured parameters for different water types.
6
1
Water type
ω0
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
are the most precise to this date, to the best of the author’s knowledge, and widely
cited in the literature. The measurements were carried in the Bahamas, San Pedro in
California, and San Diego harbor also in California, corresponding to clear, coastal,
and harbor waters, respectively.
3. Solving the radiative transfer equation: Monte Carlo simulation
The Monte Carlo name was coined by Nicholas Metropolis [7]. The modern
version of the method was invented in the late 1940s and found use early use in the
studies of neutron transport for the design of nuclear weapons [8, 9]. More recently,
Monte Carlo techniques evolved and are now used to solve problems in physical
sciences, economics, engineering, and pure mathematics.
Monte Carlo methods may be used to solve any problem that have a probabilistic
interpretation, meaning that if the probability of occurrence of each separate event
in a sequence of events is known, one can determine the probability that the entire
sequence of events will occur. By tracing the fate of millions of photons according to
statistical probabilities, a solution for the RTE is built. Each simulated photon path
is randomly distinct from the others, as determined by the probabilities of absorption and scattering in the underwater channel.
The Monte Carlo simulations of photon propagation offer a flexible yet rigorous
approach toward solving numerically the RTE. In the last decades, the exponential
improvement in computer speeds attenuated the problem of the high computational cost of Monte Carlo simulations.
In the following subsections, a detailed description for building a solution for the
RTE, step by step, will be shown.
As seen in [10], there are four main Monte Carlo methods for photon migration
in turbid media, that is, four different ways to build photon’s trajectories. The four
methods are referred as the albedo-weight method (AW), the albedo-rejection
method (AR), the absorption-scattering path length rejection method (ASPR), and
the microscopic Beer-Lambert law method (mBLL).
The AW method considers a source emitting packets of many photons and after
each interaction a fraction of the photons in the packet are absorbed, i.e., lost. In
this type of simulation, the packet is emitted with an initial weight of 1 and at each
interaction the current weight is multiplied by the albedo to account for the loss of
photons by absorption and the packet continues propagating in the scattered direction with a reduced weight. This method has in its favor that there are no wasted
computations because all photon packets eventually reach the detector. This
approach was favored by several research papers [11–16].
In this chapter, we present a simulation using the AR method; the main advantage of this method is that it mimics what happens in nature by tracking individual
photons emitted by some source. Naturally for this case, the computational time of
the absorbed photons are wasted.
It is seen in the investigation conducted by Sassaroli that the total probability of
detecting a photon is equivalent in the four methods and they have statistical
equivalence, with some differences regarding the convergence time. When comparing the AW to the AR method, as one may expect, the AR will require, in most
scenarios, more launched photons for convergence. However, it is also noted that in
the AR method, the photon is detected or lost similarly to what happens in a timeresolved experiment in which photons are detected using the time-correlated single
photon counting technique and the noise on the simulated response reproduces the
Poisson statistics that typically characterize the noise in experiments. Therefore, as
we already stated, the AR method reproduces a more physical simulation of
7
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 2.
Flow chart for the AR Monte Carlo simulation.
propagation than the AW method, where the simulated photons are energy packets
rather than physical photons. In Figure 2, we show a flow chart of the Monte Carlo
algorithm used in this simulation, each step in this chart will be explained with
further detail in the following sections.
3.1 Photon path length
The photon path length, or step size, determines how far the photon will travel
before encountering a water molecule or particle, meaning the distance a photon
travels between an absorption event and a scattering event. The path length is
directly related to the water type and, consequently, to the beam attenuation
coefficient introduced in Section 2.1, and shown for several water types in Table 1.
As stated in [4], from the derivation of the RTE, the radiance in a particular
direction ðθ; ϕÞ decays due to absorption and scattering according to:
8
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
dLðr; θ; ϕÞ
¼
dr
cðrÞLðr; θ ϕÞ,
(31)
where cðrÞ is the beam attenuation coefficient and r is the distance from a given
starting point. Integrating Eq. (31), we have:
Ðr
Lðr; θ; ϕÞ ¼ Lð0; θ; ϕÞe
0
cðr0 Þdðr0 Þ
putting it in terms of the photon path length τ ¼
Ðr
0
Lðτ; θ; ϕÞ ¼ Lð0; θ; ϕÞe τ :
,
(32)
cðr0 Þdr0 :
(33)
The exponential decay of radiance can be explained in terms of the fate of
individual photons if the probability of any photon being absorbed or scattered out
of the beam between photon path lengths τ and τ þ dτ is:
pT ðτÞdτ ¼ e τ dτ:
(34)
Mobley [3] notes that pT ðτÞ satisfies the normalization condition for a probability density function (PDF):
ð x2
x1
px ðxÞdx ¼ 1,
(35)
with x1 ¼ 0 and x2 ¼ ∞. The corresponding cumulative distribution function
(CDF) is:
PT ðτÞ ¼ 1
e τ,
(36)
selecting a random number q from a uniform distribution, 0 ≤ q ≤ 1, and solving
for τ, gives:
q ¼ PT ðτÞ ¼ 1
τ¼
ln ð1
e τ,
qÞ ¼
ln ðqÞ:
(37)
(38)
For homogenous water, where c(r) does not depend on r, where τ ¼ cr, the
photon will travel a geometric distance, r:
r¼
1
ln ðqÞ:
c
(39)
3.2 Absorption events and scattering angles
Recalling that the photon path length, or step size, is the geometric distance a
photon will travel between a scattering event and an absorption event, after determining the distance the photon will travel from origin, one must decide if the next
interaction will be an absorption or a scattering event. For this, we recall the
definition of the single scattering albedo, ωo , that is the ratio between the scattering
and the absorption coefficients. The single scattering albedo is referred sometimes
as the probability of photon survival, because this ratio is the probability that the
photon will undergo scattering instead of absorption. As we have mentioned
before, in this Monte Carlo simulation, we used the AR method and to numerically
implement this, we must select again a random number from a uniform
distribution, 0 ≤ q ≤ 1, to determine which event will take place:
9
Wireless Mesh Networks - Security, Architectures and Protocols
event ¼
absorption, for q ≤ ωo
:
scattering, for q > ωo
(40)
From this description, one may clearly see where the name “Albedo-Rejection”
comes from: if the photon is absorbed, it must be terminated, and its position is
recorded for further analysis; if the photon was scattered, the scattering angles must
be calculated, and its direction of propagation is updated.
3.2.1 The Henyen-Greenstein phase function
Unlike the FSO channel, in an underwater environment, the photons will
encounter a much larger number of particles, and thus multiple scattering events
will play a more important role in the received optical power, when compared to
other optical wireless communication channels.
One of the models proposed to describe the SPF is the Henyen-Greenstein phase
function (HG) [12, 17, 18]. The HG was originally proposed by Henyen and
Greenstein to describe the scattering of light by interstellar dust [19]:
1
1
e
β ð g; θÞ ¼ pHG ¼
4π ð1 þ g 2
g2
3
2gcosθÞ2
,
(41)
where g is referred to as the anisotropy factor and is also the mean cosine of the
scattering angle cos θ, which takes values between 1 and 1. When g ¼ 1, there
is complete backscattering; for g ¼ 0, the scattering is isotropic; and g ¼ 1 gives
complete forward scattering. Using g ¼ 0:924, the HG is a good approximation for
photon propagation in water and a good fit for the Petzold’s measurements.
To find the scattering angle, Eq. (41) must be solved for θ:
(
1
1 þ g2
cos θ ¼
2g
1
1
2 )
g2
,
g þ 2gq
(42)
where q is a uniformly distributed random number between 0 and 1. The complete method for obtaining Eq. (42) may be found elsewhere [4].
Because scattering is a three-dimensional process, besides the polar angle (θÞ
defined above, the azimuthal ðϕÞ angle must also be defined. For an unpolarized
beam, the azimuthal angle has an equal probability to take any value between 0 and
2π [3]. Selecting again a random number, the azimuthal angle ϕ can be found by:
ϕ ¼ 2πq:
(43)
Comparing the HG phase function using Eq. (41) with the Petzold’s avg. Phase
Function extracted from [3], we can see in Figure 3 that the HG is indeed an
adequate approximation for the SPF.
The HG phase function found great use in ocean optics simulations and is
heavily used in Monte Carlo simulations and other relevant channel modeling works
by several researchers [20].
3.2.2 The two-term Henyey-Greenstein phase function
Despite being an approximate fit, it can be seen in Figure 4 that the HG phase
function shows noticeable differences from the Petzold’s measurements for small
angles, θ < 20° , and large angles θ > 150° .
10
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 3.
Comparison of the HG (blue) and the Petzold’s average phase function (red).
Figure 4.
Comparison between the HG and Petzold’s average phase function (left) and large (right) scattering angles.
To overcome this limitation, Haltrin [21] proposed a two-term HeyenGreenstein (TTHG) phase function that matched Petzold’s results better at small
scattering angles and provided the elongation seen at larger angles. The TTHG is
given by:
pTTHG ðθÞ ¼ αpHG θ; g 1 þ ð1
αÞpHG θ; g 2 ,
(44)
where pHG is the HG given by Eq. (41), α is a weighting factor, and g 1 is given a
value near 1, which will make the TTHG increase more strongly at small angles. The
parameter g 2 takes a negative value that will make the TTHG increase, as the angle θ
approaches 180° . The values of α and g2 may be calculated by:
11
Wireless Mesh Networks - Security, Architectures and Protocols
g2 ¼
0:01826332g 21 þ 0:03643748g 1 3 ,
g2 1 þ g2
α¼
:
g1 þ g2 1 þ g2 g1
0:3061446 þ 1:000568g 1
(45)
(46)
Haltrin also derived a relationship that relates the scattering coefficient, b, with
the anisotropy factor, g 1 :
g1 ¼ 1
0:001247
:
b
(47)
After calculating the parameters using Eqs. (44)–(47) and computing the phase
function for clear ocean waters with b ¼ 0:037 m 1 using Eq. (44), we now compare the resulting phase function with the Petzold’s phase function for clear ocean
water. We can see from Figure 5 that the TTHG represents a better approximation
for small angles and presents the characteristic elongation at large angles. Despite
that, in Figure 6, it can also be seen that the TTHG with g 1 ¼ 0:9943 calculated
using Eq. (47), differs greatly from the Petzold’s average phase function data at
intermediate angles, 20–110° . Comparing the results for the TTHG for different
Figure 5.
Comparison of small and large angles for two phase functions, TTHG and HG, with Petzold’s data.
Figure 6.
Comparison of scattering angles for two phase functions with Petzold’s data for different values of the anisotropy
factor g 1 .
12
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
values of g 1 against the Petzold’s average phase function, we can see that
g 1 ¼ 0:9809, as proposed in [17], provides a good fit at small angles, where the
phase function has its largest values, and maintains the elongation at larger angles
characteristic of the Petzold’s phase function.
3.2.3 Haltrin empirical phase function
Despite the attempts from several researchers of developing analytical phase
functions for natural waters, we can see that the ones described in the previous
sections do not provide a good fit over the total range of θ. In [22], Kirk computed
a CDF from Petzold’s data for θ over the range of 2.5–180° in intervals of 5° . Unfortunately, his CDF did not include data for angles smaller than 2:5° , angles where, as
seen before, a lot of scattering happens. In the results presented later in Section 4, we
computed a CDF of Petzold’s data over the whole range of scattering angles and used
linear interpolation to compute θ for each random value selected. The accuracy
obtained by this method comes with the drawback of increased computation time.
In [23], Haltrin developed an empirical fit using strong regression for all
Petzold’s measurements. One of the strongest points of this model is that the Phase
Function may be obtained only by using the scattering coefficient and the single
scattering albedo as inputs:
5
n
n
βðθÞ ¼ exp ψ 1 þ ∑ ð 1Þ kn θð2Þ ,
(48)
n¼1
and the scattering phase function:
5
n
4π
n
ð
Þ
2
exp ψ 1 þ ∑ ð 1Þ kn θ
,
pð θ Þ ¼
b
n¼1
(49)
with:
1
2
ðπ
0
pðθÞ sin ðθÞdθ ¼ 1,
(50)
where b is the scattering coefficient and the remaining parameters are:
pffiffiffi
ψ ¼ 2:598 þ 17:748 b
k1 ¼ 1:188
pffiffiffi
16:722b þ 5:932b b,
0:688ω0 ,
k2 ¼ 0:1ð3:07
k3 ¼ 0:01ð4:58
k4 ¼ 0:001ð3:24
k5 ¼ 0:0001ð0:84
1:90ω0 Þ,
3:02ω0 Þ,
2:25ω0 Þ,
0:61ω0 Þ:
(51)
(52)
(53)
(54)
(55)
(56)
The strong regression given by these equations can be used for an empirical
model for the phase functions. As seen from Figure 7, this model represents a
better fit to the Petzold’s measurements than the phase functions presented before,
especially in the regions where the phase function is the strongest.
In our simulations, we found that when comparing the CDF computed by
directly interpolating the Petzold’s measurements, the CDF computed from the
empirical phase function showed some divergence. Our CDF showed that for
13
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 7.
Comparison of the empirical phase function, the TTHG and HG, against the Petzold’s measurements for small
and large angles.
harbor waters with b ¼ 1:8241 and ω0 ¼ 0:83, P θ ≤ 10° ¼ 0:64; while in the
empirical phase function, we found P θ ≤ 10° ¼ 0:84.
The above presented divergence resulted in a small difference in the final
received power; however, the empirical phase function largely overestimated the
channel bandwidth, especially in harbor waters. This is due to the smaller scattering
angles generated by this empirical phase function, causing the photons to travel in a
straighter line when compared to the case of scattering with the CDF computed
from Petzold’s data.
3.2.4 The Fournier-Forand phase function
As seen the previous sections, several analytical and empirical formulae were
used as empirical fits for ocean waters but none, specially the analytical ones
provide a perfect fit for all scattering angles.
In [24], Fournier and Forand derived an analytic expression for the phase function of a set of particles that have a hyperbolic particle size distribution, with each
particle scattering according to the anomalous diffraction theory. The FournierForand (FF) phase function is given by:
e
βFF ðθÞ ¼
1
v
½vð1 δÞ ð1 δ Þ þ ½δð1
δÞ2 δv
1 δv180
3cos2 θ 1 ,
þ
16πðδ180 1Þδv180
4πð1
v
δ Þ
vð1
θ
δÞsin
2
2
(57)
where
v¼
δ¼
4
3ð n
3
μ
,
(58)
θ
sin
:
2
2
1Þ
(59)
2
2
Being n in Eq. (59) the real index of refraction of the particles, μ in Eq. (58) the
slope parameter of the hyperbolic distribution, and δ180 in Eq. (57) is δ evaluated at
θ ¼ 180° .
14
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
It can be seen in [25] that with n ¼ 1:16 and μ ¼ 3:4586, the FF phase function
gives a good fit to the average Petzold phase function over all scattering angles. This
comparison is shown in Figure 8.
When comparing the CDF computed from the FF phase function and the CDF
from Petzold’s data for harbor waters, we get very close probabilities over all
scattering angles.
3.2.5 Sahu and Shanmugam (SS) phase function with flat back approximation
In [25], a semi-analytical model for a phase function was proposed by Sahu and
Shanmugam (SS). The SS phase function predicts the light scattering properties for
all forward angles, namely 0° ≤ θ ≤ 90° as function of the real index of refraction and
the slope parameter of the hyperbolic distribution. The phase function was divided
in two parts, one from 0.1 to 5° and the other from 5 to 90° . The authors state that
this division was necessary, since the phase function is highly nonlinear with a
change of slope around 5° . This phase function is given by:
8
<
P1 ðln ðθÞÞ2 þ P2 ðln ðθÞÞ þ P3 0:1∘ ≤ θ ≤ 5∘
e
log β ðθÞ ¼
(60)
,
: P ðln ðθÞÞ3 þ P ðln ðθÞÞ2 þ P ðln ðθÞÞ þ P 5∘ ≤ θ ≤ 90∘
1
2
3
4
where
Pm ¼ am e
x
þ bm ðxÞ þ cm ,
(61)
am ¼
dm
þ em sin ð4πyÞ þ f m y þ g m ,
y2
(62)
bm ¼
hm
þ im sin ð4πyÞ þ jm y þ km ,
y2
(63)
cm ¼
lm
þ om sin ð4πyÞ þ pm y þ qm ,
y2
(64)
x¼μ
y¼n
3,
(65)
1,
(66)
Figure 8.
Comparison of the FF and the Petzold average phase function over all scattering angles and for small scattering
angles.
15
Wireless Mesh Networks - Security, Architectures and Protocols
with m ¼ 1, 2, 3 for 0:1° ≤ θ ≤ 5° and m ¼ 1, 2, 3, 4 for 5° ≤ θ ≤ 90° . The values for
the regression coefficients dm , em , f m , g m , hm , im , jm , km , lm , pm , and qm may be found
in [26] and are not shown here for brevity.
Sahu found, in his most recent works [14, 26], that the set of values, n ¼ 1:16
and μ ¼ 3:4319 gives the best fit to Petzold average phase function. Using
Eqs. (60)–(66) and with the values mentioned above for n and μ, we compared the
obtained values with the Petzold average phase function and found that it indeed
provides an excellent fit at all scattering angles, but specially at small angles where
most of the scattering will happen. Results for both parts of the phase function,
0.1–5° and 5–90° , are presented in Figure 9.
3.2.6 Computing phase functions for Monte Carlo simulations
As seen in [27], phase functions satisfy the normalization condition:
2π
ðπ
0
e
β ðθÞsinθdθ ¼ 1,
(67)
and to determine the polar scattering angle:
q ¼ 2π
ðθ
0
e
β ðθÞsinθdθ,
(68)
where q is a uniformly distributed random number between 0 and 1.
However, given the shape of most phase functions, it is more efficient [27] to use
the equation for the phase function to build up a table of θ vs. e
β ðθÞ, compute a CDF
making use of the normalization condition, and then interpolate into the CDF the
value of θ for each random number, q, drawn. As seen in Section 3.2.3, Kirk [22]
computed a CDF for Petzold’s measurements in San Diego Harbor for angles
between 2:5 and 177:5° in intervals of 5° .
For the results presented in Table 2, we computed a CDF directly from Petzold
average phase function extracted from [3] and used the same method described in
[27] to obtain values for all angles between 0 and 90° . We chose to proceed with
numerical integration for this range, and not from 0 to 180° , because as mentioned
in Section 3.2.4, the SS phase function provides values only for forward scattering
angles.
Figure 9.
Comparison between the SS and the Petzold average phase function over the forward scattering angles.
16
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
θ (°)
P(θPETZ)
P(θFF)
P(θSS)
P(θTTHG)
0:126
0.009094
0.012355
0.004059
0.002382
0:2
0.031941
0.041314
0.025992
0.011982
0:251
0.046044
0.057930
0.039860
0.020878
0:316
0.062621
0.076605
0.056320
0.034634
0:398
0.081934
0.097344
0.075422
0.055321
0:501
0.104274
0.120254
0.097285
0.085451
0:631
0.129946
0.145647
0.122155
0.127856
0:794
0.158838
0.173561
0.149991
0.184165
1
0.191217
0.204415
0.181067
0.254854
2:512
0.352290
0.358677
0.335353
0.596932
5:012
0.498215
0.505639
0.472517
0.788187
10
0.659192
0.665537
0.630018
0.897428
30
0.891501
0.882151
0.892093
0.974829
50
0.958125
0.951993
0.960280
0.989702
70
0.986110
0.983167
0.987188
0.996168
80
0.994130
0.992771
0.994681
0.998674
85
0.997280
0.996631
0.997562
0.999169
90
1
1
1
1
Table 2.
Computed CDF for selected phase functions.
In Table 2 and Figures 10 and 11, we compare the CDF obtained for three phase
functions, the SS, FF, and the TTHG with g ¼ 0:9809 against the CDF for Petzold’s
average phase function. In Figure 11, we show the abscissa in logarithmic scale for a
better visualization of the differences between phase functions for small values of θ.
In this table, we can see that both the FF and SS phase functions are a great
approximation for the Petzold’s average phase function, while the TTHG heavily
overestimates the amount of scattering up to 10° and underestimates the amount of
scattering in intermediate angles, i.e., between 30° and 80° . As one may expect, the
effects of those deviations will be larger for waters where scattering is the main
attenuation factor.
3.3 Temporal dispersion and the underwater channel bandwidth
3.3.1 Temporal dispersion
One of the most useful information regarding the underwater channel that may
be provided by the Monte Carlo simulation, which is unavailable in most simple
models like the Beer law, is the channel temporal dispersion. As seen before in this
chapter, seawater is a dispersive medium in which light will suffer the effect of
multiple scattering events. Scattering will cause the photons to arrive at the receiver
at different time intervals causing temporal dispersion and a penalty in the channel
bandwidth. In recent years, several researchers employed the double gamma functions to model the impulse response in underwater channels [15, 28]. Despite the
accuracy of the double gamma function models, we find that the method described
17
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 10.
CDF comparison for selected phase functions.
below consists in a simpler method by means of tracking the individual time of
propagation of each photon.
By tracking the propagation distance of each photon (dprop ), the time of
propagation (TOP), may be estimated by:
top ¼
dprop
,
co
(69)
where co is the speed of light in the water given by:
co ¼
cvacuum
,
nwater
(70)
with the index of refraction of sea water at 20° C and λ ¼ 500 nm being
nwater ¼ 1:34295 [3].
For a better visualization of the impulse response, the time delay may be
normalized by the straight-line distance, or direct distance, from source to the
receiver using:
tdelay ¼ top
18
tdirect :
(71)
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 11.
Log-linear plot of selected phase functions.
Making a histogram of the received photons and distributing the received intensity in time slots defined as the bin widths, one can get a result for the channel
impulse response (CIR) of the underwater channel. The results will be presented in
Section 4.
3.3.2 Underwater channel transfer function
When using laser diodes in UOWC systems, it is desirable to modulate the light
source at very high frequencies. The CIR generated by the Monte Carlo simulation
can be converted to frequency domain by the means of a Fourier Transform algorithm, from which we can obtain the channel transfer function.
Given that the impulse response is approximated by a discrete histogram of
the time of arrival, it can be transformed to a frequency response, as proposed
in [18, 19] by means of a discrete Fourier Transform (DFT), with the usual
definition:
N 1
F½n ¼ DFT½CIR ¼ ∑ tdelay ½ke
j2π
N nk
,
(72)
k¼0
To numerically implement this, the fast Fourier transform (FFT) algorithm can
be used to convert the tdelay data into the frequency domain:
Fð f Þ ¼ abs fft tdelay ðtÞdt ,
(73)
where dt is the time step, or bin width, used when making the histogram of the
time delay data.
19
Wireless Mesh Networks - Security, Architectures and Protocols
As expected, due to Nyquist sampling theorem, the maximum frequency that
can be approximated in the frequency domain will be related to the time step by:
f max ¼
1
:
2tstep
(74)
4. Simulation
In this section, we will provide some simulation results for several emitter
parameters, underwater conditions, and different phase functions.
In the first part, we will simulate the received optical power for clear and coastal
waters for different beam divergence values and compare it with the popular BeerLambert law. It can be seen in [14] that in these water conditions, the delay spread
is so small that they can be considered as a nondispersive medium even for high
data rate communications, and for this reason the channel impulse response (CIR)
will be studied only for conditions where the medium is dispersive.
In the second part, we will study the effect of different phase functions in
Harbor waters with a perfectly collimated beam. We believe that this is an interesting approach since in these waters’, due to the higher scattering coefficient b,
scattering will be a much larger factor and the amount of scattering that occurs at a
given angle is an important factor for the received optical power and the delay
spread. We will also provide numerical results for the CIR and an estimative for the
3-dB underwater channel bandwidth. The simulation parameters for Sections 4.1
and 4.2 are listed in Tables 3 and 4, respectively.
Parameter
Value
Simulated beam
Beam divergence at 1=e2
Gaussian beam
0, 10, 20 and 30 mrad
Beam width
1 mm
Receiver diameter
0.1 m
Receiver FOV
180°
100 106
Photons per simulation
Phase function
Fournier-Forand
Table 3.
Simulation parameters for Section 4.1.
Parameter
Simulated beam
Beam divergence at 1=e2
Gaussian beam
0
Beam width
1 mm
Receiver diameter
0.1 m
Receiver FOV
180o
Photons per simulation
Phase functions
Table 4.
Simulation parameters for Section 4.2.
20
Value
100 106
Petzold average, FF, SS and HG g ¼ 0:924
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
4.1 Clear and coastal waters
4.1.1 Simulation results
In Figure 12, we compare the received power for clear ocean waters with
varying beam divergence parameters of 0, 10, 20, and 30 mrad for an emitted beam
with a FWHM width of 1 mm. For the phase function, we used the FF, since, as seen
in Section 3.2.4, it provides a good fit for Petzold’s phase function. However, as
mentioned in the “Introduction,” when discussing simulation, the impact of different phase functions, especially those which are a good approximation for the
Petzold’s average phase function, is small for the clear and coastal waters due to the
low contribution of scattered photons on the final received power.
Our simulation matches very well the Beer-Lambert law for the case of an
emitted Gaussian beam with no divergence, since one of the main assumptions of
the law is that the emitted beam is perfectly collimated, and as the beam divergence
increases the power loss becomes more accentuated.
We note that for a beam with no divergence, after 10 m, the normalized
received optical power would be 2:2 10 1 for a receiver aperture of 10 cm. Since
the beam spread after a distance z can be calculated by:
ϕ
wðzÞ ¼ 2z tan
,
2
Figure 12.
Received power for clear ocean waters with varying beam divergence parameter.
21
(75)
Wireless Mesh Networks - Security, Architectures and Protocols
we may expect a beam spread of 20 cm after 10 m for a beam divergence of
10 mrad, potentially reducing by more than half the received optical power (due to
the limited aperture of 10 cm). In fact, we can verify from Figure 12 that the
received power for a 10 mrad divergence is 8:1 10 2 , 40% of the power when
the beam was perfectly collimated. The power loss due to the beam spread can be
further confirmed by analyzing the position of the received photons along the
receiver axis of both configurations, the result being shown in Figure 13. From this
figure, we can see that for the case of a beam with a divergence of 10 mrad, the
photons are spread all over the receiver; while, in the case of a collimated beam, the
photons are concentrated along the initial FWHM width of 1 mm.
To further illustrate this, we show, in Figure 14, the geometrical losses
corresponding to beam divergence. Given that for a perfectly collimated beam in
clear ocean waters, power loss is mainly due to absorption and the additional power
loss solely due to the spread of the Gaussian beam.
For the case of coastal waters with c ¼ 0:4 m 1 , we can see from Figure 15 that
the Beer law is still a good approximation for a collimated beam, but the law starts
to underestimate the received optical power due to the increasing contribution of
scattered photons to the final received power.
As was the case for clear waters, the received power drops considerably with an
increase in the beam divergence parameter.
Here, the contribution of scattered photons to the final received power increases:
after 20 m, 25% of the received power came from scattered photons, and after 30 m,
this contribution is elevated to almost 30%. This is further illustrated in Figure 16,
where the scattering histograms for z ¼ 20 and 30 m are depicted, and in Figure 17,
where we compare the received beam for the same distances in clear and coastal
waters. For coastal waters, the effect of scattering is already noticeable at the
distance of z ¼ 10 m and becomes more accentuated at z ¼ 20 m.
Table 5 gives the received normalized power for clear ocean waters for various
distances and beam divergence parameters. In the following table (Table 6), the
same results are given for coastal waters.
4.2 Harbor waters
In this section, we compare the performance of different phase functions for
Harbor I and Harbor II waters. First, we interpolated values for the Petzold average
Figure 13.
Beam spread in clear waters for a receiver located at z = 10 m. 0 rad (left) and 10 mrad (right).
22
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 14.
Beam spread loss in clear waters for beam divergence parameters of 10, 20, and 30 mrad.
Figure 15.
Received power for coastal waters with varying beam divergence parameter.
23
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 16.
Scattering histogram for coastal waters with a perfectly collimated beam.
Figure 17.
Received beams in clear and coastal waters.
phase function and considered, as in Section 3.2.5, the benchmark for all other phase
functions. The amount of unique scattering angles generated by this interpolation
was enough for each photon to have a unique scattering angle.
24
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Distance
ϕ (mrad)
10 m
30 m
50 m
0
2:2 10
1
1:09 10
10
8:1 10
2
8:7 10
20
2:5 10
2
30
1:2 10
2
2
70 m
5:3 10
4
2:6 10
5
4
1:8 10
5
5:8 10
7
2:3 10
4
4:7 10
6
1:3 10
7
1:1 10
4
2:5 10
6
7 10
8
Table 5.
Received power for clear ocean waters.
Distance
ϕ (mrad)
10 m
20 m
30 m
40 m
4:4 10
4
8:6 10
6
1:1 10
7
2
1:15 10
4
1:4 10
6
3:6 10
8
4:1 10
3
4:2 10
5
7:2 10
7
1:9 10
8
2:1 10
3
2:3 10
5
4:5 10
7
1:2 10
8
0
2:2 10
10
1:13 10
20
30
2
Table 6.
Received power for coastal waters.
The SS phase function, seen in Section 3.2.4, predicts scattering for all forward
angles. For the angles between 90 and 180° , we first calculated the backscattering
ratio, as given in [25]:
B p ¼ P1
1
ðn
1Þ2
!P2
,
(76)
3Þ þ c m ,
(77)
where
Pm ¼ am ðμ
3Þ2 þ bm ðμ
with the values for am , bm , and cm being the ones extracted from [25].
Using the values of n ¼ 1:16 and μ ¼ 3:4319, the backscattering ratio
Bp ¼ 0:0185 is obtained, considering the scattering to be uniformly distributed over
all angles above 90° , similarly to what was done in [14].
4.2.1 Simulation results
For the received optical power for Harbor I and Harbor II waters, we can see in
Figures 18 and 19 that the results obtained when using the FF and the SS phase
functions are closely matched to the results obtained using the interpolated Petzold
average PF, meanwhile the one term HG severely underestimates the received optical
power. The main reason for this, as we have seen before, is that the HG is unable to
reproduce the high amount of scattering that occur at small angles, as is the case in
the Petzold average PF and matched by the other analyzed phase functions.
In Figure 20, the received beams are compared for Harbor I and Harbor II
waters for z ¼ 4 and 6 m. This figure clearly shows the higher scattering coefficient
of Harbor II waters and how it impacts the received beam just after 4 m.
25
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 18.
Received power for Harbor I waters for selected phase functions.
Figure 19.
Receiver power for Harbor II waters for selected phase functions.
26
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 20.
Received beams in Harbor I and Harbor II waters.
When looking at the CIR for Harbor I waters with z = 15 m, plotted in Figure 21,
we can see that the results for the Petzold, FF, and the SS phase functions show
similar behavior on the delay profile, while the CIR for the HG phase function
shows a considerable number of delayed photons.
Figure 22 plots the transfer function of the channel, from which the 3-dB
bandwidth can be extracted, for Harbor I waters corresponding to z = 15 and 20 m
using the method described in Section 3.3.2. It is possible to verify from these results
that the behavior of the SS phase function closely matches the one for the Petzold
PF for both distances and the HG underestimates the bandwidth by a large margin
in both scenarios. The FF phase function on the other hand, overestimates the
bandwidth for the case z = 15 m. The 3-dB bandwidth is limited to 300 and
90 MHz, for 15 and 20 m, respectively.
The same analysis is now applied to Harbor II waters. The impulse response is shown
in Figure 23. The interpretation is that the HG phase function exhibits a response where
the power is much more distributed over time when comparing to what is seen in the
other phase functions where most of the power is located at the instant tdirect ¼ 44 ns.
The effect of the delay spread is even more noticeable when we analyze the
transfer function. In Figure 24, this analysis is done for z = 10 and 12 m, where we
see that in these conditions, due to higher number of scattering events that each
photon undergoes, the bandwidth is lower than what we found at Harbor I waters,
even when considering smaller propagation distances. In this case, the phase
27
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 21.
Channel impulse response for Harbor I waters with z = 15 m.
Figure 22.
Transfer function for Harbor I waters. z = 15 and 20 m.
28
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 23.
Channel impulse response for Harbor II waters with z = 10 m.
Figure 24.
Transfer function for Harbor II waters. z = 10 and 12 m.
29
Wireless Mesh Networks - Security, Architectures and Protocols
Received power
τ RMS (ns)
3-dB BW (MHz)
HI z = 15 m
2:6 10
6
2.93
310
HI z = 20 m
1:8 10
7
6.18
80
HII z = 8 m
5:4 10
6
2.06
440
HII z = 10 m
7:7 10
7
3.51
185
HII z = 12 m
1:5 10
7
4.78
110
HII z = 14 m
2:8 10
8
7.87
53
Table 7.
Results for received power, rms delay, and 3-dB BW for Harbor I and Harbor II waters.
functions FF and SS give similar results to the Petzold average PF, the bandwidth
being limited to 150 and 100 MHz, for distances of 10 and 12 m, respectively.
In the following table, Table 7, we summarize the most relevant numerical
results for Harbor I and Harbor II waters as we did for the case of clear and coastal
waters in Section 4.1. For brevity, we only show the results obtained using the
Petzold average phase function; however, as noted in the results presented previously, one may expect similar numerical results when using the FF or the SS phase
function as they predict scattering angles similar to those predicted by Petzold.
4.3 Validation of simulation results
For further validation of the simulation model, we compared the values obtained
in our simulation with the ones obtained by Sahu in [14]. For this, we used the same
parameters used in his simulations, a Gaussian beam profile with FWHM beam width
of 2 mm, a beam divergence parameter of 1.5 mrad, and a receiver aperture of 10 cm.
The comparison for clear and coastal waters between our simulation and [14] is
shown in Figure 25, and we can see that our simulation closely matches the results
obtained by Sahu.
For Harbor waters, Sahu only considered Harbor II waters with c ¼ 2:19 m 1 ,
hence we could not compare the values we obtained for Harbor I waters with
c ¼ 1:1 m 1 . The results are illustrated in Figure 26, where a good agreement is seen.
Furthermore, Sahu defined a delay spread quantity which is the time period over
which the CIR falls to 20 dB below its peak. He provided numerical values for
Figure 25.
Comparison between this simulation model and [14].
30
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Figure 26.
Comparison between this simulation model and [14].
different propagation distances and the comparison against his intercept points for
these distances are depicted in Figure 27. The CIR in that figure is normalized by
the straight-line time of propagation (tdirect ) for a better visualization of the delay
spread. The black line in the figure is the 20 dB intercept point as given in [14],
and we can see that again the values obtained in the simulation are quite close to the
ones obtained by Sahu, which validate the model proposed here.
5. Conclusion
In this chapter, a simple yet powerful tool for modeling the propagation of photons in an underwater channel is presented, which provides more accurate results
than the conventional Beer-Lambert law. The algorithm presented here highlights the
fundamental processes of absorption and scattering, hopefully helping the reader to
better understand the physics of photon propagation in an underwater environment.
In the first part of the simulation, we showed how important the collimation of a
beam is and how the beam spread can cause losses up to 13 dB even after only
10 m. It was also possible to verify that in clear waters, a good signal to noise ratio is
achievable for several tens of meters in both clear and coastal waters.
In the second part of the simulation, the impact of the water turbidity on the
underwater link was assessed, where it was concluded that increasing turbidity limits
not only the received optical power, but also the maximum bandwidth of the channel.
An adequate choice of phase function is important as simpler phase functions, such as
the Henyey-Greenstein, may wrongly predict the numerical results in these conditions.
For the first time results for the received power, channel impulse response and 3-dB
bandwidth are obtained from Monte Carlo RTE simulations for different phase functions. The SS and FF are compared with the Petzold average PF and results are shown to
31
Wireless Mesh Networks - Security, Architectures and Protocols
Figure 27.
CIR for various distances in Harbor II waters, where the delay spread is highlighted.
agree well. Moreover, for validation purposes, our simulation results are further compared with those obtained by Sahu [14] and good agreement is shown to exist.
We believe that the approach presented here provides a general and flexible
technique for numerically solving the radiative transfer equation, accessible to any
reader with basic programming skills. Besides the conditions simulated in this
chapter, the simulation program may be used to simulate innumerous other conditions, like misalignment between receiver and transmitter, smaller aperture sizes,
limited FOVs, and other light sources other than Gaussian beams.
Acknowledgements
This work is financed by the ERDF—European Regional Development Fund
through the Operational Programme for Competitiveness and Internationalisation
—COMPETE 2020 Programme, and by National Funds through the FCT—
Portuguese Foundation for Science and Technology, I.P., within project
POCI-01-0145-FEDER-031971.
32
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
Author details
Rafael M.G. Kraemer1,2, Luís M. Pessoa2 and Henrique M. Salgado1,2*
1 Faculty of Engineering, University of Porto, Portugal
2 INESC TEC—INESC Technology and Science (Formerly INESC Porto), Porto,
Portugal
*Address all correspondence to: henrique.salgado@inesctec.pt
© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms
of the Creative Commons Attribution License (http://creativecommons.org/licenses/
by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,
provided the original work is properly cited.
33
Wireless Mesh Networks - Security, Architectures and Protocols
References
[1] Inherent Optical Properties. Ocean
[10] Sassaroli A, Martelli F. Equivalence
Optics Web Book [Online]. Available
from: http://www.oceanopticsbook.inf
o/view/overview_of_optical_oceanogra
phy/inherent_optical_properties
[Accessed: 28 September 2018]
of four Monte Carlo methods for photon
migration in turbid media. Journal of the
Optical Society of America A. 2012;
29(10):2110
[2] Zeng Z, Fu S, Zhang H, Dong Y,
Cheng J. A survey of underwater optical
wireless communications. IEEE
Communication Surveys and Tutorials.
2017;19(1):204-238
[3] Mobley CD. Light and Water—
Radiative Transfer in Natural Waters.
San Diego: Academic Press; 1994
[4] Leathers RA, Downes TV, Davis CO,
Mobley CD. Monte Carlo Radiative
Transfer Simulations for Ocean Optics: A
Practical Guide. Fort Belvoir, VA: Defense
Technical Information Center; 2004
[5] The Volume Scattering Function.
[11] Cox WC, Jr, Simulation, Modeling,
and Design of Underwater Optical
Communication Systems. PhD Thesis,
North Carolina State University; 2012
[12] Gabriel C, Khalighi M-A, Bourennane
S, Léon P, Rigaud V. Monte-Carlo-based
channel characterization for underwater
optical communication systems. Journal
of Optical Communications and
Networking. 2013;5(1):1
[13] Li J. Monte Carlo study on pulse
response of underwater optical channel.
Optical Engineering. 2012;51(6):066001
[14] Sahu SK, Shanmugam P. A
Ocean Optics Web Book [Online].
Available from: http://www.oceanoptic
sbook.info/view/overview_of_optical_
oceanography/the_volume_scattering_
function [Accessed: 30 September 2018]
theoretical study on the impact of
particle scattering on the channel
characteristics of underwater optical
communication system. Optics
Communications. 2018;408:3-14
[6] Petzold TJ. Volume Scattering
[15] Li Y, Leeson MS, Li X. Impulse
Functions for Selected Ocean Waters.
Fort Belvoir, VA: Defense Technical
Information Center; 1972
response modeling for underwater
optical wireless channels. Applied
Optics. 2018;57(17):4815
[7] Metropolis N, Ulam S. The Monte
[16] Vadakke-Chanat S, Shanmugam P,
Carlo method. Journal of the American
Statistical Association. 1949;44(247):
335-341
Sundarabalan B. Monte Carlo simulations
of the backscattering measurements for
associated uncertainty. Optics Express.
2018;26(16):21258
[8] Computing and the Manhattan
Project. Atomic Heritage Foundation
[Online]. Available from: https://www.
atomicheritage.org/history/computingand-manhattan-project [Accessed: 22
October 2018]
[9] Monte Carlo and the Bomb [Online].
Available from: http://jmoses.co/2014/
05/27/monte-carlo-and-the-bomb.html
[Accessed: 01 November 2018]
34
[17] Mobley CD, Sundman LK, Boss E.
Phase function effects on oceanic light
fields. Applied Optics. 2002;41(6):1035
[18] Wang L, Jacques SL. MCML—
Monte Carlo modeling of light transport
in multi-layered tissues in standard
C. 183. Computer Methods and
Programs in Biomedicine. Boston, MA:
Springer; 1995;47(2):131-146
Monte Carlo Radiative Transfer Modeling of Underwater Channel
DOI: http://dx.doi.org/10.5772/intechopen.85961
[19] Henyey L, Greenstein J. Diffuse
radiation in the galaxy. Astophysical
Journal. 1941;93:70-83
[20] Miramirkhani F, Uysal M. Visible
light communication channel modeling
for underwater environments with
blocking and shadowing. IEEE Access.
2018;6:1082-1090
[21] Haltrin VI. One-parameter twoterm Henyey-Greenstein phase function
for light scattering in seawater. Applied
Optics. 2002;41(6):1022
[22] Kirk JT. Monte Carlo Procedure for
Simulating the Penetration of Light into
Natural Waters. Australia:
Commonwealth Scientific and Industrial
Research Organization; 1981
[23] Haltrin VI. Theoretical and
empirical phase functions for Monte
Carlo calculations of light scattering in
seawater. In: Presented at the Fourth
International Conference on Remote
Sensing for Marine and Coastal
Environments. 1997;17:19
[24] Fournier GR, Luc Forand J. Analytic
phase function for ocean water. Proc.
SPIE Ocean Optics XII. 1994;2258:2258
[25] Sahu SK, Shanmugam P. Semi-
analytical modeling and
parameterization of particulates-inwater phase function for forward
angles. Optics Express. 2015;23(17):
22291
[26] Sahu SK, Shanmugam P. Scattering
phase function for particulates-in-water:
Modeling and validation. In: Presented
at the SPIE Asia-Pacific Remote Sensing;
New Delhi, India; 2016. p. 98821H
[27] Mobley CD et al. Comparison of
numerical models for computing
underwater light fields. Applied Optics.
1993;32(36):7484
[28] Tang S, Dong Y, Zhang X. Impulse
response modeling for underwater
35
wireless optical communication links.
IEEE Transactions on Communications.
2014;62(1):226-234