Optimal, Recursive and Sub-Optimal Linear Solutions to Attitude Determination from Vector Observations for GNSS/Accelerometer/Magnetometer Orientation Measurement

: The integration of the Accelerometer and Magnetometer (AM) provides continuous, stable and accurate attitude information for land-vehicle navigation without magnetic distortion and external acceleration. However, magnetic disturbance and linear acceleration strongly degrade the overall system performance. As an important complement, the Global Navigation Satellite System (GNSS) produces the heading estimates, thus it can potentially beneﬁt the AM system. Such a GNSS/AM system for attitude estimation is mathematically converted to a multi-observation vector pairs matching problem in this paper. The optimal and sub-optimal attitude determination and their time-varying recursive variants are all comprehensively investigated and discussed. The developed methods are named as the Optimal Linear Estimator of Quaternion (OLEQ), Suboptimal-OLEQ (SOLEQ) and Recursive-OLEQ (ROLEQ) for different application scenarios. The theory is established based on our previous contributions, and the multi-vector matrix multiplications are decomposed with the eigenvalue factorization. Some analytical results are proven and given, which provides the reader with a brand new viewpoint of the attitude determination and its evolution. With the derivations of the two-vector case, the n -vector case is then naturally formed. Simulations are carried out showing the advantages of the accuracy, robustness and time consumption of the proposed OLEQs, compared with representative methods. The algorithms are then implemented using the C++ programming language on the designed hardware with a GNSS module, three-axis accelerometer and three-axis magnetometer, giving an effective validation of them in real-world applications. The designed schemes have proven their fast speed and good accuracy in these veriﬁcation scenarios.


Introduction
Attitude determination (or estimation) from vector observation pairs is a significant technology in aerospace engineering and related geodetic applications [1][2][3].For instance, inertial navigation, having an important role in modern military applications, has a high demand for attitude determination accuracy for the initial alignment [4][5][6][7].A typical attitude-measuring system integrating a three-axis Accelerometer with a three-axis Magnetometer (AM) is extensively applied for real-time, continuous, stable and accurate attitude estimation for various navigation campaigns [8,9].For most applications, AM sensors are very efficient for low-cost attitude determination [10].However, the magnetometer is easily disturbed by unknown and unexpected magnetic fields in electromagnetic signal-contaminated environments.On the other hand, for large-scale region applications, the reference magnetic vector is no longer a constant vector and needs to be corrected in a timely manner with other additional heading information.Otherwise, the overall system performance will be heavily degraded.Moreover, accelerometers inevitably suffer from biases, thus leading to inaccurate attitude estimation.Therefore, auxiliary sensors are necessary to overcome such a problem.
Alternatively, the Global Navigation Satellite System (GNSS) provides precise position and velocity information [11].It has been successfully used for land, marine and aircraft attitude determination applications [12].Traditional methods integrate GNSS with inertial sensors and simultaneously estimate the orientation with synchronized position, velocity and attitude loops [13][14][15].However, this leads to the risk of degrading the accuracy of the attitude solution associated with position and velocity estimation loops.Thus, Gebre-Egziabher and Elkaim [16] proposed an independent attitude estimation loop by means of vector matching.Compared with GNSS antenna arrays, which compute highly precise solutions of baselines by using carrier-phase measurement, a single GNSS antenna is preferred in many low-cost and low-power consumption platforms.It mainly uses the simultaneous velocity vector information generated by GNSS Doppler [17,18].Indirectly, this produces important complementary orientation information, thus potentially benefiting the AM sensor system, especially for land vehicles.In addition, integrating high-rate sensors also contributes to seismogeodetic systems [19,20].
Efficient attitude estimation strategy is very crucial for a GNSS/AM integrated multi-sensor system.In essence, it can be mathematically converted into a multi-observation vector pairs matching problem.Representative methods are mainly about solutions to the famous Wahba's problem [21], which aims to obtain the optimal attitude determination results using weighted least-squares.Among these algorithms, Shuster's QUaternionESTimator(QUEST [22]), Markley's Singular Value Decomposition (SVD [23]) and Mortari's algorithm (ESOQ [24]) are the most frequently-used ones, which are mostly inspired by Davenport's q-method [25,26].Our newly-developed Fast Linear Attitude Estimator (FLAE [27]) obtains the fastest Wahba's solution, to our existing knowledge.Some other interesting approaches are proposed, as well, investigating the other aspects of this problem, e.g., Yang's analytical method, the Riemannian-manifold algorithm and Forbes' Linear Matrix Inequality (LMI) solution [28][29][30].
There are still some weight-less algorithms for multi-sensor attitude determination.They are usually used in applications like vision attitude determination, where the a priori information of weights can hardly be accurately determined.For example, using the nonlinear special orthogonal group on SO3 [31], it is able to obtain the attitude quaternion from arbitrary pairs of vectors.Via optimization approaches like the Gradient-Decent Algorithm (GDA [32]) and the Gauss-Newton Algorithm (GNA [9]), we may also compute stable solutions.Apart from these, a famous analytical method was proposed by Arun et al.where Singular Value Decomposition (SVD) is employed to calculate the attitude matrix [33].Similarly, it is then introduced in computing both the attitude and translation vector in machine vision systems [34].
For Wahba's problem, it has been shown that most existing algorithms are based on Davenport's q-method.For a long time, the attitude-solving process has been fixed in this framework, which aims to find the largest eigenvalue of the Davenport matrix K.Is it possible to seek another quite different attitude determination approach?The answer is positive, and in this paper, three novel quaternion attitude determination algorithms from pairs of vector measurements are proposed in the sense of optimal, time-recursive and sub-optimal formulations.The main contributions are listed below: 1.
The main theory is based on our previous research contributions and is extended to linear arbitrary pairs vector measurements for GNSS/AM attitude application.

