Performances of Polarization-Retrieve Imaging in Stratiﬁed Dispersion Media

: We constructed an active imaging model within 10 km of the atmosphere from the satellite to the ground based on Monte Carlo (MC) algorithm, and, because of the inhomogeneous distributions of the scattering particles in atmosphere environment, 10 km atmosphere layer was divided into ten layers in our model. The MC algorithm was used to simulate the transmission process of photons through the atmosphere. By launching lasers of linear polarization states from satellites to ground, the intensity, degree of polarization (DoP), polarization di ﬀ erence (PD), and polarization retrieve (PR) images can be obtained. The contrast of the image, peak signal to noise ratio (PSNR), and structural similarity index (SSI) were used to evaluate the imaging quality. The simulated results demonstrate that the contrast of images is degraded as the atmosphere becomes worse. However, PR imaging have a better contrast and better visibility in di ﬀ erent atmospheric conditions. Meanwhile, we found that Mueller matrix (MM) can retrieve the original images very well in a certain range of atmospheric conditions. Finally, the simulation also shows that di ﬀ erent wavelengths of light sources have di ﬀ erent penetration characteristics, and, in general, infrared light shows better performances than visible light for imaging.


Introduction
Imaging quality will be degraded for two main reasons [1]: (1) strong scattering and absorption of light by the disturbance, resulting in object information loss; and (2) noisy photons scattered in the turbid medium without the target's information. Therefore, it is important to reduce the scattering and absorption effects to improve the imaging quality for various applications, such as underwater imaging [2], remote sensing [3], haze removal [4][5][6], and biological detection [7,8]. On the one hand, it has been proposed to numerically remove the negative effect based on computer vision methods [9,10]. On the other hand, the physical methods to reconstruct image could achieve more real target's information. For the traditional imaging method based on intensity, it is difficult to have a clear image when the propagation distance is large due to the energy loss.
As an intrinsic characteristics of electromagnetic wave, polarization has attracted tremendous research interests over the past decades. Polarization has also been proven to be a powerful tool for characterizing the target. For example, Mie ellipsometry [11,12] and full radiative transfer simulations [13,14] have been used to investigate the property of the information channel composed of dust particles and gaseous medium. Polarization imaging, in particular, can be employed to obtain not only irradiance information of the object but also the polarization information, therefore improving the ability to detect and reconstruct the object information from a complex scene. It has been demonstrated that polarization imaging is an effective way to improve qualities for underwater imaging [15] and dehazing imaging [16]. It is certain that the polarization information will degenerate during light scattering, in which the image quality will degrade as well. To solve this problem, a polarization difference (PD) technology can be utilized to remove the depolarized light [17]. However, when the absorption is strong, PD technology will not be applicable. In 2001, a scattering model was proposed by Schechner, where the backward scattering light is removed by estimating the background value, thus obtaining the irradiance of the original image [18]. On this basis, they further presented an illumination compensation method to overcome degradation effects occurring in weak irradiance underwater imaging [19]. Dubreuil et al. also proposed an optical correlation-based method to improve underwater target detection [20]. In addition, polarimetric dehazing method has also been introduced to enhance the visibility of hazy images based on the orientation angle information from the Mueller matrix (MM), the angle of polarization distribution analysis, and visible and infrared image fusion technology [21][22][23][24]. Furthermore, Hu et al. achieved high quality imaging under the condition of nonuniform optical field by estimating the light intensity of background at different positions of image [25]. In simple words, polarization imaging has already shown the great ability of increasing the imaging quality through scattering medium. However, all the above methods focus on removing the scattering light. In contrast, based on the Monte Carlo (MC) algorithm for obtaining the Mueller matrix of the transmission medium [26][27][28][29][30][31], the polarization retrieve (PR) method [32,33] has been proposed for investigating the transmission characteristics of the polarization information in different environments, such as poly-dispersion system [34], layered atmosphere [35,36], and underwater [37,38]. This method can retrieve the original target information and obtain much clearer images, which demonstrates its potential applications in a complex scattering system with long transmission distance. For the case of atmosphere remote sensing imaging, most polarization imaging methods may lose their effectiveness due to the complexity of scattering process and the long transmission distance. Therefore, it is an urgent requirement to realize polarization imaging of target on the ground through the atmosphere.
In this paper, we establish a PR imaging model based on stratifying theory of atmosphere by MC algorithm. In our simulations, the whole of atmospheric environment was set to be 10 km. We utilize the active polarization imaging technique to reconstruct the target image on the ground. Compared with traditional imaging methods, we can find that the contrast of the image reconstructed by PR method has been enhanced significantly. This paper is organized as follows. In Section 2, we briefly describe the PR imaging model for stratified atmosphere by MC algorithm. In Section 3, we present the simulations of our proposed method and compare our results with the conventional imaging methods, including intensity imaging, degree of polarization (DoP) imaging, and PD imaging. We have also investigate the performance of reconstructing images for the objects with different geometric complexities. In addition, we also discuss the robustness of the initial MM to the variation of particle sizes in first layer (0-2 km) of atmosphere. Finally, we investigate the effect of different wavelengths on the imaging quality. Section 4 summarizes this work.

