Dot Product Equality Constrained Attitude Determination from Two Vector Observations: Theory and Astronautical Applications

Abstract

In this paper, the attitude determination problem from two vector observations is revisited, incorporating the redundant equality constraint obtained by the dot product of vector observations. Analytical solutions to this constrained attitude determination problem are derived. It is found out that the studied two-vector attitude determination problem by Davenport q-method under the dot product constraint has deterministic maximum eigenvalue, which leads to its advantage in error/perturbation analysis and covariance determination. The proposed dot product constrained two-vector attitude solution is applied then to solve several engineering problems. Detailed simulations on spacecrafts attitude determination indicate the efficiency of the proposed theory.

1. Introduction

Attitude determination from vector observations is usually employed in astronautical applications for state estimation and control of spacecrafts [1,2,3]. Commonly, the vector observations are acquired from vector-measurement sensors e.g., Sun sensor, Nadir sensor, magnetometer, star tracker, etc. [4,5,6]. The usage of direct attitude determination from such sensors has the advantage of compensating for biases in rate gyroscopes to cancel the long-endurance drifts inside inertial navigation results [7,8].

In the past several decades, many algorithms for optimal attitude determination have been developed [9,10]. Apart from some intuitive geometric approaches, most of the others are based on the Wahba’s problem posed in 1965, giving optimal attitude matrix estimates with normalized vector measurements along with their statistical weights [11]. There have been a lot of Wahba’s solvers developed in the past 30 years, including celebrated ones such as the QUaternion ESTimator (QUEST [12]), singular value decomposition (SVD [13]) and Euler-q [14], and recent ones such as the vector-inertia-based method by Patera [15], the Riemannian-manifold-based on by Yang [16], the optimal linear estimator of quaternion (OLEQ [17]), etc. These solver are all general ones facing the multivector cases. Actually, when there are few pairs of vector observations, the Wahba’s optimization can also be solved via numerical algorithms like gradient-descent algorithm (GDA [18]), Gauss–Newton algorithm (GNA [19]), Levenberg–Marquardt algorithm (LMA [20]), etc. Such optimization-based solvers are popular in mechatronic navigation tasks with accelerometers and magnetometers only [21].

It was found that with two vector observation pairs, we can compute the full attitude results using the methods such as TRIAD [22]. For instance, with a sun sensor and a magnetometer, the full attitude angles along with the orbit can already be determined accurately [23,24]. Based on the two-vector case, Markley has developed two approaches [25,26]. For the very first one, the attitude quaternion is employed for parameterizing the orientation transformation; and, in the second one, the optimal attitude matrix of such problem has been successfully found out. Markley’s methods are proven to be efficient in engineering practice. These two optimized approaches are all based on conventional Wahba’s framework and can not improve the essential accuracy. However, for some two-vector based attitude determination system, e.g., accelerometer–magnetometer-based one, the magnetometer’s reference vector for true-north-finding is determined by the absolute position of the object in Earth coordination system [27]. Commonly, when there is no magnetic distortion, the attitude results will be quite reliable. When adding some magnetic interferences, all Euler angles using a true-north reference vector will all be disturbed, thus generating large attitude biases. However, the accelerometer can independently determine the roll and pitch angles, which indicates that the disturbances from magnetic sensing can be separated. In the work by the authors of [28], this is achieved by introducing an equality constraint formed by vector dot product from rigid-body motion; whereas, the two-vector attitude determination system is far beyond the single accelerometer–magnetometer one, it is believed that there might be some more generalized theory inside such two-vector system with dot product equality constraint, as some system dynamics will be interesting by incorporating external equality constraints, e.g., the works by the authors of [29,30,31].

Guided by the aforementioned problem, in this paper, we present a novel attitude determination algorithm from two vector observations. The framework is still within the Wahba’s optimization, but is penalized by adding a dot product equality constraint derived from vector relations. We also show that under the proposed conditions of dot product constraint, the loss function of Wahba’s problem can be reduced to 0, i.e., the sensors can extract self information to for possible maximum separation of different sensor effects and angle observabilities. The proposed two-vector approach can be an alternative for spacecraft attitude determination backup that subjects to the minimum equipment list (MEL).

The remainder of the paper is organized as follows. Section 2 contains problem backgrounds. Section 3 includes our main results on the studied problem. Section 4 consists an astronautical application and some related experimental results. In Section 5, we draw the concluding remarks.

2. Problem Fomulation

Given two normalized vector observation pairs in the body frame b and reference frame r, one can relate them by

$$\begin{array}{c}{\mathit{b}}_1=\mathit{C}{\mathit{r}}_1\hfill \\ {\mathit{b}}_2=\mathit{C}{\mathit{r}}_2\hfill \end{array}$$

(1)

where $\mathit{C}$ in the special orthogonal group $SO\left(3\right):=\left\{\mathit{C}\right|{\mathit{C}}^T\mathit{C}=\mathit{I},det\left(\mathit{C}\right)=1\}$ is the direction cosine matrix (DCM); the sensor measurements are given by

