High-Precision Heading Determination Based on the Sun for Mars Rover

Since the American Mars Exploration Rover Opportunity landed on Mars in 2004, it has travelled more than 40 km, and headingdetermination technology based on its sun sensor has played an important role in safe driving of the rover. A high-precision heading-determinationmethod will always play a significant role in the rover’s autonomous navigation system, and the precision of themeasured heading strongly affects the navigation results. In order to improve the heading precision to the 1-arcminute level, this paper puts forward a novel calibration algorithm for solving the comparable distortion of large-field sun sensor by introducing an antisymmetric matrix.The sun sensor and inclinometer alignment model are then described in detail to maintain a high-precision horizon datum, and a strict sun image centroid-extraction algorithm combining subpixel edge detection with circle or ellipse fitting is presented. A prototype comprising a sun sensor, electronic inclinometer, and chip-scale atomic clock is developed for testing the algorithms, models, and methods presented in this paper. Three field tests were conducted in different months during 2017. The results show that the precision of the heading determination reaches 0.28–0.97 (1σ) and the centroid error of the sun image and the sun elevation are major factors that affect the heading precision.


Introduction
Spirit and Opportunity, the twin rovers of NASA's Mars Exploration Mission, landed on Mars in 2004.Some discoveries of the rovers, such as the evidence of liquid water on the surface of Mars, have been a source of excitement.Opportunity has been kept in service for 13 years and has travelled over 42.195 km, becoming a veritable marathon champion.High-precision attitude-determination technology of Mars rovers, especially the heading-determination technology, has played an important role in safe driving and the realization of scientific goals.Because there is no satellite navigation system like GPS on Mars, the rovers utilize Pancam cameras as sun sensors for heading, which can restrict the error growth of the IMU and improve the dead-reckoning accuracy [1].When the rover remains static, the sun sensor makes a 10minute tracking observation of the sun, and the quaternion estimator (QUEST) method is used to calculate the attitude and heading, the precision of which reaches 1.5 ∘ .Because the field view of sun sensor on Spirit and Opportunity is only 16 ∘ , a rotation platform with a dial is needed.The rover first rotates the sun sensor to the predicted elevation of the sun and then horizontally scans the sky to find the sun [2].Longterm searching and observation of the sun do not fulfill the real-time navigation requirement for Mars rovers, and the low-precision heading greatly affects the dead reckoning.
Several Jet Propulsion Laboratory (JPL) studies have also used the sun sensor for heading determination for Rocky 7 and FIDO field rover, but the method is similar to that of Spirit and Opportunity, and the precision in the field test is also to within a few degrees [3,4].Deans et al. bundled a fish-eye camera and an inclinometer together, and the precision of the heading determination is superior to 1 ∘ [5].Most published heading-determination methods require long-duration observations only when inclinometer data is unavailable.Enright [6,7].Yang et al. use a large-field-of-view sun sensor for heading determination, and the precision reaches 0.1123 ∘ when observing the sun for 30 minutes, which is the best reported precision so far [8].Illyas et al. present a novel algorithm for micro-planetary rover-heading determination using a low-cost sun sensor, and a large number of experiments show that the heading precision reaches 0.09 ∘ (1), which plays an important role in reducing the accumulated heading error of MEMS sensors [9].In a GPS-denied environment, visual navigation can provide accurate localization [10], but the error grows sharply with the distance travelled.Lambert et al. develop a novel approach incorporating the sun sensor and inclinometer measurements directly into the visual odometry pipeline to reduce the error growth of path estimation, and the resulting localization error is only 1.1% of a 10-km distance travelled [11].
The Mars rover-heading determination method via the inertial navigation system and star sensor has also been thoroughly researched.A novel method based on star sensors and the strap-down inertial navigation system (SINS) is put forward to accurately determine the rover's position and attitude, and the initial position and attitude determination for planetary rovers by INS/Star Sensor integration is researched [12].A high-precision SINS/star sensor deeply integrated navigation scheme is effected by He et al. [13].In recent years, the inertial navigation system, star sensor, and visual navigation system have been integrated to improve the navigation accuracy and reliability [14].However, Mars has some atmosphere, which makes stars invisible in the daytime.Furthermore, the rovers always move during the day and rest at night, which makes the star sensor difficult to apply.
It is evident that heading determination by the sun will always play a significant role in Mars exploration rovers, and the precision of the heading strongly affects the navigation results.Because the rovers move with a maximum speed of 5 cm/s on the surface of Mars and are stationary most of time, this paper only considers the problem of static heading determination.The main focus of this paper is using only one image of the sun to achieve a heading-determination precision on the order of 1 arcminute.This paper is organized as follows.First, we introduce the basic theories of heading determination using the sun, and major factors that affect the precision of the heading are analyzed.Next, the sun sensor calibration model, sun sensor, inclinometer alignment model, and sun image centroid-extraction algorithm are described in detail.After that, our prototype is introduced and three field tests are conducted to verify our models, methods, and algorithms.Finally, we give the basic conclusion according to the results of the field tests.

