Estimator-Based GPS Attitude and Angular Velocity Determination

: In this paper, the estimator-based Global Positioning System (GPS) attitude and angular velocity determination is presented. Outputs of the attitude estimator include the attitude angles and attitude rates or body angular velocities, depending on the design of estimator. Traditionally as a position, velocity and time sensor, the GPS also offers a free attitude-determination interferometer. GPS research and applications to the field of attitude determination using carrier phase or Doppler measurement has been extensively conducted. The raw attitude solution using the interferometry technique based on the least-squares approach is inherently noisy. The estimator such as the Kalman filter (KF) or extended Kalman filter (EKF) can be incorporated into the GPS interferometer, potentially providing several advantages, such as accuracy improvement, reliability enhancement, and real-time character-istics. Three estimator-based approaches are investigated for performance comparison, including (1) KF with measurement involving attitude angles only; (2) EKF with measurements based on attitude angles only; (3) EKF with measurements involving both attitude angles and body angular rates. The assistance from body mounted gyroscopes, if available, can be utilized as the measurements for further performance improvement, especially useful for the case of signal-challenged environment, such as the GPS outages. Modeling of the dynamic process involving the body angular rates and derivation of the related algorithm will be presented. Simulation results for various estimator-based approaches are conducted; performance comparison is presented for the case of GPS outages.

The remainder of this paper is organized as follows. Preliminary background on GPS attitude and angular velocity determination is reviewed in Section 2. In Section 3, attitude and angular velocity estimator using the extended Kalman filter is briefly introduced. In Section 4, illustrative examples are presented to evaluation of the effectiveness of the proposed algorithm. Conclusions are given in Section 5.

GPS Attitude and Angular Velocity Determination
Carrier phase observables in GPS include sum of range, an unknown integer ambiguity and some ranging errors, expressed as: where the parameters are defined as 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 correction; d trop : tropospheric correction; λ: carrier phase wavelength; N: carrier phase integer ambiguity; v ρ , v φ : measurement noises of code and carrier phases.
GPS interferometry has been firstly used in precise static relative positioning, and thereafter in kinematic positioning. By using carrier phase observables, the relative positioning is obtained in centimeter level provided that the integer ambiguity is resolved. In the beginning of 1990s, Van Grass et al. [15] conducted research on the GPS to the field of aircraft attitude determination using carrier phase was developed. In their work, the receiver-satellite double difference observation was employed.
The receiver-satellite double differenced observable is formed by a linear combination of four observations: where '•' could represent either Φ, r, or N, and the following operators are defined : between-receiver differencing operator; ∇(•) = (•) satellite_j − (•) satellite_i : between-satellites differencing operator. Double differenced observable eliminates or greatly reduces the satellite and receiver timing errors. The observation equation between two receivers and two satellites by combining phase data from master (denoted as '1') and slave (denoted as '2') receivers to satellites i and j can be written as For a short baseline, e.g., less than one meter between two antennas, ionospheric and tropospheric parameters become negligible. The resulting double-difference phase equation when ignoring atmospheric, satellite ephemeris, and residual clock errors is given by Attitude is defined by the rotation transformation that relates a coordinate system fixed in space to a coordinate system fixed in the body (body frame). The space coordinate system is typically defined to be a local level NED (north-east-down) frame. The baseline vector is defined as the vector between the antennas designated master and one of the slave antennas. The carrier phase differential GPS based on interferometer principles can solve for the antenna vector. The approach is based on the difference in the GPS carrier phase received at two or more nearby antennas connected to a rigid body. Multiple antenna vectors from an antenna array can be used to calculate the vehicle attitude. In general, to determine the three-dimensional attitude, three non-collinear antennas simultaneously receiving signals from two satellites are the minimum requirement. Referring to the configuration as in Fig. 1, 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: Similarly, the single differenced observable received for the same antennas from satellite j is where b is the baseline vector formed by two antennas, and e represents the line-of-sight unit vector from antennas to satellites. The receiver-satellite double difference is obtained by taking two independent single differences: When the integer ambiguity is known, the range-equivalent of Eq. (9) is ⎡ which can be expressed in matrix form ∇ r = Eb.
There have been many methods proposed for the integer ambiguity resolution, typically including ambiguity function, antenna exchange/swap, baseline rotation methods. The residuals are the differences between the actual double difference measurements and those predicted based on the leastsquares solution for b = [ b x b y b z ] T , the approximate interferometer coordinates, and the known satellite coordinates. 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.
There are several methods available for solving vehicle attitudes, typically including Euler angle method and quaternion method [12,21]. The rotation transformation matrix that relates the NED and body frames provides the information for finding the vehicle attitude: where the subscripts n and b represent the local and body frames, respectively, and we defined the notations S ≡ sin and C ≡ cos. Since R n2b is an orthonormal matrix, its inverse is easily founded to be Finally, the vehicle attitude can be obtained through the relation: To perform the calculation of attitude determination, the antennas should be set up as nonzero and are on the non-collinear vector. Once the baseline vector is determined, estimation of the coordinate frame transformationR b2n can be achieved and subsequently the transformation matrix from body to local frame can be estimated through the calculation . Once the Euler attitude angle rates are available, the body angular velocity can be found ⎡ However, the attitude angle rate directly obtained from the difference of attitude angles, = [ φ θ ψ ] T between epochs k and k − 1 with time interval t: is usually corrupted and too noisy. Utilization of the Kalman type of filters provides resolves some of the problem.

