A Bio-Inspired Polarization Sensor with High Outdoor Accuracy and Central-Symmetry Calibration Method with Integrating Sphere

A bio-inspired polarization sensor with lenses for navigation was evaluated in this study. Two new calibration methods are introduced, referred to as “central-symmetry calibration” (with an integrating sphere) and “noncontinuous calibration”. A comparison between the indoor calibration results obtained from different calibration methods shows that the two proposed calibration methods are more effective. The central-symmetry calibration method optimized the nonconstant calibration voltage deviations, caused by the off-axis feature of the integrating sphere, to be constant values which can be calibrated easily. The section algorithm proposed previously showed no experimental advantages until the central-symmetry calibration method was proposed. The outdoor experimental results indicated that the indoor calibration parameters did not perform very well in practice outdoor conditions. To establish the reason, four types of calibration parameters were analyzed using the replacement method. It can be concluded that three types can be easily calibrated or affect the sensor accuracy slightly. However, before the sensor is used outdoors every time, the last type must be replaced with the corresponding outdoor parameter, and the calculation needs a precise rotary table. This parameter, which is mainly affected by the spectrum of incident light, is the main factor determining the sensor accuracy. After calibration, the sensor reaches an indoor accuracy of ±0.009° and a static outdoor accuracy of ±0.05° under clear sky conditions. The dynamic outdoor experiment shows a ±0.5° heading deviation between the polarization sensor and the inertial navigation system with a ±0.06° angular accuracy.


Introduction
Polarization vision, the ability of some animals to detect the oscillation plane of an electric-field vector of light, has inspired scientists for decades. In 2004, Cochran reported that migrating songbirds recalibrate their internal magnetic compass daily from twilight polarization cues [1]. In 2006, Muheim reported that polarized light cues, in general, underlie the compass calibration of migratory songbirds [2]. The dynamic polarization vision in mantis shrimps was investigated in 2016 [3]. Furthermore, polarization cues also play an important role in homing ants [4]. Polarization navigation based on polarization cues (i.e., the scattered skylight) is an important application of polarization vision. Its advantages are the anti-interference and lack of error accumulations.
A mechanism for bio-inspired polarization sensors for navigation was proposed in [5], and a study of angle measurements using an optoelectronic model to emulate the polarization-sensitive compound eyes of insects was proposed in [6]. Bio-sensors for navigation [7] attempt to imitate the orthogonal sphere, but at some distance away from the port of integrating sphere source. The uniformity is 100% at the plane of port. However, it is almost impossible to position the sensor accurately, because the linear-polarizing laminated films which generate linear polarization light are located at the port, and the sensor has a complex inner structure. Therefore, the measured sensor is often placed at some distance away from the port of integrating sphere, thus decreasing uniformity. This effect leads to systematic errors. To weaken this effect, both the sphere and sensor are generally positioned coaxial. However, in practice, there is always a slight off-center distance. As a result, the section algorithm [30], proposed by us in 2012 to use the redundancy of sensor's light path when detecting the linear polarization light, shows no advantage in indoor calibration experiments. Both the problem of a slight off-center distance and the unavailability of section algorithm are solved in this paper.
In addition, in this study, a polarization bio-sensor with plano-convex lenses was constructed. Because we found that there is no lens in [5,7,9,14], the imitation of compound eyes of insects is not very good. The accuracy limit, determined by a 16-bit analog-to-digital converter (ADC), was calculated in our simulation for different degrees of polarization of incident light.
We also studied and improved the calibration method by exploring two new calibration techniques. To increase the light uniformity decreased by the slight off-center distance when using an integrating sphere, we introduce, what we named as, the "central-symmetry method". The second calibration technique is the "noncontinuous method", which is an extension of the section algorithm and serves to further improve the calibration. Subsequently, indoor calibration experiments were performed and compared, and the effect of different calibration parameters on the accuracy was analyzed. Finally, outdoor experiments were performed, and the calibration results were studied using the replacement method. Table 1 includes the main mathematical variables as below.

Structure
The sensor proposed in this paper mainly consisted of a blue filter, six linear-polarizing laminated films, six plano-convex lenses, six photodiodes, and an electronic circuit (see Figure 1a). The height of the sensor is 100 mm and the diameter is 63 mm. The focal length of the plano-convex lens is 50 mm. The diameter is 0.5 inch (12.7 mm). The photodiode is a ceramic package photodiode. Its spectral response range is 320 to 730 nm. The pixel of the photodiode is 1.3 mm × 1.3 mm. The electronic circuit mainly contains a MSP430 with four 16-bit ADCs inside. The refresh rate of the sensor is 150 Hz. There are six light paths that share the same blue filter. In each path, there is a linear-polarizing laminated film (i.e., polarizer), a plano-convex lens, and a photodiode (see Figure 1b). The plano-convex can narrow the field of view and shield stray light, such as direct sunlight. Inspired by insects, six direct-current signals, which are derived from six photodiodes, are divided into three groups. In each group, the two main polarization directions of linear-polarizing films are orthogonal, and the two current signals act as the inputs to a logarithmic amplifier, which has an exponential response to the sunlight that varies from 0.0000146 to 14.64 W/m 2 when the elevation varies from -6 • to 15 • . The six main polarization angles are 0 • (the reference direction), 90 • , 60 • , 150 • , 120 • , and 210 • (see Figure 1c). The corresponding photodiode numbers are 1, 2, 3, 4, 5, and 6.

Polarization Angle Calculation
The outputs of the three groups [7] are direct voltages. V1_ideal, V2_ideal, and V3_ideal are the ideal outputs of photodiodes 1 and 2, 3 and 4, 5 and 6, respectively. They can be described as follows: where "0.5" is the parameter of logarithmic amplifier; d is the (dimensionless) polarization degree of incident light. The output of sensor, θ, is the angle between the reference direction and the main polarization direction of incident light. The values d = 0.5 and θ = [0°, 180°] with a 0.1° interval for Equation (1) are shown in Figure 2. We can obtain three equal outputs (θ12, θ13, θ23) simultaneously. These are derived from V1_ideal and V2_ideal, V1_ideal and V3_ideal, V2_ideal and V3_ideal, respectively.

Polarization Angle Calculation
The outputs of the three groups [7] are direct voltages. V 1_ideal , V 2_ideal , and V 3_ideal are the ideal outputs of photodiodes 1 and 2, 3 and 4, 5 and 6, respectively. They can be described as follows: where "0.5" is the parameter of logarithmic amplifier; d is the (dimensionless) polarization degree of incident light. The output of sensor, θ, is the angle between the reference direction and the main polarization direction of incident light. The values d = 0.5 and θ = [0 • , 180 • ] with a 0.1 • interval for Equation (1) are shown in Figure 2.

