2.1 Altitude model

As shown in Figure 1, O and NP are assumed to be the centre and north pole of Mars, respectively, and P is assumed to be the position of the rover. An altitude h ₁ is obtained by observing a star, the zenith of which is z ₁ = 90° – h ₁. At the epoch of the observation, the star's projection point O₁ on the surface of Mars can be calculated from the star catalogue. Thus, an altitude observation is equivalent to determining a position circle on the surface of Mars with a centre O₁ and a radius of the star's zenith z ₁. If two stars are observed, two position circles can be determined, and one of the two intersections between these two position circles is the final position of the rover. Additionally, ϒ is the apparent vernal equinox on the equator referenced to define the star's apparent right ascension α.

Figure 1. Altitude model for celestial positioning

2.2 Positioning equation and solution

To solve for the position of the rover, the positioning equation is introduced as Equation (1), which shows the relationship between the longitude λ, latitude φ, apparent right ascension of the star α, apparent declination of the star δ, altitude of the star h, and prime meridian sidereal time S as follows:

(1)\begin{equation}\sin h = \sin \delta \sin \varphi + \cos \delta \cos \varphi \cos (S + \lambda - \alpha )\end{equation}

α and δ can be calculated based on the star catalogue and solar system ephemeris in combination with the observation epoch. S can be calculated based only on the observation epoch. Both the International Astronomical Union and United States Naval Observatory (USNO) have published open-source programs for the precise calculation of α, δ and S (Bangert et al., Reference Bangert, Puatua, Kaplan, Bartlett and Monet2009). Ideally, if the optical axis of the star sensor coincides with the local plumb line, the star's altitude h can be measured directly. Therefore, the altitudes of only two stars are needed to determine the position of the rover. However, because a star sensor can measure the altitude angles of dozens or hundreds of stars in one star image, there are always redundant observations. Thus, the least square method is usually used for optimal estimation of the longitude λ and latitude φ.

The accuracies of α and δ can reach the order of milliseconds, and thus the error can be ignored (Bangert et al., Reference Bangert, Puatua, Kaplan, Bartlett and Monet2009); moreover, the accuracy of S depends on the time record, and the current timing accuracy can similarly reach the order of milliseconds, so this error can also be ignored. Therefore, the celestial positioning accuracy depends mainly on the measurement accuracy of the star's altitude.

3.1 Imaging model

In this paper, a fish-eye lens and a charged couple device (CCD) are integrated as a star sensor. A star sensor with a super-large field of view was chosen to reduce the weight of the Mars rover through multiplexing technology. Due to an optical filter setup between the lens and the CCD, the star sensor can observe both the stars and the sun by removing and adding the filter. This means that no additional rotating platform is needed to search for the sun in the daytime. The detailed configuration of the star sensor is introduced in section 5.

As depicted in Figure 2, ${z_S}$ is the half-angle of view of a star; r is the pixel distance between the image centroid of a star point and the image principal point; $(x,y)$ are the coordinates of the image centroid of a star point; $({x_0},{y_0})$ are the coordinates of the image principal point; f is the focal length; $({k_1},{k_2},{k_3})$ are radial distortion parameters; and $({X_0},{Y_0},{Z_0})$ are the coordinates of the optical centre in the horizontal frame O-X_HY_HZ_H. Additionally, $(a,b,c)$ is used to describe the attitude of the star sensor in the horizontal frame. Since the fish-eye lens utilises a solid angle projection, ${z_S}$ is written, as in Zhan et al. (Reference Zhan, Zheng, Li, Wang, Zhu and Chen2020), as

(2)\begin{equation}{z_S} = 2\arcsin \left( {\frac{r}{{2f}}} \right) + {k_1}{\left( {\arcsin \left( {\frac{r}{{2f}}} \right)} \right)^2} + {k_2}{\left( {\arcsin \left( {\frac{r}{{2f}}} \right)} \right)^3} + {k_3}{\left( {\arcsin \left( {\frac{r}{{2f}}} \right)} \right)^4}\end{equation}

