# Polarization Navigation Simulation System and Skylight Compass Method Design Based upon Moment of Inertia

Authors: Huaju Liang; Hongyang Bai; Ning Liu; Xiubao Sui
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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. Then, orientation determination is transformed into solving the extreme value of MOI. This 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.


---

## Body

## 1. 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–4]. Locusts exploit the polarization pattern to extract directional information from the sky [5–7]. Dung beetles can sense polarized light to roll dung piles away along a straight path [8–10]. Drosophilae detect the electric field vector (E-vector) of polarized skylight as a compass cue [11–13].The 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 (Figure1(a)) [6, 11, 12]. In the DRA, ommatidial microvilli are arranged in different angles (Figure 1(b)) [7, 14]. This physiological structure enables insects to detect polarized light. Inspired by this physiological structure, a pixelated polarized camera is designed [15–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.Figure 1
The basis of polarized skylight navigation. (a) Polarization-sensitive dorsal rim area (DRA); (b) DRA ommatidium, the ommatidial microvilli are arranged in different angles. The microvilli of photoreceptor 1 are perpendicular to those of 2, 3, 4, 7, and 8. And microvilli of photoreceptors 5 and 6 are oriented irregularly; (c) the Rayleigh sky model, where the yellow point represents the Sun, the red line is the SM-ASM, double arrows are E-vectors, and the blue lines indicate the degree of polarization (DOP). The thicker the blue line, the greater the DOP.
(a)
(b)
(c)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]. The 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. The 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. The 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–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–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. Then, 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.


## 2. Orientation Determination Algorithm Design

As shown in Figure2, 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–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. Thus, the MOI of the polarization image can be calculated to extract the symmetry axes.Figure 2
Three rigid bodies. The three red lines are their symmetry axis, respectively, which also make the moment of inertia reach the extreme values.
(a)
(b)
(c)To simplify the description of the proposed method, a body coordinate systemobxbybzb of the polarization camera is established. As shown in Figure 3, ob is the center of the polarization image, xb and yb are aligned with the direction of the column and row of the polarization image, and zb points to the shoot direction of the polarization camera.Figure 3
Body coordinate system of the polarization camera. (a)obxbybzb is the body coordinate system of the polarization camera. The ob is the center of the polarization image, xb and yb are aligned with the direction of the column and row of the polarization image, and zb points toward the shoot direction of the polarization camera. (b) obxbyb plane of obxbybzb on the polarization image, and mx and my represent the count of column and row pixels, respectively.
(a)
(b)For a rigid body, each mass unit isdm, and the distance to a line is r. The definition of MOI is(1)I=∑r2dm.In a grayscale image, the rigid body is composed of a set of pixel pointsP′=xi,yi,i∈1,N, where N is the number of pixels in the image. For polarized skylight navigation, P′ is a circular region with ob as its origin. Each pixel point can be considered to be equivalent to each mass unit in the rigid body. The MOI of a rigid body relative to a line can be obtained by integrating the MOI of each mass unit.To ensure that all lines are considered, the equation of the axis of symmetry is defined as(2)ay+bx+c=0a≠0∨b≠0,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. Therefore, the axis of symmetry would definitely pass through the center of the polarization image, so c=0. Then, the whole MOI of a polarization image is(3)I=∑i=1NμiDi2,where Di=ayi+bxi/a2+b2 is the distance between the pixel point xi,yi 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.The MOI of a rigid body reaches the extreme values when the rigid body rotates on its symmetry axes. Therefore, the partial derivatives ofI should be calculated to extract symmetry axes.The derivative ofI with respect to a is(4)∂I∂a=2ba2+b2b2−a2∑i=1Nμixiyi+ab∑i=1Nμiyi2−xi2.The derivative ofI with respect to b is(5)∂I∂b=2aa2+b2a2−b2∑i=1Nμixiyi−ab∑i=1Nμiyi2−xi2.LetW=∑i=1Nμiyi2−xi2 and G=∑i=1Nuixiyi, when the MOI reaches the extreme values, we have ∂I/∂a=0 and ∂I/∂b=0. Thus,(6)bb2−a2G+abW=0,aa2−b2G−abW=0.For equation (6), when a=0 or b=0, G=0. When a≠0 and b≠0, a/b=W±W2+4G2/2G. So, the axes of symmetry can be determined by(7)a=0∨b=0G=0,ab=W±W2+4G22GG≠0.In short, two straight lines are obtained. And we have(8)W+W2+4G22G⋅W−W2+4G22G=−1.So, the two straight lines are perpendicular to each other. Only one of them is the SM-ASM. Then, 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.


