Estimation of Quaternion Motion for GPS-based Attitude Determination Using the Extended Kalman Filter

In this paper, the Global Positioning System (GPS) interferometer provides the preliminarily computed quaternions, which are then employed as the measurement of the extended Kalman filter (EKF) for the attitude determination system. The estimated quaternion elements from the EKF output with noticeably improved precision can be converted to the Euler angles for navigation applications. The aim of the study is twofold. Firstly, the GPS-based computed quaternion vector is utilized to avoid the singularity problem. Secondly, the quaternion estimator based on the EKF is adopted to improve the estimation accuracy. Determination of the unknown baseline vector between the antennas sits at the heart of GPS-based attitude determination. Although utilization of the carrier phase observables enables the relative positioning to achieve centimeter level accuracy, however, the quaternion elements derived from the GPS interferometer are inherently noisy. This is due to the fact that the baseline vectors estimated by the least-squares method are based on the raw double-differenced measurements. Construction of the transformation matrix is accessible according to the estimate of baseline vectors and thereafter the computed quaternion elements can be derived. Using the Euler angle method, the process becomes meaningless when the angles are at 90 where the singularity problem occurs. A good alternative is the quaternion approach, which possesses advantages over the equivalent Euler angle based transformation since they apply to all attitudes. Simulation results on the attitude estimation performance based on the proposed method will be presented and compared to the conventional method. The results presented in this paper elucidate the superiority of proposed algorithm.


Introduction
The Global Positioning System (GPS) [1][2][3] is a satellite-based navigation system that provides a user with the proper equipment access to useful and accurate positioning information anywhere on the globe. In addition to the code observable commonly used for position and velocity determination, the carrier phase is the other type of observable that can be extracted from the GPS signals. The carrier phase observables contain noise much smaller than that of code observables. Due to its higher accuracy and precision compared to code observations, carrier phase observations can be used for relative positioning in centimeter level and has been widely applied to surveying, attitude determination [4][5][6][7][8][9][10][11][12][13][14], precision approach and automatic landing.
Traditionally used in precise static relative positioning, and thereafter in kinematic positioning, the GPS offers the interferometer [5] for attitude determination by processing the carrier phase observables, which enables the relative positioning in centimeter level. The relative positioning techniques using the carrier phase differential GPS (DGPS) based on interferometer principles can be adopted to solve for the baseline vectors, defined as the vectors between the antenna, designated master and one of the slave antennas, shown as in Fig. 1. The attitude of a vehicle can be precisely determined using the GPS carrier phase observables received at two or more nearby antennas attached to the vehicle. In the beginning of 1990s, Van Grass and Braasch [4] conducted research on the GPS to the field of aircraft attitude determination using carrier phase. In their work, the receiver-satellite double differenced observable was employed. The solution of the baseline vector is the approximate interferometer coordinates, which directly influence the performance of the GPS-based attitude determination. Real-time integer ambiguity resolution techniques [5,8,9,11] and attitude determination are two main issues to be resolved for determining the vehicle attitude when applying GPS double-differenced carrier phase. Very accurate relative position estimate will be available once the integer cycle ambiguities are properly resolved. Attitude determination using GPS does not have the error accumulation, which usually happens in the inertial navigation system (INS) [2].
The rotation angles that relate a coordinate system fixed in the body (body frame) to a coordinate system fixed in space are referred to as the attitudes. The space coordinate system is typically defined to be a local level NED (north-east-down) frame, also referred to as the navigation frame. The purpose of attitude determination essentially involves calculation of the three Euler angles, namely roll, pitch, and yaw. The quaternion method [15][16][17][18][19] uses four parameters instead of nine as in the Euler angle method, by defining the generalized complex number. The quaternion method possesses advantages over the equivalent Euler angle based transformation since they apply to all attitudes with the error equations bounded by the constraint equation. Furthermore, the numerical value of each parameter always lies within the range of −1 to 1, so that the scaling problems in the computing mechanization can be easily handled.
The purpose of the Kalman filter (KF) [2,20,21] is to provide an optimal (minimum mean-square error) estimate of the system state vector. As the nonlinear version of KF, the extended Kalman filter (EKF) deals with the case governed by the nonlinear stochastic differential equations. The EKF linearizes about a trajectory that is continually updated with the state estimates. Both the KF and EKF have been widely applied to the field of navigation, such as GPS receiver position and velocity determination, attitude Figure 1: Interferometric configuration considering the plane wave approximation of a GPS interferometer [5] determination, integrated navigation system design, and the carrier-smoothed-code (CSC) processing. Utilization of the EKF as the estimator of attitude related parameters enhances the accuracy and reliability of the attitude solution.
The remainder of this paper is organized as follows. In Section 2, preliminary background on the GPS carrier phase observation model is reviewed; the computed quaternion vectors based on the GPS interferometer is introduced. Section 3 presents the modelling of quaternion dynamics for the extended Kalman filter. In Section 4, simulation experiments are carried out to evaluate the performance on estimation accuracy using the proposed method as compared to the conventional one. Two numerical examples are presented for illustration. Conclusions are given in Section 5.

