A Novel Fast Multiple-Scattering Approximate Model for Oceanographic Lidar

: An effective lidar simulator is vital for its system design and processing algorithms. However, laser transmission is a complex process due to the effects of sea surface and various interactions in seawater such as absorption, scattering, and so on. It is sophisticated and difﬁcult for multiple scattering to accurately simulate. In this study, a multiple-scattering lidar model based on multiple-forward-scattering-single-backscattering approximation for oceanic lidar was proposed. Compared with previous analytic models, this model can work without assuming a homogeneous water and ﬁxed scattering phase function. Besides, it takes consideration of lidar system and environmental parameters including receiver ﬁeld of view, different scattering phase functions, particulate sizes, stratiﬁed water, and rough sea surface. One should note that because the scattering phase function is difﬁcult to determine accurately, the simulation accuracy may be reduced in a complex oceanic environment. The Cox–Munk model used in our method simulates capillarity waves but ignores gravity waves, and the pulse stretching is not included. The wide-angle scattering occurs in the dense subsurface phytoplankton, which sometimes makes it hard to use this model. In this study, we ﬁrstly derived this method based on an analytical solution by convolving Gaussians of the forward-scattering contribution of layer dr and the energy density at R in the small-angle-scattering approximation. Then, the effects of multiple scattering and water optical properties were analyzed using the model. Meanwhile, the validation with Monte Carlo model was implemented. Their coefﬁcient of determination is beyond 0.9, the RMSE is within 0.02, the MAD is within 0.02, and the MAPD is within 8%, which indicates that our model is efﬁcient for oceanographic lidar simulation. Finally, we studied the effects of FOV, SPF, rough sea surface, stratiﬁed water, and particle size. These results can provide reference for the design of the oceanic lidar system and contribute to the processing of lidar echo signals.


Introduction
Ocean optical properties are significant indicators of the marine ecological environment [1] and their detection is vital for upper-ocean research [2]. There are several commonly used detection methods, but each has its advantages and disadvantages. There is no doubt that in situ measurements have the highest accuracy and are essential for calibration and validation of remote sensing [1], but they are limited by poor efficiency [3]. Passive ocean color remote sensing can obtain global-scale optical data efficiently, but most of the ocean color signals emanate from the first optical depths of the upper ocean [4]. Sonars can profile the vertical structure of seawater, but they have to be deployed underwater due to the nearly total reflection of the air-water interface [5]. Oceanic lidar could address these limitations due to its depth-resolved capacity, and has been widely used for plankton layer detection [6][7][8][9], ocean optical profiling [10][11][12], etc. Presently, there is an increasing demand for the development of oceanic lidar techniques to meet the needs of ocean detection [13].
An effective and accurate lidar simulation is necessary for the development of the oceanographic lidar system. The widely used single scattering lidar equation ignores the effect of multiple scattering, while it brings great errors and difficulties to signal processing [14]. Fortunately, several methods have been proposed to deal with multiple scattering. Most of them are based on standard Monte Carlo (MC) simulations [15,16], semi-analytic MC simulations [17][18][19][20], and quasi-single small-angle approximation (QSAA) analytic methods [21][22][23][24][25]. MC simulations provide a quite accurate and flexible result by using a statistical approach that consumes considerable computing power and time because of random error noise [17]. Besides, MC simulations contribute little to the solution of inversion problems [14]. Previous QSAA methods run fast and offer a robust result but are only suitable for several specific cases, such as ignoring wide-angle scattering caused by dense subsurface phytoplankton [17] and using fixed phase function [22].
Overall, an accurate model with fast speed and wider applicability is needed. In this paper, we introduce a novel, fast multiple-scattering approximate model for ocean lidar based on Gaussian laser beam. Compared with previous analytic models, our model can deal with stratified water, and it applies to any water scattering phase function (SPF). The model is in good agreement with MC simulation results, which demonstrates its efficiency. Then, the effects of multiple scattering, SPF, and stratified water were analyzed. These results can provide reference to the design of the oceanographic lidar system and contribute to the processing of lidar echo signals.
The structure of this article is as follows. In Section 2, we derive this method based on an analytical solution by convolving Gaussians of the forward-scattering contribution of layer dr and the energy density at R in the small-angle-scattering approximation. In Section 3, we present a study of the effects of multiple scattering and water optical properties. The simulation results are subsequently compared with MC simulation and measurement. In Section 4, we discuss the impacts of FOV, SPF, stratified water, wind-driven rough sea surface, and particulate size. In Section 5, we summarize our findings.

