1. Introduction

With the developing demands on the tasks like target detection, polarization information become more and more important than ever, which can provide more additional information, such as degree of polarization, angle of polarization, the image in a certain Stokes vector and so on. The interaction of light with the interface of two media will affect the concrete polarization information greatly. Therefore, it is very important to characterize the process of the reflective polarization for exploring the characteristics of the polarization distributions and revealing the features of the surface, which can be used in polarization remote sensing [1,2], target detection and recognition [3–6].

To describe the reflective process on rough surface, a bidirectional reflection distribution function (BRDF) model has been proposed based on the micro-facets theory [7]. Combined the micro-facets theory with polarization information, the Mueller matrix is introduced in the original BRDF model as a replacement of the Fresnel factor to deal with the polarization information. And many expansions and explorations have also been made in this fields later [8–10]. The simulated results of the modified polarized BRDF model has also been compared with the measured data [11], which shows that the model is suitable for the circumstances where single reflection on the surface dominates. Then, some works have been reported with the BRDF in arbitrary layers' media [12,13]. And the reflective polarization on rough surface with smooth transparent layers have also been studied [14,15]. The Monte Carlo method obtains numerical results based on repeated sampling and statistics, which makes it very flexible and effective to simulate a physics phenomenon and get the results. Several Monte Carlo methods have been provided to investigate the propagation of light in multi-layer media [16–18]. The Monte Carlo methods have also been introduced to solve the cases with polarization [19–28]. Some of them focus on calculating the polarization in an atmosphere-ocean system or turbid media [24–27]. Some of them are applied for solving the radiative transfer in the atmosphere [21–23,28].

In this paper, based on the Monte Carlo model, we make an investigation on the reflective polarization from a multi-layer surface. In Section 2, we review the specular reflection and refraction described by Fresnel’s law. In Section 3, the reflective polarization model based on the Monte Carlo algorithm is presented in detail. The structure of a multi-layer surface is introduced at first. Then, based on the geometric optics and micro-facets theory, the normal direction of each micro-facet is calculated by sampling from its normal distribution. At last, the Monte Carlo calculating processes are introduced and described carefully. In Section 4, we make a comparison of our simulated results with the reported measured data, and they agree well with each other, which demonstrate the correctness and effectiveness of our proposed simulated model. And we have also analyzed the effects of different surface structures and different incidences (different wavelength, different polarization states and so on) to the reflective polarization information. In Section 5, we make a summary and prospect of this paper.

2. Reflection and refraction of polarized light at interface

Generally, light will be reflected and refracted when it reaches to an interface of two media. Fresnel and Snell’s law can describe the reflection and refraction on smooth interface between layer 1 and layer 2 whose refraction indexes are $n_0$ and $n_1+ik$ respectively, as shown in Fig. 1.

 

Fig. 1 Reflection and refraction at interface.

Download Full Size | PDF

Each linear polarization could be decomposed into s and p polarization components, denoted by $E_s$and$E_p$, so that $E_{rs}$($E_{rp}$) and $E_{ts}$($E_{tp}$) are the reflective and refractive s (p) polarization components respectively. According to Fresnel’s law, the reflection factors$F_s$, $F_p$ and respective phase shifts ${\delta }_s$, ${\delta }_p$ can be expressed as [29,30]:

The Mueller matrix for Fresnel reflection can be written as [29,30]:

Then the Stokes vector of the reflection can be easily calculated by incident Stokes vector and the Mueller matrix of reflection.

As for the refraction, considering the criterion of energy conservation, the transmission factors are $T_s=1-F_s$, $T_p=1-F_p$. Therefore, the Mueller matrix for refraction can be expressed as [29,30]:

3. Monte Carlo model for rough-surface reflection

3.1 Multi-layer surface structure

Here, we just consider a two-dimensional multi-layer surface. There are a few assumptions for the surface in our model:

From the top to the bottom, a multi-layer surface is composed of several transparent layers and an opaque or metallic substrate. Figure 2 describes the cross section of a three-layer surface with two transparent layers (layer 1 and 2) and an opaque or metallic substrate (layer 3). The incident light will be reflected and refracted at each interfaces between layers and traced until being collected after leaving the surface or being ignored due to lacking enough energy.

 