where r can be calculated by

(3)\begin{equation}r = \sqrt {{{(x - {x_0})}^2} + {{(y - {y_0})}^2}}\end{equation}

Figure 2. Imaging model of the star sensor with a fish-eye lens

If the optical axis of the star sensor is rigorously opposite to the local gravity vector g, then ${z_S}$ becomes the zenith distance of the star. In practice, the optical axis always deviates from the gravity vector due to the random attitude of the rover, so an inclinometer is necessary for a tilt correction. Additionally, the expression of the star's azimuth A_S in the star sensor frame O-X_SY_SZ_S is given as

(4)\begin{equation}{A_S} = \begin{cases} {\textrm{atan}\dfrac{{y - {y_0}}}{{{x_0} - x}}}& {{x_0} - x > 0,y - {y_0} \ge 0}\\ {\textrm{atan}\dfrac{{y - {y_0}}}{{{x_0} - x}} + \mathrm{\pi }}& {{x_0} - x < 0}\\ {\textrm{atan}\dfrac{{y - {y_0}}}{{{x_0} - x}} + 2\mathrm{\pi }}& {{x_0} - x > 0,y - {y_0} < 0} \end{cases}\end{equation}

Additionally, when ${x_0} = x,y = {y_0}$, ${A_S} = {0^ \circ }$; when ${x_0} = x,y - {y_0} > 0$, ${A_S} = {90^ \circ }$; and when ${x_0} = x,y - {y_0} < 0$, ${A_S} = {270^ \circ }$. $({x_0},{y_0},f,{k_1},{k_2},{k_3})$ are called the inner parameters of the star sensor, which are usually obtained by an indoor calibration before the rover is launched. During the long operation of the Mars rover, the star sensor must work in a complex environment with various vibrations and fluctuations in the temperature and atmospheric pressure, all of which will inevitably change the inner parameters and introduce measurement errors into the zenith distance and altitude angle according to Equations (2) and (4), respectively. Therefore, it is particularly important to carry out research on self-calibrating the star sensor.

3.2 Self-calibration equation

Considering that all the stars are far enough from the solar system, $({X_0},{Y_0},{Z_0})$ can be set to zero, and the distances between the optical centre and the stars are taken as a unified length. Then, $({A_H},{z_H})$ is used to express the coordinates of the star in the horizontal frame O-X_HY_HZ_H, which can be calculated by a computer according to the star catalogue, solar system ephemeris, observation epoch, and initial longitude and latitude. The star vectors in the horizontal frame and the star sensor frame are

(5)\begin{align}{\textbf{{S}}_S} & = \left[ \begin{array}{l} {X_S}\\ {Y_S}\\ {Z_S} \end{array} \right] = \left[ \begin{array}{l} \sin {z_S}\cos {A_S}\\ \sin {z_S}\sin {A_S}\\ \cos {z_S} \end{array} \right]\end{align}

(6)\begin{align}{\textbf{{S}}_H} & = \left[ \begin{array}{l} {X_H}\\ {Y_H}\\ {Z_H} \end{array} \right] = \left[ \begin{array}{l} \sin {z_H}\cos {A_H}\\ \sin {z_H}\sin {A_H}\\ \cos {z_H} \end{array} \right]\end{align}

Then, a rotation matrix R is introduced to express the relationship between the horizontal frame and the star sensor frame:

(7)\begin{equation}{\textbf{{S}}_H} = \textbf{{R}} \cdot {\textbf{{S}}_S}\end{equation}

According to attitude determination theory, at least three parameters are essential to depict R. In this paper, three Rodrigues parameters $(a,b,c)$ are used instead of Euler angles and quaternions. Any rotation matrix can be decomposed into a unit matrix I and an antisymmetric matrix Q (Yao et al., Reference Yao, Han and Yang2006):

(8)\begin{equation}\textbf{{R}} = {({\rm I} - \textbf{{Q}})^{ - 1}}({\rm I} + \textbf{{Q}})\end{equation}