$$\left\{\begin{array}{c}{\mathit{b}}_1={\left(b_{x,1},b_{y,1},b_{z,1}\right)}^T\\ {\mathit{r}}_1={\left(r_{x,1},r_{y,1},r_{z,1}\right)}^T\end{array}\right.,\left\{\begin{array}{c}{\mathit{b}}_2={\left(b_{x,2},b_{y,2},b_{z,2}\right)}^T\\ {\mathit{r}}_2={\left(r_{x,2},r_{y,2},r_{z,2}\right)}^T\end{array}\right.$$

(2)

provided that

$$\left\{\begin{array}{c}{b}_{x,1}^2+b_{y,1}^2+b_{z,1}^2=1\\ r_{x,1}^2+r_{y,1}^2+r_{z,1}^2=1\\ b_{x,2}^2+b_{y,2}^2+b_{z,2}^2=1\\ r_{x,2}^2+r_{y,2}^2+r_{z,2}^2=1\end{array}\right.$$

(3)

In astronautical missions, the spacecraft is always equipped with a Sun sensor, magnetometer, horizon sensor, etc. for attitude determination (see Figure 1). Wahba’s problem aims to find out the optimal DCM of Equation (1) by the following least-square optimization.

$$\underset{\mathit{C}\in SO\left(3\right)}{argmin}\left[w{∥{\mathit{b}}_1-\mathit{C}{\mathit{r}}_1∥}^2+(1-w){∥{\mathit{b}}_2-\mathit{C}{\mathit{r}}_2∥}^2\right]$$

(4)

where w is the weight belonging to the first vector observation pair. This optimization can be reduced to finding the optimal eigenvector of the following matrix $\mathit{K}$ [32,33],

$$\mathit{K}=\left[\begin{array}{cc}\mathit{B}+{\mathit{B}}^T-tr\left(\mathit{B}\right)\mathit{I}& \mathit{z}\\ {\mathit{z}}^T& tr\left(\mathit{B}\right)\end{array}\right]$$

(5)

such that

$$\mathit{Kq}={\lambda }_{max}\mathit{q}$$

(6)

where ${\lambda }_{max}$ is the largest eigenvalue of $\mathit{K}$ while other parameters are given by

$$\begin{array}{c}\mathit{B}=w{\mathit{b}}_1^{}{\mathit{r}}_1^T+(1-w){\mathit{b}}_2^{}{\mathit{r}}_2^T\hfill \\ \mathit{z}=w{\mathit{b}}_1^{}\times {\mathit{r}}_1^{}+(1-w){\mathit{b}}_2^{}\times {\mathit{r}}_2^{}\hfill \end{array}$$

(7)

The characteristic polynomial of $\mathit{K}$ can be computed by

$$det(\mathit{K}-\lambda \mathit{I})=0$$

(8)

In fact, many closed-form formulations of this polynomial have been developed. For instance, in ESOQ2, it has been given that [34,35]

$${\lambda }^4+{\tau }_1{\lambda }^2+{\tau }_2\lambda +{\tau }_3=0$$

(9)

where

$$\begin{array}{c}{\tau }_1=-2t{r}^2\left(\mathit{B}\right)+tr\left[adj(\mathit{B}+{\mathit{B}}^T)\right]-{\mathit{z}}^T\mathit{z}\hfill \\ {\tau }_2=-tr\left[adj\left(\mathit{K}\right)\right]\hfill \\ {\tau }_3=det\left(\mathit{K}\right)\hfill \end{array}$$

(10)

where $adj$ denotes the adjoint matrix. In the following sections, we derive the analytical solution to ${\lambda }_{max}$ while generating some internal perspectives on the relationship between body and reference vectors.

3. Proposed Theory

3.1. Dot Product-Equality Constraint

Theorem 1.

Defining

$$\left(\begin{array}{ccc}{H}_{x1}& H_{x2}& H_{x3}\\ H_{y1}& H_{y2}& H_{y3}\\ H_{z1}& H_{z2}& H_{z3}\end{array}\right)={\mathit{B}}^T$$

(11)

the following matrix [36],

$$\begin{array}{ll}{\mathit{W}}_{1,1}=H_{x1}+H_{y2}+H_{z3}\hfill & {\mathit{W}}_{1,2}=-H_{y3}+H_{z2}\\ {\mathit{W}}_{1,3}=-H_{z1}+H_{x3}\hfill & {\mathit{W}}_{1,4}=-H_{x2}+H_{y1}\\ {\mathit{W}}_{2,1}=-H_{y3}+H_{z2}\hfill & {\mathit{W}}_{2,2}=H_{x1}-H_{y2}-H_{z3}\\ {\mathit{W}}_{2,3}=H_{x2}+H_{y1}\hfill & {\mathit{W}}_{2,4}=H_{x3}+H_{z1}\\ {\mathit{W}}_{3,1}=-H_{z1}+H_{x3}\hfill & {\mathit{W}}_{3,2}=H_{x2}+H_{y1}\\ {\mathit{W}}_{3,3}=H_{y2}-H_{x1}-H_{z3}\hfill & {\mathit{W}}_{3,4}=H_{y3}+H_{z2}\\ {\mathit{W}}_{4,1}=-H_{x2}+H_{y1}\hfill & {\mathit{W}}_{4,2}=H_{x3}+H_{z1}\\ {\mathit{W}}_{4,3}=H_{y3}+H_{z2}\hfill & {\mathit{W}}_{4,4}=H_{z3}-H_{y2}-H_{x1}\end{array}$$

(12)

where ${\mathit{W}}_{i,j}$ denotes the matrix entry of $\mathit{W}$ with row and column indices of i and j respectively, has the similar structure with $\mathit{K}$ and owns the same eigenvalues and associated eigenvectors.

Proof. 

We can directly obtain

$$\begin{array}{c}\mathit{B}+{\mathit{B}}^T-tr\left(\mathit{B}\right)\mathit{I}=\left(\begin{array}{ccc}{H}_{x1}-H_{y2}-H_{z3}& H_{x2}+H_{y1}& H_{x3}+H_{z1}\\ H_{x2}+H_{y1}& H_{y2}-H_{x1}-H_{z3}& H_{y3}+H_{z2}\\ H_{x3}+H_{z1}& H_{y3}+H_{z2}& H_{z3}-H_{y2}-H_{x1}\end{array}\right)\hfill \\ tr\left(\mathit{B}\right)=H_{x1}+H_{y2}+H_{z3}\hfill \\ \mathit{z}={\left({\mathit{B}}_{2,3}-{\mathit{B}}_{3,2},{\mathit{B}}_{3,1}-{\mathit{B}}_{1,3},{\mathit{B}}_{1,2}-{\mathit{B}}_{2,1}\right)}^T={\left(-H_{y3}+H_{z2},-H_{z1}+H_{x3},-H_{x2}+H_{y1}\right)}^T\hfill \end{array}$$

(13)

Comparing the above results with Equation (12), we have

$$\mathit{W}=\left[\begin{array}{cc}tr\left(\mathit{B}\right)& {\mathit{z}}^T\\ \mathit{z}& \mathit{B}+{\mathit{B}}^T-tr\left(\mathit{B}\right)\mathit{I}\end{array}\right]$$

(14)

It is evident that $\mathit{W}$ and $\mathit{K}$ has the similar structure. This completes the proof. □

Lemma 1.

${\tau }_2=0$ holds for arbitrary two pairs of vector observations.

Proof. 

The diagonal elements of $\mathit{K}$ matrix can be given by

$$\begin{array}{c}{\mathit{K}}_{1,1}=\left(b_{z,2}{r}_{z,2}+b_{y,2}{r}_{y,2}-b_{x,2}{r}_{x,2}\right)(w-1)+\left(b_{x,1}{r}_{x,1}-b_{y,1}{r}_{y,1}-b_{z,1}{r}_{z,1}\right)w\hfill \\ {\mathit{K}}_{2,2}=\left(b_{z,2}{r}_{z,2}-b_{y,2}{r}_{y,2}+b_{x,2}{r}_{x,2}\right)(w-1)+\left(-b_{x,1}{r}_{x,1}+b_{y,1}{r}_{y,1}-b_{z,1}{r}_{z,1}\right)w\hfill \\ {\mathit{K}}_{3,3}=\left(-b_{z,2}{r}_{z,2}+b_{y,2}{r}_{y,2}+b_{x,2}{r}_{x,2}\right)(w-1)+\left(-b_{x,1}{r}_{x,1}-b_{y,1}{r}_{y,1}+b_{z,1}{r}_{z,1}\right)w\hfill \\ {\mathit{K}}_{4,4}=\left(-b_{z,2}{r}_{z,2}-b_{y,2}{r}_{y,2}-b_{x,2}{r}_{x,2}\right)(w-1)+\left(b_{x,1}{r}_{x,1}+b_{y,1}{r}_{y,1}+b_{z,1}{r}_{z,1}\right)w\hfill \end{array}$$

(15)

Therefore, it is obvious that the sum of the eigenvalues to $\mathit{K}$ equals to ${\mathit{K}}_{1,1}+{\mathit{K}}_{2,2}+{\mathit{K}}_{3,3}+{\mathit{K}}_{4,4}=0$. This is why the coefficient of ${\lambda }^3$ is zero. Note that ${\tau }_2$ is de facto in the form of [36]

$${\tau }_2=8\left(\begin{array}{c}{H}_{x3}{H}_{y2}{H}_{z1}-H_{x2}{H}_{y3}{H}_{z1}-H_{x3}{H}_{y1}{H}_{z2}+\hfill \\ H_{x1}{H}_{y3}{H}_{z2}+H_{x2}{H}_{y1}{H}_{z3}-H_{x1}{H}_{y2}{H}_{z3}\hfill \end{array}\right)$$

(16)

Actually, the right part can be represented by

$${\tau }_2=-8det\left(\mathit{B}\right)$$

(17)

In view of $\mathit{B}$ for the two-vector case, we can see that it takes the following form

$$\begin{array}{c}\mathit{B}=w{\mathit{b}}_1^{}{\mathit{r}}_1^T+(1-w){\mathit{b}}_2^{}{\mathit{r}}_2^T\Rightarrow rank\left(\mathit{B}\right)\le rank\left(w{\mathit{b}}_1^{}{\mathit{r}}_1^T\right)+rank\left[(1-w){\mathit{b}}_2^{}{\mathit{r}}_2^T\right]=2\hfill \end{array}$$

(18)

indicating that $\mathit{B}$ is rank-deficient. Therefore, we arrive at

$${\tau }_2=-8det\left(\mathit{B}\right)=0$$

(19)

which completes the proof. □

With Lemma 1, we can see that the characteristic polynomial is essentially a quadratic equation of ${\lambda }^2$. Then the eigenvalues can be solved immediately via

$$\lambda =\pm \sqrt{\frac{-{\tau }_1\pm \sqrt{{\tau }_1^2-4{\tau }_3}}{2}}$$

(20)

The maximum eigenvalue should be

$${\lambda }_{max}^{}=\sqrt{\frac{-{\tau }_1+\sqrt{{\tau }_1^2-4{\tau }_3}}{2}}$$

(21)

because [36]

$${\tau }_1=-2\left(\begin{array}{c}{H}_{x1}^2+H_{x2}^2+H_{x3}^2+\hfill \\ H_{y1}^2+H_{y2}^2+H_{y3}^2+\hfill \\ H_{z1}^2+H_{z2}^2+H_{z3}^2\hfill \end{array}\right)<0$$

(22)

As the q-method is derived from the Lagrange multiplier, with increasing maximum eigenvalue, the loss function becomes even smaller. The maximum eigenvalue has an upper bound, as proven by Shuster that [12]

$${\lambda }_{{max}^{}}⩽\sum _{i=1}^n{w}_i=1$$

(23)

where $w_i$ is the weight of the i-th vector observation pair. Regarding ${\lambda }_{max}$ as a function of the vector observations, we can seek its largest point. This is equivalent to figuring out the solution to the following system,

$$\left\{\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{x,1}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{y,1}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{z,1}}=0\end{array}\right.,\left\{\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{x,2}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{y,2}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial b_{z,2}}=0\end{array},\left\{\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{x,1}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{y,1}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{z,1}}=0\end{array}\right.,\left\{\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{x,2}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{y,2}}=0\\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r_{z,2}}=0\end{array}\right.\right.$$

(24)

which is rather difficult in detail. Let us write out the algebraic form of ${\tau }_1$ and ${\tau }_3$:

$$\begin{array}{c}{\tau }_1=-2+4w-4w^2+\hfill \\ \begin{array}{ccc}& & \end{array}4w(w-1)\left(b_{x,1}{b}_{x,2}+b_{y,1}{b}_{y,2}+b_{z,1}{b}_{z,2}\right)\left(r_{x,1}{r}_{x,2}+r_{y,1}{r}_{y,2}+r_{z,1}{r}_{z,2}\right)\hfill \end{array}$$

(25a)

$$\begin{array}{c}{\tau }_3=1-4w+12w^2-16w^3+8w^4-\hfill \\ \begin{array}{ccc}& & \end{array}4w(w-1)\left(b_{x,1}{b}_{x,2}+b_{y,1}{b}_{y,2}+b_{z,1}{b}_{z,2}\right)\left(r_{x,1}{r}_{x,2}+r_{y,1}{r}_{y,2}+r_{z,1}{r}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}8w^2(w-1)\left(b_{x,1}{b}_{x,2}+b_{y,1}{b}_{y,2}+b_{z,1}{b}_{z,2}\right)\left(r_{x,1}{r}_{x,2}+r_{y,1}{r}_{y,2}+r_{z,1}{r}_{z,2}\right)-\hfill \\ \begin{array}{ccc}& & \end{array}8w^3(w-1)\left(b_{x,1}{b}_{x,2}+b_{y,1}{b}_{y,2}+b_{z,1}{b}_{z,2}\right)\left(r_{x,1}{r}_{x,2}+r_{y,1}{r}_{y,2}+r_{z,1}{r}_{z,2}\right)-\hfill \\ \begin{array}{ccc}& & \end{array}4w^2{(w-1)}^2\left(\begin{array}{c}{b}_{y,1}^2-2b_{x,1}{b}_{x,2}{b}_{y,1}{b}_{y,2}+b_{y,2}^2-2b_{y,1}^2{b}_{y,2}^2+\hfill \\ b_{z,1}^2-b_{y,2}^2{b}_{z,1}^2-2b_{x,1}{b}_{x,2}{b}_{z,1}{b}_{z,2}-2b_{y,1}{b}_{y,2}{b}_{z,1}{b}_{z,2}+\hfill \\ b_{z,2}^2-b_{y,1}^2{b}_{z,2}^2-2b_{z,1}^2{b}_{z,2}^2+r_{y,1}^2-2r_{x,1}{r}_{x,2}{r}_{y,1}{r}_{y,2}+\hfill \\ r_{y,2}^2-2r_{y,1}^2{r}_{y,2}^2+r_{z,1}^2-r_{y,2}^2{r}_{z,1}^2-2r_{x,1}{r}_{x,2}{r}_{z,1}{r}_{z,2}-\hfill \\ 2r_{y,1}{r}_{y,2}{r}_{z,1}{r}_{z,2}+r_{z,2}^2-r_{y,1}^2{r}_{z,2}^2-2r_{z,1}^2{r}_{z,2}^2\hfill \end{array}\right)\hfill \end{array}$$

(25b)

Note that we have the following simplification.

$$\begin{array}{c}{b}_{y,1}^2+b_{z,1}^2+b_{y,2}^2+b_{z,2}^2-b_{y,1}^2{b}_{y,2}^2-b_{y,2}^2{b}_{z,1}^2-b_{y,1}^2{b}_{z,2}^2-b_{z,1}^2{b}_{z,2}^2\hfill \\ =b_{y,1}^2+b_{z,1}^2+b_{y,2}^2+b_{z,2}^2-\left(b_{y,1}^2+b_{z,1}^2\right)\left(b_{y,2}^2+b_{z,2}^2\right)\hfill \\ =2-b_{x,1}^2-b_{x,2}^2-\left(1-b_{x,1}^2\right)\left(1-b_{x,2}^2\right)\hfill \\ =2-b_{x,1}^2-b_{x,2}^2-1-b_{x,1}^2{b}_{x,2}^2+b_{x,1}^2+b_{x,2}^2\hfill \\ =1-b_{x,1}^2{b}_{x,2}^2\hfill \end{array}$$

(26)

With

$$\begin{array}{c}p={\mathit{b}}_1·{\mathit{b}}_2=b_{x,1}{b}_{x,2}+b_{y,1}{b}_{y,2}+b_{z,1}{b}_{z,2}\hfill \\ q={\mathit{r}}_1·{\mathit{r}}_2=r_{x,1}{r}_{x,2}+r_{y,1}{r}_{y,2}+r_{z,1}{r}_{z,2}\hfill \end{array}$$

(27)

${\tau }_1$ and ${\tau }_3$ are simplified to

$$\begin{array}{c}{\tau }_1=-2+4w-4w^2+4w(w-1)pq\hfill \\ {\tau }_3=1-4w+12w^2-16w^3+8w^4-\hfill \\ \begin{array}{ccc}& & \end{array}4w\left(1-2w+2w^2\right)(w-1)pq-4w^2{(w-1)}^2\left(2-p^2-q^2\right)\hfill \end{array}$$

(28)

Therefore, it is shown that the eigenvalue is actually the function of p and q, with parameters ranging from

$$\begin{array}{c}p={\mathit{b}}_1·{\mathit{b}}_2\in \left[-1,1\right]\hfill \\ q={\mathit{r}}_1·{\mathit{r}}_2\in \left[-1,1\right]\hfill \\ w\in \left(0,1\right)\hfill \end{array}$$

(29)

The figures of the maximum eigenvalue Equation (21) with $w=0.5$ are depicted as follows in Figure 2.

With different w, the maximum eigenvalues are drawn in Figure 3.

From the above figures, we see that the maximum eigenvalue is a convex function with respect to p and q no matter how w changes. Based on this point, the global optimum-seeking procedure is then reduced to solve the following system.

$$\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial p}=0\hfill \\ \frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial q}=0\hfill \end{array}$$