Fig. 2 Model of multi-layer surface.

Download Full Size | PDF

3.2 Monte Carlo model for reflective polarization from multi-layer surface

To acquire the reflective polarization information from an arbitrary multi-layer surface, the Monte Carlo algorithm is employed in this paper to trace and simulate the light propagation processes effectively.

Firstly, based on micro-facets theory, each layer of a surface is assumed to be made up of many micro-facets, and the normal vector of each micro-facet can be calculated by sampling from its normal distribution. Secondly, the procedures of lights’ propagation in multi-layer media are traced repeatedly. At last, the polarization information of reflections is collected in the upper hemisphere area.

3.2.1 Normal vector calculation

Based on the micro-facet theory, each layer of surface is assumed to be made up of micro-facets which are perfectly smooth, so that those micro-facets reflect and refract light specularly. As the normal vector is crucial in the computation of specular reflection and refraction with the Fresnel’s law, a method to sample the normal vector of the micro-facets will be introduced firstly.

A rough surface located in the surface coordinate system (XYZ coordinate system) can be decomposed into lots of micro-facets as shown in Fig. 3(a). Based on the micro-facets theory, there are a lot of methods to describe a rough surface statistically. A physically-based method to sample normal direction from an anisotropic surface is employed here [31,32]. This method presents a normal direction distribution of rough surfaces. To model an anisotropic surface, two empirical but intuitive parameters $e_x$ and $e_y$ are introduced to describe the roughness along X direction and Y direction, and to control the azimuth angle dependency of normal direction.

 

Fig. 3 (a) Surface composed of the micro-facets in the XYZ coordinate system, (b) the schematics of a single micro-facet.

Download Full Size | PDF

The probability density function for the normal direction distribution can be expressed as:

In this distribution, if$e_x=e_y$, the surface roughness has no dependency on azimuth angle and the surface is isotropic, while if$e_x\ne e_y$, the surface becomes anisotropic. And the larger the parameters are, the smoother the surface will be. The normal direction $\overrightarrow{n}$ of micro-facet shown in Fig. 3(b) can be sampled from the probability distribution function in Eq. (9) by applying the following equation:

3.2.2 The transformation of polarization coordinate

Firstly, we introduce a coordinate system to trace and measure the polarized light. Light propagating in direction $\overrightarrow{d}$ can be decomposed into any two orthogonal polarization components which oscillates in directions of $\overrightarrow{s}$ and $\overrightarrow{p}$, respectively. The combination of these direction vectors $(\overrightarrow{d},\overrightarrow{s},\overrightarrow{p})$ forms the polarization coordinate system of the light.

In our simulation, we employ two special polarization coordinate system that are global polarization coordinate system and local polarization coordinate system.

Assuming there is a light propagating in XYZ coordinate system with the direction of $\overrightarrow{d}=\left(\sin \theta \cos \delta ,\sin \theta \sin \delta ,\cos \theta \right)$, in which $\theta$ and $\delta$are the zenith angle and azimuth angle respectively.

The incident light can be decomposed in any orthogonal direction and expressed with different Stokes Vector. While, in order to measure these polarized light uniformly, a coordinate system is needed. To achieve this goal, we establish a connection of the polarization coordinate system with the surface coordinate system (XYZ coordinate system). That is, the surface$\overrightarrow{Z}$direction is employed in the polarization decomposition, which leads to $\overrightarrow{s_g}=\overrightarrow{d}\times \overrightarrow{Z}$ and $\overrightarrow{p_g}=\overrightarrow{s_g}\times \overrightarrow{d}$. The combination of $(\overrightarrow{d},\overrightarrow{s_g},\overrightarrow{p_g})$ forms the global polarization coordinate system for the incident light. With this global polarization coordinate system, any polarized light in our simulation can be observed and measured clearly.

While, the global polarization coordinate system is not enough to account for the interaction of light with micro-facets. For applying the Fresnel's law, a polarized light should be decomposed into the local coordinate system determined by the normal vector $\overrightarrow{n}$ of the selected facet and the light propagating direction. Therefore, if the light reaches to a micro-facet with the normal vector of $\overrightarrow{n}$ as shown in Fig. 4, the polarized light can be decomposed as $\overrightarrow{s_l}=\overrightarrow{d}\times \overrightarrow{n}$ and $\overrightarrow{p_l}=\overrightarrow{s_l}\times \overrightarrow{d}$respectively. The obtained$(\overrightarrow{d},\overrightarrow{s_l},\overrightarrow{p_l})$is the local polarization coordinate system for the incident light.

 