Basic Theories of Heading Determination Using the Sun
In this section, we introduce the main concepts related to this paper and provide basic theories of heading determination using sun sensor, inclinometer, and local clock.Table 1 describes the main frames in this paper and gives their names, notations, and definitions.
P r im e M e r id ia n Due to the objective conditions, all our field tests are conducted on Earth, so we only discuss the headingdetermination problem in the Earth reference frame.Figure 1 provides the definitions and space relation of the important frames, and  is the rover's heading to be estimated.
We assume that the sun vectors in the E and H frames are  E and  H , respectively. H E represents the rotation matrix from frame E to frame H, which can be written by According to Figure 1,  EH can be derived by double rotations: where  is the Greenwich sidereal time calculated by the existing formula according to observation epoch, and  and  are the longitude and latitude of the rover provided by dead reckoning on Mars.Because  E can be calculated by ephemeris such as DE405, we can obtain the sun vector  H relative to the local horizontal frame by (1).Furthermore, we can determine the sun's predicted azimuth  H relative to local north according to  H . Figure 2 gives the schematic of sun azimuth measurement using a sun sensor.Suppose the sun sensor is adjusted to be horizontal and the sun sensor frame C coincides with the transition frame T. Once a sun image is captured, we can calculate  C or  T according to sun image centroid coordinates, projection model, and distortion model.Then, the rover's heading  can be calculated by It is evident that the precision of  T determines the heading results.In order to improve the measurement precision of the sun azimuth, we must strictly solve the following three problems:  (1) The sun sensor always utilizes large-field-of-view lens to avoid platform rotation, but comparable imaging distortion is introduced.An accurate distortion calibration model is needed.
(2) In practice, the sun sensor is difficult to keep absolutely horizontal, so a suitable and accurate algorithm for inclination correction must be seriously considered.
(3) Because the sun is a disk object, it always occupies a regular area on the image plane, which means a proper image-processing algorithm is needed to calculate the centroid of the sun image.Due to the large field of view, the sun sensor always suffers poor resolution.We must optimize the existing centroid-extraction algorithm to improve the centroid-extraction precision.