Polarization Angle Calculation
The outputs of the three groups [7] are direct voltages. V1_ideal, V2_ideal, and V3_ideal are the ideal outputs of photodiodes 1 and 2, 3 and 4, 5 and 6, respectively. They can be described as follows: where "0.5" is the parameter of logarithmic amplifier; d is the (dimensionless) polarization degree of incident light. The output of sensor, θ, is the angle between the reference direction and the main polarization direction of incident light. The values d = 0.5 and θ = [0°, 180°] with a 0.1° interval for Equation (1) are shown in Figure 2. We can obtain three equal outputs (θ12, θ13, θ23) simultaneously. These are derived from V1_ideal and V2_ideal, V1_ideal and V3_ideal, V2_ideal and V3_ideal, respectively.  We can obtain three equal outputs (θ 12 , θ 13 , θ 23 ) simultaneously. These are derived from V 1_ideal and V 2_ideal , V 1_ideal and V 3_ideal , V 2_ideal and V 3_ideal , respectively.
where A = 10 2V1_ideal and B = 10 2V2_ideal . In Equation (2), the entire range for θ 12 is 90 • . However, this value should be 180 • for linear polarized light. Hence, a method is needed to extend the range. Each θ 12 has a corresponding d 12 . If the sign of d 12 is negative or positive, then the output angle for V 1_ideal and V 2_ideal is θ 12 +90 • or θ 12 and the final degree of polarization is −d 12 or d 12 . In other words, the sign of d 12 is used to extend the range to 180 • . In a similar manner, we can obtain θ 13 and θ 23 .

Section Algorithm
Three groups of outputs can be obtained, because there are three equations and only two unknown variables. The section algorithm [30] is proposed to use the redundancy of the sensor's light path when detecting the linear polarization light. The basic principle is to avoid using V, for which the corresponding derivative is 0.
The derivatives are 0 when the denominators in Equation (1) are minimal. In this case, the electrical noise produces large errors. For each V, there are two bad points, the derivatives of which are 0. Then, there are six such bad points for our sensor. Therefore, we obtain six sections by averaging the range (0 • -180 • ). The middle point of each section corresponds with the bad point. For each section, V with a bad point is excluded and the other two remaining V are used to calculate the final output θ sec . In this section algorithm, when θ is in the different section, the corresponding θ is considered as the final output. The functional diagram is shown in Figure 3. θ 12 is used to determine the selected section in this paper. For instance, if we assume that θ 12 is 60 • , θ 13 is 60.1 • , and θ 23 is 60.2 • . Then, because θ 12 lies within the section [45 • , 75 • ), the final output θ sec is θ 13 (60.1 • ), and θ 23 is meaningless. More details can be seen in Sections 4.4 and 5.1.
where A = 10 2V1_ideal and B = 10 2V2_ideal . In Equation (2), the entire range for θ12 is 90°. However, this value should be 180° for linear polarized light. Hence, a method is needed to extend the range. Each θ12 has a corresponding d12. If the sign of d12 is negative or positive, then the output angle for V1_ideal and V2_ideal is θ12+90° or θ12 and the final degree of polarization is −d12 or d12. In other words, the sign of d12 is used to extend the range to 180°. In a similar manner, we can obtain θ13 and θ23.

Section Algorithm
Three groups of outputs can be obtained, because there are three equations and only two unknown variables. The section algorithm [30] is proposed to use the redundancy of the sensor's light path when detecting the linear polarization light. The basic principle is to avoid using V, for which the corresponding derivative is 0.
The derivatives are 0 when the denominators in Equation (1) are minimal. In this case, the electrical noise produces large errors. For each V, there are two bad points, the derivatives of which are 0. Then, there are six such bad points for our sensor. Therefore, we obtain six sections by averaging the range (0°-180°). The middle point of each section corresponds with the bad point. For each section, V with a bad point is excluded and the other two remaining V are used to calculate the final output θsec. In this section algorithm, when θ is in the different section, the corresponding θ is considered as the final output. The functional diagram is shown in Figure 3. θ12 is used to determine the selected section in this paper. For instance, if we assume that θ12 is 60°, θ13 is 60.1°, and θ23 is 60.2°. Then, because θ12 lies within the section [45°, 75°), the final output θsec is θ13 (60.1°), and θ23 is meaningless. More details can be seen in Sections 4.4 and 5.1.

Accuracy Limit for Different d Values
The accuracy limit for the sensor depends on the 16-bit ADC as well as the degree of polarization d. The digital output of the 16-bit AD is from −32,767 to +32,767. In our simulation, aimed at calculating the best theoretical accuracy determined by the integer effect, the voltages can be described as follows: where Vi_AD is rounded to the nearest integer by the rule known as round half up and i is 1, 2, or 3. Because Vi_AD includes both positive and negative values, the power of 2 is 15 and not 16 for the 16-bit AD. Vref is the reference voltage for the logarithmic amplifier, and its value is 600 mV in this paper.
The difference between V1_ideal and V1_AD is shown in Figure 4. The horizontal axis represents the theoretical values of linear angle with a 0.1° interval without any errors. The differences between V2_ideal and V2_AD and between V3_ideal and V3_AD are not shown because they look similar.

Accuracy Limit for Different d Values
The accuracy limit for the sensor depends on the 16-bit ADC as well as the degree of polarization d. The digital output of the 16-bit AD is from −32,767 to +32,767. In our simulation, aimed at calculating the best theoretical accuracy determined by the integer effect, the voltages can be described as follows: where V i_AD is rounded to the nearest integer by the rule known as round half up and i is 1, 2, or 3. Because V i_AD includes both positive and negative values, the power of 2 is 15 and not 16 for the 16-bit AD. V ref is the reference voltage for the logarithmic amplifier, and its value is 600 mV in this paper. The difference between V 1_ideal and V 1_AD is shown in Figure 4. The horizontal axis represents the theoretical values of linear angle with a 0.1 • interval without any errors. The differences between V 2_ideal and V 2_AD and between V 3_ideal and V 3_AD are not shown because they look similar.  The error curves were also similar. Therefore, only the error curve derived from V1_AD and V2_AD is shown in Figure 5. The maximum of Vi_AD and the maximum absolute error (accuracy limit) of θ12 that change with d are shown in Figure 6. The accuracy limit for θ13 and θ23 is identical to that for θ12. Therefore, only one curve is shown. Vref is inserted (red) into Figure 6. The corresponding d for Vref is 0.8813. This number is high enough for outdoor measurement [31]. The error curves were also similar. Therefore, only the error curve derived from V 1_AD and V 2_AD is shown in Figure 5.  The error curves were also similar. Therefore, only the error curve derived from V1_AD and V2_AD is shown in Figure 5. The maximum of Vi_AD and the maximum absolute error (accuracy limit) of θ12 that change with d are shown in Figure 6. The accuracy limit for θ13 and θ23 is identical to that for θ12. Therefore, only one curve is shown. Vref is inserted (red) into Figure 6. The corresponding d for Vref is 0.8813. This number is high enough for outdoor measurement [31]. The maximum of V i_AD and the maximum absolute error (accuracy limit) of θ 12 that change with d are shown in Figure 6. The accuracy limit for θ 13 and θ 23 is identical to that for θ 12 . Therefore, only one curve is shown. V ref is inserted (red) into Figure 6. The corresponding d for V ref is 0.8813. This number is high enough for outdoor measurement [31]. The maximum of Vi_AD and the maximum absolute error (accuracy limit) of θ12 that change with d are shown in Figure 6. The accuracy limit for θ13 and θ23 is identical to that for θ12. Therefore, only one curve is shown. Vref is inserted (red) into Figure 6. The corresponding d for Vref is 0.8813. This number is high enough for outdoor measurement [31]. Figure 6. Maximum of Vi_AD and accuracy limit for θ when d is between 0.05 and 0.95. Figure 6. Maximum of V i_AD and accuracy limit for θ when d is between 0.05 and 0.95.

