Small-angle rotation method for star tracker orientation

Ivan S. Kruzhilov

doi:10.1117/1.JRS.7.073479

19 November 2013 Small-angle rotation method for star tracker orientation

Ivan S. Kruzhilov

Author Affiliations +

Journal of Applied Remote Sensing, Vol. 7, Issue 1, 073479 (November 2013). https://doi.org/10.1117/1.JRS.7.073479

Abstract

The Wahba problem is the task of constrained optimization seeking the matrix from SO(3), which maximally converges (based on the least squares criterion) two sequences of unit vectors. The solution of this task is vital for satellite attitude determination using star trackers. An iterative method for solving the Wahba problem is proposed. Each iteration of the proposed method is reduced to sequential rotation of the vectors and solving the system of linear algebraic equations. Usage of the method implies that the corresponding vectors of both sequences are located sufficiently close to each other. Two variants of the method are proposed, having linear and quadratic convergence. The Wahba problem solution is interpreted in terms of finding the angular velocity of a system of material points, which have certain angular momentum. Taking into consideration the characteristics of state-of-the-art star trackers, one to two iterations are sufficient for finding the optimal solution using the small-angle rotation method. The primary advantage of the proposed method as compared with classical methods based on calculation of eigenvectors and singular decomposition is the simplicity of its implementation.

1. Introduction

Let us assume that matrix $R$ belongs to a special rotation group:

Eq. (1)

R \in SO (3) \Leftrightarrow R \in R^{3 \times 3} : R^{T} R = I, \det (R) = 1,

where

I

is a unit matrix with dimension 3. Let us assume that there are two sequences of unit vectors

v

and

w

in three-dimensional (3-D) space

Eq. (2)

v_{i} \in R^{3} и w_{i} \in R^{3} : | v_{i} | = 1, | w_{i} | = 1,

where

i = 1, \dots, n

, and

n

is the number of vectors in the sequence.

We shall denote the function (3) as loss function

Eq. (3)

L (R) = \frac{1}{2} \sum_{i = 1}^{n} k_{i} {| w_{i} - R v_{i} |}^{2},

where

k_{i}

are weight coefficients, such that

\sum_{i = 1}^{n} k_{i} = 1

.

The norm of vector $x$ hereinafter is understood to be the Euclidean norm

Eq. (4)

| x | = \sqrt{x^{T} x} .

It is necessary to find the matrix $R_{opt}$ from SO(3) which minimizes the loss function

Eq. (5)

R_{opt} : L (R_{opt}) = \underset{R \in SO (3)}{Min} L (R) .

This task is called the Wahba problem in honor of Grace Wahba, who has first posed it in Ref. 1. The geometrical meaning of the formulated task is as follows: finding such a rotation in 3-D space, which will maximally bring together the two systems of vectors.

The Wahba problem has important applied significance because its solution is necessary in stellar orientation systems for satellite attitude determination. Stellar orientation systems are intended for determination of satellite attitude in the inertial (geocentric) reference system, and presently they are the most precise among existing² satellite orientation systems.

General functional diagram of a star tracker is shown in Fig. 1. Light from stars is passing through the optical path of the device and is projected on a photodetector. Using special algorithms,³^,⁴ stars are segregated on the background of photodetector noise and optical distortions, and their directional cosines are determined in the star tracker reference system.

Fig. 1

Functional diagram of a star tracker. 1—stellar sky, 2—optical system, 3—electronic module with photodetector array, 4—computing module, and 5—star tracker.

If vectors $w$ are interpreted as the directional cosines of stars in the inertial geocentric stellar coordinate system and vectors $v$ are interpreted as the directional cosines of the same stars but in the star tracker reference system, then matrix $R_{opt}$ obtained as the result of solving task (5) will determine the attitude of the satellite with reference to the Earth (if the universal time is known). Therefore, hereinafter we shall name $R_{opt}$ as optimal orientation matrix. Coefficients $k_{i}$ are chosen depending on the accuracy of determination of star coordinates on the photosensitive matrix, which is usually proportional to the star energy.

The values of vectors $v$ are known not accurately because of star orientation instrument errors in the course of determining star positions. Such errors are caused by optical perversions (lens distortion and aberration), noise, and discrete structure of photosensitive matrix. The values of vectors $w$ are also determined with errors due to inaccuracy of existing star catalogues.

Let $v^{true}$ and $w^{true}$ be the true (without introduced error) directional cosines of stars in the inertial reference system and the star tracker reference system, and let $R_{true}$ be the matrix determining the rotation, which aligns vectors $v^{true}$ and $w^{true}$

Eq. (6)

w_{i}^{true} = R_{true} v_{i}^{true} .

It is possible to consider that vectors $v$ and $w$ are

Eq. (7)

v_{i} = v_{i}^{true} + ε_{i}^{v} w_{i} = w_{i}^{true} + ε_{i}^{w},

where

ε_{i}^{v} и ε_{i}^{w} \sim N (0, σ)

. If

| ε^{v} |

and

| ε^{w} |

are small, then it holds that

R_{opt} \approx R_{true}

. In practice, the error of star determination by instrument is substantially larger than the star catalogue error,

| ε^{v} | ≫ | ε^{w} |

.

There are two approaches to solve Wahba problem (5): direct minimization of loss function and exact analytical solution. Since the first formulation of the Wahba problem, numerous methods for its analytical solution have been proposed.⁵ Many of these methods seem quite simple, but actual computations come up against complications, e.g., finding the square root from an ill-conditioned matrix.

On the other hand, not much attention in the literature has been devoted to the possibility of direct minimization of the function $L (R)$ . The reason is that the application of minimization methods requires the existence of sufficiently close initial approximation and calculation of function $L (R)$ values at every step. The laboriousness of calculating the values of function $L (R)$ grows with increase of the number of vectors $n$ .

When solving the Wahba problem in connection with the task of a satellite attitude control, it is necessary to take into consideration the special features of the application domain and the resultant constraints on the character of input data. Usage of these constraints allows the creation of more efficient methods for solving the Wahba problem as compared with the generalized problem formulation. Two such constraints on the values of input data should be noted.

First, the number $n$ of stars for determining orientation usually is not over 50 because of optical characteristics of star trackers and distribution of stars over the celestial sphere. Actually, not all the stars are used for recognition, but only 10 to 15 of the brightest stars, because their coordinates have the highest accuracy.

Second, when considering the solution of the Wahba problem in connection with the task of satellite attitude control, one should take into account the fact that the errors of vectors $v$ and $w$ are small as compared with the minimal distance between two arbitrary vectors of the sequence.

Eq. (8)

arcos (v_{i} v_{j}) ≫ | ε_{i}^{v} | ≫ | ε_{i}^{w} | arcos (w_{i} w_{j}) ≫ | ε_{i}^{v} | ≫ | ε_{i}^{w} | \forall i = 1, \dots, n j = 1, \dots, n, i \neq j .

Inequalities in Eq. (8) are fulfilled since star catalogues are compiled providing a condition that a distance between the nearest stars by a few times exceeds permissible error $ε_{i}^{v}$ . If the star position error significantly exceeds the possible one (due to proton influence, pixel defects, etc.), then such stars will not be detected by identification algorithm, since they are rejected according to the criterion of mutual angular distances.