where

(9)\begin{align}\textbf{{Q}} & = \left[ {\begin{array}{*{20}{c}} 0& { - c}& { - b}\\ c& 0& { - a}\\ b& a& 0 \end{array}} \right]\end{align}

(10)\begin{align}{\rm I} & = \left[ {\begin{array}{*{20}{c}} 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1 \end{array}} \right]\end{align}

Furthermore, Equation (8) can be written as

(11)\begin{equation}({\rm I} - \textbf{{Q}}){\textbf{{S}}_H} = ({\rm I} + \textbf{{Q}}) \cdot {\textbf{{S}}_S}\end{equation}

Then, Equation (9) and Equation (10) are inserted into Equation (11) and the basic observation equation of self-calibration is obtained as follows:

(12)\begin{equation}\left[ {\begin{array}{*{20}{c}} 0& { - {Z_S} - {Z_H}}& { - {Y_S} - {Y_H}}\\ { - {Z_S} - {Z_H}}& 0& {{X_S} + {X_H}}\\ {{Y_S} + {Y_H}}& {{X_S} + {X_H}}& 0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} a\\ b\\ c \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{X_H} - {X_S}}\\ {{Y_H} - {Y_S}}\\ {{Z_H} - {Z_S}} \end{array}} \right]\end{equation}

Considering the random measurement noise of the star sensor, the error equation can be written as

(13)\begin{equation}\left[ {\begin{array}{*{20}{c}} {{v_1}}\\ {{v_2}}\\ {{v_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0& { - {Z_S} - {Z_H}}& { - {Y_S} - {Y_H}}\\ { - {Z_S} - {Z_H}}& 0& {{X_S} + {X_H}}\\ {{Y_S} + {Y_H}}& {{X_S} + {X_H}}& 0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} a\\ b\\ c \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{X_H} - {X_S}}\\ {{Y_H} - {Y_S}}\\ {{Z_H} - {Z_S}} \end{array}} \right]\end{equation}

where v ₁, v ₂ and v ₃ are the residuals caused by random measurement noise and the residual vector is expressed as $\textbf{v} = {\left[ {\begin{array}{*{20}{c}} {{v_1}}& {{v_2}}& {{v_3}} \end{array}} \right]^\textrm{T}}$. Obviously, if one star image contains n star points, 3n error equations will be obtained, and the least square method is utilised to estimate the optimal star sensor parameters $X = {({x_0},{y_0},f,{k_1},{k_2},{k_3},a,b,c)^T}$

3.3 Partial derivatives and solutions

Equation (13) is nonlinear, so it needs to be linearised to solve for the parameters iteratively. Thus, the partial derivatives of v with respect to X are deduced in detail as follows:

1. 1. The partial derivatives of $({A_S},{z_S})$ with respect to X are written as

   (14)\begin{align}& \left\{ \begin{array}{@{}l} \dfrac{{\partial {A_S}}}{{\partial {x_0}}} = \dfrac{{{y_0} - y}}{{{r^2}}}\\ \dfrac{{\partial {A_S}}}{{\partial {y_0}}} = \dfrac{{x - {x_0}}}{{{r^2}}} \end{array} \right.\end{align}
   (15)\begin{align}& \left\{ \begin{array}{@{}l} \dfrac{{\partial {z_S}}}{{\partial {x_0}}} = (2 + 2{k_1}T + 3{k_2}{T^2} + 4{k_3}{T^3})\dfrac{{\partial T}}{{\partial r}}\dfrac{{\partial r}}{{\partial {x_0}}}\\ \dfrac{{\partial {z_S}}}{{\partial {y_0}}} = (2 + 2{k_1}T + 3{k_2}{T^2} + 4{k_3}{T^3})\dfrac{{\partial T}}{{\partial r}}\dfrac{{\partial r}}{{\partial {y_0}}} \end{array} \right.\end{align}
   where
   (16)\begin{equation}\left\{ \begin{array}{l} T = \arcsin \dfrac{r}{{2f}}\\ \dfrac{{\partial T}}{{\partial r}} = \dfrac{1}{{\sqrt {4{f^2} - {r^2}} }} \end{array} \right.\end{equation}
   and
   (17)\begin{equation}\left\{ \begin{array}{l} \dfrac{{\partial r}}{{\partial {x_0}}} = \dfrac{{{x_0} - x}}{r}\\ \dfrac{{\partial r}}{{\partial {y_0}}} = \dfrac{{{y_0} - y}}{r} \end{array} \right.\end{equation}

The partial derivatives of $({A_S},{z_S})$ with respect to f are written as

(18)\begin{align}\frac{{\partial {A_S}}}{{\partial f}} & = 0\end{align}

(19)\begin{align}\frac{{\partial {z_S}}}{{\partial f}} & = (2 + 2{k_1}T + 3{k_2}{T^2} + 4{k_3}{T^3})\frac{{\partial T}}{{\partial f}}\end{align}

where

(20)\begin{equation}\frac{{\partial T}}{{\partial f}} = \frac{{ - r}}{{f\sqrt {4{f^2} - {r^2}} }}\end{equation}

The partial derivatives of $({A_S},{z_S})$ with respect to $({k_1},{k_2},{k_3})$ are written as

(21)\begin{align}& \frac{{\partial {A_S}}}{{\partial {k_1}}} = \frac{{\partial {A_S}}}{{\partial {k_2}}} = \frac{{\partial {A_S}}}{{\partial {k_3}}} = 0\end{align}

(22)\begin{align}& \left\{ \begin{array}{l} \dfrac{{\partial {z_S}}}{{\partial {k_1}}} = {T^2}\\ \dfrac{{\partial {z_S}}}{{\partial {k_2}}} = {T^3}\\ \dfrac{{\partial {z_S}}}{{\partial {k_3}}} = {T^4} \end{array} \right.\end{align}

1. 2. The partial derivatives of ${\textbf{{S}}_S}$ with respect to X are written as

   (23)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial {x_0}}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}& { - \sin {z_S}\sin {A_S}}\\ {\cos {z_S}\sin {A_S}}& {\sin {z_S}\cos {A_S}}\\ { - \sin {z_S}}& 0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {z_S}}}{{\partial {x_0}}}}\\ {\dfrac{{\partial {A_S}}}{{\partial {x_0}}}} \end{array}} \right]\end{align}
   (24)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial {y_0}}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}& { - \sin {z_S}\sin {A_S}}\\ {\cos {z_S}\sin {A_S}}& {\sin {z_S}\cos {A_S}}\\ { - \sin {z_S}}& 0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {z_S}}}{{\partial {y_0}}}}\\ {\dfrac{{\partial {A_S}}}{{\partial {y_0}}}} \end{array}} \right]\end{align}
   (25)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial f}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}\\ {\cos {z_S}\sin {A_S}}\\ { - \sin {z_S}} \end{array}} \right]\frac{{\partial {z_S}}}{{\partial f}}\end{align}
   (26)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial {k_1}}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}\\ {\cos {z_S}\sin {A_S}}\\ { - \sin {z_S}} \end{array}} \right]\frac{{\partial {z_S}}}{{\partial {k_1}}}\end{align}
   (27)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial {k_2}}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}\\ {\cos {z_S}\sin {A_S}}\\ { - \sin {z_S}} \end{array}} \right]\frac{{\partial {z_S}}}{{\partial {k_2}}}\end{align}
   (28)\begin{align}\frac{{\partial {\textbf{{S}}_S}}}{{\partial {k_3}}}& =\left[ {\begin{array}{*{20}{c}} {\cos {z_S}\cos {A_S}}\\ {\cos {z_S}\sin {A_S}}\\ { - \sin {z_S}} \end{array}} \right]\frac{{\partial {z_S}}}{{\partial {k_3}}}\end{align}