Fig. 4 Polarization coordinate system rotation.

Download Full Size | PDF

In our simulating process, because all of the incident and scattering lights are detected in global polarization coordinate system, while all of the calculations for the interaction between the incident light and the selected micro-facets are arising in the local polarization coordinate system, we should know the concrete rotating angles of${\theta }_{rot}$between the global polarization coordinate system and the local polarization coordinate system for every selected micro-facets. For example, if the polarization coordinate system of a light is needed to be transformed from the global one to the local one, it will be rotated from $\overrightarrow{s_g}$ to $\overrightarrow{s_l}$ as shown in Fig. 4, vice versa. And the rotation angle can be expressed as ${\theta }_{rot}=\text{arc}\cos \left(\overrightarrow{s_g}·\overrightarrow{s_l}/\left|\overrightarrow{s_g}\right|\left|\overrightarrow{s_l}\right|\right)$ that lies in the range of 0 to π. In fact, due to the definition range of the inverse cosine function, ${\theta }_{rot}$do not contain enough rotating direction information, i.e. rotating clockwise or anticlockwise. So it is necessary to determine the sign of ${\theta }_{rot}$, i.e. rotating ${\theta }_{rot}$or $-{\theta }_{rot}$, to make ${\theta }_{rot}$ cover the range of -π to π. Here, a new vector $\overrightarrow{m_{lg}}= \overrightarrow{s_l}\times \overrightarrow{s_g}$ is introduced to determine the sign. From the definition, $\overrightarrow{m_{lg}}$is always parallel to $\overrightarrow{d}$, and if $\overrightarrow{m_{lg}}$points to the same direction as $\overrightarrow{d}$ does, the rotation angle is ${\theta }_{rot}$, otherwise the rotation angle is $-{\theta }_{rot}$. The rotation operation in 3D coordinate system is $R{(\theta )}_{rot}$ which can be expressed in the formation of Mueller matrix $M{(\theta )}_{rot}$ [29, 30]:

3.2.3 Monte Carlo model

To acquire the reflective polarization information from the multi-layer surface, parameters, such as incident light direction, incident polarization states of $S_{in}=\left[IQUV\right]$, and the layer number of the surface, the refraction index, the normal distribution and roughness parameters $e_x$ and $e_y$ of every layers, should be set firstly. Due to the property of the Monte Carlo method, reflective polarization information is calculated statistically by repeated sampling, so the sample number of incident photons should also be set at the same time.

When light transmits from the $i_{th}$ layer to $j_{th}$ layer, we use $S_{r\_ij}$ and $S_{t\_ij}$ to label the Stokes vectors of reflected and refracted lights respectively. Specifically, if the light propagates upwards, $j=i-1$; if the light propagates downwards, then$j=i+1$. By default, lights are decomposed and measured in global polarization coordinate system. And when a light is decomposed in local polarization coordinate system, its Stokes Vector is labeled with subscript$l$. Lights in global polarization coordinate system are not labeled in our following express.

In the beginning, the air is set as layer 0. Before the incident light hit the first layer of a surface, the light is in layer 0 and will interact with layer 1, and its Stokes Vector $S_{in}$can be denoted as$S_{01}$.

We cast incident lights onto a multi-layer surface, and the lights are traced as follows:

During the tracing procedure, if the layer is not transparent, we will only calculate and trace the reflected light on it. And if light energy is lower than a threshold, the light will be abandoned.

