Semi-Analytic Monte Carlo Model for Oceanographic Lidar Systems: Lookup Table Method Used for Randomly Choosing Scattering Angles

Featured Application: A lookup table method is proposed for solving complicated or even unsolvable inverse cumulative distribution function of scattering phase function. Abstract: Monte Carlo (MC) is a signiﬁcant technique for ﬁnding the radiative transfer equation (RTE) solution. Nowadays, the Henyey-Greenstein (HG) scattering phase function (spf) has been widely used in most studies during the core procedure of randomly choosing scattering angles in oceanographic lidar MC simulations. However, the HG phase function does not work well at small or large scattering angles. Other spfs work well, e.g., Fournier-Forand phase function (FF); however, solving the cumulative distribution function (cdf) of the scattering phase function (even if possible) would result in a complicated formula. To avoid the above-mentioned problems, we present a semi-analytic MC radiative transfer model in this paper, which uses the cdf equation to build up a lookup table (LUT) of ψ vs. P Ψ ( ψ ) to determine scattering angles for various spfs (e.g., FF, Petzold measured particle phase function, and so on). Moreover, a lidar geometric model for analytically estimating the probability of photon scatter back to a remote receiver was developed; in particular, inhomogeneous layers are divided into voxels with different optical properties; therefore, it is useful for inhomogeneous water. First, the simulations between the inverse function method for HG cdf and the LUT method for FF cdf were compared. Then, multiple scattering and wind-driven sea surface condition effects were studied. Finally, we compared our simulation results with measurements of airborne lidar. The mean relative errors between simulation and measurements in inhomogeneous water are within 14% for the LUT method and within 22% for the inverse cdf (ICDF) method. The results suggest feasibility and effectiveness of our simulation model.


Introduction
Often used to solve the radiative transfer equation (RTE) of oceanographic lidar, Monte Carlo (MC) techniques refer to the algorithms that use probability theory and random numbers to simulate the fates of numerous photons propagating through sea water. The strengths of MC techniques are that they are conceptually simple, which require fewer simplifying approximations, and can be used for finding complicated 3D RTE solutions [1]. In recent years, MC simulation methods have been widely used in lidar propagation in ocean water. Peng-Wang et al. [2,3] developed a three-dimensional vector radiative transfer MC method for light propagation in atmosphere-ocean systems, which provides complete information about radiative coupling between different regions of the system. Ramella-Roman et al. [4,5] presented three MC programs of polarized light transport into scattering media; comparison between MC runs and Adding Doubling program yielded less than 1% error. Berrocal et al. [6] investigated scattering and multiple scattering of a typical laser beam in a turbid environment and obtained good agreement between experimental measurements and results from MC simulations. Gordon [7] studied multiple scattering effects by solving the RTE via MC techniques. The effective attenuation coefficient of this exponential decay is found to be strongly dependent on the parameters of the lidar system and on the optical properties of the water. Poole et al. [8] described a semi-analytic MC radiative transfer model (SALMON), which is particularly well-suited for addressing oceanographic lidar systems. It is based on the method of expected values in which an analytical estimate is made of the probability of collection by a remote receiver of scattered or emitted photons at appropriate points in the stochastically constructed underwater photon trajectory. Chen et al. [9] also developed a semi-analytic MC radiative transfer model to study laser propagation in inhomogeneous sea water within the subsurface plankton layer. Besides, there are many other MC studies for simulating photons propagating in a scattering medium [10][11][12][13].
The probability density function (pdf) for a random variable that can have values only between 0 and 1 is fundamental to MC simulation. Scattering phase function (spf) can be interpreted as a pdf: β(ψ) sin(ψ). The Henyey-Greenstein (HG) spf [14] has been widely used in many studies. The HG phase function, however, does not work well at small or large scattering angles. Though other spfs work well, e.g., Fournier-Forand phase function (FF) [15], solving cumulative distribution function (cdf) of spf P Ψ (ψ) as Equation 1 for scattering angle ψ (even if possible), namely inverse cdf (ICDF), would give a complicated formula. Besides, in almost all of these studies, they assume the water is optically vertical homogeneous; a few studies focus on simulating light propagation in inhomogeneous water. For some use cases there is a need to extend and improve existing MC methods to apply to inhomogeneous water, for instance, the ocean with subsurface chlorophyll maximum layer.
This study introduces an improved semi-analytic MC radiative transfer model based on the lookup table (LUT) method and an analytically lidar geometric model to describe photons scattering event for oceanographic lidar systems. The new model presented in this study is based on a lidar geometric model for analytically estimating the probability of photon scatter back to a remote receiver, and inhomogeneous layers are divided into voxels with different optical properties with a resolution as high as 0.11 m. The LUT method is used for solving complicated or even unsolvable inverse cumulative distribution function of various spfs. Multiple scattering and wind-driven sea surface condition effects are also studied. Finally, we compared our simulation results with measurements of airborne lidar.

