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

Authors: Yinhu Zhan; Shaojie Chen; Donghan He
Journal: Advances in Astronomy
 (2018)

---

## 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.


---

## Body

## 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. Table1 describes the main frames in this paper and gives their names, notations, and definitions.Table 1
Definition of main frames.
Name
Notation
X axis
Y axis
Z axis
Horizontal frame
H
Local north
Local west
Opposite gravity
Earth-centered inertial frame
E
Vernal equinox
Right-handed
North pole
Earth-centered fixed frame
F
Prime meridian
Right-handed
North pole
Transition frame
T
Rover heading
Right-handed
Opposite gravity
Inclinometer frame
I
X axis of itself
Y axis of itself
Right-handed
Sun sensor frame
C
Horizontal pixels
Right-handed
Optical axisDue 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. Figure1 provides the definitions and space relation of the important frames, and ω is the rover’s heading to be estimated.Figure 1
Space relation of the important frames.We assume that the sun vectors in the E and H frames areS
E and S
H, respectively. R
E
H represents the rotation matrix from frame E to frame H, which can be written by(1)
S
H
=
R
E
H
·
S
E
.According to Figure1, R
E
H can be derived by double rotations:(2)
R
E
H
=
R
Y
φ
-
9
0
∘
·
R
Z
G
+
λ
+
18
0
∘
,where G 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 S
E can be calculated by ephemeris such as DE405, we can obtain the sun vector S
H relative to the local horizontal frame by (
1
). Furthermore, we can determine the sun’s predicted azimuth A
H relative to local north according to S
H.Figure2 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 A
C or A
T according to sun image centroid coordinates, projection model, and distortion model. Then, the rover’s heading ω can be calculated by(3)
ω
=
A
H
-
A
T
.Figure 2
Schematic of sun azimuth measurement.It is evident that the precision ofA
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.


## 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′′. Figure3 shows the dome model and measurement sketch of the control points.Figure 3
10-m diameter dome and measurement sketch of control points.The purpose of sun sensor calibration is to obtain the parameters of the sun sensor, including the image principal point(
x
0
,
y
0
), focal length f, radial distortion parameters (
k
1
,
k
2
,
k
3
), rotation parameters (
a
,
b
,
c
), and translation parameters (
X
0
,
Y
0
,
Z
0
). Figure 4 depicts the geometric or physical meaning of the 12 parameters.Figure 4
Geometric or physical meaning of the 12 sun sensor parameters.(
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) are interior elements relative to the sun sensor; (
a
,
b
,
c
) and (
X
0
,
Y
0
,
Z
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(4)
θ
=
2
arcsin
⁡
r
2
f
+
Δ
θ
.Suppose(
x
,
y
) are coordinates of the image centroid of a control point. The polar distance r is represented as(5)
r
=
x
-
x
0
2
+
y
-
y
0
2
,and we adopt the polynomial model to describe the radial distortion [15]:(6)
Δ
θ
=
k
1
arcsin
⁡
r
2
f
2
+
k
2
arcsin
⁡
r
2
f
3
+
k
3
arcsin
⁡
r
2
f
4
.The azimuth of the control point in the sun sensor frame is written as(7)
A
C
=
atan
⁡
y
p
x
p
x
p
>
0
,
y
p
>
0
atan
⁡
y
p
x
p
+
π
x
p
<
0
atan
⁡
y
p
x
p
+
2
π
x
p
>
0
,
y
p
<
0
,where x
p
=
x
-
x
0 and y
p
=
y
-
y
0. Hence, the vector of the control point in the sun sensor frame is represented as(8)
S
C
=
X
C
Y
C
Z
C
=
sin
⁡
θ
cos
⁡
A
C
sin
⁡
θ
sin
⁡
A
C
cos
⁡
θand the vector of the control point in horizontal frame can be written as(9)
S
H
=
X
H
Y
H
Z
H
T
.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(10)
S
=
X
H
+
X
0
Y
H
+
Y
0
Z
H
+
Z
0
T
.
S is normalized as(11)
S
0
=
S
S
=
X
Y
Z
T
.The relationship between S
C and S
0 can be described by a rotation matrix R, as follows:(12)
S
0
=
R
·
S
C
.Equation (12) is our basic observation equation of sun sensor calibration. Unlike previous studies, we first employ the antisymmetric matrix Q instead of Euler angles to express R, 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]:(13)
R
=
I
-
Q
-
1
I
+
Q
,where I and Q are written as(14)
I
=
1
0
0
0
1
0
0
0
1
Q
=
0
-
c
-
b
c
0
-
a
b
a
0
.Then, (12) can be written as(15)
I
-
Q
·
S
0
=
I
+
Q
·
S
C
.
v
1
v
2
v
3
T represents the error vector of (15), and the error equation is expressed by(16)
V
1
V
2
V
3
=
0
-
Z
C
-
Z
-
Y
C
-
Y
-
Z
C
-
Z
0
X
C
+
X
Y
C
+
Y
X
C
+
X
0
a
b
c
-
X
-
X
C
Y
-
Y
C
Z
-
Z
C
.Equation (16) needs to be linearized before being solved. Appendix A gives the partial derivatives of V
i 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:(17)
v
=
A
·
X
-
l
,where v is the residual vector; A is coefficient matrix; X
=
a
b
c
T, which is the vector of unknown parameters; and l is the free term vector. The least-square method is used to estimate the 12 parameters:(18)
X
=
A
T
A
-
1
A
T
l
.


