Monte Carlo Radiative Transfer Modeling of Underwater Channel

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

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=1gives 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 θ:

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:

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°.

To overcome this limitation, Haltrin [21] proposed a two-term Heyen-Greenstein (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:

where pHGis the HG given by Eq. (41), αis a weighting factor, and g1is given a value near 1, which will make the TTHG increase more strongly at small angles. The parameter g2takes a negative value that will make the TTHG increase, as the angle θapproaches 180°. The values of αand g2may be calculated by:

Haltrin also derived a relationship that relates the scattering coefficient, b, with the anisotropy factor, g1:

After calculating the parameters using Eqs. (44)–(47) and computing the phase function for clear ocean waters with b=0.037m−1using 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 g1=0.9943calculated 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 values of g1against the Petzold’s average phase function, we can see that g1=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:

and the scattering phase function:

with:

where b is the scattering coefficient and the remaining parameters are:

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 harbor waters with b=1.8241and ω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 Fournier-Forand (FF) phase function is given by:

where

Being n in Eq. (59) the real index of refraction of the particles, μin Eq. (58) the slope parameter of the hyperbolic distribution, and δ180in Eq. (57) is δevaluated at θ=180°.

It can be seen in [25] that with n=1.16and μ=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:

where

with m=1,2,3for 0.1°≤θ≤5°and m=1,2,3,4for 5°≤θ≤90°. The values for the regression coefficients dm,em,fm,gm,hm,im,jm,km,lm,pm, and qmmay 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.16and μ=3.4319gives the best fit to Petzold average phase function. Using Eqs. (60)–(66) and with the values mentioned above for nand μ, 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:

and to determine the polar scattering angle:

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. β∼θ, 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.5and 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.

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.9809against 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.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 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:

where cois the speed of light in the water given by:

with the index of refraction of sea water at 20°Cand λ=500nmbeing 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:

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:

To numerically implement this, the fast Fourier transform (FFT) algorithm can be used to convert the tdelaydata into the frequency domain:

where dtis the time step, or bin width, used when making the histogram of the time delay data.

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:

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−1for a receiver aperture of 10 cm. Since the beam spread after a distance zcan be calculated by:

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.4m−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=20and 30mare 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=10mand becomes more accentuated at z=20m.

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.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=4and 6m. This figure clearly shows the higher scattering coefficient of Harbor II waters and how it impacts the received beam just after 4 m.

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=44ns.

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 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.