PolarizationNavigation Simulation System and Skylight Compass Method Design Based upon Moment of Inertia

Unpolarized sunlight becomes polarized by atmospheric scattering and produces a skylight polarization pattern in the sky, which is detected for navigation by several species of insects. Inspired by these insects, a growing number of research studies have been conducted on how to effectively determine a heading angle from polarization patterns of skylight. However, few research studies have considered that the pixels of a pixelated polarization camera can be easily disturbed by noise and numerical values among adjacent pixels are discontinuous caused by crosstalk. So, this paper proposes a skylight compass method based upon the moment of inertia (MOI). Inspired by rigid body dynamics, the MOI of a rigid body with uniform mass distribution reaches the extreme values when the rigid body rotates on its symmetry axes. So, a whole polarization image is regarded as a rigid body. *en, orientation determination is transformed into solving the extreme value of MOI. *is method makes full use of the polarization information of a whole polarization image and accordingly reduces the influence of the numerical discontinuity among adjacent pixels and measurement noise. In addition, this has been verified by numerical simulation and experiment. And the compass error of the MOI method is less than 0.44° for an actual polarization image.


Introduction
Unpolarized sunlight becomes polarized after being scattered by gas molecules, dust, aerosols, and other scatters in the atmosphere and produces polarization patterns in the sky. It has been reported that several species of insects can sense the pattern of polarization of light in the sky as a compass. Desert ants are able to explore their desert habitat hundreds of meters away and return to their nest on a straight line by detecting the polarization of skylight [1][2][3][4]. Locusts exploit the polarization pattern to extract directional information from the sky [5][6][7]. Dung beetles can sense polarized light to roll dung piles away along a straight path [8][9][10]. Drosophilae detect the electric field vector (E-vector) of polarized skylight as a compass cue [11][12][13]. e premise of polarized skylight navigation is that insects have polarization-sensitive photoreceptors. Desert ants, locusts, and other insects have a specialized group of ommatidia in the dorsal rim area (DRA) of the compound eye (Figure 1(a)) [6,11,12]. In the DRA, ommatidial microvilli are arranged in different angles (Figure 1(b)) [7,14].
is physiological structure enables insects to detect polarized light. Inspired by this physiological structure, a pixelated polarized camera is designed [15][16][17][18], which has a four-directional polarizer and can capture a four-directional polarization image in one shot. In order to take advantage of polarized light for polarization navigation, a polarized skylight navigation platform has been constructed based upon this polarized camera in this paper, which will be described in Section 4.2.
To describe skylight polarization patterns, in 1870, Rayleigh proposed a mathematical theory to explain why the sky is blue, which is then known as the Rayleigh scattering theory [19,20]. e skylight polarization pattern has been well described according to the Rayleigh theory, even when only single-scatter behavior is considered [21,22]. As shown in Figure 1(c), the Rayleigh sky model relates to the position of the Sun. e E-vector of skylight is always perpendicular to the solar vector. Specially, the E-vector in the solar meridian and antisolar meridian (SM-ASM) is always perpendicular to SM-ASM. And the E-vector in the zenith is always perpendicular to the solar azimuth. In addition, the Rayleigh sky model is symmetry about SM-ASM. e Rayleigh sky model is simple and convenient for practical applications and has become a widely used model in bio-inspired polarized skylight navigation [23].
Several attitude estimation methods have been proposed so far based on the Rayleigh sky model. Art Lompado et al. have developed a practical application that accurately determined the azimuthal pointing direction of a platform by measuring the polarization at the zenith [23][24][25][26][27][28]. Moshe Hamaoui et al. have estimated the angel of the solar azimuth based on the gradient of the angle of polarization (AOP) [29,30]. Lu et al. have presented an angle algorithm based upon Hough transform to determine orientation [31,32]. Chen Fan et al. have proposed an algorithm with the total least square to calculate the solar vector [33][34][35]. However, a more comprehensive study should consider the fact that the pixels of a pixelated polarization camera can be easily disturbed by noise and the numerical values among adjacent pixels are discontinuous caused by crosstalk [15,16]. To address this technology gap, this paper proposes a new methodology for orientation determination based on symmetry axis extraction using the moment of inertia (MOI), which makes full use of the symmetry characteristics of the whole sky to determine orientation and accordingly reaches a high level of stability.
To summarize, in order to suppress the interference of measurement noise and numerical discontinuity among pixels of the pixelated polarization camera, in this paper, a skylight compass method based on MOI is proposed, which extracts the symmetry axis of the whole skylight polarization pattern to determine orientation. en, based on the Rayleigh sky model, a simulation system is constructed to validate the MOI method. Finally, a field experiment using a polarized skylight navigation platform is carried out.