Provided that Eq. (8) is fulfilled, it is possible to calculate easily and with sufficient accuracy the initial approximation for the matrix $R_{0} \approx R_{opt}$ using the triangle attitude determination (TRIAD) method.⁶^,⁷ Using the initial approximation $R_{0}$ , it is possible to rotate vectors $v$ so that they are sufficiently close to vectors $w$ :

Eq. (9)

{\tilde{v}}_{i} = R_{0} v_{i} .

The result from this is that

Eq. (10)

arcos (w_{i} {\tilde{v}}_{i}) \approx | ε_{i}^{v} | .

Hence, when solving the Wahba problem for the task of satellite orientation, it may be assumed that Eq. (10) applies.

2. Review of Methods for Solution of the Wahba Problem

A short but recent review of methods for solving the Wahba problem is given in Ref. 8. An earlier but detailed analysis of methods for solving the Wahba problem was made by Markley and Mortari in Ref. 5. The basic results related to this problem were obtained in 1980s and 1990s;⁹^–¹³ however, due to rapid progress in space instruments development, in particular, in star sensors, this topic is still actively discussed during the last decade.⁵^–⁷^,¹⁴^–²¹

The majority of classic methods for solution of the Wahba problem are based on the usage of matrix $B$ (sometimes called the attitude profile matrix⁶), which is formed on the basis of vectors $v$ and $w$ :

Eq. (11)

B = Σ_{i = 1}^{n} k_{i} w_{i} v_{i}^{T} .

One of the first analytical solutions of the Wahba problem has been proposed by Farrell and Stuelpnagel,²² who have shown that $R_{opt}$ coincides with orthogonal matrix $W$ of polar decomposition $B = WS$ , where $S$ is a symmetrical matrix.

According to Farrell, matrix $W$ may be expressed using matrix $B$ as

Eq. (12)

∥ R_{opt} - {B ∥}_{F} = \min_{A \in SO (3)} ∥ A - {B ∥}_{F},

where

{‖ A ‖}_{F} = tr (A^{T} A)

is the Frobenius norm.

According to Ref. 22, the $W$ matrix may be expressed through the $B$ matrix as

Eq. (13)

R_{opt} = W = B {(B^{T} B)}^{- 0.5} .

Though Eq. (13) appears quite simple, its practical application is connected with considerable difficulties. During solution of the task of satellite orientation, the vectors of sequence $v$ corresponding to the directional vectors of the stars are located sufficiently close to each other (usually within the cone of 10 angular degrees or smaller), due to the size of the star tracker’s field of view (FOV). For this reason, matrix $B^{T} B$ is often an ill-conditioned matrix, and the process of computing its square root, for instance, with the Denman-Beaver iteration²³ is either diverging or has large error.

Matrix $W$ of polar decomposition may be expressed in terms of singular value decomposition (SVD) of matrix

Eq. (14)

B = {UDV}^{T}

as

Eq. (15)

W = {UV}^{T} .

Equation (14) is the basis of the method proposed by Markley.¹¹ The fast optimal attitude matrix (FOAM) and slower optimal matrix algorithm (SOMA) methods by Markley¹⁰ are also based on the idea of SVD of matrix $B$ , which in the mentioned methods is reduced to solving the system of nonlinear algebraic equations. According to the computational experiments performed by Markley,¹⁰ determination of the optimal orientation matrix using the SVD method has higher accuracy than the FOAM and the SOMA methods; however, the latter two methods are faster.

The most popular group of methods for solving Wahba problem are the methods (usually called q-method or Davenport method) calculating the quaternion of the respective orientation matrix $R_{opt}$ , as the eigenvector corresponding to maximal eigenvalue $λ_{\max}$ of matrix $K$ :

Eq. (16)

K = (\begin{matrix} S - s I & b \\ b^{T} & s \end{matrix}),

where

Eq. (17)

s = Tr (B),

Eq. (18)

S = B + B^{T},

Eq. (19)

b = \sum_{i = 1}^{n} k_{i} [w_{i}, v_{i}] .

Vector $b$ may also be expressed using components of matrix $B$

Eq. (20)

b = [\begin{matrix} B_{23} - B_{32} \\ B_{31} - B_{13} \\ B_{12} - B_{21} \end{matrix}] .

Matrix $K$ is a symmetrical matrix; therefore, efficient algorithms exist for calculation of eigenvectors of the matrix, e.g., Jacobi method or the methods based on QR decomposition with preliminary reduction to tridiagonal form.²⁴ However, implementation and debugging of such algorithms may be quite labor consuming, considering that coding often has to be performed using low-level programming languages or erasable programmable logic device (EPLD). It is necessary to keep in mind that the specified task has to be solved in real time. Besides, the capacity of spaceborne processors is significantly constrained. Therefore, despite small dimension of matrix $K$ , calculation of its eigenvectors may be laborious.

To avoid straightforward calculation of eigenvectors of matrix $K$ , Shuster and Oh⁹ have proposed the quaternion estimator (QUEST) algorithm, which gives an explicit formula for finding the eigenvector corresponding to the maximal eigenvalue $λ_{\max}$ , provided that this value is known. The maximal eigenvalue is searched for as the root of a nonlinear algebraic equation, which the authors propose to solve using Newton-Raphson method. However, such method of finding eigenvalues is the least stable method, therefore, direct search of eigenvectors of matrix $K$ gives more accurate results.

Other important results obtained by Schuster are as follows: he was the first to consider the problem from the statistical point of view and has shown that the loss function as a random value is distributed by chi-square distribution law with $2 n - 3$ degrees of freedom:

Eq. (21)

L (R_{opt}) \sim χ_{2} (2 n - 3) .

Schuster also has shown that the covariance matrix $P$ of rotational angle error equals to:

Eq. (22)

P = {(I - B)}^{- 1} .

Attitude angle errors are the values $ϕ = (ϕ_{1}, ϕ_{2}, ϕ_{3})$ , such that

Eq. (23)

R_{opt} = \exp (- [ϕ \times]) R_{true},

where

[ϕ \times]

is a skew-symmetric matrix determined by elements of vector

ϕ

.

Markley¹¹ has shown that in case of a large number of observations, the following approximate estimate of the covariance matrix will be valid

Eq. (24)

P = \frac{σ_{0}^{2}}{ζ} (κ I + {BB}^{T}),

where

κ = (1 / 2) (λ^{2} - B^{2})

,

ζ = κ λ - \det (B)

, and

λ = tr (R_{opt} B^{T})

.

The results obtained by Schuster and Markley show that the presented estimates of the Wahba problem are consistent. However, the issue of unbiasedness of the proposed estimates remains open. In order to determine what is the unbiased estimate for SO(3), at first it is necessary to clarify the concepts of mathematical expectation and the average value for the elements of this manifold. Two approaches for determination of the concept of average value for the rotation group—Euclidean and Riemannian—are given, for example, in Refs. 25 and 26.

The approaches to solve the Wahba problem that are closest to the approach proposed in this work were presented by Park and Kim,²⁷ Mortari et al.,²⁸ and Mortari.¹³ Authors of Ref. 27 were solving the task of constrained minimization similar to Wahba problem for position fixing using global positioning system. The methodology proposed by them is using the exponential relation between Lie group SO(3) and Lie algebra $so (3)$ for determining the gradient of the loss function and for sequential search of its minimum using Newton’s method and steepest descent method.

The “energy approach algorithm”¹³ proposed by Mortari has a common feature with the proposed, in this article, small-angle rotation (SAR) method: based on the closeness of vectors $v$ and $w$ , it substitutes the loss function by an approximate function—an energy function. Despite the similarity of the basic approach, the energy approach algorithm is principally different from the SAR because it is based on a computation of eigenvalues or SVD, while the SAR algorithm is based on solving the system of algebraic equations.