Central-Symmetry Method with Integrating Sphere
In the calibration process, a precise rotary table is used to provide precise angle variations, which generate the accuracy of the polarization sensor. The centers of the integrating sphere (O IS ), the polarization sensor (O S ), and the rotary table (O R ) should coincide, in ideal conditions. The central-symmetry method with the integrating sphere is proposed to improve calibration when the three centers do not coincide in the practical application. Figure 7a shows a light-source diagram, where the black solid arc above the black dashed line represents the integrating sphere with lamps. The black dashed line represents the port of sphere. The blue solid line represents the projection of photosensitive surface of sensor. L o is the length between the port and the photosensitive surface. The irradiance at the center is E o . The irradiance at the off-axis edge is E e . L is the off-axis distance between the off-axis edge and the center. Φ is the intersection angle between the two lines in Figure 7a. Figure 7b shows the eccentric features of the three centers on the photosensitive surface. O R is considered the center. O IS is at the position of E o in Figure 7a. O S is the center of circle, where the six photodetectors are located. The active photodetector area is considered as a point in the simulation. Point 1 represents the position of photodetector 1. L photo is 14 mm. The coordinate system O R -XY is built. Axis X is parallel to Line O S 1. The eccentric distances of the integrating sphere and the sensor are L IS and L S , respectively. The eccentric angles are α IS and α S , respectively. This figure is also the initial state of the anticlockwise rotation in Section 3.2. Figure 7c shows how the central-symmetry method is used in the rotation. Rotating the sensor by 180 • around O R , Points 1 and 2 results in positions 1' and 2'. The off-axis distances between Point 1, 2, 1', and 2' and O IS are L 1 , L 2 , L 1 ', and L 2 ', respectively. The eccentric distances of Points 1 (1') and 2 (2') are L R1 and L R2 , respectively. Assuming that the scattered light in Figure 1b is the light from the sphere, the positions of the six photodiodes with respect to the sphere can be visualized with the aid of Figure 1b,c. Figure 7c shows how the central-symmetry method is used in the rotation. Rotating the sensor by 180° around OR, Points 1 and 2 results in positions 1' and 2'. The off-axis distances between Point 1, 2, 1', and 2' and OIS are L1, L2, L1', and L2', respectively. The eccentric distances of Points 1 (1') and 2 (2') are LR1 and LR2, respectively. Assuming that the scattered light in Figure 1b   According to the integrating sphere manual, the irradiance (W/m 2 ) at the off-axis edge can be expressed as Ee = Eocos 4 φ. Then, the irradiances at Points 1, 2, 1' and 2' are: (4) According to the integrating sphere manual, the irradiance (W/m 2 ) at the off-axis edge can be expressed as E e = E o cos 4 ϕ. Then, the irradiances at Points 1, 2, 1' and 2' are: (4) The voltage derived from the amplifier of photodetectors 1 and 2, before rotation, is V 1br , while the voltage after rotation is V 1ar . Assuming that the spectral responsivity (A/W) and the active area size (m 2 ) for each photodiode are the same, we obtain: The mean value (V 1mr ) of two voltages is used-and not V 1br or V 1ar -for calibration using the central-symmetry method for indoor experiments.