PR Method
The polarization states of any light can be represented by the four-component Stokes vector as: where E x and E y are the amplitude of xand y-components of the electric field, respectively; I represents the total intensity; Q is the intensity difference between xand y-polarized components; U indicates the intensity difference between 45 • and 135 • linear polarizations; and V is the intensity difference between left and right circular polarizations. The following points are worth mentioning: (1) only I remains for completely unpolarized light; (2) V disappears for completely linear polarizations; and (3) When light transmit through a scattering medium, the relationship between the incident and transmitted Stokes vectors can be expressed by Equation (2): where S out and S in represent the Stokes vectors of transmitted and incident light, respectively. M is the MM expressing interacting effects of the medium system [32]. It has been widely demonstrated that the MM can describe the optical properties of the medium completely. In our simulations, the effective Mueller matrix (EMM) could be obtained by launching six different polarization states of photons as the incident lights: (1) horizontal polarization; (2) vertical polarization; (3) 45 • polarization; (4) 135 • polarization; (5) right circular polarization; and (6) left circular polarization. The corresponding effective Mueller matrix (EMM) for the single incident light can be expressed as [32,39]: where (I k , Q k , U k , V k ) with k = 1-6 is the Stokes vectors of the transmitted light from the media system under the kth (k = 1-6) incidence, such as the fourth incidence of (1, 0, −1, 0). However, in our imaging simulations, we used a plane source as the incidence, and the MC simulation need a very huge number of photons to gain reliable results, which results in a very low computational efficiency. To improve the computational efficiency, we introduce a method for obtaining the EMM for the system of incident plane source by shifting position and superposition principle in the MC model [38]. It is a typical linear problem for light transmission in the medium. Thus, the EMM of the plane incident light system can be defined as the superposition of the effective Mueller matrix (SEMM), which can be expressed as: where the M S (x,y) is the SEMM and M (i,j) (x,y) means the EMM in the position of (i, j); (x,y) expresses the two-dimensional plane of the planar light source, which is composed of the point light source array and the incident position of the point light source can be labeled as (i, j) [30,35]; and K and N are the maximum value for the positions of the incident plane source in the directions of x and y, respectively. Therefore, once we accurately obtain the SEMM of the medium in advance, the original polarization information can be reconstructed by In Equation (5), we define the original information S in recovered by S out and M S as the PR method. On the other hand, the DoP of the reconstructed image can be calculated based on PR method as Remote Sens. 2020, 12, 2895 4 of 16 DoP = Q 2 r + U 2 r + V 2 r /I r , where I r , Q r , U r and V r are the four elements of the retrieved Stokes vector. The DoP values are 0 and 1 demonstrating unpolarized and completely polarized lights, respectively. The polarization difference (PD) imaging can be expressed as PD = |I co − I cross |/|I co + I cross | in which Ico and Icross are the received optical intensities by using analyzers with orientation parallel and perpendicular to the illuminating polarization, respectively [40,41].