3. SAR Method for Wahba Problem Solution

3.1.

Basic Definitions

Let $so (3)$ be the set of skew-symmetric matrices

Eq. (25)

so (3) = {A \in R^{3 \times 3} : A + A^{T} = 0} .

The exponent of a square matrix $A$ is the square matrix equaling the sum of infinite series

Eq. (26)

\exp (A) = \sum_{k = 0}^{\infty} \frac{A^{k}}{k!} .

The properties of the matrix exponent are verified easily

Eq. (27)

\exp (A^{T}) = \exp {(A)}^{T},

Eq. (28)

\exp (A) \exp (- A) = I,

Eq. (29)

\det [\exp (A)] > 0 .

In that case, if matrix $M \in so (3)$ , then according to Eqs. (27) and (28)

Eq. (30)

\exp {(M)}^{T} \exp (M) = \exp (- M) \exp (M) = I .

Taking into consideration Eq. (29), this means that $\exp (M) \in SO (3)$ . That is, the exponent of a skew-symmetric matrix is a rotation matrix

Eq. (31)

\forall M \in so (3) \Rightarrow \exp (M) \in SO (3) .

This remarkable fact has deep connection with the theory of Lie groups and algebra.²⁹

For vector $x = (x_{1}, x_{2}, x_{3})$ , we shall denote, hereinafter, as $[x \times]$ the skew-symmetric matrix:

Eq. (32)

[x \times] = (\begin{matrix} 0 & x_{3} & - x_{2} \\ - x_{3} & 0 & x_{1} \\ x_{2} & - x_{1} & 0 \end{matrix}) .

Let $[a, b]$ be the operation of cross product, then the following relations hold true

Eq. (33)

[a, b] = [a \times] b = - [b \times] a .

3.2.

SAR Algorithm of the First Order (Angular Momentum Method)

Let us assume that a system of vectors $v$ and $w$ exists, satisfying condition (10) specified in Sec. 1. As already mentioned in Sec. 1, after determining the approximate attitude using TRIAD method,⁶^,⁷ using transformations (9), it is possible to achieve the fulfillment of condition (10).

Let $ω = (ω_{1}, ω_{2}, ω_{3})$ be the rotation vector of the coordinate system. Then by reason of expression (31), it is possible to replace $L (R) = L [\exp (ω)] = L (ω)$ , and in this case, the function $L$ may be represented in the form

Eq. (34)

L (ω) = \sum_{i = 1}^{n} k_{i} {| w_{i} - \exp ([ω \times]) v_{i} |}^{2} .

As mentioned in Sec. 1, the rotation angle $ω_{opt}$

Eq. (35)

ω_{opt} : R_{opt} = \exp ([ω_{opt} \times])

is small. Therefore, it is possible to limit the expansion into series (26) by the first several terms.

Eq. (36)

L (ω) \approx f_{1} (ω) = 0.5 \sum_{i = 1}^{n} k_{i} {| w_{i} - (I + [ω \times]) v_{i} |}^{2};

in this case

Eq. (37)

f_{1} (ω) = const + \sum_{i = 1}^{n} k_{i} (- w_{i}^{T} [ω \times] v_{i} + 0.5 v_{i}^{T} {[ω \times]}^{2} v_{i}) .

Taking into account Eq. (33)

Eq. (38)

f_{1} (ω) = const + \sum_{i = 1}^{n} k_{i} ({[w_{i}, v_{i}]}^{T} ω + 0.5 ω^{T} {[v \times]}^{2} ω),

where const means a constant independent from

ω

. Thus, the task of constrained minimization of the loss function

L (ω)

may be reduced to the task of unconstrained minimization of the quadratic form

f_{1} (ω)

.

The necessary condition for achieving the extreme point by a function is equality of its gradient to zero. Due to symmetric property of matrix ${[v \times]}^{2}$ , gradient $f_{1} (ω)$ equals to

Eq. (39)

\nabla f_{1} (ω) = \sum_{i = 1}^{n} k_{i} {[w_{i}, v_{i}]}^{T} + (\sum_{i = 1}^{n} k_{i} {[v_{i} \times]}^{2}) ω .

The search for $ω$ satisfying the equality $\nabla f_{1} (ω) = 0$ is reduced to solve the system of algebraic equations

Eq. (40)

A ω = b,

where

Eq. (41)

A = \sum_{i = 1}^{n} k_{i} {[v_{i} \times]}^{2}

and

b

is determined in Eq. (19).

It is easy to note that matrix $A$ coincides with the inertia tensor of material points determined by vectors $v_{i}$ and having masses $k_{i}$ . Hence, the solution of Eq. (40) has a vivid physical interpretation—it is necessary to find such angular velocity $ω$ for which the inertia momentum of a rigid system of material points defined by vectors $v$ and having masses $k_{i}$ will equal to the inertia momentum of the same system, if unit velocity $w_{i}$ is applied to each material point. At that, $f_{1} (ω)$ function may be regarded as the value of kinetic energy of the system of material points.

Considering the mentioned mechanical interpretation, the SAR method of the first order may be also called the “angular momentum method.”

It is worth mentioning that Mortari in Ref. 13 also approximated the loss function by a function that he regarded as the kinetic energy function-energy of $n$ compressed springs with coefficients of rigidity $k_{i}$ . In all the rest, Mortari’s approach has nothing common with the solution proposed here.

3.3.

SAR Algorithm of the Second Order

It is logical to expect from a solution of the Wahba problem that in case of permutation of the sequences of vectors $v$ and $w$ , the sought-for rotation matrix $A$ would have to be inversed. For the SAR, it means that the SAR vector $ω$ changes its sign. However, for the SAR of the first order, in case of permutation of $v$ and $w$ , the SAR vector, in general case, changes not only its sign but also the absolute value as well. The reason is that matrix $A$ (41) depends only upon vectors $v$ but does not depend upon vectors $w$ .

A modification of the SAR method invariant to permutation of vectors and giving the more accurate estimate of $ω$ than Eq. (40) is presented in this section.

Let us transform the loss function $L$ into the gain function, as was done in Ref. 6:

Eq. (42)

g (ω) = 1 - L (ω) .

Thus, the task of minimizing the loss function was transformed into the task of maximizing the function $g (ω)$ . Using uncomplicated transformations, it is possible to demonstrate that

Eq. (43)

g (ω) = \sum_{i = 1}^{n} k_{i} v_{i}^{T} \exp ([ω \times]) w_{i} .

Replacing the exponent by the first three terms of relevant series, it is possible to obtain the approximation of the function $g (ω) \approx f_{2} (ω)$

Eq. (44)

g (ω) \approx f_{2} (ω) = \sum_{i = 1}^{n} k_{i} v_{i}^{T} (i + [ω \times] + \frac{{[ω \times]}^{2}}{2}) w_{i}

or

Eq. (45)

f_{2} (ω) = const + (\sum_{i = 1}^{n} k_{i} {[w_{i}, v_{i}]}^{T}) ω + 0.5 ω^{T} (\sum_{i = 1}^{n} k_{i} [v_{i} \times] [w_{i} \times]) ω .

In that case

Eq. (46)

