Compared with underwater acoustic imaging technologies, underwater optical imaging technology has the advantages of high imaging resolution, high imaging frame frequency, good imaging quality, and simple structure^[1-3]. Because of these advantages, underwater optical imaging becomes a hotspot of research and is applied in the fields of military and civilian. However, the environment of underwater imaging is relatively complex and there are many factors affecting the quality of underwater optical imaging. Moreover, it is quite difficult to quantitatively describe the effects of these factors on imaging quality by simple mathematical formulas, which results in underwater optical imaging challenges. Among the factors that affect the quality of underwater optical imaging, the forward scattering and back scattering are the primary factors to reduce the quality of imaging^[4-5]. In order to reduce the effects of back scattering and absorption on the attenuation of the light, researchers focus on a new imaging scheme, which is so-called ghost imaging (GI) or single pixel imaging (SPI). Hereafter, both GI and SPI are referred to as SPI to avoid confusion. Because of great imaging performance, the SPI has been utilized in a variety of imaging field, such as Remote sensing^[6-7], Laser Radar^[8-9], polarization imaging^[10], Hyperspectral imaging^[11], et. al. The new imaging scheme has the advantages of high spatial resolution, strong robustness, multi-wavelength imaging, non-local imaging and high imaging range. Especially, SPI can improve the signal-noise-ratio (SNR) of imaging system^[12-13]. Therefore, it is suitable for underwater optical imaging by greatly overcoming the influence of scattering and absorption on the imaging range.

With further research of underwater SPI, there are many researches to be reported^[14-18], such as the influence of water turbidity on underwater SPI, the influence of water depth on underwater SPI and underwater polarization SPI. Up to now, underwater Fourier single pixel imaging (FSPI) has not been investigated yet and the water scattering effects of the imaging quality of underwater FSPI is lacking. The influence of forward scattering on underwater FSPI system is analyzed in the paper. Through analyzing the influence of forward scattering, it is found that the forward scattering of water can make the distribution of speckle pattern distorted, thus reducing the spatial resolution of underwater FSPI reconstructed results. In order to reduce the effect of water forward scattering on imaging quality, an underwater FSPI method based on water degradation function compensation method is proposed. In the proposed method, the model of water degradation function is assumed based on the degradation function of atmosphere turbulence. According to the spatial spectrum obtained by FSPI system, the water degradation function can be fitted. Utilizing the water degradation function, the real spatial spectrum distribution of target can be solved inversely. With the help of the inverse Fourier transform, the image of target can be reconstructed. The proposed method has high resistance to water forward scattering and higher spatial resolution than conventional FSPI. Therefore, the proposed method is suitable to apply in the FSPI system to reconstruct the image of underwater targets.

The FSPI system has the high imaging quality and imaging efficiency^[19-22]. The scheme of the FSPI system is shown in Fig.1. In the FSPI system, the laser illuminates on the spatial light modulator after being expanded by the lens. The computer pre-generates the illumination speckle patterns and sends them to the spatial light modulator to modulate the distribution of speckle. The difference of FSPI and other SPI system is that the pre-generated speckle patterns are sinusoidal structured patterns, which are designed according to the measured position of the target spatial spectrum. The modulated speckle patterns are emitted by the transmitting antenna and travel through the water. After arriving the position where the targets are located, the speckle patterns are reflected by the targets. A receiving antenna is used to collect the reflected light. The collected light is then focused on a bucket detector which has no spatial resolution. The measurements of the bucket are sent to the computational system to reconstruct the distribution of targets spatial spectrum.

The pre-generated modulation patterns in FSPI system is generally designed as:

${P_\phi }\left( {x,y;{f_x},{f_y}} \right) = a + b \cdot \cos \left( {2\pi {f_x}x + 2\pi {f_y}y + \phi } \right)$ (1)

Here, x and y represent Cartesian coordinates in the space domain. a is the average intensity of the modulation pattern, and b represent the amplitude of the modulation pattern which is equal to the contrast. ϕ is the initial phase and (f_x, f_y) is the spatial frequency. To simplify the derivation, the intensity of laser which illuminates on the spatial light modulator is uniform and the intensity is denoted as E₀. Therefore, the distribution of speckle emitted by the transmitting antenna is:

$E\left( {x,y} \right) = {P_\phi }\left( {x,y;{f_x},{f_y}} \right){E_0}$ (2)

The speckle patterns travel through the water, and illuminate the targets. The light is mainly absorbed and scattered by water in the process of propagation^[23-24]. The absorption of water is mainly reflected in the attenuation of light energy in water. The scattering includes back scattering and forward scattering. The back scattering further attenuates the energy of the modulated speckles and enhance the noise in the bucket detector. The noise is denoted as n. The effect of forward scattering on FSPI system is mainly reflected in the distortion of the modulated speckles which illuminates on the targets. The underwater environment is relatively complex. It is difficult to express the distortion of light distribution accurately by mathematical model. It is generally believed that the point spread function of water is of Gaussian-like type^[25-26]. The atmosphere turbulence model is usually utilized to approximate water degradation function. According to the degradation model of atmosphere turbulence, the degradation model of water can be expressed as:

$H\left( {{f_x},{f_y}} \right) = \exp \left[ { - k{{\left( {{f_x}^2 + {f_y}^2} \right)}^s}} \right]$ (3)

where k is blurring factor and s is turbulent constant. In addition, α₁ is the attenuation proportional coefficient, which can be used to represent the attenuation of light energy due to water absorption and back scattering. The distribution of light illuminating on the targets can be expressed as:

$\begin{split} E'\left( {x,y} \right) =& {\alpha _1}{P_\phi }\left( {x,y;{f_x},{f_y}} \right) \cdot\\ &{{E_0} * {{\mathbb F}^{ - 1}}\left[ {H\left( {{f_x},{f_y}} \right)} \right]} \end{split} $ (4)

Here, * represents the convolution calculation, and ${{\mathbb F}^{ - 1}}\left( . \right)$ represents inverse Fourier transform.

The distribution of reflectivity of targets is denoted as f(x, y). The light reflected by the targets travels through the water and is received by the receiving antenna. The bucket detector records the energy of focused light. The light in the optical path from the target to the optical receiving antenna is also be absorbed and scattered by the water, which is similar to the emitting optical path. However, the bucket detector has no spatial resolution and is not used for imaging, which leads to that the distortion of water forward scattering has no effect on the imaging quality. In the receiving optical path, only the effect of water on the attenuation of echo energy is considered. The total energy attenuation coefficient caused by the receiving optical path is denoted as α₂. Therefore, the measurement of bucket detector is:

$\begin{split} {I_\phi }\left( {{f_x},{f_y}} \right) =& n + \displaystyle\iint {{\alpha _1}{\alpha _2}{E_0}f\left( {x,y} \right)}\;\cdot \\ &{{P_\phi }\left( {x,y;{f_x},{f_y}} \right)} \;\cdot\\ &{ * {{\mathbb F}^{ - 1}}\left[ {H\left( {{f_x},{f_y}} \right)} \right]{\rm{d}} x{\rm{d}} y} \end{split} $ (5)

In FSPI system, the target spatial spectrum can be obtained by the four-step phase-shifting approach, where all four speckle patterns have the same spatial frequency (f_x, f_y), but different phases. The phase shift between the two adjacent patterns is a constant $\pi $/2. When the values of the phase ϕ are equal to 0, $\pi $/2, $\pi $ and 3 $\pi $/2, the measurements of bucket detector are denoted as I₀, $I_{{\pi}/2} $, $I_{\pi} $ and $I_{3{{\pi}/2}} $respectively. The four bucket detector measurements are utilized to estimate the Fourier spectrum value of the targets at the spatial frequency (f_x, f_y), which can be expressed as:

$C\left( {{f_x},{f_y}} \right) = \left( {{I_{{{3\pi } / 2}}} - {I_{{\pi / 2}}}} \right) + j\left( {{I_\pi } - {I_0}} \right)$ (6)

By submitting (1) and (5) into (6), (6) can be represented as:

$C\left( {{f_x},{f_y}} \right) = {\alpha _1}{\alpha _2}{E_0}H\left( {{f_x},{f_y}} \right)C'\left( {{f_x},{f_y}} \right)$ (7)