3.
Simulations and real experiments are carried out, which verify the accuracy, flexibility, robustness and time consumption of various algorithms for GNSS/AM attitude determination.Detailed comparisons with representative methods demonstrate the superiority of the proposed OLEQs.
This paper is structured as follows: Section 2 briefly introduces the GNSS-AM sensor system along with its functional and stochastic models formulated by vector pair matching.Section 3 contains the problem formulation and starts with the quaternion estimation from a single sensor observation.In Section 4, the two-vector attitude determination theory along with the n-vector theory are given accordingly.Section 5 involves the numerical examples and the real field test where comparisons between the proposed SOLEQ and other representative methods are given.Finally, Section 6 consists of concluding remarks and future work.

Fundamentals of the GNSS/AM System
For attitude determination, we require GNSS receivers and AM sensor arrays for low-cost and accurate solutions.First, considering the motion behaviours of land vehicles, a vector pair for GNSS velocity can be established as: where v denotes velocity vector; the subscript G denotes the observation source 'GNSS'; the symbols r and b represent the navigation frame (r-frame, North-East-Down (NED)) and body frame (b-frame, forward-right-down), respectively; the vector is transformed from the r-frame to the b-frame; ε b G is the Gaussian white noise with variance of R G .For the land-vehicle application, it should be pointed out that the velocities in the r-frame and the b-frame are respectively: where e denotes an Earth-Centred-Earth-Fixed (ECEF) coordinate frame (e-frame), i.e., WGS-84; C r e denotes the transformation from the e-frame to the r-frame, which can be computed according to the GNSS position and Earth ellipsoid metrics in advance.
AM sensors consist of a three-axis MEMS accelerometer and three-axis magnetometer.The accelerometer gives the specific force measurement of a rigid body, and the magnetometer provides the users with the sensed local Earth geomagnetic field according to the IGRFmodel [35].
where the subscripts A and M denote the accelerometer and magnetometer sources, respectively; z denotes the observed vector; b denotes the accelerometer bias; G r = [0, 0, −g] T , where g is the gravity, is a function of the geo-location; the normalized Earth's magnetic field reference vector M r = [cosα, 0, −sinα] T , where α is the local dip angle; µ is the linear acceleration vector, which is usually treated as an external disturbance; ε is the Gaussian white noise with the variance of R. For simplicity, two points need to be clarified that: (1) the bias term has been obtained and compensated for the accelerometer readings.(2) The linear acceleration µ is estimated by using GNSS velocity information with differentiation between two adjacent epochs.Then, the vectors and matrices for the system model in the b-frame and the r-frame can be merged into the following multi-vector pair equations: where: Neglecting the cross-correlations between sensors, the stochastic model of the system is:

Problem Formulation
The conventional Wahba's problem aims to find the optimal attitude matrix from vector observation pairs in the sense of least-squares, such that: in which C denotes the optimal Direction Cosine Matrix (DCM); are the i-th pair of normalized vector observations from the body frame b and the reference frame r, respectively; a i is the weight of the i-th sensor output, which is given by the standard deviations of the input vectors from σ 1 to σ n : provided that the variance information such as shown in Equation ( 6) is predetermined [22].This ensures the biggest variance scalar with respect to the smallest weight.Wahba's problem is solved via many representative methods [36].Many of these algorithms solve the problem via eigenvalue decompositions analytically or numerically [28,29].When there is only one pair of vector observations, Wahba's solutions fail to obtain the optimal quaternion since there will be ambiguous quaternions corresponding to the maximum eigenvalue of one [37,38].In our previous contribution [39], the continuous stable quaternion solution to an accelerometer-based attitude determination system is derived.

Quaternion from a Single Sensor Observation
Considering an attitude determination model from a pair of vector observations: this section deals with the attitude determination from a single pair of sensor observations.Note that the DCM is decomposed with quaternions q = (q 0 , q 1 , q 2 , q 3 ) [39] via: C = (P 1 q, P 2 q, P 3 q) (10) in which: q 0 q 1 −q 2 −q 3 −q 3 q 2 q 1 −q 0 q 2 q 3 q 0 q 1    , P 2 =    q 3 q 2 q 1 q 0 q 0 −q 1 q 2 −q 3 −q 1 −q 0 q 3 q 2    , P 3 =    −q 2 q 3 −q 0 q q 1 q 0 q 3 q q 0 −q 1 −q 2 q    (11) In this section, the theory is extended to an arbitrary sensor with a similar approach as in [39].Inserting Equation (10) into Equation ( 9) yields: It has been proven in [39] that: where † stands for the Moore-Penrose pseudo-inverse.In fact, another property can be deduced from the following theorem: Theorem 1. P 1 P 2 + P 2 P 1 = P 1 P 3 + P 3 P 1 = P 2 P 3 + P 3 P 2 = 0 3×3 .
Proof.This can be proven via brute-force calculation: Proof.We have: which gives Lemma 1.
Hence, Equation ( 12) can be transformed into: Since K(q) is in the form of a quaternion, K (q)D b can be expanded by: Furthermore, it can be obtained that: Then, we have: where W is given by: Therefore, the attitude determination problem is shifted to: Theorem 2. For any normalized vector observation pair D r , D b , W 2 = I, where I is the four-dimensional identity matrix.
Proof.The characteristic polynomial of W is given by: Then, with the Cayley-Hamilton theorem, we substitute W for its eigenvalues and obtain: and as W 2 − I is real symmetric, it should be 0, which completes the proof.
According to our previous contribution [39], Equation ( 21) can be treated as an iterative dynamical system: where q(n) denotes the quaternion for the n-th iteration.Furthermore, as has been proven, the discrete-time system can be converted to a corresponding continuous system.If W 2 = I, the stable solution to the continuous-time system is calculated as: where q rand denotes a randomly-chosen unit quaternion.This provides us with a new approach to obtaining the measurement quaternion from a single strapdown sensor.