For example, considering the case that there is a rough surface with 3 layers as depicted in Fig. 5, incident light labeled as $S_{in}$ propagates downwards in the air at first. When the incident light reaches to the interface between layer 0 and layer 1, light is reflected and refracted into $S_{r\_01}$ and $S_{t\_01}$. $S_{t\_01}$ will propagate downwards in layer 1, then be reflected and refracted at the next interface, and generate reflected and refracted light $S_{r\_12}$ and $S_{t\_12}$. The $S_{t\_12}$ will go upwards and be refracted and reflected at the interface between layer 1 and layer 0 into $S_{t\_01}$ and $S_{r\_10}$. The $S_{t\_12}$ will only be reflected at the interface between layer 2 and layer 3 due to that the layer 3 is opaque or metal material. The $S_{r\_23}$ will go upwards and produce $S_{r\_21}$, $S_{r\_10}$ and $S_{t\_10}$ finally. Then, two $S_{r\_10}$ and one $S_{r\_21}$ will be set as new incidences and traced as$S_{in}$. At last, $S_{r\_01}$ and all $S_{t\_10}$ that leave the top of the surface will be collected.

 

Fig. 5 Multi-layer reflections and refractions.

Download Full Size | PDF

Reflective polarization information will be collected in the upper hemisphere using a Stokes vectors matrix consisting of 91 × 360 grids with a step of 1°. Here, the surfaces are assumed to be a tiny point and regarded as the spherical center of the hemisphere. The reflected lights leaving from the top layer of the surface are collected according to their azimuth $\delta$and zenith angle$\theta$ in the corresponding grid of the matrix, which means that the lights with the same propagation direction will be counted in the same grid and their Stokes vector will be summed up. Finally, the obtained Stokes vectors matrix with 91 × 360 elements can describe the reflective polarization information from the multi-layer surfaces effectively.

4. Experiment and analyses

To prove the correctness of our Monte Carlo model, firstly, we will try to compare the simulated results by our model with the reported measured data and the results by analytical BRDF model. Then, we will also investigate the influences of the layer number, the incident wavelengths, the surface geometry (roughness) and the incident polarization states on the obtained reflective polarization information.

4.1 Comparison with the reported measured data

Firstly, the simulated results by our model have been compared with the reported measured data [33]. We simulate the reflective polarization information of an aluminum and a flat green paint surface respectively. The parameters of our model are set as follows: the polarization state of the incident light is $S_{in}=[1000]$; the light sample number is one million. The refraction index of flat green paint is $1.47+0.47i$. And in the simulation, if the intensity of the light is less than 5% of the original ($\tau =0.05$), we will stop tracing the light and ignore its contribution to the final reflective information. According to the description of the measured sample, the aluminum is sandblasted and the flat green paint is covered on foam. As an estimation of the roughness of the measured sample, the surface roughness parameters $e_x=e_y=15$ are employed in our simulation.

Figure 6 shows the variations of the obtained reflective degree of polarization (collected at the reflection angle of 60°) of the aluminum surface and flat green paint respectively, with the incident angles ranging from 40° to 80°. For the flat green paint, the obtained reflective degree of polarization reaches to its peak when the incident angle is around 55°, which agrees well with the reported measured data in experiments [33]. For the aluminum surface, as the incident angles increase from 40° to 80°, the obtained reflective degree of polarization keeps rising, which also agrees well with the experimental data. It indicates that our model is correct and effective to achieve the reflective polarization information.

 

Fig. 6 Comparisons with measured data of the aluminum surfaces (a) and flat green paint (b).

Download Full Size | PDF

4.2 Comparisons with the analytical model of BRDF

Our Monte Carlo model has also been compared with the analytical micro-facet based BRDF model included in the SCATMECH [34] that is a light scattering library and published by the NIST in 2015.

Firstly, we investigate single layered copper surface (with the refraction index of $0.23+3.46i$), and the parameters of our model are set as follows: the incident polarization state is non-polarized $S_{in}=[1000]$; the wavelength is 635nm; the incident angle is 40°; the light sample number is one million and the threshold is$\tau =0.05$ and the surface roughness parameters are $e_x=e_y=15$.

As shown in Fig. 7, the reflective angle of polarization (AoP), degree of polarization (DoP) and Stokes I, Q, U and V parameters from our model and the analytical BRDF model are listed from the left to the right. It is obvious that the reflective polarization pattern from our model matches the results from analytical BRDF model very well. In the simulation, the normalized total reflective intensity over the whole upper hemisphere space from our model is 0.8692, which satisfies the criterion of energy conservation. It should be noted that the total reflective intensity from the analytical BRDF model is bigger than 1.

 

Fig. 7 Comparison with single-layer analytical BRDF model.

Download Full Size | PDF

