High-similarity analytical model of skylight polarization pattern based on position variations of neutral points

Peifeng Pan; Xin Wang; Xin Wang; Tian Yang; Xiankun Pu; Wenli Wang; Changhao Bao; Jun Gao

doi:10.1364/OE.489534

1. Introduction

The skylight polarization pattern is an important natural property of Earth’s atmosphere. Sunlight transmitted in the atmosphere is affected by the absorption and scattering of gas molecules, aerosol particles, and cloud particles. Part of the sunlight becomes polarized light, which forms a stable polarized light pattern throughout the sky [1]. The skylight polarization pattern contains rich information regarding a target and its environment, such as surface albedo [2], aerosols [3,4], clouds [5,6], and the variation of wavelength [7–10]. Biological studies have shown that desert ants [11], locusts [12], mantis shrimp [13], honeybees [14] and other organisms use polarized light to navigate and obtain location information. Then, scholars delved more deeply into the study of the changes in sky polarization pattern under different light and sky conditions, such as situations under thick clouds [15,16], situations under moonlit clear night sky [17], situations under twilight sky [18], and situations during solar eclipse [19]. Recently, it has been widely used in navigation [20], meteorological monitoring [21] and polarization remote sensing [22,23].

Two methods have been used to characterize skylight polarization patterns. The first is to use the atmospheric radiative transfer theory to solve the vector radiative transfer equation. However, this method is relatively complicated to calculate, requires a great deal of computing power and time, and has relatively high costs. The discrete coordinate [24], successive scattering [25], spherical harmonic [26], and Monte Carlo [27] methods are examples of this approach that can be used to solve the radiation transfer equation. The second approach is an analytical model based on the distribution of the skylight polarization pattern that can fully describe the distribution, variation, and symmetry of the ideal skylight polarization pattern [15]. For example, the Rayleigh model establishes a standard analytical model by placing the position of the neutral point at the sun and anti-solar positions [28].

The neutral point is one of the most important features of a skylight polarization pattern. Its position is closely related to that of the sun, which can provide accurate positioning information for bionic polarized light navigation. In recent years, many scientists have proven that there are four neutral points in the solar meridian direction: the Brewster point below the sun, Babinet point above the sun, and the Arago points above [29] and below the anti-sun positions. The Arago point below the anti-sun position is called the second Brewster point [30]. Thus, determining the locations of two or more of these neutral points provides information on the axis of symmetry of the polarization pattern, which can be utilized for navigation purposes. In previous studies, navigation algorithms that leverage the axis of symmetry or the positions of neutral points have already been proposed [31,32]. Dennis and Berry provided an analytical model of skylight polarization patterns based on the research of these four neutral points [33]. The sun was placed in the center between the Babinet and Brewster points. Additionally, the second Brewster and Arago points were placed in a position symmetrical to the Babinet and Brewster points. To date, the similarity between these two analytical models and the measured data of skylight polarization patterns has not been sufficiently high. One possible reason for this is that there is no parameter for variations in the neutral point positions. Furthermore, the fisheye lens imaging mode must also be considered. Therefore, this study proposes a new analytical model for this situation. This study uses the polarization angle to characterize the skylight polarization pattern because the stability of the angle of polarization (AoP) is greater than that of the degree of polarization (DoP). As such, the AoP is more conducive to navigation.

It has been found from measured data that the position of the neutral point is closely related to the position of the solar elevation angle. Wang et al. initially found the offset of the split angle between the neutral point and incident direction of the sun through a Monte Carlo simulation, which provided a new idea for the improvement of the skylight polarization pattern [34]. However, this approach lacked a comparison with the measured data. Sun et al. presented an improved model of the skylight polarization pattern that introduced the principle of fisheye lens imaging, which provided a method for comparing the analytical model with measured data [35]. In this study, a new analytical model of skylight polarization patterns is proposed by considering the relation between the neutral point and position of the sun. The Monte Carlo simulation method was used to initially determine this relationship. Subsequently, the relation was corrected using actual data to establish an analytical model of the skylight polarization pattern based on the variation of the solar altitude angle. Next, the model was compared with the original analytical model, which proved that the simulation data of the analytical model modified by the proposed method were closer to the actual data. Finally, the model was compared with actual data over a longer time period, which verified the applicability of the model at different periods.

2. Method

The classical analytical model of skylight polarization patterns and its shortcomings are introduced in this section. Subsequently, the analytical model is modified using simulation data and a new analytical model for the skylight polarization pattern is proposed.

The three-dimensional coordinate system of the skylight polarization pattern is shown in Fig. 1.

Fig. 1. Three-dimensional coordinate system diagram of skylight polarization pattern. The origin O represents the observer's position, the Z axis represents the zenith direction of the observation point, the X axis represents the east direction, the Y axis represents the north direction, and S represents the position of the sun. δ represents the elevation angle of P. $\varphi $ represents the azimuth angle of P. The solar elevation angle is ${\delta _s}$, and solar azimuth is ${\varphi _s}$. $\theta $ represents the scattering angle.

Download Full Size | PDF

$P({r,\delta ,\theta } )$ on the sphere is the position observed by the observer at any point. The spatial position of the sun is denoted by $S({r,{\delta_s},{\theta_s}} )$.

2.1 Rayleigh model

