Monte Carlo Simulation with Experimental Research about Underwater Transmission and Imaging of Laser

Abstract

Attenuation of the laser beam in underwater transmission and detection due to absorption and scattering results in a rapid reduction in energy and blurring of the image. By combining the bidirectional reflectivity distribution function (BRDF) with the Monte Carlo (MC) method, a full-link underwater imaging process model was established which comprehensively investigated the influence of water quality, transmission distance and target characteristics on imaging performance. In order to describe the transmission process of the light more accurately, by adding particles with both absorption and scattering functions in the medium, the Mie scattering theory was employed to simulate the real channel. Moreover, while setting the gate width, the pre-calibrated detector response curve was employed to build a corresponding relationship between the image grayscale and the detector collection energy, aiming to simulate the working mode of the detector in the experiment. In various imaging scenarios, the maximum relative errors between the simulated images and experimental results were within 30%, which proved the correctness of the imaging simulation model and the feasibility of the imaging MC (IMC) method to evaluate the quality of whole imaging process.

1. Introduction

In recent years, with the rapid development of technologies such as underwater sensor networks and underwater vehicles, the scientific significance, economic benefits, and strategic position of ocean surveys have attracted increasing attention. Therefore, a proper light transmission simulation method suitable for underwater detection technology is urgently needed. Laser technology has been widely used in the fields of underwater target detection [1], recognition [2] and seabed mapping [3].

The early applications of laser-based underwater target detection mainly focused on the detection of seabed depth and seabed topography by airborne light detection and ranging (Lidar). Since the development of the airborne underwater Lidar system and related technologies has matured, the detection performance of underwater targets has been gradually improved. In the field of underwater detection, researchers have carried out a series of work to improve the performance of detection such as multiple scanning methods [4], high laser repetition frequency [5], photon counting [6] and polarization detection [7], etc. For the purpose of meeting the visualization requirements in many scenarios [8], laser underwater detection technology has gradually been integrated with advanced imaging technologies such as range-gated imaging technology [9], hyperspectral imaging [10], and polarization imaging [11].

Compared with the above imaging technologies, the range-gated imaging is able to capture the target at a specific distance based on controlling the relative delay between the detector gate and the laser pulse, which can filter out the influence of scattered radiation during the propagation process [12]. In addition, range-gated imaging has the advantages of a high S/N ratio, strong anti-interference ability, fast imaging speed, a simple structure and being unaffected by ambient light. It has become the most widely used underwater imaging technology against the current technical background [13,14]. Although, due to the inherent optical properties of seawater, light will be absorbed and scattered when transmitted in seawater, and the width of the time-domain impulse response caused by the multipath effect increases [2], which makes the propagating procedure difficult to solve by conventional and deterministic mathematical methods.

At present, the main analysis methods of light scattering in water are divided into three types: analytical method, experimental method and MC simulation method [15]. Among them, the analytical method needs to set strict initial conditions and boundary conditions to obtain the analytical results. Different results will lead to different calculation results. On this basis, scholars proposed the use of the phenotype radar equation for analysis and calculation. This method has a small amount of calculation and intuitive results, but it is not suitable for quantitative comparative analysis of real situations in practical applications [16]. In addition, the analytical scope of experimental methods is often limited. On the basis of adding the auxiliary verification of experimental analysis, the method of cooperating with MC has been proved to be effective [17].

The MC method, which is a category of computational method that involves the random sampling of a physical quantity, is widely used for its flexibility in many different fields. It does not need to introduce constraints, boundary conditions or artificial assumptions, and can qualitatively analyze the trend of the detection system parameters and the influence of environmental changes on the results [18]. The laws are consistent with the real situation and can predict the actual results to a certain extent [19]. Therefore, MC simulation has become the main method for analyzing and studying light scattering in various scattering media.

The MC method is mainly concentrated on the fields of wireless optical communication, underwater target detection and remote sensing, of which, research into light transmission in optical communication is the earliest, dating back to the 1960s [20]. In recent years, Zhang proposed a more general random channel model for turbid seawater, and based on the model, the path loss and impulse response of the channel were evaluated [21]. In 2017, Sahu quantified the received power loss and channel dispersion for different water types, transmission distances, divergence angles and receiving angles based on different scattering phase function models, and affirmed the accuracy and stability of the SS phase function [22]. In 2020, Yuan illustrated and improved the convergence of the MC integration model by introducing different sampling methods. The MC integration model based on partial importance sampling has higher computational efficiency in high-order scattering scenarios [23].

