Analytical solution of beam spread function for ocean light radiative transfer

Zao Xu; Dick K. P. Yue

doi:10.1364/OE.23.017966

1. Introduction

Knowledge of the response to highly collimated light beam in an anisotropic scattering medium such as turbid ocean, atmosphere and human tissue is widely used in underwater optical imaging and communication, ocean light detection and ranging (LIDAR), graphic rendering, and optical diffuse tomography. In these contexts, the spatial distribution of irradiance (usually a lateral spreading with increasing depth) in response to a narrow uni-directional light beam is defined as the beam spread function (BSF) (also called the point spread function or PSF). In a macroscopical perspective, the light propagation in a scattering medium is fully described by the (3D) radiative transfer equation (RTE). Thus, the BSF generally requires the solution of the RTE, which is a significant task. Our interest is in applications in ocean optics where, due to the highly peaked anisotropic distribution of the scattering phase function for oceanic particles, the BSF for ocean light radiative transfer is particularly difficult to solve both analytically and numerically [1].

The BSF in ocean optics has traditionally been obtained through experiments [2, 3] and empirical models [4, 5]. Direct measurements of spatially distributed irradiance in the field is challenging and expensive. Empirical models, while practically efficient in some applications, generally lack rigorous derivation and are often unable to reveal the underlying physics and apparent optical properties. Direct numerical or Monte Carlo solutions of the RTE are also available [6, 7], but the computational effort is generally large. Computational solution of the BSF directly has also been attempted. Sanchez and McCormick [8] utilize the discrete ordinate method for directional variables and finite difference in spatial variables to obtain the numerical BSF; however, the computational cost is still comparable to solving the RTE.

The analytical solution of the RTE Green’s function is fundamental and has the advantage of direct physical insights and the potential of highly computational efficiency. For the BSF, most analytical work is based on the diffusion approximation (DA) [1, 9] of the RTE. However, a diffusion approximation is mainly accurate for strong (single scattering albedo ω₀ ≳ 0.99) and/or isotropic scattering. For highly forward peaked scattering, asymptotic models [10, 11] have also been developed. These models assume that the scattering polar angle θ of radiance is highly peaked at small values so that μ = cosθ ≈ 1. Expanding the RTE in an asymptotic series for small θ obtains a sequence of advection-diffusion equations for different orders. The asymptotic model is valid when the leading term dominates which means very small optical depth and lateral spreading. Another important approach is the small-angle approximation (SAA) [12, 13]. The basic idea in SAA is again the assumption of θ ≪ 1 (μ ≈ 1), but in addition, the angle between any two ray directions is approximated by the difference of their projections in the lateral plane. With SAA, the integral operator in RTE is reduced to a convolution operator and the RTE can be solved directly in frequency domain. For a uni-directional ocean source, McLean et al [14] obtained an expression for the BSF involving two double Fourier integrals, whose quadrature still requires significant numerical effort. Further simplifications of SAA are possible (e.g.,[15]), but generally leads to appreciable loss in accuracy. In general, Swanson et al [16, 17] show that SAA is valid only for small optical depths. More recently, Markel [18] and Liemert and Kienle [19] developed a rigorous method for the 3D RTE Green function using the rotated-frame method. The resulting solution form is quite involved (containing multiple summations of Fourier integrals) and is computationally less practical.

We present a highly-robust and efficient analytical solution for the 3D BSF using a new azimuthal single scattering approximation to reduce the original 3D RTE to a standard 2D RTE, and then combine the convolution operator idea of SAA and the rotated-frame method to effectively solve for the (exact) 2D RTE. The azimuthal single scattering approximation is based on the assumption that the azimuthal dependence of the radiance is highly-peaked in the forward direction, which typically obtains in ocean optics. The 2D BSF solution [see Eq. (37)] contains only two summations, while the 3D BSF is related to the 2D BSF via a simple relationship [see Eq. (44)]. Comparing our analytical BSF results to Monte Carlo simulation of the original RTE [20], we confirm that the 2D BSF is exact with excellent comparisons between the two for all parameter ranges. For the 3D BSF, the comparisons are very satisfactory for single scattering albedo values typical in the ocean. As expected, deviations occur only for large optical depths (> O(5)) at lateral positions close to the beam, and for very large single scattering albedo (up to 0.9).