Multiple-Scattering Solution Based on QSAA
As depicted in Figure 1, a laser beam is received by a detector after multiple scattering in the water, where scattering at small angles θ i is dominant. Under the framework of the quasi-single-scattering small-angle approximation (QSAA) [26], the received power is mainly composed of the multiple forward scatterings at small angles and a single backscattering at a large angle of nearly 180 • . Figure 1. The geometry of multiple scattering of lidar propagation. 2Θ is full-angle receiver field of view (FOV); θ i is forward scattering angle; θ b is backscattering angle; Θ and D are half-angle receiver field of view and position of the "equivalent" lidar system due to refraction. For a source located at a point with a radius vector R 0 , the laser beam first propagates to a position R towards the direction n by multiple forward scattering. Then, the laser beam travels to a position R in the direction −n after one large-angle backscattering at the position R in the direction −n b . Finally, the return signals are collected by the receiver at the position R . The echo signals are given by [27]: where S src (R 0 , n 0 ) is the normalized source radiance distribution with dR 0 dn 0 S src (R 0 , n 0 ) = E 0 , G(R, n; R 0 , n 0 ) is the Green function for the radiative transfer equation, b(R) is the scattering coefficient, p b (R; −n b , n) is the scattering phase function (SPF) for the large backscattering angle (−n b , n) between vectors n and −n b , and S rec is receiver angular-spatial pattern of sensitivity.
Using the reciprocity principle that states G R , −n ; R, −n b = G R, n b ; R , n and defining the source and receiver-source radiances that describe the multiply scattered radiance at the point due to the true source S src and receiver-source S rec src , respectively, I rec src (R, n) = dR 0 dn 0 S rec src (R 0 , n 0 )G(R, n; R 0 , n 0 ) where S rec src R , n = S rec R , −n is defined as a "receiver-source" which can be referred to as a radiance distribution of a fictitious source that is equivalent to the real receiver.
Then, we can obtain the echo signal power scattered at position R of the form: Equation (4) simplifies the complex lidar problem to two simple problems of light propagation that are connected by a single backscattering event.
For the light propagation problem in the conditions of QSAA, the solutions in the Fourier space are [28]: I src (z, q, p) = S src (q, p + qz) G(z, q, p) I rec src (z, q, p) = S rec src (q, p + qz) G(z, q, p) (5) where the Fourier transforms have forms: where r and n ⊥ are the projections of vectors R and n onto the plane z = constant, respectively. For the nearly backward scattering, the phase function p b (R; −n b , n) depends on the difference |n ⊥ − n ⊥b | because of small angles of n ⊥ and n ⊥b . Here p b (z; −n b , n) ≈ p b (z; π−|n ⊥b − n ⊥ |) . Substituting n ⊥ = n ⊥b − n ⊥ , Equation (4) can be rewritten as follows using SAA variables: Integrating Equation (7) over r in the plane z, the lidar return from depth z is approximated as follows: Then, the integral from Equation (8) can be transformed into: where S(z, r, n ⊥ ) is the true receiver function S rec (z, r, n ⊥ ), and I e (z, r, n ⊥ ) is the effective radiance that is transmitted to the position (z, r) from the true source S src (0, r, n ⊥ ) in the effective medium with optical parameters defined as follows: c e (z) = 2c(z), b e (z) = 2b(z), p e (z; π− n ⊥ ) = p(z; π− n ⊥ ) (10)