#### 3.1.2. Calibration Results

Figure5 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.Figure 5
Image of the 37 control points.Figure 6
Residuals of the fish-eye sun sensor 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.Table 2
Results of the 6 interior parameters in 8 directions.
Raw
azimuth
x
0
(pixel)
y
0
(pixel)
f
(pixel)
k
1
k
2
k
3
Std
(
Δ
x)
Std
(
Δ
y)
0°
1502.292
1585.044
855.242
0.212471
-0.740583
0.877803
0.1621
0.1728
45°
1502.209
1584.992
855.142
0.221231
-0.764446
0.897239
0.1604
0.1852
90°
1502.280
1585.034
855.156
0.211819
-0.734883
0.870884
0.1788
0.1553
135°
1502.365
1585.170
855.201
0.181307
-0.653738
0.807112
0.1640
0.1481
180°
1502.293
1585.219
855.047
0.161114
-0.594148
0.755274
0.1454
0.1674
225°
1502.209
1585.248
855.101
0.160012
-0.594534
0.758023
0.1547
0.1603
270°
1502.133
1585.201
855.086
0.191343
-0.678503
0.825247
0.1763
0.1502
315°
1502.141
1585.089
855.214
0.203631
-0.713076
0.853851
0.1809
0.1521In Table2, the standard deviations of x
0, y
0, and f 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 (
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) for the 8 directions cannot be adopted as the final results because of the strong correlations between f and (
k
1
,
k
2
,
k
3
) and between (
x
0
,
y
0
) and (
a
,
b
,
c
). 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 Section3.1, (
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) are determined; hence, we use the least-square method to calculate (
a
,
b
,
c
) and (
X
0
,
Y
0
,
Z
0
), which describe the position and attitude of the sun sensor frame relative to the horizontal frame. Then, the rotation matrix R
C
H from the sun sensor frame to the horizontal frame is calculated by Appendix C, and R
C
H 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.R
C
H is written as(19)
R
C
H
=
R
Z
κ
R
Y
ψ
R
X
γ
.AppendixC gives the expressions for R
Z
(
κ
), R
Y
(
ψ
), and R
X
(
γ
). Then, (19) is expanded as(20)
R
C
H
=
cos
⁡
ψ
cos
⁡
κ
sin
⁡
γ
sin
⁡
ψ
cos
⁡
κ
+
cos
⁡
γ
sin
⁡
κ
-
cos
⁡
γ
sin
⁡
ψ
cos
⁡
κ
+
sin
⁡
γ
sin
⁡
κ
-
cos
⁡
ψ
sin
⁡
κ
-
sin
⁡
γ
sin
⁡
ψ
sin
⁡
κ
+
cos
⁡
γ
cos
⁡
κ
cos
⁡
γ
sin
⁡
ψ
sin
⁡
κ
+
sin
⁡
γ
cos
⁡
κ
sin
⁡
ψ
-
sin
⁡
γ
cos
⁡
ψ
cos
⁡
γ
cos
⁡
ψ
.Here, we provide the expressions of γ and ψ:(21)
γ
=
-
atan
⁡
R
C
H
3,2
R
H
C
3,3
ψ
=
asin
⁡
R
C
H
3,1
.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 Section3.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.Table 3
Values ofγ and ψin 8 directions.
Direction
number
Approximate
azimuth
γ
(′′)
ψ
(′′)
1
0°
1123.9
121.7
2
45°
1120.6
142.1
3
90°
1124.1
139.5
4
135°
1126.6
122.5
5
180°
1129.3
129.2
6
225°
1124.7
135.1
7
270°
1122.1
144.0
8
315°
1125.0
138.4Standard 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.