For the sake of simplicity, the partial derivatives of ${\textbf{{S}}_S}$ with respect to X are defined as

(29)\begin{equation}\frac{{\partial {\textbf{{S}}_S}}}{{\partial \textbf{{X}}}}={\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {x_0}}}}& {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {y_0}}}}& {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial f}}}& {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_1}}}}& {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_2}}}}& {\begin{array}{*{20}{c}} {\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_3}}}}& 0& 0& 0 \end{array}} \end{array}} \right]^\textrm{T}}\end{equation}

1. 3. The partial derivatives of residual vector v with respect to X are written as

   (30)\begin{align}\dfrac{{\partial \textbf{{v}}}}{{\partial \textbf{{X}}}} & = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {v_1}}}{{\partial \textbf{{X}}}}}& {\dfrac{{\partial {v_2}}}{{\partial \textbf{{X}}}}}& {\dfrac{{\partial {v_3}}}{{\partial \textbf{{X}}}}} \end{array}} \right]\notag\\ & = \left[ \begin{array}{l} \begin{array}{*{20}{c}} {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {x_0}}}}& {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {y_0}}}}& {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial f}}}& {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_1}}}}& {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_2}}}}& {{\textbf{{l}}_1}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_3}}}}& 0& { - {Z_S} - {Z_H}}& { - {Y_S} - {Y_H}} \end{array}\\ \begin{array}{*{20}{c}} {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {x_0}}}}& {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {y_0}}}}& {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial f}}}& {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_1}}}}& {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_2}}}}& {{\textbf{{l}}_2}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_3}}}}& { - {Z_S} - {Z_H}}& 0& {{X_S} + {X_H}} \end{array}\\ \begin{array}{*{20}{c}} {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {x_0}}}}& {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {y_0}}}}& {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial f}}}& {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_1}}}}& {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_2}}}}& {{\textbf{{l}}_3}\dfrac{{\partial {\textbf{{S}}_S}}}{{\partial {k_3}}}}& {{Y_S} + {Y_H}}& {{X_S} + {X_H}}& 0 \end{array} \end{array} \right]^{{\rm T}}\end{align}
   where
   (31)\begin{equation}\textbf{{L}}=\left[ {\begin{array}{*{20}{c}} {{\textbf{{l}}_1}}\\ {{\textbf{{l}}_2}}\\ {{\textbf{{l}}_3}} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} 1& { - c}& { - b}\\ c& 1& { - a}\\ b& a& 1 \end{array}} \right]\end{equation}

Once the initial parameter ${\textbf{{X}}_0}$ is determined, the partial derivatives above can be used to linearise Equation (13), and then the least square method can be adopted to solve for the optimal $\textbf{{X}}$. The initial values of $({x_0},{y_0},f,{k_1},{k_2},{k_3})$ can be obtained by indoor calibration before launching the rover. In this situation, the vectors ${\textbf{{S}}_S}$ and ${\textbf{{S}}_H}$ can be calculated according to Equations (5) and (6), respectively, and the approximate values of $(a,b,c)$ can be calculated by Equation (13) without any linearisation or iteration. Regardless of the attitude of the rover, the initial rotation parameters $(a,b,c)$ are easy to determine due to the special Rodrigues matrix.

In this paper, the Hipparcos catalogue is utilised to calculate the vector ${\textbf{{S}}_H}$, the precision of which can reach the milliarcsecond level without considering atmospheric refraction (Perryman et al., Reference Perryman, Lindegren, Kovalevsky, Høg, Bastian, Bernacca, Crézée, Donati, Grenon, Grewing, Leeuwen, Marel, Mignard, Murray, Le Poole, Schrijver, Turon, Arenou, Froeschlé and Petersen1997). Since the directional measurement accuracy of the artificial control points has difficulty reaching the milliarcsecond level, theoretically, precise self-calibration for the star sensor can be realised according to the model above.

4.1 Star image processing

