Validation of the polarized Monte Carlo model of shipborne oceanic lidar returns

Huixin He; Huixin He; Qi Liu; Qi Liu; Junwu Tang; Junwu Tang; Peizhi Zhu; Peizhi Zhu; Shuguo Chen; Xiaoquan Song; Xiaoquan Song; Songhua Wu; Songhua Wu

doi:10.1364/OE.511445

1. Introduction

Ocean lidar has a better penetration into the seawater, up to three times of detection depth than that of passive ocean color sensors [1,2]. It has been applied to detect optical scattering layers [3–6], estimate bio-optical properties [7–9], and observe the diel vertical migration (DVM) of ocean animals [10]. Based on depth-resolved lidar data, the net primary production (NPP) estimation can be improved up to 54% [11]. Since the first underwater polarization study by Waterman [12], researchers have proved that polarimetry can help to understand the optical properties of oceanic particles [13–19]. The results of the polarization measurements suggest that it is possible to distinguish various types of phytoplankton using their polarization characteristics [20–22]. Polarization observation can motivate the development of retrieval algorithms for oceanic particles [23,24]. Although many polarization measurements of seawater have been conducted, the Mueller matrix proposed by Voss & Fry is currently the only publicly available measured matrix [15]. This matrix was obtained as an average of many observations from 10° to 160° in the Gulf of Mexico, Atlantic, and Pacific Oceans. The matrix was parameterized by Kokhanovsky [25], and it can be applied to theoretical polarized radiative transfer.

A variety of methods have been developed for the lidar signals, such as the analytical model [26–28] and Monte Carlo (MC) simulation [6,29–31]. The formula of the analytical model is complex and consists of a variety of independent variables. The direct solution of the equation relies on assuming a variety of ideal conditions, while the approximate solution needs to simplify the transmission equation [32,33]. MC techniques use probability theory and random numbers to simulate the numerous light rays propagating through a medium [34]. This method has been widely applied in scalar and vector radiative transfer for several decades [33,35–40]. Meanwhile, many variance reduction techniques have been adopted to improve the accuracy of simulation results. For example, the strategy of estimating results has changed from traditional photon number statistics to photon weight statistics and then to local estimate statistics (or directional estimate). For a given photon number, the simulated results of the local estimate have less variance than the basic photon weight statistical method [33,41]. Other methods include Russian roulette [42], detector directional importance sampling (DDIS) [43], phase function forward truncation (PFFT) [44], etc. All these techniques can accelerate the convergence rate and improve numerical performance. Some MC models of radiative transfer have been extended to simulate the lidar signal [39,45]. Combined with the measured signals of atmospheric lidar, the polarized Monte Carlo (PMC) techniques have been employed to discriminate the differences among water clouds, ice clouds, and snow [46–48].

Several MC models specific to oceanic lidar have also been proposed [31,49–52]. Limited by the available particle Mueller matrices, these models can only simulate the signal profiles and have not been further applied to discriminate the particles. Most models focus on the normal incidence or scalar signal simulation, which is not capable of simulating the non-normal incidence of polarization signals in three-dimensional space [31,52]. Besides, polarized model validation relies on comparisons of measured signals and depolarization ratios [51]. Such a validation scheme has some limitations. The particle Mueller matrix, which affects the scattering direction and the depolarization ratio, is a crucial input variable for simulation. The single scattering depolarization ratio of the Voss & Fry Mueller matrix is 0.1173 in theory [15,25,53]. However, if the matrix is not applicable to the study area, the objective differences between measurement and simulation will increase remarkably. So, such validation will not be sufficient. Compared with airborne or spaceborne lidar systems, shipborne lidars have no advantage in observation range, and the measurements are easily disturbed by sea surface conditions. Nevertheless, the shipborne platform can simultaneously observe a large amount of seawater optical profile information. This is important to interpret the relationship between particle properties and the polarization lidar signal.

Previous research has introduced our dual-wavelength, variable-FOV, and polarized oceanic lidar LOOP [54]. In this work, we develop a PMC model for shipborne lidar and aim at model validation. The model is verified by the Rayleigh scattering test case, lidar equation, and in-situ LOOP measurements. Rayleigh scattering is used to validate the polarization theoretical calculations for MC. The lidar equation and LOOP measurements are applied to validate the simulations of single scattering signals and total scattering signals, respectively. The simulation of the underwater light field is also illustrated. The average relative errors of the simulation results demonstrate that this polarized model is accurate and effective.

2. Methodology

This PMC model uses the meridian plane method for tracking the status of polarized light [Fig. 1(a)] [55,56], the coordinate system of lidar simulation is shown in Fig. 1(b). To improve accuracy and efficiency, we use the local estimate technique in signal estimation. Pool based on the local estimate method, proposed a semi-analytic Monte Carlo radiative transfer model (SALMON), which has been widely applied in oceanographic lidar simulation [30,31,50,52]. For the simulation research of the atmosphere, the local estimate technique is used in combination with the variance reduction methods [41,43,57]. Our current studies neglect the influence of the atmosphere, and the detection geometric distance is greater than 0 m. Therefore, neither the problems of the singular term [Eq. (23), the first term], nor the more precise variance reduction methods need to be considered.

Fig. 1. (a) Flow chart of the PMC model, (b) global coordinate system of lidar simulation in this study, and (c) geometry of polarization scattering in the meridian plane method of PMC model.

Download Full Size | PDF

2.1 Polarized Monte Carlo model

Usually, the incident light can be described by the Jones vector or Stokes vector. However, the Stokes vector can also represent partially polarized light. The Stokes vector is defined as [55,58]:

(1)$$\textbf{S} = \left[{\begin{array}{c} I\\ Q\\ U\\ V \end{array}} \right] = \sqrt {\varepsilon /{\mu _\textrm{m}}} \left[{\begin{array}{c} {{E_\parallel }E_\parallel^\ast{+} {E_ \bot }E_ \bot^\ast }\\ {{E_\parallel }E_\parallel^\ast{-} {E_ \bot }E_ \bot^\ast }\\ {{E_\parallel }E_ \bot^\ast{+} {E_ \bot }E_\parallel^\ast }\\ {i({E_\parallel }E_ \bot^\ast{-} {E_ \bot }E_\parallel^\ast )} \end{array}} \right],$$

where $\sqrt {\varepsilon /{\mu _\textrm{m}}}$ is a constant factor, E_|| and E_⊥ are the components of the electric field vector parallel and perpendicular to the reference plane, respectively. The initial Stokes vector for the polarized laser beam is S₀ = [1,1,0,0].

The initial direction of the photon is defined by the incident zenith angle θ₀ and the incident azimuth angle φ₀ [Fig. 1(b)]. The optical depth τ that a photon can propagate is sampled by:

(2)$$\tau = \ln ( - q),$$

here q is a random number that is uniformly distributed over (0,1). The geometric distance Dist is calculated as:

(3)$$\left\{ \begin{array}{l} \tau \textrm{ = }\sum\limits_i {{\tau_i} = \sum\limits_i {{c_i}{d_i}} } \\ Dist = \sum\limits_i {{d_i}} \end{array} \right.,$$

where i indicates the layer that the photon will travel through, d_i and c_i are the geometric distance and attenuation coefficient in each layer, respectively. The position of the photon P is updated by the geometric distance Dist and direction unit vector ξ_in [59]:

(4)$$\textbf{P} = \textbf{P} + {\boldsymbol{\mathrm{\xi}}_{in}} \cdot Dist,$$

(5)$${\boldsymbol{\mathrm{\xi}}_{in}}\textrm{ = }\left\{ \begin{array}{l} {\mu_x} = \sin {\theta_1}\cos {\varphi_1}\\ {\mu_y} = \sin {\theta_1}\sin {\varphi_1}\\ {\mu_z} = \cos {\theta_1} \end{array} \right.,(0 \le {\theta _1} \le \pi ,0 \le {\varphi _1} \le 2\pi ),$$

where θ₁ and φ₁ are the polar angle and azimuth angle for the incident light [Fig. 1(c)]. When a scattering event occurs, the photon’s weight w is updated by multiplying the single scattering albedo ω₀,

(6)$$w = w{\omega _0}.$$

And the Stokes vector becomes [33,39,60,61]:

(7)$${\textbf{S}_{sca}} = \frac{\textbf{Z}}{{{M_{11}}}}{\textbf{S}_{in}}\textrm{ = }\frac{{\textbf{R}({\gamma _2})\textbf{M}({\theta _{sca}})\textbf{R}({\gamma _1})}}{{{M_{11}}({\theta _{sca}})}}{\textbf{S}_{in}},$$

where S_sca and S_in are incident and scattering Stokes, and S_sca will be the incident Stokes vector S_in for the next scattering event. Z is the scattering matrix, M(θ_sca) is the scattering phase matrix (Mueller matrix), M₁₁(θ_sca) is the (1,1) element of the matrix, and R(γ) is the rotation matrix:

(8)$$\textbf{R}(\gamma ) = \left[{\begin{array}{cccc} 1&0&0&0\\ 0&{\cos (2\gamma )}&{ - \sin (2\gamma )}&0\\ 0&{\sin (2\gamma )}&{\cos (2\gamma )}&0\\ 0&0&0&1 \end{array}} \right].$$

The rotation of polarization scattering is in Fig. 1(c). There are two schemes for obtaining the scattering polar angle θ_sca, rotation angles γ₁, γ₂, and scattering direction unit vector ξ_sca. One is the scalar sample method [33,39,62], and the other is the rejection method [56,60,63]. Previous research has introduced these sampling methods [64]. This work applies the scalar method for angle sampling. The scalar method is based on the importance sampling theory [33,65], and it is normalized by M₁₁ to remove the bias [Eq. (7)].

First, use M₁₁ as a probability distribution function (PDF) for sampling scattering angle θ_sca. The typical ocean particle scattering matrix proposed by Voss and Fry [15] is:

(9)$${\boldsymbol{\tilde{\textbf{M}}}_{VF}}({\theta _{sca}}) = \left[ \begin{array}{cccc} 1&{a_1}&0&0\\ {a_1}&{a_2}&0&0\\ 0&0&{a_3}&0\\ 0&0&0&{a_3} \end{array}\right],$$

(10)$$\left\{ \begin{array}{l} {a_1} = \frac{{ - p({{90}^ \circ }){{\sin }^2}{\theta_{sca}}}}{{1 + p({{90}^ \circ }){{\cos }^2}{\theta_{sca}}}}\\ {a_2} = \frac{{p({{90}^ \circ })[1 + {{\cos }^2}({\theta_{sca}} - {\theta_0})] + \zeta \exp ( - \alpha {\theta_{sca}})}}{{1 + p({{90}^ \circ }){{\cos }^2}({\theta_{sca}} - {\theta_0}) + \zeta \exp ( - \alpha {\theta_{sca}})}}\\ {a_3} = \frac{{2p({{90}^ \circ })\cos {\theta_{sca}} + \zeta \exp ( - \alpha {\theta_{sca}})}}{{1 + p({{90}^ \circ }){{\cos }^2}{\theta_{sca}} + \zeta \exp ( - \alpha {\theta_{sca}})}} \end{array} \right..$$

The elements a₁, a₂, and a₃ are calculated by Kokhanovsky [25]. This matrix is a reduced Mueller matrix, which has normalized the elements by M₁₁. So, we adopt the Petzold average-particle phase function as M₁₁ and sampling angles [35,66]. The scattering azimuth angle φ_sca is randomly selected from 0 to 2π. The scattering direction unit vector ξ_sca = [μ_x_′, μ_y_′, μ_z_′] is calculated by:

(11)$$\left[\begin{array}{l} \mu_x^{\prime}\\ \mu_y^{\prime}\\ \mu_z^{\prime} \end{array} \right] = \left\{ \begin{array}{l} \left[ {\begin{array}{ccc} {\frac{{{\mu_x}{\mu_z}}}{{\sqrt {1 - \mu_z^2} }}}&{\frac{{ - {\mu_y}}}{{\sqrt {1 - \mu_z^2} }}}&{{\mu_x}}\\ {\frac{{{\mu_y}{\mu_z}}}{{\sqrt {1 - \mu_z^2} }}}&{\frac{{{\mu_x}}}{{\sqrt {1 - \mu_z^2} }}}&{{\mu_y}}\\ { - \sqrt {1 - \mu_z^2} }&0&{{\mu_z}} \end{array}} \right]\left[{\begin{array}{c} {\sqrt {1 - {{\cos }^2}{\theta_{sca}}} \cos {\varphi_{sca}}}\\ {\sqrt {1 - {{\cos }^2}{\theta_{sca}}} \sin {\varphi_{sca}}}\\ {\cos {\theta_{sca}}} \end{array}} \right],({\cos^2}{\theta_{sca}} < 1)\\ sign({\mu_z})\left[{\begin{array}{c} {\sqrt {1 - {{\cos }^2}{\theta_{sca}}} \cos {\varphi_{sca}}}\\ {\sqrt {1 - {{\cos }^2}{\theta_{sca}}} \sin {\varphi_{sca}}}\\ {\cos {\theta_{sca}}} \end{array}} \right],({\cos^2}{\theta_{sca}} \approx 1) \end{array} \right..$$

The scattering direction unit vector ξ_sca will be used as the new direction unit vector ξ_in for photon moving. Finally, obtain the rotation angles γ₁ and γ₂ as follows [62]:

(12)$$\left\{ \begin{array}{l} {\gamma_\textrm{1}}\textrm{ = arccos[(}{\boldsymbol{\mathrm{\xi}}_{in}} \times {\textbf{e}_\textbf{z}}\textrm{)} \cdot \textrm{(}{\boldsymbol{\mathrm{\xi}}_{in}} \times {\boldsymbol{\mathrm{\xi}}_{sca}}\textrm{)]}\\ {\gamma_2}\textrm{ = arccos[(}{\boldsymbol{\mathrm{\xi}}_{sca}} \times {\textbf{e}_\textbf{z}}\textrm{)} \cdot \textrm{(}{\boldsymbol{\mathrm{\xi}}_{sca}} \times {\boldsymbol{\mathrm{\xi}}_{in}}\textrm{)]} \end{array} \right.,$$

where e_z is the unit vector of the z-axis. After N_S scattering events, the Stokes vector is expressed as:

(13)$$\textbf{S}_{sca}^{Ns} = \left( {\prod\limits_{i = 1}^{Ns} {\frac{{{\textbf{Z}_i}}}{{{{({M_{11}})}_i}}}} } \right){\textbf{S}_\textbf{0}}.$$

When the photon hits the ocean surface, reflection or refraction needs to be considered. The new direction unit vector is updated by the reflected or transmitted angles [59]. The reflectance matrix R_f and transmission matrix T_f are:

(14)$${\textbf{R}_\textbf{f}} = \left[ {\begin{array}{cccc} {\frac{1}{2}({r_\parallel }r_\parallel^ \ast{+} {r_ \bot }r_ \bot^ \ast )}&{\frac{1}{2}({r_\parallel }r_\parallel^ \ast{-} {r_ \bot }r_ \bot^ \ast )}&0&0\\ {\frac{1}{2}({r_\parallel }r_\parallel^ \ast{-} {r_ \bot }r_ \bot^ \ast )}&{\frac{1}{2}({r_\parallel }r_\parallel^ \ast{+} {r_ \bot }r_ \bot^ \ast )}&0&0\\ 0&0&{\Re \{ {r_\parallel }r_ \bot^ \ast \} }&{\Im \{ {r_\parallel }r_ \bot^ \ast \} }\\ 0&0&{ - \Im \{ {r_\parallel }r_ \bot^ \ast \} }&{\Re \{ {r_\parallel }r_ \bot^ \ast \} } \end{array}} \right],$$

(15)$${\textbf{T}_\textbf{f}} = {f_T}\left[ {\begin{array}{cccc} {\frac{1}{2}({t_\parallel }t_\parallel^ \ast{+} {t_ \bot }t_ \bot^ \ast )}&{\frac{1}{2}({t_\parallel }t_\parallel^ \ast{-} {t_ \bot }t_ \bot^ \ast )}&0&0\\ {\frac{1}{2}({t_\parallel }t_\parallel^ \ast{-} {t_ \bot }t_ \bot^ \ast )}&{\frac{1}{2}({t_\parallel }t_\parallel^ \ast{+} {t_ \bot }t_ \bot^ \ast )}&0&0\\ 0&0&{\Re \{ {t_\parallel }t_ \bot^ \ast \} }&{\Im \{ {t_\parallel }t_ \bot^ \ast \} }\\ 0&0&{ - \Im \{ {t_\parallel }t_ \bot^ \ast \} }&{\Re \{ {t_\parallel }t_ \bot^ \ast \} } \end{array}} \right].$$

The elements of the reflectance matrix R_f and transmission matrix T_f are calculated as follows:

(16)$$\left\{ \begin{array}{l} {r_\parallel }r_\parallel^ \ast{=} {\left( {\frac{{\cos {\theta_{in}} - ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}{{\cos {\theta_{in}} + ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}} \right)^2}\\ {r_ \bot }r_ \bot^ \ast{=} {\left( {\frac{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} - \cos {\theta_{tra}}}}{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} + \cos {\theta_{tra}}}}} \right)^2}\\ \Re \{ {r_\parallel }r_ \bot^ \ast \} \textrm{ = }\left( {\frac{{\cos {\theta_{in}} - ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}{{\cos {\theta_{in}} + ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}} \right)\left( {\frac{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} - \cos {\theta_{tra}}}}{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} + \cos {\theta_{tra}}}}} \right)\\ \Im \{ {r_\parallel }r_ \bot^ \ast \} \textrm{ = 0} \end{array} \right.,$$

(17)$$\left\{ \begin{array}{l} {t_\parallel }t_\parallel^ \ast{=} {\left( {\frac{{2({n_{in}}/{n_{tra}})\cos {\theta_{in}}}}{{\cos {\theta_{in}} + ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}} \right)^2}\\ {t_ \bot }t_ \bot^ \ast{=} {\left( {\frac{{2({n_{in}}/{n_{tra}})\cos {\theta_{in}}}}{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} + \cos {\theta_{tra}}}}} \right)^2}\\ \Re \{ {t_\parallel }t_ \bot^ \ast \} \textrm{ = }\left( {\frac{{2({n_{in}}/{n_{tra}})\cos {\theta_{in}}}}{{\cos {\theta_{in}} + ({n_{in}}/{n_{tra}})\cos {\theta_{tra}}}}} \right)\left( {\frac{{2({n_{in}}/{n_{tra}})\cos {\theta_{in}}}}{{({n_{in}}/{n_{tra}})\cos {\theta_{in}} + \cos {\theta_{tra}}}}} \right)\\ \Im \{ {t_\parallel }t_ \bot^ \ast \} \textrm{ = 0} \end{array} \right.,$$

For the total internal reflection,

(18)$$\left\{ \begin{array}{l} {r_\parallel }r_\parallel^ \ast{=} \textrm{1}\\ {r_ \bot }r_ \bot^ \ast{=} \textrm{1}\\ \Re \{ {r_\parallel }r_ \bot^ \ast \} \textrm{ = }\frac{{\textrm{2}{{\sin }^\textrm{4}}{\theta_{in}}}}{{1 - [1 + {{({n_{tra}}/{n_{in}})}^2}]{{\cos }^2}{\theta_{in}}}} - 1\\ \Im \{ {r_\parallel }r_ \bot^ \ast \} \textrm{ = }\frac{{\textrm{2}\cos {\theta_{in}}{{\sin }^\textrm{2}}{\theta_{in}}\sqrt {{{\sin }^\textrm{2}}{\theta_{in}} - {{({n_{tra}}/{n_{in}})}^2}} }}{{1 - [1 + {{({n_{tra}}/{n_{in}})}^2}]{{\cos }^2}{\theta_{in}}}} \end{array} \right.,$$

In simulation, we use the transmission matrix for irradiance [50,58], so the f_T factor is:

(19)$${f_T} = \frac{{{n_{tra}}\cos {\theta _{tra}}}}{{{n_{in}}\cos {\theta _{in}}}},$$

where n_in and n_tra are the refractive indexes of the incident and transmitted medium, respectively. The reflection coefficient R_f(1,1) will determine whether the photon is reflected or transmitted. If a random number q < R_f(1,1), the photon is reflected; otherwise, it is transmitted. The Stokes vector is updated by Eq. (7), with the matrix M replaced by matrix R_f or T_f.

2.2 Estimation of detected signal profiles

For airborne or spaceborne lidar, when the incident zenith angle of the laser beam is small, it can be modeled as a normal incidence. If a photon gets Ns scattering events, and its position P_Ns is in the detector’s field of view (FOV), the photon’s signal contribution of this position will be estimated before the new scattering event [Figs. 1(a) and 2(a)]. For a photon in seawater, this condition is expressed as [52,67]:

(20)$$\left\{ \begin{array}{l} {\textbf{P}_{Ns}}{(x)^2} + {\textbf{P}_{Ns}}{(y)^2} \le {A_r}(z)\\ {A_r}(z) = {[H \cdot \tan (\frac{{{\theta _{FOV}}}}{2}) + {R_{tele}}+\frac{{\left| {{\textbf{P}_{Ns}}(z)} \right| \cdot \sin (\frac{{{\theta_{FOV}}}}{2})}}{{\sqrt {n_w^2 - {{\sin }^2}(\frac{{{\theta_{FOV}}}}{2})} }}]^2} \end{array} \right.,$$

where θ_FOV is the FOV in the unit of rad, R_tele is the radius of the telescope, and n_w is the refractive index of seawater.

Fig. 2. Observation geometry of PMC model simulating lidar signals, (a) nadir observation, (b) oblique observation (for shipborne lidar). The green dot P_Ns represents the position of the photon at the Ns^th scattering event, the green line $\boldsymbol{\mathrm{\xi}}_{in}^{Ns}$ shows the photon’s incident direction in this scattering event. The blue line is the estimated trajectory of the photon signal, the dotted line $\boldsymbol{\mathrm{\xi}}_{up}^{Ns}$ is the underwater estimated trajectory, and the solid line $\boldsymbol{\mathrm{\xi}}_{det}^{Ns}$ represents the atmospheric estimated trajectory. The blue dot $\textbf{P}_{surf}^{Ns}$ indicates the sea surface interaction position for this scattering event.

Download Full Size | PDF

Shipborne lidar will emit the laser beam at a tilt angle to avoid strong reflected signals from the sea surface. With a conical field-of-view, an elliptical patch of sea surface is viewed by the detector. The oblique observation geometry is shown in Fig. 2(b). If the photon’s interaction position with the surface $\textbf{P}_{surf}^{Ns}$ is in this ellipse area, the photon can be detected. Such a relationship can be expressed as [67,68]:

(21)$$\left\{ \begin{array}{l} \frac{{\textbf{P}_{surf}^{Ns}{{(x)}^2}}}{{R_a^2}} + \frac{{\textbf{P}_{surf}^{Ns}{{(y)}^2}}}{{R_b^2}} \le 1\\ {R_a} = [H \cdot \tan (\frac{{{\theta_{FOV}}}}{2})/\cos {\theta_v} + {R_{tele}}]/\cos {\theta_v}\\ {R_b} = H \cdot \tan (\frac{{{\theta_{FOV}}}}{2})/\cos {\theta_v} + {R_{tele}} \end{array} \right.,$$

where θ_v is oblique observation angle, the R_a and R_b are the long axis and short axis of the ellipse, respectively. The detector’s position is P_det = [-Htan(θ_v),0,H]. For the simplified calculation, the detector’s approximate position P′_det = [-Htan(θ_v),0,n_wH] is used for obtaining the upwelling direction unit vector $\boldsymbol{\mathrm{\xi}}_{up}^{Ns}$ and the sea surface interaction position $\textbf{P}_{surf}^{Ns}$ [28,57]:

(22)$$\left\{ \begin{array}{l} \boldsymbol{\mathrm{\xi}}_{up}^{Ns}\textrm{ = }\frac{{\textbf{P}_{det}^{\prime} - {\textbf{P}_{Ns}}}}{{|{\textbf{P}_{det}^{\prime} - {\textbf{P}_{Ns}}} |}}\\ \textbf{P}_{surf}^{Ns}\textrm{ = }{\textbf{P}_{Ns}}\textrm{ + }\left|{\frac{{{\textbf{P}_{Ns}}(z)}}{{\boldsymbol{\mathrm{\xi}}_{up}^{Ns}(z)}}} \right|\boldsymbol{\mathrm{\xi}}_{up}^{Ns} \end{array} \right..$$

If $\textbf{P}_{surf}^{Ns}$ satisfies Eq. (21), the underwater signal contribution of this photon in the Ns^th scattering event is estimated as [39,50,61]:

(23)$$\begin{array}{r} {\textbf{S}_D}({z_{det}},{\textbf{P}_{Ns}}) = \frac{{{A_{tele}}}}{{{{|{\textbf{P}{^{\prime}_{det}} - {\textbf{P}_{Ns}}} |}^\textrm{2}}}}\left( {\prod\limits_{i = 1}^{Ns} {{\omega_0}({\textbf{P}_i})} } \right){\textbf{Z}_{Est}}({\textbf{P}_{Ns}})\exp [ - {\tau _{atm}}(\textbf{P}_{surf}^{Ns},{\textbf{P}_{det}})]\\ \cdot \exp [ - {\tau _{wat}}({\textbf{P}_{Ns}},\textbf{P}_{surf}^{Ns})] \times \frac{{\textbf{Z}({\textbf{P}_{Ns - 1}})}}{{{M_{11}}({\textbf{P}_{Ns - 1}})}} \cdots \frac{{\textbf{Z}({\textbf{P}_\textbf{1}})}}{{{M_{11}}({\textbf{P}_\textbf{1}})}}{\textbf{S}_\textbf{0}}, \end{array}$$

(24)$${\textbf{Z}_{Est}}({\textbf{P}_{Ns}}) = \textbf{R}(\gamma _2^{\prime}){\textbf{T}_\textbf{f}}({\theta _{in}},{\theta _{tra}})\textbf{R}(\gamma _1^{\prime}) \cdot \textbf{R}({\gamma _2})\textbf{M}(\theta {^{\prime}_{sca}})\textbf{R}({\gamma _1}),$$

(25)$${z_{det}} = \frac{1}{2}[|{\textbf{P}_{surf}^{Ns} - {\textbf{P}_{Ns}}} |\textrm{ + }\sum\limits_{i = 1}^{Ns} {Dist} ({\textbf{P}_i})],$$

where A_tele is the detected area of the telescope, θ′_sca is the scattering angle for estimation [Fig. 2(a)], τ_atm and τ_wat are the optical depth in the atmosphere and the seawater, respectively. Here we neglect the atmospheric attenuation, so τ_atm = 0. For homogeneous water, τ_wat =cd_wat, c is the beam attenuation coefficient of the water medium (Section 2.3), and d_wat is the distance from P_Ns to $\textbf{P}_{surf}^{Ns}$. For inhomogeneous water, τ_wat=Σc(i)d_wat(i), c(i) are the beam attenuation coefficients of water medium layers, and d_wat(i) are the layer distances from P_Ns to $\textbf{P}_{surf}^{Ns}$. z_det is the estimated water depth of the return signal, and it mainly depends on the cumulative geometric distance of N_S scattering events that the photon has traveled in water. It should be noted that if the scattering event N_S =1, the signal contribution will belong to single scattering signals. Sum all the photons’ single scattering contributions, and the simulated single scattering results are:

(26)$${\textbf{S}_{single}}({z_{det}}) = \frac{1}{{{N_P}}}\sum\limits_{i = 1}^{{N_P}} {{\textbf{S}_D}(i,{\textbf{P}_{Ns = 1}},{z_{det}})} .$$

where N_P is total photon number of simulation. Sum all the photons’ scattering contributions, the simulated total signal results are:

(27)$${\textbf{S}_{total}}({z_{det}}) = \frac{1}{{{N_P}}}\sum\limits_{i = 1}^{{N_P}} {\sum\limits_{Ns = 1}^\infty {{\textbf{S}_D}(i,{\textbf{P}_{Ns}},{z_{det}})} } .$$

The parallel-channel and cross-channel signal results are calculated by the I and Q components of S_single or S_total:

(28)$$\left\{ \begin{array}{l} P{P_{weight}}({z_{det}}) = {I_f}({z_{det}}) + {Q_f}({z_{det}})\\ C{P_{weight}}({z_{det}}) = {I_f}({z_{det}}) - {Q_f}({z_{det}}) \end{array} \right..$$

After normalizing the results by the sea surface point detection signal, the depolarization ratio can be expressed by the normalized parallel-channel and cross-channel signals:

(29)$$Dp({z_{det}}) = \frac{{C{P_N}({z_{det}})}}{{P{P_N}({z_{det}})}}.$$

The relative error used in this work is defined as follows:

(30)$$\delta \textrm{ = }\frac{{|{{X_S}\textrm{ - }{X_R}} |}}{{{X_R}}} \times 100\%.$$

To facilitate the error comparison, this study takes positive values of the absolute errors for statistics. X_S is the simulated value, and X_R is the measured or reference value.

2.3 IOPs model of seawater

The inherent optical properties (IOPs) of seawater are used as input variables in the simulation. The total seawater absorption coefficient a_total and scattering coefficient b_total can be expressed as:

(31)$$\left\{ \begin{array}{l} {a_{total}} = {a_p} + {a_g} + {a_w}\\ {b_{total}} = {b_p} + {b_w} \end{array} \right.,$$

where the absorption and scattering coefficients of pure seawater are a_w(532 nm) = 0.043m⁻¹ and b_w(532 nm) = 0.00223m⁻¹, respectively [69,70]. a_p and a_g are the absorption coefficients of phytoplankton and colored dissolved organic matter (CDOM). b_p is the scattering coefficients of phytoplankton. These parameters can be calculated by the in-situ measured chlorophyll a concentration [Chl_a]. Specific in-situ calculations refer to previous research [54].

3. Model validation

Comparing the simulation model results with the lidar measured signals is generally used for model validation. Due to the limitations of the oceanic particle Mueller matrix, there may be some objective differences between measurement and simulation. To ensure the correctness of the simulation model, we first use the Rayleigh scattering case for theoretical polarization calculation validation and then compare the simulated signals with the lidar equation results and the shipborne lidar measurements.

3.1 Rayleigh scattering for vectorization validation

The Rayleigh scattering test case is common in vector radiation transfer (VRT) theoretical calculations [62,71,72]. To validate the vectorization (polarization), such as Stokes calculation and reference plane rotation, we simulate Rayleigh scattering in atmospheres and compare it with the benchmark results of Natraj [73]. The conditions of this test case are described as follows: set a single layer of scattering molecules (non-absorbing) as the medium. The optical thickness of the layer is 0.5, and the bottom albedo is 0. The cosine of the solar zenith angle is 0.92, the cosine of the viewing zenith angle is 0.4, and the relative azimuth angles of the viewing direction are from 0° to 180° with an increment of 30°. The Mueller matrix of Rayleigh scattering is defined as follows:

(32)$${\textbf{M}_{\textbf{Ray}}}\textrm{ = }\frac{\Delta }{{4\pi }}\left[ {\begin{array}{cccc} {\frac{\textrm{3}}{\textrm{4}}( \textrm{1 + co}{\textrm{s}^2}\theta ) }&{\frac{{ - 3}}{4}{{\sin }^2}\theta }&0&0\\ {\frac{{ - 3}}{4}{{\sin }^2}\theta }&{\frac{\textrm{3}}{\textrm{4}}( \textrm{1 + co}{\textrm{s}^2}\theta ) }&0&0\\ 0&0&{\frac{3}{2}\cos \theta }&0\\ 0&0&0&{\Delta ^{\prime}\frac{3}{2}\cos \theta } \end{array}} \right] + (1 - \Delta )\left[ {\begin{array}{cccc} 1&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right],$$

where Δ=(1-δ)/(1+δ/2) and Δ′=(1-2δ)/(1-δ). Here the depolarization factor δ is 0, so Δ=Δ′=1; and this matrix can be further simplified. The (1,1) element of matrix M_Ray is the Rayleigh phase function ${\tilde{\beta }_{R\textrm{ay}}}(\theta )$, which satisfies $\textrm{2}\pi \int_\textrm{0}^\pi {{{\tilde{\beta }}_{R\textrm{ay}}}(\theta )} \sin \theta d\theta = 1$. This case displays horizontal invariance, so for a photon getting through N scattering events, the estimation for the detector can be simplified as [33]:

(33)$$\begin{array}{r} {\textbf{S}_D}({z_D},{\boldsymbol{\mathrm{\xi}}_D}) = \frac{{|{{\textbf{e}_\textbf{z}} \cdot {\boldsymbol{\mathrm{\xi}}_1}} |}}{{|{{\textbf{e}_\textbf{z}} \cdot {\boldsymbol{\mathrm{\xi}}_D}} |}}\sum\limits_{N = 1}^\infty {\left( {\prod\limits_{i = 1}^N {{\omega_0}({z_i})} } \right)} \textbf{Z}({z_N};{\boldsymbol{\mathrm{\xi}}_N} \to {\boldsymbol{\mathrm{\xi}}_D})\exp [ - \tau ({z_N},{z_D})]\\ \times \frac{{\textbf{Z}({z_{N - 1}};{\boldsymbol{\mathrm{\xi}}_{N - 1}} \to {\boldsymbol{\mathrm{\xi}}_N})}}{{{M_{11}}({z_{N - 1}};{\boldsymbol{\mathrm{\xi}}_{N - 1}} \to {\boldsymbol{\mathrm{\xi}}_N})}} \cdots \frac{{\textbf{Z}({z_1};{\boldsymbol{\mathrm{\xi}}_\textrm{1}} \to {\boldsymbol{\mathrm{\xi}}_2})}}{{{M_{11}}({z_1};{\boldsymbol{\mathrm{\xi}}_1} \to {\boldsymbol{\mathrm{\xi}}_2})}}{\textbf{S}_\textbf{0}}, \end{array}$$

where S₀ = [1,0,0,0] is the incident Stokes, ξ₁ is the incident direction, ξ_D is the viewing direction, ξ_N is the scattering direction, z_D is the position of the detector, and z_i is the scattering position. For the non-absorbing layer, the single scattering albedo ω₀ = 1.0. The matrix Z = R(γ₂)M_Ray(θ)R(γ₁), M₁₁ is the (1,1) element of matrix M_Ray.

The result of the comparison is shown in Fig. 3. There is a good agreement between PMC and the Rayleigh scattering benchmark results. Using 10⁷ photons simulation, the relative errors of the Stokes components are less than 0.07%.

Fig. 3. Comparing PMC simulation and Rayleigh scattering benchmark results.

Download Full Size | PDF

3.2 Validation with lidar equation and in-situ lidar returns

The lidar equation is often applied to simulate lidar signal profiles, but it cannot simulate multiple scattering or polarization. PMC simulation solves both of these problems. Neglecting the effects of multiple scattering, the oblique observation lidar signals can be described by the lidar equation [74,75]:

(34)$${N_{\textrm{signal}}}(z)\textrm{ = }\frac{{{E_\textrm{0}}}}{{hv}} \cdot \frac{{{A_{tele}}{{\cos }^2}{\theta _v}}}{{{{({n_w}H + z)}^2}}}O(z){T_{\textrm{atm}}}^2T_{surf}^2{\eta _{oe}}{\eta _{qe}}\Delta z{\beta _\pi }(z)\exp [\frac{{ - 2\int_0^z {{K_{lidar}}(z^{\prime})dz^{\prime}} }}{{\cos {\theta _w}}}],$$

where N_signal(z) is the detector received signal, E₀ is the transmitted pulse energy, v is the frequency of the laser, h is the Planck constant, A_tele is the detected area of the telescope, O(z) is the lidar overlap function, T_atm and T_surf are the transmittance of the atmosphere and sea surface, η_oe and η_qe are the optical efficiency and quantum efficiency of the system, H is the distance from lidar system to the surface, z is the transmitted length in water, Δz is the vertical resolution. The zenith angle of laser beam in the seawater θ_w is calculated by the observation angle θ_v using Snell law. β_π(z) is the volume scattering coefficient at the angle π:

(35)$${\beta _\pi }(z) = {\tilde{\beta }_w}(\pi ){b_w} + {\tilde{\beta }_p}(\pi ){b_p}(z),$$

${\tilde{\beta }_\textrm{w}}$ and ${\tilde{\beta }_p}$ are the phase functions of seawater and particles, respectively. The lidar attenuation coefficient [76]:

(36)$${K_{lidar}}(z) \approx {K_d}(z) + [{c_{total}}(z) - {K_d}(z)]\exp [{ - 0.85{c_{total}}(z)D(z)} ],$$

D is the lidar spot diameter of FOV, and the diffuse attenuation coefficient K_d (z) = 0.0452 + 0.0474 [Chl_a]^0.67 [77,78]. The total attenuation coefficient of seawater c_total = a_total + b_total is calculated based on the concentration profiles of Chl_a (Section 2.3).

The shipborne lidar LOOP carried out a field experiment in the Western Pacific from October 23, 2019, to December 10, 2019 (Fig. 4). The trial area is an oligotrophic ocean that belongs to the western end of the North Equatorial Circulation and the southern birthplace of the Kuroshio [54]. A series of in-situ instruments were deployed to provide ocean optical parameter profiles. The Chl_a concentration was measured by the Scattering Meter ECO BB-9 (WET Labs Inc.) [54]. We select eight stations for validation. Table 1 displays the parameters for the lidar system.

Fig. 4. (a) Shipborne oceanic lidar LOOP, (b) the distribution of the selected stations in the Western Pacific

Download Full Size | PDF

Table 1. The system parameters of LOOP

View Table

The Chl_a concentration of this trial area observed by MODIS was perennially less than 0.06 mg/m³ [54]. This low concentration led to some inaccurate in-situ measurements of the chlorophyll profiles, such as negative values, the chlorophyll data need to be corrected before simulation. Considering the small FOV of shipborne lidar, it is inferred that the signals from the upper water body (z < 30 m) are mainly based on single scattering effects. Therefore, the results of the lidar equation should be highly consistent with the LOOP measurement signal (z < 30 m), and the lidar equation can be used to determine the benchmark concentrations of chlorophyll. Meanwhile, the corrected profile also refers to the chlorophyll concentrations of BGC Argo in the trial area to ensure the rationality of the correction [79]. Figure 5 shows the concentration profiles of chlorophyll after correction and fitting.

Fig. 5. Chlorophyll concentration profiles from the field experiments.

Download Full Size | PDF

According to Section 2.3, the seawater parameters used in the PMC simulation can be calculated based on these chlorophyll profiles. The simulation lidar system parameters are shown in Table 1. Considering the beam divergence is less than 0.5 mrad, use a collimated beam for PMC simulation. The simulated photon number is 10⁷. To compare with scalar lidar equation results, calculate the scalar signals of PMC and LOOP based on their parallel-channel and cross-channel signals, respectively. All the signals have been normalized and shown in Figs. 6 and 7.

Fig. 6. Comparing lidar equation signals (lidarEq) with single scattering scalar signals of PMC (PMC_SS) and LOOP’s scalar signals (LOOP_S). The left sub-graph is the signal comparison and the right sub-graph shows the relative errors of lidar equation with LOOP_S or PMC_SS, respectively. The sampling interval of error statistics is 2 m, the arrows represent that the error value is more than 50%, and MRE is the mean relative error.

Download Full Size | PDF

Fig. 7. Comparing simulated total scattering signals of PMC (PP_PMC, CP_PMC) with measured LOOP signals (PP_LOOP, CP_LOOP). The left sub-graph is the signal comparison and the right is the relative errors of PMC and LOOP, respectively. PP is the parallel-channel signal, and CP is the cross-channel signal, the sampling interval of error statistics is 2 m, the arrows represent that the error value is more than 50%, and MRE is the mean relative error.

Download Full Size | PDF

Figure 6 shows that the single scattering scalar signals of PMC (PMC_SS) are roughly consistent with the lidar equation returns (lidarEq), the mean relative errors (MRE) of these stations range from 6.23% to 30.94%. However, the MREs of the lidar equation results and LOOP’s scalar signals (LOOP_S) are from 44.22% to 138.58%, and as the depth increases, the relative errors gradually increase; some relative errors below 40 m are even larger than 50.00%. These imply that using the lidar equation to simulate the measured signal is insufficient, it cannot represent the effect of multiple scattering signals as the depth increases. These also verify the theory that the simulated signal from the lidar equation is mainly based on single-scattering effect. The total scattering signal simulations of PMC are in good agreement with LOOP, with the maximal MREs of parallel-channel and cross-channel being 33.29% and 22.37%, respectively (Fig. 7). The simulated depolarization ratio is close to the theoretical single scattering depolarization ratio Dp = 0.1173, and the maximal MRE of the depolarization ratio is 24.13% (Fig. 8). Those results prove the validity of the PMC model for the estimation of lidar signals. Based on the analysis above, we infer that the measured signals of LOOP mainly come from the single scattering of the particles, and the influence of multiple scattering signals becomes significant as the depth increases.

Fig. 8. Comparison of depolarization ratio between PMC and LOOP, the left sub-graph is the depolarization ratio, the right sub-graph is the relative error, the sampling interval of error statistics is 2 m, and MRE is the mean relative error.

Download Full Size | PDF

However, there are some unsatisfactory results in PMC simulation. For simulated total signals and depolarization ratios (Figs. 7 and 8), some stations’ relative errors at 0-10 m or 30-50 m are significant. For the upper water body, the increase in error may be related to the influence of solar background light, sea waves, seawater bubbles, etc. These factors are not considered in the current simulation results. In addition, sea waves and bubbles can also affect the measurement of seawater optical parameters. For deep water, seawater particles may be the main factor. As the depth increases, the chlorophyll concentration changes, and so do the seawater optical parameters. At present, the Voss & Fry Mueller matrix is the only available measured particle Mueller matrix. It may be inappropriate to simulate the detection of the inhomogeneous seawater profile by only using one Mueller matrix.

To observe the underwater light field of photons, we also build a light field recorder to count the photon numbers in the simulation. This recorder has 206 × 206 × 252 sub-recorders, and each of them has a size of 0.2 × 0.2 × 0.2 m. The recorder result is first normalized by the simulation photon number, then multiplied by the number of photons emitted (calculated by laser beam energy). Figure 9 describes the propagation process of the laser beam. As the depth of penetration increases, the center of the beam gradually deviates from the center of the x-y plane, and the energy of the beam gradually spreads around.

Fig. 9. The underwater light field simulation, the color bar represents the recorded photon numbers.

Download Full Size | PDF

4. Discussion and conclusions

Simulation models are extremely significant for interpreting polarized lidar signals. In this study, a polarized Monte Carlo model for shipborne lidar is developed. This model is based on the meridian plane method and the local estimate technique. We use the Rayleigh scattering test case for basic theoretical polarization validation, and the relative errors are less than 0.07%. The maximal mean relative error of simulated single scattering scalar signals and lidar equation results is 30.94%. The maximal MREs of total scattering signals and LOOP returns in parallel-channel and cross-channel are 33.29% and 22.37%, respectively, and the maximal MRE of the depolarization ratio is 24.13%. These results demonstrate that this simulation model is effective. Combined with the comparison of the signals, and depolarization ratios, it can be inferred that the signals of the shipborne lidar LOOP mainly come from the single scattering of particles. This also indicates the PMC model helps understand the signal contributions from polarized multiple scattering and single scattering. The underwater light field of the oblique incidence laser beam is also displayed, which illustrates the phenomenon of beam diffusion in the water.

In addition, the analysis also shows that it is difficult to fully verify the model by only comparing the simulated signal with the measured signal. Many factors can influence the simulation results, such as the measurements of chlorophyll profiles, the hydrosol model of seawater, and the rough sea surface influence. For the oceanic lidar with a large FOV, the relative error between the simulation and the measured signal will be more significant, which may cover up the theoretical calculation error of the simulation model itself. Therefore, it is necessary to introduce other effective validation methods. The refined validations, including more environmental influence factors, will be performed in future research.

Funding

Laoshan Laboratory (LSKJ202201202); National Key Research and Development Program of China (2022YFB3901705).

Acknowledgments

The authors appreciate the BGC-Argo program. We are grateful to Xiaogang Xing from the Second Institute of Oceanography, Ministry of Natural Resources, China, for discussing the chlorophyll profiles of BGC-Argo.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. M. Ottaviani, R. Foster, A. Gilerson, et al., “Airborne and shipborne polarimetric measurements over open ocean and coastal waters: Intercomparisons and implications for spaceborne observations,” Remote Sens. Environ. 206, 375–390 (2018). [CrossRef]

2. C. Jamet, A. Ibrahim, Z. Ahmad, et al., “Going Beyond Standard Ocean Color Observations: Lidar and Polarimetry,” Front. Mar. Sci. 6(251), 1 (2019). [CrossRef]

3. A. P. Vasilkov, Y. A. Goldin, B. A. Gureev, et al., “Airborne polarized lidar detection of scattering layers in the ocean,” Appl. Opt. 40(24), 4353–4364 (2001). [CrossRef]

4. J. H. Churnside and P. L. Donaghay, “Thin scattering layers observed by airborne lidar,” ICES J. Mar. Sci. 66(4), 778–789 (2009). [CrossRef]

5. M. J. Behrenfeld, Y. Hu, R. T. O’Malley, et al., “Annual boom-bust cycles of polar phytoplankton biomass revealed by space-based lidar,” Nat. Geosci. 10(2), 118–122 (2017). [CrossRef]

6. G. M. Krekov, M. M. Krekova, and V. S. Shamanaev, “Laser sensing of a subsurface oceanic layer II Polarization characteristics of signals,” Appl. Opt. 37(9), 1589–1595 (1998). [CrossRef]

7. B. L. Collister, R. C. Zimmerman, C. I. Sukenik, et al., “Remote sensing of optical characteristics and particle distributions of the upper ocean using shipboard lidar,” Remote Sens. Environ. 215, 85–96 (2018). [CrossRef]

8. M. J. Behrenfeld, Y. Hu, C. A. Hostetler, et al., “Space-based lidar measurements of global ocean carbon stocks,” Geophys. Res. Lett. 40(16), 4355–4360 (2013). [CrossRef]

9. J. H. Churnside, J. W. Hair, C. A. Hostetler, et al., “Ocean backscatter profiling using high-spectral-resolution lidar and a perturbation retrieval,” Remote Sens. 10(12), 2003 (2018). [CrossRef]

10. M. J. Behrenfeld, P. Gaube, A. Della Penna, et al., “Global satellite-observed daily vertical migrations of ocean animals,” Nature 576(7786), 257–261 (2019). [CrossRef]

11. J. A. Schulien, M. J. Behrenfeld, J. W. Hair, et al., “Vertically- resolved phytoplankton carbon and net primary production from a high spectral resolution lidar,” Opt. Express 25(12), 13577–13587 (2017). [CrossRef]

12. T. H. Waterman, “Polarization Patterns in Submarine Illumination,” Science 120(3127), 927–932 (1954). [CrossRef]

13. A. Ivanoff, N. Jerlov, and T. H. Waterman, “A comparative study of irradiance, beam transmittance and scattering in the sea near Bermuda,” Limnol. Oceanogr. 6(2), 129–148 (1961). [CrossRef]

14. G. W. Kattawar, G. N. Plass, and J. A. Guinn, “Monte Carlo Calculations of the Polarization of Radiation in the Earth’s Atmosphere-Ocean System,” J. Phys. Oceanogr. 3(4), 353–372 (1973). [CrossRef]

15. K. J. Voss and E. S. Fry, “Measurement of the Mueller matrix for ocean water,” Appl. Opt. 23(23), 4427–4439 (1984). [CrossRef]

16. K. J. Voss and N. Souaidia, “POLRADS: polarization radiance distribution measurement system,” Opt. Express 18(19), 19672–19680 (2010). [CrossRef]

17. Y. Zhao, C. Poulin, D. McKee, et al., “A closure study of ocean inherent optical properties using flow cytometry measurements,” J. Quant. Spectrosc. Radiat. Transf. 241, 106730 (2020). [CrossRef]

18. M. Chami and M. D. Platel, “Sensitivity of the retrieval of the inherent optical properties of marine particles in coastal waters to the directional variations and the polarization of the reflectance,” J. Geophys. Res. 112(C5), C05037 (2007). [CrossRef]

19. J. Chowdhary, B. Cairns, F. Waquet, et al., “Sensitivity of multiangle, multispectral polarimetric remote sensing over open oceans to water-leaving radiance: Analyses of RSP data acquired during the MILAGRO campaign,” Remote Sens. Environ. 118, 284–308 (2012). [CrossRef]

20. Ø Svensen, J. J. Stamnes, M. Kildemo, et al., “Mueller matrix measurements of algae with different shape and size distributions,” Appl. Opt. 50(26), 5149–5157 (2011). [CrossRef]

21. H. Volten, J. F. De Haan, J. W. Hovenier, et al., “Laboratory measurements of angular distributions of light scattered by phytoplankton and silt,” Limnol. Oceanogr. 43(6), 1180–1197 (1998). [CrossRef]

22. Y. Wang, J. Dai, R. Liao, et al., “Characterization of physiological states of the suspended marine microalgae using polarized light scattering,” Appl. Opt. 59(5), 1307–1312 (2020). [CrossRef]

23. M. E. Zugger, A. Messmer, T. J. Kane, et al., “Optical scattering properties of phytoplankton: Measurements and comparison of various species at scattering angles between 1° and 170°,” Limnol. Oceanogr. 53(1), 381–386 (2008). [CrossRef]

24. D. Koestner, D. Stramski, and R. A. Reynolds, “Characterization of suspended particulate matter in contrasting coastal marine environments with angle-resolved polarized light scattering measurements,” Appl. Opt. 60(36), 11161–11179 (2021). [CrossRef]

25. A. A. Kokhanovsky, “Parameterization of the Mueller matrix of oceanic waters,” J. Geophys. Res. Ocean. 108(C6), 10–13 (2003). [CrossRef]

26. Y. Zhou, W. Chen, X. Cui, et al., “Validation of the analytical model of oceanic lidar returns: Comparisons with Monte Carlo simulations and experimental results,” Remote Sens. 11(16), 1870 (2019). [CrossRef]

27. R. E. Walker and J. W. McLean, “Lidar equations for turbid media with pulse stretching,” Appl. Opt. 38(12), 2384–2397 (1999). [CrossRef]

28. Y. I. Kopilevich and A. G. Surkov, “Mathematical modeling of the input signals of oceanological lidars,” J. Opt. Technol. 75(5), 321–326 (2008). [CrossRef]

29. P. Bruscaglioni, G. Zaccanti, S. del Bianco, et al., “Monte Carlo for multiple scattering and nonspherical particles,” Proc. SPIE 5237, 223–227 (2004). [CrossRef]

30. Q. Liu, D. Liu, J. Bai, et al., “Relationship between the effective attenuation coefficient of spaceborne lidar signal and the IOPs of seawater,” Opt. Express 26(23), 30278–30291 (2018). [CrossRef]

31. L. R. Poole, D. D. Venable, and J. W. Campbell, “Semianalytic Monte Carlo radiative transfer model for oceanographic lidar systems,” Appl. Opt. 20(20), 3653–3656 (1981). [CrossRef]

32. G. Zhou, C. Li, D. Zhang, et al., “Overview of Underwater Transmission Characteristics of Oceanic LiDAR,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 14, 8144–8159 (2021). [CrossRef]

33. H. R. Gordon, Physical Principles of Ocean Color Remote Sensing (University of Miami, 2019), Chap.2. [CrossRef]

34. C. D. Mobley, The Oceanic Optics Book (International Ocean Colour Coordinating Group (IOCCG), 2022), Appendix E.

35. C. D. Mobley, B. Gentili, H. R. Gordon, et al., “Comparison of numerical models for computing underwater light fields,” Appl. Opt. 32(36), 7484–7504 (1993). [CrossRef]

36. J. T. O. Kirk, “Estimation of the scattering coefficient of natural waters using underwater irradiance measurements,” Mar. Freshw. Res. 32(4), 533–539 (1981). [CrossRef]

37. C. Emde, R. Buras, and B. Mayer, “ALIS: An efficient method to compute high spectral resolution polarized solar radiances using the Monte Carlo approach,” J. Quant. Spectrosc. Radiat. Transf. 112(10), 1622–1631 (2011). [CrossRef]

38. D. Zawada, G. Franssens, R. Loughman, et al., “Systematic comparison of vectorial spherical radiative transfer models in limb scattering geometry,” Atmos. Meas. Tech. 14(5), 3953–3972 (2021). [CrossRef]

39. P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution to the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47(8), 1037–1047 (2008). [CrossRef]

40. D. Ramon, F. Steinmetz, D. Jolivet, et al., “Modeling polarized radiative transfer in the ocean-atmosphere system with the GPU-accelerated SMART-G Monte Carlo code,” J. Quant. Spectrosc. Radiat. Transf. 222–223, 89–107 (2019). [CrossRef]

41. B. Mayer, “Radiative transfer in the cloudy atmosphere,” EPJ Web Conf. 1, 75–99 (2009). [CrossRef]

42. H. Iwabuchi, “Efficient Monte Carlo methods for radiative transfer modeling,” J. Atmos. Sci. 63(9), 2324–2339 (2006). [CrossRef]

43. R. Buras and B. Mayer, “Efficient unbiased variance reduction techniques for Monte Carlo simulations of radiative transfer in cloudy atmospheres: The solution,” J. Quant. Spectrosc. Radiat. Transf. 112(3), 434–447 (2011). [CrossRef]

44. H. Iwabuchi and T. Suzuki, “Fast and accurate radiance calculations using truncation approximation for anisotropic scattering phase functions,” J. Quant. Spectrosc. Radiat. Transf. 110(17), 1926–1939 (2009). [CrossRef]

45. B. Mayer, F. Faure, L. Bugliaro, et al., “Realistic Simulations of Earthcare Observations,” in Proceedings of the 2nd ESA EarthCARE Workshop, 10–12 June, 2009.

46. Y. You, G. W. Kattawar, P. Yang, et al., “Sensitivity of depolarized lidar signals to cloud and aerosol particle properties,” J. Quant. Spectrosc. Radiat. Transf. 100(1-3), 470–482 (2006). [CrossRef]

47. Y. Hu, X. Lu, X. Zeng, et al., “Deriving Snow Depth From ICESat-2 Lidar Multiple Scattering Measurements,” Front. Remote Sens. 3, 855159 (2022). [CrossRef]

48. C. Jimenez, A. Ansmann, R. Engelmann, et al., “The dual-field-of-view polarization lidar technique: A new concept in monitoring aerosol effects in liquid-water clouds - Theoretical framework,” Atmos. Chem. Phys. 20(23), 15247–15263 (2020). [CrossRef]

49. P. Bruscaglioni, G. Zaccanti, and Q. Wei, “Transmission of a pulsed polarized light beam through thick turbid media: numerical results,” Appl. Opt. 32(30), 6142–6150 (1993). [CrossRef]

50. P. G. Stegmann, B. Sun, J. Ding, et al., “Study of the effects of phytoplankton morphology and vertical profile on lidar attenuated backscatter and depolarization ratio,” J. Quant. Spectrosc. Radiat. Transf. 225, 1–15 (2019). [CrossRef]

51. Q. Liu, X. Cui, W. Chen, et al., “A semianalytic Monte Carlo radiative transfer model for polarized oceanic lidar: Experiment-based comparisons and multiple scattering effects analyses,” J. Quant. Spectrosc. Radiat. Transf. 237, 106638 (2019). [CrossRef]

52. P. Chen, D. Pan, Z. Mao, et al., “Semi-analytic Monte Carlo radiative transfer model of laser propagation in inhomogeneous sea water within subsurface plankton layer,” Opt. Laser Technol. 111, 1–5 (2019). [CrossRef]

53. Q. Liu, X. Cui, C. Jamet, et al., “A semianalytic monte carlo simulator for spaceborne oceanic lidar: Framework and preliminary results,” Remote Sens. 12(17), 1–11 (2020). [CrossRef]

54. Q. Liu, S. Wu, B. Liu, et al., “Shipborne variable-FOV, dual-wavelength, polarized ocean lidar: design and measurements in the Western Pacific,” Opt. Express 30(6), 8927–8948 (2022). [CrossRef]

55. S. Chandrasekhar, Radiative Transfer (Dover Publications, 1960), Chap.1.

56. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef]

57. G. I. Marchuk, G. A. Mikhailov, M. A. Nazaraliev, et al., The Monte Carlo Methods in Atmospheric Optics (Springer Berlin, 1980), Chap.5.

58. C. D. Mobley, HydroPol Mathematical Documentation:Invariant Imbedding Theory for the Vector Radiative Transfer Equation (2014).

59. R. A. Leathers, C. O. D. Trijntje Valerie Downes, and C. D. Mobley, Monte Carlo Radiative Transfer Simulations for Ocean Optics : A Practical Guide (2004).

60. H. H. Tynes, G. W. Kattawar, E. P. Zege, et al., “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40(3), 400–412 (2001). [CrossRef]

61. Y. Hu, D. Winker, P. Yang, et al., “Identification of cloud phase from PICASSO-CENA lidar depolarization: A multiple scattering sensitivity study,” J. Quant. Spectrosc. Radiat. Transf. 70(4-6), 569–579 (2001). [CrossRef]

62. C. Emde, R. Buras, B. Mayer, et al., “The impact of aerosols on polarized sky radiance: Model development, validation, and applications,” Atmos. Chem. Phys. 10(2), 383–396 (2010). [CrossRef]

63. S. Bartel and A. H. Hielscher, “Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,” Appl. Opt. 39(10), 1580–1588 (2000). [CrossRef]

64. H. He, M. Shi, J. Tang, et al., “Scattering direction sampling methods for polarized Monte Carlo simulation of oceanic lidar,” Appl. Opt. 62(23), 6253–6263 (2023). [CrossRef]

65. R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo Method (John Wiley & Sons, 2016), Chap.5.

66. T. J. Petzold, Volume Scattering Functions for Selected Ocean Waters (1972).

67. T. Yin, N. Lauret, and J. P. Gastellu-Etchegorry, “Simulation of satellite, airborne and terrestrial LiDAR with DART (II): ALS and TLS multi-pulse acquisitions, photon counting, and solar noise,” Remote Sens. Environ. 184, 454–468 (2016). [CrossRef]

68. J. P. Gastellu-Etchegorry, T. Yin, N. Lauret, et al., “Simulation of satellite, airborne and terrestrial LiDAR with DART (I): Waveform simulation with quasi-Monte Carlo ray tracing,” Remote Sens. Environ. 184, 418–435 (2016). [CrossRef]

69. Z. Lee, J. Wei, K. Voss, et al., “Hyperspectral absorption coefficient of “pure” seawater in the range of 350–550 nm inverted from remote sensing reflectance,” Appl. Opt. 54(3), 546–558 (2015). [CrossRef]

70. X. Zhang and L. Hu, “Scattering by pure seawater at high salinity,” Opt. Express 17(15), 12685–12691 (2009). [CrossRef]

71. J. Chowdhary, P. W. Zhai, F. Xu, et al., “Testbed results for scalar and vector radiative transfer computations of light in atmosphere-ocean systems,” J. Quant. Spectrosc. Radiat. Transf. 242, 106717 (2020). [CrossRef]

72. C. Emde, V. Barlakas, C. Cornet, et al., “IPRT polarized radiative transfer model intercomparison project - Phase A,” J. Quant. Spectrosc. Radiat. Transf. 164, 8–36 (2015). [CrossRef]

73. V. Natraj, K. F. Li, and Y. L. Yung, “Rayleigh scattering in planetary atmospheres: Corrected tables through accurate computation of X and y functions,” Astrophys. J. 691(2), 1909–1920 (2009). [CrossRef]

74. C. A. Hostetler, M. J. Behrenfeld, Y. Hu, et al., “Spaceborne lidar in the study of marine systems,” Ann. Rev. Mar. Sci. 10(1), 121–147 (2018). [CrossRef]

75. J. H. Churnside, “Review of profiling oceanographic lidar (erratum),” Opt. Eng. 56(7), 079802 (2017). [CrossRef]

76. J. H. Churnside, “Review of profiling oceanographic lidar,” Opt. Eng. 53(5), 051405 (2013). [CrossRef]

77. J. H. Churnside, J. M. Sullivan, and M. S. Twardowski, “Lidar extinction-to-backscatter ratio of the ocean,” Opt. Express 22(15), 18698–18706 (2014). [CrossRef]

78. A. Morel and S. Maritorena, “Bio-optical properties of oceanic waters: A reappraisal,” J. Geophys. Res. 106(C4), 7163–7180 (2001). [CrossRef]

79. H. C. Bittig, T. L. Maurer, J. N. Plant, et al., “A BGC-Argo Guide: Planning, Deployment, Data Handling and Usage Mission Considerations or the Global,” Front. Mar. Sci. 6, 502 (2019). [CrossRef]

Wavelength (nm)	532	Aperture diameter (m)	0.06
Pluse energy (mJ)	200	FOV (mrad)	60
Beam divergence (mrad)	≤0.5	Sampling frequency (ns)	1
Plane Height (m)	13	Oblique observation angle (°)	15

Validation of the polarized Monte Carlo model of shipborne oceanic lidar returns

Abstract

1. Introduction

2. Methodology

2.1 Polarized Monte Carlo model

2.2 Estimation of detected signal profiles

2.3 IOPs model of seawater

3. Model validation

3.1 Rayleigh scattering for vectorization validation

3.2 Validation with lidar equation and in-situ lidar returns

4. Discussion and conclusions

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (9)

Tables (1)

Equations (36)

Optics Express