### 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:(
1
)
The Sobel operator is used to detect the pixel-level edge points(
x
,
y
) of the sun image.
(
2
)
For each pixel-level edge point, the Zernike moment is used to detect the subpixel edge point(
x
′
,
y
′
). 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(
x
′
,
y
′
), the center of circle is fitted and the value of the cost function J
c is calculated.(22)
J
c
x
c
,
y
c
,
r
=
∑
i
=
1
n
x
i
′
-
x
c
2
+
y
i
′
-
y
c
2
-
r
2
2Figure 7
Two classic sun images and their edge-detection results.
(a)
A circular sun image and its edge-detection results
(b)
An elliptical sun image and its edge-detection resultswherev
i is residual. Suppose the initial values of the unknown parameters are X
0
=
x
c
0
y
c
0
r
0
T; thus, (23) is linearized as(23)
v
i
=
a
i
δ
x
c
+
b
i
δ
y
c
+
c
i
δ
r
-
l
iand(24)
a
i
=
2
x
c
0
-
2
x
i
′
b
i
=
2
y
c
0
-
2
y
i
′
c
i
=
r
0
2
-
x
i
′
-
x
c
0
2
-
y
i
′
-
y
c
0
2
l
i
=
r
0
2
-
x
i
′
-
x
c
0
2
-
y
i
′
-
y
c
0
2
.The error equations of all the edge points can be written in the form of a matrix, as follows:(25)
v
=
A
·
δ
X
-
l
,where v is the residual vector, A is the coefficient matrix, δ
X is the correction vector corresponding to X
0, and l is the free term vector. The expanded forms of the four matrices and vectors are(26)
v
=
v
1
v
2
…
v
n
T
,
(27)
A
=
a
1
b
1
c
1
⋮
⋮
⋮
a
n
b
n
c
n
,
(28)
δ
X
=
δ
x
c
δ
y
c
δ
r
T
,
(29)
l
=
l
1
l
2
…
l
n
T
.We assume all the edge points have the same weight, and thenδ
X is calculated by [23](30)
δ
X
=
A
T
A
-
1
A
T
l
.The unknown parameters are updated by X
=
X
0
+
δ
X, and iterations are needed until the absolute values of δ
X are smaller than 0.001 pixels. Then, J
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 functionJ
e is calculated.(31)
J
e
A
,
B
,
C
,
D
,
E
=
∑
i
=
1
n
A
x
i
′
2
+
B
x
i
′
y
i
′
+
C
y
i
′
2
+
D
x
i
′
+
E
y
i
′
+
1
2A, B, C, D, andE are the basic parameters of the ellipse that need to be estimated, and the ellipse’s center (
x
e
,
y
e
) can be calculated by(32)
x
e
=
B
E
-
2
C
D
4
A
C
-
B
2
y
e
=
B
D
-
2
A
E
4
A
C
-
B
2
.The error equation of each edge point is represented as(33)
v
i
=
A
x
i
′
2
+
B
x
i
′
y
i
′
+
C
y
i
′
2
+
D
x
i
′
+
E
y
i
′
+
1
.Equation (33) is linear; therefore, the coefficient matrix A and free item vector l are written as(34)
A
=
x
1
′
2
x
1
′
y
1
′
y
1
′
2
x
1
′
y
1
′
⋮
⋮
⋮
⋮
⋮
x
n
′
2
x
n
′
y
n
′
y
n
′
2
x
n
′
y
n
′
,
(35)
l
=
-
1
-
1
…
-
1
T
.Suppose the vector of unknown parameters is X
=
A
B
C
D
E
T, which can be calculated by X
=
(
A
T
A
)
-
1
A
T
l. Then, x
e and y
e are derived by (32), and we use (31) to calculate J
e.(
5
)
The shape of the sun image is judged by comparingJ
c to J
e. If J
c
<
J
e, (
x
c
,
y
c
) should be adopted. Otherwise, (
x
e
,
y
e
) should be adopted.
(
6
)
The precision is estimated. The root-mean-square error (RMSE) of the centroids can be calculated by(36)
R
M
S
E
=
m
x
2
+
m
y
2
.m
x and m
y are the RMSE of centroid coordinates (
x
c
,
y
c
) or (
x
e
,
y
e
), which can be calculated according to the residual vector v and coefficient matrixA.Figure8 is the flow chart of our centroid-extraction algorithm for the sun image.Figure 8
Flow chart of centroid extraction of 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:(
1
)
Extract the centroid of a sun image using the method mentioned in Section3.3 and calculate the vector S
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 byR
X
(
γ
) and R
Y
(
ψ
); then calculate the vector of the sun in the transition frame S
T by(37)
S
T
=
R
Y
ψ
R
X
γ
S
C
=
X
T
Y
T
Z
T
.The azimuthA
T of the sun relative to the X axis of the new coordinate system can be calculated by S
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 A
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. Figure9 shows the flow chart of our heading-determination algorithm, where the contents of the red boxes are the key algorithms presented in Section 3.Figure 9
Flow chart of heading-determination algorithm.