\nabla f_{2} (ω) = \sum_{i = 1}^{n} k_{i} {[w_{i}, v_{i}]}^{T} + 0.5 (\sum_{i = 1}^{n} k_{i} ([v_{i} \times] [w_{i} \times] + [{w_{i} \times]}^{T} {[v_{i} \times]}^{T})) ω .

Finding $ω$ satisfying the equality $\nabla f_{2} (ω) = 0$ is reduced to solving the system of algebraic equations

Eq. (47)

A ω = b,

where

Eq. (48)

A = 0.5 \sum_{i = 1}^{n} k_{i} ([v_{i} \times] [w_{i} \times] + {[w_{i} \times]}^{T} {[v_{i} \times]}^{T}) .

Let us assume $v_{i} = (v_{i}^{1}, v_{i}^{2}, v_{i}^{3})$ and $w_{i} = (w_{i}^{1}, w_{i}^{2}, w_{i}^{3})$ , then matrix $A$ looks as follows:

Eq. (49)

A = \frac{1}{2} \sum_{i = 1}^{n} k_{i} [\begin{matrix} - 2 (v_{i}^{3} w_{i}^{3} + v_{i}^{2} w_{i}^{2}) & v_{i}^{1} w_{i}^{2} + v_{i}^{2} w_{i}^{1} & v_{i}^{1} w_{i}^{3} + v_{i}^{3} w_{i}^{1} \\ v_{i}^{1} w_{i}^{2} + v_{i}^{2} w_{i}^{1} & - 2 (v_{i}^{3} w_{i}^{3} + v_{i}^{1} w_{i}^{1}) & v_{i}^{2} w_{i}^{3} + v_{i}^{3} w_{i}^{2} \\ v_{i}^{1} w_{i}^{3} + v_{i}^{3} w_{i}^{1} & v_{i}^{2} w_{i}^{3} + v_{i}^{3} w_{i}^{2} & - 2 (v_{i}^{2} w_{i}^{2} + v_{i}^{1} w_{i}^{1}) \end{matrix}] .

Matrix $A$ is connected with attitude profile matrix $B$ by the following equation:

Eq. (50)

A = 0.5 S - s I,

where scalar

s

and matrix

S

are expressed through attitude profile matrix

B

, as shown in Eqs. (17) and (18).

3.4.

Matrix Corresponds to the Small-Rotation Vector

When the SAR vector $ω$ is found, the sought-for approximation of matrix $R_{opt}$ is obtained as

Eq. (51)

R_{opt} = \exp ([ω \times]) .

There are many methods for calculating matrix exponent,³⁰ but for the case of a skew-symmetric matrix, the simplest way is to use the Rodriguez formula:²⁹

Eq. (52)

\exp (ω) = I + \sin (Θ) [μ \times] + [1 - \cos (Θ)] {[μ \times]}^{2},

where

Θ = | ω |

,

μ = ω / Θ

.

Instead of using the matrix form for setting the rotation in $R^{3}$ , it is preferable to use the quaternion presentation, which is more convenient as it requires less memory space and allows reduction of the number of operations necessary for vector rotation. Calculation of the orientation quaternion corresponding to rotation $ω$ is easier than calculation of the corresponding matrix:

Eq. (53)

Q = [\cos (Θ), \sin (Θ) μ] .

Though usage of quaternions is more convenient from the computational point of view than usage of matrices, we shall continue using the matrix notations in this work for the sake of convenience of theoretical calculations. If required, the matrix notations may be easily substituted by quaternion analogs.

4. Evaluation of Solution Accuracy

It is necessary to determine the criteria for evaluation of the accuracy of the obtained solution. Different approaches to this issue are possible.

In technical applications, the error of star trackers is measured most often by three Euler angles by which matrix $R_{1}$ must be additionally rotated in order to be aligned with matrix $R_{2}$ . Euler angles set the sequential rotation around axes OX, OY, OZ and are commonly used in mechanics and computer graphics and are called roll, pitch, and yaw.

For star sensors, axis OZ traditionally coincides with the optical axis of the device; therefore, the largest orientation errors are achieved for the third angle yaw (rotation around line of sight), whereas the roll and pitch errors are usually equaling to each other.

Euler angles are convenient for practical usage. Usually, a satellite is always facing the Earth with one of its sides (where scanning equipment or radio transmitter is located). The location of the star tracker within the satellite is known; therefore, the orientation error measured in Euler angles determines the error of directing the satellite on the Earth.

However, Euler angles are extremely inconvenient for theoretical evaluations, as they are not invariant to the change of coordinate system. Besides, it is desirable to express the combined error by one number rather than by three numbers.

For evaluation of a difference of rotation matrices, usually two metrics are used-Riemannian and Euclidean. In the first case, the set of rotation matrices SO(3) is regarded as Riemannian metric manifold with metric²⁵

Eq. (54)

d_{R} (R_{1}, R_{2}) = \frac{1}{\sqrt{2}} {‖ \log (R_{1}^{T} R_{2}) ‖}_{F} .

The Euclidean metric is based on the Frobenius norm

Eq. (55)

d_{E} (R_{1}, R_{2}) = \frac{1}{\sqrt{2}} {‖ R_{1} - R_{2} ‖}_{F} .

The value $d_{R}$ varies between 0 and $π$ , whereas $d_{E}$ varies between 0 and 2. The set of all rotation matrices $SO (3)$ is homeomorphic to a 3-D sphere in four-dimensional space, and each matrix may be interpreted as a point on this sphere. In such an interpretation, $d_{R}$ may be regarded as angle, whereas $d_{E}$ may be regarded as chord between respective points on the sphere.

The Riemannian distance $d (R_{1}, R_{2})$ defines the minimal absolute value of the angle by which the coordinate system $R_{1}$ must be rotated around an arbitrary axis in order to align it with the coordinate system $R_{2}$ .

If the rotation is performed around one and the same axis for the angles $φ_{1}$ and $φ_{2}$ , then for the matrices $R_{1} (φ_{1})$ and $R_{2} (φ_{2})$ , the value of Riemannian metric $d_{R} (R_{1}, R_{2})$ equals to $| φ_{1} - φ_{2} |$ , naturally complying with our perception of orientation error.

Considering that accuracy in optics and in astronomy is usually measured by angular minutes (or angular seconds), the value $d$ determining accuracy will be expressed hereinafter in angular minutes.

It is worth mentioning that in case of close matrices, Euclidean metric and Riemannian metric are close to each other

Eq. (56)

d_{E} (R_{1}, R_{2}) \approx d_{R} (R_{1}, R_{2}) .

Based on geometrical interpretation of metrics, relation (56) means equality of arc length and chord length for small angles. As star trackers are high-precision instruments with accuracy of the order of angular seconds, it may be assumed that relation (56) holds for them. In order to highlight this fact, hereinafter, we shall denote the metric $d_{R}$ simply as $d$ .

5. Step-by-Step Form of the Algorithm

Input data for the algorithm are:

$R_{0}$ —initial approximation of orientation matrix;
$w$ , $v$ —sequences of unit 3-D vectors of dimensionality $n$ ;
$ε$ —required accuracy of solution (in radians).

Output data:

$R$ —optimal orientation matrix.

Intermediate data:

$R_{ω}$ —orientation matrix.

All matrices and vectors have the dimension equal to 3.

The general scheme of the SAR is as follows. Gain function $g (ω)$ is approximated by quadratic form $f_{2} (ω)$ . By solving the system of linear Eq. (42), the minimum of function $f_{2} (ω)$ is found. When the respective rotation angles $ω$ are found which minimize the function $f_{2} (ω)$ , the rotation matrix $R$ is calculated on its basis. Vectors $v$ are rotated by matrix $R$ and approach to $w$ .