Simulation Method
The schematic of the system model is shown in Figure 1a. Actually, the atmosphere is inhomogeneous along vertical direction. To simulate the real case, 10-km atmosphere is evenly divided into 10 layers, and each layer can be thought of as a uniformly dispersed system. The more layers you divide, the more accurate results you can obtain in theory. For simplicity, only ten layers are divided here. In fact, more layers would burden the simulation calculation. In addition, the size of whole scattering system is modeled to be 10 km × 10 km × 10 km. A linearly polarized illumination source (S = (1,1,0,0) T ) with central wavelength of 1536 nm is normally launched from a satellite to the ground, and a detector is placed in the same side with the light source, where both sizes are 1 × 1 m 2 . In the simulations, we set the ground as a scatter without total absorption to make most of the photons that reach the ground able to be reflected back. The target is placed in the center of the ground facing the light source; the size is also set to 1 × 1 m 2 . The scattering characteristics of the target can be represented by the corresponding MM. In the preliminary simulation, we found that the energy received by the detector is about 3-5% of the energy emitted in different weather conditions. Considering that scattered light can carry enough target information, we emitted 1 × 10 8 photons. Finally, the reflected light can be received on the satellite. layers you divide, the more accurate results you can obtain in theory. For simplicity, only ten layers are divided here. In fact, more layers would burden the simulation calculation. In addition, the size of whole scattering system is modeled to be 10 km × 10 km × 10 km. A linearly polarized illumination source (S = (1,1,0,0) T ) with central wavelength of 1536 nm is normally launched from a satellite to the ground, and a detector is placed in the same side with the light source, where both sizes are 1 × 1 m 2 .
In the simulations, we set the ground as a scatter without total absorption to make most of the photons that reach the ground able to be reflected back. The target is placed in the center of the ground facing the light source; the size is also set to 1 × 1 m 2 . The scattering characteristics of the target can be represented by the corresponding MM. In the preliminary simulation, we found that the energy received by the detector is about 3-5% of the energy emitted in different weather conditions. Considering that scattered light can carry enough target information, we emitted 1 × 10 8 photons. Finally, the reflected light can be received on the satellite.  Figure 1b shows the imaging object with three different sub-objects. The first part is steel sheets with different geometric shapes, including the letters of "HFUT", figure "-1945-", and a building model. They have high reflectivity and low depolarization characteristics. The second part is made up of marble with both high reflectivity and high depolarization characteristics. The third part is the bottom wood, which has high depolarization characteristics but low reflectivity. The corresponding Mueller matrices of steel, marble, and wood are summarized in Table 1 [42]. When the illumination beam is aligned with the line of sight or the objects are purely depolarizing, the Mueller matrix of most materials will be diagonal [43][44][45]. It was suggested that the Muller matrices for steel, marble, and wood are independent of the wavelength. In our simulations, the scattering medium system is modeled with randomly distributed absorbing particles based on Mie scattering theory, whose size and dispersion properties are taken into consideration. In a real atmospheric model, it should be noted that the molecular scattering and absorption are not taken into account [46]. We mainly consider the influence of aerosol particles which contains 70% water soluble and 30% dust, whose refractive indices are 1.51-0.023i and 1.4-  Figure 1b shows the imaging object with three different sub-objects. The first part is steel sheets with different geometric shapes, including the letters of "HFUT", figure "-1945-", and a building model. They have high reflectivity and low depolarization characteristics. The second part is made up of marble with both high reflectivity and high depolarization characteristics. The third part is the bottom wood, which has high depolarization characteristics but low reflectivity. The corresponding Mueller matrices of steel, marble, and wood are summarized in Table 1 [42]. When the illumination beam is aligned with the line of sight or the objects are purely depolarizing, the Mueller matrix of most materials will be diagonal [43][44][45]. It was suggested that the Muller matrices for steel, marble, and wood are independent of the wavelength. In our simulations, the scattering medium system is modeled with randomly distributed absorbing particles based on Mie scattering theory, whose size and dispersion properties are taken into consideration. In a real atmospheric model, it should be noted that the molecular scattering and absorption are not taken into account [46]. We mainly consider the influence of aerosol particles which contains 70% water soluble and 30% dust, whose refractive indices are 1.51-0.023i and 1.4-0.008i, respectively, at wavelength of λ = 1536 nm [47]. Thus, their reference mixture ratio is 7/3. Correspondingly, the collision probability between the incident photons and water-soluble particles is 0.7, and the collision probability of between the incident photons and dust particles is 0.3. The relative refractive index of air is 1. For simplicity in calculations, those particles are roughly equivalent to spherical particles. In addition, the increase in height is accompanied by decrease in both number density and size of particles. Table 2 summarizes the parameters of the clear/cloudy atmosphere used in this paper, in which PND is particle number density. Intuitively, PND decreases with the increase of altitude. It should be described that, when the height is above 5 km, the parameters of clear and cloudy weather are nearly the same. When the altitude is below 5 km, in each layer, for the clear and cloudy weather, the PNDs will be different, which will make the scattering coefficient and absorption coefficient be different. Thus, there are two values for each layer below 5 km in Table 2. r eff represents the equivalent radius of particles. The scattering coefficient and absorbing coefficient of the medium can be calculated in each layer by Mie scattering theory. All of these parameters are obtained from the works in [47][48][49]. Note that we choose a mid-latitude location of N31 • , E117 • (China) as an example. As is known, most previous works simply focus on homogeneous media or only discuss the medium, but here the model is also applied in the real free space system. Thus, a modified MC program is developed to apply in stratified dispersion media. When light travels through a homogeneous medium, its mean free path (mfp) can be calculated as [26,36]: where is ξ a pseudo-random number between 0 and 1 and µ e is the extinction coefficient of the medium equaling to the sum of scattering coefficient and absorption coefficient. The change of photon energy during scattering process can be characterized by Lambert-Beer law as: where t is the transmittance and I s and I 0 are the intensities of scattered light and incident light, respectively. By comparing Equations (6) and (7), we conclude that t is numerically equal to ξ. τ is the Remote Sens. 2020, 12, 2895 6 of 16 optical thicknesses, which is the product of the extinction coefficient and the transmission distance.
Once that the random number ξ is chosen, the transmittance can be determined as well. As shown in Figure 2, the total energy loss is equal to the sum of that in each layer. Obviously, t is equal to the product of the transmittances of corresponding sublayers and can be expressed as: Remote Sens. 2020, 12, x 6 of 16 Combining this result with Equation (7), we can get where q, j, and s are distance parameters describing the position parameters of the next colliding particles. Then, we need to obtain the position of the photon. As shown in Figure 2, the single photon is scattered at position p0 and the next scattered position is p1. We suppose that the coordinate of p0 is (x0, y0, z0), and the propagation direction is (cosα, cosβ, cosγ), where α, β, and γ are the angles between u and the x, y, and z axes respectively. The extinction coefficient of each layer can be calculated by the basic parameters of the medium. By generating a random number, a reference distance can be obtained by Equation (6). Thus, when the initial position and direction were determined, the iterative algorithm can be utilized to solve the photon pathlength by Equation (9). Then, we consider three situations for understanding the calculation fully.
A. cosγ > 0. The photon is downlink scattering. j and s are If wn+1 ≤ j, the point p1 locates in the (n + 1)th layer. The next scattering still happens in this sublayer. Clearly, the photon free path is, Otherwise, the next scattering happens in another layer. In conjunction with Equations (9) and (10), q holds: If s > q, m adds 1, and then continue with the next iteration. Else the iteration is over, and the pathlength is B. cosγ < 0. The photon is uplink scattering. Then, Equation (9) should be rewritten as Combining this result with Equation (7), we can get where q, j, and s are distance parameters describing the position parameters of the next colliding particles. Then, we need to obtain the position of the photon. As shown in Figure 2, the single photon is scattered at position p 0 and the next scattered position is p 1 . We suppose that the coordinate of p 0 is (x 0 , y 0 , z 0 ), and the propagation direction → u is (cosα, cosβ, cosγ), where α, β, and γ are the angles between u and the x, y, and z axes respectively. The extinction coefficient of each layer can be calculated by the basic parameters of the medium. By generating a random number, a reference distance can be obtained by Equation (6). Thus, when the initial position and direction were determined, the iterative algorithm can be utilized to solve the photon pathlength by Equation (9). Then, we consider three situations for understanding the calculation fully.
A. cosγ > 0. The photon is downlink scattering. j and s are If w n+1 ≤ j, the point p 1 locates in the (n + 1)th layer. The next scattering still happens in this sublayer. Clearly, the photon free path is, Otherwise, the next scattering happens in another layer. In conjunction with Equations (9) and (10), q holds: Remote Sens. 2020, 12, 2895 7 of 16 If s > q, m adds 1, and then continue with the next iteration. Else the iteration is over, and the pathlength is B. cosγ < 0. The photon is uplink scattering. Then, Equation (9) should be rewritten as In this case, j and s are Finally, the remaining calculation process of pathlength has the same form as Equation (11).
The propagation direction is perpendicular to the z-axis. Thereby, the photon free path has the same form as Equation (9).
When the scattering position is determined, the next scattering direction of the photon is determined by the scattering phase function. For the incident light with a Stokes vector S 0 = [I 0 , Q 0 , U 0 , V 0 ] T , the phase function P(θ, ϕ) can be expressed as where s 11 (θ) and s 12 (θ) are two elements of particle's single scattering matrix and θ and ϕ represent scattering angle and azimuth angle, respectively [50]. Then, we calculate the MM of the scattering system, and, based on the obtained MM, the PR method has been used for reconstructing the target's images. In general, the MM of a system is an intrinsic property and independent on incident light source. However, we know that, in the atmosphere, the uplink and downlink of light are not reversible processes, and the corresponding MM of uplink and downlink are also different. Thus, to obtain the MM for whole of the optical link, we need calculate two MMs for downlink (from satellite to ground) and uplink (from ground to satellite), respectively. In fact, we can send a laser from satellite to ground for determining the downlink MM. To reconstruct the image of the target, we need to get the reflected light from the target surface from the receiving light. Based on PR method, the corresponding MM is the MM of uplink. We placed a laser with the same wavelength of 1536 nm near the target, and light travels perpendicular to the ground. Here, we choose the incident light source as the point light source, in which, to ensure the accuracy of simulation, 1 × 10 8 photons are emitted. The basic parameters of the scattering medium are demonstrated in Table 2, from which we can obtain the EMM. We know that the transmission of light in the scattering system is a typical linear problem. Because of the symmetric and linear characteristics, the superposition principle is consistent. By utilizing EMM of the point light source, we can then use the shifting position and superposition principle to calculate the SEMM of the plane light source for polarization imaging [38]. Figure 3 shows the calculated Mueller matrices of the proposed atmosphere model in clear atmosphere by using MC method, in which all elements are normalized to the m 11 element to compensate for the radial decay of intensity. Each of the images is displayed with its own color map to enhance the image contrast. Generally speaking, m 22 and m 33 show the depolarization characteristics of the scattering system to the linearly polarized light, and m 44 shows the depolarization characteristics of the scattering system to the circularly polarized light. We can see that some of the elements in the matrix actually have some symmetry, and the patterns are identical to those works reported previously [51,52], which are correct and valid to some extents.
to enhance the image contrast. Generally speaking, m22 and m33 show the depolarization characteristics of the scattering system to the linearly polarized light, and m44 shows the depolarization characteristics of the scattering system to the circularly polarized light. We can see that some of the elements in the matrix actually have some symmetry, and the patterns are identical to those works reported previously [51,52], which are correct and valid to some extents. Based on the layered scattering theory and MC algorithm, we used the traditional method and the PR method to image in the atmospheric environment. The proposed scheme is described in Figure  4, which contains three steps: Based on the layered scattering theory and MC algorithm, we used the traditional method and the PR method to image in the atmospheric environment. The proposed scheme is described in Figure 4, which contains three steps: Remote Sens. 2020, 12, x 8 of 16 Step 1: Traditional polarization imaging methods. In the first step, based on the layered scattering theory and the MC method, we can obtain the Stokes vectors (Sout) of output light. Therefore, the DoP = + + / , where Iout, Qout, Uout, and Vout are the four Stokes elements of the Sout. The polarization difference (PD) imaging can be concluded with Ico and Icross.
Step 2: Imaging based on PR method. The second step aims to retrieve the original image with higher quality. Indeed, the traditional imaging method loses its effect when the cloudy medium is very heavy. The PR imaging method can complete retrieve the original image in theory. Then, we compare the imaging quality under different methods through the corresponding image-quality indicators: peak signal to noise ratio (PSNR) and structural similarity index (SSI).
Step 3: Discussion. Finally, in the third step, to better analyze the PR imaging model in atmospheric scattering medium, we consider the following points: (1) effects of different scattering environments on imaging quality; (2) the penetration characteristics of light at different wavelengths in the atmosphere; and (3) MM's applicability with changing scattering particle sizes.