Attitude and Angular Velocity Estimator Using the Extended Kalman Filter
The EKF has been successfully applied in the GPS receiver position/velocity determination, and the integrated GPS/INS navigation system design. In this paper, the KF and EKF will be employed as the estimator for attitude, attitude rates, and body angular velocities to enhance the accuracy and reliability of the attitude solution.

Implementation of the KF and EKF as the State Estimator
Consider a dynamical system whose state is described by a linear, vector differential equation. The process model and measurement model are represented as The Kalman filter equations are summarized in Tab. 1.
k Consider a dynamical system whose state is described by a nonlinear, vector difference equation. Assuming the process to be estimated and the associated measurement relationship may be written in the form: where the state vector x k ∈ n , process noise vector w k ∈ n , measurement vector z k ∈ m , and measurement noise vector v k ∈ 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 δ ik stands for the Kronecker delta function: Without external aiding such as an inertial sensor to provide a reference trajectory, the process dynamics represent the total observer dynamics, as shown in Fig. 2. Implementation algorithm for the discrete EKF equations is provided in Tab. 2.

Modeling of Attitude Dynamics
The appropriate description of the dynamic process model depends on the type of dynamics encountered in a given application.
where the linear approximation equations for system and measurement matrices are obtained through the relations

.1 Linear System Model Using the KF as the Estimator
A three-state filter model involving the state vector x = = φ θ ψ T is suitable for a stationary receiver where the attitude states are described as the random walk. This model is equivalent to the P (Position) model in GPS positioning. If the observer is moving nearly constant angular rate, the model corresponding to the double integrator will be more suitable. When the velocity cannot reasonably be modeled as constant, then acceleration states can be added.
In addition to a pure random walk, an alternative model using the body angular acceleration for the process model can be described as This model is reasonable since accelerations may not be always constant, but are correlated over short time intervals. Based on Eq. (18), the dynamic process model is given by Furthermore, they can also be modeled as the Gauss-Markov process: T , the CAA filter model is given by

Nonlinear System Model Using the EKF as the Estimator
If the body angular velocity information is required, an alternative dynamic model can be applied such that the body angular velocity (p, q, r) estimate can be provided directly via the filter. This design provides a convenient way when the system is assisted by the body mounted rate gyroscopes, which is especially useful in case the GPS attitude information is not available.
(1) CAV model for the EKF The state vector is composed of 6 states: x = φ θ ψ p q r T . The attitude rates and the body angular velocity are related bẏ φ = p + qsin(φ)tan(θ) + r cos(φ) tan(θ) θ = q cos(φ) − r sin(φ) ψ = qsin(φ)sec(θ) + r cos(φ) sec(θ) The above system of equations can be written in matrix form ⎡ ⎣φ θψ By combining the dynamic models of the body angular velocities:ṗ = u p ;q = u q ;ṙ = u r , we have the following dynamic process model One of the alternative models for the dynamic of body angular rate is given by:ṗ = − 1 τp p + u p ; q = − 1 τq q + u q ;ṙ = − 1 τr r + u r . In such case, we have Similarly, the Gauss-Markov can also be employed for the case of:ṗ = −β p p + 2σ 2 p β p u p ; q = −β q q + 2σ 2 q β q u q ;ṙ = −β r r + 2σ 2 r β r u r . The Jacobian matrix F = ∂f ∂x can be derived to be and the corresponding discrete-time model is given by φ k = exp(F · t).
(2) CAA model for the EKF The state vector for the CAA model involves 9 states: x = φ θ ψ p q rṗqṙ T , where the models for angular accelerations can be represented as:p = u p ;q = u q ;r = u r , yielding the dynamic model Similarly, the alternative model that can be considered is given byp = − 1 τpṗ + u p ;q = − 1 τqq + u q ; r = − 1 τrṙ + u r , leading to the process model: As mentioned, the Gauss-Markov process can also be employed when necessary: p = −β pṗ + 2σ 2 p β p u p ;q = −β qq + 2σ 2 q β q u q ;r = −β rṙ + 2σ 2 r β r u r