Orientation Determination Algorithm Design
As shown in Figure 2, in rigid body dynamics, the MOI of a rigid body with uniform mass distribution reaches the extreme values when the rigid body rotates on its symmetry axes [36][37][38]. So, inspired by this theory, a whole polarization image can be regarded as a rigid body with uniform mass distribution and each pixel point of polarization image can be considered to be equivalent to each mass unit in the rigid body. us, the MOI of the polarization image can be calculated to extract the symmetry axes.
To simplify the description of the proposed method, a body coordinate system o b x b y b z b of the polarization camera is established. As shown in Figure 3, o b is the center of the polarization image, x b and y b are aligned with the direction of the column and row of the polarization image, and z b points to the shoot direction of the polarization camera.
For a rigid body, each mass unit is dm, and the distance to a line is r. e definition of MOI is In a grayscale image, the rigid body is composed of a set of pixel points P ′ � (x i , y i ), i ∈ [1, N] , where N is the number of pixels in the image. For polarized skylight navigation, P ′ is a circular region with o b as its origin. Each pixel point can be considered to be equivalent to each mass unit in the rigid body. e MOI of a rigid body relative to a  To ensure that all lines are considered, the equation of the axis of symmetry is defined as where a, b, and c are the coefficients of the axis of symmetry. When we use the polarized skylight navigation platform to navigate, the polarization camera points toward the zenith. erefore, the axis of symmetry would definitely pass through the center of the polarization image, so c � 0. en, the whole MOI of a polarization image is where is the distance between the pixel point (x i , y i ) and line ay + bx � 0 and μ i represents the mass density. For polarized skylight navigation, μ i is the absolute value of AOP, whose reference direction is the local meridian.
e MOI of a rigid body reaches the extreme values when the rigid body rotates on its symmetry axes. erefore, the partial derivatives of I should be calculated to extract symmetry axes. e derivative of I with respect to a is For equation (6), when /2G. So, the axes of symmetry can be determined by In short, two straight lines are obtained. And we have So, the two straight lines are perpendicular to each other. Only one of them is the SM-ASM. en, according to the polarization E-vector on SM-ASM perpendicular to SM-ASM, the final axis of symmetry can be determined. Above all, SM-ASM is determined based on MOI.

Polarization Navigation Simulation System
To validate the feasibility of the MOI algorithm, a simulation system is constructed, in which the polarization image can be captured by a hypothetical imager based on the Rayleigh sky model. Figure 4, ox g y g z g is the geographic coordinate system. e east is labeled x g , the north is labeled y g , and the up is labeled z g . e Rayleigh sky model predicts the degree of polarization (DOP) [24] as

Rayleigh Sky Model. As shown in
where DOP max denotes the maximum DOP of the whole sky and equals 1 for an ideal (single scattering) sky, c P is the angular distance between the Sun and the observation point P. e views of DOP are shown in Figures 5(a) and 5(b). e Rayleigh sky model predicts the AOP [31] as where θ S is the solar zenith angle, θ P is the zenith angle of P, φ S is the solar azimuth angle, and φ p is the zenith angle of P in the geography coordinate. e views of AOP are shown in Figures 5(c) and 5(d).
AOP is the angle between the polarization E-vector and reference direction. So, the E-vector is calculated, and then AOP can be determined according to the reference direction. e E-vector at the point P is where In equations (13) and (14), V gP �� �→ and H gP ���→ are tangential directions of the meridian and weft at the point P, respectively. e superscript T represents matrix transpose.
Substituting equations (10) and (11) into equation (12), components of E gP ��→ can be given by e views of the E-vector are shown in Figures 5(e) and 5(f ). Figure 4: e geographic coordinate system. ox g y g z g is the geographic coordinate system. e yellow point represents the Sun, and the pink point represents the observation point P. SM-ASM is the solar meridian and antisolar meridian corresponding to the red line. c P is the scattering angle of P. θ S is the solar zenith angle. h S is the solar elevation angle. θ P is the zenith angle of P. φ S is the solar azimuth angle. φ P is the azimuth angle of P.