Once a star image is captured, the median filter is used to reduce the white noise, the eight connected domain algorithm is used to determine the edge of the star point, and the grey centroid method is used to calculate the centre of the star point. Finally, an improved quadrilateral matching algorithm is utilised to identify the star.

4.2 Star sensor calibration

After identifying the star, abundant observations ${\textbf{{S}}_S}$ and ${\textbf{{S}}_H}$ are obtained, and the self-calibration model proposed in section 3 is applied to solve for the star sensor parameters $({x_0},{y_0},f,{k_1},{k_2},{k_3})$ .

4.3 Tilt and atmospheric refraction correction

In engineering practice, it is almost impossible to adjust the optical axis of the star sensor to be consistent with the plumb line, so a tilt correction is needed to correct the observing vector ${\textbf{{S}}_S}$ according to the output values ${\theta _\textrm{X}}$ and ${\theta _\textrm{Y}}$ of the inclinometer. The corrected altitude $h^{\prime}$ is obtained as follows:

(32)\begin{align}{\textbf{{S}}^{\prime}_S} & = \textbf{{R}}_\textrm{I}^\textrm{H} \cdot {\textbf{{S}}_S}\end{align}

(33)\begin{align}h^{\prime} & = asin ({\textbf{{S}}^{\prime}_S}(3))\end{align}

where

(34)\begin{equation}\textbf{{R}}_\textrm{I}^\textrm{H} = \left[ {\begin{array}{*{20}{c}} {\cos \beta }& 0& { - \sin \beta }\\ {\sin \gamma \sin \beta }& {\cos \gamma }& {\sin \gamma \cos \beta }\\ {\cos \gamma \sin \beta }& { - \sin \gamma }& {\cos \gamma \cos \beta } \end{array}} \right]\end{equation}

and, as in Yang et al. (Reference Yang, Xie and Liu2014),