Due to the relative simple form(s) of the 2D and 3D analytical solution(s), the present BSF is very amenable to numerical evaluation, and more efficient than direct Monte Carlo simulation especially for single scattering albedo that is not small. This approach provides great potentials for practical applications of both forward and inverse ocean optics problems, particularly in areas of optical oceanography includes predicting the LIDAR returns and the light scattering in the ocean turbulence. The method is espeecially suitable to deal with light scattering by turbulence of a turbid ocean, which is usually obtained by solving 3D RTE with temporally/spatially varying inherent optical properties. Considering the fact that the BSF is equivalent to the Green’s function of RTE, one can immediately obtain the scattered light field by the turbid ocean turbulence using the first-order approximation and the present analytic approach.

2. Methods

2.1. Variables reduction of three-dimensional radiative transfer equation

In general, the beam spread function (BSF) is the spatial distributions of the light intensity in response to a narrow incident beam. In this work, we define the BSF as the scalar irradiance which is the function of the lateral spatial variables (x, y) and depth variable z for an uni-directional point source inside a homogenous three-dimensional horizontally-unbounded ocean column of large depth (Fig. 1). Here all position variables are normalized by the beam attenuation coefficient to represent the optical distances. We use E⁽³⁾(x, y, z) to represent the scalar irradiance, where for later reference (and in contrast to the two-dimensional case), the superscript denotes the irradiance is in three dimensions (3D). We assume the uni-directional light source is placed at the origin illuminates in positive z-direction. The scalar irradiance is ax-isymmetric with respect to the z-axis, and it is convenient to consider the problem in cylindrical coordinates E⁽³⁾(x, y, z) = E⁽³⁾(r, z).

Fig. 1 Axisymmetric geometry of the problem: uni-directional point source illuminating in positive z-axis under the azimuthal single scattering approximation.

Download Full Size | PDF

The scalar irradiance is determined by the radiance L(r, z, θ, ϕ − ϕ_r); L is the function of both the spatial variables (r, ϕ_r, z), and directional variables (θ, ϕ); θ is the polar angle and ϕ is the azimuthal angle. Introducing a new azimuthal variable ψ = ϕ − ϕ_r, the scalar irradiance can be defined as:

E^{(3)} (r, z) = \int_{0}^{2 π} \int_{0}^{π} L (r, z, θ, ψ) sin θ d θ d ψ .

The radiance L is a function of five independent variables and is the solution of the 3D radiative transfer equation (RTE):