However, the research into underwater target detection based on the MC method started relatively late. On the basis of Wang’s research [24], in 2005, the first systematic research into underwater target detection was carried out using standard MC; three different MC programs were proposed that keep track of the status of polarization after every scattering event were described. Based on the research of polarization, a backscattered light coherence investigation was added, an electric field MC (EMC) approach for directly simulating coherent backscattering (CBS) of coherent or partially coherent polarized light was proposed by Xu [25]. To further improve the EMC model, a transmission matrix-based electric field MC (TEMC) method was introduced to study the propagation characteristics of Bessel beams with different orbital angular momentum (OAM) in turbid media in which beam transmission patterns, intensity attenuation, and the degree of polarization (DOP) through turbid media of varying thickness were analyzed [26].Based on a semi-analytic MC technique, in 2019, Liu built an MC model to simulate Lidar signals based on the simultaneous inherent in situ optical properties of seawater. The Lidar measurements and the MC simulation can provide both the Lidar signals and the retrieved Lidar attenuation coefficient α. In 2020, Chen proposed an Oceanic Lidar Emulator (OLE), which is an improved semi-analytic MC-based Lidar simulator with the capability of dealing with the physics of light propagating through a wind-driven rough air-sea interface and into the ocean, being scattered by subsurface phytoplankton and reflected from the sea bottom, and returning to the Lidar receiver.

As an important part of underwater detection, the reflection characteristic of the target is a very key factor affecting the quality of detection and imaging. To describe the underwater imaging process more accurately, with the basic principle of MC simulation, it is an effective method to organically combine BRDF based on the microsurface principle. Some scholars have carried out related research work; in 2006, Huang designed the random process between photons and crop canopy, using field measured BRDF data to validate MC model. The results showed that the simulated BRDF with MC model was in good accordance with the field data [27]. In 2015, Wang proposed a wetland aquatic vegetation canopy BRDF-MC model, which involved determining track photon value, and performing an aquatic vegetation total weight calculating process corresponding to track the photon value. The method enabled the realization of the cap layer structure parameter-obtaining process and the BRDF-MC model-establishing process in a precise manner [28]. Recently, model-optimization algorithms based on the two theories have been successively proposed and applied in the fields of radiation uncertainty classification [29] and 3D canopy modeling [30].

In summary, the scope of MC mainly focuses on the simulation of transmission attenuation, and most of the research relating to the combination of MC and BRDF stays in the field of space remote sensing, while its application in underwater detection is rare. In our research, the IMC method was established to simulate the underwater target range-gated detection performance of different targets at various ranges and media. The structural similarity (SSIM), mean structural similarity (MSSIM) and the peak signal-to-noise ratio (PSNR) were calculated to evaluate the imaging quality, and the trend curves of the detection signals with the water parameters were also obtained. Finally, the experiment results indicated the applicability of the proposed transport imaging model under different water types and target characteristics.

2. Materials and Methods

The proposed IMC method comprised five main sections: initialization, movement, reflection, gathering and scattering. Figure 1 shows the flow chart of the IMC method and the specific algorithm process is as follows:

Initialization: we assume that a beam of Gaussian laser emerges from the XOY plane of Cartesian coordinates, and the limited size beam propagates along the z-axis from the plane with z = 0. The divergence angle ${\theta }_0$ and the radius of spot $r$ follow the Gaussian distribution, while the azimuths of location $\phi$ and the azimuths of propagation ${\phi }_0$ are uniformly distributed over $\left[0, 2\pi \right]$. Hence, the initial coordinates and the initial direction of sampled packets can be obtained as:

$$\{\begin{array}{c}x=r\cos \left(\phi \right)\\ y=r\sin \left(\phi \right)\\ z=0\end{array}$$

(1)

$$\{\begin{array}{c}{u}_x=\sin \left({\theta }_0\right)\cos \left({\phi }_0\right)\\ u_y=\sin \left({\theta }_0\right)\sin \left({\phi }_0\right)\\ u_z=\cos \left({\theta }_0\right)\end{array}$$

(2)

Considering the refraction of beam on the surface of the emission window, the direction of refracted light could be determined according to the Law of Snell.

Movement: photon packets will travel a random distance $l_R$ along the direction of transmission while the number of scattering is less than the preset maximum, which is called a free path. The sampling function of it is:

$$l_R=-\ln \left(\varsigma \right)/c$$

(3)

where $\varsigma$ is uniformly distributed over $\left[0, 1\right]$ and $c$ is the attenuation coefficient of scattering media, which is the sum of absorption coefficient $a$ and scattering coefficient $b$. To describe the inherent optical properties of medium more precisely, we solve for the coefficients $a$, $b$ and $c$ according to the Mie scattering theory [31,32].

The extinction efficiency factor $Q_e$ and the scattering efficiency factor $Q_s$ represent the ratio of extinction and scattering of the energy incident on the geometric cross-section of the particle, respectively.

$$\{\begin{array}{l}{Q}_e\left(x,m\right)=\frac{2}{x^2}{\displaystyle \sum _{n=1}^{\infty }\left(2n+1\right)}\mathrm{Re}\left[a_n+b_n\right]\\ Q_s\left(x,m\right)=\frac{2}{x^2}{\displaystyle \sum _{n=1}^{\infty }\left(2n+1\right)}\left({\left|a_n\right|}^2+{\left|b_n\right|}^2\right)\end{array}$$

(4)

Factors $a_n$ and $b_n$ are defined by:

$$\{\begin{array}{l}{a}_n\left(x,m\right)=\frac{{\psi }_n^{\prime }\left(mx\right){\psi }_n\left(x\right)-m{\psi }_n\left(mx\right){\psi }_n^{\prime }\left(x\right)}{{\psi }_n^{\prime }\left(mx\right){\xi }_n\left(x\right)-m{\psi }_n\left(mx\right){\xi }_n^{\prime }\left(x\right)}\\ b_n\left(x,m\right)=\frac{m{\psi }_n^{\prime }\left(mx\right){\psi }_n\left(x\right)-{\psi }_n\left(mx\right){\psi }_n^{\prime }\left(x\right)}{m{\psi }_n^{\prime }\left(mx\right){\xi }_n\left(x\right)-{\psi }_n\left(mx\right){\xi }_n^{\prime }\left(x\right)}\end{array}$$

(5)

where the size parameter $x=2\pi r/\lambda$ represents the ratio of the particle perimeter to the incident wavelength and $m$ is the refractive index. ${\xi }_n\left(x\right)$ and ${\psi }_n\left(x\right)$ are the Riccati-Bessel function described by spherical Bessel Function of the first kind $J_{n+1/2}$ and half-integer orders Hankel functions of the second kind $H_{n+1/2}^{\left(2\right)}$, respectively.

$$\{\begin{array}{l}{\psi }_n\left(x\right)=\sqrt{\pi x/2}{J}_{n+1/2}\left(x\right)\\ {\xi }_n\left(x\right)=\sqrt{\pi x/2}{H}_{n+1/2}^{\left(2\right)}\left(x\right)\end{array}$$

(6)

The extinction and scattering efficiency factor of a single scattering particle to photons are obtained, but the key to characterizing the extinction and scattering specialty of the media are the extinction and scattering cross section within the unit volume of the medium, namely attenuation coefficient $c$ and scattering coefficient $b$. Suppose the particle concentration is uniformly distributed, and given the number of particles per unit volume $\rho$ and the mean radius of particles $r$, the coefficients could be expressed as:

$$\{\begin{array}{l}c=\rho Q_e\pi r^2\\ b=\rho Q_s\pi r^2\end{array}$$

(7)

Therefore, the free path $l_R$ could be determined by attenuation coefficient c. Furthermore, updates to the transmission direction, position and time-of-flight of photon packets are accompanied by each movement.

Reflection: the process of which is closely related to the material properties of the target surface and the direction of light incidence [33]. In order to describe the reflected energy distribution of the target, the bidirectional reflectance distribution function (BRDF) model proposed by Cook and Torrance is used in this study to divide the target surface into microfacets satisfying the Gaussian distribution, and each micro element follows a Fresnel reflection formula [34,35]. The reflection $f_s$ is expressed as:

$$f_s=\frac{F\left({\theta }_i,{\theta }_r\right)G\left({\theta }_i,{\theta }_r,\Delta \phi \right)}{\pi \cos {\theta }_i\cos {\theta }_r}$$

(8)

$F\left({\theta }_i,{\theta }_r\right)$ is the Fresnel reflection coefficient, and $G\left({\theta }_i,{\theta }_r,\mathsf{\Delta }\phi \right)$ represents the geometric attenuation factor that the reflected light is attenuated due to the mutual occlusion of the microfacets, including blocking the incident light and the outgoing light, the following equation is used:

$$G\left({\theta }_i,{\theta }_r,\Delta \phi \right)=\min \left(1;\frac{2\cos {\theta }_N\cos {\theta }_r}{\cos \beta };\frac{2\cos {\theta }_N\cos {\theta }_i}{\cos \beta }\right)$$

(9)

$\mathsf{\Delta }\mathsf{\phi }$ is the angle between the incident light and the reflected light at the microfacet, $\beta$ is the angle between the incident light and the normal to the micro surface, ${\theta }_i$, ${\theta }_r$ and ${\theta }_N$ are the incident angle, the reflection angle and the angle between the normal direction of the microfacet and the normal direction of the target surface, respectively. The diffuse reflection part $f_d$ is expressed as:

$$f_d={\rho }_d{G}_d$$

(10)

${\rho }_d$ is the light intensity ratio between diffuse and full Lambertian surface scattering. Here, $R_i$ is the Fresnel reflectivity of the surface, which can be calculated as $R_i=\left(r_s+r_p\right)/2$. $K_2$ is the Fresnel reflection of diffuse light on the material surface and is usually assumed to be approximately equal to $R_i$. $R_{\infty }$ is the reflectivity coefficient when the surface thickness of the object is seen as infinite, which can be seen as a constant for a given material at a fixed incident wavelength. $G_d$ represents the angular distribution of the Lambert diffuse scattering intensity on the surface of the target [36].

$${\rho }_d=\left(1-R_i\right)\frac{\left(1-K_2\right)R_{\infty }}{1-K_2{R}_{\infty }}$$

(11)

$$G_d=\frac{\cos {\theta }_r}{\pi }$$

(12)

Thus, the residual weight of photon packets in the direction of reflection can be determined according to Equations (8) and (10). The condition for identifying photon packets reach target is whether the projection of the distance on the z-axis exceeds the position of target. Once the condition is met, it means that the target was incident in the last random step, and the current position needs to be corrected to make it return to the target plane.

$$\{\begin{array}{c}X=X^{\prime }-l_{Rx}(Z^{\prime }-L)/l_{Rz}\hfill \\ Y=Y^{\prime }-l_{Ry}(Z^{\prime }-L)/l_{Rz}\hfill \\ Z=L\end{array}$$

(13)

in which $X^{\prime },Y^{\prime }$ and $Z^{\prime }$ are the original coordinates while $X,Y$ and $Z$ are the corrections in the target plane, whose distance from origin is $L$. $l_{Rx}$, $l_{Ry}$ and $l_{Rz}$ are the projection of propagation direction on $X,Y$ and $Z$ axis.

Gathering: the receiver is placed at $\left(l_x,l_y,l_z\right)$ with a specific aperture, focal length and field of view (FOV), and the weights of photon packets are accumulated while the following conditions are met:

$$\{\begin{array}{c}x\le l_x\\ y\le l_y\\ z\le l_z\\ u_z^{\prime }\ge \cos \left(\mathrm{FOV}\right)\end{array}$$

(14)

Subsequently, the image is obtained according to the calibrated response curve, and the SSIM, MSSIM and PSNR indicators are used as the evaluation standard for it [37,38]. SSIM is used to measure the similarities between two given images in three dimensions: brightness, contrast and structure. The specific form of SSIM index is:

$$\mathrm{SSIM}\left(x,y\right)=\frac{\left(2{\mu }_x{\mu }_y+C_1\right)\left(2{\sigma }_{xy}+C_2\right)}{\left({\mu }_x^2+{\mu }_y^2+C_1\right)\left({\sigma }_x^2+{\sigma }_y^2+C_2\right)}$$

(15)

${\mu }_x$, ${\mu }_y$, ${\sigma }_x^2$, ${\sigma }_y^2$ and ${\sigma }_x{\sigma }_y$ are the average value, variance and covariance of two nonnegative image signals $x$ and $y$, respectively. $C_1$ and $C_2$ are constants to maintain stability with $C_1={\left(0.01L\right)}^2$ and $C_2={\left(0.03L\right)}^2$, where $L$ is the dynamic range of the pixel values. Simultaneously, the MSSIM index is used to evaluate the overall image quality from local image features:

$$\mathrm{MSSIM}\left(X,Y\right)=\frac{1}{MN}{\displaystyle \sum _{i=1}^M{\displaystyle \sum _{j=1}^N\mathrm{SSIM}\left(x_i,y_j\right)}}$$

(16)

where $X$ and $Y$ are the reference and the distorted images; $x_j$ and $y_j$ are the image contents at the $i$th row and $j$th column local window; $M$, $N$ is the number of local windows of the image. Here, we use a symmetric Gaussian filter matrix of $11\ast 11$ with a standard deviation of 1.5 as a weighting window to filter the image, which is equivalent to sliding the window on the distorted image, so the MSSIM can be obtained by averaging the matrix mapped by the local SSIM index. PSNR is used to measure the quality of the degraded image:

$$\mathrm{PSNR}=10{\log }_{10}\left(\frac{L^2}{\mathrm{MSE}}\right)$$

(17)

$$\mathrm{MSE}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(x_i-y_i\right)}^2}$$

(18)

where MSE is mean squared error, $N$ is the total number of pixels in the image. Generally, the larger the SSIM, MSSIM and PSNR, the higher the quality of the degraded image and the more similar the two images are.

Scattering: photon packets collide after the free path $l_R$ with the scattering angle θ and the azimuth angle $\gamma$. The transmission direction after collision converted into the global coordinate system is expressed as $\left(u_x^{\prime },u_y^{\prime },u_z^{\prime }\right)$, denoted by the original direction $\left(u_x,u_y,u_z\right)$:

$$\left[\begin{array}{c}{u}_x^{\prime }\\ u_y^{\prime }\\ u_z^{\prime }\end{array}\right]=\left[\begin{array}{ccc}{u}_x{u}_z/\sqrt{1-u_z^2}& -u_y/\sqrt{1-u_z^2}& u_x\\ u_y{u}_z/\sqrt{1-u_z^2}& u_x/\sqrt{1-u_z^2}& u_y\\ -\sqrt{1-u_z^2}& 0& u_z\end{array}\right]\left[\begin{array}{c}\sin \theta \cos \gamma \\ \sin \theta \sin \gamma \\ \cos \theta \end{array}\right]$$

(19)

Scattering is also accompanied by the loss of energy: $E_s=\omega E_{pre}$ where the albedo $\omega =b/c$, $E_s$ is the energy after collision while $E_{pre}$ is the energy before it.

The bidirectional transmission of photons is a cyclic process with the termination conditions:

(1)

The photon packet returns to the detector within the field of view;

(2)

The photon packet reaches the target plane but does not hit the target area;

(3)

The photon packet energy is less than the threshold, namely, it is lost during the transmission process;

(4)

The number of photon packet scattering exceeds the preset maximum number.

The initialized photon packet participates in the calculation of the modules in turn, once the termination conditions are met, the tracking is interrupted and the next photon packet is emitted.

3. Results

3.1. Simulation

3.1.1. Unidirectional Transmission

In the simulation experiment, the influence of different beam parameters on the propagation result was investigated firstly by the forward transmission part of the IMC method. The medium in the simulation was set to Jerlov II type water, which is an oceanic water type and mostly free of extra colored dissolved organic matter from terrigenous sources, resulting in small signal attenuation, and the detailed values of $c$ and $\omega$ were $0.42 {\mathrm{m}}^{-1}$ and $0.89$ [39] The scattering angle θ was determined by the Henyey-Greenstein (H-G) scattering phase function [40], where $\gamma$ satisfies the uniform distribution of $\left[0,2\pi \right]$. To ensure that in underwater laser communication or detection applications, the laser transmitted through the channel can not only illuminate the core area, but also maximize the energy at the target or receiver, the impact of divergence angle and spot radius in the initial on spatial and temporal variation characteristics after unidirectional transmission was analyzed emphatically in the simulation. The target was located at different distances, transmission results of spatial broadening and temporal broadening on the target plane are shown in the figure below where the pulse width (full width at half maximum (FWHM)) was 4.5 ns, the divergence angles were 20, 60, 100 and 200 mrad, with the initial radius of 0.5, 2.5 and 10 mm (FWHM), respectively. The number of simulated photon packets was 1 × 10⁸ and each packet contained 5.35 × 10⁵ photons (each photon has an energy of 3.74 × 10¹⁹ $J$ at a wavelength of 532 nm, thus, the total energy can be calculated as 1 × 10⁸ × 5.35 × 10⁵ × 3.74 × 10⁻¹⁹ $J$ = 20 × 10⁻⁶ $J$).

As the transmission length increased, the spot and pulse width of all beams were broadened, but it is worth noting that when the divergence angle was the same, the effect of the spatial broadening of lasers with different initial radius at the same length was extremely weak and could be ignored. While the changes of divergence angle bring contrasting results: when the distance was relatively short, the different divergence angles had significant differences in the effect of spatial broadening, the larger the divergence angle, the more obvious the broadening. With the increase in transmission distance, the results of smaller divergence angles gradually tended to be the same, such as 20 mrad (solid blue line) and 60 mrad (dashed blue line); however, larger divergence angles within a certain range maintained more distinct broadening than smaller ones, as shown by 200 mrad (densely dashdotted blue line) and 20 mrad in Figure 2. Conversely, as the red line shown in Figure 2, the temporal broadening effect of the laser was almost unaffected by the divergence angle and the initial radius. For simplicity, only the temporal broadening results of 20 and 200 mrad are shown here. Consequently, it was necessary to comprehensively consider the application scenario and experimental conditions to select the most appropriate initial parameters to ensure the best results.

3.1.2. Bidirectional Transmission

In order to explore the effect of transmission distance on the laser detection capability, a complex underwater simulation environment was constructed. Relevant study has shown that the minimum particle size in the surface of the South China Sea is about 3 μm (Zheng et al.) [41], so we chose calcium carbonate powder with a relative refractive index of 1.59 and the size of 3 μm as scattering particles. We controlled the detection distance to be 1 m in our simulation and changed the attenuation coefficient of the medium by adjusting the quantity of the added particles, that is, adjusting the attenuation lengths (AL). In order to describe the corresponding attenuation coefficient and scattering angle distribution during transmission more realistically, the properties of the calcium carbonate powder and the attenuation efficiency factor $Q_e$ and scattering efficiency factor $Q_s$ represented by Equation (4) were combined to obtain the corresponding attenuation coefficient through the calculation of Equation (7).

In the underwater scattering medium set as above, the IMC method was executed to simulate the underwater detection of artificial targets: “Diver” and “Double Slit”. They were binary images with two regions, the image part (A) and the background part (B). The normal of the microfacet follows a Gaussian distribution and can be obtained by the Marsaglia-Bray algorithm, where the expectation is 0, and the standard deviation is related to the surface roughness. The approximation of the silver plated part A is the specular reflection $f_s$ with roughness of 0.001 and part B is based on the diffuse reflection $f_d$ [36]. Since we adopted the coaxial structure of transceiver, and the coordinate transformation of photons in the scattering process was always kept in the global coordinate system, the image coordinate system and pixel coordinate system of the detector were also established under it, so the imaging results of the ideal optical system in the object space could be directly described by the position of the photon. Targets with side length of 5 cm × 5 cm were placed in the water tank 1 m away from the incident window, and the image resolution was 200 × 200. To ensure they could be illuminated, the laser properties were confirmed by forward transmission simulation of IMC, where the divergence angle was 200 mrad and the initial spot diameter was 1 mm.

3.1.3. Optimization of IMC Method

During the process of underwater transmission, the backscattering light of particles return the imaging system and introduce noise, especially when imaging long-distance targets in complex media. In the IMC method, the detector and laser are coaxially placed, as the power of illumination increases, so does the intensity of backscattering, which results in laser active imaging systems being unable to detect objects effectively. The detected normalized intensities of the diffuse planar target under AL of 0.42 and 1.42 by IMC model before optimization are shown in Figure 3a,c. It can be seen from the echo timing, in which the black dashed rectangular area was the response of the detector to backscattering, while the red one was the reflected signal of the target. It indicates that with the increase in AL, the influence of the backscattering is enhanced and part of the reflected signals from targets are submerged in the invalid scattering signals, which greatly limits the underwater detection distance.

With the development of a precise delay in time domain and ultra-fast shutter, the employment of image intensifiers in range-gated technology ensured that the imaging system only received signals in the vicinity of the target and would not be affected by scattering signals on the transmission path [9]. It is an effective method to overcome the backscattering of medium, which can improve the contrast and S/N ratio of the detected image. Therefore, referring to the principle of the technology, the echo signals within the gating range were marked effective in the optimized IMC simulation, whereas the photons with channel information, in other words, the backscattered part, was marked as invalid photons. The normalized intensities depicted in Figure 3b,d were obtained by counting the response of the detector to the effective photons exclusively. Comparing the simulation results before and after optimization, the influence of the backscattered photons to the detector was removed, which described the underwater detection process of the range-gated system appropriately.

The detector was divided into grids in the IMC simulation to separately count the photon energy falling in each bin, and the response curve of the camera was used to establish an illuminance-grayscale lookup table, so the simulation image could be obtained. Original and optimized simulated images of “Diver” and “Double Slit” with AL equals 0.42 and 1.42 are shown in Figure 4. With the increase in the AL, the photon energy returning to the detector decreased, and obvious background noise appeared in the non-target area with the loss of image details. The main reason for the phenomena is that the absorption and scattering of the water lead to the attenuation of the laser energy during the transmission process, which, accompanied by the multipath effects, reduces the received photon energy, lowers the contrast and degrades the image quality.

To further explore the relationship between the detection effect and water quality, and to verify the reliability of the simulation, we built a range-gated imaging system. While obtaining the images corresponding to the simulated situation, imaging experiments in the water with continuously changing ALs were also conducted. Moreover, the influence of the transmission distance on the photon energy returning to the detector were studied.

3.2. Experimental Setup

On the basis of simulation, the selected water tank made of high transparency glass is 1.5 m in length, 0.42 m in width and 0.4 m in water depth. The single pulse energy is 20 μJ, the beam is expanded into the water tank and then travels 1 m to hit targets. The reflected beam returns to the detector, which consists of an image intensifier and a complementary metal oxide semiconductor (CMOS) camera with a maximum resolution of 1024 × 1024 and a frame frequency of 1000. The image intensifier has a gate width of 3 ns and includes a double-layer micro channel plate (MCP) to enable the maximum gain. However, the image intensifier has a power limit on the internal power supply to protect MCP and the phosphor screen, so once the detection light intensity exceeds the preset threshold, the image enhancement effect is invalid. Considering this issue, there was about 20 cm of spatial separation between the ICMOS and the laser to minimize the overlap between them and to ensure that the image intensifier was not affected by the reflected light from the glass pool. After being triggered by the signal generator, the synchronous output signal obtained by the Si-PIN detector inside the laser was used as the external trigger of the signal generator, and the delay of the timing sequence provided by the generator to the intensifier was adjusted based on it. When the reflected signal from the target was about to reach the detector, the gate of the intensifier was opened and the camera exposed to record information afterwards. However, the gate was closed for the rest of the time and the camera could not work, which overcame the strong backscattering in the water. The specific parameters were consistent with those in simulation and the experimental setup with timing diagram of synchronization signals is depicted as Figure 5.

Considering the insoluble properties of calcium carbonate, a pump was applied until the medium was uniform and two targets placed in the designated position were sequentially detected; the results are shown in Figure 6. As AL increased, the energy returning to the detector decreased due to absorption and scattering; it is displayed in the images that the overall gray value became smaller, more noises appeared around the edges and became blurred. The same phenomenon as the simulation was also observed in the experimental results.

Taking the image at AL of 0.42 as the reference, we calculated the image quality evaluation indicators of the unoptimized and optimized simulation and experimental results with AL of 1.42, including structural similarity (SSIM), mean structural similarity (MSSIM) and peak signal-to-noise ratio (PSNR). The statistical results are shown in

Table 1. Comparison of image quality evaluation indicators between simulation and experimental results with the change of AL.

Target	Source	SSIM	MSSIM	PSNR/dB
Diver	Unoptimized simulation	0.494	0.999	29.71
	Optimized simulation	0.527	0.999	30.05
	Experiment	0.534	0.992	27.04
Double slit	Unoptimized simulation	0.503	0.999	23.72
	Optimized simulation	0.53	0.998	24.44
	Experiment	0.553	0.993	22.87

. It can be seen that the optimized results not only outperformed the unoptimized simulation results on SSIM and PSNR, but also the optimized simulation results were closer to the experimental results. At the same time, it is also intuitively seen from the figures that the increase in the AL led to obvious degradation of the image quality. Notably, since the Gaussian filter matrix is used to smooth and denoise the image in the process of solving the MSSIM index, the SSIM was directly solved on the original image. It led to the MSSIM index of the processed distorted image being closer to 1 than the SSIM index, that was closer to the original image.

Furthermore, we used the nickel alloy targets with a roughness of 0.01 and 0.1, respectively, and increased the width of the target seam to investigate the effect of roughness on the simulation and experimental results while AL was 0.42. The images are shown in Figure 7 and the statistical results with the roughness of 0.01 as a reference are shown in

Table 2. Comparison of image quality evaluation indicators between simulation and experimental results with the change of roughness.

Target	Source	SSIM	MSSIM	PSNR/dB
Double Slit *	Simulation	0.318	0.998	15.31
	Experiment	0.367	0.989	14.07

* The width of the target seam is widened.

.

Figure 7 shows that due to the change of the roughness of the target, the fluctuation of the microsurface elements on the target surface increased, so that the photons reflected at the target became divergent, reducing the brightness and contrast of the image. Since the transmission and reception were coaxial in the simulation, the background noise would be stronger than the experimental results, but the longitudinal comparison results in Table 2 prove the applicability of the IMC method to the simulation of different roughness. Finally, in order to quantitatively investigate the impact of water quality on imaging, we continued to carry out experiments under different ALs and supplemented with simulation analysis.

4. Discussion

It is worth noting that, for the purpose of ensuring the visibility of all images in the experiments of this paper, the gain of the intensifier was adjusted according to the attenuation coefficient of the medium. The calibrated relative magnification of the image intensifier at different gains should be used to convert grayscale images into the energy of the photosensitive surface before comparison. On the other hand, since the weight of each photon packet was assigned to 1 at the initialization of the simulation, the corresponding energy could be obtained by counting the received effective photon packets and accumulating them. Comparing the simulation and experimental results of the two targets in different ALs, it can be seen in Figure 8 that as AL increased, the energy declined significantly. For convenience, the energy in each case was normalized to the value at AL of 0.42. Because the transceiver was separated to protect the image intensifier in the experimental setup, the overlap factor ($\mathrm{O}\left(\mathrm{z}\right)$) in the light detection and ranging (LiDAR) equation was less than 1, and considering the measurement error and systematic error in the experiment, theoretically the received energy in the simulation should be higher than the value in the experiment. In fact, this phenomenon is also observed in Figure 8, and it maintained a slight discrimination under different water conditions, which means that the results of IMC method are basically consistent with the experimental results.

There are several studies on the underwater transmission MC model, among them Chen et al. [42] developed a semi-analytic MC radiative transfer model to study laser propagation in sea water and obtained the chlorophyll α concentration profile; the lookup table method has been shown to be appropriate for complex scattering phase function. Compared with the lidar experimental data of the Sanya Bay in 2018, the relative maximum error was about 40%, while the experiment proved that the relative maximum error of our IMC simulation results did not exceed 30%. In general, the results indicate that the IMC method is effective for researching laser transmission, target detection and imaging in complex underwater medium, and significantly improves the modeling accuracy.

5. Conclusions

Based on the MC method, a simulation method IMC was established for underwater imaging. Firstly, the influence of laser divergence angle and initial spot radius on the transmission effect was clarified by analyzing the spatiotemporal broadening effect of the laser with different initialization parameters after unidirectional transmission. Then, we used the Mie theory to relate the scattering medium in the simulation to the real situation. In order to collect the reflected signal carrying the target information more effectively, the IMC method was optimized to a range-gated simulation model and performed an imaging simulation of specific targets. It showed that the optimized method reduced the influence of backscattering caused by the medium on detection. It is an effective method to overcome the backscattering of a medium, which can improve the contrast and signal-to-noise ratio of the detected image. Furthermore, the influence of roughness on imaging results were analyzed, and the applicability of IMC to targets with different roughness was verified. Finally, the applicability of the IMC method was demonstrated by well-matched quantitative comparison results under continuously changing AL, where the maximum error between the simulation and experimental results was less than 30%.

In summary, the proposed IMC method establishes a channel simulation model based on the Mie theory, connects the scattering environment in the simulation with the real situation, and uses the time of flight (TOF) of photon to realize gated reception. Supplemented by modeling the detection process to simulate the whole imaging link, it can be considered that our proposed IMC simulation method has the ability to describe the transmission process of the Gaussian light field in complex underwater environments and the imaging process for two-dimensional simple targets. Therefore, the unidirectional transmission process of the IMC method could determine the optimal beam parameters in the pre-design period of the underwater communication system, and the complete bidirectional transmission IMC process lays a research foundation for the establishment of the imaging detection model for complex underwater targets.