Determining Photon-Free Path Length in Inhomogeneous Water
Detailed MC simulation for light propagating in sea water requires random determination of photon path length, photon scattering angle, and reflection or transmission at sea surface or bottom boundaries [16]. Most existing MC simulation methods assume that sea water is optically homogeneous, while true sea water is sometimes optically vertical inhomogeneous. It should be noted that "inhomogeneous" water in this study mainly means that vertical optical properties of the water are inhomogeneous. In order to simulate optically inhomogeneous sea water, inhomogeneous layers are divided into voxels with different optical properties. We allow α(z) and b(z) vary with depth z in this study, where α(z) and b(z) are sea water absorption and scattering coefficient at the depth of z, respectively. Therefore, a series of photon path lengths at different depths are generated as follows: where is a series of random numbers, which follow uniform distribution over the interval [0, 1], and c(z) is sea water attenuation coefficient, which is the sum of α(z) and b(z). In the case of homogeneous water, the average distance that a photon travels is 1 c .

Determining Scattering Angle: Lookup Table Method
The pdf is fundamental to MC simulation for determining the scattering angle for each photon collision. The spf β(ψ) can be interpreted as a pdf in the formula of β(ψ) sin(ψ). The HG spf has been widely used in many MC studies. However, the HG phase function does not work well at small or large scattering angles. As shown in Figure 1, the black curve is the Petzold's average spf measurement [17], which comprises the most carefully made and widely cited scattering measurements. The red curve shows the Fournier-Forand phase function with the same backscatter fraction as the measured phase function. The blue curve presents the HG phase function; we can see that it does not fit well at small angles (<20 • , especially bad for angles of less than 1 • ) or large angles (>140 • ). There is a deviation for HG spf among the curves in angles between 140 and 180 degrees because oceanic particles are quite different in their physical properties from interstellar dust, which is initially widely used. Table 1 shows the curve-fitting statistics of the two spfs compared with Petzold's measurements. Note that the HG parameter g is 0.924, which already gives the best fit to the Petzold's measurements, but the curve still does not work well.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 13 the depth of z , respectively. Therefore, a series of photon path lengths at different depths are generated as follows: where ℜ is a series of random numbers, which follow uniform distribution over the interval [0,1] , and (z) c is sea water attenuation coefficient, which is the sum of (z) α and (z) b . In the case of homogeneous water, the average distance that a photon travels is 1 c .

