Abstract
Since the American Mars Exploration Rover Opportunity landed on Mars in 2004, it has travelled more than 40 km, and heading-determination technology based on its sun sensor has played an important role in safe driving of the rover. A high-precision heading-determination method will always play a significant role in the rover’s autonomous navigation system, and the precision of the measured 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.
1. 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 10-minute 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]. Long-term 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 et al. and Furgale et al. use a digital sun sensor with 140° field of view and an inclinometer called HMR-3000 for heading, and practical tests on Earth indicate that the precision reaches 1° [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.
2. 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.
Due to the objective conditions, all our field tests are conducted on Earth, so we only discuss the heading-determination 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 and , respectively. represents the rotation matrix from frame E to frame H, which can be written by
According to Figure 1, 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 can be calculated by ephemeris such as DE405, we can obtain the sun vector relative to the local horizontal frame by . Furthermore, we can determine the sun’s predicted azimuth relative to local north according to .
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 or 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 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.
3. Algorithms
3.1. Sun Sensor Calibration Model
3.1.1. 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 10-m 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 high-precision 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 , focal length , radial distortion parameters , rotation parameters , and translation parameters . Figure 4 depicts the geometric or physical meaning of the 12 parameters.
, , and are interior elements relative to the sun sensor; and 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 asand we adopt the polynomial model to describe the radial distortion [15]:The azimuth of the control point in the sun sensor frame is written aswhere and . Hence, the vector of the control point in the sun sensor frame is represented asand the vector of the control point in horizontal frame can be written asWe 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 asThe relationship between and 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 and are written asThen, (12) can be written as 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; , 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.169 pixels 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 , , and are ±0.081 pixels, ±0.097 pixels, and ±0.068 pixels, respectively, which indicate high consistency of the calibration results in 8 directions. However, the mean values of , , and for the 8 directions cannot be adopted as the final results because of the strong correlations between and and between 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.
3.2. 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 high-precision 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, , , and are determined; hence, we use the least-square method to calculate and , which describe the position and attitude of the sun sensor frame relative to the horizontal frame. Then, the rotation matrix from the sun sensor frame to the horizontal frame is calculated by Appendix C, and can be described by 3-step rotations around the , , and axes in order:Rotate an angle of around the axis.Rotate an angle of around the axis.Rotate an angle of around the axis.
is written as
Appendix C gives the expressions for , , and . Then, (19) is expanded asHere, we provide the expressions of and :
When the outputs of the inclinometer are adjusted to zero, we rotate the sun sensor frame by steps and ; 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 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.
3.3. 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:The Sobel operator is used to detect the pixel-level edge points of the sun image.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.According to the subpixel edge points , the center of circle is fitted and the value of the cost function is calculated.
(a) A circular sun image and its edge-detection results
(b) An elliptical sun image and its edge-detection results
where is residual. Suppose the initial values of the unknown parameters are ; thus, (23) is linearized asandThe 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 , 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 , and iterations are needed until the absolute values of are smaller than 0.001 pixels. Then, can be calculated according to (22).According to the subpixel edge points, the ellipse’s center is fitted and the value of the cost function 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 can be calculated by
The error equation of each edge point is represented asEquation (33) is linear; therefore, the coefficient matrix and free item vector are written asSuppose the vector of unknown parameters is , which can be calculated by . Then, and are derived by (32), and we use (31) to calculate .The shape of the sun image is judged by comparing to . If , should be adopted. Otherwise, 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 or , 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.
3.4. 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:Extract the centroid of a sun image using the method mentioned in Section 3.3 and calculate the vector of the sun with respect to the sun sensor frame according to the intrinsic parameters of the sun sensor.Rotate the camera coordinate system by and ; then calculate the vector of the sun in the transition frame by
The azimuth of the sun relative to the X axis of the new coordinate system can be calculated by according to . We then use the Naval Observatory Vector Astronomy Software version 3.0 (NOVAS3.0) and DE405 ephemeris to calculate the sun’s predicted azimuth . The observation epoch in UTC is provided by local clock.Calculate the heading according to 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.
4. Field Tests and Results
In order to test our algorithms and methods for high-precision 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,056-by-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.(4)Estimate the heading-determination precision.
4.2. 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.
4.3. Results
Because the real heading is difficult to obtain, we consider the mean value of the measured heading as a reference for precision estimation. Figure 12 depicts the time series of the heading errors of the three tests.
In the first test, the heading errors fluctuate greatly with an amplitude of 2.5 arcminutes. However, in the second and third tests, the heading errors become quite small, and the amplitude is only 1 arcminute. The standard deviations of the three tests are ±0.97′, ±0.30′, and ±0.28′, respectively. It is evident that the second and third tests achieve more precise heading results. It appears that the elevation of the sun has a great impact on the heading determination. Because the sun elevation has a strong correlation with the sun imaging position, Figure 13 directly shows the trajectories of the sun image centroids of the three tests on the image plane.
On July 20, the trajectory of the sun image centroids was close to the principal point with a mean distance of 220.6 pixels. However, the trajectories on Oct 14 and Nov 11 were 642.8 pixels and 785.8 pixels away from the principal point, respectively. Obviously, radial errors of the sun image centroids have no effect on the sun azimuth measurement, whereas it is linearly influenced by tangential errors. Because the tangential error of the sun image centroid is a small quantity, the error of the sun azimuth measurement can be expressed as follows: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 . Supposing that , , , and , we can calculate the measurement errors of the sun azimuth. The results are , , and , 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.
5. 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. When the outputs of the electronic inclinometer are not strictly adjusted to zero, we need a proper model to modify it. 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. Dust on the sun sensor, resulting from dust storms on Mars, and clouds obscuring the view of the sun may produce irregular 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 centroid-extraction precision.
Appendix
A. Partial Derivatives of Vi with respect to the 12 Parameters
Partial derivatives of and with respect to and arewhere and , , , and are written asPartial derivatives of and with respect to arePartial derivatives of and with respect to () arePartial derivatives of with respect to and arePartial derivatives of with respect to and () arePartial derivatives of with respect to , , and arewhere . 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
C. Zernike 7×7 Models
We use Zernike 7×7 models deduced by Gao et al. as follows [24]:
For each pixel-level edge point obtained by Sobel operator, we use , , and for convolution operations. Then we get three important Zernike moments as , , and , and rotation angle is calculated byand the length from center point to subpixel edge point is obtained by
Hence we can get the coordinates of the subpixel edge point as follows:
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (NSFC) under Grant nos. 41704006 and 41504018 and funded by State Key Laboratory of Geo-Information Engineering under Grant no. SKLGIE2016-Z-2-1.