(30)

The partial differentiation can be computed by

$$\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r}=-\frac{\partial {\tau }_1}{\partial r}+\frac{1}{2\sqrt{{\tau }_1^2-4{\tau }_3}}\left(2{\tau }_1\frac{\partial {\tau }_1}{\partial r}-4\frac{\partial {\tau }_3}{\partial r}\right)$$

(31)

where r is the differentiation independent variables:

$$r=p,q$$

(32)

By which the specific derivatives are obtained by

$$\left\{\begin{array}{c}\frac{\partial {\tau }_1}{\partial p}=4w(w-1)q\\ \frac{\partial {\tau }_1}{\partial q}=4w(w-1)p\end{array}\right.,\left\{\begin{array}{c}\frac{\partial {\tau }_3}{\partial p}=4w\left(1-2w+2w^2\right)(w-1)q+8w^2{(w-1)}^{2}p\\ \frac{\partial {\tau }_3}{\partial q}=4w\left(1-2w+2w^2\right)(w-1)p+8w^2{(w-1)}^{2}q\end{array}\right.$$

(33)

Inserting these equalities into Equation (30) and invoking

$$\begin{array}{c}\frac{\partial \left(2{\lambda }_{max}^2\right)}{\partial r}=0\hfill \\ \Rightarrow \frac{\partial {\tau }_1}{\partial r}=\frac{1}{\sqrt{{\tau }_1^2-4{\tau }_3}}\left({\tau }_1\frac{\partial {\tau }_1}{\partial r}-2\frac{\partial {\tau }_3}{\partial r}\right)\hfill \\ \Rightarrow \left(\frac{\sqrt{{\tau }_1^2-4{\tau }_3}}{2}-\frac{{\tau }_1}{2}\right)\frac{\partial {\tau }_1}{\partial r}=\frac{\partial {\tau }_3}{\partial r}\hfill \end{array}$$

(34)

we obtain

$$p=q$$

(35)

at which, we also have

$$\sqrt{{\tau }_1^2-4{\tau }_3}=4w-4w^2+4w(w-1)pq$$

(36)

leading to

$${\lambda }_{max}=1$$

(37)

Reconsidering the forms of p and q, we find out

$$p=q\Rightarrow {\mathit{b}}_1·{\mathit{b}}_2={\mathit{r}}_1·{\mathit{r}}_2$$

(38)

which can actually be derived from another aspect by

$${\mathit{b}}_1·{\mathit{b}}_2={\mathit{b}}_1^T{\mathit{b}}_2={\mathit{r}}_1^T{\mathit{C}}^T\mathit{C}{\mathit{r}}_2={\mathit{r}}_1^T{\mathit{r}}_2={\mathit{r}}_1·{\mathit{r}}_2$$

(39)

This finding reflects that the equality constraint Equation (38) can achieve least value of Wahba’s loss function, such that [32]

$$\begin{array}{c}L\left(\mathit{C}\right)=w{∥{\mathit{b}}_1-\mathit{C}{\mathit{r}}_1∥}^2+(1-w){∥{\mathit{b}}_2-\mathit{C}{\mathit{r}}_2∥}^2\hfill \\ =1-tr\left(\mathit{C}{\mathit{B}}^T\right)\hfill \\ =1-{\mathit{q}}^T\mathit{Kq}\hfill \\ =1-{\mathit{q}}^T{\lambda }_{max}\mathit{q}\hfill \\ =0\hfill \end{array}$$

(40)

The occurrence of ${\lambda }_{max}=1$ has in fact been shown by Markley [26]. However, we need to note that the attitude matrix employed for attitude representation in the work by the authors of [26] is not as flexible as the quaternion used in this paper, as a quaternion only requires normalization rather than orthonormalization for an attitude matrix [37].

3.2. Quaternion Solution

As the loss function is zero, we can regard the optimization as a properly solved one. That is, at this time, the information from each sensor has been maximumly extracted to reduce the loss function value; and, geometrically, every sensor has fully contributed to all those angles that they can observe. This means the optimal quaternion is finally computed via

$$\left(\mathit{K}-\mathit{I}\right)\mathit{q}=\mathit{0}$$

(41)

generating the following row-echelon form

$$\mathit{K}-\mathit{I}\Rightarrow \left(\begin{array}{cccc}1& & & \alpha /N\\ & 1& & \beta /N\\ & & 1& \gamma /N\\ & & & 0\end{array}\right)$$

(42)

The symbolic solutions to this function can be obtained by

$$\begin{array}{c}N=1-b_{y,1}{r}_{y,1}-b_{y,2}{r}_{y,2}+b_{z,1}{r}_{z,1}+b_{z,2}{r}_{z,2}\hfill \\ \begin{array}{ccc}& & \end{array}+\left(b_{x,1}{r}_{y,1}-b_{x,1}{r}_{x,1}\right)r_{y,2}^2+q^2\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{x,2}{r}_{y,1}\left(b_{y,2}{r}_{x,1}-b_{x,2}{r}_{y,1}+b_{x,2}{b}_{y,1}-b_{x,1}{b}_{y,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{y,2}{r}_{y,1}\left(-b_{x,1}{r}_{y,2}+b_{x,1}{r}_{x,2}+b_{y,2}{r}_{y,1}+b_{y,1}{r}_{y,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{x,1}{r}_{y,2}\left(b_{y,1}{r}_{x,2}+b_{x,2}{r}_{y,1}-b_{x,2}{b}_{y,1}+b_{x,1}{b}_{y,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{x,2}{r}_{z,1}\left(-b_{z,2}{r}_{x,1}-b_{x,2}{r}_{z,1}-b_{x,2}{b}_{z,1}+b_{x,1}{b}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{x,1}{r}_{z,2}\left(-b_{z,1}{r}_{x,2}-b_{x,1}{r}_{z,2}+b_{x,2}{b}_{z,1}-b_{x,1}{b}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{z,1}{r}_{z,2}\left(b_{x,2}{r}_{x,1}+b_{x,1}{r}_{x,2}-b_{z,2}{r}_{z,1}-b_{z,1}{r}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{y,1}{r}_{z,2}\left(-b_{z,1}{r}_{y,2}+b_{y,2}{r}_{z,1}+b_{y,2}{b}_{z,1}-b_{y,1}{b}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{r}_{y,2}{r}_{z,1}\left(-b_{z,2}{r}_{y,1}+b_{y,1}{r}_{z,2}-b_{y,2}{b}_{z,1}+b_{y,1}{b}_{z,2}\right)\hfill \end{array}$$

(43a)

$$\begin{array}{c}\alpha =b_{z,2}{r}_{x,1}{r}_{x,2}{r}_{y,1}+b_{x,2}{b}_{y,1}{r}_{x,2}{r}_{z,1}-b_{y,2}{r}_{x,1}{r}_{x,2}{r}_{z,1}-b_{x,2}{b}_{y,1}{r}_{x,1}{r}_{z,2}-\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{y,1}{r}_{x,1}{r}_{x,2}{r}_{z,2}+b_{x,1}\left(b_{z,2}{r}_{x,2}{r}_{y,1}-b_{z,2}{r}_{x,1}{r}_{y,2}-b_{y,2}{r}_{x,2}{r}_{z,1}+b_{y,2}{r}_{x,1}{r}_{z,2}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{z,1}\left[r_{y,2}\left(b_{x,2}{r}_{x,1}+r_{x,2}{r}_{x,1}+r_{z,1}{r}_{z,2}\right)+r_{y,1}\left(-b_{x,2}{r}_{x,2}+r_{y,2}^2-1\right)\right]+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{z,2}{r}_{y,2}{r}_{y,1}^2-b_{y,2}{r}_{y,2}{r}_{y,1}{r}_{z,1}-b_{y,1}{r}_{y,2}{r}_{y,1}{r}_{z,2}+b_{z,2}{r}_{y,1}{r}_{z,1}{r}_{z,2}-\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{y,1}{r}_{z,1}{r}_{z,2}^2-b_{z,2}{r}_{y,2}+b_{y,1}{r}_{z,1}-b_{y,2}{r}_{z,1}^2{r}_{z,2}+b_{y,2}{r}_{z,2}\hfill \end{array}$$

(43b)

$$\begin{array}{c}\beta =b_{z,2}{r}_{x,2}{r}_{y,1}^2-b_{z,2}{r}_{x,1}{r}_{y,2}{r}_{y,1}-b_{z,1}{r}_{x,2}{r}_{y,2}{r}_{y,1}+b_{x,2}{r}_{y,2}{r}_{y,1}{r}_{z,1}+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{x,1}{r}_{y,2}{r}_{y,1}{r}_{z,2}+b_{z,1}{r}_{x,1}{r}_{y,2}^2+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{y,2}\left(-b_{z,1}{r}_{x,2}{r}_{y,1}+b_{x,1}{r}_{y,1}{r}_{z,2}+b_{z,1}{r}_{x,1}{r}_{y,2}-b_{x,1}{r}_{y,2}{r}_{z,1}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{y,1}\left(b_{z,2}{r}_{x,2}{r}_{y,1}-b_{x,2}{r}_{y,1}{r}_{z,2}-b_{z,2}{r}_{x,1}{r}_{y,2}+b_{x,2}{r}_{y,2}{r}_{z,1}\right)+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{z,2}{r}_{x,2}{r}_{z,1}^2+b_{z,1}{r}_{x,1}{r}_{z,2}^2+b_{x,1}{r}_{z,1}{r}_{z,2}^2-b_{x,1}{r}_{z,1}+b_{x,2}{r}_{x,1}{r}_{x,2}{r}_{z,1}+\hfill \\ \begin{array}{ccc}& & \end{array}{b}_{x,2}{r}_{z,1}^2{r}_{z,2}-b_{x,2}{r}_{z,2}+b_{x,1}{r}_{x,1}{r}_{x,2}{r}_{z,2}-b_{z,2}{r}_{x,1}{r}_{z,1}{r}_{z,2}-b_{z,1}{r}_{x,2}{r}_{z,1}{r}_{z,2}\hfill \end{array}$$

(43c)

$$\gamma =-\left\{\begin{array}{c}{b}_{y,2}{r}_{x,2}{r}_{z,1}^2+b_{y,1}{r}_{x,1}{r}_{z,2}^2+b_{y,2}{b}_{z,1}{r}_{x,2}{r}_{z,1}-b_{y,1}{b}_{z,2}{r}_{x,2}{r}_{z,1}-\hfill \\ b_{y,2}{b}_{z,1}{r}_{x,1}{r}_{z,2}+b_{y,1}{b}_{z,2}{r}_{x,1}{r}_{z,2}-b_{y,2}{r}_{x,1}{r}_{z,1}{r}_{z,2}-b_{y,1}{r}_{x,2}{r}_{z,1}{r}_{z,2}+\hfill \\ b_{x,2}\left[r_{y,1}{r}_{z,2}\left(b_{z,1}+r_{z,1}\right)+r_{y,2}\left(-b_{z,1}{r}_{z,1}+r_{y,1}^2-1\right)+r_{x,1}{r}_{x,2}{r}_{y,1}\right]+\hfill \\ b_{x,1}\left[r_{y,2}\left(r_{z,1}{b}_{z,2}+r_{z,1}{r}_{z,2}+r_{x,1}{r}_{x,2}\right)+r_{y,1}\left(-b_{z,2}{r}_{z,2}+r_{y,2}^2-1\right)\right]+\hfill \\ b_{y,2}{r}_{x,2}{r}_{y,1}^2-b_{y,2}{r}_{x,1}{r}_{y,2}{r}_{y,1}-b_{y,1}{r}_{x,2}{r}_{y,2}{r}_{y,1}+b_{y,1}{r}_{x,1}{r}_{y,2}^2\hfill \end{array}\right\}$$

(43d)

The final quaternion solution is obtained by

$$\mathit{q}=\frac{1}{\sqrt{{\alpha }^2+{\beta }^2+{\gamma }^2+N^2}}{\left(\alpha ,\beta ,\gamma ,-N\right)}^T$$

(44)

We can see that the obtained quaternion is weight-free. A well-known two-vector quaternion solution is given by Markley [25], such that

$$\mathit{q}=\left\{\begin{array}{c}\frac{1}{2\sqrt{\gamma \left(\gamma +\alpha \right)\left(1+{\mathit{b}}_3·{\mathit{r}}_3\right)}}\left[\begin{array}{c}\left(\gamma +\alpha \right)\left({\mathit{b}}_3\times {\mathit{r}}_3\right)+\beta \left({\mathit{b}}_3+{\mathit{r}}_3\right)\\ \left(\gamma +\alpha \right)\left(1+{\mathit{b}}_3·{\mathit{r}}_3\right)\end{array}\right]\\ \frac{1}{2\sqrt{\gamma \left(\gamma -\alpha \right)\left(1+{\mathit{b}}_3·{\mathit{r}}_3\right)}}\left[\begin{array}{c}\beta \left({\mathit{b}}_3\times {\mathit{r}}_3\right)+\left(\gamma -\alpha \right)\left({\mathit{b}}_3+{\mathit{r}}_3\right)\\ \beta \left(1+{\mathit{b}}_3·{\mathit{r}}_3\right)\end{array}\right]\end{array}\begin{array}{cc}for& \alpha ⩾0\\ for& \alpha ⩽0\end{array}\right.$$

(45)

where

$$\begin{array}{c}{\mathit{b}}_3=\frac{{\mathit{b}}_1\times {\mathit{b}}_2}{∥{\mathit{b}}_1\times {\mathit{b}}_2∥},{\mathit{r}}_3=\frac{{\mathit{r}}_1\times {\mathit{r}}_2}{∥{\mathit{r}}_1\times {\mathit{r}}_2∥}\hfill \\ \alpha =\left(1+{\mathit{b}}_3·{\mathit{r}}_3\right)\left[w{\mathit{b}}_1·{\mathit{r}}_1+(1-w){\mathit{b}}_2·{\mathit{r}}_2\right]+\hfill \\ \begin{array}{cc}& \end{array}\left({\mathit{b}}_3\times {\mathit{r}}_3\right)·\left[w{\mathit{b}}_1\times {\mathit{r}}_1+(1-w){\mathit{b}}_2\times {\mathit{r}}_2\right]\hfill \\ \beta =\left({\mathit{b}}_3+{\mathit{r}}_3\right)·\left[w{\mathit{b}}_1\times {\mathit{r}}_1+(1-w){\mathit{b}}_2\times {\mathit{r}}_2\right]\hfill \\ \gamma =\sqrt{{\alpha }^2+{\beta }^2}\hfill \end{array}$$

(46)

Markley’s solution has the same accuracy as QUEST. However, if we put the proposed equality $p=q$ into this solution, it is difficult to draw the conclusion that the quaternion is independent of w. The is because Markley’s solution relies on the normalized cross products ${\mathit{b}}_3$ and ${\mathit{r}}_3$. It is the nonlinearity of the normalization that makes the further simplification harder. The above results can also be computed by introducing the Lagrange multiplier [38] and in combination with the GDA [18]; as, in such an optimization, there will be two constraints: (1) quaternion unity and (2) dot product equality, and because the solving process will belong to the type of interior-point solvers for nonlinear programming [39], deriving the closed-form solution for such a problem seems much harder than the presented approach.

3.3. Error and Covariance Anlysis

Equation (1) is an ideal model for two-vector attitude determination without noises. The noise-corrupted model is expressed as follows

$$\begin{array}{c}{\mathit{b}}_1+{\mathit{e}}_1=\mathit{C}({\mathit{r}}_1+{\mathit{\epsilon }}_1)\hfill \\ {\mathit{b}}_2+{\mathit{e}}_2=\mathit{C}({\mathit{r}}_2+{\mathit{\epsilon }}_2)\hfill \end{array}$$

(47)

where ${\mathit{e}}_1,{\mathit{e}}_2,{\mathit{\epsilon }}_1,$ and ${\mathit{\epsilon }}_2$ are independent white Gaussian noise items with covariances of ${\mathsf{\Sigma }}_{{\mathit{b}}_1},{\mathsf{\Sigma }}_{{\mathit{b}}_2},{\mathsf{\Sigma }}_{{\mathit{r}}_1},$ and ${\mathsf{\Sigma }}_{{\mathit{r}}_2}$, respectively. Assume that for the true value $\mathit{K}$ we have a corresponding true quaternion $\mathit{q}$ with eigenvalue 1, such that

$$\mathit{Kq}=\mathit{q}$$

(48)

With noise perturbation, $\mathit{K}$ is interfered with by $\delta \mathit{K}$, also generating an eigenvector perturbation $\delta \mathit{q}$, say

$$\left(\mathit{K}+\delta \mathit{K}\right)\left(\mathit{q}+\delta \mathit{q}\right)=\mathit{q}+\delta \mathit{q}$$

(49)

Expanding it yields

$$\begin{array}{c}\mathit{Kq}+\delta \mathit{Kq}+\mathit{K}\delta \mathit{q}+\delta \mathit{K}\delta \mathit{q}=\mathit{q}+\delta \mathit{q}\hfill \\ \Rightarrow \delta \mathit{Kq}+\mathit{K}\delta \mathit{q}+\delta \mathit{K}\delta \mathit{q}=\delta \mathit{q}\hfill \\ \Rightarrow -\delta \mathit{Kq}=(\mathit{K}+\delta \mathit{K}-\mathit{I})\delta \mathit{q}\hfill \end{array}$$

(50)

As we also have $det(\mathit{K}+\delta \mathit{K}-\mathit{I})=0$, the quaternion perturbation is derived by

$$-{(\mathit{K}+\delta \mathit{K}-\mathit{I})}^{†}\delta \mathit{Kq}=\delta \mathit{q}$$

(51)

which produces the following quaternion difference [40],

$$\delta \mathit{q}={\left(\mathit{I}-\mathit{K}\right)}^{†}\delta \mathit{K}\mathit{q}$$

(52)

with

$$\begin{array}{c}\delta \mathit{K}=\left[\begin{array}{cc}\delta \kappa & \delta {\mathit{z}}^T\\ \delta \mathit{z}& \delta \mathit{B}+\delta {\mathit{B}}^T-\delta \kappa \mathit{I}\end{array}\right]\hfill \\ \delta \mathit{B}=w\left({\mathit{b}}_1{\mathit{\epsilon }}_1^T+{\mathit{e}}_1{\mathit{r}}_1^T\right)+\left(1-w\right)\left({\mathit{b}}_2{\mathit{\epsilon }}_2^T+{\mathit{e}}_2{\mathit{r}}_2^T\right)\hfill \\ \delta \mathit{z}=w\left({\mathit{b}}_1\times {\mathit{\epsilon }}_1^{}+{\mathit{e}}_1\times {\mathit{r}}_1^{}\right)+\left(1-w\right)\left({\mathit{b}}_2\times {\mathit{\epsilon }}_2^{}+{\mathit{e}}_2\times {\mathit{r}}_2^{}\right)\hfill \\ \delta \kappa =w\left({\mathit{\epsilon }}_1^T{\mathit{b}}_1+{\mathit{r}}_1^T{\mathit{e}}_1\right)+\left(1-w\right)\left({\mathit{\epsilon }}_2^T{\mathit{b}}_2+{\mathit{r}}_2^T{\mathit{e}}_2\right)\hfill \end{array}$$

(53)

The above error analysis is the standard approach for the Wahba’s solution using Davenport q-method. This error-analysis approach is presented to cope with the restriction of the unknown eigenvalues associated with the Davenport q-method. However, such results are also related to the selection of the weight w. Apparently, seen from Solution Equation (43), one may easily see that the obtained quaternion is free of w. This type of quaternion owns closed forms without high-order or highly nonlinear terms. Therefore for the proposed method, the error analysis can be performed easily using

$$\delta \mathit{q}=\frac{1}{\sqrt{{\alpha }^2+{\beta }^2+{\gamma }^2+N^2}}{\left(\delta \alpha ,\delta \beta ,\delta \gamma ,-\delta N\right)}^T+\delta \left(\frac{1}{\sqrt{{\alpha }^2+{\beta }^2+{\gamma }^2+N^2}}\right){\left(\alpha ,\beta ,\gamma ,-N\right)}^T$$

(54)

with

$$\delta \left(\frac{1}{\sqrt{{\alpha }^2+{\beta }^2+{\gamma }^2+N^2}}\right)=-{\left({\alpha }^2+{\beta }^2+{\gamma }^2+N^2\right)}^{-\frac{3}{2}}\left(\alpha \delta \alpha +\beta \delta \beta +\gamma \delta \gamma +N\delta N\right)$$

(55)

and $\delta \alpha ,\delta \beta ,\delta \gamma ,$ and $\delta N$ can be computed directly from Equation (43).

Then, the covariance of the derived quaternion using first-order approximation is also straightforward, given by

$${\mathsf{\Sigma }}_{\mathit{q}}=\mathit{J}\left(\begin{array}{cccc}{\mathsf{\Sigma }}_{{\mathit{b}}_1}& & & \\ & {\mathsf{\Sigma }}_{{\mathit{b}}_2}& & \\ & & {\mathsf{\Sigma }}_{{\mathit{r}}_1}& \\ & & & {\mathsf{\Sigma }}_{{\mathit{r}}_2}\end{array}\right){\mathit{J}}^T$$

(56)

where the Jacobian $\mathit{J}$ is

$$\mathit{J}=\frac{\partial \mathit{q}}{\partial {\left({\mathit{b}}_1^T,{\mathit{b}}_2^T,{\mathit{r}}_1^T,{\mathit{r}}_2^T\right)}^T}$$

(57)

4. Applications: Attitude Determination from Horizon Sensor and Another Generalized Sensor

The derived theory can be applied to those problems in which some of the reference vectors are hard to be determined. In this section, we choose the system configuration with a horizon sensor and another generalized sensor for two-vector attitude determination. The horizon sensor can be either an infrared Nadir detector [41] or a camera-based horizon finder [42]. It estimates the horizontal Euler angles of the mounting object by tracking Nadir or analyzes the geometry of the camera-captured horizon. Also, we need to note that in near-ground static applications, accelerometers can be used for reconstruction of roll and pitch information and can also be treated as a kind of horizon sensor.

There are many other types of sensors for a horizon sensor to accompany. For instance, the attitude determination by means of accelerometer and magnetometer can be used for pedestrian navigation. The infrared Nadir detector together with Sun sensor or magnetometer can also provide accurate enough attitude determination results. Take the accelerometer–magnetometer case for example, the magnetic reference vector is usually time-varying with the position of the rigid body. This is because the Earth geomagnetic field can be described by the following potential [43].

$$\mathcal{U}=A\sum _{n=1}^k\sum _{m=1}^n{\left(\frac{R_e}{A}\right)}^{n+1}{P}_n^m(cos\theta )(g_n^{m}cosm\phi +h_n^{m}sinm\phi )$$

(58)

where A denotes the distance of the object to the Earth center; $R_e$ is the Earth radius; $\theta$ and $\phi$ are co-latitude and longitude, respectively; $P_n^m$ is the associated Legendre function of degree m and order n; and $g_n^m$ and $h_n^m$ are Gaussian coefficients from global satellite geomagnetic measurements and are released in annual IGRF models. Then, the magnetic vector in the geodetic frame is the function of $A,\theta$, and $\phi$. The most common way to obtain such reference information is by looking up the table or directly computing the vector by spheric harmonic functions. However, when there is no source of position measurements (such as low-cost indoor smart wearables), it is very hard to obtain such information. A compromising way is to use the magnetic north rather than the true north provided by IGRF. Assuming that the accelerometer and magnetometer have their normalized vector measurements pairs as follows,

$$\left\{\begin{array}{c}{\mathit{A}}^b={\left(a_x,a_y,a_z\right)}^T\\ {\mathit{A}}^r={\left(0,0,1\right)}^T\end{array}\right.,\left\{\begin{array}{c}{\mathit{M}}^b={\left(m_x,m_y,m_z\right)}^T\\ {\mathit{M}}^r={\left(m_N,0,m_D\right)}^T\end{array}\right.$$

(59)

we can directly construct the magnetic reference vector by

$${\mathit{A}}^b·{\mathit{M}}^b={\mathit{A}}^r·{\mathit{M}}^r\Rightarrow m_D=a_x{m}_x+a_y{m}_y+a_z{m}_z$$

(60)

whereas another component can be determined by

$$m_N=cos\kappa =\sqrt{1-m_D^2}$$

(61)

where $\kappa$ is the local dip angle [44]. From the derived results shown in last section, the quaternion calculated now only obtains roll, pitch from accelerometer, and yaw from magnetometer, respectively. The magnetometer does not influence the estimation of other two angles and neither does the accelerometer [28].

A record of accelerometer/magnetometer data along with reference angles are logged from the 3DM-GX3-25 attitude and heading reference system (AHRS). This AHRS module is attached to the data collection computer via the serial port with baudrate of 921600. The sampling frequency is set to 500 Hz. The motion has been generated using a human-operated inertial stabilized gimbal so that the motion in all axes will be produced. The motivation of using the accelerometer and magnetometer as an experimental object is that we would like to theoretically explain why attitude estimation can be conducted in a reference-free manner as presented in previous literatures [21,28,45,46,47]. Using MATLAB r2016b, we evaluate the attitude determination results in Figure 4. The related loss functions are plotted in Figure 5. The results in Figure 5 magnify the nonintuitive errors presented in Figure 4.

In this experiment, the weights between accelerometer and magnetometer are both $0.5$. The local magnetic reference vector of FLAE is set to ${\mathit{M}}^r={(0.9800,0,-0.1989)}^T$ in the Northern Hemisphere. It is noticed from the figure that the roll and pitch angle determination has more observable differences than that of the proposed method between the true angles. This is because the roll and pitch angles used herein are not fully determined by the accelerometer but influenced by magnetometer as well. For the conventional Wahba’s problem, the two sensors both contribute to the observabilities of three Euler angles roll, pitch, and yaw. However, using the proposed method, we can separate these observabilities. Apart from the errors in the local magnetic reference vector during movement, the AHRS has been somewhat distorted by outer unknown electromagnetic disturbances. In such occasion, the non-yaw angles would contain evident differences with the reference values. If the proposed method behaves, the angles fit those from the reference instrument very well. The internal reason has been introduced before i.e., the roll and pitch angles are fully measured from the accelerometer while the yaw is only determined by the magnetometer. The loss function results also reflect this point. Coinciding with previous findings, the loss function values from the proposed method vary at the basic scale of $1\times {10}^{-27}$ which is tiny enough to be considered as zero.

While for the more general case in which there are one horizon sensor and another sensor, using the results in Equation (43), we can give the analytical solution to such attitude determination problem by

$$\begin{array}{c}\mathit{q}=\frac{1}{\sqrt{{\tilde{q}}_0^2+{\tilde{q}}_1^2+{\tilde{q}}_2^2+{\tilde{q}}_3^2}}{\left({\tilde{q}}_0,{\tilde{q}}_1,{\tilde{q}}_2,{\tilde{q}}_3\right)}^T\hfill \\ {\tilde{q}}_0=f_x\left(g_y{h}_x-g_x{h}_y\right)+f_y\left(g_z-f_z{h}_z\right)+\left(f_z^2-1\right)h_y\hfill \\ {\tilde{q}}_1=f_x\left(g_x{h}_z-g_z{h}_x\right)+f_y\left(g_z{h}_y-g_y{h}_z\right)-f_x{g}_x+f_y{g}_y+f_z\left(f_x{h}_x-f_y{h}_y\right)+\left(f_z^2-1\right)\left(1-h_z\right)\hfill \\ {\tilde{q}}_2=f_y\left[h_x\left(g_z-f_z\right)-g_x\left(h_z-1\right)\right]+f_x\left[h_y\left(g_z-f_z\right)-g_y\left(h_z-1\right)\right]\hfill \\ {\tilde{q}}_3=h_x\left(f_y{g}_y+f_z^2-1\right)-f_y{g}_x{h}_y+f_x\left(f_z{h}_z-g_z\right)\hfill \end{array}$$

(62)

where

$$\begin{array}{c}{\mathit{h}}^b={(h_x,h_y,h_z)}^T\\ {\mathit{h}}^r={(0,0,1)}^T\\ {\mathit{g}}^b={(g_x,g_y,g_z)}^T\\ {\mathit{g}}^r={(f_x,f_y,f_z)}^T\end{array}$$

(63)

and ${\mathit{h}}^b$ and ${\mathit{h}}^r$ are vector observation pair for the horizon sensor in b and r frames, respectively, whereas ${\mathit{g}}^b$ and ${\mathit{g}}^r$ represent the same meaning for another generalized sensor.

Now, we simulate a mission for a satellite. Four sensors, including a star tracker, a horizon sensor, a Sun sensor, and a magnetometer, are employed (see Figure 6). Only the horizon sensor and Sun sensor are employed to compute the attitude in this section; the other sensors are used for generating high-resolution attitude reference information. All of the sensors presented in this simulation study are transformed into the J2000 Earth-Centered Interval (ECI) frame for attitude computation. The Sun vector can be computed straightforward according to historical astronomy ephemeris [48]. Given the Julian time $T_{\mathrm{JD}}$ [49]:

$$\begin{array}{c}\hfill \begin{array}{c}{T}_{\mathrm{JD}}=367\mathrm{year}-\hfill \\ ⌊\begin{array}{c}7\left(\mathrm{year}+\frac{⌊\frac{\mathrm{month}+9}{12}⌋}{4}\right)+⌊\frac{275\mathrm{month}}{9}⌋\hfill \\ +\mathrm{day}+1721013.5+\frac{\mathrm{hour}}{24}+\frac{\mathrm{minute}}{1440}+\frac{\mathrm{second}}{86400}\hfill \end{array}⌋\hfill \end{array}\end{array}$$

(64)

where $\lfloor \cdots \rfloor$ is the downward floor operator, thus the normalized Sun vector ${\mathit{g}}^r$ can be determined by

$$\begin{array}{c}\hfill \begin{array}{c}{U}_{T1}=\frac{T_{\mathrm{JD}}-2451545}{36525}\hfill \\ M_{sun}=357.5277233+35999.05034U_{T1}\hfill \\ {\gamma }_{ecliptic}=280.4606184+36000.77005361U_{T1}+\hfill \\ 1.914666471sin\left(\frac{\pi M_{sun}}{180}\right)+0.918994643sin\left(\frac{2\pi M_{sun}}{180}\right)\hfill \\ e_{sun}=23.439291-0.0130042U_{T1}\hfill \\ s_x=cos\left(\frac{\pi {\gamma }_{ecliptic}}{180}\right)\hfill \\ s_y=cos\left(\frac{\pi e_{sun}}{180}\right)sin\left(\frac{\pi {\gamma }_{ecliptic}}{180}\right)\hfill \\ s_z=sin\left(\frac{\pi e_{sun}}{180}\right)sin\left(\frac{\pi {\gamma }_{ecliptic}}{180}\right)\hfill \\ {\mathit{g}}^r=\frac{1}{\sqrt{s_x^2+s_y^2+s_x^2}}{(s_x,s_y,s_z)}^T\hfill \end{array}\end{array}$$

(65)

where $M_{sun}$ denotes the Sun mean anomaly, ${\gamma }_{ecliptic}$ is the ecliptic longitude of the Sun, and $e_{sun}$ stands for the eccentricity of the Sun orbit.

The satellite motion is simulated by creating an orbit around the Earth with eccentricity of $0.15$ and semimajor axis of $9966.14\mathrm{km}$. The satellite attitude is generated with Nadir alignment subject to orbit normal constraint of $3\mathrm{deg}$. The spacecraft attitude has been simulated through the following equation of angular momentum,

$$\mathcal{I}\dot{\mathit{\omega }}+\dot{\mathit{h}}+\mathit{\omega }\times (\mathcal{I}\mathit{\omega }+\mathit{h})=\mathit{u}$$

(66)

where $\mathcal{I}$ denotes the inertia of spacecraft, $\mathit{\omega }$ is the angular rate, and $\mathit{h}$ stands for the angular momentum driven by the control input $\mathit{u}$. The attitude model has been constructed via the 4-th Runge–Kutta method, and the control method used herein is a Propotional-Derivative (PD)-type control algorithm. The star tracker in this simulation points exactly to the Hipparcos 2 24608 star. The outer geomagnetic environment has the perturbation based on the Olson-Pfitzer field. The IGRF model for the magnetometer measurement in the local geodetic frame is expressed as follows.

$$\begin{array}{c}{\left(B_A,B_{\theta },B_{\phi }\right)}^T=-\nabla \mathcal{U}\hfill \\ B_A={\displaystyle \sum _{n=1}^k}{\left(\frac{R_e}{A}\right)}^{n+2}(n+1){\displaystyle \sum _{m=0}^n}\left(g_n^{m}cosm\phi +h_n^{m}sinm\phi \right)P_n^m(cos\theta )\hfill \\ B_{\theta }=-{\displaystyle \sum _{n=1}^k}{\left(\frac{R_e}{A}\right)}^{n+2}{\displaystyle \sum _{m=0}^n}\left(g_n^{m}cosm\phi +h_n^{m}sinm\phi \right)\frac{\partial P_n^m(cos\theta )}{\partial \theta }\hfill \\ B_{\phi }=-\frac{1}{sin\theta }{\displaystyle \sum _{n=1}^k}{\left(\frac{R_e}{A}\right)}^{n+2}{\displaystyle \sum _{m=0}^n}m\left(-g_n^{m}sinm\phi +h_n^{m}cosm\phi \right)P_n^m(cos\theta )\hfill \end{array}$$

(67)

The Sun sensor not always work, as there may be an area of shade due to the Earth’s positioning during flight that prohibits Sun detection. The simulation is run for 1 single day from 1 March 2019 04:00:00.000 to 2 March 2019 04:00:00.000. The available solar luminance is indicated in Figure 7. The Sun sensor has accuracy of $20\mathrm{arcsec}$ and for horizon sensor that would be $1.5\mathrm{arcmin}$. The sensor noises are simulated to be additive and subject to the Gaussian distribution. When the spacecraft is in the Earth shade area, we do not compute the attitude since only one vector observation pair does not have enough observability for three-axis attitude determination.

The raw measurements of the four attitude sensors are presented in Figure 8 where the measured values of Sun sensor, star tracker, and the horizon sensor are presented in the normalized form.

Using the representative QUEST solver for conventional Wahba’s problem, we show the comparison with the proposed method. Figure 9 and Figure 10 show the results from QUEST and the proposed algorithm, respectively.

The results show that both methods are efficient for attitude determination. However, since the proposed method separates the effects of the horizon sensor and Sun sensor, the attitude error covariance has a smaller scale than QUEST. The quaternion error responses of QUEST and the proposed algorithm are shown in Figure 11. In Figure 11, the $3\sigma$ bounds are also presented using the square roots of the diagonal elements of the computed quaternion error covariances. Corresponding Euler angle errors are shown in Figure 12. The QUEST has the root-mean-squared error (RMSE) attitude accuracy of $5.12\mathrm{deg}$ on roll, $2.34\mathrm{deg}$ on pitch, and $4.06\mathrm{deg}$ on yaw. The proposed method has advantages on roll and pitch and has RMSE attitude accuracy of $2.93\mathrm{deg}$ on roll, $1.01\mathrm{deg}$ on pitch, and $3.82\mathrm{deg}$ on yaw. The proposed method can make the attitude propagation from sensors independently and thus when emergency occurs the engineers may be also easier to separate the faults and will find out the sensor failures much faster; that is, if the horizon sensor has exceptions, only roll and pitch will be affected, whereas the Sun sensor encounters faults only if yaw can be influenced.

5. Conclusions

In this paper, we revisit the problem of two-vector attitude determination, which is based on the theory of Markley [25]. A novel dot product equality constraint is proposed to further support the conventional optimization. The quaternion solutions and covariance analysis are also obtained for implementation on chips. By virtue of the equality constraint, the derived quaternion result can separate sensor effects for independent attitude angle observabilities. The derived quaternion solution is not as neat as that provided in the work by the authors of [26], but it provides another perspective without matrix orthonormalization. An application on the attitude determination from the horizon sensor and another generalized sensor illustrates that the proposed method is effective and accurate. It has also been tested smaller covariance scales in the presented simulation study for spacecraft attitude determination. Future works should be devoted to invoke more vector observation pairs for generalized attitude determination tasks. Also, according to the independent angle separation functionality of the proposed work, it may be used for flexible sensor calibration in further studies.