Determining Scattering Angle: Lookup Table Method
The pdf is fundamental to MC simulation for determining the scattering angle for each photon collision. The spf ( ) β ψ  can be interpreted as a pdf in the formula of ( )sin( ) The HG spf has been widely used in many MC studies. However, the HG phase function does not work well at small or large scattering angles. As shown in Figure 1, the black curve is the Petzold's average spf measurement [17], which comprises the most carefully made and widely cited scattering measurements. The red curve shows the Fournier-Forand phase function with the same backscatter fraction as the measured phase function. The blue curve presents the HG phase function; we can see that it does not fit well at small angles (< 20°, especially bad for angles of less than 1°) or large angles (> 140°). There is a deviation for HG spf among the curves in angles between 140 and 180 degrees because oceanic particles are quite different in their physical properties from interstellar dust, which is initially widely used. Table 1 shows the curve-fitting statistics of the two spfs compared with Petzold's measurements. Note that the HG parameter g is 0.924, which already gives the best fit to the Petzold's measurements, but the curve still does not work well.
Here, y i is the difference between spf-calculated and measured data,ŷ i is the mean value of y i , y i is the relative difference andŷ i is the mean relative difference.
Though the Fournier-Forand phase function could give the best fit to the Petzold's measurements of average particle phase function as shown in Figure 1, its cdf (Equation (3)) [18] for ψ would give a complicated formula. Solving the cdf equation of FF for ψ seems impossible. It is numerically more efficient to use Equation (4) to build up a lookup table of ψ vs. P Ψ (ψ) to determine scattering angle.
Here, n is the real part of refraction of the particles, µ is the slope parameter of the hyperbolic distribution, and δ υ 180 is δ at ψ = 180 • . In this study, we first obtain the cdfs for various spfs using the Kaplan-Meier non-parametric method [19], and then build a LUT with two columns of ψ vs. P Ψ (ψ). The Ψ angle increments in the LUT table are dense and incremental from 0.02 • to 5 • in a range between 0.1 • and 10 • ; when it reaches 10 • , the increments keep constant at 5 • . This is mainly because the spfs vary significantly in the small angles, while vary gently in the larger angles. Since the table is discrete, the resulting values are linearly interpolated between the two consecutive values. Petzold's measurements mentioned above show that linear interpolation is good enough to obtain the results between two discrete values. When a random number is generated, we find the nearest P Ψ (ψ); the corresponding angle ψ is the scattering angle we need. This is illustrated in Figure 2, which shows a Fourier-Forand cdf for ψ vs. P Ψ (ψ) as tabulated using Equation (4). The blue arrows show how drawing a value of P Ψ (ψ) = 0.694 leads to a scattering angle of ψ about 10 degree, and drawing a value of P Ψ (ψ) = 0.827 leads to a scattering angle of ψ about 20 degrees.
where P Ψ (ψ min ) and P Ψ (ψ max ) are cdf values at the minimum and maximum angles, respectively. P Ψ (ψ i ) is the cdf of the angle at row i.

Semi-Analytic MC Model
We developed an improved semi-analytic MC simulation model for analytically estimating the probability of photon scatter back to a remote receiver (i) E at collision point i, which is also appropriate for inhomogeneous water: where p( ) θ is the spf value at angle θ ; ΔΩ is the small solid angle of collision photon received by a remote detector; r A means the detector aperture; H is lidar flight height; i z is the distance from point i to sea surface, which equals the photon's depth; (i) c is sea water using the beam attenuation coefficient, and d i is the thickness of each divided layer. Because the solid angle of collision photon received by a remote detector is so small, the cosine value of the actual backscattering line is similar to the line perpendicular to the air/water interface. Therefore, each layer distance d i is vertical. The analytical lidar geometric model for describing the photon's scattering event is shown in Figure 3. Before using the equation, we should determine whether the photon's position is within the receiver field of view (FOV) first.
where FOV Ω is the receiver FOV solid angle; z A is the area at depth z i seen by the receiver; i x , i y , and i z are the photon's position in Cartesian coordinates.

Semi-Analytic MC Model
We developed an improved semi-analytic MC simulation model for analytically estimating the probability of photon scatter back to a remote receiver E(i) at collision point i, which is also appropriate for inhomogeneous water: where p(θ) is the spf value at angle θ; ∆Ω is the small solid angle of collision photon received by a remote detector; A r means the detector aperture; H is lidar flight height; z i is the distance from point i to sea surface, which equals the photon's depth; c(i) is sea water using the beam attenuation coefficient, and d i is the thickness of each divided layer. Because the solid angle of collision photon received by a remote detector is so small, the cosine value of the actual backscattering line is similar to the line perpendicular to the air/water interface. Therefore, each layer distance d i is vertical. The analytical lidar geometric model for describing the photon's scattering event is shown in Figure 3. Before using the equation, we should determine whether the photon's position is within the receiver field of view (FOV) first.
where Ω FOV is the receiver FOV solid angle; A z is the area at depth z i seen by the receiver; x i , y i , and z i are the photon's position in Cartesian coordinates. At each photon scattering event, a photon is scattered into a new trajectory and the scattering angle θ is determined by the LUT method described in the previous section. Meanwhile, the photon's weight w(i) is multiplied by the single scattering albedo w(i). If a photon's weight drops under the preset threshold of 10 -4 , a roulette procedure is triggered and employed. In this study, divided layer vertical resolution can reach up to 0.11 m.
At each photon scattering event, a photon is scattered into a new trajectory and the scattering angle θ is determined by the LUT method described in the previous section. Meanwhile, the photon's weight w(i) is multiplied by the single scattering albedo w(i) . If a photon's weight drops under the preset threshold of 10 -4 , a roulette procedure is triggered and employed. In this study, divided layer vertical resolution can reach up to 0.11 m.