Results and Discussions
We simulated the whole imaging process, from which we could obtain the corresponding images by the relative imaging methods. From Table 2, we can get the scattering coefficient and absorption coefficient under the clear atmosphere and calculate the extinction coefficient in each layer, and then multiply the extinction coefficient by the corresponding transmission distance to get the optical thickness of 1τ = 1.8. The light intensity transmittance (T1) in this environment can be calculated from Equation (7) as 16.63%. Figure 5 shows the imaging results from intensity (I components of Stokes vector), DoP, PD, and PR methods through scattering medium with different concentrations. From left to right, the optical thicknesses, which is the product of the extinction coefficient and the transmission distance, were set as 1τ, 3τ, 6τ, and 9τ, respectively, by changing the extinction coefficients (changing the number density and radius of particles), which demonstrates the turbidity of the medium systems. In the four different scenarios in Figure 5, the clear weather is set as a reference standard. In clear weather, the scattering coefficient and absorption coefficient of 1τ are increased by three, six, and nine times to obtain other three different scenes. 3τ, 6τ, and 9τ equal Step 1: Traditional polarization imaging methods. In the first step, based on the layered scattering theory and the MC method, we can obtain the Stokes vectors (S out ) of output light. Therefore, the DoP = Q 2 out + U 2 out + V 2 out /I 2 out , where I out , Q out , U out , and V out are the four Stokes elements of the S out . The polarization difference (PD) imaging can be concluded with Ico and Icross.
Step 2: Imaging based on PR method. The second step aims to retrieve the original image with higher quality. Indeed, the traditional imaging method loses its effect when the cloudy medium is very heavy. The PR imaging method can complete retrieve the original image in theory. Then, we compare the imaging quality under different methods through the corresponding image-quality indicators: peak signal to noise ratio (PSNR) and structural similarity index (SSI).
Step 3: Discussion. Finally, in the third step, to better analyze the PR imaging model in atmospheric scattering medium, we consider the following points: (1) effects of different scattering environments on imaging quality; (2) the penetration characteristics of light at different wavelengths in the atmosphere; and (3) MM's applicability with changing scattering particle sizes.