The blue color of the sky results from molecular scattering and absorption of short-wavelength light by ozone in the atmosphere. It is assumed that the scattered particles in the atmosphere are homogeneous isotropic spheres that are small compared to the wavelength of the incident light. Therefore, all the scattering characteristics of the fixed component in the Earth's atmosphere are within the range of the Rayleigh scattering theory. The skylight polarization pattern based on the Rayleigh scattering theory is called the Rayleigh scattering model [28]. In the Rayleigh model, the scattering angle θ is expressed as:

(1)$$cos \; \theta \textrm{ = cos}({cos \; \delta cos \; {\delta_s} + \sin \; \delta \sin \; {\delta_s}cos ({\varphi - {\varphi_s}} )} ). $$

The angle between the Electric field vector (E-vector) at a given measurement point and the local meridian (the line connecting the zenith to the sun) is defined as the angle of polarization (AoP). In the Rayleigh model it is expressed as:

(2)$$AoP = arccos \left( {\frac{{sin ({{\varphi_s} - \varphi } )}}{{sin \theta }} \times \left( {\frac{\pi }{2} - {\delta_s}} \right)} \right). $$

The traditional Rayleigh single-scattering model considers that the degrees of polarization of light in the directions of the sun and anti-sun are zero. However, the existence of multiple scattered sunlight splits the points (with a degree of polarization close to zero) into two. Therefore, this model is not accurate in expressing the skylight polarization pattern in the case of multiple scattering.

2.2 Singularity model

The Rayleigh scattering formula of particles was first considered to explain the reason for neutral-point splitting, which is expressed as:

(3)$$F(\theta )= \varDelta \left[ {\begin{array}{{cccc}} {\frac{3}{4}({1 + {{cos }^2}\theta } )}&{\frac{3}{4}{{sin }^2}\theta }&0&0\\ {\frac{3}{4}{{sin }^2}\theta }&{\frac{3}{4}({1 + {{cos }^2}\theta } )}&0&0\\ 0&0&{\frac{3}{2}{{cos }^2}\theta }&0\\ 0&0&0&{\varDelta^{\prime}\frac{3}{2}{{cos }^2}\theta } \end{array}} \right] + ({1 - \varDelta } )\left[ {\begin{array}{{cccc}} 1&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{array}} \right]$$

where ${\Delta } = ({1 - \tau } )/({1 + \tau /2} )$, ${\Delta ^{\prime}} = ({1 - 2\tau } )/({1 - \tau } )$, and $\tau $ is the depolarization coefficient of the particle. The Stokes vector of sunlight was set to ${I_{sun}} = [{I,0,0,0} ]$ to obtain the Stokes vector after single scattering, which is expressed as:

(4)$$\left\{ {\begin{aligned} &{Q = \frac{3}{4}{{sin }^2}\theta \cdot I \cdot \Delta }\\& {U = 0} \end{aligned}} \right..$$

It should be noted that this is the only Rayleigh scattering model obtained by single scattering. However, multiple scattering still occurred even though the atmosphere was dominated by single scattering. As such, only the secondary scattering process was calculated owing to computational complexity. It was assumed that the incident light Stokes vector before each last scattering was ${I_{sun}} = [{I,\eta I,0,0} ]$, where $\eta ={-} c\frac{3}{4}{\sin ^2}\theta \cdot {\Delta }$. c represents the proportion of secondary scattering in the total atmospheric polarization scattering. The degree of polarization at the neutral point was 0 based on the characteristics of the neutral point of the skylight polarization pattern. This can be obtained by substituting Eq. (3) with:

(5)$$Q = \Delta \frac{3}{4}{sin ^2}\theta I - \Delta \frac{3}{4}({1 + {{cos }^2}\theta } )\eta I = {sin ^2}\theta -{-} \frac{3}{4}c({1 + {{cos }^2}\theta } ){sin ^2}\theta \cdot \Delta = 1 -{-} \frac{3}{4}c({1 + {{cos }^2}\theta } )\Delta = 0$$

Two scattering angle positions that are symmetrical about the incident light can be obtained by solving Eq. (5). Subsequently, the sun was placed at the center of the Babinet and Brewster points using this characteristic in the singularity theoretical model established by Dennis and Berry [33]. The state of the polarized light in this singularity theory model is complex and can be expressed as:

(6)$$\omega (\xi )= {({{E_x} + i{E_y}} )^2} = E_x^2 - E_y^2 + 2i{E_x}{E_y} = |{\omega (\xi )} |{e^{2i\gamma (\xi )}}, $$

where $|{\omega (\xi )} |$ represents the DoP and $\gamma (\xi )$ represents the AoP. $\xi $ in the level coordinate system can be expressed as:

(7)$$\xi = r{e^{2i - \left( {\frac{1}{2}\varphi } \right)}}, $$

where

(8)$$r = ({1 - tan ({{\alpha_s}/2} )} )/({1 + tan ({{\alpha_s}/2} )} ). $$

The angular separation between two adjacent neutral points was ${\Delta }\delta = 4{\tan ^{ - 1}}(A )$, where A is a constant used to describe the projection of the two singularities. The positions of the four neutral points are expressed as:

(9)$${\xi _ + } = i\frac{{({{y_s} + A} )}}{{({1 - A{y_s}} )}},{\xi _ - } = i\frac{{({{y_s} - A} )}}{{({1 - A{y_s}} )}}, - 1/\xi _ + ^\ast , - 1/\xi _ - ^\ast , $$

where ${y_s} = r \cdot \cos {\delta _s}$ represents the distance from the projection of the sun’s position in the XOY coordinate to the origin. The following was derived from these four neutral points:

(10)$$\omega (\xi )\textrm{ = } - \frac{{4({\xi - {\xi_ + }} )({\xi - {\xi_ - }} )({\xi \textrm{ + 1/}\xi_ +^\ast } )({\xi + \xi_ -^\ast } )}}{{{{({1 + {r^2}} )}^2}|{{\xi_ + } + \textrm{1/}\xi_ +^\ast } ||{{\xi_ - } + \textrm{1/}\xi_ -^\ast } |}}. $$

2.3 Extension to the singularity model

The singularity theory model established by Dennis and Berry defaulted to a fixed angle between two adjacent neutral points, with the sun at the center of the angle of the two neutral points. However, it was found that the position of each neutral point changed with a shift in the solar position for the actual measurements and simulations. Specifically, each neutral point was on the meridian, and the angle between the neutral point and sun changed with the solar elevation angle. The angle between the Babinet point, Brewster point, and sun is described as ${\delta _ + },\; \; {\delta _ - }$, as shown in Fig. 2. The Arago, Babinet, Brewster, and second Brewster points were all symmetrical around the origin O.

Fig. 2. Neutral point diagram of the skylight polarization pattern.

Download Full Size | PDF

The Babinet, Brewster, second Brewster, and Arago points can be obtained using:

(11)$${\mu _ + } = i\frac{{({{y_s} + {A_ + }} )}}{{({1 - {A_ + }{y_s}} )}},{\mu _ - } = i\frac{{({{y_s} - {A_ - }} )}}{{({1 - {A_ - }{y_s}} )}}, - 1/\mu _ + ^\ast , - 1/\mu _ - ^\ast , $$

where ${A_ + } = \tan ({{\delta_ + }/2} )$ and ${A_ - } = \tan ({{\delta_ - }/2} )$. In addition, the model places the solar position on the Y-axis without considering changes in the solar azimuth. Therefore, the X-axis is considered as the zero-direction position and the solar azimuth change parameter is introduced. This model can be expressed as:

(12)$${\mu _ + } = \frac{{({{y_s} + {A_ + }} )}}{{({1 - {A_ + }{y_s}} )}}({cos ({{\varphi_s}} )+ i\;sin ({{\varphi_s}} )} ),{\mu _ - } = \frac{{({{y_s} - {A_ - }} )}}{{({1 - {A_ - }{y_s}} )}}({cos ({{\varphi_s}} )+ i\;sin ({{\varphi_s}} )} ), - 1/\mu _ + ^\ast , - 1/\mu _ - ^\ast$$

Then $\omega (\mu )$ turns to:

(13)$$\omega (\mu )\textrm{ = } - \frac{{4({\mu - {\mu_ + }} )({\mu - {\mu_ - }} )({\mu \textrm{ + 1/}\mu_ +^\ast } )({\mu + \mu_ -^\ast } )}}{{{{({1 + {r^2}} )}^2}|{{\mu_ + } + \textrm{1/}\mu_ +^\ast } ||{{\mu_ - } + \textrm{1/}\mu_ -^\ast } |}}. $$

DoP and AoP can be derived as:

(14)$$\textrm{DoP = }\left|{ - \frac{{4({\mu - {\mu_ + }} )({\mu - {\mu_ - }} )({\mu \textrm{ + 1/}\mu_ +^\ast } )({\mu + \mu_ -^\ast } )}}{{{{({1 + {r^2}} )}^2}|{{\mu_ + } + \textrm{1/}\mu_ +^\ast } ||{{\mu_ - } + \textrm{1/}\mu_ -^\ast } |}}} \right|, $$

(15)$$\textrm{AoP = }\frac{{\ln ({\omega (\mu )/|{\omega (\mu )} |} )}}{{2i}}. $$

Here, the change parameter of the neutral point is added to expand the standard singularity model further. Next, the change rule for the neutral point must be determined to improve the model.

2.4 Simulation experiment based on Monte Carlo method

The Monte Carlo method can be used to simulate the skylight polarization pattern at different solar elevation angles and determine the relation between the neutral point of the skylight polarization pattern and solar elevation angle [27]. The Monte Carlo method achieves high accuracy according to the simulation results of Ramella-Roman et al. [36]. The Monte Carlo method was used in this study to simulate the distribution of skylight polarization patterns at different altitude angles. In particular, the flattened atmosphere model was used that assumes all particles are Rayleigh particles. The height angle step was set to 1° and the azimuth angle was set to 0°. The wavelength of the incident light was set to 450 nm. The Earth’s curvature was taken into account through the use of the Three-dimension RTE solver, MYSTIC. Originally developed as a forward tracing method for the calculation of irradiances and radiances in plane-parallel atmospheres, MYSTIC has since been extended to incorporate fully spherical geometry and a backward tracing mode. This enables calculations to be performed in a 1D spherical model atmosphere. The skylight polarization pattern simulated using the Monte Carlo method is shown in Fig. 3.

Fig. 3. Skylight polarization patterns at solar elevation angles of 60°, 45°, and 30°: (a) the AoP of the skylight polarization pattern, (b) DoP, (c) contour distribution of AoP, and (d) contour distribution of DoP. S represents the position of the sun.

Download Full Size | PDF

The position of the neutral point of the skylight polarization pattern clearly changed with a shift in the solar position, as seen in Fig. 3. Additionally, the distance between the neutral point and sun changed with a shift in the solar elevation angle. It is known that the position of the neutral point has the lowest degree of polarization in the skylight polarization pattern. Therefore, the exact position of the neutral point can be determined using the center position of the contour line with the smallest degree of polarization. The position of the neutral point clearly belongs to the intersection position of the $\infty $-type image of the polarization angle of the skylight polarization pattern, facilitating the accurate location of the neutral point. The ${\delta _ + }$ and ${\delta _ - }$ of the neutral points were calculated at the Babinet and Brewster points at the solar elevation angle ranges of 30° to 60°. The relation between the position of the neutral point and solar elevation angle in the simulation experiment is shown in Fig. 4.

Fig. 4. Relations between the positions of the neutral point and solar elevation angle in the simulation experiment. The left and right pictures illustrate the relations between ${\delta _ + }$ and ${\delta _ - }$ with the solar elevation angle, respectively.

Download Full Size | PDF

This set of simulation experiments demonstrated that the neutral point position of the skylight polarization pattern was linearly related to the solar elevation angle. The Babinet and Brewster points relative to the sun (${\delta _ + }$ and ${\delta _ - }$) can be expressed as functions of the solar elevation angle ${\delta _s}$, which is expressed as:

(16)$$\left\{ {\begin{array}{{c}} {{A_ + } = tan [{{\delta_ + }({{\delta_s}} )/2} ]}\\ {{A_ - } = tan [{{\delta_ - }({{\delta_s}} )/2} ]} \end{array}} \right.. $$

The linear relation between the neutral point and solar elevation angle was initially determined through a preliminary Monte Carlo simulation, which laid the foundation for the study of the measured data. However, the simulated skylight polarization pattern still had a large gap compared to the actual pattern because the simulation experiment could not fully simulate the measured environment. Therefore, a measured experiment was used to improve and correct the analytical model of the skylight polarization pattern.

2.5 Polar coordinate calculation in camera images

The polar coordinate system was used in the analytical model to represent the positions of the sun S and observation point P. However, the results obtained from the fisheye lens used in the actual measurement experiment were in the form of a matrix. Therefore, it was necessary to convert the images captured by the camera from rectangular to polar coordinates. To achieve improved consistency between our simulation model and the measured data, our previous work evaluated the impacts of various types of fisheye lenses. In Ref. [35], the potential variability of the appropriate imaging mode for a fisheye lens was explored. Therefore, according to the optical design of the fisheye lens used in this current skylight polarization pattern acquisition setup, we made the compensation based on the equidistance imaging mode. In the coordinate system shown in Fig. 2, the azimuth of the observation point equidistant imaging mode is expressed as:

(17)$$\varphi \textrm{ = arctan}({\textrm{y/x}} )\qquad \quad\textrm{}({x \ge 0,y \ge 0} ).$$

Subsequently, the following can be obtained because the equidistant imaging mode can be expressed as ${y_0} = f\omega $:

(18)$$\delta = \frac{\pi }{2}\sqrt {{x^2} + {y^2}} . $$

According to the above equations, the θ and φ corresponding to each pixel point on the image taken by the wide-angle lens of the equidistant imaging mode can be calculated to establish the correspondence between the pixel points in the camera and the three-dimensional coordinate system of the skylight polarization pattern. In addition, the output data in the analytical model can be converted to the camera coordinate system. The zenith point was at the center of the matrix (n/2, n/2) because measured and simulation data were collected in an n × n square matrix in the comparison experiment. According to the MatLab matrix program:

(19)$$\delta = \omega \sqrt {{{({x - n/2} )}^2} + {{({y - n/2} )}^2}},$$

and

(20)$$\begin{cases} \varphi = \arctan \left( {\frac{{x - n/2}}{{n/2 - y}}} \right)\qquad\qquad({x \ge n/2,y\, \le n/2} )\\ \varphi ={-} \arctan \left( {\frac{{x - n/2}}{{y - n/2}}} \right)\qquad\qquad({x \ge n/2,y\, \ge n/2} )\\ \varphi ={-} \arctan \left( {\frac{{n/2 - x}}{{y - n/2}}} \right)\, - \pi \quad({x \le n/2,y\, \ge n/2} )\\ \varphi = \pi - \arctan \left( {\frac{{n/2 - x}}{{n/2 - y}}} \right)\qquad({x \le n/2,y\, \le n/2} )\end{cases},$$

where $\omega $ represents the field of the camera view. The matrix coordinate can subsequently be substituted from Eq. (19) into Eq. (20). In this manner, a comparison between the measured and simulation data can be realized. The size of the picture for data collected in the measured experiment was $4001 \times 4001\; \textrm{p}$, thus the value of n was 4001. The flowchart of the proposed model can be illustrated as Fig. 5.

Fig. 5. The flowchart of the proposed model.

Download Full Size | PDF

3. Experiments and discussion

3.1 Processing of measured data

The skylight polarization pattern imaging system used in this study is shown in Fig. 6.

Fig. 6. Acquisition device of skylight polarization pattern measured data.

Download Full Size | PDF

This device was mainly composed of a Nikon D850 camera, Sigma 8 mm F/3.5 fisheye lens, polarizer, and wireless shutter. The fisheye lens was used for this purpose as it was capable of collecting sky polarization information over a wide field of view. Thus the location of the neutral points can be determined by measuring their altitude and azimuth. The setup has been well-established in studies of neutral point observation, as described in previous works [29,31]. In addition, a specialized apparatus was designed for installing the polarizer in front of the fisheye lens, which enables manual rotation of the polarizer to obtain polarization information in different directions, and clearly indicates the polarizer's rotation angle. The duration required for rotating the polarizer to acquire the four necessary images is approximately 6 seconds. In general, a sun occulter should be set to avoid ghost spots. However, this study aimed at establishing the relationship between the variations in the neutral point position and the solar altitude angle. It is worth noting that the determination of the neutral points in the sky is crucial for deriving this relationship. The use of a sun occulter would result in the obscuration of the neutral points. In light of this, the sun occulter was not included in our measurement setup. The field of this device was 140°, and it could collect the light intensity in the four directions of 0°, 45°, 90°, and 135° by manually rotating the linear polarizer. These directions could subsequently be used to calculate the Stokes vector of the skylight polarization pattern, $S = [{I,Q,U,V} ]$, and the value of V could be neglected, which is expressed as [37]:

(21)$$\left\{ {\begin{aligned}& {I = \frac{1}{2}({{I_0} + {I_{45}} + {I_{90}} + {I_{135}}} )}\\& {Q = {I_0} - {I_{90}}}\\& {U = {I_{\textrm{45}}} - {I_{135}}}\\& {V = 0} \end{aligned}} \right..$$

The AoP and DoP can be expressed as:

(22)$$\left\{ {\begin{array}{{c}} {AoP = \frac{1}{2}{{tan }^{ - 1}}\frac{Q}{U}}\\ {DoP = \sqrt {\frac{{{Q^2} + {U^2}}}{{{I^2}}}} } \end{array}} \right.. $$

To present AoP as the angle between the E-vector at the given measurement point and the local meridian:

(23)$$AoP({x,y} )= AoP({x,y} )- {\tan ^{ - 1}}[{({y - n/2} )/({n/2 - x} )} ], $$

then

(24)$$AoP({x,y} )= \left\{ {\begin{array}{{c}} {AoP({x,y} )+ 180\; \; \; ({AOP({x,y} )< - 90} )}\\ {AoP({x,y} )- 180\; \; \; \; \; \; \; ({AOP({x,y} )> 90} )} \end{array}} \right.. $$

The experiment was performed in the gymnasium of the Feicuihu campus at the Hefei University of Technology (longitude east $117^\circ 17^{\prime}43^{\prime\prime}$, latitude north $31^\circ 50^{\prime}49^{\prime\prime}$). The experiment was conducted from 13:50 to 19:00 on August 5, 2022. The weather was bright and cloudless. The data were measured every ten minutes. Some of the measurement results are shown in Fig. 7.

Fig. 7. Measurement results from August 5, 2022: (a) the original image, (b) AoP of measured data, and (c) DoP of measured data.

Download Full Size | PDF

From the above measurement results, the AoP and DoP diagrams of the skylight polarization pattern can be used to determine the neutral points [31,32]. The matrix coordinates of the neutral points can be transformed into polar coordinates by calculating Eqs. (19) and (20), as shown in Fig. 8.

Fig. 8. Relations between the positions of the neutral points and solar elevation angle in the measured experiment. The left and right pictures represent the relations between ${\delta _ + }$ and ${\delta _ - }$ with the solar elevation angle, respectively.

Download Full Size | PDF

${\delta _ + }$ decreased in the measured data with an increase in the solar altitude angle in the range of 0° to 64°. ${\delta _ - }$ increased with an increase in the solar altitude angle from 0° to 27°. However, ${\delta _ + }$ decreased with an increase in the solar altitude angle in the range of 27° to 64°. ${\delta _ + }$ and δ − had a linear relationship with the solar elevation angle according to the Monte Carlo simulation results. Therefore, ${\delta _ + }$ and ${\delta _ - }$ were linearly fitted with the solar elevation angle. We attributed the lack of a good linear relationship at lower solar elevation angles to the fact that, when the sun is at a low angle, the amount of light entering the measurement device form the sun is greatly affected by the stray light from the ground, leading to disparities between the measured data and theoretical predictions. And it was common to this dataset. The relationship between the neutral point position and solar elevation angle can be expressed as:

(25)$$\left\{ {\begin{array}{{c}} {{\delta_ - }({{\delta_s}} )\textrm{ = 56}\textrm{.84} - \textrm{0}\textrm{.25} \times {\delta_s}\qquad\textrm{}({{\delta_s} \ge \textrm{27}} )}\\ {{\delta_ - }({{\delta_s}} )\textrm{ = 37}\textrm{.34 + 0}\textrm{.49} \times {\delta_s}\qquad\textrm{}({{\delta_s} \le \textrm{27}} )} \end{array}} \right.,$$

and

(26)$${\delta _ + }({{\delta_s}} )\textrm{ = 42}\textrm{.53} - \textrm{0}\textrm{.56} \times {\delta _s}. $$

The direct relationship of the angle between the neutral point and sun and the solar altitude angle under a clear sky in the summer is presented in Eqs. (25) and (26). Based on this characteristic, the analytical model of the skylight polarization pattern of the sun can fit at any position under a clear sky. It is worth noting that we have investigated the influence of wavelength in a previous paper [38], and constructed a model to incorporate the wavelength factor as a variable. However, in this study, results from measurements indicated that the influence of wavelength on the accuracy of the relationship between neutral point positions and solar altitude angle simulated by the proposed model is not significant, and thus it was not included here.

3.2 Comparison of models

The size of the image collected in the measured experiment was $4001 \times 4001$ p, that is, the value of n in Eq. (19) and Eq. (20) was 4001. The method proposed by Gabor was used to compare the simulated and measured data to verify the validity of the model [37]. First, $|{Ao{p_1}({i,j} )- Ao{p_2}({i,j} )} |= \varDelta Aop({i,j} )$ was calculated based on a point of identical positions in the measured and simulation data. These points were considered as similar pairwise data when $\varDelta Aop({i,j} )\le 5$. ${N_{Aop}}$ was used to denote the number of similar pairwise data. ${\varepsilon _{Aop}} = {N_{Aop}}/N$ was defined as an index of similarity between the simulation and measured data. The location and shooting time were used to calculate the position of the sun [39]. The Rayleigh, standard singularity, and proposed models were simulated and compared according to this method, as shown in Fig. 9.

Fig. 9. Comparison between the AoPs in the measured skylight polarization pattern and AoPs of the: (a) measured data, (b) proposed model, (c) Rayleigh model, and (d) singularity model.

Download Full Size | PDF

The difference between each model and the measured data can be described with $\varDelta AoP$ images, as shown in Fig. 10.

Fig. 10. Difference in skylight polarization patterns between the measured data and the: (a) proposed model, (b) Rayleigh model, and (c) singularity model.

Download Full Size | PDF

Nine sets of experiments were performed to compare the measured and simulated data. The observation times of the measured data were 15:10, 15:30, 16:00, 16:30, 17:00, 17:30, 18:00, 18:30, and 19:00. The similarity between the simulation and measured data was obviously improved in the proposed model, as shown in Figs. 9 and 10. The similarities between the simulation and measured data in various analytical models are listed in Table 1 and illustrated in Fig. 11.

Fig. 11. Similarity between simulation and measured data in different models.

Download Full Size | PDF

Table 1. Similarity between simulation and measured data in different models

View Table

According to the comparison of the three models, the similarity between the proposed skylight polarization pattern model and the actual measurement was higher than those of the traditional Rayleigh and standard singular point models. This clearly proves the advantages of the proposed model.

3.3 Comparison of models over different time periods

The accuracy of the model for different periods was also tested. The data from the morning of August 5, 2022, were compared using the simulation data of this model, as shown in Fig. 12.

Fig. 12. Comparison between the morning simulation and measured data for: (a) the AoP of measured data, (b) AoP of the proposed model, and (c) difference between the proposed model and measured data.

Download Full Size | PDF

Then, the measured data from 15:30 to 18:00 were compared with the simulation data of the model, as shown in Fig. 13.

Fig. 13. Comparison between simulation and measured data: (a) the AoP of measured data, (b) AoP of the proposed model, and (c) difference between the proposed model and measured data.

Download Full Size | PDF

Furthermore, data were collected over a longer time range, including June, July, and August, and compared with those of the proposed model to calculate their similarity. The data were collected on June 29, July 24, and August 21. The statistical results are as presented in Fig. 14.

Fig. 14. Similarity statistics between simulated data from the proposed model and measured data in different months.

Download Full Size | PDF

All measured data were taken from the gymnasium of the Feicuihu campus of the Hefei University of Technology (longitude east 117°17'43'‘, latitude north 31°50'49'‘). The proposed model can be used in many periods with high similarity, as shown in Figs. 12 to 14. This proves that the proposed model achieved strong universality. It can be noticed that there was better similarity for August than June and July. We consider that the similarities between the measured data and simulation results for August, June and July did not have an essential difference. It is worth noting that the measurement of the sky polarization pattern is susceptible to the weather conditions. Although all the experiments were conducted on clear days, it was not feasible to guarantee that all weather conditions would be identical. The more favorable weather conditions on August 21 compared to June 29 and July 24 likely contributed to these outcomes.

4. Conclusion

This study proposed a high similarity analytical model of a skylight polarization pattern based on position variations of the neutral points. This model greatly improved the similarity between simulation and measurement data to reflect the actual situation more accurately. First, it was confirmed through a Monte Carlo simulation that there were linear relations between the positions of the neutral points and the solar altitude angle. Second, a function was built for the neutral point and solar altitude angle using a large number of measured data. A high similarity analytical model for the skylight polarization pattern was achieved by introducing this function and the influence of the imaging mode of the fisheye camera. The proposed model was proven to attain a higher similarity than the traditional method by comparing various models with the measured data. Finally, the proposed analysis model was proven to have universality based on its high similarity when comparing the measured data for three consecutive months. This study provides a more accurate reference for the application of skylight polarization patterns in practical scenarios.

Funding

National Natural Science Foundation of China (62171178, 61801161, 61971177).

Acknowledgments

The authors thank Daqian Wang for the valuable discussions and experimental data collection.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data supporting the findings of this study in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. G. S. Smith, “The polarization of skylight: An example from nature,” Am. J. Phys. 75(1), 25–35 (2007). [CrossRef]

2. A. R. Dahlberg, N. J. Pust, and J. A. Shaw, “Effects of surface reflectance on skylight polarization measurements at the Mauna Loa Observatory,” Opt. Express 19(17), 16008–16021 (2011). [CrossRef]

3. W. G. Egan, “Optical sky polarimetry and photometry at Mauna Loa Observatory affected by stratospheric dust from El Chichon,” Appl. Opt. 23(7), 1013–1020 (1984). [CrossRef]

4. A. Kreuter, C. Emde, and M. Blumthaler, “Measuring the influence of aerosols and albedo on sky polarization,” Atmos. Res. 98(2-4), 363–367 (2010). [CrossRef]

5. N. J. Pust and J. A. Shaw, “Digital all-sky polarization imaging of partly cloudy skies,” Appl. Opt. 47(34), H190–H198 (2008). [CrossRef]

6. L. M. Eshelman, M. J. Tauc, and J. A. Shaw, “All-sky polarization imaging of cloud thermodynamic phase,” Opt. Express 27(3), 3528–3541 (2019). [CrossRef]

7. L. M. Eshelman and J. A. Shaw, “The VIS-SWIR spectrum of skylight polarization,” Appl. Opt. 57(27), 7974–7986 (2018). [CrossRef]

8. R. L. Lee and O. R. Samudio, “Spectral polarization of clear and hazy coastal skies,” Appl. Opt. 51(31), 7499–7508 (2012). [CrossRef]

9. N. J. Pust and J. A. Shaw, “Wavelength dependence of the degree of polarization in cloud-free skies: simulations of real environments,” Opt. Express 20(14), 15559–15568 (2012). [CrossRef]

10. G. V. Harten, F. Snik, J. H. H. Rietjens, J. M. Smit, and C. U. Keller, “Spectral line polarimetry with a channeled polarimeter,” Appl. Opt. 53(19), 4187–4194 (2014). [CrossRef]

11. R. Wehner and S. Wehner, “Path integration in desert ants. Approaching a long-standing puzzle in insect navigation,” Ital. J. Zool. 20(3), 309–331 (1986).

12. M. Bech, U. Homberg, and K. Pfeiffer, “Receptive fields of locust brain neurons are matched to polarization patterns of the sky,” Curr. Biol. 24(18), 2124–2129 (2014). [CrossRef]

13. M. Garcia, T. Davis, S. Blair, N. Cui, and V. Gruev, “Bioinspired polarization imager with high dynamic range,” Optica 5(10), 1240–1246 (2018). [CrossRef]

14. M. V. Srinivasan, “Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics,” Physiol. Rev. 91(2), 413–460 (2011). [CrossRef]

15. R. Hegedüs, S. Åkesson, and G. Horváth, “Polarization patterns of thick clouds: overcast skies have distribution of the angle of polarization similar to that of clear skies,” J. Opt. Soc. Am. A 24(8), 2347–2356 (2007). [CrossRef]

16. I. Pomozi, G. Horváth, and R. Wehner, “How the clear-sky angle of polarization pattern continues underneath clouds: full-sky measurements and implications for animal orientation,” J. Exp. Biol. 204(17), 2933–2942 (2001). [CrossRef]

17. J. Gál, G. Horváth, A. Barta, and R. Wehner, “Polarization of the moonlit clear night sky measured by full-sky imaging polarimetry at full Moon: Comparison of the polarization of moonlit and sunlit skies,” J. Geophys. Res. 106(D19), 22647–22653 (2001). [CrossRef]

18. T. W. Cronin, E. J Warrant, and B. Greiner, “Polarization patterns of the twilight sky,” Polarization Science and Remote Sensing II. SPIE 5888, 58880R (2005). [CrossRef]

19. J. A. Shaw, L. M. Eshelman, M. J. Tauc, and G. E. Shaw, “Temporal evolution of sky polarization during solar eclipse totality,” Proc. SPIE 11132, 111320C (2019). [CrossRef]

20. H. Lu, K. Zhao, Z. You, and K. Huang, “Angle algorithm based on Hough transform for imaging polarization navigation sensor,” Opt. Express 23(6), 7248–7262 (2015). [CrossRef]

21. G. Zhang, V. N. Mahale, B. J. Putnam, Y. Qi, Q. Cao, A. D. Byrd, and R. J. Doviak, “Current status and future challenges of weather radar polarimetry: Bridging the gap between radar meteorology/hydrology/engineering and numerical weather prediction,” Adv. Atmos. Sci. 36(6), 571–588 (2019). [CrossRef]

22. L. Yan, Y. Li, V. Chandrasekar, H. Mortimer, J. Peltoniemi, and Y. Lin, “General review of optical polarization remote sensing,” Int. J. Remote Sens. 41(13), 4853–4864 (2020). [CrossRef]

23. T. Tanikawa, K. Masuda, and H. Ishimoto, “Spectral degree of linear polarization and neutral points of polarization in snow and ice surfaces,” J. Quant. Spectrosc. Radiat. Transfer 273, 107845 (2021). [CrossRef]

24. F. M. Schulz and K. Stamnes, “Angular distribution of the Stokes vector in a plane-parallel, vertically inhomogeneous medium in the vector discrete ordinate radiative transfer (VDISORT) model,” J. Quant. Spectrosc. Radiat. Transfer 65(4), 609–620 (2000). [CrossRef]

25. Q. Min and M. Duan, “A successive order of scattering model for solving vector radiative transfer in the atmosphere,” J. Quant. Spectrosc. Radiat. Transfer 87(3-4), 243–259 (2004). [CrossRef]

26. K. F. Evans, “The spherical harmonics discrete ordinate method for three-dimensional atmospheric radiative transfer,” J. Atmos. Sci. 55(3), 429–446 (1998). [CrossRef]

27. F. F. M. de Mul, “Monte-Carlo simulations of light scattering in turbid media,” Handbook of Coherent-Domain Optical Methods593–661 (Springer, 2013).

28. A. Wilkie, C. Ulbricht, R. F. Tobler, G. Zotti, and W. Purgathofer, “An analytical model for skylight polarization,” in 15th Eurographics Symposium on Rendering (pp. 387–398) (2004).

29. J. Gál, G. Horváth, V. B. Meyer-Rochow, and R. Wehner, “Polarization patterns of the summer sky and its neutral points measured by full–sky imaging polarimetry in Finnish Lapland north of the Arctic Circle,” Proc. R. Soc. London, Ser. A 457(2010), 1385–1399 (2001). [CrossRef]

30. G. Horváth, B. Bernáth, B. Suhai, A. Barta, and R. Wehner, “First observation of the fourth neutral polarization point in the atmosphere,” J. Opt. Soc. Am. A 19(10), 2085–2099 (2002). [CrossRef]

31. Z. Fan, X. Wang, H. Jin, C. Wang, N. Pan, and D. Hua, “Neutral point detection using the AOP of polarized skylight patterns,” Opt. Express 29(4), 5665–5676 (2021). [CrossRef]

32. H. Zhao, W. Xu, Y. Zhang, X. Li, H. Zhang, J. Xuan, and B. Jia, “Polarization patterns under different sky conditions and a navigation method based on the symmetry of the AOP map of skylight,” Opt. Express 26(22), 28589–28603 (2018). [CrossRef]

33. M. V. Berry, M. R. Dennis, and R. L. Lee, “Polarization singularities in the clear sky,” New J. Phys. 6(1), 162 (2004). [CrossRef]

34. W. Jue, Q. Jianqiang, and G. Lei, “An equivalent incident light model for multiple scattering and polarization neutral point of skylight,” J. Opt. (Bristol, U. K.) 22(7), 075603 (2020). [CrossRef]

35. S. Sun, J. Gao, D. Wang, T. Yang, and X. Wang, “Improved models of imaging of skylight polarization through a fisheye lens,” Sensors 19(22), 4844 (2019). [CrossRef]

36. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef]

37. T. Yang, X. Wang, X. Pu, Z. Shi, S. Sun, and J. Gao, “Adaptive method for estimating information from a polarized skylight,” Appl. Opt. 60(30), 9504–9511 (2021). [CrossRef]

38. X. Wang, J. Gao, Z. Fan, and N. W. Roberts, “An analytical model for the celestial distribution of polarized light, accounting for polarization singularities, wavelength and atmospheric turbidity,” J. Opt. (Bristol, U. K.) 18(6), 065601 (2016). [CrossRef]

39. I. Reda and A. Andreas, “Solar position algorithm for solar radiation applications,” J. Sol. Energy 76(5), 577–589 (2004). [CrossRef]

Time	Solar altitude angle (°)	Solar azimuth angle (°)	Proposed model	Rayleigh model	Singularity model
15:10	48.29	-14.84	54.81%	30.72%	23.83%
15:30	44.09	-11.54	62.07%	31.07%	22.84%
16:00	37.74	-7.13	73.34%	37.86%	27.96%
16:30	31.37	-3.13	80.14%	41.36%	33.54%
17:00	25.00	0.62	81.33%	45.62%	42.84%
17:30	18.68	4.26	82.51%	46.44%	53.24%
18:00	12.43	7.88	84.29%	48.35%	64.71%
18:30	6.30	11.58	84.25%	52.75%	73.95%
19:00	0.51	15.42	86.10%	55.03%	76.53%

High-similarity analytical model of skylight polarization pattern based on position variations of neutral points

Abstract

1. Introduction

2. Method

2.1 Rayleigh model

2.2 Singularity model

2.3 Extension to the singularity model

2.4 Simulation experiment based on Monte Carlo method

2.5 Polar coordinate calculation in camera images

3. Experiments and discussion

3.1 Processing of measured data

3.2 Comparison of models

3.3 Comparison of models over different time periods

4. Conclusion

Funding

Acknowledgments

Disclosures

Data availability

References

Data availability

Cited By

Figures (14)

Tables (1)

Equations (26)

Optics Express