Boundary Consideration
As a photon reaches the sea surface, an estimation for total transmission t t of ocean interface due to specific viewing angle and sea wave conditions of rough sea surface [20,21] is given: Actual measurements of bidirectional reflectance distribution functions (BRDFs) of ocean bottom materials like sand or sea grass canopies have rarely been made. Because of the lack of measurements and models of the BRDF for actual ocean bottom materials, it is usually assumed that a bottom is a Lambertian reflecting surface [22]. In this study, the sea bottom is simply treated as a Lambertian surface, and for that reason lidar bathymetry simulation is not useful for our purposes.

Boundary Consideration
As a photon reaches the sea surface, an estimation for total transmission t t of ocean interface due to specific viewing angle and sea wave conditions of rough sea surface [20,21] is given: where t p is means the part of transmission under calm sea surface calculated with Snell's law; and t r is the part of transmission in a wind-driven roughened surface and can be obtained as a function of wind speed u: Actual measurements of bidirectional reflectance distribution functions (BRDFs) of ocean bottom materials like sand or sea grass canopies have rarely been made. Because of the lack of measurements and models of the BRDF for actual ocean bottom materials, it is usually assumed that a bottom is a Lambertian reflecting surface [22]. In this study, the sea bottom is simply treated as a Lambertian surface, and for that reason lidar bathymetry simulation is not useful for our purposes. Its reflectivity is the same in all directions. A reflection event is forced by multiplying the photon's weight by sea bottom albedo.

Simulation Comparison between Inverse Function for HG cdf and LUT Method for FF cdf
The inverse equation of HG cdf for scatter angle ψ is Figure 4 shows the simulation comparison between the inverse cdf method for HG and the LUT method for FF. Simulations were conducted using a total of 10,000,000 photons for a given detector height of H = 300 m, FOV = 50 mrad, detector aperture of A = 0.06 m 2 , sea bottom depth of d = 40 m, vertical layer resolution of 0.11 m, water absorption coefficient of α = 0.04 m −1 , scattering coefficient of b = 0.06 m −1 , and sea bottom albedo of ρ b = 0.03. The given n and µ for FF spf were 1.10 and 3.62, respectively. The HG parameter was set to 0.924. It took 129 s and 140 s for the LUT method and ICDF method simulations, respectively. It shows that the simulation curve using the LUT method of FF is larger and more stable than using the ICDF method of HG under the same input condition. It may be because the forward scattering of HG spf is much less than that of the Petzold's measurements (see Figure 1). Besides, the back scattering is also less; therefore, the p(θ) value in Equation (2) is less. We can see that the black curve by LUT method looks smoother than the red curve by the inverse function method. The adjusted R-square (N = 498) between the two curves in Figure 1 is 0.948. The residual sum of squares is 3.2 × 10 −8 . The RMSE is 8.03 × 10 −6 . The min value, max value, mean value, and standard deviation of the red curve are 1.54 × 0 −8 , 7.70 × 0 −4 , 9.82 × 0 −5 , and 1.73 × 0 −4 , respectively. The min value, max value, mean value, and standard deviation of the black curve are 1.11 × 10 −9 , 1.63 × 10 −4 , 2.48 × 10 −5 , and 3.53 × 10 −5 , respectively (see Table 2). The signal intensity by LUT method appears higher than ICDF method: This is because the ICDF method to calculate the backscattering probability return to receiver was to obtain a possible certain backscattering angle, which may be less than the small receive angle range, while the LUT method to calculate the backscattering probability return to receiver was to obtain the integrals between two discrete values, which is a little larger than the small receive angle range. The standard deviation of the red curve by ICDF method is more than 5 times larger than the black curve by LUT method. It looks like there is more noise and disturbance in the red curve by ICDF method under the same preset parameters. The results also indicate the LUT method is more stable than the inverse function method. We also find that the sea bottom curve based on the LUT method for FF cdf has less pulse widening effect than the inverse function method for HG cdf.    Figure 1 by the two methods.   Figure 1 by the two methods.