## 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′′. Figure3 shows the dome model and measurement sketch of the control points.Figure 3
10-m diameter dome and measurement sketch of control points.The purpose of sun sensor calibration is to obtain the parameters of the sun sensor, including the image principal point(
x
0
,
y
0
), focal length f, radial distortion parameters (
k
1
,
k
2
,
k
3
), rotation parameters (
a
,
b
,
c
), and translation parameters (
X
0
,
Y
0
,
Z
0
). Figure 4 depicts the geometric or physical meaning of the 12 parameters.Figure 4
Geometric or physical meaning of the 12 sun sensor parameters.(
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) are interior elements relative to the sun sensor; (
a
,
b
,
c
) and (
X
0
,
Y
0
,
Z
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(4)
θ
=
2
arcsin
⁡
r
2
f
+
Δ
θ
.Suppose(
x
,
y
) are coordinates of the image centroid of a control point. The polar distance r is represented as(5)
r
=
x
-
x
0
2
+
y
-
y
0
2
,and we adopt the polynomial model to describe the radial distortion [15]:(6)
Δ
θ
=
k
1
arcsin
⁡
r
2
f
2
+
k
2
arcsin
⁡
r
2
f
3
+
k
3
arcsin
⁡
r
2
f
4
.The azimuth of the control point in the sun sensor frame is written as(7)
A
C
=
atan
⁡
y
p
x
p
x
p
>
0
,
y
p
>
0
atan
⁡
y
p
x
p
+
π
x
p
<
0
atan
⁡
y
p
x
p
+
2
π
x
p
>
0
,
y
p
<
0
,where x
p
=
x
-
x
0 and y
p
=
y
-
y
0. Hence, the vector of the control point in the sun sensor frame is represented as(8)
S
C
=
X
C
Y
C
Z
C
=
sin
⁡
θ
cos
⁡
A
C
sin
⁡
θ
sin
⁡
A
C
cos
⁡
θand the vector of the control point in horizontal frame can be written as(9)
S
H
=
X
H
Y
H
Z
H
T
.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(10)
S
=
X
H
+
X
0
Y
H
+
Y
0
Z
H
+
Z
0
T
.
S is normalized as(11)
S
0
=
S
S
=
X
Y
Z
T
.The relationship between S
C and S
0 can be described by a rotation matrix R, as follows:(12)
S
0
=
R
·
S
C
.Equation (12) is our basic observation equation of sun sensor calibration. Unlike previous studies, we first employ the antisymmetric matrix Q instead of Euler angles to express R, 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]:(13)
R
=
I
-
Q
-
1
I
+
Q
,where I and Q are written as(14)
I
=
1
0
0
0
1
0
0
0
1
Q
=
0
-
c
-
b
c
0
-
a
b
a
0
.Then, (12) can be written as(15)
I
-
Q
·
S
0
=
I
+
Q
·
S
C
.
v
1
v
2
v
3
T represents the error vector of (15), and the error equation is expressed by(16)
V
1
V
2
V
3
=
0
-
Z
C
-
Z
-
Y
C
-
Y
-
Z
C
-
Z
0
X
C
+
X
Y
C
+
Y
X
C
+
X
0
a
b
c
-
X
-
X
C
Y
-
Y
C
Z
-
Z
C
.Equation (16) needs to be linearized before being solved. Appendix A gives the partial derivatives of V
i 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:(17)
v
=
A
·
X
-
l
,where v is the residual vector; A is coefficient matrix; X
=
a
b
c
T, which is the vector of unknown parameters; and l is the free term vector. The least-square method is used to estimate the 12 parameters:(18)
X
=
A
T
A
-
1
A
T
l
.