(35)\begin{equation}\left\{ \begin{array}{l} \gamma ={-} {\theta_\textrm{Y}}\\ \beta = asin \dfrac{{\sin {\theta_\textrm{X}}}}{{\cos \gamma }} \end{array} \right.\end{equation}

In addition, the atmosphere bends light, so a refraction model (related mainly to the temperature, the pressure and the star's altitude) is needed. In this paper, the common model provided by Bangert et al. (Reference Bangert, Puatua, Kaplan, Bartlett and Monet2009) is utilised to calculate the refraction at the Earth's surface. For the Mars rover, a simple refraction model is sufficient because of the extremely thin atmosphere around Mars.

4.4 Quasi-position determination

Due to the inevitable inconsistency between the star sensor frame and the inclinometer frame, it is impossible to correct totally the inclination of the star sensor, which means that the corrected optical axis is still slightly inconsistent with the plumb line. Hence, the position determined by $h^{\prime}$ according to Equation (1) is called the quasi-position $({\lambda _i},{\varphi _i})$ .

4.5 Accurate position determination

Zhan et al. (Reference Zhan, Zheng, Li, Wang, Zhu and Chen2020) numerically proved that the quasi-positions are distributed along a circle on the surface of the planet, while the real position is located at the centre of the circle. Therefore, the spherical circle fitting algorithm is utilised to estimate the accurate position $(\lambda ,\varphi )$ and the tilt correction error $\theta$ as follows (Zhan et al., Reference Zhan, Zheng and Zhang2016):

(36)\begin{equation}\cos (\theta ) = \sin \varphi \sin {\varphi _i} + \cos \varphi \cos {\varphi _i}\cos ({\lambda _i} - \lambda )\end{equation}

Therefore, it is necessary to rotate the platform in different directions and capture several star images to determine the unknown parameters $(\lambda ,\varphi ,\theta )$. In the field test, eight uniform directions were set to fit the accurate position.

4.6 Iterative calculation

After calculating the accurate position, the new star sensor parameter X is used to update the initial parameters, which is beneficial for star image processing; the position coordinates $(\lambda ,\varphi )$ are used to update the initial position to calculate more accurate horizontal coordinates for the star, which may improve the accuracy of the star sensor self-calibration. The iterative calculation terminates when the coordinates $(\lambda ,\varphi )$ change little, for example, by only 0⋅01 arcseconds. In summary, the proposed ACP algorithm can realise celestial positioning autonomously regardless of how many days the rover travels.

5.1 Experimental conditions

The field test data were collected using a prototype with a star sensor and an inclinometer, as shown in Figure 4 (Zhan et al., Reference Zhan, Zheng, Li, Wang, Zhu and Chen2020). The star sensor is composed of a high-resolution CCD and a fish-eye lens. The CCD is an ALT-16M made by Apogee Company with a 4,096-by-4,096 array of 9 μm square pixels, while the fish-eye lens is an AF DX made by Nikon with a 16 mm focal length. Based on the integration of the fish-eye lens with the ALT-16M CCD, the star sensor has a 120° field of view. As mentioned in section 3.1, a star sensor with a super-large field of view was chosen to reduce the weight of the Mars rover through multiplexing technology by adding or removing a median filter between the lens and the CCD. The super-large field-of-view sun sensor easily captures sun images for determining the azimuth without a rotating platform, which is beneficial for reducing the weight of the Mars rover.

Figure 4. Prototype used for data collection

The inclinometer is a Leica Nivel230 with a resolution of 0⋅2′′. The sensor data are processed by an ARK-1122F embedded computer with an Intel Atom dual-core processor. Additionally, the high-accuracy time base is maintained by a Microsemi SA.45s chip-scale atomic clock with a minor aging rate of 3⁻¹⁰ s, which corrects the embedded computer's clock error before performing celestial positioning. The time stamp of each image is recorded by the embedded computer with an error of the order of milliseconds. According to Equation (2), the time observation error is a systematic error of S impacting the astronomical longitude by approximately 0⋅015′′ per millisecond. Since the radius of Mars is only half that of Earth, the systematic positioning error is equivalent to approximately 20 cm, which can be ignored.

To validate the performance of the ACP algorithm, three field tests were conducted on three nights with clear weather or thin clouds in Central China. On the nights of 23 July, 7 November and 13 November 2017, 36, 24 and 24 positioning tests were conducted, respectively. In each positioning test, the prototype was rotated into eight different directions with an angle interval of 45°; that is, that every eight quasi-positions were utilised to fit the accurate position. The real position of the experimental site used to evaluate the accuracy was obtained by a Y/JGT-01 astrogeodetic system based on a high-accuracy image total station. The accuracy of the real position is greater than ±6 m (Zhang, Reference Zhang2009; Shi et al., Reference Shi, Zhang, Du, Li, Ye, Zhang, Chen, Li, Ma, Lin and Mi2020), which satisfies the validation to some extent.

5.2 Self-calibration results

Although hundreds of star images were obtained each night, only one star image is needed to self-calibrate the star sensor. In this paper, the first star image of each positioning test is adopted to calibrate the star sensor. The initial values of the star sensor parameters on 23 July are provided by the indoor calibration, the initial values of the star sensor parameters on 7 November are provided by the self-calibration results from 23 July, and the initial values of the star sensor parameters on 13 November are provided by the self-calibration results from 7 November, thereby simulating the real situation of long-term travel for the rover without indoor calibration. Relevant studies about the indoor calibration can be found in Zhan et al. (Reference Zhan, Chen and He2018b), for which a 10 m dome with 37 simulated stars was used as the control point.

Figure 5 shows the first star image from the field test on 23 July 2017, which contains 1,112 star points that were successfully recognised after image processing. To ensure the reliability of extracting the centroids of the star points, the limiting magnitude of each star is set to 6⋅5. Figure 6 depicts the positions of the star points on the image plane, revealing an almost perfect uniform distribution with over 1,000 star points. For the indoor calibration, a two-axis rotary table and a star simulator are needed, and hundreds of images are captured by adjusting the attitude of the star sensor to obtain a uniform star distribution (Sun et al., Reference Sun, Xing and You2013; Wei et al., Reference Wei, Zhang, Fan, Jiang and Li2014; Xiong et al., Reference Xiong, Wei, Zhang and Jiang2015). Compared with the indoor calibration, the self-calibration obtains thousands of star points by imaging the night sky only one time.

Figure 5. Typical star image captured during the field tests

Figure 6. Distributions of the star points from the first star images during the three field tests

Figure 7 depicts the pixel residuals with respect to the self-calibration of each night. In total, the standard deviations of the residuals in the horizontal and vertical axes are smaller than ±0⋅39 pixels, which indicates the reliability of the self-calibration. However, the distributions of the residuals are not uniform, as the residuals on the vertical axis are larger than those on the horizontal axis. The authors hypothesise that the optical axis is not perpendicular to the CCD plane due to the manual assembly of the CCD camera and lens. In other words, a small angle should exist between the optical axis and the normal direction to the CCD plane, and the optical axis should lean towards the y-axis of the CCD plane, affecting the symmetry of the radial distortion of the imaging model. A more precise imaging model that considers the misalignment between the optical axis and the CCD plane's normal direction will be constructed in the authors’ future research.

Figure 7. Self-calibration residuals on each of three nights

Table 1 describes the final self-calibration results of the three field tests, including the initial values of the star sensor parameters provided by the indoor calibration. Table 1 indicates that ${x_0}$ and ${y_0}$ vary from six to 14 pixels, while f varies by approximately two pixels in the three field tests. The radial distortion parameters ${k_1}$, ${k_2}$ and ${k_3}$ also vary among the three field tests.

Table 1. Self-calibration results of the three field tests

Notably, the parameters of the star sensor possibly change daily. The variations in the parameters are caused mainly by temperature fluctuations, vibrations during exploration and so on. Therefore, self-calibration is necessary to reduce the systematic error of the parameters to realise high-accuracy celestial positioning.

5.3 Positioning results

For comparison, the position is solved with two algorithms. For the first algorithm, all three field tests utilise the star sensor parameters obtained by the indoor calibration to solve the position, which is recorded as the noncalibrated celestial position (‘No calibration’). Then, the ACP algorithm proposed in this paper is used to solve the position. Figure 8 gives the positioning errors of the two algorithms.

Figure 8. Positioning errors of the two algorithms

Furthermore, the standard deviations and the mean values of the positioning results are calculated, as shown in Table 2. When the star sensor is not self-calibrated, the mean positioning error reaches 35⋅8 m, and the standard deviation varies from ±28⋅6 m to ±45⋅7 m in the longitude and latitude directions. When the star sensor is self-calibrated, the mean positioning error is less than 15⋅2 m, and the standard deviation is less than 39⋅7 m. Figure 7 and Table 2 clearly indicate that the ACP algorithm improves both the average error and the standard deviation of the positioning results.

Table 2. Statistics of the positioning results (m)

Note: Ave = average; Std = standard deviation.

Systematic error still exists in the positioning results, however, even when the ACP algorithm is used. There may be two causes of this error. First, as analysed in section 5.2, a small angle should exist between the optical axis and the normal direction to the CCD plane, which might produce the systematic observation error in the star's altitude. Second, the star sensor frame is largely misaligned with the inclinometer frame, introducing a systematic tilt correction error. A deeper investigation will be performed to reduce the systematic errors in the authors’ future research.

It should be noted that the self-calibration and celestial positioning results are obtained on clear nights on Earth, which means the star points are always uniformly distributed on the image plane. However, cloudy weather will definitely reduce the number of star points, and as a consequence, the positioning error will likely increase. To date, land and sea-based celestial navigation has typically been affected by the weather conditions. Fortunately, the density of Mars’ atmosphere is only 1% of that of Earth's, so the Mars rover will operate under better weather conditions for celestial navigation than exist on Earth.