Our model is then compared to a multi-layer analytical model included in the SCATMECH. We use the same parameters of the incident light, and employ a sample of two-layered surface. And the refraction index for each layer is 1.5 and $0.23+3.46i$ (copper) from the top to the bottom. As depicted in Fig. 8, the patterns of reflective polarization obtained by our model are still similar with that of analytical BRDF model. And our model still satisfies the criterion of energy conservation, while the analytical model still have some problems in some reflection angles which also results in the inaccuracy of the distribution of degree of polarization.

 

Fig. 8 Comparison of the results between our proposed mode and the analytical BRDF model for the multi-layer sample.

Download Full Size | PDF

4.3 Analyses and discussions

In this section, with our proposed Monte Carlo model, we will investigate the influences of the layer number, the incident wavelengths, the surface geometry and the incident polarization states on the obtained reflective polarization information.

4.3.1 The influence of the layer number

One advantage of our Monte Carlo model is that it could deal with a surface with arbitrary layers. Firstly, we will discuss the influence of the layer numbers on the reflective polarization information.

In this experiment, we use three kinds of surfaces and all the materials of the last layer of these surfaces are the copper (refraction index:$1.31+1.85i$). Surface one is a three-layered surface, with the refraction indexes of 1.3, 1.4 and $1.31+1.85i$ from the top to the bottom. The surface two has five layers with refraction indexes of 1.3, 1.4, 1.5, 1.6 and $1.31+1.85i$ from the top to the bottom. And the surface three has seven layers with refraction indexes of 1.3, 1.4, 1.5, 1.6, 1.7, 1.8 and $1.31+1.85i$ from the top to the bottom.

All dielectric layers in the three samples have the same roughness parameters of$e_x=e_y=30$, and the substrates of the three samples are smoother with roughness parameters of $e_x=e_y=40$. And the incident polarization state is$S_{in}=[1000]$; the incident wavelength is 340nm; the incident light sample number is three million; the incident angle is 40° and the threshold is set as$\tau =0.05$.

The corresponding results (the reflective polarization information) from the samples with three, five and seven layers are listed in Fig. 9 respectively. It is clear that with increasing the number of layers, all Stokes parameters (I, Q, U and V), become weaker, while all of the angle of polarization appear to be similar. It can be concluded that with the increment of the layer, the degree of polarization of the obtained light and the total intensity continue decreasing, while the angle of polarization properties are kept well, which could be useful in the target detection.

 

Fig. 9 The reflective polarization information from 3 layers, 6 layers and 8 layers samples.

Download Full Size | PDF

4.3.2 Influence of surface geometry

To study the influence of surface roughness on the reflection, we make a comparison using our Monte Carlo model with different surface roughness parameters. In order to show the variation clearly, we employ a single-layer copper surface (without coatings) as the sample, and show the concrete simulating results about the Stokes parameter I and degree of polarization only. The incident angle, the incident wavelength and the threshold are set as 40°, 340nm and $\tau =0.05$ respectively.

For an anisotropic surface labeled with surface A, the roughness are set as $e_x=30$, $e_y=15$respectively. As contrasts, surface B is isotropic surface with $e_x=e_y=30$, and surface C is anisotropic with $e_x=45$, $e_y=30$accordingly. As shown in Fig. 10, from the left to the right, we have obtained the Stokes parameter of I and degree of polarization from surface A, B and C respectively. It is clear that, compared to the isotropic surface B, the reflection pattern of surface A are stretched in the horizental direction, because as $e_y$ decreases from 30 to 15, the surface becomes more diffuse horizontally. While as $e_x$ arise from 30 to 45, the vertical direction of surface C becomes smoother than that of surface B, which results in a more compressed reflection area in surface C.

 

Fig. 10 Reflection distributions of the samples with difference surface parameters.

Download Full Size | PDF

4.3.3 Effect of the incident wavelength

The effect of the incident wavelength is discussed in detail with employments of two kinds of surfaces which are made of copper and aluminum respectively. Parameters are set as follows: the polarization state of the incident light is$S_{in}=\left[1000\right]$, the incident angle is 40°, the incident light sample number is three million and the threshold is$\tau =0.05$. The top layer of the two surface are both covered by a dielectric layer (refraction index:1.5), and the substrate is copper or aluminum. All of the layers have the same roughness parameters of $e_x=e_y=10$.