Simulation of the Central-Symmetry Method
The performance of the central-symmetry method, namely the voltage deviation between V mr and V ideal , was studied in detail during the simulation. During the anticlockwise rotation in Figure 7c, the polarization angle, namely the rotational angle (θ R ) of the precise rotary table, kept varying and had a big influence on the voltage deviation. However, the polarization degree (d) had no influence on the voltage deviation, according to Equation (5). Four parameters (L IS , L S , α IS , and α S ) were investigated, and L o was set to 20 mm.
First, the influences of the two eccentric centers (O IS and O S in Figure 7b) on the voltage were analyzed separately. Figure 8a shows the separate influence of O IS , when O S and O R coincide. The voltage deviations between V 1ar and V 1_ideal , and V 1br and V 1_ideal can be considered approximate sine curves. It is easy to conclude that the amplitude depends on the eccentric distance (L IS ) and the phase depends on the eccentric angle (α IS ). The other two groups of voltage deviations were omitted because of the similarity. The voltage deviation between V 1mr and V 1_ideal was close to 0, which proved the validity of the central-symmetry method preliminarily. Figure 8b shows the separate influence of O S , when O IS and O R coincide. The three deviations between V mr and V ideal seemed constant. Under this condition (L IS = 0 mm), V ar was the same as V br , according to Equation (5) . Therefore, they were both the same as V mr . The other two groups of deviations were omitted. phase depends on the eccentric angle (α IS ). The other two groups of voltage deviations were omitted because of the similarity. The voltage deviation between V1mr and V1_ideal was close to 0, which proved the validity of the central-symmetry method preliminarily. Figure 8b shows the separate influence of OS, when OIS and OR coincide. The three deviations between Vmr and Videal seemed constant. Under this condition (LIS = 0 mm), Var was the same as Vbr, according to Equation (5). Therefore, they were both the same as Vmr. The other two groups of deviations were omitted. Second, the integrated influence of the two eccentric centers on the voltage is shown in Figure 9. It was found that Figure 9 was the superposition of Figure 8a,b. The influences of the two eccentric centers were not stuck to each other. On the contrary, they were independent of each other. Second, the integrated influence of the two eccentric centers on the voltage is shown in Figure 9. It was found that Figure 9 was the superposition of Figure 8a,b. The influences of the two eccentric centers were not stuck to each other. On the contrary, they were independent of each other. Third, the analytical expression was used to explain the geometric patterns in Figures 8 and 9, and the general regular was proposed. The deviations between V1mr and V1_ideal in the four calculation steps are shown as Equation (6) Figure 7c, the off-axis distances (L1, L2, L'1, and L'2) can be calculated using the coordinates of the moving points (1 and 2) and the fixed points (OIS) directly, because the moving points can be calculated easily with the constant eccentric distances (LR1 and LR2). Therefore, the second intermediate result was obtained. An approximate assumption in Equation (7) was proposed to obtain the third intermediate result, because LIS was chosen to be 0.02 mm, LR2 was 14 mm, and Lo was 20 mm. Finally, the simplification of the deviation between V1mr and V1_ideal, based on the assumption, was found without θR or α IS . Hence, the deviation was constant, as the solid lines shown in Figure 8a,b and Figure 9. In addition, it is 0 if LR1 is equal to LR2, as shown in Figure 8a.  Third, the analytical expression was used to explain the geometric patterns in Figures 8 and 9, and the general regular was proposed. The deviations between V 1mr and V 1_ideal in the four calculation steps are shown as Equation (6) Figure 7c, the off-axis distances (L 1 , L 2 , L' 1 , and L' 2 ) can be calculated using the coordinates of the moving points (1 and 2) and the fixed points (O IS ) directly, because the moving points can be calculated easily with the constant eccentric distances (L R1 and L R2 ). Therefore, the second intermediate result was obtained. An approximate assumption in Equation (7) was proposed to obtain the third intermediate result, because L IS was chosen to be 0.02 mm, L R2 was 14 mm, and L o was 20 mm. Finally, the simplification of the deviation between V 1mr and V 1_ideal, based on the assumption, was found without θ R or α IS . Hence, the deviation was constant, as the solid lines shown in Figure 8a,b and Figure 9. In addition, it is 0 if L R1 is equal to L R2 , as shown in Figure 8a.
If L 2 Based on Equation (6), there was a general regular: The method optimized the nonconstant voltage deviations generated by the off-axis distance of the integrated sphere during rotation. Constant voltage deviations, which can be calibrated easily, were the final outputs of the central-symmetry method during rotation.
In addition, L R1 and L R2 in Equation (6) were expanded further. Then we obtained: A new assumption, as shown in Equation (9), was proposed to investigate the derivation of Equation (8) when α S is the unknown parameter. As a result, when α S was 240 • , Equation (8) reached the maximum. The three deviations between V mr and V ideal are as shown in Figure 10 when the resolution of α S is 1 • .
Sensors 2019, 19, x FOR PEER REVIEW 11 of 28 A new assumption, as shown in Equation (9), was proposed to investigate the derivation of Equation (8) when α S is the unknown parameter. As a result, when α S was 240°, Equation (8) reached the maximum. The three deviations between Vmr and Videal are as shown in Figure 10 when the resolution of α S is 1°. Figure 10. The influence of α S on the three voltage deviations after using the central-symmetry method.
Fourth, the standard deviations of voltage deviations were used to evaluate the central-symmetry method, which does not change mathematical expectations. During a rotation of 360° with 1° intervals, the number of samples is 360. Standard deviations were investigated when LS and LIS were both between 0.001 and 1 mm. The standard deviations, before and after the central-symmetry method is used, are defined as BSD and ASD, respectively. ASD, the percent of ASD to BSD, the expectation, and the deviation percent of the expectation are shown in Figure 11, when LS is 1 mm and LIS is between 0.001 and 1 mm, and in Figure 12 when LS is between 0.001 and 1 mm and LIS is 1 mm, respectively. Only the first one of the three voltage deviations was chosen to be studied because of the similarity between them. α S was set to be 240° because the corresponding Fourth, the standard deviations of voltage deviations were used to evaluate the central-symmetry method, which does not change mathematical expectations. During a rotation of 360 • with 1 • intervals, the number of samples is 360. Standard deviations were investigated when L S and L IS were both between 0.001 and 1 mm. The standard deviations, before and after the central-symmetry method is used, are defined as BSD and ASD, respectively. ASD, the percent of ASD to BSD, the expectation, and the deviation percent of the expectation are shown in Figure 11, when L S is 1 mm and L IS is between 0.001 and 1 mm, and in Figure 12 when L S is between 0.001 and 1 mm and L IS is 1 mm, respectively. Only the first one of the three voltage deviations was chosen to be studied because of the similarity between them. α S was set to be 240 • because the corresponding voltage deviation was the maximum, as shown in Figure 10. Fourth, the standard deviations of voltage deviations were used to evaluate the central-symmetry method, which does not change mathematical expectations. During a rotation of 360° with 1° intervals, the number of samples is 360. Standard deviations were investigated when LS and LIS were both between 0.001 and 1 mm. The standard deviations, before and after the central-symmetry method is used, are defined as BSD and ASD, respectively. ASD, the percent of ASD to BSD, the expectation, and the deviation percent of the expectation are shown in Figure 11, when LS is 1 mm and LIS is between 0.001 and 1 mm, and in Figure 12 when LS is between 0.001 and 1 mm and LIS is 1 mm, respectively. Only the first one of the three voltage deviations was chosen to be studied because of the similarity between them. α S was set to be 240° because the corresponding voltage deviation was the maximum, as shown in Figure 10.  In Figure 12a, ASD and the percent vary very slightly and can be approximately considered as a constant value. This also happens when LIS is another value between 0.001 and 1 mm. The smaller the LIS, the more precise is this phenomenon. Therefore, it is concluded that LS hardly makes an effect on ASD and the deviation percent when LS and LIS are between 0.001 and 1 mm. Therefore, as shown in Figure 11a, the standard deviations have been largely reduced by the central-symmetry method in the simulation. This result supports the general regular derived from Equation (6).
A similar conclusion to the ASD is made that LIS hardly affects expectations, based on the red curve in Figure 11b. Therefore, Figure 12b shows the influence of LS on the expectations of voltage deviation when LS and LIS are between 0.001 and 1 mm. The expectations will be handled in the following calibration part.
The expectations remain unchanged when the method is used. During the 360° rotation, the expectations of Vbr and Var are the same, because Vbr and Var contain the same elements, which are exactly one period (360°) of a periodic function. The difference is that there is a 180° phase deviation. Then, the expectations of Vbr-Videal and Var-Videal are the same. Therefore, the expectations of Vmr-Videal are the same as Vbr-Videal or Var-Videal. In other words, the central-symmetry method does not change the expectations of voltage deviation.
To further verify Equation (6), the "constant" values calculated using Equation (6) are added to In Figure 12a, ASD and the percent vary very slightly and can be approximately considered as a constant value. This also happens when L IS is another value between 0.001 and 1 mm. The smaller the L IS , the more precise is this phenomenon. Therefore, it is concluded that L S hardly makes an effect on ASD and the deviation percent when L S and L IS are between 0.001 and 1 mm. Therefore, as shown in Figure 11a, the standard deviations have been largely reduced by the central-symmetry method in the simulation. This result supports the general regular derived from Equation (6).
A similar conclusion to the ASD is made that L IS hardly affects expectations, based on the red curve in Figure 11b. Therefore, Figure 12b shows the influence of L S on the expectations of voltage deviation when L S and L IS are between 0.001 and 1 mm. The expectations will be handled in the following calibration part.
The expectations remain unchanged when the method is used. During the 360 • rotation, the expectations of V br and V ar are the same, because V br and V ar contain the same elements, which are exactly one period (360 • ) of a periodic function. The difference is that there is a 180 • phase deviation. Then, the expectations of V br -V ideal and V ar -V ideal are the same. Therefore, the expectations of V mr -V ideal are the same as V br -V ideal or V ar -V ideal . In other words, the central-symmetry method does not change the expectations of voltage deviation.
To further verify Equation (6), the "constant" values calculated using Equation (6) are added to both Figures 11b and 12b. To be exact, the deviation percent between the "constant" values (CV) and the expectations (ME), namely (CV-ME)/ME × 100%, was investigated. The deviation percentages are below 0.23% in Figures 11b and 12b, which verify Equation (6).
At the beginning of this part, L o was set to 20 mm, which is close to the actual value. If its value meets the requirements of Equation (7), the above conclusion of the method can still be made. Because L R2 itself has already met the requirements of Equation (7), L o is allowed to be any real number. We suggest that it should be [0, 40] mm.