These operations are repeated until rotation angle $ω$ becomes negligible in the framework of the task being solved.

The steps of iterative algorithm for solving Wahba problem using the SAR method:

1. Determining the initial approximation $R_{0}$ using TRIAD method.
2. $R_{ω} = R = R_{0}$ .
3. $v_{i} = R_{ω} v_{i}$ , $i = 1, \dots, n$ .
4. Determining, on the basis of formula (42), the parameters of the system of linear algebraic equations $A$ and $b$ using data $w$ and $v$ .
5. Solving the system of linear algebraic equations $A ω = b$ .
6. According to formula (52), obtaining rotation matrix $R_{ω}$ from vector $ω$ .
7. $R = R_{ω} R$ .
8. If the condition ${‖ R_{ω} ‖}_{F}^{2} < ε^{2}$ is satisfied, then go to step 9, else to step 3.
9. End of algorithm.

If the rotation is set by quaternion, the algorithm will be similar.

All the steps of the SAR algorithm are sufficiently simple for coding including low-level programming languages. The most source-consuming among these steps are extraction of square root and solving system of linear equations. Binary arithmetic has fast and effective methods for calculating square root. Implementation of Gauss method is also rather simple even using low-level programming languages. Moreover, solution of 3-D system may be easily implemented directly using Cramer formula (for such dimensions, it has a negligible influence on execution time and accuracy). At the same time, finding eigen numbers, even by means of the simplest Jacobi method, requires nontrivial logic which is not simple for implementation using assembler or EPLD without saying about more complicated methods.²⁴ Methods of singular matrix decomposition are also complicated for implementation.

The first-order SAR algorithm was implemented in software of 329K star tracker,³¹ which successfully operates on board the data relay satellite Luch-5A since December 2011 and on board the satellite Luch-5B since November 2012. Implementation of the algorithm using digital signal processor (DSP) has shown that it is simple for coding and program debugging.

6. Simulation Results

6.1.

Convergence of the Method

The SAR of the first order and the second order is obtained using approximation of matrix exponent by its decomposition into exponential series up to the first and the second power, respectively [Eqs. (36) and (44)]. Thus, approximation by means of loss functions $f_{1} (ω)$ and $f_{2} (ω)$ has, respectively, an error of orders $ω$ and $ω^{2}$ . Since minimum of functions $f_{1} (ω)$ and $f_{2} (ω)$ is seeking exactly at each step according to formulae (40) and (47), so the final error of function $L (ω)$ minimum will also have order of $ω$ and $ω^{2}$ .

In order to test the accuracy and convergence of the proposed method, its modeling in MATLAB system has been performed. Data of Tables 1 Table 2 Table 3–4 are obtained by means of averaging not less than 100,000 experiments.

Table 1

Average error of estimate Ropt using the singular decomposition method and the small-angle rotation (SAR) method of the first order in angular minutes.

Method	Step zero	First step	Second step	Third step	Fourth step	Fifth step
Singular value decomposition (SVD), $Δ_{1}$	26.2657745
SAR, $Δ_{2}$	99.3452572	26.2961055	26.2657045	26.2657761	26.2657744	26.2657745
Difference, $Δ$	73.1	0.0303	$- 6.99 \times 10^{- 5}$	$1.65 \times 10^{- 6}$	$3.83 \times 10^{- 8}$	$3.77 \times 10^{- 10}$

Table 2

Value of loss function L(R) for the estimate using the singular decomposition method and the SAR method of the first order in angular minutes.

Method	Step zero	First step	Second step	Third step	Fourth step	Fifth step
SVD, $Δ_{1}$	2.9383489
SAR, $Δ_{2}$	4.523	2.9395197	2.9383489	2.9383489	2.9383489	2.9383489
Difference, $Δ$	1.58	0.00117	$2.37 \times 10^{- 8}$	$1.22 \times 10^{- 11}$	$7.99 \times 10^{- 15}$	$- 2.6 \times 10^{- 15}$

Table 3

Average error of estimate Ropt using the singular decomposition method and the SAR method of the second order in angular minutes.

Method	Step zero	First step	Second step	Third step	Fourth step
SVD, $Δ_{1}$	26.0686937
SAR, $Δ_{2}$	98.7480572	26.0750189	26.0686937	26.0686937	26.0686937
Difference, $Δ$	72.68	0.0063252	$- 5.5 \times 10^{- 10}$	$1.83 \times 10^{- 12}$	$1.78 \times 10^{- 12}$

Table 4

Value of loss function L(R) for the estimate using the singular decomposition method and the SAR method of the second order in angular minutes.

Method	Step zero	First step	Second step	Third step	Fourth step
SVD, $Δ_{1}$	2.9331031
SAR, $Δ_{2}$	4.5330924	2.9331278	2.9331031	2.9331031	2.9331031
Difference, $Δ$	1.6	$2.5 \times 10^{- 5}$	$2.2 \times 10^{- 15}$	$- 4.4 \times 10^{- 16}$	$- 4.4 \times 10^{- 16}$

Simulation was performed as follows. Rotation matrix $R_{true}$ was set in random manner. $n$ unit vectors $w_{i}$ were uniformly placed in a random manner in the FOV (spherical square $2 W$ ) of the star sensor. Some errors were added in coordinates of each vector to $x$ and $y$ components of each vector. Vectors $v_{i} = R_{true} w_{i} + ε$ were determined with subsequent normalization $v_{i} = v_{i} / | v_{i} |$ , where $ε \sim N (0, σ)$ —normally distributed random values. In this case, it was not allowed that star angular distances to be less than 10 SD $σ$ of the added error. Such restriction is quite reasonable since the adjacent stars are not included usually at the guide star catalogues.

Using the SAR and singular decomposition method, estimates of the matrix $R_{opt}$ were obtained for the sequence $w$ and $v$ : $R_{opt 1} и R_{opt 2}$ . According to formula (49), the errors $Δ_{1} = d (R_{opt 1}, R_{true})$ and $Δ_{2} = d (R_{opt 2}, R_{true})$ were calculated.

The mentioned sequence of operations was performed for two SAR methods: the first and second order.

The values $Δ_{1}$ and $Δ_{2}$ and their difference $Δ = Δ_{2} - Δ_{1}$ for the SAR methods of the first and second order are given in Tables 1 and 3. The values of square root from the loss function $L$ are given in Tables 2 and 4.

The following parameters were used during modeling: the side of the square FOV $2 W = 20$ angular degrees, number of stars $n = 15$ , SD $σ = 10.0$ arc min. It should be noted that the SD value used during the simulation is much larger than actual values. This was done in order to increase the number of steps executed by the algorithm before reaching the optimal point. Solution by TRIAD method was used as the initial approximation.

As may be seen from Tables 1 Table 2 Table 3–4, the singular decomposition method provides more accurate estimate than the SAR though the difference becomes negligibly small after several iterations. Since the exact value $R_{opt}$ is not known, the estimate $R_{opt}$ may be approximated (with the accuracy of singular decomposition by method³²) by the estimate $R_{opt}$ of the SVD method. In this case, the third row of the table may be interpreted as the error of the SAR method.

