A Novel Skylight Orientation Sensor for Autonomous Navigation

—The angle of the polarization and the degree of polarization (DOP) of the scattered skylight are symmetrically distributed concerning the solar meridian. Based on the symmetry of the skylight polarization distribution pattern, this paper proposes a novel skylight orientation sensor consists of a camera, an S-waveplate, and a linear polarizer. The skylight orientation sensor is using the image polarization encoding capability of the S-waveplate and the linear polarizer to convert the skylight polarization information into the image’s symmetry axis extraction, which has the advantages of no resolution loss and instantaneous field of view error. The symmetry axis in the image is consistent with the solar meridian. Therefore, the angle between the solar meridian and the skylight orientation sensor reference axis can be obtained without calculating the polarization information, which is also beneficial for real-time performance. The angle measurement accuracy and uncertainty of the skylight orientation sensor are verified by numerical simulation and outdoor experiments. The results demonstrate that the skylight orientation sensor has good application potential in autonomous navigation.


I. Introduction
UNLIGHT enters the atmosphere and is scattered by atmospheric particles to form partially polarized light. The angle of the polarization (AOP) and degree of polarization (DOP) of the scattered partially polarized light, known as the skylight polarization distribution pattern, are symmetrically distributed concerning the solar meridian [1]. The skylight polarization distribution pattern is also related to the position of the sun. Therefore, it can be used for autonomous navigation. Many creatures, such as ants, crickets, beetles, etc, have been found to use the skylight polarization distribution pattern for navigation [2]- [9]. As a new navigation technology type, bionic polarized light navigation has increasingly attracted researchers' attention, given its unique advantages, such as antielectromagnetic interference and no error accumulation. Up to now, researchers have developed different kinds of polarized light navigation sensors, mainly classified into "point-source" polarization navigation sensors and imaging polarization navigation sensors. The "point-source" polarization navigation sensor [10]- [21] has the advantages of a simple structure, low costs, high precision, and good real-time performance. The simple structure is also conducive to its integration [15], [18]. However, the "point-source" polarization navigation sensor is easily affected by inclement weather, such as clouds, fog, and haze. Researchers developed the imaging polarization navigation sensor [22]- [34] to improve the anti-jamming capability of the sensor.
The imaging polarization navigation sensor mainly includes multi-camera [22], [24], [25], [27], [28], division-of-time [30]- [32], and division-of-focal-plane [26], [33]- [35] polarization imaging sensors. Compared with multi-camera and division-oftime polarization navigation sensors, the division-of-focalplane sensor has the advantages of integration and real-time. However, the division-of-focal-plane sensor regards the adjacent pixels as a superpixel, which leads to resolution loss and instantaneous field of view error. To avoid the abovementioned shortcomings, this paper proposes a novel skylight orientation sensor consists of a camera, an S-waveplate, and a linear polarizer. Similar to the application of magnetic compass in the earth's magnetic field: sensing the orientation of magnetic field used for navigation, the skylight orientation sensor is used in skylight field: using the image polarization encoding capability of the S-waveplate and the linear polarizer to convert the skylight polarization information into the image's symmetry axis extraction. The symmetry axis in the image is consistent with the solar meridian. Therefore, the angle between the solar meridian and the skylight orientation sensor reference axis can be obtained without calculating the polarization information which is also beneficial to real-time performance. We can further obtain the azimuth of the skylight orientation sensor used for autonomous navigation by inquiring the astronomical almanac to get the azimuth of the solar meridian.
In the following paper, we first analyze the skylight polarization pattern based on the Rayleigh scattering theory and S mathematically deduce the lightwave modulation principle of the proposed skylight orientation sensor. Then, a numerical simulation is used to analyze the angle measurement accuracy and uncertainty of the skylight orientation sensor in the case of considering ambient light noise, and the time cost is also analyzed. Finally, we set up an experimental system and carry out outdoor experiments to further verify the feasibility of the skylight orientation sensor. The numerical simulation and outdoor experiments results indicate that the skylight orientation sensor has good application potential in autonomous navigation.

A. Skylight Polarization Pattern
Sunlight is scattered by atmospheric particles to form a partially polarized skylight. In clear weather conditions, the atmospheric particles' size is much smaller than the wavelength of sunlight, and Rayleigh scattering of sunlight occurs [36]. The skylight polarization distribution pattern based on Rayleigh scattering is described by establishing a horizontal coordinate system. As shown in Fig. 1, the O, P, S, and Z points represent the observer, observed point, sun, and zenith point, respectively. A p and h p represent the azimuth and altitude of the observed point respectively. As and hs represent the azimuth and altitude of the sun respectively. is the scattering angle (the angle between the direction of incident light and the direction of observation), and is the polarization azimuth of the scattered light at the point (the angle between the E vector vibration direction of point and the meridian passing through point). (1) tan . ( It can be seen from the calculation formulas that when the scattering angle is 90°, the degree of polarization of the scattered light is 1, forming fully polarized light. Further, when the scattering angle is other angles, a partially polarized light is formed. In Fig. 1, the direction and thickness of the blue dotted line represent the electric vector's direction of the polarized light and the degree of polarization, respectively. Fig. 2 shows the schematic diagram corresponding to the skylight orientation sensor. The skylight orientation sensor comprises a fisheye lens, a filter, an S-waveplate, a linear polarizer, a relay lens, and a mono camera. This paper used a combination of filter, S-waveplate, and linear polarizer, called lightwave modulator. The lightwave modulator cannot be placed directly on the photosensitive plane of the camera. Therefore, a relay lens was used here. The lightwave modulator is placed in the object plane of the relay lens. The relay lens's object plane coincides with the image plane of the fisheye lens, while the image plane of the relay lens coincides with the photosensitive plane of the camera.