### 3.1.2. Calibration Results

Figure5 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.Figure 5
Image of the 37 control points.Figure 6
Residuals of the fish-eye sun sensor 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.Table 2
Results of the 6 interior parameters in 8 directions.
Raw
azimuth
x
0
(pixel)
y
0
(pixel)
f
(pixel)
k
1
k
2
k
3
Std
(
Δ
x)
Std
(
Δ
y)
0°
1502.292
1585.044
855.242
0.212471
-0.740583
0.877803
0.1621
0.1728
45°
1502.209
1584.992
855.142
0.221231
-0.764446
0.897239
0.1604
0.1852
90°
1502.280
1585.034
855.156
0.211819
-0.734883
0.870884
0.1788
0.1553
135°
1502.365
1585.170
855.201
0.181307
-0.653738
0.807112
0.1640
0.1481
180°
1502.293
1585.219
855.047
0.161114
-0.594148
0.755274
0.1454
0.1674
225°
1502.209
1585.248
855.101
0.160012
-0.594534
0.758023
0.1547
0.1603
270°
1502.133
1585.201
855.086
0.191343
-0.678503
0.825247
0.1763
0.1502
315°
1502.141
1585.089
855.214
0.203631
-0.713076
0.853851
0.1809
0.1521In Table2, the standard deviations of x
0, y
0, and f 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 (
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) for the 8 directions cannot be adopted as the final results because of the strong correlations between f and (
k
1
,
k
2
,
k
3
) and between (
x
0
,
y
0
) and (
a
,
b
,
c
). 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.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′′. Figure3 shows the dome model and measurement sketch of the control points.Figure 3
10-m diameter dome and measurement sketch of control points.The purpose of sun sensor calibration is to obtain the parameters of the sun sensor, including the image principal point(
x
0
,
y
0
), focal length f, radial distortion parameters (
k
1
,
k
2
,
k
3
), rotation parameters (
a
,
b
,
c
), and translation parameters (
X
0
,
Y
0
,
Z
0
). Figure 4 depicts the geometric or physical meaning of the 12 parameters.Figure 4
Geometric or physical meaning of the 12 sun sensor parameters.(
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) are interior elements relative to the sun sensor; (
a
,
b
,
c
) and (
X
0
,
Y
0
,
Z
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(4)
θ
=
2
arcsin
⁡
r
2
f
+
Δ
θ
.Suppose(
x
,
y
) are coordinates of the image centroid of a control point. The polar distance r is represented as(5)
r
=
x
-
x
0
2
+
y
-
y
0
2
,and we adopt the polynomial model to describe the radial distortion [15]:(6)
Δ
θ
=
k
1
arcsin
⁡
r
2
f
2
+
k
2
arcsin
⁡
r
2
f
3
+
k
3
arcsin
⁡
r
2
f
4
.The azimuth of the control point in the sun sensor frame is written as(7)
A
C
=
atan
⁡
y
p
x
p
x
p
>
0
,
y
p
>
0
atan
⁡
y
p
x
p
+
π
x
p
<
0
atan
⁡
y
p
x
p
+
2
π
x
p
>
0
,
y
p
<
0
,where x
p
=
x
-
x
0 and y
p
=
y
-
y
0. Hence, the vector of the control point in the sun sensor frame is represented as(8)
S
C
=
X
C
Y
C
Z
C
=
sin
⁡
θ
cos
⁡
A
C
sin
⁡
θ
sin
⁡
A
C
cos
⁡
θand the vector of the control point in horizontal frame can be written as(9)
S
H
=
X
H
Y
H
Z
H
T
.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(10)
S
=
X
H
+
X
0
Y
H
+
Y
0
Z
H
+
Z
0
T
.
S is normalized as(11)
S
0
=
S
S
=
X
Y
Z
T
.The relationship between S
C and S
0 can be described by a rotation matrix R, as follows:(12)
S
0
=
R
·
S
C
.Equation (12) is our basic observation equation of sun sensor calibration. Unlike previous studies, we first employ the antisymmetric matrix Q instead of Euler angles to express R, 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]:(13)
R
=
I
-
Q
-
1
I
+
Q
,where I and Q are written as(14)
I
=
1
0
0
0
1
0
0
0
1
Q
=
0
-
c
-
b
c
0
-
a
b
a
0
.Then, (12) can be written as(15)
I
-
Q
·
S
0
=
I
+
Q
·
S
C
.
v
1
v
2
v
3
T represents the error vector of (15), and the error equation is expressed by(16)
V
1
V
2
V
3
=
0
-
Z
C
-
Z
-
Y
C
-
Y
-
Z
C
-
Z
0
X
C
+
X
Y
C
+
Y
X
C
+
X
0
a
b
c
-
X
-
X
C
Y
-
Y
C
Z
-
Z
C
.Equation (16) needs to be linearized before being solved. Appendix A gives the partial derivatives of V
i 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:(17)
v
=
A
·
X
-
l
,where v is the residual vector; A is coefficient matrix; X
=
a
b
c
T, which is the vector of unknown parameters; and l is the free term vector. The least-square method is used to estimate the 12 parameters:(18)
X
=
A
T
A
-
1
A
T
l
.


