SkyPole: a geolocation algorithm based on polarized vision without using astronomical ephemerides

Global Navigation Satellite Systems (GNSS) are widely used due to their easy access to outdoor GNSS signals and their spatial precision. However, such systems are sensitive to jamming and spoofing. Simple and robust navigation strategies can be found in animals deprived by essence of any GNSS system. Studies have shown that animals like bees or ants utilize the sky’s polarization pattern for navigation. We recently proposed a method inspired by migratory birds, which calibrate their magnetic compass through the celestial rotation of night stars or the daytime polarization pattern. By considering the temporal properties of the sky’s polarization pattern as a relevant navigation information, we developed a bio-inspired method to find the geographical north bearing and the observer’s latitude, requiring only skylight polarization observations during day. To reduce the noise susceptibility of our method, we added a pre-processing step using a least square method based on skylight polarization models, and a segmentation process based on a convolutional autoencoder neural network, trained with simulated data.


INTRODUCTION
Global Navigation Satellite Systems (GNSS) are widely used due to their easy access to outdoor GNSS signals and their spatial precision.However, such systems are sensitive to jamming and spoofing.As an alternative, bio-inspired methods have been developed.Studies have shown that some insects, particularly desert ants, utilize the sky's polarization pattern for navigation. 1,2 he use of solar ephemeris combined with the estimation of the Sun's position through the polarization pattern enables direct geographical positioning. 3Nevertheless, as far as we know, animals do not have access to these astronomical ephemerides, and the utilization of the polarization pattern as a reference for their navigation remains poorly understood.We recently proposed an alternative method inspired by migratory birds, 4 which calibrate their magnetic compass through the celestial rotation of night stars or the daytime polarization pattern. 5,6 s suggested by Brines, 7 we considered the temporal properties of the sky's polarization pattern as a relevant navigation information.We developed a bio-inspired algorithm to find the geographical north bearing and the observer's latitude, requiring only skylight degree of linear polarization (described in the following section) observations by day, and provided by a commercial polarimetric camera.No additional information, such as ephemerides or time was needed.However, the proposed algorithm remains sensitive to noise.We therefore added a pre-processing step using a least square method based on skylight polarization models, and a segmentation process based on a convolutional autoencoder neural network, trained with simulated data.We will first present our method, then our results, and we will finally discuss about the level of performance of our algorithm.

MATERIALS AND METHODS
As the Earth rotates around its North-South Axis, the sky rotates uniformly around a point called the celestial pole.This can be observed during night, stars rotating around the celestial pole, or during day with the Sun.However, the Sun being a single point in the sky, and at a large angular distance from the celestial pole, estimating the celestial pole's position from the Sun's rotation may result in large errors, or may not be possible when the Sun is not visible.An alternative is to use the skylight polarization pattern, linked to the Sun's position in the sky.When sun rays enter the atmosphere, they interact with particles which scatter light in all directions.Scattered light from the atmosphere features polarization properties (angle and degree) depending on the position of the particle relative to the Sun, and the direction of the scattered light.From those scattered rays, a pattern of polarization appears in the sky.This phenomenon can be described by using the Rayleigh single scattering model. 8The skylight polarization, which is mostly linearly polarized, can be described by two parameters, the direction of polarization and the degree of linear polarization (DoLP).In the following, we will only consider the DoLP, which is defined as the ratio between the polarized light intensity and the total light intensity.From Rayleigh single scattering model, the skylight's DoLP is given by: 8 Where γ is the scattering angle, defined in fig 1, and DoLP max the maximum degree of linear polarization.In, 4 we described how the sky DoLP pattern varies in time.We demonstrated that when the Sun moves from a position to another, the sky DoLP pattern remains unchanged on two circles which we called radial and plane invariances.We also shown that the celestial pole is located on the radial invariance, and is the the only point in the sky that has a constant DoLP value during a whole day.To estimate the position of the celestial pole, the SkyPole algorithm segments the DoLP invariances from pairs of sky DoLP images, and induces the position of the celestial pole from the intersection of all the invariances.However, when applying this technique, multiple intersections of invariances can appear, adding ambiguity to the result.To overcome this issue, the radial invariance can be distinguished from the plane invariance.In order to isolate the radial invariance, we trained a U-Net neural network, 9, 10 taking as input two DoLP images, and returning a segmentation image of the radial invariance.To simplify the training task, we trained our network on simulated images only, using the OpenSky simulator. 11Radial invariance segmented images were generated using the analytical equation of this set (see supplemental information 4 ).Real images where pre-processed using a least square method, in order to match a simulated image of the DoLP, based on Rayleigh single scattering model, to the real image.The method is described in figure 1.