Analytic Model Based on Convolutions of Gaussian Energy Density Functions
Equation (9) gives the general solution of return lidar signals, but it is still hard to use in practice. Here we derive a more practical model based on a Gaussian dependence of the transmitted laser beam on the divergence angle.
As shown in Figure 2a, for a source located at a point with a radius vector R 0 , the number of incident photos dN in the increment area rdrdϕ on a slab of thickness c∆t/2 at a distance R without scattering is given by [29]: where N 0 is the number of photos in the emitted pulse, θ 0 is half of the divergence of the laser beam, τ(R) = R 0 c(z )dz is the optical depth between the laser and the slab, and c(z ) is the beam attenuation coefficient as the function of z.
Photons arrive at the slab c∆t/2 after a small-angle scattering event that occurs on the layer dz 1 at a distance z 1 from the slab c∆t/2 (as shown in Figure 2c). If the scattering phase function is a Gaussian function of the scattering angle, the angle-spatial distribution of incident photons on the slab can be calculated as the convolution of the Gaussian laser beam with the Gaussian phase function, which has: where γ(z) is the scattering energy in the forward peak of the scattering phase function at z and is γ(z) = πΘ 2 s p(0, z 1 )/4π, p(θ, z 1 )/4π is the scattering phase function at z 1 , Θ 2 s is the mean square angular width of the scattering phase function.
Then, consider photons incident on the slab after additional scatterings at distances z 2 , z 3 , . . . , z m between z 1 and the slab. The number of incident photos dN m in the increment area rdrdϕ on the slab after m scattering can be calculated as the convolution of the Gaussian laser beam with m Gaussian phase functions, which has: The beam backscattered from the area rdrdϕ in the slab c∆t/2 is collected by the receiver after an additional scattering at the layer dz m+1 at distance z m+1 is expressed as follows ( Figure 2b): where A r is the area of the telescope, 4π is the average backscatter phase function near 180 • for nth-order scattering, b(R) is the scattering coefficient, θ is the angle between the beam path and the radius vector, and φ is the angle of the radius vector.
Photons are collected by the receiver after additional scatterings at z m+2 , z m+3 , . . . , z n−1 between z m+1 and the telescope, as shown in Figure 2d. The angular distribution of the photons can be calculated as the convolution of the Gaussian phase function, which has: After integrating over positions of the increment area rdrdϕ, the angular distributions of photons collected by a telescope can be expressed as follows: where 0 < ρ < ρ t , and ρ t is half of the receiver FOV. Then, integrate over the FOV and divide by the signal scattering, which has: dP n (R) Finally, the ratio between the signal power of nth-order scattering and the signalscattering can be calculated by integrating over dz as follows:

Sea Surface Modeling
The Cox and Munk sea surface probability distribution function [30,31] is used here to approximate wind-driven rough sea surface modeling: where µ n = cos(θ n ), θ n is the polar angle, and ϕ n is the azimuth angle for the wave facet normal vector → n . σ 2 is the variance as a function of the wind speed W (m/s) above the sea surface.

Effects of Multiple Scattering
Multiple scattering has a great influence on the lidar echo signals, so the effects of multiple scattering need to be determined for accurate lidar signal processing. To address the multiple scattering during lidar transmission, we simulated the ratios of double, three-order, and four-order scattering to single scattering, and their signal intensity is simulated, as well ( Figure 3). The parameters of the lidar system used are H = 300 m, FOV = 10 mrad, and the receiver aperture Ar = 0.09 m 2 . The water optical properties of clear ocean shown in Table 1 are used here. As shown in Figure 3, when the laser beam enters the water, the echo signal is dominated by single scattering, while the contribution of multiple scattering is negligible. However, multiple scattering signal increases with the propagation of laser beam in the water. In particular, when the depth exceeds 27 m, the double scattering signal begins to be more than single scattering signal, and at this depth range, the dominating backscattered signal received by the lidar FOV is multiple scattering. However, the total signal magnitude decreases as the depth increases because of water attenuation. Some backward signal exists out of the FOV after huge amounts of multiple scattering, which causes the decreasing of the signal magnitude as well. The difference between the total echo signals and the single scattering signal is shown in the result, and it increases with the propagation of the laser beam, which indicates that the multiple scattering effect is not negligible for lidar inversion.