## 3.1.2. Calibration Results

Figure5 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.Figure 5
Image of the 37 control points.Figure 6
Residuals of the fish-eye sun sensor 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.Table 2
Results of the 6 interior parameters in 8 directions.
Raw
azimuth
x
0
(pixel)
y
0
(pixel)
f
(pixel)
k
1
k
2
k
3
Std
(
Δ
x)
Std
(
Δ
y)
0°
1502.292
1585.044
855.242
0.212471
-0.740583
0.877803
0.1621
0.1728
45°
1502.209
1584.992
855.142
0.221231
-0.764446
0.897239
0.1604
0.1852
90°
1502.280
1585.034
855.156
0.211819
-0.734883
0.870884
0.1788
0.1553
135°
1502.365
1585.170
855.201
0.181307
-0.653738
0.807112
0.1640
0.1481
180°
1502.293
1585.219
855.047
0.161114
-0.594148
0.755274
0.1454
0.1674
225°
1502.209
1585.248
855.101
0.160012
-0.594534
0.758023
0.1547
0.1603
270°
1502.133
1585.201
855.086
0.191343
-0.678503
0.825247
0.1763
0.1502
315°
1502.141
1585.089
855.214
0.203631
-0.713076
0.853851
0.1809
0.1521In Table2, the standard deviations of x
0, y
0, and f 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 (
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) for the 8 directions cannot be adopted as the final results because of the strong correlations between f and (
k
1
,
k
2
,
k
3
) and between (
x
0
,
y
0
) and (
a
,
b
,
c
). 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 Section3.1, (
x
0
,
y
0
), f, and (
k
1
,
k
2
,
k
3
) are determined; hence, we use the least-square method to calculate (
a
,
b
,
c
) and (
X
0
,
Y
0
,
Z
0
), which describe the position and attitude of the sun sensor frame relative to the horizontal frame. Then, the rotation matrix R
C
H from the sun sensor frame to the horizontal frame is calculated by Appendix C, and R
C
H 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.R
C
H is written as(19)
R
C
H
=
R
Z
κ
R
Y
ψ
R
X
γ
.AppendixC gives the expressions for R
Z
(
κ
), R
Y
(
ψ
), and R
X
(
γ
). Then, (19) is expanded as(20)
R
C
H
=
cos
⁡
ψ
cos
⁡
κ
sin
⁡
γ
sin
⁡
ψ
cos
⁡
κ
+
cos
⁡
γ
sin
⁡
κ
-
cos
⁡
γ
sin
⁡
ψ
cos
⁡
κ
+
sin
⁡
γ
sin
⁡
κ
-
cos
⁡
ψ
sin
⁡
κ
-
sin
⁡
γ
sin
⁡
ψ
sin
⁡
κ
+
cos
⁡
γ
cos
⁡
κ
cos
⁡
γ
sin
⁡
ψ
sin
⁡
κ
+
sin
⁡
γ
cos
⁡
κ
sin
⁡
ψ
-
sin
⁡
γ
cos
⁡
ψ
cos
⁡
γ
cos
⁡
ψ
.Here, we provide the expressions of γ and ψ:(21)
γ
=
-
atan
⁡
R
C
H
3,2
R
H
C
3,3
ψ
=
asin
⁡
R
C
H
3,1
.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 Section3.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.Table 3
Values ofγ and ψin 8 directions.
Direction
number
Approximate
azimuth
γ
(′′)
ψ
(′′)
1
0°
1123.9
121.7
2
45°
1120.6
142.1
3
90°
1124.1
139.5
4
135°
1126.6
122.5
5
180°
1129.3
129.2
6
225°
1124.7
135.1
7
270°
1122.1
144.0
8
315°
1125.0
138.4Standard 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.