Optimal Linear Estimator of the Quaternion
With Equation ( 9), it is possible to list all the single rotation equations together as follows: which can in turn be converted to: This system is rewritten as: Noticing that: hence the pre-multiplication leads to: It should be then noted that Equation ( 28) is free of noises.In this case, the maximum eigenvalue of n ∑ i=1 a i W i in engineering practice is very close to the noise-free theoretical result of one.Corresponding to Equation (24), based on Brouwer's fixed-point theorem, it is possible to recursively obtain the normalized optimal quaternion by rotating a randomly given initial quaternion over and over again to infinity.This is something similar to the power method of the numerical eigenvalue factorization, but can be accelerated with geometric series-like iterations as follows: where R denotes the rotation operator over the Hamilton space H.In fact, the above iterations can hardly be achieved when the maximum eigenvalue approaches one.The reason is that at this time, the power R 2 approaches I, as well.A more robust way is shown to solve this problem.Considering both sides of Equation ( 28), we find out that the right side is in fact the mixture of solutions to single vector observation pairs.As mentioned in Equation ( 25), a stable and continuous solution to each single equation can be done by pre-multiplying 1 2 the quaternion evolutes are from the mixture of the single optimal solution.Here, the rotation operator is revised as: This equals the least-square of the set of pre-computed single rotated quaternions, which is definitely faster and more robust than rotation from a randomly given initial quaternion.
The covariance of the calculated optimal quaternion is equal to that of QUEST under the condition of least-square optimality: where Q = (q 1 , q 2 , q 3 ) is the vector part of the quaternion.
[q] defines the following matrix:

Variant 1: Recursive-OLEQ
We have seen from the above formulations that for each epoch, the vector observations are batch processed thoroughly.When used in aerospace electronic systems, the measured vector observation pairs in neighbouring time epochs are usually continuous since they are always smoothed by the sum filters and Low-Pass Filters (LPF).Therefore, with this consideration, the attitude quaternion can be propagated from the last estimated one using the rotation operator described before.In this way, the quaternions are recursively computed with much fewer computations, and the accuracy is maintained.Furthermore, for highly-reliable attitude determination systems, high-precision gyroscopes are usually employed.This provides us with a second-stage acceleration scheme, inspired by the conventional recursive algorithms like filter QUEST, REQUEST, etc. [40], by which we may first rotate the estimated quaternion in the last time epoch with the zero-order angular transition matrix by: where: in which ∆t is the sampling time and [Ω×] composed by the angular rate from the gyroscope, which is detailed in many classical literature works [41].After this, even a single rotation by R would then be very accurate.The one-step covariance matrix of the obtained quaternion is calculated by: 4.2.Variant 2: SOLEQ

Two-Vector Case
When there are two pairs of vector observations, regardless of the weights of the respective sensors, it is natural that the attitude quaternion may be computed via: where W 1 and W 2 are for the first and second sensors.Equation ( 25) is reasonable since we have: However, for the two-vector case, one can write: This reflects that the accurate quaternion can be recursively computed with: which exits when the Euler distance q k − q k−1 is less than one predetermined threshold η.Note that the above initialization procedure is equivalent to the following process: where q rand is a randomly chosen quaternion; j is chosen to make sure q j − q j−1 < η.Let us define the following transformation operators: Theorem 3.For the two-vector attitude determination case, the steady-state evolution in Equation ( 43) is not affected by the transformation operators' order, such that: Proof.The integrated transformation can be computed by: and accordingly These generate: and: Finally, we have: This proves that the mixed steady-state transformation (AB) j can be achieved by independent transformations from A or B, which completes the proof.
Following this theorem, the problem is to compute the power A j .In fact, A is formed with 1 2 (W 1 + I) and 1  2 (W 2 + I).Their respective eigenvalue decompositions can be given by: where V and D are constituted by eigenvectors and eigenvalues respectively as (W i + I) is real symmetric [42].Since W 1 and W 2 are in the same form, the eigenvalue matrices are equal to each other, i.e., Then, A can be rewritten as: Identically, we have: Combining Equations ( 53) and (54), it is obtained that: Note that: Equation ( 55) is simplified as: Here, we define: Actually, it can be decomposed by: An interesting fact is that the eigenvalue matrix S can be analytically calculated and is given by: where diag(•) represents the diagonal matrix.This further yields H to be a matrix with the form of: Then, U is computed by: where u 2 = u 3 .Using this, we have: U can be decomposed with eigenvalue decomposition as well, such that: which can be analytically given by: with: Letting: U j is finally computed by: Therefore: The required computation of V i is computed as follows: where Ṽi (x, y) stands for the element of Ṽi in the x-th row and the y-th column, the details of which are given by Equation (71): in which the factor is computed by: Related information can also be acquired from [27].It should be noted that: Thus, the Gram-Schmidt orthogonalization is applied to Ṽi , enabling V i V i = V i V i = I [42].A typical commitment to achieve this is to implement the QR decomposition [43], such that: where R denotes an invertible upper triangular matrix.If q rand = (1, 0, 0, 0), the suboptimal quaternion is equal to the normalized first column of (AB) j .