Hypothetical Imager.
To capture the AOP and DOP images, DOP and AOP of each pixel are captured. In the body coordinate, the vector of the pixel P ′ is where D x and D y are the pixel sizes of the chargecoupled device (CCD) or complementary metal oxide semiconductor (CMOS) array, η x and η y represent the image has η x × η y pixels, and m x and m y represent the count of the column and row pixels of the pixel P ′ . Suppose the attitude of the polarization camera is given, and then the rotation matrix from the body to the geographic coordinate system is where φ represents the yaw angle, α represents the pitch angle, and β represents the roll angle.
According to the definition of the body coordinate, the shooting direction of the polarization camera is along z b → � (0, 0, 1) T . us, in the geographic coordinate system, the shooting direction of the pixel P ′ is where f is the focal length of the polarization camera. So, the azimuth angle φ P ′ and the zenith angle θ P ′ of each where V gP ′ ��� �→ (1, 1), V gP ′ ��� �→ (2, 1), and V gP ′ ��� �→ (3, 1) are the first, second, and third components of V gP ′ ��� �→ . Noting that φ P ′ needs to make quadrant judgments, then the scattering angle c P ′ of P ′ can be obtained by Substituting equation (21) into equation (9), the DOP of each pixel can be obtained, and a hypothetical DOP image is shown in Figure 6(a). And substituting equations (19)- (21) into equation (15) e shooting direction of the polarization camera is along z b , and the AOP reference direction is along y b . So, we have where E bP ′ ���→ (1, 1) and E bP ′ ���→ (2, 1) are the first and second components of E bP ′ ���→ . Figure 6(b) shows a hypothetical AOP image. When we determine orientation, the local meridian is the reference direction of AOP. e AOPLM is defined as equation (24), whose reference direction is the local meridian.
where ξ is the angular distance between axis y b and the local meridian. When the polarization camera points toward at the zenith, Figure 6(c) shows a hypothetical AOPLM image captured by the hypothetical imager.

Simulation and Experiment
Numerical simulation and experiment are carried out to verify the reliability and robustness of the proposed MOI algorithm. Moreover, the results are compared with other two classical methods. e first method is to measure the polarization at the zenith, which is named as PATZ [23][24][25][26][27][28]. e second method is an angle approach based upon Hough transform (HT) [31,32].

Simulation.
e procedure of simulation is shown in Figure 7. As described in Section 3, a hypothetical polarization camera is constructed based on the Rayleigh sky model. en, polarization images are captured by the hypothetical camera. After that, the orientation is determined by skylight compass methods. In addition, we investigate orientation estimation error at different solar elevation angles and two thousand random trials have been performed in each case.
We are particularly interested in three issues: (1) Effectiveness of the orientation determination method based on MOI. (2) Sensitivity to measurement noise.
(3) Sensitivity to numerical value discontinuity among adjacent pixels. Table 1 shows the parameters used in this simulation.

Algorithm Validation.
To test the feasibility of the proposed MOI algorithm, the polarization camera noise is set to 0. As shown in Figure 8, the bar chart compares the maximum errors, variances, and average errors of three different orientation determination methods under perfect data at different solar elevation angles. e bar charts (Figure 8) illustrates that the maximum error, variance, and average error of the PATZ method are always less than 10 −4°, 10 −14 , and 10 −5°u nder perfect data, respectively, which are the smallest compared with other two methods. Next, the HTmethod is compared with the designed MOI method. e maximum error and average error of the MOI method are always less than 0.01°and always smaller than those of HT at each solar elevation angle. In addition, the variance of the MOI method is smaller than that of the HT method except for a larger variance at h S � 50°.
Overall, the maximum errors of these three methods are always less than 0.1°, the average errors are always less than 0.1°, and the variances are always less than 10 − 6 in ideal situation. So, these three methods meet the requirement of navigation where the polarization camera noise is set to zero. In addition, the maximum error and average error of the PATZ method are the smallest, the second is MOI, and the third is HT.

Sensitivity to Measurement Noise.
To investigate the sensitivity of methods to measurement noise, Gaussian noise is added to the polarization image. e mean of Gaussian noise is set to zero, and the variance of Gaussian white noise is set to 1%(Max(AOPLM) − Min(AOPLM)). Figure 9 shows the results of sensitivity to measurement noise at h S � 0°, 10°, 20°,30°, 40°, and 50°.
For sensitivity to measurement noise, the maximum error of the PATZ method is always greater than 1°, the average error of the PATZ method is always greater than 0.1°, and the variance of the PATZ method is always greater than 10 − 6 . e results demonstrate that the POTZ method is very sensitive to Gaussian white noise. Compared to no noise situation, the error of the HT method has increased to some extent. e maximum error and average error of the HT are always greater than 0.1°and 0.01°, respectively, but less than those of PATZ. What can be clearly seen from Figure 9 is that the maximum error of the MOI method is always less than 0.011°, the variance of the MOI method is always less than 10 −7 , and the average error of the MOI method is always less than 0.01°at different solar elevation angles.
Overall, the most striking result to emerge from Figure 9 is that the stability of the MOI method is the highest for Gaussian white noise, the second is HT method and the third is PATZ.