## 3. 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.


### 3.1. Rayleigh Sky Model

As shown in Figure4, oxgygzg is the geographic coordinate system. The east is labeled xg, the north is labeled yg, and the up is labeled zg.Figure 4
The geographic coordinate system.oxgygzg is the geographic coordinate system. The 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. γP is the scattering angle of P. θS is the solar zenith angle. hS 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.The Rayleigh sky model predicts the degree of polarization (DOP) [24] as(9)DOP=DOPmaxsin2γP1+cos2γP,where DOPmax denotes the maximum DOP of the whole sky and equals 1 for an ideal (single scattering) sky, γP is the angular distance between the Sun and the observation point P. The views of DOP are shown in Figures 5(a) and 5(b).Figure 5
Hypothetical skylight polarization patterns. (a) Three-dimensional (3D) view of DOP; (b) zenith view of DOP; (c) 3D view of AOP; (d) zenith view of AOP; (e) 3D view of the E-vector; (f) zenith view of the E-vector.
(a)
(b)
(c)
(d)
(e)
(f)The Rayleigh sky model predicts the AOP [31] as(10)cosAOP=sinφS−φPsinγPsinθS,(11)sinAOP=sinθPcosθS−cosφS−φPcosθPsinθSsinγP,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. The 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.The E-vector at the pointP is(12)EgP⟶=VgP⟶cosAOP+HgP⟶sinAOP,where(13)VgP⟶=−cosθPsinφP,−cosθPcosφP,sinθPT,(14)HgP⟶=−cosφP,sinφP,0T.In equations (13) and (14), VgP⟶ and HgP⟶ are tangential directions of the meridian and weft at the point P, respectively. The superscript T represents matrix transpose.Substituting equations (10) and (11) into equation (12), components of EgP⟶ can be given by(15)EgP⟶=cosφScosθPsinθS−cosφPsinθPcosθS/sinγP−sinφScosθPsinθS+sinφPsinθPcosθS/sinγPsinφS−φPsinθSsinθP/sinγP.The views of the E-vector are shown in Figures5(e) and 5(f).


### 3.2. 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 pixelP′ is(16)VbP′⟶=Dx×mx−ηx+12,Dy×−my+ηy+12,0T,where Dx and Dy 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 mx and my 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(17)Cbg=cosβcosφ+sinβsinαsinφcosαsinφsinβcosφ−cosβsinαsinφ−cosβsinφ+sinβsinαcosφcosαcosφ−sinβsinφ−cosβsinαcosφ−sinβcosαsinαcosβcosα,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 alongzb⟶=0,0,1T. Thus, in the geographic coordinate system, the shooting direction of the pixel P′ is(18)VgP′⟶=CbgVbP′⟶+f0,0,1T,where f is the focal length of the polarization camera.So, the azimuth angleφP′ and the zenith angle θP′ of each pixel P′ are(19)φP′=arctanVgP′⟶1,1VgP′⟶2,1,(20)θP′=arccosVgP′⟶3,1VgP′⟶,where VgP′⟶1,1, VgP′⟶2,1, and VgP′⟶3,1 are the first, second, and third components of VgP′⟶. Noting that φP′ needs to make quadrant judgments, then the scattering angle γP′ of P′ can be obtained by(21)γP′=arccoscosθP′cosθS+sinθP′sinθScosφS−φP′.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), the polarization E-vector EgP′⟶ of the pixel P′ in the geographic coordinate system can be obtained. So, the polarization E-vector EbP′⟶ of P′ in the body coordinate system is(22)EbP′⟶=CgbEgP′⟶.Figure 6
Polarization images captured by the hypothetical imager. (a) Hypothetical DOP image; (b) hypothetical AOP image; (c) hypothetical AOPLM image.
(a)
(b)
(c)The shooting direction of the polarization camera is alongzb, and the AOP reference direction is along yb. So, we have(23)AOP=arctanEbP′⟶1,1EbP′⟶2,1,where EbP′⟶1,1 and EbP′⟶2,1 are the first and second components of EbP′⟶. Figure 6(b) shows a hypothetical AOP image. When we determine orientation, the local meridian is the reference direction of AOP. The AOPLM is defined as equation (24), whose reference direction is the local meridian.(24)AOPLM=AOP−ξ,where ξ is the angular distance between axis yb and the local meridian. When the polarization camera points toward at the zenith,(25)ξ=arctan−my+ηy+1/2mx−ηx+1/2.Figure6(c) shows a hypothetical AOPLM image captured by the hypothetical imager.