Noncontinuous Method
Not all voltages determined using the central-symmetry method were used in calibration parameter estimation. More precisely, voltages with small derivatives (dV/dθ) were rejected, and only the remaining values were used to calculate the calibration parameter. The derivative is given by Equation (10). The voltage (V 1_ideal ) and its derivative are shown in Figure 13 when d is 0.5. Generally, when voltages are close to 0, the corresponding derivatives are close to their extreme values. When voltages are close to their extreme values, the corresponding derivatives are close to 0. The mathematical regular between voltages and their derivatives is similar to a sine function. Note that at least six points remained; otherwise, a poor result will be obtained. More details are shown in Section 4.4.

Calibration Parameters
The incident skylight of the polarization sensor includes two components: unpolarized light and totally polarized light. The two respective gains [32] can be described as follows: where gup represents the gain of unpolarized light; gtp is the gain of totally polarized light. τM is the transmittance when the reference angle and main polarization angle of incident light are parallel. τm is the transmittance when the two directions are orthogonal. Then, the irradiance at one photodiode and the output of one photodiode are as follows:

Calibration Parameters
The incident skylight of the polarization sensor includes two components: unpolarized light and totally polarized light. The two respective gains [32] can be described as follows: where g up represents the gain of unpolarized light; g tp is the gain of totally polarized light. τ M is the transmittance when the reference angle and main polarization angle of incident light are parallel. τ m is the transmittance when the two directions are orthogonal. Then, the irradiance at one photodiode and the output of one photodiode are as follows: where E in is the irradiance of incident light. τ f is the transmittance of the corresponding blue filter. τ c is the transmittance of the corresponding plano-convex lens. s r is the spectral responsivity of photodiode. A r is the active area size of photodiode. Then, V i can be formulated as follows: where j is 2i -1; h is equal to 2i. k v1 is the deviation parameter of reference voltage of logarithmic amplifier. The installation angles of polarizer 0 • and polarizer 90 • are α 1 and α 2 , respectively. k ci is the constant value generated by the integrating sphere method, namely the expectation in Figure 12b. k di is the coefficient of the non-ideal polarizer. Then, Equation (1) can be written as follows: Additionally, k i is: To obtain the final output θ, Equation (15) should be transformed as follows.
By eliminating the unknown number d and expanding all cosine functions, we obtain:

Decoupling Method for Calibration Parameters
The iterative least-squares estimation method was used to obtain the calibration parameters. There are four types of calibration parameters (k, k v , k d , α) (see Equation (15)). Among them, k, k v , and α are easy to obtain. However, the exact value for k d cannot be determined. Clearly, d and (1 + k d ) are coupled, because there are seven unknown numbers (d and k d1-6 ) but only six equations (see Equation (19)). Therefore, the decoupling method is necessary. The coupled parts are as follows: where d cm is obtained using the iterative least-squares estimation method. The value (d a ) for d was determined by the authors and is not an exact value obtained by an algorithm of decoupling method. Note that d a cannot be zero. Then, we can obtain k d .
Then, Equation (15) becomes: The unknown variable in Equation (21) can be considered as θ and d/d a . The inaccurate d a is a part of the unknown variables and does not affect other calibration parameters. Therefore, it can be concluded that the decoupling method does not affect θ but d.

Indoor Results
The indoor calibration experiments were performed and the performances of different methods in Section 3 were analyzed by three different types of calibrations, as shown in Figure 14. Before the least squares were used to calculate the calibration parameters, the three calibrations were distinguished by arrows with different colors. After the least squares, they shared the same process.
The first type of calibration, which is indicated in Section 4.2 and is shown by the blue arrows and the black arrows, only included the section algorithm. The voltages were sent to the least squares directly. The same process shown by the black arrows was a combination of Sections 2.3, 3.4, 3.5 and 4.5. Section 4.5 indicates "Accuracy". Section 2.3 indicates "Section". Sections 3.4 and 3.5 indicate the remaining part. The central symmetry method was inserted into the second type of calibration in Section 4.3, shown by the red arrows. The central symmetry method and the noncontinuous method were both inserted into the third type of calibration in Section 4.4, shown by the green arrows.

Calibration Device
The calibration tools, as shown in Figure 15, include an integrating sphere and a precise rotary table. The integrating sphere provides uniform light, which is unpolarized. A linear polarizer is fixed at the port of the sphere. Then, the light from the sphere becomes linear polarization light after passing through the polarizer. Note that the lamps in the sphere need at least 15 min before using, otherwise the light intensity might not be stable. The model of the sphere is JF-500 (Hefei Xingyue Luminous Technology Applications Institute, Hefei, China). Its diameter is 500 mm. The diameter of its port, where light comes out, is 100 mm. The rotary table, equipped with an optical incremental encoder, has a 1″ accuracy (0.000277°). The central symmetry method was inserted into the second type of calibration in Section 4.3, shown by the red arrows. The central symmetry method and the noncontinuous method were both inserted into the third type of calibration in Section 4.4, shown by the green arrows.

Calibration Device
The calibration tools, as shown in Figure 15, include an integrating sphere and a precise rotary table. The integrating sphere provides uniform light, which is unpolarized. A linear polarizer is fixed at the port of the sphere. Then, the light from the sphere becomes linear polarization light after passing through the polarizer. Note that the lamps in the sphere need at least 15 min before using, otherwise the light intensity might not be stable. The model of the sphere is JF-500 (Hefei Xingyue Luminous Technology Applications Institute, Hefei, China). Its diameter is 500 mm. The diameter of its port, where light comes out, is 100 mm. The rotary table, equipped with an optical incremental encoder, has a 1" accuracy (0.000277 • ).

Calibration Device
The calibration tools, as shown in Figure 15, include an integrating sphere and a precise rotary table. The integrating sphere provides uniform light, which is unpolarized. A linear polarizer is fixed at the port of the sphere. Then, the light from the sphere becomes linear polarization light after passing through the polarizer. Note that the lamps in the sphere need at least 15 min before using, otherwise the light intensity might not be stable. The model of the sphere is JF-500 (Hefei Xingyue Luminous Technology Applications Institute, Hefei, China). Its diameter is 500 mm. The diameter of its port, where light comes out, is 100 mm. The rotary table, equipped with an optical incremental encoder, has a 1″ accuracy (0.000277°).

Section-Only Algorithm
The error curves for θ 12 , θ 13 , θ 23 , and θ sec after calibration, without using the integrating sphere central-symmetry method but by using the section algorithm, are shown in Figure 16.

Section-Only Algorithm
The error curves for θ12, θ13, θ23, and θsec after calibration, without using the integrating sphere central-symmetry method but by using the section algorithm, are shown in Figure 16. The 57 (3 × 19) points (Vmr) were obtained using the central-symmetry method (see Figure 17).

Adding the (Integrating-Sphere) Central-Symmetry Method
The 57 (3 × 19) points (V mr ) were obtained using the central-symmetry method (see Figure 17). Figure 16. Error curves when only the section algorithm is used.