Illustrative Examples
Simulation was conducted for evaluating the estimation performance of GPS-based attitude and angular velocity determination. Computer codes were developed using the Matlab ® software with assistance from the commercial software Satellite Navigation (SatNav) Toolbox by GPSoft LLC [24] employed for generation of the GPS satellite orbits/positions and thereafter, the satellite pseudorange, and carrier phase observables, required for simulation. The attitude dynamics for simulation is as follows: p = π 15 sin π 60 t ; q = π 20 cos π 150 t ; r = 3π 100 cos π 300 t + 0.03 with initial attitude φ 0 = θ 0 = ψ 0 = 0. The attitude rates and body angular velocities in the simulation are shown in Fig. 3.
Three algorithms are involved for testing their performance, including (1) KF with measurement involving attitude angles only; (2) the first algorithm using EKF (herein referred to as the EKF1), for which the measurements involving attitude angles only; (3) the second algorithm using EKF (herein referred to as the EKF2), for which the measurements involving both attitude angles and body angular rates (possibly from gyroscopes).
Tab. 3 summarizes the three architectures employed in this paper for comparison evaluation.
To avoid the filter divergence, it should be noticed that setting the process noise PSD to zero is not acceptable, since it can be very detrimental to the filtering performance, especially if the vehicle is maneuvering rapidly. As a result, when the system reaches steady-state condition, the steady-state gain won't become zero and, subsequently, the filter is able to follow the time-varying attitude dynamics. A small value is admissible in practical design. However, one still needs to find the appropriate values that meet the requirement. The CAV model is employed for investigation since the dynamic for the illustrative is moderate. For the CAV model, the noise strengths for the process and measurement models are selected as:   Fig. 4 presents estimate of attitude angles without signal outage and Fig. 5 presents those with signal outage. It can be seen that the results from both architectures are equivalent for the case that no outage occurs while the results from EKF1 provides improved performance in some situation, but not so noticeably. Figs. 6 and 7 show the performance comparison for EKF1 and EKF2 without signal outage. Fig. 6 shows estimate of attitude angles and Fig. 7 provides estimate of body angular rates. It can be seen that the EKF2 approach, with more measurement information available, outperforms the EKF1 approach noticeably. Notice that the solution provided by EKF2 is very close to the true states for the high quality rate gyros. Comparison of the error statistics for EKF1 and EKF2 regarding the attitude angles and body angular rates, respectively, are summarized in Tabs. 4 and 5, respectively. In addition, Fig. 8 shows estimate of attitude angles with signal outage for KF and EKF2. The assistance from body mounted gyroscopes, if available, can be utilized as the measurements for further performance improvement, especially useful for the case of GPS outages. Incorporation of the information of angular velocities into the EKF as the measurement not only achieves better accuracy but also catches the attitude dynamics well.  The above examples has demonstrated the attitude solutions obtained by KF and EKF based on measurements involving attitude angles only vs. those by EKF based on measurements involving both attitude angles and body angular velocities, possibly available from gyroscopes. In the case of signal outage, the body angular velocity is especially helpful from the view point of continuity of service and accuracy improvement.  The raw attitude solutions provided by the interferometry technique based on the least-squares approach are inherently noisy. The performance can be improved noticeably by using the Kalman family of filters. Incorporating the filter into the GPS interferometer effectively improves the performance in estimation performance. However, the case for GPS signal outage leads to performance degradation and even divergence. The proposed design provides the opportunity when body angular rates information from other sensor such as the gyroscopes is available for better resistance to the signal-challenged environment.
Details on dynamic models, including the CAV model as the second-order filter, and CAA model as the third-order filter, respectively, have been discussed. Numerical simulation has been carried out using the Kalman filter and extended Kalman filter, using the CAV model. Results involved include the estimate of attitude angles, attitude rates, and body angular velocities, for the cases with and without signal outages. Examples shown includes the attitude solutions obtained by KF and EKF based on measurements involving attitude angles only vs. those by EKF based on measurements involving both attitude angles and body angular velocities, possibly from gyroscopes. In the case of signal outage, the body angular velocity, possibly from gyroscopes is noticeably helpful from the view point of continuity of service and accuracy improvement.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.