Error Equation.
Because there is no regular control point on the surface of Mars, the sun sensor is usually calibrated in the laboratory before launch.We build a 10m diameter dome, on the surface of which we uniformly install 37 super-bright optical fibers as control points.The 3D coordinates of the control points are surveyed by highprecision total station Leica TDRA6000, which has an orientation measurement precision of 0.5  .Figure 3 shows the dome model and measurement sketch of the control points.
The purpose of sun sensor calibration is to obtain the parameters of the sun sensor, including the image principal point ( 0 ,  0 ), focal length , radial distortion parameters ( 1 ,  2 ,  3 ), rotation parameters (, , ), and translation parameters ( 0 ,  0 ,  0 ). Figure 4 depicts the geometric or physical meaning of the 12 parameters.
( 0 ,  0 ), , and ( 1 ,  2 ,  3 ) are interior elements relative to the sun sensor; (, , ) and ( 0 ,  0 ,  0 ) are exterior orientation elements relative to local horizontal frame; and Δ is the error of the half angle of view caused by radial distortion.Because the sun sensor utilizes a solid-angle projection, the real half angle of view is represented as Suppose (, ) are coordinates of the image centroid of a control point.The polar distance  is represented as and we adopt the polynomial model to describe the radial distortion [15]: )) 3 The azimuth of the control point in the sun sensor frame is written as where   =  −  0 and   =  −  0 .Hence, the vector of the control point in the sun sensor frame is represented as and the vector of the control point in horizontal frame can be written as We move the horizontal frame to make its origin coincide with that of the sun sensor frame, and the vector of the control point in the new frame becomes is normalized as The relationship between  C and  0 can be described by a rotation matrix , as follows: Equation ( 12) is our basic observation equation of sun sensor calibration.Unlike previous studies, we first employ the antisymmetric matrix  instead of Euler angles to express , which is beneficial for reducing the amount of calculation through reduced use of trigonometric sines and cosines.It has been proven that all the rotation matrices can be expressed by an antisymmetric matrix as follows [16]: where I and  are written as Then, ( 12) can be written as T represents the error vector of ( 15), and the error equation is expressed by Equation ( 16) needs to be linearized before being solved.Appendix A gives the partial derivatives of   with respect to the 12 parameters, which can be used to linearize (16).Because there are 37 control points, 37×3 linearized error equations can be obtained.We express the error equations in the form of a vector: where  is the residual vector;  is coefficient matrix;  = [  ] T , which is the vector of unknown parameters; and  is the free term vector.The least-square method is used to estimate the 12 parameters: 3.1.2.Calibration Results. Figure 5 is an image of the 37 control points from our sun sensor introduced in Section 4.1.
After threshold operation, we use the gray centroid method to detect the image centroids of the control points [17].Figure 6 depicts the residuals of the 37 image points after calibration.The standard deviations of residuals along x and y axes are ±0.169pixels and ±0.185 pixels, respectively.The calibration precision is similar to the method based on the 2D or 3D calibration board [4,18].As the minimum altitude of control point is about 15 ∘ , the effective calibration field of the sun sensor reaches 150 ∘ .Later, we put the sun sensor near the center of the dome and rotate it in 8 directions in order.In each direction, the sun sensor captures an image of the 37 control points for calibration.The results of the 6 interior parameters of the sun sensor in each direction are listed in Table 2.
In Table 2, the standard deviations of  0 ,  0 , and  are ±0.081pixels, ±0.097 pixels, and ±0.068 pixels, respectively, which indicate high consistency of the calibration results in 8 directions.However, the mean values of ( 0 ,  0 ), , and ( 1 ,  2 ,  3 ) for the 8 directions cannot be adopted as the final results because of the strong correlations between  and ( 1 ,  2 ,  3 ) and between ( 0 ,  0 ) and (, , ).The standard deviation of the residuals, which should be as small as possible, is a good criterion to choose the best set of parameters.Direction 4 achieves the smallest standard deviation of the residuals, so the fourth group of interior parameters should be the final result.

Sun Sensor and Inclinometer Alignment Model.
Prior studies generally use the sun as a control point to determine the relationship between the sun sensor and inclinometer, which must wait for the sun to move across several long traces and obtain low-precision results due to having only one control point in the sky [6,7,19].However, the 10-m diameter dome with 37 control points provides ideal conditions for the sun sensor and inclinometer calibration.Because the coordinates of the 37 control points are surveyed by the highprecision total station, they contain local gravity information, which can be used to determine the relationship between the sun sensor and inclinometer.
In Section 3.1, ( 0 ,  0 ), , and ( 1 ,  2 ,  3 ) are determined; hence, we use the least-square method to calculate (, , ) and ( 0 ,  0 ,  0 ), which describe the position and attitude of the sun sensor frame relative to the horizontal frame.Then, the rotation matrix  H C from the sun sensor frame to the horizontal frame is calculated by Appendix C, and  H C can be described by 3-step rotations around the X C , Y C , and Z C axes in order: (1) Rotate an angle of  around the X C axis.
(2) Rotate an angle of  around the Y C axis.
(3) Rotate an angle of  around the Z C axis.