Adding the (Integrating-Sphere) Central-Symmetry Method
The 57 (3 × 19) points (Vmr) were obtained using the central-symmetry method (see Figure 17). The performance of the central-symmetry method is shown in Figure 18. Vicf are the outputs of Equation (21) when the inputs are the angle provided by the precise rotary table. The calibration parameters in Equation (21) are the indoor parameters as shown in Table 2. Vicf are considered as ideal values, corresponding to Vi_ideal in Figure 9, in this performance analysis. The three vertical lines with red, blue, and green color are used to indicate three intersection points (see more details in the discussion). The error curves for θ12, θ13, θ23, and θsec, after calibration, using both the central-symmetry method and section algorithm are shown in Figure 19a. The performance of the central-symmetry method is shown in Figure 18. V icf are the outputs of Equation (21) when the inputs are the angle provided by the precise rotary table. The calibration parameters in Equation (21) are the indoor parameters as shown in Table 2. V icf are considered as ideal values, corresponding to V i_ideal in Figure 9, in this performance analysis. The three vertical lines with red, blue, and green color are used to indicate three intersection points (see more details in the discussion). The error curves for θ 12 , θ 13 , θ 23 , and θ sec , after calibration, using both the central-symmetry method and section algorithm are shown in Figure 19a.  Figure 18. Experimental data analysis when the central-symmetry method is also used. Figure 18. Experimental data analysis when the central-symmetry method is also used. Figure 18. Experimental data analysis when the central-symmetry method is also used. Figure 19. (a) Error curves when the central-symmetry method is also used; (b) error curves when the noncontinuous method is added.

Adding the Noncontinuous Method
The points (Vnc), marked as triangles, were used with the iterative least-squares estimation to obtain the calibration parameters (see Figure 17). The error curves for θ12, θ13, θ23, and θsec after calibration, obtained using the central-symmetry method, noncontinuous method, and section algorithm are shown in Figure 19b. Note that the start point of the  Figures 16-19. According to the actual angles, how the section algorithm works can be seen.

Comparison of Three Calibration Methods
To compare the calibration results, the half value of difference between the maximum and minimum in the error curve is defined as the accuracy of the sensor. For example, the maximum and minimum of black dashed error curve (Error_θsec) in Figure 19

Adding the Noncontinuous Method
The points (V nc ), marked as triangles, were used with the iterative least-squares estimation to obtain the calibration parameters (see Figure 17). The error curves for θ 12 , θ 13 , θ 23 , and θ sec after calibration, obtained using the central-symmetry method, noncontinuous method, and section algorithm are shown in Figure 19b. Note that the start point of the  Figures 16-19. According to the actual angles, how the section algorithm works can be seen.

Comparison of Three Calibration Methods
To compare the calibration results, the half value of difference between the maximum and minimum in the error curve is defined as the accuracy of the sensor. For example, the maximum and minimum of black dashed error curve (Error_θ sec ) in Figure 19

Analysis of Four Calibration Parameters
The calibration parameters derived from Section 4.4 are considered as the indoor calibration parameters because of the highest accuracy. They are shown in the third line in Table 2. The theoretical calibration parameters are shown in the second line.

Analysis of Four Calibration Parameters
The calibration parameters derived from Section 4.4 are considered as the indoor calibration parameters because of the highest accuracy. They are shown in the third line in Table 2. The theoretical calibration parameters are shown in the second line.  Figure 21 shows the accuracies, when parts of the indoor calibration parameters shown in Table 2 were replaced with the corresponding theoretical calibration parameters. For example, the fifth column "k d k v -" represents the accuracy when indoor k d and k v were replaced with the corresponding theoretical calibration parameters (0 and 1), respectively. Furthermore, "-" was used to distinguish the comments between different columns. Considering the accuracy limit in Figure 6, all accuracies are correct to three decimal places when the unit is 1 • . For example, the original accuracy of the indoor calibration was ±0.0088 • , and the accuracy limit was about ±0.001 • when d was~0.7. As a result, the accuracy was ±0.009 • (see the last row in Table 2). The unit for accuracy is 0.001 • in Figure 21. all accuracies are correct to three decimal places when the unit is 1°. For example, the original accuracy of the indoor calibration was ±0.0088°, and the accuracy limit was about ±0.001° when d was ~0.7. As a result, the accuracy was ±0.009° (see the last row in Table 2). The unit for accuracy is 0.001° in Figure 21.