Effects of Water Optical Property
To analyze the effects of different seawater optical properties, echo signals in three typical types of ocean water (Table 1) are simulated based on the same lidar system parameters as in Section 3.1. As shown in Figure 4, with the increases of water turbidity, the ratio of multiple scattering to single scattering becomes greater (Figure 4a-c). Meanwhile, due to the relatively high attenuation of the LiDAR signal with range, the lidar echo signal intensity in turbid water decreases faster, resulting in quite shallow penetration depth. Besides, the strong multiple scattering results in a more complex signal inversion. Therefore, a high signal-to-noise performance and an accurate inversion algorithm are needed for the detection of turbid water. The existent lidar inversion approaches and algorithms were often based on the solution form of the lidar equation with single scattering hypothesis. While multiple scattering leads to the estimated lidar attenuation coefficient value falling between the diffuse attenuation coefficient and the beam attenuation coefficient, the contribution of multiple scattering is considered as a hindrance for these algorithms. We still need to build a relationship between the estimated lidar attenuation coefficient and diffuse attenuation coefficient under different FOVs after solving the lidar equation with single scattering hypothesis. Hence, an improved multiple scattering lidar equation could be constructed and solved for higher accuracy based on a deeper understanding and a more accurate simulation of multiple scattering. The single-scattering lidar equation S = Aβ(π)exp(−2K lidar z) is then transformed to S = Aβ(π)exp(−2K lidar z) * F(z) [22]. F(z) represents multiple scattering effect function, which could be solved by the accurate simulation under different lidar conditions. We then could obtain the relationship between the estimated lidar attenuation coefficient and diffuse attenuation coefficient expressed as K lidar = K d − 1 2z lnF(z). Besides, accurate simulation of multiple scattering could be helpful for interpreting the depolarization introduced by multiple scattering, which in turn helps in the inversion of particles backscattering based on depolarization. In addition, simulation training data could be produced based on simulating under different water environment parameter conditions, then a signal matching algorithm (SMA) could be applied in the future, based on finding the proximate simulation signal compared to the measurement signal. This algorithm has the advantage of obtaining multiple water parameters simultaneously. A theoretical assessment of the stability of the algorithms to uncertainties could also be carried out. In addition, a tentative signal-to-noise (S/N) scheme for ocean lidar design which, comprehensively considering the inversion accuracy requirement, hardware cost, technical complexity, and so on, could be proposed based on a comprehensive lidar simulation system.

Validation under Extreme Condition
To validate the model, we simulated the results in extreme condition with a small FOV of 0.1 mrad. The used lidar system parameters are the same as those in Section 3.1 and the water type is the clear ocean. As shown in Figure 5, the multiple scattering signals are far below the single scattering, so the total scattering results meet with the single scattering results very well. The good agreement indicates the validity of the model.

Validation Using MC Model
Then, we used a semi-analytic MC model [17] to validate our simulation. Figure 6 represents the comparison results between our analytical model and MC simulation with several FOVs of 0.1 and 10 mrad under the same lidar system parameter settings as in Section 3.1. The MC model used a total of 10 million photons. The results show that our analytic results are in good agreement with the MC results, indicating the effectiveness of our model. The comparative statistical analysis between them is shown in Table 2. The coefficient of determination (R 2 ) is beyond 0.9, the mean absolute difference (MAD) is within 0.02, the root mean square error (RMSE) is within 0.02, and the mean absolute percent difference (MAPD) is within 8%, which indicates that our model works well for oceanographic lidar simulation.

Validation with In Situ Measurements
In situ measurement is used to validate the model. The data were collected at the station of 109.8164 • E, 18.3144 • N off Wuzhihzou Island on 30 September 2017, with the airborne oceanographic lidar designed by the Shanghai Institute of Optics and Fine Mechanics (SIOM), Chinese Academy of Sciences, Shanghai, China. The wavelength of 532 nm is used for oceanic detection, and the FOV is 6 mrad for shallow waters [7,10]. The water is quite clear, and the beam attenuation coefficient is less than 0.08 m −1 . The chlorophyll concentration is less than 0.17 mg/m 3 and is nearly vertical homogeneous. The comparison result is shown in Figure 7. As the R 2 is 0.97, RMSE is 0.004, MAD is 0.004, and MAPD is 8.9%, our simulation result is in good agreement with the measurement, which indicates that our model works well for real ocean conditions.