Results and Discussions
We simulated the whole imaging process, from which we could obtain the corresponding images by the relative imaging methods. From Table 2, we can get the scattering coefficient and absorption coefficient under the clear atmosphere and calculate the extinction coefficient in each layer, and then multiply the extinction coefficient by the corresponding transmission distance to get the optical thickness of 1τ = 1.8. The light intensity transmittance (T1) in this environment can be calculated from Equation (7) as 16.63%. Figure 5 shows the imaging results from intensity (I components of Stokes vector), DoP, PD, and PR methods through scattering medium with different concentrations. From left to right, the optical thicknesses, which is the product of the extinction coefficient and the transmission Remote Sens. 2020, 12, 2895 9 of 16 distance, were set as 1τ, 3τ, 6τ, and 9τ, respectively, by changing the extinction coefficients (changing the number density and radius of particles), which demonstrates the turbidity of the medium systems. In the four different scenarios in Figure 5, the clear weather is set as a reference standard. In clear weather, the scattering coefficient and absorption coefficient of 1τ are increased by three, six, and nine times to obtain other three different scenes. 3τ, 6τ, and 9τ equal 3.7, 6.9, and 10.1 respectively. Thus, the other three different light transmittances under cloudy weather are T2 = 2.532%, T3 = 0.1012%, and T4 = 0.0039%, respectively. In each row, Figure 5a-d presents imaging results with four different imaging methods in different atmosphere conditions (with the optical thicknesses of 1τ, 3τ, 6τ, and 9τ, respectively). The first column shows different imaging results in a clear day (with the optical thickness of τ). It could be seen that the intensity imaging identifies the target with a bad performance even on a clear atmosphere. However, it should be noted that the intensity imaging is covered with a thin "white mist" due to effect of backscattered light without target information. In contrast, the DoP, PD, and PR methods can reconstruct the image without obvious distortions except for that of the larger optical thickness of 9τ.
Remote Sens. 2020, 12, x 9 of 16 When the optical thickness of the medium increases gradually, the intensity image is covered by more backscattered light, from which the object information cannot be obtained absolutely. Although the results of the other three methods were slightly reduced, they can still get the outline of the intermediate target. In the fourth column of Figure 5 (with the optical thickness of 9τ), I, DoP, and PD methods have lost their effects, but PR method can still get a relatively clear pattern, which is consistent with the theoretical analysis. Moreover, absorption of medium has been considered in the medium, and, due to the longer scattering distance, there will be too many scattering times, in which most of the polarization photons from the target would die out. Thus, compared with DOP method, PD method has no significant improvement in image visibility, because only a few depolarization photons have been eliminated based on the PD method.
To analyze the advantages and disadvantages of these four methods numerically, we calculated the contrasts of the 16 images in Figure 5, and the formula for calculating contrast can be expressed as: C = |IT − IB|/|IT + IB|, where IT and IB express the intensities of the targets and the background, respectively. In addition, the reconstructed effects of PR methods for targets with different geometric complexities are also discussed. Thus, we define the model of the main building in the school badge of Hefei University of Technology as Target 1, the English letters "HFUT" and the figure "-1945-" as Target 2, and the circular marble section as background. IT is calculated from the mean intensity values of Targets 1 and 2, respectively and IB is calculated from the average strength value of the background. Obviously, the geometry of Target 1 is more complex.
In Figure 6a,b, the x-coordinate represents different scattering environments with increasing optical thickness of 1τ, 3τ, 6τ, and 9τ, respectively, and the y-coordinate represents the contrast of the images. We can find that the image contrast retrieved by PR method has a small improvement when the scattering is small, but, when the medium becomes very cloudy, such as 6τ and 9τ, the PR method improves significantly. Comparing Figure 6a,b, we can find the contrast of Target 2 is higher than When the optical thickness of the medium increases gradually, the intensity image is covered by more backscattered light, from which the object information cannot be obtained absolutely. Although the results of the other three methods were slightly reduced, they can still get the outline of the intermediate target. In the fourth column of Figure 5 (with the optical thickness of 9τ), I, DoP, and PD methods have lost their effects, but PR method can still get a relatively clear pattern, which is consistent with the theoretical analysis. Moreover, absorption of medium has been considered in the medium, and, due to the longer scattering distance, there will be too many scattering times, in which most of the polarization photons from the target would die out. Thus, compared with DOP method, PD method has no significant improvement in image visibility, because only a few depolarization photons have been eliminated based on the PD method.
To analyze the advantages and disadvantages of these four methods numerically, we calculated the contrasts of the 16 images in Figure 5, and the formula for calculating contrast can be expressed as: C = |I T − I B |/|I T + I B |, where I T and I B express the intensities of the targets and the background, respectively. In addition, the reconstructed effects of PR methods for targets with different geometric complexities are also discussed. Thus, we define the model of the main building in the school badge of Hefei University of Technology as Target 1, the English letters "HFUT" and the figure "-1945-" as Target 2, and the circular marble section as background. I T is calculated from the mean intensity values of Targets 1 and 2, respectively and I B is calculated from the average strength value of the background. Obviously, the geometry of Target 1 is more complex.
In Figure 6a,b, the x-coordinate represents different scattering environments with increasing optical thickness of 1τ, 3τ, 6τ, and 9τ, respectively, and the y-coordinate represents the contrast of the images. We can find that the image contrast retrieved by PR method has a small improvement when the scattering is small, but, when the medium becomes very cloudy, such as 6τ and 9τ, the PR method improves significantly. Comparing Figure 6a,b, we can find the contrast of Target 2 is higher than that of Target 1 for the four different methods, which shows that the complexity of the target geometry structure also affects the contrast of the image.
as: C = |IT − IB|/|IT + IB|, where IT and IB express the intensities of the targets and the background, respectively. In addition, the reconstructed effects of PR methods for targets with different geometric complexities are also discussed. Thus, we define the model of the main building in the school badge of Hefei University of Technology as Target 1, the English letters "HFUT" and the figure "-1945-" as Target 2, and the circular marble section as background. IT is calculated from the mean intensity values of Targets 1 and 2, respectively and IB is calculated from the average strength value of the background. Obviously, the geometry of Target 1 is more complex.
In Figure 6a,b, the x-coordinate represents different scattering environments with increasing optical thickness of 1τ, 3τ, 6τ, and 9τ, respectively, and the y-coordinate represents the contrast of the images. We can find that the image contrast retrieved by PR method has a small improvement when the scattering is small, but, when the medium becomes very cloudy, such as 6τ and 9τ, the PR method improves significantly. Comparing Figure 6a,b, we can find the contrast of Target 2 is higher than that of Target 1 for the four different methods, which shows that the complexity of the target geometry structure also affects the contrast of the image. According to the atmospheric distribution characteristics, there are more scattered particles and aerosols near the ground; once away from the ground, the atmosphere would become thinner. Thus, it will result in a negative correlation between atmospheric extinction coefficient and altitude. On the other hand, especially within 0-2 km of the Earth's surface, the particle number density (PND) and radius of aerosol are always changes with weather conditions. In general, when the temperature or the wind speed is higher, the Brownian motion of particles will be intensified, and the diffusion ability caused by Brownian motion of particles will also increase. Of course, under different atmospheric environment conditions, the MM of atmospheric environment is also different. However, the atmosphere is a dynamic environment; thus, when using PR method for atmospheric remote sensing imaging, we cannot calculate the changing environment's MM in time. Thus, it is important for the real applications to determine the MM's recovery range when environmental parameters changes in 0-2 km atmosphere. In other words, we investigated whether the original MM can still retrieve the original image well when environmental parameters changes based on the PR method. When the particle size increases, the extinction coefficient of the medium will increase correspondingly. The changing particle size was used to simulate the change of atmospheric environment, and the initial conditions were set corresponding to clear atmosphere. The particle radius increases by 0.1-0.5 times correspond to particle size of clear atmosphere in 0-2 km atmospheres above the ground. The simulation results are demonstrated in Figure 7, in which simulation parameters from 1d to 1.5d mean the particles' sizes are increased from 1.0 to 1.5 times. Figure 7a-f shows that the contrast and visibility of image gradually decrease, in which Figure 7a is for reference, corresponding to the retrieved image based on PR method under sunny day. Obviously, the air is clear and the MM is proper for retrieving, so the retrieved image is the best. As the air becomes more turbid, we still use the original MM for restoration; in Figure 7b-d, imaging quality gradually decreases, and we can still distinguish the outline of the clear target. However, in Figure 7d, the contrast has dropped dramatically, and, in Figure 7e, the imaging quality has become very poor. As depicted in Figure 7f, the target cannot be identified completely. The results show the original MM has a suitable range for image restoration in different atmospheric environments.
Obviously, the air is clear and the MM is proper for retrieving, so the retrieved image is the best. As the air becomes more turbid, we still use the original MM for restoration; in Figure 7b-d, imaging quality gradually decreases, and we can still distinguish the outline of the clear target. However, in Figure 7d, the contrast has dropped dramatically, and, in Figure 7e, the imaging quality has become very poor. As depicted in Figure 7f, the target cannot be identified completely. The results show the original MM has a suitable range for image restoration in different atmospheric environments. Similarly, we calculated the contrast of images, as shown in Figure 8. The x-axis shows increased times of the particle size. Black and red solid lines are the contrast of Targets 1 and 2 in the original image restored by MM, which is the reference value. The results show that the image contrast of Target 1 was always higher than that of Target 2. However, with increasing the radius of particle, the contrasts of both Targets 1 and 2 continue to decline. We can find that, when the increased size factor is less than or equal to 0.3, the contrasts of Targets 1 and 2 decrease very slowly. When the increased Similarly, we calculated the contrast of images, as shown in Figure 8. The x-axis shows increased times of the particle size. Black and red solid lines are the contrast of Targets 1 and 2 in the original image restored by MM, which is the reference value. The results show that the image contrast of Target 1 was always higher than that of Target 2. However, with increasing the radius of particle, the contrasts of both Targets 1 and 2 continue to decline. We can find that, when the increased size factor is less than or equal to 0.3, the contrasts of Targets 1 and 2 decrease very slowly. When the increased times are bigger than 0.3, both contrasts drop sharply. When it is bigger than 0.45, the contrasts of the image become very small, and, as shown in Figure 7f, the target completely disappears from the speckle map.
Remote Sens. 2020, 12, x 11 of 16 times are bigger than 0.3, both contrasts drop sharply. When it is bigger than 0.45, the contrasts of the image become very small, and, as shown in Figure 7f, the target completely disappears from the speckle map. The imaging effect in the scattering medium is often closely related to the incident wavelength of light. In water, for example, to avoid absorbing too much, use blue-green light as incident light. In the atmospheric environment, compared with the visible light, the infrared wavelength usually has better penetration characteristics. Figure 9 shows the spectral characteristics of PR-based imaging under different weather conditions (clear in first row and cloudy in second row); we obtained the images in visible band (633 nm), near infrared (1536 nm), and mid-infrared (4 μm), respectively, and calculated the corresponding image contrasts. Comparing the two groups of images (clear and cloudy), it can be found that different imaging bands have different performances, in which visible light with an incident wavelength of 633 nm is no longer able to recognize the target's information at all, but the incident light in infrared band can still present the original target information well, and the contrast is increased about 3-4 times. Figure 9b,c shows that the contrasts between bands of 1536 nm and 4 μm are almost the same. As depicted in Figure 9e,f, the imaging performances are slightly worse at the same band than in clear weather (Figure 9b,c), which can be attributed to the cloudier atmosphere. Moreover, we can also find that the imaging performances is also influenced by the wavelength in the infrared band, and the contrast of Figure 9f is approximately 1.47 times that of The imaging effect in the scattering medium is often closely related to the incident wavelength of light. In water, for example, to avoid absorbing too much, use blue-green light as incident light. In the atmospheric environment, compared with the visible light, the infrared wavelength usually has better penetration characteristics. Figure 9 shows the spectral characteristics of PR-based imaging under different weather conditions (clear in first row and cloudy in second row); we obtained the images in visible band (633 nm), near infrared (1536 nm), and mid-infrared (4 µm), respectively, and calculated the corresponding image contrasts. Comparing the two groups of images (clear and cloudy), it can be found that different imaging bands have different performances, in which visible light with an incident wavelength of 633 nm is no longer able to recognize the target's information at all, but the incident light in infrared band can still present the original target information well, and the contrast is increased about 3-4 times. Figure 9b,c shows that the contrasts between bands of 1536 nm and 4 µm are almost the same. As depicted in Figure 9e,f, the imaging performances are slightly worse at the same band than in clear weather (Figure 9b,c), which can be attributed to the cloudier atmosphere. Moreover, we can also find that the imaging performances is also influenced by the wavelength in the infrared band, and the contrast of Figure 9f is approximately 1.47 times that of Figure 9e. We also note that the difference will become more pronounced when the medium becomes more turbid.
under different weather conditions (clear in first row and cloudy in second row); we obtained the images in visible band (633 nm), near infrared (1536 nm), and mid-infrared (4 μm), respectively, and calculated the corresponding image contrasts. Comparing the two groups of images (clear and cloudy), it can be found that different imaging bands have different performances, in which visible light with an incident wavelength of 633 nm is no longer able to recognize the target's information at all, but the incident light in infrared band can still present the original target information well, and the contrast is increased about 3-4 times. Figure 9b,c shows that the contrasts between bands of 1536 nm and 4 μm are almost the same. As depicted in Figure 9e,f, the imaging performances are slightly worse at the same band than in clear weather (Figure 9b,c), which can be attributed to the cloudier atmosphere. Moreover, we can also find that the imaging performances is also influenced by the wavelength in the infrared band, and the contrast of Figure 9f is approximately 1.47 times that of Figure 9e. We also note that the difference will become more pronounced when the medium becomes more turbid.  In addition, to further verify the superiority and effectiveness of our method, we used various evaluation criteria to evaluate the imaging quantify, including PSNR and SSI [53]. A higher value of PSNR or SSI indicates higher image quality. SSI values range 0-1, and the larger is the value, the smaller is the image distortion. The unit of PSNR is dB, and the larger is the PSNR value, the lower is distortion, which can be expressed as: where n is the number of bits per pixel, and the general grayscale image is taken as 8, that is, the grayscale of pixel is 256. For two m × n monochrome images X and Y, if one is the noise approximation of the other, then MSE is the square of the error between the estimated value and the reference value. Here, we define the 1τ image as the reference image because the 1τ image is very clear. Similarly, SSI can be expressed as: where C 1 and C 2 are constants; µ X and µ Y represent the mean values of the image X and Y, respectively; σ X and σ Y represent the variances of image X and Y, respectively; and σ XY represents the covariance of image X and Y. SSI measures image similarity from brightness, contrast, and structure. Tables 3-5 show values of PSNR and SSI for different images obtained in different conditions, respectively. Table 3. Quantitative comparison of recovered images for Figure 5.   Table 5. Quantitative comparison of recovered images for Figure 9.   Table 3 demonstrates the quantitative comparison of different recovered images that is obtained in different environmental conditions with different optical thicknesses. When the optical thickness of the medium is 1τ, we can see that the PSNR value of the recovered image based on the PR method is about six times that of the intensity image, and, in this low concentration medium, the SSI results show that the recovered image is basically the same as the original. When the optical thickness reaches 9τ, it still maintains a high value, even though other methods are unable to image. Figure 7 explores the MM's retrieving range. In Table 4, PSNR and SSI numerical display that, when the particle size is greater than 1.3 times, the quality of the image drops so dramatically that the MM is no longer useful. Such results are of great significance for us to explore the applicable range of MM in different environments.
PSNR and SSI of Figure 9, as displayed in Table 5, also demonstrate the universality and reliability of our method and theory. At the wavelength of 633 nm, both the PSNR and SSI of images show relative smaller values in both clear and cloudy atmospheres. When the incident light lies in the infrared bands, and the values of PSNR increases nearly eight times compared to that under clear atmosphere. Similarly, the SSI values of the images is nearly double. When the atmospheric condition is cloudy, the trend that the imaging quality is proportional to the wavelength will be more obvious, as depicted in Table 5.

Conclusions
In this paper, we apply polarization imaging theory to real turbid environments (from the satellite to the ground) based on stratification model and PR theory. The simulation results show that PR imaging method can be used for object discrimination in turbid atmosphere medium with enhanced imaging contrasts. Compared with intensity imaging, DoP, and PD imaging methods, PR imaging shows superior power to distinguish the targets and its background. Through calculating the contrast of images, the proposed PR imaging method can further reduce the influence of the turbid atmosphere medium, and, in a highly turbid atmosphere environment, the imaging based on PR method still has higher contrast and visibility than those of intensity, DoP, and PD imaging methods. Meanwhile, when atmosphere scattering medium becomes more turbid in 0-2 km, the target can also be retrieved to a certain extent based on PR method with the original MM. Then, we demonstrated that the PR imaging method has higher contrast and visibility for a target with lower geometric complexity. Finally, the influence of the incident wavelengths to the retrieving performance was also investigated, and the results show a longer wavelength is better for transmitting in the real environment. In addition, we also introduced image evaluation indices, such as PSNR and SSI, to evaluate the images quantitatively, which shows the reliability of our proposed method. These results are significant for active imaging in remote sensing atmospheric environments.