𝑅 H
C is written as Appendix C gives the expressions for  Z (),  Y (), and  X ().Then, (19)   Here, we provide the expressions of  and : When the outputs of the inclinometer are adjusted to zero, we rotate the sun sensor frame by steps (1) and (2); then the Z axis of the new frame opposes local gravity.In other words, if the outputs of the inclinometer are adjusted to zero and the angles  and  are known, the relationship between sun sensor and inclinometer can be determined.
In the experiment of Section 3.1.2,we always adjust the legs of the prototype introduced in Section 4.1 to keep the outputs of the inclinometer close to zero (smaller than 1  in practice).Table 3 lists the results of  and  in each direction.
Standard deviations of  and  are ±2.6  and ±8.6  , respectively.The mean values ( = 1124.5 ,  = 134.1  ) are adopted as the final parameters.Table 3 indicates the high-precision relationship determination between the sun sensor and inclinometer, which is the basis of high-precision heading determination.

Sun Image Centroid-Extraction Algorithm.
Because large-field sun sensors suffer from poor angular resolution, we must improve the precision of the sun image centroid extraction as much as possible.The gray centroid method is widely used for the sun image centroid extraction, but the precision is not high when the sun image is irregular [17,20].Cui et al. consider the sun image as a circle and use the Sobel operator to detect the pixel-level edge; thus, the centroid is achieved by circle fitting of the edge points [21].Yang et al. adopt the Zernike moment to obtain the subpixel edge points and make the precision of the circular sun image centroid reach 0.07 pixel, which is obviously better than that of the gray centroid algorithm [22].However, for the general projection models, such as the pinhole and solid-angle models, when the sun is away from the boresight direction, the shape of the sun image is more similar to an ellipse than a circle.Our algorithm for sun image centroid extraction is similar to [22], but the difference is that both circle and ellipse sun image are considered, and a reasonable criterion for shape judgment is put forward.Our algorithm is realized by 6 steps: (1) The Sobel operator is used to detect the pixel-level edge points (, ) of the sun image.
(2) For each pixel-level edge point, the Zernike moment is used to detect the subpixel edge point (  ,   ).Appendix C gives the computation method in detail.Figures 7(a) and 7(b) show a practical circular and elliptical sun image, respectively, from our fish-eye sun sensor and its edge-detection results.Obviously, the Zernike moment produces a smoother edge than the Sobel operator, which is beneficial for improving the centroid fitting precision.(3) According to the subpixel edge points (  ,   ), the center of circle is fitted and the value of the cost function  c is calculated.
where V  is residual.Suppose the initial values of the unknown parameters are  0 = [ c0  c0  0 ] T ; thus, ( 23) is linearized as and The error equations of all the edge points can be written in the form of a matrix, as follows: where  is the residual vector,  is the coefficient matrix,  is the correction vector corresponding to  0 , and  is the free term vector.The expanded forms of the four matrices and vectors are We assume all the edge points have the same weight, and then  is calculated by [23] The unknown parameters are updated by  =  0 + , and iterations are needed until the absolute values of  are smaller than 0.001 pixels.Then,  c can be calculated according to (22).
(4) According to the subpixel edge points, the ellipse's center is fitted and the value of the cost function  e is calculated.A, B, C, D, and E are the basic parameters of the ellipse that need to be estimated, and the ellipse's center ( e ,  e ) can be calculated by The error equation of each edge point is represented as Equation ( 33) is linear; therefore, the coefficient matrix  and free item vector  are written as Suppose the vector of unknown parameters is  = [    ] T , which can be calculated by  = ( T ) −  (5) The shape of the sun image is judged by comparing  c to  e .If  c <  e , ( c ,  c ) should be adopted.Otherwise, ( e ,  e ) should be adopted. ( The precision is estimated.The root-mean-square error (RMSE) of the centroids can be calculated by and   are the RMSE of centroid coordinates ( c ,  c ) or ( e ,  e ), which can be calculated according to the residual vector  and coefficient matrix .
Figure 8 is the flow chart of our centroid-extraction algorithm for the sun image.

Heading-Determination Algorithm.
Without loss of generality, we consider the X axis of the prototype to be aligned with that of the charge-coupled device (CCD), and the dip angles of the electronic inclinometer are always adjusted to zero in the field test.Then, we take 4 steps to calculate the heading: (1) Extract the centroid of a sun image using the method mentioned in Section 3.3 and calculate the vector  C of the sun with respect to the sun sensor frame according to the intrinsic parameters of the sun sensor.
(2) Rotate the camera coordinate system by  X () and  Y (); then calculate the vector of the sun in the transition frame  T by The azimuth  T of the sun relative to the X axis of the new coordinate system can be calculated by  T according to (7).We then use the Naval Observatory Vector Astronomy Software version 3.0 (NOVAS3.0)and DE405 ephemeris to calculate the sun's predicted azimuth  H .The observation epoch in UTC is provided by local clock.
(3) Calculate the heading  according to (3)and estimate the precision.
It is clear that only one sun image is needed to calculate the heading in our algorithm, which is beneficial for reducing the data-processing time.Figure 9 shows the flow chart of our heading-determination algorithm, where the contents of the red boxes are the key algorithms presented in Section 3.

Field Tests and Results
In order to test our algorithms and methods for highprecision heading determination, three field tests were conducted in central China's Henan province in 2017.The details of the tests are presented in this section, including the hardware configuration, test conditions, and test results.
4.1.Hardware Configuration.Our sun sensor is composed of a fish-eye lens, a CCD, and a filter.The fish-eye lens is AF DX made by Nikon, with 10.5-mm focal length and 180 ∘ field of view.The CCD is Alta U9000 made by Apogee, with a 3,056by-3,056 array of 12-micron square pixels and a 7 square-inch and 3-inch thick body.The filter is mounted between the lens and CCD to weaken the light of the sun.The sun sensor can capture images of the sun without searching, which means that no rotating platform is needed.
The Leica Nivel230 electronic inclinometer is chosen to obtain the dip angles relative to local horizontal plane, and the precision is up to 1 arcsecond with a measurement range of -3.78 arcminutes to +3.78 arcminutes, a weight of 700 g, and a 3.5-square-inch and 2.7-inch thick body.We choose the SA.45s chip-scale atomic clock made by Microsemi with a month aging rate of 3 −10 s, low power consumption of 120 mW, and volume of 17 cc.The ARK-1122F embedded computer with a dual Atom-core processor is used for processing data including sun images, dip angles, and time information.
We develop a prototype by integrating the sun sensor, inclinometer, and chip-scale atomic clock.Figure 10 shows a photograph of our prototype.
Each test was conducted by the following steps: (1) Put the prototype on a pillar and adjust the legs of the prototype to make the outputs of the inclinometer close to zero.
(2) Keep the prototype static and continuously shoot images of the sun.During shooting, the inclinometer data and time data are collected together.
(3) Use the heading-determination algorithm presented in this paper to calculate the heading.

Test Conditions.
Detailed conditions of the three tests are presented in Table 4, including the date and time, observation times, sun elevation fluctuation, and weather.The three tests were conducted in different months in order to analyze the impact of the sun elevation on the heading determination.On July 20, the maximum elevation of the sun was up to 75.7 ∘ , which means the sun imaging positions were close to the principal point on the image plane.On Oct 14 and Nov 11, the maximum elevation of the sun was 46.7 ∘ and 37.7 ∘ , respectively.This means that the sun imaging positions were away from the principal point on the image plane.Figure 11 shows the time series of the predicted elevation of the sun.Additionally, during the first and third tests, thin clouds or fog sheltered the sensor from the sun, which slightly impacts the observation of the sun.where   is the tangential error of sun image centroid and  is the polar distance of sun image centroid relative to the principal point.The RMSE series of the sun image centroids, which are calculated by fitting the residuals of the subpixel edge points of the sun image, are depicted in Figure 14.
Figure 14 shows that the RMSE of the sun image centroids are similar in the three tests, which indicates that the sun elevation is irrelevant to the sun image centroid extraction.The RMSE of several sun image centroids in the first and third tests are up to 0.1 pixels, probably due to the cloudy or foggy weather.However, the mean RMSE is about 0.065 pixels, which indicates the superiority of our sun image centroid-extraction algorithm.The tangential mean RMSE can be estimated by 0.065/ √ 2 = 0.046pixels.Supposing that   = 0.046 pixels,  1 = 220.6 pixels,  2 = 642.8pixels, and  3 = 785.8pixels, we can calculate the measurement errors of the sun azimuth.The results are  1 C = 0.72  ,  2 C = 0.25  , and  3 C = 0.20  , which essentially coincide with the standard deviations of the heading error series in Figure 12.
Additionally, the heading error series of the three tests are not totally random, and some trend in the variation is  obvious.This is mainly caused by the residual distortion after sun sensor calibration, alignment error of the sun sensor and inclinometer, and periodic change in the weather.As the minimum imaging period of the sun sensor is 10.8 s and the heading calculation time is less than 0.1 s, we conclude that just one sun image can produce a 1-arcminute heading using our prototype.

Conclusions
This paper attempts to improve the heading-determination precision to 1-arcminute level for Mars rovers using only one sun image.Algorithms for the sun sensor calibration, sun sensor and inclinometer alignment, sun image centroid extraction, and heading determination are presented in detail.A prototype is developed, and the results of three ground-based field tests indicate that the precision of heading reaches 0.28-0.97 (1), which is the best reported precision for heading determination in Mars rovers so far.We not only improve the heading precision to the 1-arcminute level, but also reduce the observation time for the sun from 10 to 30 minutes to about 10 seconds.
However, some questions still need to be studied in the future.(1) When the outputs of the electronic inclinometer are not strictly adjusted to zero, we need a proper model to modify it.(2) Because the rover works on Mars for several or tens of years, we want to utilize the stars as control points to calibrate the sun sensor and inclinometer.sun images, which have an impact on the precision of the heading.We are trying to use the robust estimation method to adjust the weight of the edge points to improve the centroidextraction precision.
Partial derivatives of  C and  with respect to   ( = 1, 2, 3) are Partial derivatives of  C with respect to  0 and  0 are where  = √( H +  0 ) 2 + ( H +  0 ) 2 + ( H +  0 ) 2 .According to equations above, it is easy to obtain the partial derivatives of   with respect to 12 fish-eye sun sensor parameters.

B. Basic Rotation Matrix
Rotation matrix around X, Y, and Z axis can be represented by

Figure 1 :
Figure 1: Space relation of the important frames.

Figure 3 :Figure 4 :
Figure 3: 10-m diameter dome and measurement sketch of control points.

Figure 5 :
Figure 5: Image of the 37 control points.

Figure 6 :
Figure 6: Residuals of the fish-eye sun sensor calibration.
A circular sun image and its edge-detection results An elliptical sun image and its edge-detection results

Figure 7 :
Figure 7: Two classic sun images and their edge-detection results.

Figure 8 :
Figure 8: Flow chart of centroid extraction of the sun image.

Figure 11 :
Figure 11: Times series of the predicted elevations of the sun.

Figure 12 :
Figure 12: Time series of the heading errors.

Figure 14 :
Figure 14: Mean square error series of sun image centroids.

Table 1 :
Definition of main frames.

Table 2 :
Results of the 6 interior parameters in 8 directions.

Table 3 :
Values of  and  in 8 directions.
9)Partial derivatives of  C with respect to  and   ( = 1, 2, 3) are Partial derivatives of  0 with respect to  0 ,  0 , and  0 are