Multiple Scattering and Wind-Driven Roughened Sea Surface Effect
Influences of multiple scattering and wind-driven roughened sea surface conditions on simulation results were then analyzed. We separated the first order of scattering, second order of scattering, third order of scattering, and more than third order of scattering parts in the total scattering signal, as shown in Figure 5a. The black line is the total signal, which is the sum of the single, double, triple and more than triple scattering event curves. Figure 5a shows that the single scattering signal decreases as depth increases, while multiple scattering signal increases first and then decreases. This is because, as the photon travels forward, multiple scattering occurs more frequently in water. Besides, the signal is the deepest depth the photon attained prior to scattering back to the surface and being detected. It also shows that the single, double, triple and greater than triple scattering events curves cross at about 20 m. The multiple scattering occurs more frequently as the water scattering coefficient becomes larger. When the water scattering coefficient is large, the depth of curves cross decreases, and the depth increases when the water scattering coefficient is small. The first order of scattering in logarithmic form decreases linearly with the increase of depth, which agrees with Mobley's conclusion that light falls exponentially in water [1], and it indicates that our simulation method is correct and feasible. The contribution of high-order scattering signal increases with depth, until it dominates at a certain depth. Figure 5b shows the slope of return signal, namely, lidar effective attenuation coefficient (K sys ) falls between beam attenuation coefficient c and absorption coefficient α. The system attenuation coefficient, K sys , lies somewhere between α and c, and is dependent on the relative values of the transmitter beam spot size, transmitter divergence, receiver aperture, and receiver field of view. As Figure 5b shows, we can see that the cyan curve using a narrow FOV fits well with the blue curve of c and the black curve using a large FOV fits well with the dark green curve of α. It approaches c with a narrow FOV and approaches α with a large FOV. K sys approaches with a narrow field of view c since almost all of the light that had been scattered was not detected by the receiver; thus, the decay rate of the backscatter signal is the sum of α and b. Since the field of view was large enough, K sys approaches α since the only light lost at the receiver was what had been absorbed. However, for a very large beam and large receiver aperture and field of view, K sys approaches α + b b since almost all of the forward propagating photons remain in the receiver FOV until absorbed or scattered in the backward direction. However, when the field of view was large enough, K sys approaches α, since the only light lost at the receiver was what had been absorbed. These results agree with previous studies [23][24][25]. They indicate that the configuration of multiple FOV applied for lidar system is useful for optical properties' inversion in various types of water. Figure 5c,d show the influence of wind-driven roughened sea surface on simulation signals. The incident angle in Figure 5c,d is the lidar incident angle. The reason why it spans to about 90 degrees in the "From air to sea" plot, but only to about 48 degrees in the "From sea to air" plot, is due to the critical angle for total internal reflection. We can see that the lidar transmittance coefficients remain unchanged when wind speed is less than 9 m/s, but it decreases significantly when the wind speed is greater than 9 m/s. It indicates that both wind-driven roughened sea surface and multiple scattering have significant effects on simulation results, which should not be ignored.
to about 90 degrees in the "From air to sea" plot, but only to about 48 degrees in the "From sea to air" plot, is due to the critical angle for total internal reflection. We can see that the lidar transmittance coefficients remain unchanged when wind speed is less than 9 m/s, but it decreases significantly when the wind speed is greater than 9 m/s. It indicates that both wind-driven roughened sea surface and multiple scattering have significant effects on simulation results, which should not be ignored.

Simulated-Signal vs. Airborne Lidar Measured-Signal
Subsequently, simulation results (from both LUT and ICDF methods) were compared with airborne lidar measurements. The experimental airborne oceanic lidar was developed by the Shanghai Institute of Optics and Fine Mechanics (SIOM), Chinese Academy of Sciences (Figure 6a). The transmitting system utilized a 532 nm frequency-doubled Nd:YAG laser. Its pulse duration is 1.5 ns, and pulse repetition rate (PRF) is 1 kHz. The detector telescope's diameter is 200 mm and FOV is 40 mrad. The layer thickness employed in our simulation is 0.11 m; the phase function lookup table was kept unchanged during a simulation process. The profile optical properties of sea water