GPS-Based Computed Quaternion Vectors Based on the GPS Interferometer
In a GPS interferometer, the receiver-satellite double-differenced carrier phase observable has been commonly utilized to solve for the antenna baseline vector. The GPS carrier phase observables representing sum of range, an unknown integer ambiguity and some ranging errors -can be represented by: where the parameters involved in Eq. (1) are defined as follows. r: True range between a satellite and receiver; c: Speed of light; dt: offset of the satellite clock from GPS time; dT : Offset of the receiver clock from GPS time; d ion : Ionospheric error; d trop : Tropospheric error; : Carrier phase wavelength; N : Carrier phase integer ambiguity; v q , v ' : Measurement noises of code and carrier phases, respectively.

Formulation of the Transformation Using the Baseline Vectors
The receiver-satellite double differencing operator is defined asrDðÁÞ ¼ DrðÁÞ ¼ rðÁÞ À DðÁÞ, where DðÁÞ ¼ ðÁÞ 1 À ðÁÞ 2 denotes the between receiver single differencing operator for receivers 1 and 2, and rðÁÞ ¼ ðÁÞ i À ðÁÞ j denotes the between satellites single differencing operator for satellites i and j. Referring to the configuration as in Fig. 2 [4,14], when using the carrier phase signal from satellite i, the between-receiver single-differenced observable is a linear combination of two phase observables received by two antennas where the effects of errors associated with the satellites are greatly reduced. Similarly, for satellite j, we have where b is the baseline vector formed by two antennas, and e represents the line-of-sight unit vector from antennas to a satellite. Taking two independent single-differenced observables leads to the receiversatellite double-differenced observable: which eliminates or greatly reduces the satellite and receiver timing errors.
The signals received from n satellites by one GPS interferometer provide n À 1 independent double differences. When the integer ambiguity parameter (rDN ij 12 ) is resolved, the range-based equivalent of Eq. (4) is depicted as follows: which can be expressed in matrix formrDr ¼ Gb. The baseline vector can be obtained by the least-squares approachb ¼ ðG T GÞ À1 G T rDr. The solution of the baseline vector b ¼ ½ b x b y b z T is the approximate interferometer coordinates, which directly influences the performance of the GPS-based attitude determination.
Referring to Fig. 3, there are two body-frame-mounted non-collinear baseline vectors formed: b 1 ¼ S A À M and b 2 ¼ S B À M, where the master antenna (M) position is located at ½ 0 0 0 T in the body frame, while the other two slave antennas S A and S B are at ½ 1 0 0 T ' and ½ cos c sin c 0 T ', respectively. The parameter ' is the baseline length parameter used to adjust the length, and c is the angles between two baseline vectors which are adjustable for design flexibility. The accuracy of the attitude measurement depends on the baseline to noise ratio, and is also a function of antenna placement and GPS satellite geometry. The coordinate transformation matrix from body to local frameR b2n can be formed once the baseline vector is determined through the following calculation.
Since the conventional least-squares approach for baseline vector estimation is inherently noisy, incorporation of the Kalman filter into the GPS interferometer for performance improvement is accessible.