RESULTS
We tested our method on data recorded with a division of focal plane polarimetric camera (PHX050S-QC from Lucid Vision Labs, Sony IMX250MYR) equipped with a 185°fisheye lens (FE185C57HA-1, Fujinon), which was situated on the roof of a laboratory (Institut des Neurosciences de la Timone, INT UMR7289) in Marseille, France (43.2870135°N, 5.4034385°E). 12The pose of the camera (3D orientation) was calibrated by means of a method based on measured position of the sun in the camera's frame compared with a theoretical position of the sun provided by ephemeris.Processed data was acquired during the whole month of August 2022, each day from 8:00am to 6:00pm. 13Only data with a DoLP max higher than 0.2 was processed.Since we were located on the Earth north hemisphere, we detected the position of the north celestial pole (NCP).Each estimation of the NCP was made using four images taken with time intervals of 100 minutes.Once the celestial pole's position was estimated, we computed the True North and the camera's latitude from the coordinates of the celestial pole.We obtained a mean absolute error of 2.6194 degrees in True North and of 1.8059 degrees in latitude of the camera.Results are displayed in figure 2.
We also tested our algorithm on data simulated using Rayleigh single scattering model and a camera model equivalent to the one used for experimentation.Data was simulated with random dates between the 1st January 2015 and the 1st January 2018.For each date, a second and a third date were randomly selected in a range varying between 0.5 and 2h from the first and second date respectively.For each set of three dates, a random latitude and longitude on Earth was generated (also on the south hemisphere, but by estimating the position of the south In the first row are the DoLP patterns taken at 3 different times (more images will improve the results).Least square method was applied to match simulated images to the real images, giving the second row.The radial invariances were then segmented by using a U-Net neural network, the results being presented in the third row.Lastly, the radial invariances images were summed, and the NCP was then located at the intersection between the radial invariances, which has the highest value in terms of pixel intensity.
Proc. of SPIE Vol.13050 1305007-3 Figure 2. North Celestial Pole (NCP) coordinates estimation error.NCP coordinates were computed from experimental data with the proposed algorithm (Fig. 1d).δsun is the solar declination in degrees.∆αNCP and ∆θNCP are the azimuth and altitude errors, respectively, in degrees, of the NCP measured with respect to the ground truth values.nmeas is the number of measurements for each error interval of the histogram.The zoomed in histogram corresponds to NCP coordinates values with absolute error lower than 5 and 10 degrees for azimuth and latitude respectively.It represents 162 measurements out of 226 in azimuth, and 217 measurements out of 226 in elevation.For each estimation, we computed our algorithm using four images taken with time intervals of 100 minutes celestial pole).For each date, latitude and longitude, the Sun's position was obtained from solar ephemeris, using Python libraries suncalc and PyAstronomy. 14Then, Rayleigh single scattering model was simulated using the OpenSky simulator. 11For each set of three images, we computed the position of the celestial pole using the proposed algorithm, without the least square step, which was skipped.We observe that some errors occurs near the poles or the equator.Near the poles, the celestial pole is close to the zenith.Since our simulated camera is placed horizontally, we obtain high errors due to the singularity in azimuth of the camera near the zenith.Near the equator, since the celestial pole is at the horizon, we may sometimes detect the opposite intersection of radial invariances leading to 180°errors in azimuth.We therefore decided to quantitatively present our results before and after removing outliers with Matlab function rmoutliers (default parameters).We obtained 1879 measurements after removing outliers from a total of 2078 measurements.We analysed our results with and without outliers using Matlab rmoutliers function.We finally obtained a mean absolute error (MAE) of 1.0273 degrees in azimuth (3.5528 degrees with outliers), and 0.6752 degrees in azimuth (1.0173 degrees with outliers).Results are plotted on figure 3 4. DISCUSSION In this study, we presented a new algorithm to implement the SkyPole method.This algorithm has the advantage of removing the ambiguity in the celestial pole estimation due to multiple intersections of plane and radial invariances, when using for instance only 3 images, or during short periods of time.Therefore, contrary to the SkyPole algorithm, presented in, 4 where at least four DoLP images were required to estimate the position of NCP without ambiguity, three DoLP images is now sufficient.Another advantage of this algorithm is that, by using a neural network, the algorithm needs no tuning of its parameters, contrary to the SkyPole algorithm, 4 where a threshold has to be tuned, ideally as a function of the time interval between two images, to obtain a suitable width of the radial and plane invariances.Moreover, the least square step of our algorithm is robust to noise such as slighly cloudy skies, allowing to compute the SkyPole method even when clouds are located near the celestial pole.However, some issues remain.First, when using the least square method as a prepossessing step in our algorithm, the presence of buildings in the field of view, modifying the polarization properties of light, or overcast skies causing the overall DoLP to drop, highly reduce the accuracy of the algorithm.Here, we only processed images where DoLP max was higher than 0.2, removing data with overcast sky.Second, when the time intervals between three or more images, used to compute at least two radial invariances, are too short, the radial invariances overlap on a large region, reducing the celestial pole's estimation accuracy.To get precise results, the width of the radial invariances simulated to train the neural network should be as little as possible.