B. Skylight orientation sensor Lightwave Modulation Principle
S-waveplate is a kind of fused quartz glass irradiated by femtosecond laser beam, which forms a superstructure wave plate with birefringence index periodic nanostructure inside. It can convert the linearly polarized beam into the radial polarized beam or angular polarized beam [37], [38]. The structural schematic diagram of the S-waveplate is given in Fig. 3. The Swaveplate is equivalent to half-wave plates with different orientations in different radial orientations. The geometric center O of the S-waveplate is the origin. The X-axis's positive direction is the zero-degree reference direction, and the angle in the clockwise direction is positive. In the radial direction of the S-waveplate with the azimuth angle , there is a half-wave with azimuth 2 ⁄ . The Mueller matrix of the half-wave with azimuth angle is written as Then we obtain the Mueller matrix of S-waveplate: The Mueller matrix of the linear polarizer with an azimuth is defined as The normalized Stokes vector of the polarized beam with azimuth , degree of polarization P, and ellipticity can be given by After modulated by S-waveplate and linear polarizer, the output beam can be expressed as follows: (7) where is the Stokes vector of the normally incident polarized beam, is the Mueller matrix of S-waveplate, is the Mueller matrix of a linear polarizer, and is the Stokes vector of the modulated beam. In this paper, the transmission axis of the linear polarizer is consistent with the zero-degree reference direction of the S-wave plate, i.e., 0. When the skylight orientation sensor is used in the skylight field, the skylight is first imaged in the object plane of the relay lens through the fisheye lens. It is then imaged in the photosensitive plane of the camera after the modulation of the lightwave modulator and the transmission of the relay lens shown in Fig. 2. The obtained image has a regular light intensity distribution pattern symmetric to the solar meridian. Additionally, the skylight orientation sensor is composed of independent optical devices. Therefore, the effects of sky polarization distribution under different weather conditions on skylight orientation sensor angle's measurement accuracy can be studied by replacing optical devices with different operating wavelengths.