Simulated-Signal vs. Airborne Lidar Measured-Signal
Subsequently, simulation results (from both LUT and ICDF methods) were compared with airborne lidar measurements. The experimental airborne oceanic lidar was developed by the Shanghai Institute of Optics and Fine Mechanics (SIOM), Chinese Academy of Sciences (Figure 6a). The transmitting system utilized a 532 nm frequency-doubled Nd:YAG laser. Its pulse duration is 1.5 ns, and pulse repetition rate (PRF) is 1 kHz. The detector telescope's diameter is 200 mm and FOV is 40 mrad. The layer thickness employed in our simulation is 0.11 m; the phase function lookup table was kept unchanged during a simulation process. The profile optical properties of sea water measured at the station (109 • 43.808 E, 18 • 2.205 N) in Sanya Bay in the South China Sea on 12 March 2018 were input into the model, as shown in Figure 6b. The green line is the chlorophyll a concentration profile, and the black line represents water attenuation coefficient profile. We can see that a subsurface chlorophyll maximum (SCM) layer exists at about 30 m. Clearly, most current MC simulation methods based on the assumption of optically vertical homogeneous are no longer suitable. Figure 6c shows the comparison between simulation signal based on our method with airborne lidar measured-signal. The simulations were conducted using a total of 10,000,000 photons and actual lidar detector height is H = 307 m. We find that the simulated curve using the LUT method (blue line) agrees well with the lidar measured-signal (red line). The ICDF method appears to overestimate the actual signal, which is because the ICDF method used to calculate the backscattering probability return to receiver obtained a possible certain backscattering angle, which is less than small receive angle range. The LUT method appears to underestimate the actual signal at greater depths and also at shallower depths, crossing at about 20 m. This is because the LUT method used to calculate the backscattering probability return to receiver obtained the integrals between two discrete values, which is a little larger than small receive angle range at deeper depths. It may be better to set the LUT increments according to the solid receiver angle. Figure 6d shows the relative errors between simulation and measurements for the LUT method within 14% mean relative error and for the ICDF method within 22% mean relative error. We can see that the relative error was less than 10% for the LUT method and less than 20% for the ICDF method within 20 m. Both methods have the maximum relative errors when the depth is between 25 and 30 m, where there was a chlorophyll maximum scattering layer. It may be due to the fact that the spfs in our study are no longer applicable when reaching the chlorophyll maximum scattering layer. Figure 7 shows the regression plots of simulated-signal vs. airborne lidar measured-signal. Figure 7a,b presents the regression plots of ICDF and LUT retrieved data vs. measured data, respectively. Figure 7c,d presents the residual plots of ICDF and LUT, respectively. It shows that the adjusted R-square between simulation and measurements (N = 347) for LUT method (R 2 = 0.996) is larger than that by ICDF method (R 2 = 0.996). The data retrieved from both methods have significant correlations with measured data, although the RMSE by LUT method (RMSE = 0.016) is relatively smaller than that by ICDF method (RMSE = 0.052). As Figure 7c,d shows, the LUT method has relatively smaller residual than that by ICDF method in almost all data. In general, the results indicate that the LUT method is effective for solving complicated or even unsolvable inverse cumulative distribution function of various scattering phase functions, and it significantly improves modeling accuracy.

Conclusions
In this study, we introduced an improved semi-analytic MC radiative transfer model based on the LUT method and a lidar geometric model for scattering event, which is also applicable to inhomogeneous water. Comparison of simulation and measurements indicates that our method is effective for solving complicated or even unsolvable inverse cumulative distribution function of various scattering phase functions and for significantly improving modeling accuracy. This model also considers multiple scattering, wind-driven roughened sea surface and water optical inhomogeneous effects. Thus, fundamental radiative transfer problems in oceanographic lidar and effects of a host of environmental factors on lidar system performance can be addressed. We plan to apply our method to in-situ measured spfs, especially subsurface phytoplankton scattering spfs. In addition, simulation results will be used to study the relationships between lidar effective coefficients and water bio-optical properties for quantificational applications. More comparisons between measurements and simulations will be carried out to improve our method.