## 3.1. Rayleigh Sky Model

As shown in Figure4, oxgygzg is the geographic coordinate system. The east is labeled xg, the north is labeled yg, and the up is labeled zg.Figure 4
The geographic coordinate system.oxgygzg is the geographic coordinate system. The 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. γP is the scattering angle of P. θS is the solar zenith angle. hS 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.The Rayleigh sky model predicts the degree of polarization (DOP) [24] as(9)DOP=DOPmaxsin2γP1+cos2γP,where DOPmax denotes the maximum DOP of the whole sky and equals 1 for an ideal (single scattering) sky, γP is the angular distance between the Sun and the observation point P. The views of DOP are shown in Figures 5(a) and 5(b).Figure 5
Hypothetical skylight polarization patterns. (a) Three-dimensional (3D) view of DOP; (b) zenith view of DOP; (c) 3D view of AOP; (d) zenith view of AOP; (e) 3D view of the E-vector; (f) zenith view of the E-vector.
(a)
(b)
(c)
(d)
(e)
(f)The Rayleigh sky model predicts the AOP [31] as(10)cosAOP=sinφS−φPsinγPsinθS,(11)sinAOP=sinθPcosθS−cosφS−φPcosθPsinθSsinγP,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. The 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.The E-vector at the pointP is(12)EgP⟶=VgP⟶cosAOP+HgP⟶sinAOP,where(13)VgP⟶=−cosθPsinφP,−cosθPcosφP,sinθPT,(14)HgP⟶=−cosφP,sinφP,0T.In equations (13) and (14), VgP⟶ and HgP⟶ are tangential directions of the meridian and weft at the point P, respectively. The superscript T represents matrix transpose.Substituting equations (10) and (11) into equation (12), components of EgP⟶ can be given by(15)EgP⟶=cosφScosθPsinθS−cosφPsinθPcosθS/sinγP−sinφScosθPsinθS+sinφPsinθPcosθS/sinγPsinφS−φPsinθSsinθP/sinγP.The views of the E-vector are shown in Figures5(e) and 5(f).