III. SIMULATION
Based on the single Rayleigh scattering model and the light wave modulation principle of the skylight orientation sensor, we numerically simulate the skylight polarization distribution mode and the skylight orientation sensor imaging map. Fig 4(a, b) shows the top view of the simulation calculation: the angle of the polarization (AOP) and the degree of polarization (DOP) of the scattered skylight polarization distribution pattern with a solar altitude angle of 20 deg and an azimuth angle of 0 deg. The AOP and DOP are symmetrically distributed concerning the solar meridian. Besides, the scattered sunlight is mainly linearly polarized light, that is, =0. To further simplify the calculation, we assume that the intensity of scattered sunlight at any point in the sky is 1. Therefore, the scattered sunlight can be expressed by is approximate to the angle between the electric vector at point P in the projection plane of the skylight body relative to the due south direction. Substituting (8) into (7) yields the Stokes vector of the skylight after being modulated by the lightwave modulator. Fig. 4(c) shows the light intensity distribution of the modulated skylight. We added AOP random noise with a mean of 0 and a standard deviation of 0.005π and DOP random noise with a mean of 0 and a standard deviation of 0.005 to consider realistic factors, such as the influence of ambient light. Fig. 4(d) shows the modulated light intensity distribution after considering ambient light noise.
In particular, the electric vector direction at the position of the solar meridian is perpendicular to the solar meridian direction. In the camera coordinate system, the azimuth of the half-wave plate in the S-waveplate is 2 ⁄ , while the solar meridian azimuth is . The electric vector azimuth is always 90 deg at the solar meridian position after modulated by the Swave plate, which is perpendicular to the linear polarizer's transmission axis. Therefore, the gray value is the smallest in the solar meridian direction. According to the intensity image, we do not need to calculate the polarization information. We can obtain the angle between the solar meridian and the skylight orientation sensor reference axis by the following steps. First, the center point of the intensity image was determined, and the center point was passed by the solar meridian. Following this, the image was cropped symmetrically around this point to remove the edge portion. Because the edge portion cannot be well modulated due to the vignetting effect, removing these edge portions can reduce the image's size, which is conducive to improving processing speed and measurement accuracy, respectively. Second, the image was divided into two parts according to a predefined threshold, and the two dark parts were used as the area of interest to calculate its angle by the curve fitting method. After that, by inquiring the astronomical almanac to get the azimuth of the solar meridian, we can further obtain the azimuth of the celestial light compass.
We used numerical simulation to obtain modulated light intensity images with a solar altitude angle of 20 deg and azimuth angle 0,10, … ,180 deg, considering ambient light noise (AOP random noise with a mean of 0 and a standard deviation of 0.005π and DOP random noise with a mean of 0 and a standard deviation of 0.005) to verify the accuracy and uncertainty of the method mentioned above. Following this, we calculated the azimuth of the solar meridian. Fig. 5 shows the simulation results of the calculated angle and calculated angle error. The solid red line is the real angle, the blue pentagram line is the calculated angle, and the solid green line is the calculated angle error. The angle error of the simulation result where x i represents the calculated angle, y i is the actual angle, and N is the number of the measurement sample. The sample standard deviation of the angle measurement is 0.093. In addition to the accuracy and uncertainty of angle measurement, the time cost is another important indicator that should be paid attention to. The traditional skylight orientation sensors need to first obtain four light intensity image (I 0 , I 45 , I 90 , I 135 ), then calculate the polarization information, and the angle between the sun meridian and the carrier axis is finally obtained by some algorithms, such as Hough transform [29], symmetry axis scanning [1] or curve fitting [22]. Our skylight orientation sensor only needs to obtain a light intensity image. The symmetry axis in the image is consistent with the solar meridian. Therefore, the angle between the solar meridian and the skylight orientation sensor reference axis can be obtained without calculating the polarization information, which is beneficial for real-time performance . Fig 6 (a) and (b) show the flowchart of angle calculation about the traditional skylight orientation sensor and our proposed skylight orientation sensor. Because the polarization angle image obtained by numerical calculation and the light intensity image obtained by our sensor are both symmetrical about the solar meridian. Therefore, in the angle calculation step, we use the same method (curve fitting) to obtain the azimuth of the solar meridian. Fig 6 (a) and (b) take about 127 ms and 78 ms to process a 1000*1000-pixel image by MATLAB R2017b on a desktop (Intel Core i5-7400 @3 GHz, 8 GB memory), respectively. In the future practical application, we think that using our sensor can improve the response speed by about 39%.

IV. EXPERIMENT AND RESULT
To further verify the outdoor angle measurement accuracy of the skylight orientation sensor, we set up an experimental system to collect outdoor experimental data. The skylight orientation sensor was first calibrated by integrating sphere in the laboratory. The performance parameters of the main components are given in Table 1. The experimental system is shown in Fig. 7. As shown, the precision turntable was placed horizontally on the ground plane. The skylight orientation sensor was fixed on a precision turntable so that it faces the apex of the sky. The skylight orientation sensor collected an image for each 10 deg rotation of the precision turntable. A total of 37 images were collected as the precision turntable rotated for 360 deg ， which takes less than 2 minutes. After acquiring the  image, the angle between the solar meridian and the skylight orientation sensor reference axis was calculated. We repeated the outdoor experiment five times. In every experiment, the first measured angle was taken as the true value. The reference angle was the sum of the first measured angle and the turntable's relative rotation angle. Fig 8 (a) shows the results of the first outdoor experiment. The solid red line is the reference angle, and the blue pentagram line is the measured angle. In Fig  8 (b), Ai, i = 1, 2, 3, 4,5 denotes the measured angle error of the ith outdoor experiment. The measured angle error is the difference between the measured angle and the reference angle.
In the process of collecting outdoor experimental data, the selfrotation of the sky polarized light field with time and the deviation between the body axis of the skylight orientation sensor and the direction facing the zenith will affect the accuracy of the calculation results. As shown, the measured angle errors have periodicity and good consistency. Therefore, we use the average measured angle error of Ai, i = 1, 2, 3 to compensate for A4 and A5. The maximum measured angle error after compensation of A4 and A5 is less than 0.3 deg shown in Fig 7 (b). The sample standard deviations of A4 and A5 after compensation are 0.165 and 0.171, respectively.

V. CONCLUSION
In conclusion, a novel skylight orientation sensor is designed and built. Compared with the division-of-focal-plane polarization imaging sensors, the novel skylight orientation sensor has the advantages of no loss of resolution and instantaneous field of view error. The skylight orientation sensor is using the image polarization encoding capability of the S-waveplate and the linear polarizer to convert the skylight polarization information into the extraction of the symmetry axis in the image, which is verified by numerical simulation and outdoor experiments. Under the conditions of clear weather, the compass is used to image the skylight field and the measured angle error after compensation is less than 0.3 deg. The experimental results show that the skylight orientation sensor has good application potential in autonomous navigation.
This skylight orientation sensor needs further research and improvement. First, the angle measurement algorithm under complex weather conditions must be studied to improve the skylight orientation sensor's adaptability and robustness. Second, the skylight orientation sensor needs to be optimized for integrated fabrication and real-time angle measurement.