Open Access
Add to BibSonomy
Abstract
A continuous field Monte Carlo radiative transfer model with an improved semianalytic approach is developed to study laser propagation in an inhomogeneous dust environment. In the proposed model, the photon step size can vary with the mass concentration of the dust environment. Additionally, the scattering properties of the dust particles are calculated with the T-matrix method and the T-matrix scattering phase function is applied to the Monte Carlo simulation with a rejection method. Using this model, the influences of the particle sizes and shapes on the backscattering properties are studied. Finally, the laser echoes simulated by our proposed model are compared with those of traditional Monte Carlo method and experimental results. Different mass concentration distributions indeed influence the simulated laser echo. The simulated results (of our proposed model) agree well with the measured data, demonstrating the effectiveness and accuracy of our approach for inhomogeneous media.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Light propagation in inhomogeneous media is a complicated phenomenon. The multiple scattering of light causes dramatic alterations in the intensity, propagation direction, polarization, and coherence of light [1].
Because light transport through scattering media is difficult to solve due to the disadvantages of analytical solutions [2–9] for the radiative transfer of light within a bounded medium [10], numerical solutions of the radiative transfer equations and Monte Carlo (MC) method [11–13] have been widely used and applied to the optical scattering of scattering media [1]. Compared with the numerical solutions, the MC method requires fewer simplified approximations, and it can be extended to more complex media with relative ease to study laser radiative transfer [14]. In the MC method, the scattering is described by the scattering coefficient and the scattering angle is typically drawn from the inverse cumulative distribution of the phase function. The Henyey-Greenstein (H-G) scattering phase function [15] is the most commonly used since it offers an analytical inverse of the cumulative distribution function (CDF), which is convenient for sampling the scattering angles with a random number drawn from a uniform distribution. However, the H-G phase function is well known to underestimate large-angle backward scattering observed in turbid media [16]. Consequently, other phase functions have been investigated, such as the modified Henyey-Greenstein (MHG) [17,18], Gegenbauer kernel (G-K) [19,20], Mie [21], and Mie fractal [22] phase functions. However, the phase functions mentioned above mainly describe spherical particles, which are usually based upon a perceived lack of practical alternatives in the MC method [23]. In fact, the scattering properties of non-spherical particles can significantly differ from those of volume-equivalent or surface-equivalent spheres. In recent years, the electromagnetic scattering properties of non-spherical particles have been of great interest in a broad range of applications. At present, the T-matrix method is a powerful exact technique for computing light scattering with non-spherical particles based on numerically solving Maxwell’s equations [24]. Unlike the H-G phase function, the T-matrix method does not offer the analytical inverse of the CDF. In order to use this type of phase function in the MC simulation, a numerical sample method has been proposed [25,26]. Lu Bai et al. [27] used the T-matrix method to calculate the scattering phase function of non-spherical particles with a single particle size and they applied the sample scattering phase function method to some typical non-spherical particles in the UV band. However, they did not further use the results in the radiative propagation-related studies. Many researchers have used T-matrix method to model and analyze the optical properties of complex particles [28]. However, no one has used it in combination with the MC method for echo research, as far as the authors know.
Additionally, in the MC simulations for scattering media, such as smoke and dust, the scattering coefficient and the extinction coefficient are always assumed to be constant throughout the motion of the photons [29–32]. However, in an actual situation, homogeneous media rarely exist. An MC model of light transport in multi-layered media (MCML) has been developed [33] by dividing inhomogeneous media into several layers with various optical properties, and this model is now widely used for ocean detection, photomedicine and photobiology [34]. The optical properties of homogeneous media such as smoke, dust, and soot environments are always changing continuously, and the media cannot be divided simply. Bohu Liu et al. [35] used the dynamic changes of the smoke aerosol scattering coefficient and extinction coefficient with the smoke concentrations based on the Gaussian plume model in the MC simulation. However, they still adopted the H-G function to calculate the scattering angles of sphere particles. Furthermore, despite the changing extinction coefficient, the photon step size was based on the fixed extinction coefficient at the current position, ignoring the parameter changes on the path. Additionally, the consistency of the plume model and the actual smoke situation was not verified and the precision of the model was unknown.
The goal of this study is to develop an improved Monte Carlo radiative transfer model of laser propagation based on the T-matrix method in a continuous concentration field. In Section 2, the proposed MC model is presented in detail. The new step size method for polydisperse particles with the extinction coefficient changing with the mass concentration is introduced first. Then, based on the semianalytic approach, the corresponding emission and reception model is established. The estimated value for the continuous field is given. Finally, the sample method based on the T-matrix phase function is described. In Section 3, we discuss the influences of different particle sizes and shapes on the scattering properties, and then analyze the effects of different particle sizes and shapes on the laser echo with the sample method based on the scattering properties obtained with the T-matrix method. Section 4 describes how a laboratory is designed and established to carry out mass concentration measurements and laser echo experiments in a dust environment for the purpose of validation of the continuous field Monte Carlo radiative transfer model. Then we make a comparison of the simulated results with the measured data, and excellent agreement is observed between the proposed MC model and the experimental results. The effects of different mass concentration distributions on the laser echo calculated with the proposed MC method are also analyzed. The differences between the results of different mass concentration distributions are obvious, illustrating the necessity of the proposed model.
2. Theory and modeling
When light encounters particles in a dust environment, various physical phenomena occur, such as refraction, absorption, reflection, and diffraction, as schematically shown in Fig. 1. Interaction with a particle may change the direction in which a photon travels, which is collectively known as scattering. In an inhomogeneous medium, more attention should be paid to the influence of multiple scattering. Therefore, we use the MC method to model the scattering process of a laser system in a dust environment, which converts the light propagation problem to the interaction between photons and randomly distributed particles. The strength of the scattering strongly depends on the particle sizes and shapes of the light for a certain wavelength.
2.1 Photon step sizes in a dust environment
The step size $s(0 \le s < \infty )$ of the photon packet is calculated based on a sampling of the probability distribution for the photon’s free path [36]. We first consider an infinite inhomogeneous medium. According to the definition of the absorption coefficient ${\mu _\textrm{a}}$ and the scattering coefficient ${\mu _\textrm{s}}$, the probability of interaction with the medium per unit path length in the interval $(s^{\prime},s^{\prime} + \textrm{d}s^{\prime})$ is
In multi-layered media, the photon packet may experience free flight over changing extinction coefficients before an interaction occurs. In this case, Eq. (2) becomes [33]
In an actual dust environment, the mass concentration usually changes with the position in the atmosphere. As a result, a photon may experience free flight over dust aerosol with different mass concentrations before an interaction occurs, as shown in Fig. 2.
In this research, the continuous mass concentration is measured directly; that is, the length of the photon can be divided into countless layers ($i \to \infty$), and then $\Delta {s_i} \to 0$, i.e. $\Delta {s_i} = \textrm{d}s$. In this case, the counterpart of Eq. (3) becomes
Where $P\{ s \ge 0\} = 1$ and ${\mu _\textrm{t}}(s)$ is the extinction coefficient for the current position. Therefore, Eq. (5) can be rearranged as
Then the probability density function of the free path ${s_{\textrm{sum}}}$ is
In the MC method, the generally non-uniform probability density function $p(\chi )$ that defines the distribution of $\chi$ over the interval $(a,b)$ can be sampled by solving the following equation [33]
where $\xi$ is a random variable generated by a computer, which is uniformly distributed over the interval (0, 1). $p({s_{\textrm{sum}}})$ can be substituted into Eq. (8) to yieldDue to the symmetry of $\xi$, we can get the sampling equation
If $U(s)$ is the primitive function of ${\mu _\textrm{t}}(s)$, then Eq. (10) will become
where $U({s_{\textrm{sum}}})$ and $U(0)$ are dependent on the corresponding coordinate position. When the extinction coefficient ${\mu _\textrm{t}}(s)$ and the random number $\xi$ are given, the step size ${s_{\textrm{sum}}}$ can be determined by solving Eq. (11). An extinction coefficient can be calculated by the relationship [37] where $N(r)$ is the particle size distribution, ${C_{\textrm{ext}}}(r)$ is the extinction cross section, and $[{r_{\min }},{r_{\max }}]$ is the particle size range. For aerosols composed of particles of the same size, it is easy to relate the different methods of characterizing the concentration. For polydisperse aerosols, there is no simple relationship between the mass and number concentrations. Hence, the particle size distribution function is normalized to unity as followsTherefore, $n(r)$ can be assumed to remain unchanged per unit volume of gas. Then for arbitrary region
where $K(x,y,z)$ represents the scale factor in a certain position. The total volume of particles in a volume of gas is [38]In this way, the mass concentration $\rho (x,y,z)$ can be written as
Therefore, the extinction coefficient in Eq. (12) is related to the mass concentration by
For a photon located at $({x_0},{y_0},{z_0})$ travelling for the step size ${s_{\textrm{sum}}}$ in the direction $(u,v,w)$, the new coordinates $(x,y,z)$ are given by [35]
Substituting Eq. (19) into Eq. (18) can produce the relationship between the extinction coefficient and the photon step size, that is, ${\mu _\textrm{t}}(s)$, and the step size ${s_{\textrm{sum}}}$ can then be solved for. Therefore, the relationship between the step sizes and the mass concentration can finally be determined in order to analyze multiple scattering in an inhomogeneous dust environment.
2.2 Emission and reception model
The emission and reception model is shown in Fig. 3. In our MC program, a Cartesian coordinate system is established with the receiving center of the detector as the origin in order to facilitate the calculation, and the axes of the field view of the emitter and the detector are parallel. Hence, the displacement of the optical axis is ${d_0}$ and the beam emission direction is along the z axis, where ${\varphi _\textrm{e}}$ and ${\varphi _\textrm{r}}$ are the transmitting/receiving half angles of the optical field, A is the receiving optics aperture area, and ${t_\textrm{p}}$ is the pulse width of the laser.
In this research, we use an improved semianalytic approach originally proposed by Poople et al. [41]. For multi-layered media with various optical properties, the estimated value ${E_m}$ of the fraction of photons collected by the receiver upon a scattering event at the collision point is given by [42]
While in inhomogeneous media, the extinction coefficient is constantly changing with the distance ${R_m}$. Hence, according to Eq. (5), the estimated value in this research can be written as
After a semianalytic sensing procedure, the photon continues to transmit along the original orientation, and the weight becomes [43]
When the photon’s weight drops below the threshold ${\omega _\textrm{e}}$ (${\omega ^{\prime}_m} < {\omega _\textrm{e}}$) after one semianalytic reception, the photon is terminated. The detected photon delay time ${t_\textrm{p}}$ is generated by three parts, including the launch delay, scattering delay, and receiving delay, which is expressed by
where ${t_\textrm{e}}$ is the launch delay, ${s_\textrm{P}}$ is the total path length in the dust environment, and c is the velocity of light. The simulated echo signal can be recovered by counting all of the times of the detected photons.2.3 Sample method of scattering angles
The choice of scattering angles is a fundamental problem in the MC method, which is based on the scattering phase function and a rejection method. Therefore, we use the discrete data for the scattering phase function obtained with the T-matrix method and a rejection method. The best advantage is that we can get a more precise distribution of scattering angles even for polydisperse non-spherical particles.
The sample method may be described as follows:
- (1) An incremental table of the discrete data of the scattering phase function ${a_1}(\theta )$ versus the scattering angle is created.
- (2) These phase function are unitary, that is [27] where n is the amount of discrete data.
- (3) A random number $\xi \in (0,1)$ is generated.
- (4) Based on the rejection method, the unknown $m$ ($m \in (1,n)$) that can be used to obtain the minimum value of $|{{P_m}(\theta ) - \xi } |$ is found. Then the new scattering angle ${\theta _m}$ and the corresponding scattering phase ${a_1}({\theta _m})$ can be determined.
3. Simulation results and discussion
3.1 Influence of particle sizes and shapes on scattering properties
In the real dust environment, the particle size is variable and cannot easily be expressed by the single particle size hypothesis. Fortunately, many plausible size distributions can be adequately represented by just two parameters, i.e., the effective radius ${r_{\textrm{eff}}}$ and the effective variance ${v_{\textrm{eff}}}$, defined by [39]
To analyze the influence ${r_{\textrm{eff}}}$ and ${v_{\textrm{eff}}}$ have on the size distribution, Eq. (26–28) are used to obtain the corresponding ${\sigma _r}$ and $\bar{r}$ in Eq. (29). Then the certain size distribution is confirmed, as shown in Fig. 4. Figure 4(a) demonstrates the broadening of the size distribution with the increasing ${v_{\textrm{eff}}}$ while the effective radius ${r_{\textrm{eff}}}$ stays constant. The effective radius mainly influences the proportion of particle size. With the increase of ${r_{\textrm{eff}}}$ and constant ${v_{\textrm{eff}}}$, the number of particles with a larger size also increases, as shown in Fig. 4(b).
Fig. 4. Log-normal particle size distributions with different (a) ${v_{\textrm{eff}}}$ and (b)${r_{\textrm{eff}}}$.
Similarly, the scattering phase functions are calculated using the varying effective radius and effective variance with the T-matrix method. The particles in the calculation are oblate spheroids with aspect ratios of $\varepsilon = 2$ ($\varepsilon$ is the ratio of the horizontal semi-axis to the rotational semi-axis), as shown in Fig. 5. The effective variance ${v_{\textrm{eff}}}$ has little influence on the forward scattering region (<80°) and the increase of ${v_{\textrm{eff}}}$ will reduce the backscattering (see Fig. 5(a)). However, it is difficult to summarize the effect of effective radius ${r_{\textrm{eff}}}$ from Fig. 5(b).
Fig. 5. Scattering phase functions with different (a) ${v_{\textrm{eff}}}$ and (b) ${r_{\textrm{eff}}}$.
Additionally, the shapes of particles influence the scattering characteristics. We mainly study the effects of spheroids, which are formed by rotating an ellipse about the minor (oblate spheroid) or major (prolate spheroid) axes (Fig. 6). The shape and the size of a spheroid can be conveniently specified by the axial ratio $\varepsilon$ and the radius of a sphere having the same volume, ${r_v}$. The $\varepsilon$ is greater than 1 for oblate spheroids, smaller than 1 for prolate spheroids, and equal to 1 for spheres.
As shown in Fig. 7(a)-(c), the scattering phase functions and the scattering intensities of polydisperse spheroid particles and corresponding monodisperse spheroid particles with the equal-volume-radius are calculated with the T-matrix method. The H-G method is also used to calculate the scattering phase function of the monodisperse particles to visually demonstrate the limitations of the H-G phase function. Furthermore, the scattering phase functions of polydisperse spheroid particles with different $\varepsilon$ are shown in Fig. 7(d). The other parameters in the calculation are as follows. The incident wavelength is 905 nm, and the equal-volume-radius of the spheroids and the effective radius of the polydisperse particles are both 1.1 µm. The effective variance ${v_{\textrm{eff}}}$ of the polydisperse particles is 0.1, and the complex refractive index is 1.52 + 0.008i [45].
Fig. 7. (a) Scattering phase function and scattering intensity of sphere particles; (b) Scattering phase function and scattering intensity of oblate particles; (c) Scattering phase function and scattering intensity of prolate particles; (d) Scattering phase function with different $\varepsilon$.
From Fig. 7(a)-(c), the obvious difference between the scattering phase functions calculated with the T-matrix method and the H-G phase function can be observed. Additionally, compared with the monodisperse sphere, the scattering phase function of the polydisperse sphere is smoother (Fig. 7(a)), as a result of the interference and resonance effects for particles of a single size [46]. However, these effects are inconspicuous for the spheroids in Fig. 7(b)-(c). Moreover, the scattering phase functions of oblate and prolate spheroids significantly differ from those of volume equivalent spheres, as shown in Fig. 7(d). Because of the dominance of the diffraction contribution to the nearly direct forward scattering (< 20°), it is the region least sensitive to particle asphericity. The latter (> 20°) is determined by the average area of the particle geometrical cross section, which is determined by the particle shape. In particular, for the nearly direct backscattering (> 150°), which has a great influence on the laser echo, the degree of sphericity plays a decisive role. In other words, the closer the spheroid particle is to a spherical shape, the greater the backscattering is. Furthermore, the scattering phase functions of oblate and prolate particles whose $\varepsilon$ are the inverse of each other have similar trends.
To further analyze the effects of the sizes and shapes on the backscattering property, the average backscattering cross section per particle $\left\langle {{C_b}} \right\rangle$ and the backscattering efficiency factor ${Q_b}$ [40] of the above particles are calculated as
Table 1. Backscattering properties
Similar to the direct backscattering (${a_1}({180^ \circ })$) shown in Fig. 5(a) and Fig. 7, the average backscattering cross section per particle $\left\langle {{C_\textrm{b}}} \right\rangle$ and backscattering efficiency factor ${Q_b}$ both decrease with the increasing ${v_{\textrm{eff}}}$ and the increasing asphericity. However, the opposite trend appears in these two parameters with the changing ${r_{\textrm{eff}}}$ when the ${v_{\textrm{eff}}}$ is constant. This variability makes the analysis of backscattering intensity difficult. Therefore, close attention should be paid to this issue in future research.
3.2 Simulation signal echo with different particle parameters
The relevant simulation parameters are listed in Table 2. These parameters are consistent with the experimental system settings.
Table 2. Simulation parameters
3.2.1 Effect of particle sizes on simulated signal echo
For a given radius range and mass concentration, we only change the effective variance ${v_{\textrm{eff}}}$ and the effective radius ${r_{\textrm{eff}}}$ of the particles (oblate with aspect ratio 2) in the MC program. Additionally, the simulation results are normalized according to the result of oblate particles with ${r_{\textrm{eff}}} = 1.1$µm, ${v_{\textrm{eff}}} = 0.1$ for intuitive comparison, as shown in Fig. 8(a)-(b). Similar to Fig. 5(a), the increase of ${v_{\textrm{eff}}}$ will reduce the backscattering and then decrease the signal echo. Separately, when ${v_{\textrm{eff}}} > 0.2$, there is no longer a remarkable difference between the signal echoes, as can be seen from Fig. 8(a). However, unlike the irregular conclusion of Fig. 5(b), with the increase of ${r_{\textrm{eff}}}$, the signal echo is significantly reduced, as shown in the trend of ${Q_b}$ in Table 1. Therefore, we can summarize that the effective radius ${r_{\textrm{eff}}}$ of dust particles has a greater effect than the effective variance ${v_{\textrm{eff}}}$ on the backscattering. Therefore, this indicates that the smaller particle size and the more compact particle size range for the same type of dust environment will be more likely to return an echo, obscuring the target signal in a laser detection system.
Fig. 8. Normalized simulation signal echo with varying (a) ${v_{\textrm{eff}}}$ and (b) ${r_{\textrm{eff}}}$.
3.2.2 Effect of particle shapes on simulated signal echo
The influence of the particle shapes on simulated results is then analyzed. Simulations are conducted for a given size distribution, and the other parameters are the same as those given in Table 2. The normalized simulation signal echoes according to the oblate particle with $\varepsilon \textrm{ = }2$ are shown in Fig. 9. The conclusion is exactly the same as the scattering phase functions (> 150°) of Fig. 7(d), which determines the backscattering intensity to some degree. By comparing the relevant results (Fig. 7(d) and Table 1), it can be found that for different spheroid shapes, the signal echo is also affected greatly by the degree of sphericity. However, the comparison between the oblate spheroid and the prolate spheroid with the reciprocal $\varepsilon$ is difficult to summarize.
4. Experiment and validation
4.1 Dust environment laboratory
For a relatively stable dust environment, a dust environment laboratory with three functional zones: control room, test room, and settling chamber (Fig. 10) is designed and established to validate the MC model in this research.
We create the required dust environment (mainly NH4Cl) in the test area, and we observe and record in the control room. Furthermore, we can adjust the length of the test area as needed. The test room is a long narrow cube with a length of $L = 12\textrm{m}$, a width of $W = \textrm{2m}$, a height of $H = \textrm{2m}$ and ${R_z} = 6\textrm{m}$ in this case (Fig. 11).
4.2 Mass concentration measurement
4.2.1 Experiment set up and procedure
The measurement of the mass concentration of the dust environment in the laboratory is implemented using the real-time dust monitor (CEL-712 Microdust Pro from CASELLA) whose extensive range is 0.001 mg/m3 to 250 g/m3. Additionally, there are three fixed concentration measuring instruments MODEL 2030 (named as CM1, CM2, and CM3) monitoring the change of the mass concentration along the y axis all of the time.
The dust is produced incessantly for about five minutes in all of the experiments. We perform the mass concentration measurements four times with the same conditions (namely, temperature, humidity, and windless room). For convenience, we number and name the four measurements Dust 1, Dust 2, Dust 3, and Dust 4 in this work. In the first three measurements (Dust 1, Dust 2, and Dust 3), the Microdust Pro moves steadily along the light path (z axis) several times at specific moments after the man-made dust environment is relatively stable. While in Dust 4, the Microdust Pro is fixed at point A (z = 0) as shown in Fig. 12, and the laser echo power is recorded at the same time. In the measurement, the Microdust Pro records the data every second, and the CMs record every three seconds.
Fig. 12. Experiment set up of mass concentration measurement (a) inside the laboratory and (b) outside the laboratory.
4.2.2 Results of mass concentration measurement
The raw data of the mass concentration measurement are shown below. It can be found from Fig. 13(a)-(d) that the dust environments become stable in the enclosed lab after 900 s, then decrease smoothly because of sedimentation and adsorption. Due to the enclosed space and the same conditions, the dust environments in the experiments display a good consistency despite the complexity and the uncertainty of the two phase flow.
Fig. 13. Mass concentration data measured by (a) Microdust Pro along the z-axis and (b) CM1 and (c) CM2 and (d) CM3.
To investigate the agreement of the dust environments, we calculate the relative errors (> 1000 s) of each concentration measuring equipment in different dust environments. The results and the maximum errors are shown in Fig. 14.
The average errors of the equipment are all less than 10%, as listed in Table 3. It can be seen that the relative errors for Microdust Pro are relatively larger than others because of the slight flow caused by the equipment movement. In general, the average errors are acceptable for the dust environment. Therefore, it can be considered that we obtain consistent dust environments in this laboratory for the same conditions.
Table 3. Average relative errors for mass concentration
As can be seen from Fig. 13, the mass concentrations keep dropping after the dust environments become stable. Therefore, it is difficult to measure the spatial mass concentration at one moment. However, we find that the downtrends of the measured mass concentrations in this work are similar (Fig. 15). Additionally, the decrement can be ignored in a short time. For instance, the moving time for Microdust Pro once in the measurement procedure is 75 s, the corresponding reduction is within 5% according to Fig. 15. With the decreasing concentration the reduction value will also go down. As a result, the mass concentration variation during the movement of Microdust Pro is ignored.
Fig. 15. Normalized mass concentration measurement data in (a) Dust 1 and (b) Dust 2 and (c) Dust 3 and (d) Dust 4.
To further analyze the law of dust concentration distribution in the test room, we process and compare the data measured with the moving Microdust Pro in Dust 1 to 3, as illustrated in Fig. 16. The results shown in Fig. 16 confirm our previous conclusion that dust environments in this enclosed laboratory for the same conditions display a good consistency. Furthermore, it can be found that when the mass concentrations are in the range of 200 mg/m3 to 400 mg/m3, the mass concentrations go up first and then go down along the z axis, and the trends of the mass concentration in three moments are similar through translation along the vertical axis. Clearly, the existing MC method based on the assumption of optically homogeneous media is no longer applicable in this case.
Since the change trends of the curves are not complicated, we use polynomial fitting to make better use of Eq. (18), that is
where $\rho (0)$ is the mass concentration of Point A. The dust concentrations near the center line (z-axis in this paper) are quite similar for the long narrow room [47] and for small particles (< 4.5 µm) the influence of gravity is negligible [48] in the enclosed space. Therefore, we assume that the mass concentrations across the cross section (x-y plane) are the same in the simulation below. Thus, Eq. (11) becomesThen the photon step size in a continuous field can be acquired through the measured mass concentrations along the photon path.
4.3 Comparison and validation
4.3.1 Laser echo experiment results
In Dust 4, we simultaneously conduct a laser echo experiment (Fig. 17). In this experiment, the laser system employs a Nanostack Pulsed Laser Diode at 905 nm (SPL PL90_3 from OSRAM) and a 5 mm Silicon PIN Photodiode (PD333-3B/H0/L2 from EVERLIGHT) is used to receive the laser echo. The laser echo power caused by the dust scattering is then recorded.
Because of the threshold of the emitter, when the mass concentration of Point A is less than 260 mg/m3, it is difficult to record effective echo results. As a consequence, we only obtain eight groups of valid data during the stable dust environment analyzed in Section 4.2.2. The experimental results are listed in Table 4. The experiment aims to analyze the influence of the dust mass concentration on the laser echo, so light-absorbing materials are used at the end of the test area and the laser system is placed in the dust environment. Therefore, the results in Table 4 are only related to the mass concentration. Likewise, the signal amplitude goes down with the decreasing mass concentration of Point A.
Table 4. Experimental results
4.3.2 Experimental and simulated echo powers
In this section, the laser echo signal is simulated with oblate particles with ${r_{\textrm{eff}}} = 1.1$ µm, ${v_{\textrm{eff}}} = 0.1$ and $\varepsilon = 2$ [49]. The MC simulations are calculated with the Dell Precision Tower 5810, and the average run time of the simulations is 280s. The simulations are conducted using 107 photons. In order to analyze the influence of the mass concentration distribution on the backscattering, the varying concentration (Eq. 32) and the fixed concentration (Point A) with the proposed MC model are calculated and compared. Figure 18 shows the signal echo normalized by the result of $\rho (0)\textrm{ = }400\textrm{mg/}{\textrm{m}^3}$ with varying concentration. The received signal increases as the mass concentration increases in both situations (varying concentration and fixed concentration). In addition, the results for the fixed concentration are all smaller than those for the varying concentration, and this effectively explains the phenomenon found in the laser echo experiment that the same mass concentration of Point A at different stages of dust diffusion can have different backscattering amplitudes. Therefore, in the study of multiple scattering in inhomogeneous media, the influence of the changing optical properties along the light path should be considered.
Fig. 18. Comparison of normalized signal echo between varying concentration and fixed concentration with different $\rho (0)$. (a) $\rho (0)\textrm{ = }200\textrm{mg/}{\textrm{m}^3}$; (b) $\rho (0)\textrm{ = }250\textrm{mg/}{\textrm{m}^3}$;(c) $\rho (0)\textrm{ = }300\textrm{mg/}{\textrm{m}^3}$;(d) $\rho (0)\textrm{ = }400\textrm{mg/}{\textrm{m}^3}$.
To validate the accuracy of the proposed model, a comparison among the proposed model, the widely used traditional MC model (with the H-G phase function) [15,35,50–55], and the experiment is performed. The simulations are all calculated for both the varying concentration and the fixed concentration. The waveform obtained in the experiment is a photoelectric converted voltage signal, which is affected by the transmittance of the optical system, the bias voltage of the detector, and other factors, so it is unrealistic to compare the experimental and simulated waveforms directly. Therefore, we calculate and compare the normalized peak amplitudes of the simulations and the experiment, as shown in Fig. 19.
Fig. 19. Comparison of normalized peak amplitudes between simulated results and experimental results.
In order to show the difference numerically, the maximum relative errors between the simulations and the experiment in Fig. 19 are listed in Table 5. As can be seen from Fig. 19 and Table 5, compared with the traditional MC model, the curves of the proposed model are closer to the experiment results, whether for the varying concentrations or the fixed concentration. In the traditional MC model, the assumption of spherical particles and the use of the H-G scattering phase function will lead to large errors in the study of light transmission in the dust environment. In contrast, using the proposed MC method with the T-matrix method can obtain accurate results. In Table 5, the maximum relative error of the normalized peak amplitude is within 0.8% for the proposed model with varying mass concentration, which is smaller than the value (3.5%) for the proposed model with fixed mass concentration. Therefore, the consideration of the changing optical properties with the mass concentration improves the precision greatly of the simulation, indicating the effectiveness and accuracy of the continuous field MC model for the inhomogeneous dust environment.
Table 5. Maximum relative errors between simulation and experiment
Because of the relatively stable dust environment in this enclosed room, the difference between the varying concentration and the fixed mass concentration is not obvious. To further analyze the effects of the mass concentration distributions, some typical functions are calculated with the proposed model, as listed in Table 6. The four distributions below have the same average mass concentration.
Table 6. Laser echo with different mass concentration distributions
The simulated results in Table 6 further indicate that different mass concentration distributions have a great influence on the light transmission in the inhomogeneous environment. The laser echoes all increase with the increasing mass concentration, which means different mass concentrations do not affect this rule. Table 6 further illustrates the necessity of our proposed model.
5. Conclusion and summary
In this study, a Monte Carlo radiative transfer model for a continuous concentration field with an improved semianalytic approach is proposed. The method considers the influences of varying extinction coefficients with continuously changing mass concentrations on the simulation. Additionally, we use the discrete data of the scattering phase function obtained with the T-matrix method and a rejection method to sample scattering angles in the MC program. Then the effects of the fundamental particle characteristics such as particle sizes and particle shapes on the scattering properties and laser echo are analyzed and discussed. The influences of the mass concentration distribution are also analyzed with our model in this paper. The comparison of the simulation results and the actual measurements verifies the effectiveness and accuracy of our proposed approach for an inhomogeneous dust environment. This research is significant for the techniques of laser detection and measurement.
Funding
Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX20_0315); Central University Special Funding for Basic Scientific Research (30918012201); National Natural Science Foundation of China (51709147).
Acknowledgments
All authors sincerely thank M. I. Mishchenko and L. D. Travis for the multiple spheroid T-matrix code.
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. X. Zhu, L. Lu, Z. Cao, B. Zeng, and M. Xu, “Transmission matrix-based Electric field Monte Carlo study and experimental validation of the propagation characteristics of Bessel beams in turbid media,” Opt. Lett. 43(19), 4835–4838 (2018). [CrossRef]
2. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20(1), 92–98 (2003). [CrossRef]
3. E. Akkermans and G. Montambaux, “Mesoscopic physics of photons,” J. Opt. Soc. Am. B 21(1), 101–112 (2004). [CrossRef]
4. F. Martelli, S. D. Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation through Biological Tissue and other Diffusive Media: Theory, Solutions and Software (SPIE, 2009).
5. A. Isimaru, Wave propagation and scattering in random media (Academic, 1978).
6. M. I. Mishchenko, “Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation,” Appl. Opt. 39(6), 1026–1031 (2000). [CrossRef]
7. R. R. Naraghi, S. Sukhov, J. J. Sáenz, and A. Dogariu, “Near-Field Effects in Mesoscopic Light Transport,” Phys. Rev. Lett. 115(20), 203903 (2015). [CrossRef]
8. B. Saleh, Introduction to subsurface imaging (Cambridge University, 2011).
9. S. Sukhov, D. Haefner, and A. Dogariu, “Coupled dipole method for modeling optical properties of large-scale random media,” Phys. Rev. E 77(6), 066709 (2008). [CrossRef]
10. M. Xu, W. Cai, M. Lax, and R. Alfano, “Photon migration in turbid media using a cumulant approximation to radiative transfer,” Phys. Rev. E 65(6), 066609 (2002). [CrossRef]
11. J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31(30), 6535–6546 (1992). [CrossRef]
12. X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt. 7(3), 279–290 (2002). [CrossRef]
13. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef]
14. L. R. Bissonnette, “Lidar and multiple scattering,” in Lidar. Springer Series in Optical Sciences (Springer, 2005), vol. 102.
15. T. Binzoni, T. S. Leung, A. H. Gandjbakhche, D. Rüfenacht, and D. T. Delpy, “The use of the Henyey-Greenstein phase function in Monte Carlo simulations in biomedical optics,” Phys. Med. Biol. 51(17), N313–N322 (2006). [CrossRef]
16. F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source-detector separations close to one transport mean free path,” J. Opt. Soc. Am. A 16(12), 2935–2945 (1999). [CrossRef]
17. J. R. Mourant, J. Boyer, A. H. Hielscher, and I. J. Bigio, “Influence of the scattering phase function on light transport measurements in turbid media performed with small source-detector separations,” Opt. Lett. 21(7), 546–548 (1996). [CrossRef]
18. F. Vaudelle, J. P. L’huillier, and M. L. Askoura, “Light source distribution and scattering phase function influence light transport in diffuse multi-layered media,” Opt. Commun. 392, 268–281 (2017). [CrossRef]
19. L. O. Reynolds and N. J. McCormick, “Approximate two-parameter phase function for light scattering,” J. Opt. Soc. Am. 70(10), 1206–1212 (1980). [CrossRef]
20. A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H. J. Schwarsmaier, “Influence of the scattering phase function approximation on the optical properties of blood determined from the integrating sphere measurements,” J. Biomed. Opt. 4(1), 47–53 (1999). [CrossRef]
21. C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles (John Wiley & Sons, 2008).
22. B. Gélébart, E. Tinet, J. M. Tualle, and S. Avrillier, “Phase function simulation in tissue phantoms: a fractal approach,” Pure Appl. Opt. 5(4), 377–388 (1996). [CrossRef]
23. S. Dap, D. Lacroix, R. Hugon, and J. Bougdira, “Retrieving particle size and density from extinction measurement in dusty plasma, Monte Carlo inversion and Ray-tracing comparison,” [J],” J. Quant. Spectrosc. Radiat. Transfer 128, 18–26 (2013). [CrossRef]
24. M. I. Mishchenko, “Light scattering by size-shape distributions of randomly oriented axially symmetric particles of a size comparable to a wavelength,” Appl. Opt. 32(24), 4652–4666 (1993). [CrossRef]
25. D. Toublanc, “Henyey-Greenstein and Mie phase functions in Monte Carlo radiative transfer computations,” Appl. Opt. 35(18), 3270–3274 (1996). [CrossRef]
26. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Lookup table-based sampling of the phase function for Monte Carlo simulations of light propagation in turbid media,” Biomed. Opt. Express 8(3), 1895–1910 (2017). [CrossRef]
27. L. Bai, Z. Wu, and H. Li, “Sample scattering phase function method of non-spherical particles in UV band,” IEEE, 4067–4070 (2011).
28. M.I. Mishchenko, N.T. Zakharova, N.G. Khlebtsov, G. Videen, and T. Wriedt, “Comprehensive thematic T-matrix reference database: A 240–246 update,” J. Quant. Spectrosc. Radiat. Transfer 202, 240–246 (2017). [CrossRef]
29. L. Ma and R. K. Hanson, “Measurement of aerosol size distribution functions by wavelength-multiplexed laser extinction,” Appl. Phys. B 81(4), 567–576 (2005). [CrossRef]
30. L. Mei, P. Lundin, S. Andersson-Engels, S. Svanberg, and G. Somesfalean, “Characterization and validation of the frequency-modulated continuous-wave technique for assessment of photon migration in solid scattering media,” Appl. Phys. B 109(3), 467–475 (2012). [CrossRef]
31. D. W. Illig, W. D. Jemison, and L. J. Mullen, “Independent component analysis for enhancement of an FMCW optical ranging technique in turbid waters,” Appl. Opt. 55(31), C25–C33 (2016). [CrossRef]
32. H. Wang and X. Sun, “Effect of laser beam on depolarization characteristics of backscattered light in discrete random media,” Opt. Commun. 470, 125971 (2020). [CrossRef]
33. L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
34. L. Zhu, X. Wang, B. Wang, and Y. Wu, “Semi-analytical Monte Carlo simulation for time-resolved light propagating in multilayered turbid media,” J. Mod. Opt. 66(11), 1–9 (2019). [CrossRef]
35. B. Liu, C. Song, and Y. Duan, “The characteristics simulation of FMCW laser backscattering signals,” Opt. Rev. 25(2), 197–204 (2018). [CrossRef]
36. S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in tissues,” in Optical-Thermal response of laser-irradiated tissue (Springer, 1995), pp. 73–100.
37. M. Grabner and V. Kvicera, “The wavelength dependent model of extinction in fog and haze for free space optical communication,” Opt. Express 19(4), 3379–3386 (2011). [CrossRef]
38. E. J. McCartney, Optics of the atmosphere: scattering by molecules and particles (John Wiley & Sons, 1976).
39. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, “Capabilities and Limitations of a Current FORTRAN Implementation of the T-Matrix Method for Randomly Oriented, Rotationally Symmetric Scatterers,” J. Quant.Spectrosc. Radiat. Transfer 60(3), 309–324 (1998). [CrossRef]
40. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).
41. 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]
42. P. Chen, D. Pan, Z. Mao, and H. Liu, “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]
43. Q. Liu, X. Cui, W. Chen, C. Liu, J. Bai, Y. Zhang, Y. Zhou, Z. Liu, P. Xu, H. Che, and D. Liu, “A semianalytic Monte Carlo radiative transfer model for polarized oceanic lidar: experiment-based comparisons and multiple scattering effects analyses,” J. Quant. Spectrosc. Radiat. Transfer 237, 106638 (2019). [CrossRef]
44. J. F. Henry, A. Gonzalez, and L. K. Peters, “Dynamics of NH4Cl particle nucleation and growth at 253-296 K,” Aerosol Sci. Tech. 2(3), 321–339 (1982). [CrossRef]
45. Radiation Commission, and IAMAP, “A preliminary cloudless standard atmosphere for radiation computation,” World Meteorological Organization Geneva, WCP-112, WMO/TD-No.24 (1986).
46. J. E. Hansen and L. D. Travis, “Light Scattering in Planetary Atmospheres,” Space Sci. Rev. 16(4), 527–610 (1974). [CrossRef]
47. L. L. L. Waduge, S. Zigan, L. E. Stone, A. Belaidi, and P. García-Triñanesa, “Predicting concentrations of fine particles in enclosed vessels using a camera based system and CFD simulations,” Process Saf. Environ. 105, 262–273 (2017). [CrossRef]
48. F. Chen, C. M. Simon, and A. C. K. Lai, “Modeling particle distribution and deposition in indoor environments with a new drift–flux mode,” Atmospheric Environ. 40(2), 357–367 (2006). [CrossRef]
49. S. C. Hill, A. C. Hill, and P. W. Barber, “Light scattering by size/shape distributions of soil particles and spheroids,” Appl. Opt. 23(7), 1025–1031 (1984). [CrossRef]
50. S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo Model of Light Propagation in Tissue,” Proc. SPIE. IS 5, 1030509 (1989). [CrossRef]
51. P. Chervet, C. Lavigne, A. Roblin, and P. Bruscaglioni, “Effects of aerosol scattering phase function formulation on point-spread-function calculations,” Appl. Opt. 41(30), 6489–6498 (2002). [CrossRef]
52. K. I. Gjerstad, J. J. Stamnes, B. Hamre, J. K. Lotsberg, B. Yan, and K. Stamnes, “Monte Carlo and discrete-ordinate simulations of irradiances in the coupled atmosphere-ocean system,” Appl. Opt. 42(15), 2609–2622 (2003). [CrossRef]
53. X. Sun, Y. Han, and X. Shi, “Monte Carlo simulation of backscattering by a melting layer of precipitation,” Acta Phys. Sin. 26(16), 21258 (2018). [CrossRef]
54. M. Premuda, E. Palazzi, F. Ravegnani, D. Bortoli, S. Masieri, and G. Giovanelli, “MOCRA: a Monte Carlo code for the simulation of radiative transfer in the atmosphere,” Opt. Express 20(7), 7973–7993 (2012). [CrossRef]
55. H. Wang, D. Liu, Z. Song, N. Yang, and C. Yang, “Monte Carlo simulation of laser transmission through a smoke screen,” Chin,” J. Comput. Phys. 29(2), 022002 (2017). [CrossRef]