Sensitivity to Numerical Value Discontinuity.
To investigate the sensitivity of orientation determination methods to numerical value discontinuity among adjacent pixels, not only Gaussian noise but also uniformly distributed noise is added to the polarization image. e characteristics of Gaussian noise are the same as those in Section 4.1.2. e uniformly distributed noise is a uniform distribution from 0 to 2% (Max(AOPLM) − Min(AOPLM)). Moreover, in order to simulate the discontinuity of pixel values due to crosstalk, as shown in Figure 10(a), adjacent pixels have different positive and negative characteristics of uniform distribution noise. And as shown in Figures 10(b) and 10(c), the values among adjacent pixels are obviously discontinuous.
As shown in Figure 11, the maximum error and average error of the PATZ method are always great compared with the other two methods. And Figure 11 also reveals that the maximum error of the HT method rises slightly with solar elevation, whereas that of the MOI method remains steady. Compared to only Gaussian noise situation, the error of the HT method is increased very much and even greater than 0.1°at h S � 50°. However, the maximum error of the MOI method is always less than 0.011°, variance of the MOI method is always less than 10 −7 , and the average error of the MOI method is always less than 0.01°which are consistent with only Gaussian noise situation. Overall, the most striking result to emerge from Figure 11 is that the stability   of the MOI method is the highest for numerical value discontinuity among adjacent pixels.
In summary, compared with PATZ and HT methods, the proposed MOI method can suppress the interference of measurement noise and numerical discontinuity. is result may be explained by the fact that the PATZ method only uses the polarization information at the zenith and the HT method only uses polarization information on SM-ASM. However, the MOI makes full use of the symmetry characteristics of the whole sky to determine orientation, so the MOI method is able to reduce the influence of the numerical discontinuity among adjacent pixels and measurement noise and has the highest reliability and robustness.

Experiment.
To verify the MOI method against realworld imagery, a field experiment was carried out. As shown in Figure 12, a polarized skylight navigation platform was constructed to capture the skylight polarization pattern. In this navigation platform, a pixelated polarization navigation camera was mounted on a tripod and pointed toward the zenith, which can capture a four-directional polarization image in one shot. And the parameters of the polarization camera are the same as in Table 1. In addition, a dual-antenna Global Positioning System (GPS) was used to determine the north as a benchmark. A data processing computer controlled the acquisition of the polarization image, compensated Sun motion and determined orientation. Field experiments were performed at the Nanjing University of Science and Technology (32°01 ′ 36 ′ , N; 118°51 ′ 13 ′ , E) on January 18, 2019, from 10 : 41 am to 5 : 43pm. e meteorological conditions were stable, with a clear sky and little wind. During over seven hours experiment, the solar elevation angle rose from 32.8°to a maximum of 37.3°a nd finally dropped to -3.9°over time. e solar azimuth angle changed from 135.9°to 249.4°, in which 482 image sets were captured using the polarized skylight navigation platform. e procedure of the experiment is shown in Figure 13. Firstly, polarization images are captured by the polarized skylight navigation platform. en, each polarization image is splited into four intensity images. After that, the DOP image and AOP image are calculated. Moreover, the AOP image is transformed into the AOPLM image by equation (24). Finally, the skylight compass method is designed and applied to determine orientation.
As shown in Figure 14, in actual, polarization images captured by pixelated polarization camera have significant measurement noise and numerical discontinuity among adjacent pixels. e experimental results are shown in Figure 15 and Table 2. For polarization images, the maximum error of the MOI is 0.44°, the variance of the MOI is 0.040, and the average error of MOI is 0.242. Compared with the other two methods, this method has the highest accuracy and best stability. e results of the experiment confirm that the MOI method is able to determine orientation against real-world polarization imagery captured by the pixelated polarization camera and reduce the influence of the numerical discontinuity among adjacent pixels and measurement noise.

Conclusions
e main result of this paper is to propose a new skylight compass algorithm based on MOI, which makes full use of the symmetry characteristics of the whole sky to determine orientation. So, the MOI method can effectively suppress the interference of numerical discontinuity among pixels and measurement noise and has a higher level of reliability, repeatability, and robustness. is solution is well suited to practical application, such as initial alignment for aircrafts, ships, vehicles, and other moving carriers.
Simulation has validated the effectiveness of the proposed MOI algorithm. en, the influence of measurement noise and numerical discontinuity among adjacent pixels was investigated, and the results showed that this algorithm can effectively acquire the orientation on the condition of small measurement noise and numerical discontinuity. In addition, this method was compared with other two classical methods. e compared result highlights the stability of the MOI method. Finally, a field experiment based on the polarized navigation platform was also undertaken and the compass error of the MOI method is less than 0.44°for the actual polarization image under clear sky.
e experimental results demonstrate that the MOI method is able to determinate orientation for the polarization image captured by the pixelated polarization camera in practice.
Although the impact of measurement noise on orientation determination has been studied in detail in this paper, polarized skylight navigation may not only be affected by measurement noises but also by weather conditions. So, the impact of different weather conditions on orientation determination would be the focus of our future research.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.