\begin{array}{l} cos θ \frac{\partial L}{\partial z} + sin θ (cos ψ \frac{\partial L}{\partial r} - \frac{sin ψ}{r} \frac{\partial L}{\partial ψ}) + L & = ω_{0} \int_{0}^{2 π} \int_{0}^{2 π} L (θ^{'}, ψ^{'}) Φ (θ^{'}, ψ^{'}; θ, ψ) sin θ^{'} d θ^{'} d ψ^{'} \\ + \frac{1}{π r} δ (r) δ (z) δ (θ) δ (ψ), \end{array}

where δ(·) is the Dirac delta function; Φ(θ′, ψ′; θ, ψ) is the phase function with arbitrary shape. So far, no assumption on the phase function and the solution has been made yet. The closed-form solution of Eq. (2) is generally very difficult to obtain [18, 19].

Now assumption can be made in the context of the typical ocean environment (e.g., case 1 water [6]) where light scattering tends to be highly forward peaked. The typical phase functions, such as the Petzold [21], and the Henyey-Greenstein (H-G) [22] with the asymmetric parameter g > 0.9, are highly peaked in the forward direction. Under the geometry of an uni-directional beam that emanates from a point source into the ocean, the scattering of light tends to spread more energy in the polar direction than in the azimuthal direction; the lateral spreading of a (narrow) incident beam is small relative to the forward direction. Therefore, a simplifying assumption is possible based on the highly peaked scattering property and is called ’azimuthal single scattering approximation’. In this approximation, we assume that in the azimuthal direction the scattering is relatively less important than in the polar direction, and therefore only single scattering occurs in the azimuthal direction. This assumption can be expressed mathematically by deliberately separating out azimuthal angle ψ in the radiance L. The dependence of L on azimuthal angle ψ is thus nontrivial only for |ψ| ≪ 1, especially for locations not close to the incident beam. According to the characteristics of single scattering [23], the radiance L now can be written as:

L (r, z, θ, ψ) = \hat{L} (r, z, θ) δ (ψ) .

Figure 1 demonstrates the geometry of the azimuthal single scattering approximation. It is seen that under such approximation, scattering of light is assumed to take place within the plane OABO_r which is identical to the plane ψ = 0 and π.

For forward-peaked scattering, the radiance is small except near θ = 0, and we write sin(θ′) ≈ sin(θ′ − θ) and cos(θ′) ≈ cos(θ′ − θ). We replace L in Eq. (2) with Eq. (3) and perform the integration over ψ to obtain:

cos θ \frac{\partial \hat{L}}{\partial z} + sin θ (\frac{\partial \hat{L}}{\partial r} + \frac{\hat{L}}{r}) + L = 2 π ω_{0} \int_{0}^{π} \hat{L} (θ^{'}) p (θ^{'} - θ) sin θ^{'} d θ^{'} + \frac{1}{π r} δ (r) δ (z) δ (θ) .

Note that the 3D phase function Φ in Eq. (2) is now replaced by a 2D phase function p(θ′ − θ). The function p can be expressed explicitly by expanding Φ in spherical harmonics, and upon integration over ψ, to yield:

p (θ^{'} - θ) = 2 π | sin (θ^{'} - θ) | \sum_{n = 0}^{\infty} g_{n} P_{n} [cos (θ^{'} - θ)],

where P_n(x) is the Legendre polynomial and the coefficient g_n is given by:

g_{n} = \frac{2 n + 1}{2} \int_{- 1}^{1} Φ (x) P_{n} (x) d x .

By defining the azimuthal-averaged radiance I(r, z, θ):

I (r, z, θ) = π r \hat{L} (r, z, θ),

the Eq. (4) can be written in two-dimensional form:

cos θ \frac{\partial I (r, z, θ)}{\partial z} + sin θ \frac{\partial I (r, z, θ)}{\partial r} + I (r, z, θ) = ω_{0} \int_{0}^{π} p (θ - θ^{'}) I (r, z, θ^{'}) d θ^{'} + δ (r) δ (z) δ (θ) .

We change the variable r to x and extending the original range of r from (0, ∞) to (−∞, ∞); at the same time, we extend range of the polar angle θ from (0, π) and [−π, π]. We are able to obtain the standard two-dimensional RTE with uni-directional point source:

cos θ \frac{\partial I (x, z, θ)}{\partial z} + sin θ \frac{\partial I (x, z, θ)}{\partial x} + I (x, z, θ) = ω_{0} \int_{0}^{2 π} p (θ - θ^{'}) I (x, z, θ^{'}) d θ^{'} + δ (x) δ (z) δ (θ) .

We define scalar irradiance for the two-dimensional RTE:

E^{(2)} (x, z) = \int_{0}^{2 π} I (x, z, θ) d θ .

The relation between the original 3D scalar irradiance E⁽³⁾(r, z) and the 2D E⁽²⁾(x, z) is simply:

E^{(3)} (r, z) = \frac{1}{π r} E^{(2)} (| x |, z) .

Note that Eq. (11) is new and is an immediate result of the single scattering approximation in azimuthal direction we introduce. Thus, for conditions such as typical ocean where the scattering is highly focussed in the forward direction, the 3D BSF can be obtained directly and explicitly from the 2D BSF.

2.2. Analytical solution of two-dimensional scalar irradiance

We solve the standard 2D RTE (9) by using Fourier spectral method [24] in combination with the rotated frame method [18, 19]. We expand the 2D radiance I(x, z, θ) in terms of Fourier integrals:

I (x, z, θ) = \frac{1}{{(2 π)}^{2}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \hat{I} (k_{x}, k_{z}, θ) e^{- i (x k_{x} + z k_{z})} d k_{x} d k_{z},

where the coefficients are:

\hat{I} (k_{x}, k_{z}, θ) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} I (x, z, θ) e^{i (x k_{x} + z k_{z})} d x d z,

k_x=k cosθ_k, k_z=k sinθ_k. Equation (9) becomes:

(i k cos (θ - θ_{k}) +) \hat{I} (k, θ_{k}, θ) = ω_{0} \int_{0}^{2 π} p (θ - θ^{'}) \hat{I} (k, θ_{k}, θ^{'}) d θ^{'} + S (k, θ_{k}, θ) .

We define a new variable α = θ − θ_k and write the radiance as

\hat{I} (k, θ_{k}, θ) = \hat{I} (k, θ_{k}, α) = \sum_{m = - \infty}^{\infty} I_{m} (k, θ_{k}) e^{i m α}, and

I_{m} (k, θ_{k}) = \frac{1}{2 π} \int_{0}^{2 π} \hat{I} (k, θ_{k}, α) e^{- i m α} d α .

Similarly, the phase function p can be written as

p (θ - θ^{'}) = \frac{1}{2 π} \sum_{m = - \infty}^{\infty} p_{m} e^{i m (θ - θ^{'})}, and

p_{m} = \int_{0}^{2 π} p (θ - θ^{'}) e^{- i m (θ - θ^{'})} d θ^{'} .

The source function in the spectral domain becomes

S_{m} (k, θ_{k}) = \frac{1}{2 π} \int_{0}^{2 π} δ (θ) e^{- i m α} d α = \frac{1}{2 π} e^{i m θ_{k}} .

Substituting I and p into Eq. (14), we obtain

i \frac{k}{2} (I_{m - 1} + I_{m + 1}) + (1 - ω_{0} p_{m}) I_{m} = S_{m} .

Truncation at some finite number of modes −N < m < N, we obtain a algebraic linear system for the unknown modal amplitudes:

i \frac{k}{2} A \bar{X} + B \bar{X} = \bar{S},

\bar{X} = {[I_{- N} \dots I_{0} \dots I_{N}]}^{T} = {[x_{1} \dots x_{N} \dots x_{2 N + 1}]}^{T},

\bar{S} = {[S_{- N} \dots S_{0} \dots S_{N}]}^{T} = {[s_{1} \dots s_{N} \dots s_{2 N + 1}]}^{T},

s_{ℓ} = S_{ℓ - N - 1} = \frac{1}{2 π} e^{i (ℓ - N - 1) θ_{k}} .

Accordingly, A and B are symmetric tri-diagonal and diagonal matrices, respectively, with elements:

A_{n, m} = δ_{n, m + 1} + δ_{n + 1, m}, and B_{n, m} = β_{n} δ_{n, m},

where δ_n,m is the Kronecker delta function and β_n=1 − ω₀p_n.

To solve X̄ analytically, we symmetrize the matrix A by defining

C = B^{- 1 / 2}, and D = CAC .

We obtain, upon formal inversion:

\bar{X} = C {(i \frac{k}{2} D + I)}^{- 1} C \bar{S},

where I is the identity matrix and D is a symmetric matrix with elements:

D_{n, m} = \frac{1}{\sqrt{β_{n} β_{m}}} (δ_{n, m + 1} + δ_{n + 1, m}) .

D is Hermitian and can be decomposed as

D = V Σ V^{- 1},

where the eigenmatrix V has real elements V_nm, and Σ is diagonal with (real) eigenvalue λ_m.

Solving Eq. (27), we obtain the solution X̄:

x_{n} = \sum_{ℓ = 1}^{2 N + 1} \sum_{m = 1}^{2 N + 1} \frac{V_{n m} V_{ℓ m}}{\sqrt{β_{n} β_{ℓ}}} \frac{1}{i k λ_{m} / 2 + 1} s_{ℓ} .

Accordingly,

\hat{I} (k, θ_{k}, θ) = \sum_{ℓ = 1}^{2 N + 1} \sum_{m = 1}^{2 N + 1} \frac{V_{N + 1, m} V_{ℓ m}}{\sqrt{β_{N + 1} β ℓ}} \frac{1}{i k λ_{m} / 2 + 1} e^{i (ℓ - N - 1) θ_{k}} .

Applying Eq. (10), we have

E^{(2)} (k, θ_{k}) = \int_{0}^{2 π} \hat{I} (k, θ_{k}, θ) d θ .

The 2D scalar irradiance in physical space E⁽²⁾(x, z) can be obtained after inverse Fourier transform in cylindrical coordinates (ρ, θ_ρ): z=ρ cosθ_ρ, x=ρ sinθ_ρ. We obtain:

E^{(2)} (x, z) = E^{(2)} (ρ, θ_{ρ}) = \frac{1}{{(2 π)}^{2}} \int_{0}^{2 π} \int_{0}^{\infty} k E^{(2)} (k, θ_{k}) e^{i k ρ cos (θ_{k} - θ_{r})} d k d θ_{k} .

To perform the θ_k integral in Eq. (33), we use the identity

e^{i k ρ cos (θ_{k} - θ_{ρ})} = \sum_{n = - \infty}^{\infty} i^{n} J_{n} (k ρ) e^{i n (θ_{k} - θ_{ρ})},

where J_n(x) is the bessel function of the first kind, to obatin

E^{(2)} (ρ, θ_{ρ}) = \frac{1}{2 π} \sum_{ℓ = 1}^{2 N + 1} \sum_{m = 1}^{2 N + 1} i^{ℓ - N - 1} \frac{V_{n m} V_{ℓ m}}{\sqrt{β_{N + 1} β ℓ}} e^{i (ℓ - N - 1) θ_{ρ}} \int_{0}^{\infty} \frac{k}{i k / 2 λ_{m} + 1} J_{ℓ - N - 1} (k ρ) d k .

To simplify this result, we invoke the properties of the eigenvalues and eigenmatrix of D:

\begin{array}{l} V_{2 N + 2 - n, m} = {(- 1)}^{m - 1} V_{n m}, V_{n, 2 N + 2 - m} = {(- 1)}^{n - 1} V_{n m}, V_{2 N + 2 - n, 2 N + 2 - m} = {(- 1)}^{m + n} V_{n m}, \\ β_{m} = β_{2 N + 2 - m}, λ_{m} = - λ_{2 N + 2 - m}, λ_{N + 1} = 0 . \end{array}

Considering the real part of the results, we obtain the exact 2D scalar irradiance expressed in the form:

\begin{array}{l} E^{(2)} (ρ, θ_{ρ}) & = \frac{4}{π} \sum_{m = 1}^{N} \frac{V_{N + 1, m}^{2}}{β_{N + 1} λ_{m}^{2}} K_{0} (\frac{2 ρ}{λ_{m}}) \\ + \frac{1}{2 π} \sum_{n = 1}^{N} \sum_{m = 1}^{N} \frac{V_{N + 1, m} V_{n + N + 1, m}}{\sqrt{β_{N + 1} β_{n + N + 1}}} {[1 + {(- 1)}^{n}] [1 + {(- 1)}^{n + m - 1}] cos [n (θ_{ρ} + π / 2)] H_{n}^{R} (ρ, λ_{2 m + 1}) \\ + \frac{2}{λ_{m}} [1 + {(- 1)}^{n + 1}] [1 + {(- 1)}^{n + m}] sin [n (θ_{ρ} + π / 2)] [\frac{1}{ρ} - H_{n}^{I} (ρ, λ_{2 m + 1})]}, \end{array}

where K_n(x) is the modified Bessel function of the second kind and

\begin{array}{l} H_{n}^{R} (ρ, λ) & = \int_{0}^{\infty} \frac{k J_{n} (k ρ)}{1 + k^{2} λ^{2} / 4} d k \\ = \frac{4}{λ^{2} n} {{}_{1}F}^{2} (1; 1 + \frac{n}{2}, 1 - \frac{n}{2}; \frac{ρ^{2}}{λ^{2}}) - \frac{2 π}{λ^{2}} I_{n} (\frac{2 ρ}{λ}) csc (\frac{n π}{2}), \end{array}

\begin{array}{l} H_{n}^{I} (ρ, λ) & = \int_{0}^{\infty} \frac{J_{n} (k ρ)}{1 + k^{2} λ^{2} / 4} d k \\ = \frac{4 ρ}{λ^{2} (n^{2} - 1)} {{}_{1}F}^{2} (1; \frac{3}{2} - \frac{n}{2}, \frac{3}{2} + \frac{n}{2}; \frac{ρ^{2}}{λ^{2}}) - \frac{π}{λ} I_{n} (\frac{2 ρ}{λ}) sec (\frac{n π}{2}) . \end{array}

Here ₁F²(x₁; x₂, x₃; x₄) is the hypergeometric function and I_n(x) is the modified Bessel function of the first kind.

Following the same approach, we can derive the series form of the solution for the ballistic source. The only differences between the complete solution and ballistic solution are the eigenvalues and eigenvectors. The eigenvalues $λ_{m}^{b}$ and eigenvectors $V_{n m}^{b}$ for ballistic solution can be obtained analytically:

λ_{m}^{b} = 2 cos (\frac{m π}{2 N + 1}),

and

V_{n m}^{b} = \frac{1}{\sqrt{N + 1}} sin (\frac{n m π}{2 N + 2}) .

Therefore, the final solution of the scattered light fields for 2D BSF from an unidirectional beam source is

E_{s}^{(2)} (ρ, θ_{ρ}) = E^{(2)} (ρ, θ_{ρ}) - E_{b}^{(2)} (ρ, θ_{ρ}) .

Change the variable of (ρ, θ_ρ) to (x, z), and we obtain the expression of the beam spread function obtained from a standard two-dimensional radiative transfer.

E_{s}^{(2)} (x, z) = E^{(2)} (\sqrt{x^{2} + z^{2}}, arctan \frac{x}{z}) - E_{b}^{(2)} (\sqrt{x^{2} + z^{2}}, arctan \frac{x}{z}) .

It is worth mentioning that $H_{n}^{R} (ρ, λ)$ and $H_{n}^{I} (ρ, λ)$ are solutions after inverse Fourier transforms for angular variables. Though they have explicit forms, they must be treated carefully. The two functions are both substractions of two terms. The two functions are non-singular but both terms in the function are singular for certain n, e.g., when n is even both ₁F² and csc(nπ/2) in $H_{n}^{R}$ are singular; whilst n is odd ₁F² and sec(nπ/2) in $H_{n}^{I}$ are singular. The accurate approach to evaluate these two functions is to expand the singular terms at n; the singular terms are crossed out and the non-singular terms in the expansion can be calculated easily. A less accurate but simpler approach is to take the limit of n numerically using the fact that both integrals are continuous with n.

To summarize the method to solve the two-dimensional standard RTE analytically, the radiance is expanded to Fourier series for both spatial and angular variables. The eigenvalues and eigenvectors must be calculated and stored in advance. The scalar irradiance as the function of depth and lateral position is obtained by the summation of terms containing the eigenvalues and elements of eigenvectors.

The solution of 3D beam spread function of light radiative transfer in an infinite ocean is obtained based on the solution of 2D RTE through Eq. (11) and (43); therefore, the 3D BSF for the scattered light fields can be expressed as:

E_{s}^{(3)} (r, z) = \frac{1}{π r} [E^{(2)} (\sqrt{r^{2} + z^{2}}, arctan \frac{r}{z}) - E_{b}^{(2)} (\sqrt{r^{2} + z^{2}}, arctan \frac{r}{z})] .

3. Results

3.1. Convergence evaluation

The validation of the present theory is made by comparisons with Monte Carlo radiative transfer simulation in two and three dimensions for arbitrary scattering phase function [25, 26]. The codes are parallelized for high-performance computing platforms. All the results of BSF obtained by the Monte Carlo simulation for comparisons are performed by launching O(10¹¹) photons.

Before validation of the (exact) 2D BSF solution [Eq. (37)], we first evaluate the convergence of Eq. (37) with increasing Fourier modes N. For simplicity, we utilize the 2D Henyey-Greenstein (H-G) phase function [22] in the evaluations:

p (Φ) = p (θ - θ^{'}) = \frac{1}{2 π} \frac{1 - g^{2}}{1 + g^{2} - 2 g cos (θ - θ^{'})}, g = 0.924,

where Φ = θ − θ′ is scattering angle and g is the asymmetric parameter.

The dependence of the solution on N is investigated with N varying from 40 to 800. The relative error (RE) of the scattered scalar irradiance is defined as

RE = \frac{| E_{s, N}^{(2)} - E_{s, \infty}^{(2)} |}{| E_{s, \infty}^{(2)} |},

where

E_{s, \infty}^{(2)}

is the best estimate of the true value using Richardson extrapolation on the calculated data. Figure 2 shows RE for three different single scattering albedos ω₀ at different horizontal positions x (at a given depth). Exponential convergence RE∼ e^−κN is clearly obtained. For given x, κ is independent of ω₀ (although RE is somewhat greater for smaller ω₀). With increasing x (at a given depth), κ is also somewhat greater.

Fig. 2 Evaluation of the convergence of the 2D analytical solution. Relative error (RE) as the function of N for ω₀=0.1, 0.5 and 0.9 at lateral position (a) x = 0.1, (b) x = 0.2, and (c) x = 0.3.

Download Full Size | PDF

It is worth noting that the convergence rate is related to the delta function angular source of the BSF. In many practical applications, the source angular spread is not infinitesimal, and one might approximate the delta function with a generalized delta function, say, of the form:

S (x, z, θ) = δ (x) δ (z) \frac{1}{\sqrt{π} ε} exp (- \frac{θ^{2}}{2 ε^{2}}) .

As ε →0, the original delta function BSF is recovered. For small ε > 0, the converge of Eq. (37) is improved. This is illustrated in Fig. 3 comparing the delta function 2D BSF with that using Eq. (37), for different ε. The deviation between the finite ε and ε=0 results are primarily for smaller x, as expected. As ε is decreased, the result approaches the ε=0 result smoothly; while the number of Fourier modes N required also increases.

Fig. 3 2D scalar irradiance solution at optical depth z = 2 under generalize delta function source with ε=0.06, 0.03, 0.015 and 0.

Download Full Size | PDF

3.2. Validation of 2D analytical BSF

We now validate the present 2D theory for ocean radiative transfer BSF by comparing the closed-form analytical solution with Monte Carlo simulations. We assume a homogeneous ocean with phase function as given earlier [Eq.(45)]. The only two factors that affect the lateral profiles of irradiance are the single scattering albedo ω₀ and optical depth z. Figure 4 shows the comparisons for broad ranges of ω₀ and z. For fixed optical depth z, Fig. 4(a), shows excellent agreement between analytical and simulation results for different ω₀ over all values of x. Similarly, for fixed ω₀ and varying z [Fig. 4(b)] the comparison is excellent.

Fig. 4 Comparisons of the lateral 2D scalar irradiance $E_{s}^{(2)}$ between Monte Carlo simulations and analytical solutions for (a) optical depth z=2, for different single scattering albedos ω₀=0.1, 0.5 and 0.9; (b) ω₀=0.5 for z=1, 2, and 3. All analytical BSFs are obtained using N=320 modes.

Download Full Size | PDF

3.3. Validation of 3D analytical BSF

We now compare the 3D BSF given by the (approximate) analytical form [Eq. (44)] against (3D) Monte Carlo result. We again assume homogeneous ocean with 3D H-G phase function [22]:

Φ (Θ) = \frac{1}{4 π} \frac{1 - g^{2}}{{(1 + g^{2} - 2 g cos Θ)}^{3 / 2}}, g = 0.924,

where Θ is the scattering angle. Figure 5 shows the comparisons for different single scattering albedos ω₀ and optical depths z. The comparisons are overall very satisfactory, showing the general validity of the azimuthal single scattering approximation we propose. For smaller optical depth, the results are especially good (for all lateral distance r) where multiple azimuthal scattering is expected to be limited. For large z, the deviations are primarily for small r, where multiple azimuthal scattering is expected to be more important. Comparing the results for different ω₀, the deviations are somewhat increased for increasing ω₀ (for larger z at small r), again as expected, since greater ω₀ is associated with enhanced multiple azimuthal scattering.

Fig. 5 Comparisons of the lateral 3D scalar irradiance $E_{s}^{(3)}$ between Monte Carlo simulations and analytical solutions for single scattering albedos ω₀=0.1, 0.5 and 0.9 at (a) z=1.05, (b) z=2.89, and (c) z=5.

Download Full Size | PDF

4. Conclusion

The main contribution of the present work is the development of the analytical 3D BSF for homogenous ocean environment using a new single scattering approximation for the azimuthal direction (only). This reduces the 3D beam spread problem to that in 2D. We obtain the exact analytical 2D BSF by combining the Fourier spectral approach and the rotated-frame method [18, 19], and obtaining for the first time an explicit expression in terms of sums of Fourier modes expressed in special functions (without inverse Fourier integrals). The present solution (involving a double Fourier sum) is a substantial improvement over existing small-angle approximations (containing four nested Fourier integrals [14]) in terms of computational simplicity and efficacy.

We compare the new 2D and 3D analytical BSFs against direct Monte Carlo simulations. The results are excellent over broad ranges of optical depths and single scattering albedos. The overall good comparisons for the 3D BSF support the general validity of the azimuthal single scattering approximation, with (small) deviations limited, as expected, to large optical depths and small lateral distances for extremely turbid water (single scattering albedos greater that ∼0.9).

The analytical BSFs are found to be computationally efficient with exponential convergence (error ∼ e^−κN) with increasing number of Fourier modes N (compared with ∼ N^−1/2 convergence rate of standard Monte Carlo simulation using N realizations).

Our results obtained under the assumption of azimuthal single scattering should be valid in the ocean environment where scattering of light is highly forward-peaked. Comparisons and validation of these theoretical results against direct ocean measurements would be very desirable.

Acknowledgments

This study has been supported by the Office of Naval Research (No. N00014-13-1-0352).

References and links

1. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967).

2. S. Q. Duntley, “Underwater lighting by submerged lasers and incandescent sources,” Tech. Rep., Scripps Institution of Oceanography, San Diego, Calif. (1971).

3. K. J. Voss and A. L. Chapin, “Measurement of the point spread function in the ocean,” Appl. Opt. 29, 3638–3642 (1990). [CrossRef] [PubMed]

4. K. J. Voss, “Simple empirical model of oceanic point spread function,” Appl. Opt. 30, 2647–2651 (1991). [CrossRef] [PubMed]

5. W. Hou, D. J. Gray, A. D. Weidemann, and R. A. Arnone, “Comparison and validation of point spread models for imaging in natural waters,” Opt. Express. 16(13), 9958–9965 (2008). [CrossRef] [PubMed]

6. C. D. Mobley, Light and Water: Radiative Transfer In Natural Waters (Academic, 1994).

7. G. N. Plass and G. W. Kattawar, “Monte Carlo calculations of radiative transfer in the earth atmosphere-ocean system: I. flux in the atmosphere and ocean,” J. Phys. Oceanog 2(13), 139–145 (1972). [CrossRef]

8. R. Sanchez and N. J. McCormick, “Analytic beam spread function for ocean optics applications,” Appl. Opt. 41(30), 6276–6288 (2002). [CrossRef] [PubMed]

9. E. P. Zege, A. P. Ivanov, and I. L. Katsev, Image Transfer Through a Scattering Medium (Springer-Verlag, 1991). [CrossRef]

10. G. C. Pomraning, A. K. Prinja, and J. W. VanDenburg, “An asymptotic model for the spreding of a collimated beam,” Nucl. Sci. Eng. 112, 347–360 (1992).

11. C. Börgers and E. W. Larsen, “Asymptotic derivation of the Fermi pencil-beam approximation,” Nucl. Sci. Eng. 123, 343–357 (1996).

12. W. H. Wells, “Theory of small angle scattering,” in Optics of the Sea (Interface and In-Water Transmission and Imaging), P. Halley, ed., Lecture Series no. 61 (North Atlantic Treaty Organization, Advisory Group for Aerospace Research and Development, Neuilly-sur-Seine, France, 1973).

13. R. E. Walker, Marine Light Field Statistics (Wiley, 1994).

14. J. W. McLean, J. D. Freeman, and R. E. Walker, “Beam spread function with time dispersion,” Appl. Opt. 37(21), 4701–4711 (1998). [CrossRef]

15. W. Wei, X. Zhang, Y. Chao, and X. Zhou, “An analytical model of the power spatial distribution for underwater optical wireless communication,” Optica Applicata XLII(1), 157–166 (2002).

16. N. L. Swanson, V. M. Gehman, B. D. Billard, and T. L. Gennaro, “Limits of the small-angle approximation to the radiative transport equation,” J. Opt. Soc. Am. A 18, 385–391 (2001). [CrossRef]

17. N. L. Swanson, B. D. Billard, V. M. Gehman, and T. L. Gennaro, “Application of the small-angle approximation to ocean water types,” Appl. Opt. 40, 3608–3613 (2001). [CrossRef]

18. V. A. Markel, “Modified spherical harmonics method for solving the radiative transport equation,” Waves Random Media 14(1), L13–L19 (2004). [CrossRef]

19. A. Liemert and A. Kienle, “Green’s functions for the two-dimensional radiative transfer equation in bounded media,” J. Phys. A: Math. Theor. 45, 175201 (2012). [CrossRef]

20. Z. Xu and D. K. P. Yue, “Monte Carlo radiative transfer simulation for the near-ocean-surface high-resolution downwelling irradiance statistics,” Opt. Eng. 53, 051408 (2014). doi: [CrossRef]

21. T. J. Petzold, “Volume scattering functions for selected ocean water,” Tech. Rep., Scripps Institution of Oceanography, San Diego, Calif. (1972).

22. L. C. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]

23. L. Florescu, J. C. Schotland, and V. A. Markel, “Single-scattering optical tomography,” Phy. Rev. E 79, 036607 (2009). [CrossRef]

24. J. P. Boyd, Chebyshev and Fourier Spectral Methods (DOVER Publications, 2000).

25. Z. Xu, D. K. P. Yue, L. Shen, and K. J. Voss, “Patterns and statistics of in-water polarization under conditions of linear and nonlinear ocean surface waves,” J. Geophys. Res. 116, C00H12 (2011). doi: [CrossRef]

26. Z. Xu, X. Guo, L. Shen, and D. K. P. Yue, “Radiative transfer in ocean turbulence and its effect on underwater light field,” J. Geophys. Res. 117, C00H18 (2012)., doi: [CrossRef]

Analytical solution of beam spread function for ocean light radiative transfer

Abstract

1. Introduction

2. Methods

2.1. Variables reduction of three-dimensional radiative transfer equation

2.2. Analytical solution of two-dimensional scalar irradiance

3. Results

3.1. Convergence evaluation

3.2. Validation of 2D analytical BSF

3.3. Validation of 3D analytical BSF

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Equations (48)

Optics Express