Transformation Matrix Involving the Euler Angles
In certain applications, the body angular velocities: p, q, r(which can be measured, for example, by body mounted rate gyros) information is required. They can be determined from the Euler rates: _ f, _ h, _ w. The relationship between the body angular velocities and the Euler rates can be written as: Once the Euler angle rates are available with sufficiently good accuracy, the body angular velocity can be obtained. Taking the inverse transformation for Eq. (7), we have the Euler rates in terms of the body angular velocities: The direction cosine matrix (DCM) R b2n is employed as a transition matrix to describe the transformation of a vector quantity defined in the body frame (denoted by 'b') to the geodetic frame (denoted by 'n'), f n ¼ R b2n f b . There are two common types of transformation approaches available for solving vehicle attitudes, typically including Euler angle method and quaternion method. The transformation matrix related b frame relative to n frame can be constructed in terms of the Euler's angles or the quaternion parameters.
where the subscripts n and b represent the local and body frames, respectively. Since R b2n is an orthonormal matrix, its inverse can obtained through its transpose: In Eq. (9), the notations S ðÁÞ sinðÁÞ and C ðÁÞ cosðÁÞ are defined. Since the vehicle attitude is defined by the angles between the NED frame and body frame, therefore, the rotation transformation matrix that relates the body and NED frames provides the information for finding the vehicle attitude [2,4]. The vehicle attitude can be obtained through the calculation: Figure 3: Antenna differential position vector geometry Using the Euler angle method, the singularity problem occurs when the angles are at 90 , the process becomes meaningless.

Transformation Matrix Involving the Quaternion Elements
A quaternion is a four-dimensional extension to complex numbers, containing four real parameters. The first is considered a scalar and the other three vector components in three-dimensional space: where a constant equation exists of the form q 2 1 þ q 2 2 þ q 2 3 þ q 2 4 ¼ 1. The transformation matrix from body frame to navigation frame axes in terms of the unit quaternion parameters is given by where q 1~q4 are the quaternion vector components. The relations between the quaternion q and the two vectors f b and f n are nonlinear.
Using the transformation matrix to obtain quaternion parameters can be done through By comparing Eq. (13) with Eq. (9), conversion of the quaternion parameters to Euler angles can be implemented through the following relationships: See [2] for more details on the quaternion method. The proposed algorithm for implementing the computation of the quaternion vector derived from the baseline vectors based on the GPS interferometer to be employed as the measurement of the EKF is provided in Fig. 4.
The interferometer offers the preliminarily computed quaternion vector using the GPS doubledifferenced carrier phase observables. The implementation procedure is highlighted as the following steps: (1) Determining the baseline vector using the receiver-satellite double-differenced carrier phase observables:b e 1 andb e 2 from the GPS interferometer; (2) Construct the transformation matrix according to Eq. (6); (3) Compute the quaternion elements.

Modelling of the Quaternion Dynamics for the Extended Kalman Filter
Consider a dynamical system whose state is described by a nonlinear, vector differential equation. The process model and measurement model are represented as It is assumed that f and h are known functions, uðtÞ and vðtÞ are two white-noise processes mutually independent with zero means and: where dðt À sÞ is the Dirac delta function, E ½Á represents expectation, and superscript "T" denotes matrix transpose. The nonlinear filtering deals with the case governed by the nonlinear stochastic difference equations. Assuming the process to be estimated and the associated measurement relationship may be written in the form: where the state vector x k 2 < n , process noise vector w k 2 < n , measurement vector z k 2 < m , and measurement noise vector v k 2 < m . The vectors w k and v k are zero mean Gaussian white sequences having zero cross-correlation with each other: where Q k is the process noise covariance matrix, R k is the measurement noise covariance matrix. The symbol d ik stands for the Kronecker delta function:

The Extended Kalman Filter
The nonlinear process model can be linearized along the currently estimated trajectory where the influence of the perturbations on the trajectory can be represented by a Taylor series expansion about the nominal trajectory.
The actual trajectory is the sum of the estimated trajectory and the small increment: x ¼x þ dx. The corresponding difference equation by converting the continuous-time model into a discrete-time model is given by When working with incremental state variables, the linearized measurement equation presented to the EKF is dz k ¼ z k À hðx À k ; kÞ ¼ H k dx k þ v k rather than the total measurement (nonlinear) z k . Consider the incremental estimate update equation at step k in which the measurement residual is given by: z k À hðx À k ; kÞ À H k dx À k , and the predictive measurement is the sum of hðx À k ; kÞ and H k dx À k . The residual involves the noisy measurement minus the predictive measurement based on the corrected trajectory rather than the nominal one.
Addingx À k on both sides of Eq. (20)x À k þ dx k ¼x À k þ dx À k þ K k ½z k Àẑ À k , we have the update equation: which shows that in an EKF it is accessible to keep track of the total estimates rather than the incremental ones. Oncex À k is determined, the predictive measurementẑ À k can be formed as hðx À k ; kÞ, and the measurement residual at k þ 1 is formed as the difference ðz k Àẑ À k Þ. Projection ofx k tox À kþ1 can be done through the nonlinear dynamicsx À kþ1 ¼ fðx k ; kÞ. Without the external aiding such as an inertial sensor to provide a reference trajectory, the process dynamics represent the total observer dynamics, as shown in Fig. 5. Implementation algorithm for the discrete EKF equations is provided in Tab. 1.

Reference
Estimate of states Figure 5: Configuration of the EKF

The Extended Kalman Filter for Quaternion Estimation
It can be shown that the quaternion elements are propagated according to the differential equation: The differential equations for describing the propagation of quaternion elements are given by (1) Evaluate the predicted state vector through the process model x À kþ1 ¼ fðx k ; kÞ (2) Propagate the predicted error covariance matrix kÞ, and the linear approximation equations for system and measurement matrices are obtained through the relations where the product terms in the parentheses are introduced by quaternion product between the angular rate and the quaternion, which can be written in matrix form Therefore, the quaternion vector is propagated according to the differential equation where q ¼ ½ q 1 q 2 q 3 q 4 T denotes the quaternion vector, and x b nb ¼ ½ p q r T describing the vector of body angular velocities. The symbol Ωðx b nb Þ ¼ ðx b nb ÂÞ represents the skew symmetric matrix with components of x nb in the body frame: The differential equations for describing the quaternion vector can be represented by where Wðx b nb Þ ¼ 0 Àp Àq Àr p 0 r Àq q Àr 0 p r q Àp 0 and Eq. (26) may also be written as Augmented by the propagation of the body angular rates, described by the random walk models, i.e., _ p ¼ u p ; _ q ¼ u q ; _ r ¼ u r , the differential equation has the form: _ x b nb ¼ ½ _ p _ q _ r T ¼ ½ u p u q u r T . The resulting state vector consists of seven states, in which the first four components are the quaternion elements, and the other three components are the body angular velocities. The corresponding Jacobian matrix can be shown to be The corresponding discrete state transition matrix is given by Φ k ¼ £ À1 ½ðsI À FÞ À1 ¼ expðFDtÞ ffi I þ FDt, where Dt is the step size.
Improved accuracy is accessible for the attitude solutions based on the EKF using the preliminarily computed quaternion vector from the GPS interferometer as the measurement. The measurement equations in this case are linear and are much simpler, and they are the four quaternion components i.e., where v i is the white noise measurement. The measurement model written in matrix form is given by With the GPS-based computed quaternions available as the measurement, the measurement matrix H k and noise v k can be expressed as respectively. The measurement model H k is a 4 Â 7 matrix. Although the measurement equations are linear, an EKF is still required since the process model is nonlinear. The attitude estimation is implemented based on the block diagram shown as in Fig. 6. Alternative models for body angular velocities in the process model.
or the first-order Gauss-Markov process