CONCLUSION
In this study, we presented a new way of processing the SkyPole method.The proposed algorithm allows us to remove an eventual ambiguity due to multiple intersections between the plane and radial DoLP time invariances.Moreover, by using a neural network, the threshold used for the segmentation of the radial invariance is tuned automatically, meaning the algorithm's parameters should not be modified manually by its user.Next, the proposed algorithm can help to reduce the acquisition time of the SkyPole method.Finally, simulated data shows that our method can be used either on the north hemisphere or on the south hemisphere, by detecting the north or south celestial pole respectively.

Figure 1 .
Figure 1.(a) Scattering angle γ, azimuth αP of a point P and altitude θS of the sun S. Parameters are presented in the ENU frame, namely East, North, and Up frame, centered on the observer O.The colored pattern stands for the skylight DoLP as described in equation 1. Dark blue corresponds to near zero DoLP and yellow, to maximum DoLP values (1 in theory, less in reality).(b) The trajectory of the sun in the ENU frame, centered on an observer O located at latitude ϕ.NCP is the North Celestial Pole.The sun moves on a plane perpendicular to the observer-NCP vector.(c) Invariance axis on the celestial sphere.Comparison between simulated and analytical sets of solutions.Green circle is the radial invariance circle, red circle is the plane invariance circle computed from analytical calculations (see supplemental information 4 ).Colored half sphere is the simulated absolute difference between two DoLP patterns associated with the Sun's positions S1 and S2 at two different times.Dark blue corresponds to near zero values.Red dot is the NCP.(d) Method for finding the NCP based on the skylight's DoLP pattern.In the first row are the DoLP patterns taken at 3 different times (more images will improve the results).Least square method was applied to match simulated images to the real images, giving the second row.The radial invariances were then segmented by using a U-Net neural network, the results being presented in the third row.Lastly, the radial invariances images were summed, and the NCP was then located at the intersection between the radial invariances, which has the highest value in terms of pixel intensity.

Figure 3 .
Figure 3. Celestial Pole (CP) coordinates estimation error (north celestial pole for positive latitudes, and south celestial pole for negative latitudes).CP coordinates were computed from simulated data with the proposed algorithm (Fig.1d), by skipping the least square step.Φcam is the camera latitude.αCPand θCP are the azimuth and altitude, respectively, in degrees, of the measured and groundtruth CP. ∆αCP and ∆θCP are the azimuth and altitude errors, respectively, in degrees, of the CP measured with respect to the ground truth values.nmeas is the number of measurements for each error interval of the histogram.The zoomed in histograms corresponds to CP coordinates values with absolute error lower than 5 degrees.It represents 1926 measurements out of 2078 in azimuth, and 2058 measurements out of 2078 in elevation.