In order to better understand the convergence behavior of the SAR, Figs. 2 and 3 show the dependence of the difference between accuracies of the methods from the number of iterations (in logarithmic scale). Every line in Figs. 2 and 3 corresponds to a single experiment (differing from Tables 1 Table 2 Table 3–4 whose data are obtained by means of averaging).

Fig. 2

Dependence of the difference between accuracies of the small-angle rotation (SAR) method of the first order and the singular value decomposition (SVD) method $Δ = | Δ_{1} - Δ_{2} |$ . OX axis, number of iterations of SAR method. OY axis, difference between accuracies of the two methods in angular minutes (logarithmic scale).

Fig. 3

Dependence of the difference between accuracies of the SAR method of the second order and the SVD method $Δ = | Δ_{1} - Δ_{2} |$ . OX axis, number of iterations of SAR method. OY axis, difference between accuracies of the two methods in angular minutes (logarithmic scale).

The simulation results confirm the conclusion about the method’s convergence rate. As may be seen from Fig. 2 and Table 1, SAR of the first order has linear convergence. Based on Fig. 3 and Table 3, the conclusion may be made that SAR of the second order has super-linear convergence. Graphs in Fig. 3 have form which differs from parabola because the difference between the SAR and the SVD method is stabilized starting from the third step. This is caused by the error of number representation with the type “double” and the error of the singular decomposition algorithm implemented by LAPACK software package.³²

At any rate, experimental data are indicative of the fact that the second-order method is by far more effective than that of the first order; therefore, the former is more expedient to be used for calculations. Given this, its computational complexity is just very little higher than that of the first order method. $A$ matrix calculation for the second-order method will require more multiplication operations than it will for the first-order method. Otherwise, computational complexity of the algorithms is identical.

Simulation shows that if the error of initial approximation of orientation is of the order of 1 angular minute, then after one iteration of the SAR, the achieved resulting accuracy is of the order of $10^{- 5}$ angular minute, abundantly meeting the requirements to the errors of state-of-the-art star trackers. In case of necessity, two iterations of the method may be performed. Performing more than three iterations of the SAR makes no sense because in this case, the accuracy of the method exceeds the accuracy of data presentation using the type double.

6.2.

Method Stability against Single Errors

An important issue in star tracker attitude estimation is the presence of inaccurate data. Radiation, pixel failures, etc., can cause one or more of the star vectors to have noise values that are significantly larger than others. Methods such as the SVD method and Davenport’s q-method are very robust to this, while faster methods such as QUEST have more problems dealing with this.

Simulation was performed to verify stability of the SAR against errors of such kind. The errors of 60 arc sec, i.e., significantly larger than usual errors for stars, were added to directing vectors of two stars from 15. Since SAR accuracy significantly depends on accuracy of initial approximation, the most interesting for simulation is the case when two vectors having maximum errors are chosen to estimate initial approximation. The simulation data are presented in Table 5. The data of Table 5 were obtained by means of averaging not less than 50,000 experiments.

Table 5

Comparison of SVD and SAR algorithms’ stability against single errors (errors of SD=60 arc sec are added to the first two vectors).

$σ$ (arc sec)	1	5	10	15
SVD (arc sec)	55.66742	57.10473	61.02592	67.56212
SAR, 1 step	55.69210	57.10473	61.02592	67.56212
SAR1–SVD	$1.3 \times 10^{- 6}$	$- 1.9 \times 10^{- 9}$	$- 3.5 \times 10^{- 9}$	$8.0 \times 10^{- 10}$
SAR, 2 steps	55.66742	57.10473	61.02592	67.56212
SAR2–SVD	$- 1.4 \times 10^{- 10}$	$- 1.9 \times 10^{- 9}$	$- 3.6 \times 10^{- 9}$	$8.7 \times 10^{- 10}$

Taking into account the characteristics of modern star trackers (with the FOV of angular radius $W \approx 10 \deg$ ), SD for star position error is not able to exceed 1 arc min since else it will not be identified by the star identification algorithm according to criterion of mutual angular distances. Proceeding from this estimation, errors of 60 arc sec were added to the first two stars’ positions at the simulation. The SD $σ$ of errors for other stars is presented in Table 5 along with the simulation results.

The given Table 5 shows that the SAR algorithm is stable against random single errors and possesses, concerning this criterion, properties close to SVD algorithm. The difference in accuracy for SVD and SAR methods is comparable with the data presentation error.

7. Computational Aspects

7.1.

Time of Algorithm’s Execution

The purpose of this section is to compare the computational complications of the SAR algorithm and the classical algorithms of solving the Wahba problem based either on SVD or on calculation of eigenvectors.

The components requiring the most computational efforts at each iteration of the SAR are rotation of $n$ vectors with subsequent computation of matrix $A$ and evaluation of rotation vector $ω$ . From the computational point of view, these tasks are reduced to sequential multiplication of $n$ vectors by square matrix or quaternions and subsequent solution of the system of linear algebraic equations.

Although finding the eigenvectors is similar to solving the system of linear equations and is formally characterized by asymptotic³³ complexity $O (m^{3})$ , the actual number of operations required for finding the eigenvectors is much higher. The reason is that the Jacobi method, similar to the methods of matrix reduction to standard forms and QR decomposition,²⁴ requires numerous calculations of square root. Besides, the methods of finding the eigenvectors are iterative, i.e., require numerous recalculations until the desired accuracy is achieved. On the other hand, solving a system of linear equations requires a finite number of steps.

According to Ref. 33, the following statement is valid.

Statement:

The Gauss method for a system of equations with dimensionality $m$ requires

Eq. (57)

\frac{m}{3} (m^{2} + 3 m - 1),

multiplication and division operations.

Based on Eq. (57), solution of the system of equations with dimensionality $3 \times 3$ using the Gauss method requires only 17 multiplication operations.

It is obvious that computational complexity of all the algorithms mentioned in the review is $O (n)$ depending on the number of stars, since these algorithms are related to calculating attitude profile matrix $B$ [Eq. (11)]. Moreover, besides the calculation of the matrix $B$ rotation of $n$ vectors is required for SAR algorithm. Due to additional computational consumption for the vectors’ rotation, SAR algorithm has linear coefficient, which is larger than SVD and q-method have.

It is also obvious that time of the algorithm operation linearly depends on quantity of iterations.

To compare algorithmic complexity of the proposed algorithms with state-of-the-art algorithms, evaluation of their execution time in PC was performed. SVD, q-method, QUEST, and SAR of the second order for one and two iterations were chosen for the comparison. The algorithms were implemented using MATLAB 7/0 and a computer with the processor Intel Core Duo 2.00 GHz. Time of the algorithms’ execution versus number of stars and the number of iterations are presented in Figs. 4 and 5.

Fig. 4

Time of the second order SAR algorithms’ execution versus number of stars. The digits denote characteristics of the following algorithms: 1—quaternion estimator (QUEST), 2—q-method, 3—SVD, 4—SAR one iteration, and 5—SAR two iterations.

Fig. 5

Time of the second order SAR algorithm’s execution versus number of iterations. The digits denote characteristics of the following algorithms: 1—SAR, 2—SVD, 3—q-method (Davenport), and 4—QUEST.

Interpretive program MATLAB is known for its slow operation, so the absolute values of execution time should be regarded very accurately. Real time of algorithm’s execution depends significantly on the method of its implementation using one or another processor. For comparison, time of SAR algorithm’s execution at a single iteration for 10 to 15 stars using 40-MHz DSP processor is equal to 10 μs approximately. For modern space class DSP, this time will be tens of milliseconds.