Effects of Multiple Scattering Contribution in Different FOVs
The results in Section 3.1 have shown the critical impact of multiple scattering on lidar return signals. Further analysis is needed to determine multiple scattering contributions in lidar echo signals with different FOVs to guide the lidar system design and inversion algorithms. Therefore, we simulated the multiple scattering under different FOVs in different seawater types based on the same lidar system parameters setting, as shown in Section 3.1. The water types used here are clear ocean and coastal water. Figure 8 shows the results with FOVs of 10, 1, and 0.1 mrad in clear water. Figure 8a-d represent the double scattering, three-order scattering, four-order scattering, and the total echo signals, respectively. For a very small FOV of 0.1 mrad, the corresponding multiple scattering is so small that the total echo signal is mainly dominated by signal scattering. As the FOV increases, the multiple scattering increases, which leads to the result that the total echo signal decays more slowly with a larger FOV.   Figure 8, we find out the dependence of multiple scattering on water turbidity. For clear water, large FOV is good for obtaining multiple scattering signals due to less multiple scattering in clear ocean. However, smaller FOV is better for the single scattering signal for coastal water because of the stronger multiple scattering. Therefore, the water inherent optical properties should be taken into account when it comes to the FOV design for a lidar system. A multiple-FOV or adjustable-FOV lidar system configuration will be a good choice for signal inversion.

Effects of SPF
To analyze the effects of different SPFs, we simulated the lidar echo signal using HG (g = 0.924), FF, TTHG, and Petzold SPFs with the same lidar system parameters as in Section 3.1. As shown in Figure 10b, FF simulation results agree with Petzold results very well as the result of a great match of their SPFs (red and purple lines in Figure 10a). However, the results of HG are slightly smaller. It is largely because HG SPF has a much smaller backscatter at a large scattering angle, as shown in Figure 10a, which is caused by the sharp forward peaks of the seawater phase function [26]. Relatively, the results of TTHG are much higher than the others due to its higher backscatter, as shown in Figure 10a. Therefore, the real oceanic SPF should be taken into account during lidar simulation.

Effects of Wind Speed
Wind exerts a subtle influence on lidar echo by changing the roughness of the sea surface. To determine effects of wind speed, we computed lidar beam reflectivity and transmissivity with different wind speeds, considering four cases when laser beam propagates through the sea surface: laser beam reflected from the air towards air (RAA); beam reflected from the water towards water (RWW); beam transmitted from air into water (TAW); beam transmitted from water into air (TWA). The results are shown in Figure 10. When the laser beam propagates from the air (Figure 11a,b), the light reflectivity changes slowly within a small laser incident angle of 30 • . Then, it increases rapidly with the increase of laser incident angle. Meanwhile, the transmissivity is, oppositely, on the decrease. The wind tends to improve the transmission of the laser. In other words, higher wind speed leads to greater transmissivity and less reflectivity. When a laser beam propagates from water (Figure 11c,d), the light reflectivity increases, while the transmissivity is on the decrease with the increase of the incident angle. When the incident angle exceeds the critical angle of total reflection, all light is reflected from water to water totally, and the critical angle increases with the increase of wind speed. The result indicates that the impact of wind-driven rough sea surface on laser transmission through the air-sea interface can be neglected in small incident angles (Figure 11b) to simplify the simulation. Although the Cox and Munk model has been widely used to simulate rough sea surface in ocean color remote sensing [17,34,35], it is obtained for wind speeds of 1-14 m/s and cannot produce wave elevations [36]. Cox-Munk simulates capillarity waves but ignores gravity waves. A more accurate sea surface model considering the shadowing effect [37] will be modified for oceanic lidar simulation in the future. As the Cox and Munk model is for isotropic sea surface, other anisotropic surfaces [38] are used to modify our model. One notes that because the scattering phase function is difficult to determine accurately, the simulation accuracy may be reduced in a complex oceanic environment. The wide-angle scattering occurs in the dense subsurface phytoplankton, which sometimes makes it hard to use this model, and the pulse stretching is included in this model.