## 3.2. 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 pixelP′ is(16)VbP′⟶=Dx×mx−ηx+12,Dy×−my+ηy+12,0T,where Dx and Dy 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 mx and my 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(17)Cbg=cosβcosφ+sinβsinαsinφcosαsinφsinβcosφ−cosβsinαsinφ−cosβsinφ+sinβsinαcosφcosαcosφ−sinβsinφ−cosβsinαcosφ−sinβcosαsinαcosβcosα,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 alongzb⟶=0,0,1T. Thus, in the geographic coordinate system, the shooting direction of the pixel P′ is(18)VgP′⟶=CbgVbP′⟶+f0,0,1T,where f is the focal length of the polarization camera.So, the azimuth angleφP′ and the zenith angle θP′ of each pixel P′ are(19)φP′=arctanVgP′⟶1,1VgP′⟶2,1,(20)θP′=arccosVgP′⟶3,1VgP′⟶,where VgP′⟶1,1, VgP′⟶2,1, and VgP′⟶3,1 are the first, second, and third components of VgP′⟶. Noting that φP′ needs to make quadrant judgments, then the scattering angle γP′ of P′ can be obtained by(21)γP′=arccoscosθP′cosθS+sinθP′sinθScosφS−φP′.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), the polarization E-vector EgP′⟶ of the pixel P′ in the geographic coordinate system can be obtained. So, the polarization E-vector EbP′⟶ of P′ in the body coordinate system is(22)EbP′⟶=CgbEgP′⟶.Figure 6
Polarization images captured by the hypothetical imager. (a) Hypothetical DOP image; (b) hypothetical AOP image; (c) hypothetical AOPLM image.
(a)
(b)
(c)The shooting direction of the polarization camera is alongzb, and the AOP reference direction is along yb. So, we have(23)AOP=arctanEbP′⟶1,1EbP′⟶2,1,where EbP′⟶1,1 and EbP′⟶2,1 are the first and second components of EbP′⟶. Figure 6(b) shows a hypothetical AOP image. When we determine orientation, the local meridian is the reference direction of AOP. The AOPLM is defined as equation (24), whose reference direction is the local meridian.(24)AOPLM=AOP−ξ,where ξ is the angular distance between axis yb and the local meridian. When the polarization camera points toward at the zenith,(25)ξ=arctan−my+ηy+1/2mx−ηx+1/2.Figure6(c) shows a hypothetical AOPLM image captured by the hypothetical imager.


## 4. 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. The first method is to measure the polarization at the zenith, which is named as PATZ [23–28]. The second method is an angle approach based upon Hough transform (HT) [31, 32].


### 4.1. Simulation

The procedure of simulation is shown in Figure7. As described in Section 3, a hypothetical polarization camera is constructed based on the Rayleigh sky model. Then, 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.Figure 7
Procedure of the simulation. (a) DOP of the Rayleigh sky model; (b) AOP of the Rayleigh sky model; (c) hypothetical DOP image; (d) hypothetical AOPLM image; (e) the result of orientation determination shown on the hypothetical DOP image; (f) the result of orientation determination shown on the hypothetical AOPLM image. The yellow point represents the Sun and the white dotted line represents the result of orientation determination.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.Table1 shows the parameters used in this simulation.Table 1
Simulation parameters.
Symbol
Value
Units
Description
D
x
3.45
Μm
Pixel column size of CMOS or CCD
D
y
3.45
Μm
Pixel row size of CMOS or CCD
η
x
2048
Pixel
Number of pixels in column direction
η
y
2448
Pixel
Number of pixels in row direction
f
8
Mm
Focus length of polarization camera
D
O
P
max
1
—
Maximum DOP of sky


#### 4.1.1. Algorithm Validation

To test the feasibility of the proposed MOI algorithm, the polarization camera noise is set to 0. As shown in Figure8, 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.Figure 8
Statistical analysis of the results for algorithm validation. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)The bar charts (Figure8) illustrates that the maximum error, variance, and average error of the PATZ method are always less than 10−4°, 10−14, and 10−5° under perfect data, respectively, which are the smallest compared with other two methods. Next, the HT method is compared with the designed MOI method. The 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 hS=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 than10−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.


#### 4.1.2. Sensitivity to Measurement Noise

To investigate the sensitivity of methods to measurement noise, Gaussian noise is added to the polarization image. The mean of Gaussian noise is set to zero, and the variance of Gaussian white noise is set to 1%(Max(AOPLM) − Min(AOPLM)). Figure9 shows the results of sensitivity to measurement noise at hS=0°, 10°, 20°,30°, 40°, and 50°.Figure 9
Statistical analysis of the results for sensitivity to measurement error. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)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 than10−6. The 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. The 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 Figure9 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.