It is seen from Fig. 4 that the methods SVD, q-method, and QUEST have the same linear coefficient of increase depending on the number of stars. The SAR method has larger linear coefficient comparing with SVD, the q-method, and QUEST and therefore its execution time will rise a greater extent at increasing number of stars. However, as was mentioned above, not all stars are used for recognition but usually only 10 to 15 of the brightest stars because their coordinates have the highest accuracy.

The fact, that SAR algorithms’ execution time rises a greater extent at increasing number of stars comparing with SVD, q-method, and QUEST is a fundamental property of the algorithms. At the same time, the constant constituent of the algorithm’s execution time (the graph value for two stars) may significantly vary depending on implementation of the algorithm using one or another processor.

The SAR algorithm’s execution time versus quantity of the method iterations is presented in Fig. 5. The number of stars in FOV is 10. Execution time for the algorithms SVD, the q-method, and QUEST (shown by horizontal lines) is given for comparison. As should be expected, the dependence has linear form. For quantity of iterations more than two, time of algorithms’ execution for SAR is essentially larger than for SVD, q-method, and QUEST. However, as was mentioned in the preceding section, to meet current accuracy requirements, the method’s one or two iterations are quite sufficient.

7.2.

Conditionality of the Equation System Solution

Unlike the mentioned methods of finding the eigenvectors, solving a system of linear equations using the Gauss method does not require the calculation of square root, and that is the advantage of the SAR. However, employing other method of solution than the Gauss method may necessitate calculation of square root. For example, since matrix $A$ is a symmetric matrix, the search for system may be solved using the Cholesky method, which, in general case, is more stable than the Gauss method.

To compare stability of the Gauss method and the Cholesky method, a computational experiment was performed. Both methods were programmed in C++ language using Visual Studio 2010. Numbers were represented with double precision. For device FOV angular radius ranging from 1 to 10 deg, matrices $A$ were generated, and SAR vectors $ω$ were calculated by the Gauss method and the Cholesky method. The computational experiment showed that the difference between solutions obtained by the two methods is within $10^{- 12}$ arc min, which is substantially smaller then the accuracy requirements raised to the present-day star trackers. Therefore, application of the Gauss method instead of the Cholesky method in the context of the task being solved is quite reasonable.

The largest difference in the accuracy of the Gauss method and the Cholesky method was observed for the smallest FOV. The reason is that the closer the stars (and the corresponding vectors) are located to each other, the more degenerated the task of finding the minimum of $L (ω)$ function. The following fact holds that the smaller the FOV of the star tracker, the more elongated function $L$ is.

The problem of degeneracy of the search task is illustrated in Fig. 6 showing the dependence of the loss function $L$ from parameters $ω_{2}$ and $ω_{3}$ . In order to enable the presentation of the dependence in 3-D space, it is supposed that $ω_{1}$ value is known precisely. Figure 6 distinctly shows that diagram $L$ is elongated in shape, and its level lines are similar to ellipses. In optimization theory,³⁴ such functions are named as “ravine functions.”

Fig. 6

Dependence of loss function $L$ from parameters $ω_{2}$ and $ω_{3}$ .

If $L$ function is approximated by a quadratic form, as was done in Eqs. (36) and (44), then the ratio of maximum and minimum singular values of matrix $A$ of the quadratic form may serve as the measure of elongation of $L$ . If matrix $A$ is positively defined, then the square root of the above mentioned ratio is the ratio of the lengths of maximum and minimum semiaxes of the ellipsoid determined by the corresponding quadratic form. The ratio of maximum and minimum singular values is the matrix condition number $cond (A)$ for Euclidean norm. Therefore, the size of the device FOV is connected with the elongation of $L$ function, which, in its turn, is connected with the condition number of matrix $A$ and, respectively, with accuracy of finding $ω$ .

Based on the mathematic simulation, a following observation was made, in case of disposition of stars in a circular FOV of angular radius $W$ , the condition number of matrix $A$ is proportional to squared cotangent of the FOV angular radius $W$ of the star tracker.

The author did not find the proof of this fact, but its truth is confirmed by simulation results for the several particular cases. During simulation, the stars were arranged in the FOV of the device as shown in Fig. 7. In total, 3 possible relative positions of stars were considered, for the number of stars $n = 2, 3, 4$ . Simulation results are shown in Fig. 8. The indicated dependence is valid also for other quantity and mutual positions of the stars.

Fig. 7

Relative disposition of stars in field of view (FOV) of angular radius $W$ .

Fig. 8

Dependence of the condition number of matrix $A$ from the FOV’s angle of the star tracker. The dots denote the condition numbers computed experimentally, the line denotes the approximation by the function $k^{*} \tan {(W)}^{- 2}$ . The numbers 1, 2, and 3 correspond to three different dispositions of stars according to Fig. 7.

It is worth mentioning that the problem of degradation of conditionality when the FOV is decreasing is not only the specific problem of the SAR but is also typical for the Wahba problem in general.

8. Conclusions

The SAR method proposed in this work is an iterative sequential method for solving the Wahba problem, having linear and quadratic convergence. For the data represented by double precision, the optimal solution is achieved after two iterations, with an error of initial approximation many times higher than actually possible. For actual data, one or two iterations of the method are sufficient. Each iteration of the method is reduced to sequential multiplication of $n$ vectors by a matrix or to multiplication of quaternions and solving of the system of linear algebraic equations with dimensionality 3. The method is stable against single errors of individual vectors.

The primary advantage of the proposed method as compared with classical methods based on calculation of eigenvectors and singular decomposition is the simplicity of implementation.

The SAR method of the first order was implemented in software of wide-angle star tracker 329K. The trackers of this type successfully operate already more than 1 year on board of the Russian data relay geostationary satellites Luch-5A and Luch-5B.

9. Future Works

It would be useful to determine the theoretical dependence of conditionality of the task being solved from size of the FOV, which, as experimentally shown in this work, is proportional to the squared cotangent of the FOV’s angle of the star tracker.

Acknowledgments

This paper is written in the context of developing at JSC «SPE «Geofizika-Cosmos» the algorithms for multiple-head APS-based star tracker 348K. I would like to thank Andrey Zakharov of the Sternberg Astronomical Institute (Moscow State University) who encouraged me to conduct this research and Olga Amosova of Moscow Power Engineering Institute for her advice. Besides, I would like to thank Antonina Ivannikova for her support.

References

1.

G. Wahba, “A least-squares estimate of satellite attitude,” SIAM Rev., 7 (3), 409 (1965). http://dx.doi.org/10.1137/1007077 SIREAD 0036-1445 Google Scholar

2.

K. M. Huffman, “Designing star trackers to meet micro-satellite requirements,” Massachusetts Institute of Technology, (2006). Google Scholar

3.

W. ZhangW. QuanL. Guo, “Blurred star image processing for star sensors under dynamic conditions,” Sensors, 12 (5), 6712 –6726 (2012). http://dx.doi.org/10.3390/s120506712 SNSRES 0746-9462 Google Scholar

4.

I. Kruzhilov, “Minimization of point light source coordinates determination error on photo detectors,” J. Opt. Commun., 32 (4), 201 –204 (2011). http://dx.doi.org/10.1515/JOC.2011.037 JOCODG 0173-4911 Google Scholar

5.

F. L. MarkleyD. Mortari, “New developments in quaternion estimation from vector observations,” Adv. Astronaut. Sci., 106 373 –393 (2000). ADASA9 0065-3438 Google Scholar