## 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:(
1
)
The Sobel operator is used to detect the pixel-level edge points(
x
,
y
) of the sun image.
(
2
)
For each pixel-level edge point, the Zernike moment is used to detect the subpixel edge point(
x
′
,
y
′
). 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(
x
′
,
y
′
), the center of circle is fitted and the value of the cost function J
c is calculated.(22)
J
c
x
c
,
y
c
,
r
=
∑
i
=
1
n
x
i
′
-
x
c
2
+
y
i
′
-
y
c
2
-
r
2
2Figure 7
Two classic sun images and their edge-detection results.
(a)
A circular sun image and its edge-detection results
(b)
An elliptical sun image and its edge-detection resultswherev
i is residual. Suppose the initial values of the unknown parameters are X
0
=
x
c
0
y
c
0
r
0
T; thus, (23) is linearized as(23)
v
i
=
a
i
δ
x
c
+
b
i
δ
y
c
+
c
i
δ
r
-
l
iand(24)
a
i
=
2
x
c
0
-
2
x
i
′
b
i
=
2
y
c
0
-
2
y
i
′
c
i
=
r
0
2
-
x
i
′
-
x
c
0
2
-
y
i
′
-
y
c
0
2
l
i
=
r
0
2
-
x
i
′
-
x
c
0
2
-
y
i
′
-
y
c
0
2
.The error equations of all the edge points can be written in the form of a matrix, as follows:(25)
v
=
A
·
δ
X
-
l
,where v is the residual vector, A is the coefficient matrix, δ
X is the correction vector corresponding to X
0, and l is the free term vector. The expanded forms of the four matrices and vectors are(26)
v
=
v
1
v
2
…
v
n
T
,
(27)
A
=
a
1
b
1
c
1
⋮
⋮
⋮
a
n
b
n
c
n
,
(28)
δ
X
=
δ
x
c
δ
y
c
δ
r
T
,
(29)
l
=
l
1
l
2
…
l
n
T
.We assume all the edge points have the same weight, and thenδ
X is calculated by [23](30)
δ
X
=
A
T
A
-
1
A
T
l
.The unknown parameters are updated by X
=
X
0
+
δ
X, and iterations are needed until the absolute values of δ
X are smaller than 0.001 pixels. Then, J
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 functionJ
e is calculated.(31)
J
e
A
,
B
,
C
,
D
,
E
=
∑
i
=
1
n
A
x
i
′
2
+
B
x
i
′
y
i
′
+
C
y
i
′
2
+
D
x
i
′
+
E
y
i
′
+
1
2A, B, C, D, andE are the basic parameters of the ellipse that need to be estimated, and the ellipse’s center (
x
e
,
y
e
) can be calculated by(32)
x
e
=
B
E
-
2
C
D
4
A
C
-
B
2
y
e
=
B
D
-
2
A
E
4
A
C
-
B
2
.The error equation of each edge point is represented as(33)
v
i
=
A
x
i
′
2
+
B
x
i
′
y
i
′
+
C
y
i
′
2
+
D
x
i
′
+
E
y
i
′
+
1
.Equation (33) is linear; therefore, the coefficient matrix A and free item vector l are written as(34)
A
=
x
1
′
2
x
1
′
y
1
′
y
1
′
2
x
1
′
y
1
′
⋮
⋮
⋮
⋮
⋮
x
n
′
2
x
n
′
y
n
′
y
n
′
2
x
n
′
y
n
′
,
(35)
l
=
-
1
-
1
…
-
1
T
.Suppose the vector of unknown parameters is X
=
A
B
C
D
E
T, which can be calculated by X
=
(
A
T
A
)
-
1
A
T
l. Then, x
e and y
e are derived by (32), and we use (31) to calculate J
e.(
5
)
The shape of the sun image is judged by comparingJ
c to J
e. If J
c
<
J
e, (
x
c
,
y
c
) should be adopted. Otherwise, (
x
e
,
y
e
) should be adopted.
(
6
)
The precision is estimated. The root-mean-square error (RMSE) of the centroids can be calculated by(36)
R
M
S
E
=
m
x
2
+
m
y
2
.m
x and m
y are the RMSE of centroid coordinates (
x
c
,
y
c
) or (
x
e
,
y
e
), which can be calculated according to the residual vector v and coefficient matrixA.Figure8 is the flow chart of our centroid-extraction algorithm for the sun image.Figure 8
Flow chart of centroid extraction of 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:(
1
)
Extract the centroid of a sun image using the method mentioned in Section3.3 and calculate the vector S
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 byR
X
(
γ
) and R
Y
(
ψ
); then calculate the vector of the sun in the transition frame S
T by(37)
S
T
=
R
Y
ψ
R
X
γ
S
C
=
X
T
Y
T
Z
T
.The azimuthA
T of the sun relative to the X axis of the new coordinate system can be calculated by S
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 A
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. Figure9 shows the flow chart of our heading-determination algorithm, where the contents of the red boxes are the key algorithms presented in Section 3.Figure 9
Flow chart of heading-determination algorithm.