Effects of Stratified Water
We analyzed the effects of stratified water that are simulated by using the chlorophyll Gaussian profile [39]. The vertical structure of chlorophyll concentration affects the echo signal by changing the water optical properties, including attenuation and scattering coefficients. The chlorophyll profile Chl(z) mg·m 3 is expressed by Gaussian function as follows: where C 0 is chlorophyll concentration mg·m 3 at the sea surface, C 1 is the constant background phytoplankton biomass mg·m 3 , z max is the depth where the maximum Gaussian chlorophyll value is reached, and σ c is the width of the chlorophyll peak. The parameters used here are C 0 = 1 mg·m 3 , C 1 = 0.5 mg·m 3 , z max = 10 m, σ c = 5 m. Figure 12 shows the simulation results of stratified water. Compared with results of homogeneous water, which represent a linear decay in the logarithm coordinates as shown in Figure 3b, the echo signal in stratified water has a bulge at the depth of 10 m where the maximum Gaussian chlorophyll value is reached. Therefore, the subsurface chlorophyll maximum layer can be detected by using the bulge of lidar measurements [12,40]. Besides, the multiple scattering occurs more frequently, as shown in Figure 12a. Considering the worldwide stratification in the real marine environment, stratified water is supposed to be considered in lidar simulation.  Figure 13 represents the ratio of total echo signal to single scattering with different particle radiuses. The lidar system parameters are the same as those in Section 3.1. The results show that particle sizes have a vital influence on multiple scattering of lidar. With a very small particle radius (0.1 µm), the ratio of total echo signal magnitude to single scattering is close to 1 at all depths, which means that the return signal is mainly dominated by single scattering, and multiple scattering does not increase with the propagation of the laser beam. When particle radius is greater, multiple scattering increases and plays a major role in the total signal. Considering the complexity of particulates in the ocean, the real particle radius should be taken into account during lidar simulation. Besides, this feature can be used for the retrieval of particulate sizes in water.

Conclusions
We introduced a novel, fast multiple-scattering approximate model for ocean lidar. Preliminary results indicate that it works well for oceanic lidar simulation of beam transmission in the air-sea coupled systems with a wind-driven rough sea surface. Compared with previous analytic models, our model can deal with stratified water, and it can be applied to arbitrary SPFs. However, because the scattering phase function is difficult to determine accurately in complex turbid oceanic water, the wide-angle scattering occurring in the strong scattering layer also makes it hard to use this model, and the pulse stretching is not included in this model. The results show that multiple scattering signal increases with the propagation of the laser beam in the water. The ratio of multiple scattering to single scattering becomes greater with the increase of water turbidity. However, the total signal intensity decreases with the increase of the depth because of water attenuation and decreases faster, especially in turbid water. Thus, the penetration in turbid water is rather shallow. Besides, the strong multiple scattering results in a more complex signal inversion. A good agreement between our model and the MC simulation demonstrates its validation. Their coefficient of determination is beyond 0.9, the RMSE is within 0.02, the MAD is within 0.02, and the MAPD is within 8%, which indicates that our model works well for oceanographic lidar simulation. Besides, the validation with measurement also demonstrates its efficiency in real oceanic water. Different FOVs were used to estimate their impact on multiple scattering. For a small FOV of 0.1 mrad, multiple scattering is so small that the total echo signal is mainly dominated by single scattering. For larger FOVs, the multiple scattering is more critical for the total signal. Simulation results with different SPFs show that the widely used HG SPF has poor performance at small or large scattering angles. Therefore, using the real oceanic SPF for the lidar simulation will improve its accuracy. A corresponding bulge appears in the simulation with stratified water, and it can be used for SCML detection. Different particle radiuses were applied to analyze their influence on lidar transmission. For a very small particle radius (0.1 µm), the multiple scattering is so small that it is negligible; thus, we can only consider the single scattering. For a larger particle radius, the total signal is dominated by multiple scattering.
Overall, multifactor influences including stratified water, SPF, and particular size were analyzed in our simulation, which shows that they are not negligible. These results are expected to provide a reference to the design of the oceanographic lidar system and contribute to the processing of lidar echo signals.