n-Vector Case
Corresponding to the above notations and derivations, the n-vector case's transformation operators are defined by: Defining: we have: Then: Accordingly, the normalized first column of (AB) j constitutes the attitude quaternion.

The Effect of Power Order
In this sub-section, we will show that the selection of j is in fact not influential with respect to the final result at all.Letting: we have: Here, we also have: and: Therefore, with increasing iteration numbers, the item multiplied by λ j U,1 gradually vanishes in the results.The limiting result of AB j turns out to be: Additionally, the quaternion solution is not about which column of AB j , and the result is the normalization of the following vector:

Discussion of OLEQs
The three derived OLEQs can be used on different occasions.The OLEQ incorporates the weights so that the determination results are optimal in the sense of least-squares.With the aid of the gyroscope, the ROLEQ can achieve faster and smoother estimates.The meaning of the proposed SOLEQ is that it has a very simple linear expression that may generate short and concise analytic results for certain sensor combinations.Furthermore, when required in applications where the weights can hardly be accurately determined, e.g., vision attitude determination, the SOLEQ could be an alternative choice, as well.The attitude determination results of the three OLEQs are evaluated in the following experimental section.

Simulation: Common Cases
In this sub-section, the sensor observations are simulated with random reference vectors and true DCM in: where ε is the noise item, which are not supposed to have any correlation with each other and be subject to Gaussian distribution.The reference vectors and the standard deviations of noise items are selected according to the classical test samples by Markley (see Table 1), where the reference DCM C true is:

Case Reference Vectors
Noise Standard Deviations Using the simulated samples, the mean attitude Root Mean Squared Errors (RMSEs) in Euler angles are evaluated with our proposed algorithms OLEQ, SOLEQ and representative algorithms including QUEST and FLAE [27], which are shown in Tables 2-4.Table 5 contains the computed average Wahba's loss function values for different cases and algorithms.These algorithms are executed with MATLAB r2016 software (r2016, MathWorks, Natick, USA) on a PC 10,000 times with each data sample.We first observe the attitude RMSEs.From the results of OLEQ, QUEST and FLAE, it is noticeable that they have a similar attitude determination accuracy.Combining the same statistics in Table 5, the proposed OLEQ is well verified for its optimality.From the presented results, we see that SOLEQ has larger attitude errors and loss function values compared to the other optimal methods.The proposed SOLEQ is sub-optimal as it actually approximates the attitude estimator where the weights are ignored.Therefore, this simulation scenario has validated the correctness and optimality of the proposed OLEQ and SOLEQ.

Simulation: Extreme Cases
Conventional Wahba's solutions face a dilemma when exposed to some critical pairs of vector observations.For instance, Markley and Mortari give an example where the sensors are configured by Equation (88) [37].In such a scenario, the root of the characteristic polynomial of the Davenport matrix cannot be easily obtained by Newton iterations.The internal reason is given by Cheng [44] showing that it results in numerical loss according to the specific CPU word storage length.A flexible transformation of the characteristic polynomial is proposed then to significantly boost the convergence.As such configurations indeed happen in engineering practice, there is a necessity to evaluate the proposed schemes by comparisons with representative solvers.With similar simulation techniques as previously mentioned, the vectors are simulated with given reference vectors and standard deviations by rotation of C true .Here, the QUEST algorithm is revised to Cheng's form.First, we mainly compare the two iterative methods QUEST and OLEQ because in our existing paper [27], QUEST and FLAE have been proven to have a very similar behaviour when faced with this extreme case.Here, the iteration stops when the Euclidean norm of the neighbouring attitude quaternion difference is less than 1 × 10 −8 .For QUEST, the maximum iteration number is set to 50.The obtained results are depicted in Figure 1.We can see that the supervised QUEST can obtain accurate quaternion solutions within several iterations.Actually, before Cheng's improvement, QUEST may exceed the maximum iteration number from time to time.The proposed OLEQ, however, shows better performance dealing with this extreme case.Furthermore, the final mean loss function values of the two algorithms are computed as 4.9890 × 10 −11 and 2.8391 × 10 −10 , which reveals that the proposed OLEQ can not only obtain faster solutions, but leads to smaller loss function values, compared with supervised QUEST.As is known to everyone, QUEST is the most representative Wahba's solution using Davenport's q-method.Many other algorithms like ESOQand FOAMactually have the same performance as QUEST.Therefore, in this way, OLEQ is proven to be faster and much more robust than the whole class of algorithms based on Davenport's q-method.This also shows that the presented novel attitude evolution method owns brand new abilities.