#### 4.1.3. 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. The characteristics of Gaussian noise are the same as those in Section4.1.2. The 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.Figure 10
Numerical discontinuity among adjacent pixels. (a) Positive and negative characteristics of uniform distribution noise; (b) the hypothetical AOPLM image with noise; (c) local magnification of the hypothetical AOPLM image.As shown in Figure11, 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 hS=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.Figure 11
Statistical analysis of the results for sensitivity to numerical value discontinuity. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)In summary, compared with PATZ and HT methods, the proposed MOI method can suppress the interference of measurement noise and numerical discontinuity. This 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.


### 4.2. Experiment

To verify the MOI method against real-world imagery, a field experiment was carried out. As shown in Figure12, 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.Figure 12
Polarized skylight navigation platform.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. The 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° and finally dropped to -3.9° over time. The solar azimuth angle changed from 135.9° to 249.4°, in which 482 image sets were captured using the polarized skylight navigation platform.The procedure of the experiment is shown in Figure13. Firstly, polarization images are captured by the polarized skylight navigation platform. Then, 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.Figure 13
Procedure of the experiment. (a) Polarization imagery with four different directional polarizing filters; (b) four different directional light intensity images; (c) DOP image; (d) AOP image; (e) AOPLM image; (f) the result of orientation determination.As shown in Figure14, in actual, polarization images captured by pixelated polarization camera have significant measurement noise and numerical discontinuity among adjacent pixels. The 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. The 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.Figure 14
Measurement noise and numerical discontinuity among adjacent pixels of the actual polarization image. (a) DOP image; (b) AOP image; (c) AOPLM image; (d) local magnification of the DOP image; (e) local magnification of the AOP image; (f) local magnification of the AOPLM image.Figure 15
The result of experiment for orientation determination.Table 2
The statistical result of experiment.
Method
Max error (deg)
Variance
Average error (deg)
PATZ
−1.11
0.115
0.263
HT
−0.55
0.060
0.266
MOI
−0.44
0.040
0.242


## 4.1. Simulation

The procedure of simulation is shown in Figure7. As described in Section 3, a hypothetical polarization camera is constructed based on the Rayleigh sky model. Then, 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.Figure 7
Procedure of the simulation. (a) DOP of the Rayleigh sky model; (b) AOP of the Rayleigh sky model; (c) hypothetical DOP image; (d) hypothetical AOPLM image; (e) the result of orientation determination shown on the hypothetical DOP image; (f) the result of orientation determination shown on the hypothetical AOPLM image. The yellow point represents the Sun and the white dotted line represents the result of orientation determination.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.Table1 shows the parameters used in this simulation.Table 1
Simulation parameters.
Symbol
Value
Units
Description
D
x
3.45
Μm
Pixel column size of CMOS or CCD
D
y
3.45
Μm
Pixel row size of CMOS or CCD
η
x
2048
Pixel
Number of pixels in column direction
η
y
2448
Pixel
Number of pixels in row direction
f
8
Mm
Focus length of polarization camera
D
O
P
max
1
—
Maximum DOP of sky


### 4.1.1. Algorithm Validation

To test the feasibility of the proposed MOI algorithm, the polarization camera noise is set to 0. As shown in Figure8, 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.Figure 8
Statistical analysis of the results for algorithm validation. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)The bar charts (Figure8) illustrates that the maximum error, variance, and average error of the PATZ method are always less than 10−4°, 10−14, and 10−5° under perfect data, respectively, which are the smallest compared with other two methods. Next, the HT method is compared with the designed MOI method. The 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 hS=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 than10−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.


### 4.1.2. Sensitivity to Measurement Noise

To investigate the sensitivity of methods to measurement noise, Gaussian noise is added to the polarization image. The mean of Gaussian noise is set to zero, and the variance of Gaussian white noise is set to 1%(Max(AOPLM) − Min(AOPLM)). Figure9 shows the results of sensitivity to measurement noise at hS=0°, 10°, 20°,30°, 40°, and 50°.Figure 9
Statistical analysis of the results for sensitivity to measurement error. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)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 than10−6. The 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. The 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 Figure9 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.