The refraction indexes of copper are $1.31+1.85i$, $0.95+2.57i$ and $0.21+3.86i$ at the wavelength of 340nm, 550nm and 670nm respectively. As shown in Fig. 11, we have obtained the reflective polarization patterns (angle of polarization, degree of polarization, Stokes I, Q, U, V) at the incident wavelengths of 340nm, 550nm and 670nm from the top row to the bottom. It is obvious that as the wavelength increases, the angle of polarization and Stokes parameter of I become stronger, while the values of the degree of polarization and Stokes parameters of Q, U and V become weaker.

 

Fig. 11 Polarization reflection distribution of copper under three wavelengths.

Download Full Size | PDF

The refraction indexes of aluminum at the wavelength of 340nm, 550nm and 670nm are $0.35+4.11i$,$0.96+6.70i$ and $1.60+8.01i$ separately. The obtained results are listed in Fig. 12 as Fig. 11 does. Unlike the copper surface, the obtained results of the aluminum under different incident wavelengths do not show many differences. The simulated reflective intensities in these two experiments are both consistent with the copper and aluminum spectral reflectance curve. And according to color theory, this is the reason that the copper and aluminum look reddish brown and shining grey [35]. From this simulation results, it can be concluded that the reflective polarization from copper surface is more sensitive to the variation of wavelength than the reflection from aluminum surface. And surfaces with different materials have different reflective polarization spectral property.

 

Fig. 12 Polarization reflection distribution of aluminum under three wavelengths.

Download Full Size | PDF

4.3.4 Effect of polarization state of incidence

At last, we analyze the impact of the incident polarization state on reflective polarization. A two layer isotropic surface is considered here, which is a copper surface (refraction index:$1.31+1.85i$) with paint layer (refraction index:1.9). All of the layers have the same roughness parameters of $e_x=e_y=10$. The incident wavelength and the incident angle is set as 340nm and 40°, the threshold is set as $\tau =0.05$. Five incident polarization states, non-polarized light of $S_{in1}=\left[1000\right]$, linear polarized light of $S_{in2}=\left[1100\right]$, 45 degree linear polarized light of $S_{in3}=\left[1010\right]$, right and left circular polarized light of $S_{in4}=\left[1001\right]$and $S_{in5}=\left[100-1\right]$, are selected as the incidenct light respectively for our simulations.

The results from the incidence of $S_{in1}$ to $S_{in5}$are showed in Fig. 13. Compared to the results of unpolarized light $S_{in1}$, the reflected lights from polarized light show much stronger polarization and higher Stokes V parameter value. The angle of polarization patterns are different with variations of polarization states of incident lights. And the angle of polarization and Stokes Q,U,V parameters of the reflected lights from $S_{in4}=\left[1001\right]$ and $S_{in5}=\left[100-1\right]$are completely reversed with each other. Overall, it is apparent that reflective polarization varies with the incident polarization states, and if the incident light is polarized, there will be more polarized lights in the obtained reflective patterns, which will be a significant guide for the target detection and remote sensing.

 

Fig. 13 The polarization reflection of the incidence with different polarization state.

Download Full Size | PDF

5. Conclusion

Polarization is of crucial importance in a variety of applications, such as target detection and remote sensing. In this paper, we have proposed a Monte Carlo approach to acquire the reflective polarization information from an arbitrary multi-layer surface. The simulated experiments show that the results of our model matches well with the reported measured data and are more precise than the results from analytical BRDF model in SCATMECH. Our model not only satisfy the criterion of energy conservation, but also is suitable for the surface with arbitrary number of layers and profile. We also have also analyzed the influences of surface layer number, the surface roughness, the incident wavelengths and the polarization states on the reflective polarization distributions from the sample surfaces. The simulation results demonstrate that (1) as the layer number increases, the reflected light get more depolarized and the reflective intensity decrease, while the angle of polarization property is well kept; (2) the surface geometry change the shape of the reflective polarization pattern; (3) the surface with different material has different reflective polarization spectral property; (4) reflective polarization varies with the incident polarization states. All of these results will be very significant for the techniques of the target detection and remote sensing.