Static Outdoor Experiments
Two groups of static outdoor experiments were performed on the roof (121°32'38''E,

Static Outdoor Experiments
Two groups of static outdoor experiments were performed on the roof (121 • 32'38"E, 38 • 52'45"N) of the School of Mechanical Engineering at Dalian University of Technology on November 16, 2018 (see Figure 22). The sky was clear and the polarization degree was high (from 0.65 to 0.71). The rotary table followed the same process as the indoor experiments. The end time of the first group was 16:11:28 (Beijing time). The end time of the second group was 16:43:08 (Beijing time). Each experiment lasted 108 s for a 360 • range, and the data were calibrated. Note that both the noncontinuous method and the section algorithm were used in the calibration. However, the central-symmetry method was not used because the light source was the light originating from the sky and not from the integrating sphere. The outdoor calibration parameters are shown in the fourth (outdoor1) and fifth columns (outdoor2) of Table 2. The error curves for outdoor2 experiment after calibration are shown in Figure 23. The polarization angle error, the polarization degree, and the standard deviation of 40 samples are shown in Figure 23a-c, respectively. The 40 samples at the start point in Figure 23a-c are shown in Figure 23d. The degree error was ±0.003 • when outdoor2 d a (0.71) in Table 2 was considered the theoretical value. Most standard deviations were smaller than 0.01 • . The error curves for outdoor1 were omitted, because they were very similar. The polarization angle of sky (the reference in this calibration) was assumed to be constant during a short time, such as 54 s for a 180 • range, to ensure that the way to calculate the accuracy is the same as the manner in [7,9,10,12].  Moreover, the outdoor calibration parameters, when the slight variation of the reference was compensated, are shown in the sixth (sun1) and the seventh column (sun2) in Table 2. According to the single-scattering Rayleigh model, the variation of the reference equals the variation of the solar azimuth when the sensor is measuring the zenith.
Note that, similar to Figure    Moreover, the outdoor calibration parameters, when the slight variation of the reference was compensated, are shown in the sixth (sun1) and the seventh column (sun2) in Table 2. According to the single-scattering Rayleigh model, the variation of the reference equals the variation of the solar azimuth when the sensor is measuring the zenith.
Note that, similar to Figure  The same replacement method as in Section 4.6 was used to analyze the four calibration parameters. For example, the fifth column "kvkd" of Figure 24a is the accuracy, when outdoor2 kd and kv were replaced with the corresponding theoretical values (0 and 1), indoor or outdoor1 calibration parameters, respectively. The other parameters (α and k) were still the same for outdoor2. The results are divided into two parts, as shown in Figures 24a and 25a, with respect to the calibration parameter k (i.e., ki), because the accuracy was lower when k was replaced. For all the replacements, all the measured data were the same as those in the experiment at 16:43:08. In Figures 24a and 25a, The same replacement method as in Section 4.6 was used to analyze the four calibration parameters. For example, the fifth column "k v k d " of Figure 24a is the accuracy, when outdoor2 k d and k v were replaced with the corresponding theoretical values (0 and 1), indoor or outdoor1 calibration parameters, respectively. The other parameters (α and k) were still the same for outdoor2. The results are divided into two parts, as shown in Figures 24a and 25a, with respect to the calibration parameter k (i.e., k i ), because the accuracy was lower when k was replaced. For all the replacements, all the measured data were the same as those in the experiment at 16:43:08. In Figures 24a and 25a, the red, green, and blue columns represent that outdoor2 calibration parameters were replaced with the theoretical values, indoor, and outdoor1 calibration parameters, respectively. In Figures 24b and 25b, the red, green, and blue columns represent that sun2 calibration parameters were replaced with the theoretical values, indoor, and sun1 calibration parameters, respectively.
Sensors 2019, 19, x FOR PEER REVIEW 21 of 28 the red, green, and blue columns represent that outdoor2 calibration parameters were replaced with the theoretical values, indoor, and outdoor1 calibration parameters, respectively. In Figures 24b and  25b, the red, green, and blue columns represent that sun2 calibration parameters were replaced with the theoretical values, indoor, and sun1 calibration parameters, respectively.

Dynamic Outdoor Experiments
The dynamic outdoor experiment was performed on the small square roof nearby the stadium of Dalian University of Technology at 19:01 on 1 July 2019 (Beijing time) (see Figure 26). The sky was clear and the polarization degree was high. The polarization sensor and the inertial navigation system with a high accuracy were fixed on the shelf, which has two layers. The shelf was fixed on the wheel barrow, which also has two layers. The path of the dynamic outdoor experiment was a "Z" path. Figure 27 shows the path offered by the inertial navigation system. Figure 28 shows both headings of the polarization sensor and the inertial navigation system, and the partial data is enlarged to show more details. Figure 29 shows their deviations. The polarization degree was about 0.61 during the dynamic outdoor experiment.

Dynamic Outdoor Experiments
The dynamic outdoor experiment was performed on the small square roof nearby the stadium of Dalian University of Technology at 19:01 on 1 July 2019 (Beijing time) (see Figure 26). The sky was clear and the polarization degree was high. The polarization sensor and the inertial navigation system with a high accuracy were fixed on the shelf, which has two layers. The shelf was fixed on the wheel barrow, which also has two layers. The path of the dynamic outdoor experiment was a "Z" path. Figure 27 shows the path offered by the inertial navigation system. Figure 28 shows both headings of the polarization sensor and the inertial navigation system, and the partial data is enlarged to show more details. Figure 29 shows their deviations. The polarization degree was about 0.61 during the dynamic outdoor experiment.          . Heading deviations between the polarization sensor and the inertial navigation system. The inertial navigation system is SPAN-CPT, the price of which is 250,000 RMB or $36,700. SPAN is synchronized position attitude navigation and CPT is a specific type of SPAN. SPAN-CPT brings together two different but complementary positioning and navigation systems: Global Navigation Satellite System (GNSS) and an inertial navigation system (INS). By combining the best aspects of GNSS and INS into one system, SPAN technology offers a solution that is more accurate Figure 29. Heading deviations between the polarization sensor and the inertial navigation system. The inertial navigation system is SPAN-CPT, the price of which is 250,000 RMB or $36,700. SPAN is synchronized position attitude navigation and CPT is a specific type of SPAN. SPAN-CPT brings together two different but complementary positioning and navigation systems: Global Navigation Satellite System (GNSS) and an inertial navigation system (INS). By combining the best aspects of GNSS and INS into one system, SPAN technology offers a solution that is more accurate and reliable than either GNSS or INS alone. The heading accuracy can be ±0.06 • if a good initial calibration is finished. The Y-axis points forward through the front of the vehicle, in the direction the vehicle is moving.
The way to obtain the head using the polarization sensor is the single-scattering Rayleigh model. The e-vector in the body coordinate system and the solar vector in the navigation coordinate system or the horizontal coordinate system are: E b = (cos θ, sin θ, 0) T S n = (cosh s cos As, cosh s sin As, cos hs) T , (22) where As is the solar azimuth and hs is the solar elevating angle. They can be calculated when the geographic position (latitude and longitude) and the time are known. The solar vector in the body coordinate system is: The navigation coordinate system of SPAN-CPT has the definition: Z-axis pointing up; Y-axis pointing north; X-axis pointing east. Therefore, we can obtain: where γ is the roll angle. λ is the pitch angle. H is the heading. Both the roll angle and the pitch angle are offered by SPAN-CPT. Based on the single-scattering Rayleigh model, the e-vector is perpendicular to the solar vector. Then, we can obtain: The heading (H) can be calculated by Equations (22)- (25) because there is only one unknown number.

Discussion
The integer effect caused by the 16-bit ADC was considered in Figures 4 and 5 in our simulation derived from Equation (3). Note that the max error is ±0.02 • when d is equal to 0.05. In practice, electrical noise can strongly reduce the accuracy because the max V is very small (21.73 mV). In other words, the signal-to-noise ratio decreases. Methods to maintain the accuracy as high as when d is large (e.g., 0.5) when d is small will be investigated in a future study, where we will focus on improving the electronic circuit. Figure 18 shows the deviations between V imr and V icf , V iar and V icf , V ibr and V icf . Overall, the experimental results are similar to the simulation results. In particular, the three red curves support the simulation in Figure 9 well. Figure 9 shows one intersection point in the range (0 • -180 • ). The same intersection points are also observed and indicated by vertical lines. The angles of blue, red, and green vertical lines are about 32 • , 98 • , and 150 • , respectively. The deviations are 66 • , 52 • , and 62 • , respectively. This regular is also similar to the relative angular relationship (0 • , 60 • , and 120 • ) of the polarizers in Figure 1c. Therefore, the central-symmetry method performs well.
In addition, k ci in Equation (14), which is not solved by the central-symmetry method, is handled well by the calibration equation (namely, Equation (15)). When L S is 0.02 mm, k ci is between −0.5 and 0.5 mV, as shown in Figure 8. k i is between −10 and 10 mV, as shown in Table 2. Therefore, k ci is a small part of k i . The value of L S depends on the mechanical processing accuracy of the threads that fix the sensor to the rotary table. Generally, the modern mechanical processing accuracy is better than 0.01 mm. Even though L S reaches 1 mm, the central-symmetry method still offers a low ASD, as shown in Figure 12a.
The only differences between the three calibration processes, as shown in Figure 14, are the voltages used in the least-square estimation method. Note that if the number of elements of V inc in Figure 17 is too low, the estimation method cannot provide good calibration parameters, and the accuracy may be very low.
Equation (21) shows that independent of d a , the accuracy remains the same. Therefore, use of the decoupling method to obtain the calibration parameter k d does not change the accuracy. It just changes the degree of polarization slightly. In this paper, d a is slightly larger than d c to ensure that k d is negative, consistent with its physical significance, and the degree of polarization after calibration approaches the real value as close as possible. Figure 20 shows that the values of accuracies of θ 12 , θ 13 , θ 23 , and θ sec decrease from left to right. Therefore, the central-symmetry method (with the integrating sphere) and the noncontinuous method play important roles in the calibration. In the first group of columns, the accuracies of θ 12 , θ 13 , θ 23 , and θ sec are similar, indicating that the section algorithm does not perform well without the central-symmetry method for indoor calibration. The central-symmetry method can be used to calibrate different types of linear polarization devices, not only the sensor used in this paper.
The theoretical-indoor replacement to process the indoor data is shown in Figure 21. When only one type of calibration parameters was replaced, the accuracies can be divided into three groups. k v and k d are in the good group. α is in the modest group. k is in the bad group. In addition, the accuracy of "k v α" or "k d α" column is the same as for "α" column. Similarly, the accuracy of "kk v " or "kk d " column is the same as that for "k" column. Therefore, it can be concluded that k has the strongest effect on accuracy, whereas k v or k d has a minimal effect with regard to the polarization sensor used in this paper.
The "all" column shows that, without calibration, the sensor has a poor initial accuracy (±1.708 • ). When theoretical-outdoor2 replacement was used to process the outdoor data, see Figures 24a and 25a, the same conclusion can be made as for theoretical-indoor replacement.
In addition, when Figures 24a and 25a were analyzed, it was found that the indoor calibration parameters did not perform very well (a ±0.85 • accuracy) in practice outdoor implications (see the green column "all" in Figure 25a). Further analysis was done to study this problem. After analyzing the good group, it was found that the accuracy of "k v k d " column is lower than for "k v " and "k d " columns in Figures 24a and 25a, except for indoor-outdoor2 replacement. This indicates that k v and k d are closely connected and should be processed together.
For the modest group, the calibration parameters are shown in Table 2. Here, α for indoor, outdoor1, and outdoor2 are similar. In addition, the blue and green columns of "α" in Figure 24a show a similar accuracy (±0.036 • ), slightly larger than outdoor2 calibration accuracy (±0.014 • ). Therefore, the directions of linear polarizers were well calibrated for both indoor and outdoor.
For the bad group, the green columns in Figure 25a show a similar accuracy (±0.850 • ), indicating that indoor k is significantly different from outdoor2 k. We suggest that the indoor light source is also significantly different from the scattered skylight. In addition, the blue columns in Figure 25a show a similar accuracy, indicating that outdoor1 k is similar to outdoor2 k. Both Equations (14) and (15) show that k is mainly affected by the spectrum of incident light. Therefore, the spectrum of outdoor light (the scattered skylight) varies very slightly.
Not only the effects of the three groups but also the full effect of the four types of parameters was analyzed. The indoor-outdoor2 replacement values in Figures 24a and 25a indicate that the indoor calibration parameters cannot be used to process outdoor data directly, even though the results are better than theoretical-outdoor2 replacement values. The blue columns in Figures 24a and 25a show the highest accuracy, and the outdoor1 and outdoor2 parameters in Table 2 are very similar. Therefore, it can be concluded that the outdoor calibration parameters can be used to process the subsequently outdoor data, and a high accuracy can be obtained.
As with the analysis of Figures 24a and 25a, in which the sun drift is not compensated, the same conclusion can be made when the sun drift is compensated, because Figure 24a,b are very similar and Figure 25a,b are also very similar. The similarity is an open question. During 54 s, the total sun drifts of "sun1" and "sun2" are 0.14746 • and 0.13792 • , respectively. If the compensation is not very precise, the accuracy after calibration will be decreased. For example, in the "sun2" calibration, when only the start point and the end point are compensated precisely by the total sun drift (0.13792 • ) and the medial points are compensated by linear inserting other than direct compensating, the accuracy is ±0.051 • , which is much bigger than ±0.012 • in Table 2.
Because the outdoor calibration and the indoor calibration obtain similar accuracies as shown in Table 2, we conclude that a clear sky and a high polarization degree is the best condition under which the sensor is used in practice application.
In the outdoor calibration process, every time before the sensor is used outdoors it needs a precise rotary table to calculate four types of parameters. This complex process should be simplified. In particular, when k d and k v are the theoretical parameters (0 and 1), α is the indoor parameter, and k is the outdoor parameter. The static accuracy of outdoor2 data is ±0.050 • . The calculation process of the accuracy (±0.050 • ) is like the blue column of "all" in Figure 25a, with a combination of the theoretical, indoor, and outdoor calibration parameters as the input calibration parameters. Hence, to obtain more accurate outdoor data, we suggest that k d and k v should be the theoretical parameters; α should be the indoor parameter; and k must be the outdoor calibration parameter. This indicates that outdoor calibration is necessary to obtain only k before the measurement, which makes it easier to remove the precise rotary table in the next research work. As a result, ±0.05 • is considered as the actual static accuracy when the polarization bio-sensor is used outdoors. It is better than the outdoor accuracy (between ±0.2 • and ±0.42 • ) of the previous bio-sensors. Therefore, both the plano-convex lens and the calibration achieve substantial improvements.
The dynamic outdoor experiment shows a ±0.5 • heading deviation between the polarization sensor and the inertial navigation system. The data in the two red rectangles in Figure 29 are bigger than ±0.05 • . The reason might be data synchronous deviation. The polarization sensor's frequency is decreased to 70 Hz to avoid loss of data in the transmission process. The inertial navigation system's frequency is 10 Hz. The big angular velocities (about 30 • /s) of the two turns in the experiment may cause big heading deviations when there is a little synchronous deviation, such as 0.1 s.
The deviation above (±0.5 • ) between the two heading measuring systems is bigger than the static experiments, even though the k has been replaced with the outdoor calibration parameter that is offered by outdoor calibration. There are many reasons. First, the observing direction of the polarization sensor, which is the light path as shown in Figure 1b, should coincide with the Z-axis of the inertial navigation system. The reference coordinate system deviation between the polarization sensor and the inertial navigation system needs a good calibration. Second, the polarization sensor should not be single. Many sensors should be used together when measuring different E-vectors, such as the device in the reference [20] which has five sensors. Finally, the wheel barrow, as shown in Figure 26, does not offer good shock absorption, thus the vibration is obvious.
In a future study, the blue filter will be replaced with a narrower band filter. Further optimization will be performed using the results described in this paper, which makes it possible to extend the effective range of calibration parameters while maintaining a high accuracy. New dynamic outdoor experiments, which are designed and calibrated well, will be performed.

Conclusions
A navigation bio-sensor that uses plano-convex lenses and two new calibration methods are described in this paper. The two methods are the central-symmetry calibration method (with an integrating sphere) and the noncontinuous calibration method. After calibration, the bio-sensor shows ±0.009 • indoor accuracy and ±0.018 • outdoor accuracy for clear-sky conditions. If more bio-sensors are used together, an even better accuracy might be obtained. In addition, the bio-sensor has a higher refresh rate (150 Hz) than the 30 Hz rate of conventional imaging polarimetry. The plano-convex lens is an essential component in this bio-sensor. The central-symmetry calibration method prominently decreased the standard deviation of voltage deviations during calibration and generated a constant value that could be easily calibrated. As a result, the redundancy of the sensor's light path has advantages when detecting the linear polarization. The calibration parameters that were mainly determined by the spectrum of incident light dominate the accuracy, according to the analysis of calibration parameters. Before the sensor is used outdoors every time, the indoor calibration parameters k must be replaced with the corresponding outdoor parameters, and a static accuracy of ±0.05 • can be regarded as the actual static outdoor accuracy. The dynamic outdoor experiment shows a ±0.5 • heading deviation between the polarization sensor and the inertial navigation system.