Experiment: Accelerometer-Magnetometer Case
In this sub-section, we conduct an experiment where the accelerometer-magnetometer combination is adopted.Such a sensor combination is extensively applied nowadays in low-cost attitude estimation schemes.The accelerometer is pre-calibrated using six-direction bias cancelling, while the magnetometer is calibrated using the method proposed by Y .Wu et al. [45].
The hardware is constituted by a battery, a high-end Attitude and Heading Reference System (AHRS) with a high precision internal accelerometer, magnetometer and gyroscope, a transmitter for remote data transmission and a micro-controller for the implementation of the algorithm using the C++ language on the FreeRTOS.With the designed hardware platform shown in Figure 2, we collect a dataset with 10,000 samples.
The main purpose of this sub-section is to validate the performances of the proposed OLEQ, SOLEQ and ROLEQ since the AHRS has angular rate readouts.The compared results with the reference angles from representative methods are obtained (see Figures 3-5).Note that here, the weights between the accelerometer and magnetometer for Wahba's solution are chosen as 0.63 and 0.37 according to their respective noise characteristics.Yet, the local magnetometer's reference vector is calculated as M r = (0.60311, 0, −0.79766) in Wuxi, Jiangsu Province, China.In Figure 3, the reference angles, QUEST solutions and SOLEQ solutions have been presented.The results indicate that the proposed sub-optimal estimator can estimate the attitude angles with similar macroscopic accuracy.Detailed attitude errors are plotted in Figures 4 and 5.It is noticed that in general, these algorithms have the same errors with respect to the reference.Figure 5 shows that the overall attitude accuracy of ROLEQ is slightly smaller than the others.This is because it is first processed using the angular rate data, which essentially plays the role of a smoothing procedure.results using QUEST and SOLEQ.

Experiment: GNSS-Aided Attitude Determination for Land Vehicles
The GNSS receiver is widely employed in the attitude determination tasks for land and unmanned aerial vehicles.In this experiment, we use a designed rover (see Figure 6) to validate the feasibility of the proposed algorithm for GNSS attitude determination.
The rover is equipped with the aforementioned navigation computer and employs an external ublox M8N GNSS module with a serial-com connection to the board at the sampling frequency of 5 Hz.This rover is controlled by a handheld 2.4-GHz transmitter, and the onboard Pixhawk autopilot generates PID controlling commands to the servos and motors according to internal measurements.In this experiment, the rover is ran on a playground of UESTC, and we pick up one period of data in which the GNSS velocity is valid.In the data history, the magnetometer was distorted by outer unknown electromagnetic and ferromagnetic fields.Furthermore, during the execution process, raw measurements from the gyroscope and accelerometer are also logged with a speed of 250 Hz.The raw data are shown in Figure 7.
According to the sensor noise characteristics, the weights of the accelerometer-magnetometer combination are 0.9 and 0.1, respectively; while for the accelerometer-magnetometer-GNSS combination, the weights are chosen as 0.474, 0.05 and 0.474, respectively.The reference vector of the magnetometer is determined by the initial GNSS positioning results with the IGRF model [35].By making use of algorithms including QUEST, FLAE, OLEQ and SOLEQ, the computation results are summarized in Figures 8 and 9.
We specifically add the GNSS velocity norm under the Euler angle results to illustrate the influence of the velocity scalar on the attitude determination results.In principle, when the vehicle is not moving with relative discriminative velocities, the GNSS receiver cannot give accurate speed estimates.Therefore, it is shown in the initial stage of the attitude determination results where the GNSS takes part, in that the determination accuracy of the yaw angle is seemingly very poor.As the velocity increases, the accuracy is improved very quickly accordingly.The accelerometer-magnetometer combination is largely distorted by magnetic disturbances.The integrated results of roll, pitch and, moreover, the yaw angles are influenced, generating very obvious differences with reference angles.With the aid of GNSS velocity, the corresponding attitude determination accuracy is not damaged because Wahba's solution balances the sensor inputs by the weights.It is observed that the OLEQ is validated to have almost the same accuracy for normal sensors in the aforementioned section, and in this section, such behaviour holds, as well.The SOLEQ, however, does not employ the weights and therefore produces relatively worse estimates, but for GNSS case, it is still better than the accelerometer-magnetometer estimates.The RMSEs of these scenarios are shown in Table 6.The results demonstrate to us the validity of the proposed algorithms, especially OLEQ.