### 4.1.3. 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. The characteristics of Gaussian noise are the same as those in Section4.1.2. The 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.Figure 10
Numerical discontinuity among adjacent pixels. (a) Positive and negative characteristics of uniform distribution noise; (b) the hypothetical AOPLM image with noise; (c) local magnification of the hypothetical AOPLM image.As shown in Figure11, 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 hS=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.Figure 11
Statistical analysis of the results for sensitivity to numerical value discontinuity. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)In summary, compared with PATZ and HT methods, the proposed MOI method can suppress the interference of measurement noise and numerical discontinuity. This 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.


## 4.1.1. Algorithm Validation

To test the feasibility of the proposed MOI algorithm, the polarization camera noise is set to 0. As shown in Figure8, 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.Figure 8
Statistical analysis of the results for algorithm validation. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)The bar charts (Figure8) illustrates that the maximum error, variance, and average error of the PATZ method are always less than 10−4°, 10−14, and 10−5° under perfect data, respectively, which are the smallest compared with other two methods. Next, the HT method is compared with the designed MOI method. The 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 hS=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 than10−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.


## 4.1.2. Sensitivity to Measurement Noise

To investigate the sensitivity of methods to measurement noise, Gaussian noise is added to the polarization image. The mean of Gaussian noise is set to zero, and the variance of Gaussian white noise is set to 1%(Max(AOPLM) − Min(AOPLM)). Figure9 shows the results of sensitivity to measurement noise at hS=0°, 10°, 20°,30°, 40°, and 50°.Figure 9
Statistical analysis of the results for sensitivity to measurement error. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)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 than10−6. The 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. The 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 Figure9 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.


## 4.1.3. 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. The characteristics of Gaussian noise are the same as those in Section4.1.2. The 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.Figure 10
Numerical discontinuity among adjacent pixels. (a) Positive and negative characteristics of uniform distribution noise; (b) the hypothetical AOPLM image with noise; (c) local magnification of the hypothetical AOPLM image.As shown in Figure11, 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 hS=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.Figure 11
Statistical analysis of the results for sensitivity to numerical value discontinuity. (a) Comparison of maximum errors; (b) comparison of variances; (c) comparison of average errors.
(a)
(b)
(c)In summary, compared with PATZ and HT methods, the proposed MOI method can suppress the interference of measurement noise and numerical discontinuity. This 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.


## 4.2. Experiment

To verify the MOI method against real-world imagery, a field experiment was carried out. As shown in Figure12, 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.Figure 12
Polarized skylight navigation platform.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. The 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° and finally dropped to -3.9° over time. The solar azimuth angle changed from 135.9° to 249.4°, in which 482 image sets were captured using the polarized skylight navigation platform.The procedure of the experiment is shown in Figure13. Firstly, polarization images are captured by the polarized skylight navigation platform. Then, 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.Figure 13
Procedure of the experiment. (a) Polarization imagery with four different directional polarizing filters; (b) four different directional light intensity images; (c) DOP image; (d) AOP image; (e) AOPLM image; (f) the result of orientation determination.As shown in Figure14, in actual, polarization images captured by pixelated polarization camera have significant measurement noise and numerical discontinuity among adjacent pixels. The 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. The 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.Figure 14
Measurement noise and numerical discontinuity among adjacent pixels of the actual polarization image. (a) DOP image; (b) AOP image; (c) AOPLM image; (d) local magnification of the DOP image; (e) local magnification of the AOP image; (f) local magnification of the AOPLM image.Figure 15
The result of experiment for orientation determination.Table 2
The statistical result of experiment.
Method
Max error (deg)
Variance
Average error (deg)
PATZ
−1.11
0.115
0.263
HT
−0.55
0.060
0.266
MOI
−0.44
0.040
0.242


## 5. Conclusions

The 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. This 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. Then, 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. The 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. The 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.



---

*Source: 4081269-2020-05-13.xml*