Illustrative Examples and Performance Evaluation
Simulation study has been carried out to evaluate the performance of the proposed approach in comparison with the conventional method for GPS-based attitude determination. Two illustrative examples were implemented through numerical experiments. Computer codes were developed using the Matlab® software. The commercial software Satellite Navigation (SatNav) Toolbox by GPSoft LLC (2003) [22] was employed for generation of the GPS satellite orbits/positions and thereafter, the satellite pseudorange, and carrier phase observables, required for simulation. Furthermore, the Inertial Navigation System Toolbox (2007) [23] was employed for generation of the Euler attitude angles of the threedimensional flight.
The error sources corrupting the GPS carrier phase observables include ionospheric error, tropospheric delay, receiver thermal noise and multipath errors. The variances of the receiver noise are assumed to be 1 m 2 . Most of the receiver-independent common errors can be corrected through the differential correction while the multipath and receiver thermal noise cannot be eliminated. The antenna geometry is set up as in Fig. 2, with the baseline length variable ' equal to 1 meter and c 90 degrees. ; where m ¼ 4 is the number of measurements and equals the number of quaternion elements. After resolving the integer ambiguity, the least-squares approach is utilized for constructing the baseline vector. Fig. 7 presents the body angular rates p, q and r for Example 1. There were 8 visible satellites available in the clear open sky during the time period of simulation. Fig. 8 provides the skyplot at the final time.
The proposed method in which Kalman filtering is used provides estimate of the quaternion vector with noticeable accuracy improvement. It should be noticed that it is risky to set the process noise variance zero to avoid filter divergence. As a result, when the system reaches steady state, the steady-state gain will not approach zero and, subsequently, the filter is able to catch the time-varying attitude dynamics. One still needs to find the suitable values to meet the specific design/mission requirement. The estimation accuracy To confirm the correctness of the solutions, the estimated quaternion vectors were examined, as shown in Fig. 9. Comparison of the attitude solutions is shown in Fig. 10; Tab. 2 summarizes the error statistic of attitude solutions for both the conventional and proposed methods. It can be seen that estimation accuracy by the proposed method has been remarkably improved.
(b) Illustrative Example 2-the three-dimensional flight. The scenario for the second example is a threedimensional flight. The 3D flight path is generated with the Inertial Navigation System toolbox, shown as in Fig. 11. Tab. 3 provides the description of the motion for the flight path. The body angular velocities as time progresses for this example are depicted in Fig. 12. In this example, the duration for simulation is 110 seconds with 8 visible satellites. The covariance matrices for the process and noise models, respectively, are given by 1e À 5 1e À 5 1e À 3 0:1 0:1 10 The standard EKF is sensitive to the changes of the process and measurement models, thus yielding poor performance. Furthermore, the EKF framework does not possess capability to deal with the non-Gaussian measurement errors/outliers. The EKF associates the contaminated measurements with an increase in the measurement covariance, causing the reformulated error covariance. This modified covariance inflation is known to cause an increase of error in the state estimation. Performance based on the EKF will degrade when the noise strength is varying and/or the actual distribution deviates from the assumed Gaussian model. To further improve the performance of the EKF, an adaptive mechanism or a robust technique can be incorporated for performance improvement. Due to appropriate tuning, the adaptive EKF (AEKF) exhibits robust behavior and therefore outperforms the standard EKF when the time-varying dynamic process and measurement models are involved.

Conclusions
A novel GPS-based attitude determination method has been presented. The quaternion vector derived from GPS interferometer has been used as the measurement of the EKF in the attitude determination system. Although utilization of the carrier phase observables enables the relative positioning to achieve centimeter level accuracy, nevertheless, the quaternion elements derived from the GPS interferometer are inherently noisy. This is essentially due to the fact that the baseline vectors estimated by the least-squares method are based on the noisy double-differenced measurements.
The preliminarily computed quaternion elements from the GPS interferometer is employed as the measurement of the EKF based quaternion estimator. The model based approach using the EKF is adopted for estimating the quaternion elements, which can then be converted to the Euler angles. Results show that, by incorporating the EKF into the GPS interferometer, the errors in the solutions of the baseline vectors, and thereafter the quaternion elements have been remarkably mitigated and the estimation accuracy of the attitude solutions has been noticeably improved. The proposed method provides several advantages, such as accuracy improvement, reliability enhancement, and real-time characteristics.
Since the EKF is sensitive to the changes of the process and measurement models, when implementing the EKF approach, poor knowledge of the noise statistics may seriously degrade the estimation performance, and even provoke the filter divergence. For achieving improved estimation accuracy, the designers are required to possess the complete a priori knowledge on both the dynamic process and measurement models. In the two illustrated examples, both the process and measurement noise parameters remained unchanged based on the stationary noise assumption. In the cases of high dynamic or multipath contaminated environments, the noise parameters in the two models need to be -properly tuned. Further investigation can be carried out using the AEKF as the noise-adaptive filter for tuning the noise covariance matrices in the high dynamic or multipath environments. Conflicts of Interest: The author declares that he has no conflicts of interest to report regarding the present study.