## 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. Figure10 shows a photograph of our prototype.Figure 10
Actual picture of the 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 Table4, including the date and time, observation times, sun elevation fluctuation, and weather.Table 4
Test conditions.
Date
Beijing Time
Observations
Sun Elevation
Weather
July 20, 2017
12:05–13:01
252
74.3–75.7°
Thin Clouds
Oct 14, 2017
11:42–13:12
250
44.3–46.7°
Sunny
Nov 11, 2017
11:59–13:13
200
35.6–37.7°
Thin FogThe 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. Figure11 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.Figure 11
Times series of the predicted elevations 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. Figure12 depicts the time series of the heading errors of the three tests.Figure 12
Time series of the heading errors.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, Figure13 directly shows the trajectories of the sun image centroids of the three tests on the image plane.Figure 13
The trajectories of sun image centroids.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:(38)
e
A
C
=
e
c
r
,where e
c is the tangential error of sun image centroid and r 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
Mean square error series of sun image centroids.Figure14 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.046
p
i
x
e
l
s. Supposing that e
c
=
0.046
p
i
x
e
l
s, r
1
=
220.6
p
i
x
e
l
s, r
2
=
642.8
p
i
x
e
l
s, and r
3
=
785.8
p
i
x
e
l
s, we can calculate the measurement errors of the sun azimuth. The results are e
A
1
C
=
0.72
′, e
A
2
C
=
0.25
′, and e
A
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.


## 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. Figure10 shows a photograph of our prototype.Figure 10
Actual picture of the 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 Table4, including the date and time, observation times, sun elevation fluctuation, and weather.Table 4
Test conditions.
Date
Beijing Time
Observations
Sun Elevation
Weather
July 20, 2017
12:05–13:01
252
74.3–75.7°
Thin Clouds
Oct 14, 2017
11:42–13:12
250
44.3–46.7°
Sunny
Nov 11, 2017
11:59–13:13
200
35.6–37.7°
Thin FogThe 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. Figure11 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.Figure 11
Times series of the predicted elevations 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. Figure12 depicts the time series of the heading errors of the three tests.Figure 12
Time series of the heading errors.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, Figure13 directly shows the trajectories of the sun image centroids of the three tests on the image plane.Figure 13
The trajectories of sun image centroids.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:(38)
e
A
C
=
e
c
r
,where e
c is the tangential error of sun image centroid and r 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
Mean square error series of sun image centroids.Figure14 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.046
p
i
x
e
l
s. Supposing that e
c
=
0.046
p
i
x
e
l
s, r
1
=
220.6
p
i
x
e
l
s, r
2
=
642.8
p
i
x
e
l
s, and r
3
=
785.8
p
i
x
e
l
s, we can calculate the measurement errors of the sun azimuth. The results are e
A
1
C
=
0.72
′, e
A
2
C
=
0.25
′, and e
A
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.


## 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.(
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. (
3
) 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.



---

*Source: 1493954-2018-08-01.xml*