Computation Time
From another point of view, the time consumption of various algorithms needs to be investigated.The time consumption is assessed on the embedded platform with the C++ programming language to ensure fairness.A rough evaluation is done with two pairs of vector observations in a few samples, which shows direct time consumption results (see Figure 10).As shown in the figure, for two pairs of vector observations, the three proposed algorithms' computation times are between QUEST and FLAE.However, from the expressions of the algorithms presented before, the number of vector observations is influential with respect to the final time consumption.Hence, with the simulation samples, each algorithm is again tested 20,000-times with different vector observation numbers.The time consumption is averaged and plotted in Figure 11.The results indicate that the algorithms are all linear, owning time complexities of O(n).QUEST, OLEQ and ROLEQ join at the vector observation number of 20.For common tasks, such a number covers most amounts of sensors.Although FLAE has the least time consumption, it can not overcome the drawbacks of extreme cases as well as OLEQ.That is to say, the proposed algorithms can replace the original algorithms for faster and more robust attitude determination.

Conclusions
This paper revisits the attitude determination from vector observations for a GNSS/AM case study.Novel linear algorithms are designed to obtain accurate attitude estimates in the sense of least-squares.Handling it in this manner, the computed quaternion is identical or suboptimal with respect to the conventional Wahba's solutions including QUEST and FLAE.We have found out that: 1.
Numerical simulations exhibit that the proposed OLEQs have a similar accuracy with representative solvers.

2.
It is also demonstrated that faced with extreme cases, OLEQs show much more robustness with less computation iterations.

3.
The computation speeds of OLEQs are tested, revealing that they belong to computationally-efficient algorithms.4.
Moreover, a real vehicular experiment of a GNSS/AM system is designed and conducted showing the effectiveness of the proposed OLEQs in real-world embedded applications.
The MATLAB source codes of the proposed algorithms are uploaded as open-source files at https://github.com/zarathustr/OLEQs.The presented approach provides the reader with a brand new viewpoint of attitude evolution and hopefully will benefit related multi-sensor attitude determination applications.

Figure 1 .
Figure 1.Iteration numbers of QUEST and OLEQ when faced with an extreme case.

Figure 2 .
Figure 2. Designed hardware for the implementation of the algorithm.

Figure 4 .
Figure 4. Experiment results of OLEQ, SOLEQ and representative algorithms using sampled data and different algorithms.

Figure 5 .
Figure 5. Experiment results of Recursive OLEQ (ROLEQ) and representative algorithms using sampled data and different algorithms.

Figure 6 .
Figure 6.The designed multi-functional rover for the validation of the proposed algorithms.

Figure 7 .
Figure 7. Raw sensor measurements from the memory logging.

Figure 8 .
Figure 8. Attitude determination from the accelerometer, magnetometer and velocity output of the GNSS receiver, by means of QUEST and FLAE.

Figure 9 .
Figure 9. Attitude determination from the accelerometer, magnetometer and velocity output of the GNSS receiver, by means of OLEQ and SOLEQ.

Figure 10 .
Figure 10.Time consumption of various algorithms.

Figure 11 .
Figure 11.Time consumption of algorithms with respect to the number of vector observations.

Table 6 .
RMSEs of the experiment (in deg).