$ \begin{split} &{C'\left( {{f_x},{f_y}} \right) = \displaystyle\iint {f\left( {x,y} \right)}}\cdot \\ &{\exp \left( { - j2\pi \left( {{f_x}x + {f_y}y} \right)} \right){\rm{d}} x{\rm{d}} y} \end{split} $ (8)

C'(f_x, f_y) represents the true spatial frequency of unknown targets at the position (f_x, f_y) of spatial spectrum. Under the influence of water back scattering and forward scattering, the obtained spatial frequency value is equal to the product of the true spatial frequency value of unknown targets, attenuation coefficient and water degradation function. Moreover, the influence of water degradation function in spatial domain is mainly reflected in the blurring and unclear image. Therefore, the effects of the water degradation function need to be eliminated to obtain the clear image.

Due to the relative complexity of water, it is rather difficult to measure or estimate the parameters k and s in the water degradation function. The way to decrease the influence of water forward scattering is equivalent to an ill-posed problem of blind image restoration. In order to solve the ill-posed problem, the water degradation function compensation method is implemented to deburr the reconstruction image of FSPI in the spectrum domain. By submitting (3) into (7), (7) can be simplified as:

$ \begin{split} C\left( {{f_x},{f_y}} \right) =& {\alpha _1}{\alpha _2}{E_0}C'\left( {{f_x},{f_y}} \right)\cdot\\ &{\exp \left[ { - k{{\left( {{f_x}^2 + {f_y}^2} \right)}^s}} \right]} \end{split}$ (9)

Suppose σ is the maximum value of spatial spectrum, which can be obtained by the FSPI system. The left and right sides of (9) are divided by σ for normalization. Ĉ (f_x, f_y) and Ĉ'(f_x, f_y) are used to represent the normalized spatial spectrum of the blurred image and the true image respectively. Taking logarithmic operation on the left and right sides of (9), it can be expressed as:

$\begin{split} \ln \hat C\left( {{f_x},{f_y}} \right) =& \ln {\alpha _1}{\alpha _2}{E_0} + \ln \hat C'\left( {{f_x},{f_y}} \right) -\\ &{ k{{\left( {{f_x}^2 + {f_y}^2} \right)}^s}} \end{split} $ (10)

According to the symmetry of Gaussian-like functions, any line passing through the origin in the spatial spectrum satisfies (10). It is reasonable to set the value of f_y as 0 and represent lnα₁α₂E₀+lnĈ'(f_x, f_y) with a constant A. (10) is reformulated as:

$\ln \hat C\left( {{f_x},0} \right) = A - k{\left( {{f_x}^2} \right)^s}$ (11)

Based on the spatial spectrum obtained by the FSPI system, the unknown parameters k and s can be calculated by the least squares fitting. The approximate water degradation function H (f_x, f_y) can be obtained when the parameters k and s are known. The spatial spectrum is divided by the water degradation function H (f_x, f_y) to calculate the true spatial spectrum of targets. Combining the inverse Fourier transform, the image of underwater targets can be finally reconstructed.

In order to further prove the effectiveness of the proposed method, experiments of the underwater FSPI system are carried out. The experiment system is shown in Fig.2. Figure 2(a) is the photo of integrated experiment system. In the experiment system, the wavelength of laser is 532±1 nm and a DMD is used to modulate the spatial distribution of laser. The DMD is an ultra-fast V-7001 of ViALUX. The modulation frequency of the DMD is set as 100 Hz. Because DMD can only perform binary modulation, the Floyd-Steinberg error diffusion dithering method is used to convert the sinusoidal structured patterns from grayscale to binary. In the Floyd-Steinberg error diffusion dithering method, the bicubic interpolation coefficient of the upsampling process is set to be 4. Therefore, in order to achieve the reconstruction results with the resolution of 256×256, the effective active pixel number of the DMD is 1 024×1 024. The detector is chosen as H11523-110-NF and it is connected to the data acquisition card, whose model is PCIe-9852. The digital reflected light signal is sent to the computational system to reconstruct the image of targets. The CPU of computational system is Intel(R) Core(TM) i5-8500 CPU with the work frequency of 3 GHz and the memory of 8 GB. Figure 2(b) and (c) are the target and detection environment of the experiment respectively. As Fig.2(b) shows, the satellite model is chosen as the target whose material is ABS plastic. The length and diameter of the model is 150 mm and 52 mm. The satellite model is placed in a sink, whose size is 95 cm×36 cm×40 cm. The sink is filled with pure water and the target is flooded by water. Because the FSPI is a kind of active illumination imaging system, in order to avoid interference from ambient background light and simulate the dark environment in deep water detection as well, the experiment is performed in the dark environment.

In order to verify that the proposed method can realize underwater imaging, the target is imaged with the proposed method in the pure water. Considering the sparsity of the target in the spatial spectrum domain, on the premise of ensuring the quality of imaging reconstruction, the down-sampling measurement can be taken when the FSPI imaging scheme is used. The down-sampling measurement can reduce the number of measurement and improve the imaging efficiency of FSPI system. The experiment results of pure water environment is shown in Fig.3. In Fig.3, different rows represent the reconstruction results with different sample rate, which are chosen as 100%, 50%, 25% and 12.5% respectively. The reconstruction results in the first column are based on the conventional FSPI and the reconstruction results in the second column are based on the FSPI with the proposed method. It can be seen from Fig.3 that the sample rate can influence the reconstructed results of the FSPI. When the sample rate is high, the reconstructed result of the FSPI has high reconstruction quality and the spatial resolution is great. The details and the edge information of the target can be completely reconstructed. When the sample rate decreases, the quality of reconstructed image turns bad. The reconstructed image will have low spatial resolution and get blurring. Moreover, there is ringing effect appearing in the reconstructed results. For instance, the lower half of the target surface has a screw-like structure. The reconstructed results with high sample rate have clearly screw-like structure, which can hardly be reconstructed with low sample rate. In addition, although the laser has a high permeability in the pure water, scattering still influences the quality of the FSPI results. Compared with the conventional FSPI results, the results of the proposed method have better visual effect. The details and edge of the conventional FSPI results are more blurring. Especially, under the high sample rate condition, the screw-like structure parts of reconstructed results with the proposed method are much clearer than those with conventional FSPI. It can be concluded that the blurring caused by the water forward scattering can be eliminated in certain degree with the proposed method.

Furthermore, the turbidity of water is increased to verify that the proposed method not only has the ability to eliminate the image blur caused by forward scattering in the pure water, but also is effective in turbid water. The way to increase the turbidity is to add milk to the water. In the experiment, the amount of milk added to water is 5 mL, 15 mL and 20 mL. Meanwhile, the target in the turbid water is imaged through differential single pixel imaging (DSPI), conventional FSPI and the proposed method respectively. And the numbers of measurements with different reconstruction methods are all 32768 corresponding to 50% sample rate. The experiment results are shown in Fig.4. In Fig.4, different rows represent the reconstructed results under different water turbidity condition. The rows from the top to the bottom are 5 mL, 15 mL and 20 mL respectively. Different columns represent the different reconstruction method, and the columns from the left to the right are DSPI, FSPI and the proposed method.

It can be seen from the Fig.4 that the turbidity of water has a great influence on the reconstructed results. The more turbidity the water is, the worse quality of reconstructed results with different method is. The DSPI method has the worst reconstruction quality among the three methods. Under each condition, all the reconstructed results of DSPI contain a lot of noise. The DSPI method cannot avoid blurring caused by water forward scattering. The turbidity of water has the greatest influence on the DSPI results. When the amount of added milk is 20 mL, the DSPI method can hardly reconstruct the details and edge profile information of the target. The FSPI has a higher reconstruction ability than the DSPI. The reconstruction results of FSPI have less imaging noise than those of DSPI. However, the reconstruction results of FSPI are still affected by the turbidity of the water seriously. With the increasing turbidity of water, the forward scattering of water increases gradually, and the results of the FSPI gradually become blurring. From the point of visual effect, the surface of the reconstructed image is like to be covered with a layer of thin gauze. The visual effect of the reconstructed image is blurred, which is mainly caused by the impurities of water. Due to the scheme of the proposed method which is based on FSPI, the proposed method has similar ability to the conventional FSPI in the aspect of anti-image noise. Meanwhile, there is a step in the proposed method to estimate the forward scattering of water. The influence of water degradation function caused by forward scattering is removed partially in the frequency domain. As a result, it can be clearly seen that the proposed method has a clearer reconstruction result compared with the previous two methods. Especially, the proposed method does a better job of detail recovery and edge preserving than the conventional FSPI. When the amount of added milk reaches 10 mL, the reconstruction results of the conventional FSPI have already become blurred. But the proposed method can still reconstruct the target clearly and the details and edges of the result can be well presented.

In order to more intuitively show the high-efficiency reconstruction of the proposed method, the gray values in the line 163 (the light blue line) of the reconstruction results with various methods are drawn into a two-dimensional curve, as shown in the Fig.5. The black curve in Fig.5 is the gray values in the line 163 of the original target image. The green one is the two-dimensional curve of DSPI, the blue one is the two-dimensional curve of the conventional FSPI, and the red one is the two-dimensional curve of the proposed method. It can be seen from Fig.5 that the reconstruction error of DSPI is relatively large, and the noise almost drowns out the useful signals. There is less noise in the reconstructed results of the conventional FSPI and the proposed method. However, the turbid water decrease the reconstruction quality of the two methods, which is reflected in the blurred edge of the reconstruction results. The curves in the area of the oval dotted line in Fig.5 are corresponding to the area of target edge. The target edge areas are enlarged and shown in the right side of Fig.5. Because the target edge area of original is the clearest, the curves of gray value in the original target image are the steepest. For the influence of water forward scattering, the reconstructed results of the conventional FSPI and the proposed method get blurry. Thus, the gradients of the gray value with the conventional FSPI and the proposed method at the edge are both smaller than that of the original target image, and the curves (the red and blue ones) are gentle. With the increase of water turbidity, the water forward scattering gets more serious and the curves get gentler. Because the proposed method compensates for the influence of forward scattering, the curve at the edge is less steep than the conventional FSPI. Therefore, the red curve at the edge is steeper than the blue one.

In order to quantitatively compare the influence of water turbidity on the quality of reconstruction results of various methods, the point spread function (PSF) of reconstructed results are calculated based on the slant edge method^[27]. The process of calculating the PSF is shown in Fig.6. First of all, a region which contains the desired edge is selected, assuming that the selected region is circled out with the red box in the first image of Fig.6. Sobel operator is used for edge detection in the selected region. The equation of the line corresponding to the edge can be fitted according to the edge. The distances between the whole pixels in the selected region and the fitted line are calculated. The distances and gray values of each pixel are stored into the set of X and Y. According to the data in the set of X and Y, the edge spread function (ESP) can be obtained and the Fermi function is used for ESP fitting. The reason to adopt Fermi function is that it can be well approximate the ESP, and it can be differentiable everywhere, which is conducive to the differentiation operation when finding the line spread function (LSP) later. In order to fit the ESP more accurately, a summation of three Fermi functions is adopted for fitting in practice.

The summation of three Fermi functions is shown as:

${\rm{F}} \left( x \right) = K + \sum\limits_{i = 1}^3 {\dfrac{{{a_i}}}{{\exp \left( {\dfrac{{x - {b_i}}}{{{c_i}}}} \right) + 1}}} $ (12)

where K, a_i, b_i and c_i are the fitting coefficients, which can be calculated according to the data in set X and Y. Considering the influence of noise on the fitting results and the accuracy improvement of the fitting process, especially the image noise is serious in the DSPI, the set of X and Y can be expanded by the mean of cubic spline interpolation. The accuracy of the fitting can be improved by increasing the data set. Then, according to the ESF, the LSF is solved by the differential. After differentiating (12), it can be obtained:

${\rm{LSF}} \left( x \right) = \dfrac{{{\rm{d}} {\rm{F}} \left( x \right)}}{{{\rm{d}} x}} = \sum\limits_{i = 1}^3 {\dfrac{{ - {a_i}\exp \left( {\dfrac{{x - {b_i}}}{{{c_i}}}} \right)}}{{{c_i}{{\left[ {\exp \left( {\dfrac{{x - {b_i}}}{{{c_i}}}} \right) + 1} \right]}^2}}}} $ (13)

Because of the isotropic distribution, the PSF can be expressed as:

${\rm{PSF}} \left( x \right) = {\rm{LSF}} \left( x \right) \times {\rm{LSF}} \left( y \right)$ (14)

In order to guarantee the fairness of PSF calculation, the calculation region of each method should be the same when calculating the PSF with slant edge method. The calculated PSF of reconstructed results with different methods are shown in Fig.7.

The different columns in Fig.7 represent the results with different reconstruction methods. The used methods from the left to the right are DSPI, FSPI and the proposed method respectively. Different rows represent the different imaging environments, and the added milk amounts are 20 mL, 10 mL and 5 mL respectively. As the calculated PSF shows in the Fig.7, the turbidity of water affects the reconstruction results, which is embodied in the broadening of the point diffusion function of the reconstruction results. When the amount of added milk is relatively small, the turbidity of water is relatively low. The full width at half maximum (FWHM) of the PSF reconstructed by various methods are relatively narrow. With the increase of the amount of added milk, the water gradually becomes turbid. All the FWHMs of the PSF with various methods become wider, and the imaging results become blurry. On the other hand, in the case of the same turbidity of water, the PSF of the reconstruction results with the DSPI is the widest, and the image gets the most blurry and the resolution is relatively low. The PSF of the reconstruction results with the conventional FSPI is slightly better than that of the DSPI, especially when the water is relatively turbid. For instance, when 20 mL milk is added, the FWHM of the PSF of FSPI is significantly smaller than that of the DSPI. The reconstructed result is clearer than that of the DSPI. The proposed method has the best imaging effect among the three methods. The FWHM of the PSF of the reconstructed results is the narrowest of all under the same turbidity degree condition. The resolution of the proposed method is the highest. Therefore, it can be concluded that forward scattering causes the least blurring to the proposed method.

In addition, the FWHM of PSF in the Fig.7 is calculated quantitatively. In Fig.7, the total number of pixels in the x direction is equaled to the total number of pixels in the y direction. Hence, this pixel number is denoted as N. For some certain PSF figure, the pixel value of different pixel coordinates (x, y) is P. Supposing that the maximum value of the PSF figure is P_max and the corresponding coordinate is (x_max, y_max). Therefore, all the point on the line of y= y_max, which satisfy the condition of P>P_max/2, can be found. The number of these points can be expressed as N_x. According to the isotropic distribution of PSF, the FWHM of PSF is defined as:

$FWHM = \dfrac{{{N_x}}}{N}$ (15)

Based on (15), the FWHM of PSF with various methods under different turbidity conditions can be calculated. The calculated results is shown in Tab.1. The same conclusion can be drawn from Tab.1 as in Fig.7. From the quantitative calculation results, it can be obviously seen that the proposed method has high quality imaging performance than the other methods.

In order to reduce the influence of forward scattering of turbid water on SPI, an FSPI based on the water degradation function compensation method is proposed. In this scheme, the measurements of the FSPI are processed in the spatial spectrum domain, and the water degradation function is estimated. Utilizing the estimated water degradation function, the eliminating forward scattering is realized in the spatial spectrum domain. The experiment is also carried out based on the proposed method. In the experiment, different amounts of milk are added in the pure water, and the target image is reconstructed by the DSPI, conventional FSPI and our method respectively. Through comparison, it can be found that the proposed method can greatly improve the quality of reconstructed results. In addition, the PSFs of the reconstructed results with various methods are quantitatively calculated based on slant edge method. Comparing the FWHM of PSF with different methods, it can be seen that the reconstructed results of the proposed method have the narrowest PSF. Thus, the image is the clearest. The quantitative calculation is consistent with the qualitative visual effect. It is proved that the proposed method can effectively reduce the influence of forward scattering of water and obtain a clearest reconstructed image.