6.

J. A. Tappe, “Development of star tracker system for accurate estimation of spacecraft attitude,” Naval Postgraduate School, (2009). Google Scholar

7.

S. TanyginM. D. Shuster, “The many triad algorithms,” Adv. Astronaut. Sci., 127 81 –99 (2007). ADASA9 0065-3438 Google Scholar

8.

Z. Menget al., “A brief survey of the deterministic solution for satellite attitude estimation,” Proc. SPIE, 7129 712928 (2008). http://dx.doi.org/10.1117/12.807456 PSISDG 0277-786X Google Scholar

9.

M. D. ShusterS. D. Oh, “Three-axis attitude determination from vector observations,” J. Guid. Control, 4 (1), 70 –77 (1981). http://dx.doi.org/10.2514/3.19717 JGCODS 0162-3192 Google Scholar

10.

F. L. Markley, “Attitude determination using vector observations: a fast optimal matrix algorithm,” J. Astronaut. Sci., 41 (2), 261 –280 (1993). JALSA6 0021-9142 Google Scholar

11.

F. L. Markley, “Attitude determination using vector observations and the singular value decomposition,” J. Astronaut. Sci., 36 (3), 245 –258 (1988). JALSA6 0021-9142 Google Scholar

12.

M. D. Shuster, “A survey of attitude representations,” J. Astronaut. Sci., 41 (4), 439 –517 (1993). JALSA6 0021-9142 Google Scholar

13.

D. Mortari, “Energy approach algorithm for attitude determination from Vector Observations,” J. Astronaut. Sci., 45 (1), 41 –55 (1997). JALSA6 0021-9142 Google Scholar

14.

F. L. Markley, “Fast quaternion attitude estimation from two vector measurements,” J. Guid. Control Dyn., 25 (2), 411 –414 (2002). http://dx.doi.org/10.2514/2.4897 JGCODS 0731-5090 Google Scholar

15.

Y. ChengM. D. Shuster, “Robustness and accuracy of the QUEST algorithm,” Adv. Astronaut. Sci., 127 41 –61 (2007). ADASA9 0065-3438 Google Scholar

16.

F. L. MarkleyM. Mortari, “Quaternion attitude estimation using vector measurements,” J. Astronaut. Sci., 48 (2), 359 –380 (2000). JALSA6 0021-9142 Google Scholar

17.

D. Mortari, “Second estimator for the optimal quaternion,” J. Guid. Control Dyn., 23 (5), 885 –888 (2000). http://dx.doi.org/10.2514/2.4618 JGCODS 0731-5090 Google Scholar

18.

T. DelabieJ. De SchutterB. Vandenbussche, “Highly efficient attitude-estimation algorithm for star trackers using optimal image matching,” J. Guid. Control Dynam., 36 (6), 1672 –1680 (2013). http://dx.doi.org/10.2514/1.61082 Google Scholar

19.

M. D. Shuster, “The QUEST for better attitudes,” J. Astronaut. Sci., 54 (3), 657 –683 (2006). http://dx.doi.org/10.1007/BF03256511 JALSA6 0021-9142 Google Scholar

20.

D. Choukrounet al., “Optimal-REQUEST algorithm for attitude determination,” in Proc. AIAA Guid. Navig. Control Conf., (2001). Google Scholar

21.

J. C. HinksM. L. Psiaki, “Solution strategies for an extension of Wahba’s problem to a spinning spacecraft,” J. Guid. Control Dyn., 34 (6), 1734 –1745 (2011). http://dx.doi.org/10.2514/1.53530 JGCODS 0731-5090 Google Scholar

22.

J. L. FarrellJ. C. Stuelpnagel, “A least squares estimate of spacecraft attitude,” SIAM Rev., 8 (3), 384 –386 (1966). http://dx.doi.org/10.1137/1008080 SIREAD 0036-1445 Google Scholar

23.

E. D. DenmanA. N. Beavers, “The matrix sign function and computations in systems,” Appl. Math. Comput., 2 (1), 63 –94 (1976). http://dx.doi.org/10.1016/0096-3003(76)90020-5 AMHCBQ 0096-3003 Google Scholar

24.

D. S. Watkins, “The QR algorithm revisited,” SIAM Rev., 50 (1), 133 –145 (2008). http://dx.doi.org/10.1137/060659454 SIREAD 0036-1445 Google Scholar

25.

M. Moakher, “Means and averaging in the group of rotations,” SIAM J. Matrix Anal. Appl., 24 (1), 1 –16 (2002). http://dx.doi.org/10.1137/S0895479801383877 SJMAEL 0895-4798 Google Scholar

26.

I. SharfA. WolfM. B. Rubin, “Arithmetic and geometric solutions for average rigid-body rotation,” Mech. Mach. Theory, 45 1239 –1251 (2010). http://dx.doi.org/10.1016/j.mechmachtheory.2010.05.002 MHMTAS 0094-114X Google Scholar

27.

F. C. ParkJ. Kim, “Geometric descent algorithms for attitude determination using the global positioning system,” J. Guid. Control Dyn., 23 (1), 26 –33 (2000). http://dx.doi.org/10.2514/2.4516 JGCODS 0731-5090 Google Scholar

28.

D. MortariL. F. MarkleyP. Singla, “Optimal linear attitude estimator,” J. Guid. Control Dyn., 30 (6), 1619 –1627 (2007). http://dx.doi.org/10.2514/1.29568 JGCODS 0731-5090 Google Scholar

29.

J. GallierD. Xu, “Computing exponentials of skew-symmetric matrices and logarithms of orthogonal matrices,” Int. J. Rob. Autom., 17 (4), 10 –20 (2002). IJAUED Google Scholar

30.

C. MolerC. van Loanz, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” SIAM Rev., 45 (1), 3 –49 (2003). http://dx.doi.org/10.1137/S00361445024180 SIREAD 0036-1445 Google Scholar

31.

I. Kruzhilovet al., “Flight experience of 329K star tracker,” Proc. SPIE, 8889 888920 (2013). http://dx.doi.org/10.1117/12.2028301 Google Scholar

32.

E. Andersonet al., LAPACK User’s Guide, 3rd Ed.SIAM, Philadelphia (1999). Google Scholar

33.

J. D. Hoffman, Numeric Methods for Engineers and Scientists, 47 Marcel Dekker, New York (2001). Google Scholar

34.

I. V. Sergienko, Methods of Optimization and Systems Analysis for Problems of Transcomputational Complexity, 42 Springer Science+Business Media, New York (2012). Google Scholar

Biography

Ivan S. Kruzhilov received his MSc degree from the Moscow Power Engineering Institute in 2007. From 2005 to 2006, he was a scholarship holder of Siemens AG and studied at the Technical University Ilmenau. He finished his PhD degree in 2010 from the Moscow Power Engineering Institute. He took part in development of star tracker for Russian satellites “Luch-5A” and “Luch-5B” and multiple-head star tracker 348K developed at the Research and Production Enterprise Geofizika-Cosmos. He is a regular member of SPIE. His research interests are in applied mathematics and mathematical simulation.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ivan S. Kruzhilov "Small-angle rotation method for star tracker orientation," Journal of Applied Remote Sensing 7(1), 073479 (19 November 2013). https://doi.org/10.1117/1.JRS.7.073479

Published: 19 November 2013

Access the abstract

JOURNAL ARTICLE
20 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY