1. INTRODUCTION
Since the appearance of the classic paper by Wahba (Reference Wahba1965) a number of improvements to the attitude estimation method have been reported in the literature. This problem was re-considered by Farrell (Reference Farrell1970) who used a Kalman filter along with a Sun sensor and magnetometer sensor data to achieve a ten-fold increase in accuracy of the estimate over the accuracy of the measurement. The problem of the attitude determination may be posed as a purely kinematic problem if one uses a set of gyros to measure the components of the angular velocity vector and the problem reduces to the case of optimum mixing of rate and attitude measurements. Potter and Vander Velde (Reference Potter and Vander Velde1968) adopted such an approach to achieve a hundred-fold increase in accuracy of the estimate over the accuracy of the measurements. An alternate possibility is to determine the angular velocity components from the spacecraft's Eulerian dynamics, but such an approach no longer remains a purely kinematic approach. Then there was the issue of the representation of the attitudes either by the three-parameter varieties of Euler angles or equivalently by the four-parameter representations such as quaternions. In spite of some significant drawbacks due to the need for quaternion normalisation (Lefferts, Markley and Schuster, Reference Lefferts, Markley and Schuster1982), quaternions have emerged as the standard method (Lerner, Reference Lerner, Wertz and Reidel1978) for representing the attitude of a space vehicle.
The problem of quaternion-based attitude estimation is patently nonlinear and most early efforts focussed on the Extended Kalman Filter (EKF). The EKF-based approaches are still exceptionally popular and some recent applications of it are illustrated by Marins et. al. (Reference Marins, Yun, Bachmann, McGhee and Zyda2001) and de Ruiter et. al. (Reference de Ruiter and Damaren2002). However the need to include higher order nonlinear effects in the quaternion prediction and the state update phases as well as in the error covariance prediction and update was soon recognised and a number of extensions have appeared in the literature. The first class of these extensions account for the second-order effects (see for example Vathsal, Reference Vathsal1987) and these employ the mean and covariance updating method outlined in Sage and Melsa (Reference Sage and Melsa1971) and Jazwinski (Reference Jazwinski1970). The EKF employs the Jacobians of the nonlinear transformations of the process states and this has spawned two other extensions to it. The first of these, termed the nonlinear predictive filter, is based on including higher order derivatives of the nonlinear transformations of the process states. Typical applications to the spacecraft attitude problem are outlined by Crassidis and Markley (Reference Crassidis and Markley1997a, Reference Crassidis and Markley1997b). The need to avoid the determination of the Jacobians has led to the development of the Unscented Kalman Filter (UKF) based on the unscented transformation by Julier and Uhlmann (Reference Julier and Uhlmann2000). It must be emphasized though that the UKF is accurate at best to third order and in most practical situations to second order. A recent application of the method to the attitude estimation problem is discussed by Vandyke et. al. (Reference VanDyke, Schwartz and Hall2004).
While most of the early applications were based on using gyroscopes, star trackers, magnetometers and gravity sensors, the recent availability of interferometry measurements of carrier phase differences arising from navigation satellite transmissions has led to cheaper and more reliable alternative sensors for spacecraft attitude determination. Because these sensors operate without the need to rely on the pseudo-random number based codes, they are often referred to as the codeless satellite navigation sensors. They operate by measuring the phase difference in the carrier signals received at two antennas located at opposite ends of a baseline vector. Given the baseline vector d and navigation satellite's sight line unit vector r as well as the carrier signal wavelength, λ, the differential phase measurement is given by:
where N is the number of integer wavelengths occurring between the two antennas. Since the antenna's baseline vector is a function of the spacecraft's attitude, the differential phase measurement Δφ may also be construed to be a nonlinear function of the spacecraft's attitude. A number of studies of the spacecraft attitude determination problem have appeared based on the differential measurement idea. The use of these codeless satellite navigation sensors coupled with the quaternion representations of the attitude, opens up the possibility of determining the large attitudes without any restriction on their magnitudes, over a large time frame. However most of these (see for example Fujikawa and Zimbelman, Reference Fujikawa and Zimbelman1995 and Huang and Juang, Reference Huang and Juang1997) have addressed the problem where the angular velocities are determined by the Eulerian dynamics. Kingston and Beard (Reference Kingston and Beard2004) have also employed the codeless satellite navigation sensors for the attitude determination but their formulations are not exclusively based on the quaternion. Thus the codeless satellite navigation sensors are the only sensors used with no prospect of mixing and the attitudes are not estimated purely from the kinematic considerations. A mixing-type attitude determination filter purely from the kinematics is proposed by Creamer (Reference Creamer1996) but experience with using this approach did not indicate a long enough prediction window for long term autonomous operations. This was identified to be due to the use of the EKF approach which is known to diverge under certain conditions. Crassidis and Markley (Reference Crassidis and Markley2003) have proposed a method of applying the UKF to the problem of spacecraft estimation by extending the multiplicative EKF, using equivalent quaternion multiplications leading to rotation additions rather than component vector additions. They transform the quaternion to modified Rodrigues parameters which need not be explicitly constrained and define the mean and covariance in this space. The transformation and inverse involve the application of trigonometric and inverse trigonometric functions which cannot be adequately approximated by second and third order functions uniformly for all rotations. Although the method is novel, it is not always applicable particularly when disturbance models are naturally additive. A similar method was also proposed by Kraft (Reference Kraft2003). With a view to improving the accuracy and speed of estimation Cheon and Kim (Reference Cheon and Kim2007) have proposed a novel unscented filter in a unit quaternion space where the sigma vectors maintain normalization during the update stage. However they implicitly assume that the entire quaternion normalisation error is due to the magnitude of the rotation which could be restrictive. Moreover in many satellite applications the speed of prediction is not as important as the ability to accurately predict the attitude over a long time frame.
In this paper a new predictor-corrector form of the unscented Kalman filter is proposed to maintain the normalization of the unscented mean quaternion estimate in the presence of additive disturbances. The three-axis attitude of a spacecraft is estimated by integrating the output of three rate gyros corrupted by drifts and zero mean additive white Gaussian noise. An attitude sensor in the form of a satellite navigation carrier phase differential sensor measures the spacecraft attitudes at regular time intervals, and this information is utilized by a nonlinear but optimal estimator, i.e., an UKF, to update both the spacecraft attitude and the gyro drift. This attitude measurement is also corrupted by zero mean noise. Attitude estimates between samples are provided by a kinematic model of the spacecraft's angular motion, relating the aircraft's angular velocity rate vector to the attitude quaternion. The estimator is shown to converge to steady-state operation, while the means and standard deviations or variances of the attitude and gyro drift estimation errors just prior to, and also just after an attitude update are obtained recursively by applying the unscented transformation twice; first to predict the updated estimate and then to correct it to maintain its normalization. The application of the unscented transformation a second time can be done by an almost trivial extension of the UKF algorithm. Thus an attitude determination system that is based on three vector measurements of non-zero, non-colinear vectors is demonstrated. The approach that is taken involves a triad of rate gyros the outputs of which are fused with an ultra-short baseline differential satellite navigation three-axis attitude determination system, without the need to integrate the attitude dynamics. Even though the rate gyros used exhibited poor long term bias stability, the satellite navigation attitude system was used to update the estimate of the gyro biases continuously and obtain attitude estimates over a relatively long term, immaterial of the attitude dynamics at the expense of speed of prediction. Moreover the entire filter is defined in the quaternion space without any need for intermediate transformations to satisfy the quaternion normalization constraint.
2. ATTITUDE KINEMATICS PROCESS MODELLING
The quaternion vector update equation relating the quaternion rate to the angular velocity vector ω and the quaternion normalisation relation are given by well known differential equations. When the angular velocity vector is measured by rate gyros and corrupted by biases and noise, the dynamics may be expressed in terms of the attitude quaternion q as (Vathsal, Reference Vathsal1987):
![{ \dot{\bf q}} \equals \left[ {\matrix{ {{ \dot{\bf q}}} \cr {{ \dot{\bf b}}} \cr} } \right] \equals \left[ {\matrix{ {{1 \over 2}\bmOmega \lpar {\bf u}\rpar } \tab { \minus {1 \over 2}\rmGamma \lpar {\bf q}\rpar } \cr 0 \tab 0 \cr} } \right]\left[ {\matrix{ {\bf q} \cr {\bf b} \cr} } \right] \plus \left[ {\matrix{ { \hskip-1\minus {1 \over 2}\rmGamma \lpar {\bf q}\rpar } \tab 0 \cr 0 \tab {\bf I} \cr} } \right]\left[ {\matrix{ {\bmeta _{\setnum{1}} } \cr {\bmeta _{\setnum{2}} } \cr} } \right] \equiv {\bf A}\lpar {\bf q}\rpar {\bf x} \plus {\bf B}\lpar {\bf q}\rpar \bmeta \comma](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqn2.gif?pub-status=live){kind=small}
where,
![\bmOmega \lpar \bmomega \rpar \equals \left[ {\matrix{ 0 \tab {\omega _{\setnum{3}} } \tab { \minus \omega _{\setnum{2}} } \tab {\omega _{\setnum{1}} } \cr { \minus \omega _{\setnum{3}} } \tab 0 \tab {\omega _{\setnum{1}} } \tab {\omega _{\setnum{2}} } \cr {\omega _{\setnum{2}} } \tab { \minus \omega _{\setnum{1}} } \tab 0 \tab {\omega _{\setnum{3}} } \cr { \minus \omega _{\setnum{1}} } \tab { \minus \omega _{\setnum{2}} } \tab { \minus \omega _{\setnum{3}} } \tab 0 \cr} } \right]\comma \;\bmomega \equals \left[ {\matrix{ {\omega _{\setnum{1}} } \cr {\omega _{\setnum{2}} } \cr {\omega _{\setnum{3}} } \cr} } \right]\comma \;\rmGamma \lpar {\bf q}\rpar \equals \left[ {\matrix{ {q_{\setnum{4}} } \tab { \minus q_{\setnum{3}} } \tab {q_{\setnum{2}} } \cr {q_{\setnum{3}} } \tab {q_{\setnum{4}} } \tab { \minus q_{\setnum{1}} } \cr { \minus q_{\setnum{2}} } \tab {q_{\setnum{1}} } \tab {q_{\setnum{4}} } \cr { \minus q_{\setnum{1}} } \tab { \minus q_{\setnum{2}} } \tab { \minus q_{\setnum{3}} } \cr} } \right]\comma](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqn3.gif?pub-status=live)
and b is a rate-gyro drift rate bias, η1 is a Gaussian white noise process corrupting the angular velocity measurements, while the bias rate is assumed to be driven by another independent Gaussian white noise process, η2.
We also assume that the quaternion rate satisfies the constraint:
To satisfy the constraints (4) and (5) we assume a time step equal to Δt and discretise the Equation (2) as:
The evaluation of F is done in accordance with the approach recommended by Markley (Reference Markley, Wertz and Reidel1978) which is equivalent to Equation (6) and given by the pair of equations:
where,
![\eqalign{ \tab {\bf F}_{\bf \setnum{1}} \lpar \bmomega \comma {\bf b}\rpar \equals \left[ {\matrix{ {\alpha _{k} {\bf I}} \tab {\beta _{k} \lpar \bmomega \plus {\bf b}\rpar } \cr { \minus \beta _{k} \lpar \bmomega \plus {\bf b}\rpar ^{\bf T} } \tab {\alpha _{k} } \cr} } \right]\comma \;\;{\bf G}_{\bf \setnum{1}} \lpar {\bf q}\rpar \equals \minus {{\rmDelta t} \over 2}\rmGamma \lpar {\bf q}\rpar \left[ {\matrix{ {\bf I} \tab {\bf 0} \cr} } \right]\comma \;{\bf G}_{\bf \setnum{2}} \equals \rmDelta t\left[ {\matrix{ {\bf 0} \tab {\bf I} \cr} } \right]\comma \cr \tab \alpha _{k} \equals \cos \left( {{{\rmDelta t} \over 2}\left\Vert {\bmomega \plus {\bf b}} \right\Vert} \right)\comma \;\beta _{k} \equals {{\sin \left( {{{\rmDelta t} \over 2}\left\Vert {\bmomega \plus {\bf b}} \right\Vert} \right)} \mathord{\left/ {\vphantom {{\sin \left( {{{\rmDelta t} \over 2}\left\Vert {\bmomega \plus {\bf b}} \right\Vert} \right)} {\left\Vert {\bmomega \plus {\bf b}} \right\Vert}}} \right. \kern-\nulldelimiterspace} {\left\Vert {\bmomega \plus {\bf b}} \right\Vert}}. \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqnU1.gif?pub-status=live)
The last column of the matrix F1 may be interpreted as change in the quaternion components due to the angular velocity vector giving the updates to the quaternion given an initial rotation magnitude or equivalently the last component of the quaternion. Thus if the quaternion normalization error is assumed to be entirely due to the last component, as is done by Cheon and Kim (Reference Cheon and Kim2007), the corrections to the quaternion components due to an angular rate error correction or a multiplicative error quaternion qe+, can be computed from the last column of the matrix F1. The quaternion update can then be done by quaternion multiplication defined by:
The use of Equation (9) involves a nonlinear transformation of the random variable representing the noise vector. The noise in the rate gyro measurements could have also been transformed to the modified Rodrigues parameters, to maintain the post-update quaternion normalization. This would lead to a filter similar to that of Crassidis and Markley (Reference Crassidis and Markley2003) provided the quaternion measurement error is also modelled as a multiplicative quaternion. A Gaussian distributed vector random variable, when transformed to an equivalent quaternion or a rotation need not remain Gaussian distributed. Hence we choose to include the rate gyro noise vector as an additive Gaussian distributed vector to demonstrate an alternate approach to unscented filtering. An interesting feature of this formulation is that mean quaternion update continues to be normalized if it is initially normalized. The matrices G1(q) and G2 are only needed for purposes of updating the covariance and for evaluating the sigma points. The mean quaternion at the sigma points therefore continues to be normalized if the initial sigma point vector is normalized. However the state estimate does not preserve the normalization as these are evaluated as weighted linear combinations of the means at the sigma points. For this reason we propose a predictor-corrector implementation of the filter (this is in addition to the already existing propagate and update structure of the Kalman filter), which is based on a process and measurement predictor model given by:
The post-estimate process and output corrector model is given by:
Equations (10 and 11) are used to implement an UKF while Equations (12 and 13) are used to define a second unscented transformation to renormalize the estimated mean and modify the corresponding covariance of the estimates. The implementation of the unscented transformation and the UKF are discussed in section 4.
3. CODELESS SATELLITE NAVIGATION CARRIER PHASE MEASUREMENT MODEL
At the outset we assume we are able to construct codeless satellite navigation receiver antenna pairs located at the two ends of a graphite-epoxy bar which is less than the wavelength λ of the carrier signal. When this is not the case it is first required that the ambiguity in N is resolved. The process of estimating the correct carrier phase integer ambiguity is called satellite navigation ambiguity resolution. In order to reach centimetre level positioning accuracies, this ambiguity term must be determined correctly. When multiple satellites are available, the residual sensitivity matrix approach proposed by Hatch and Sharpe (Reference Hatch and Sharpe2001) provides a method for estimating N. The most successful methods of ambiguity resolution are based on using pseudo-satellite (pseudolite) carrier wave transmitters (Pervan, Cohen and Parkinson, Reference Pervan, Cohen and Parkinson1994) to complement measurements from available satellites. When it can be assumed that we are able to construct codeless satellite navigation receiver antenna pairs located at the two ends of a graphite-epoxy bar which is less than the wavelength λ of the carrier signal, it is possible in principle to assume that N equals a known integer or zero in Equation (1) which may be expressed as,
The GPS signals are transmitted on two radio frequencies in the UHF band. These frequencies are referred to as L1 and L2, with:
where f 0=10·23 MHz, is a common frequency. The GPS signal is a phase-modulated signal using bi-phase shift keying (BPSK). The phase change rate is often referred as the chip rate. For the case of the highest frequency carrier which is the L1 carrier signal (1575·42 MHz), the distance between the two antennas located at the two ends of the graphite-epoxy bar should be less than 0·19 m. The distance between the two antennas would be sufficient to operate even at the L2 carrier (1227·6 MHz) frequency. In practice these frequencies ought to be Doppler shifted. It is in fact feasible to construct such an antenna by locating two micro-strip patch antennas at the two ends of the bar although the field of view is considerably limited. Yet the two-channel receiver can in principle be constructed and the differential phase angle extracted by employing the fast Fourier Transform (FFT) algorithm implemented on suitable digital signal processing hardware.
The GLONASS satellite navigation system offers many features in common with the GPS satellite navigation system. It is also based on 24 orbiting satellites forming the space segment (21 operational satellites with 3 in-orbit spares) but will use only 3 orbital planes separated by 120° of longitude and with equal spacing between satellites of 45° within the plane. The GLONASS orbits are roughly circular orbits with an inclination of about 64·8°, a semi-major axis of 25440 Km and a period of 11h 15m 44s. Similar to GPS, GLONASS also transmits two spread-spectrum signals in the L-band at around the same power levels (−160 dBW at L1, −163 dBW at L2). However individual GLONASS satellites are distinguished by a dedicated radio frequency channel rather than spread-spectrum code. In GLONASS a single code of length 511 bits repeating every 1ms is used. Information is encoded differentially in an RZ (return to zero) format with a final data rate of 50 baud.
Radio-frequency carriers used by GLONASS occupy channels within the L-band ranging from 1240–1260 MHz and 1597–1617 MHz, the channel spacing being 9/16 (or 0·5625) MHz at the higher frequencies and 7/16 (or 0·4375) MHz at the lower frequencies. The carrier frequencies themselves are also multiples of channel spacing and the number of planned channels is 24. GLONASS L1 transmission carrier frequencies (f n, in megahertz) and channel numbers (C n) are related by the expression:
A similar formula relates the GLONASS L2 transmission carrier frequencies to the channel numbers where the frequencies are in the ratio 7/9.
This is quite different from the GPS system which uses the same frequency for all satellites and differentiates one satellite from another by individual Gold codes, a form of code division multiplex (CDM); the GLONASS system uses frequency division multiple access (FDMA) to distinguish between satellites. This difference between the two systems has a major impact in designing codeless receivers capable of joint operation. However receivers capable of measuring carrier phase differential have been designed and built for both the GLONASS and GPS systems. In fact a receiver with an antennas separation distance of 0·18m would be adequate for the relation (14) to hold.
The change in attitude of the spacecraft over a period of time could be observed by comparing the current measured phase differential with the initial phase differential measured at some initial reference time. Thus this difference in the measured phase differential could be expressed as:
where rB is the navigation satellite's sight line vector at the current time and r0 is the navigation satellite's sight line vector at the initial reference time. The navigation satellite's sight line vector rB could be expressed in terms of the satellite's body coordinates. However since the body attitude may be defined in terms of the quaternion, the transformation relating the estimate of current sight line vector 
in the orbiting coordinates to the current sight line vector rB in body coordinates may be expressed in terms of the quaternion components. Hence:
where,
![{\bf T}\lpar {\bf q}\rpar \equals \left[ {\matrix{ {q_{\setnum{4}}^{\setnum{2}} \plus q_{\setnum{1}}^{\setnum{2}} \minus q_{\setnum{2}}^{\setnum{2}} \minus q_{\setnum{3}}^{\setnum{2}} } \tab {2\lpar q_{\setnum{1}} q_{\setnum{2}} \plus q_{\setnum{3}} q_{\setnum{4}} \rpar } \tab {2\lpar q_{\setnum{1}} q_{\setnum{3}} \minus q_{\setnum{2}} q_{\setnum{4}} \rpar } \cr {2\lpar q_{\setnum{1}} q_{\setnum{2}} \minus q_{\setnum{3}} q_{\setnum{4}} \rpar } \tab {q_{\setnum{4}}^{\setnum{2}} \minus q_{\setnum{1}}^{\setnum{2}} \plus q_{\setnum{2}}^{\setnum{2}} \minus q_{\setnum{3}}^{\setnum{2}} } \tab {2\lpar q_{\setnum{2}} q_{\setnum{3}} \plus q_{\setnum{1}} q_{\setnum{4}} \rpar } \cr {2\lpar q_{\setnum{1}} q_{\setnum{3}} \plus q_{\setnum{2}} q_{\setnum{4}} \rpar } \tab {2\lpar q_{\setnum{2}} q_{\setnum{3}} \minus q_{\setnum{1}} q_{\setnum{4}} \rpar } \tab {q_{\setnum{4}}^{\setnum{2}} \minus q_{\setnum{1}}^{\setnum{2}} \minus q_{\setnum{2}}^{\setnum{2}} \plus q_{\setnum{3}}^{\setnum{2}} } \cr} } \right].](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqn19.gif?pub-status=live)
An estimate of current sight line vector in the orbiting coordinates can generally obtained by an independent Kalman filter or by employing an algorithm such as NORAD's SDP4, SDP8 or SGP4 methods (Hoots et. al. Reference Hoots, Schumacher and Glover2004). It therefore follows that the difference in the measured phase differential could be expressed as:
and using Equation (19) we may write, T(q)−I as:
![{\bf T}\lpar {\bf q}\rpar \minus {\bf I} \equiv \rmDelta {\bf T}\lpar {\bf q}\rpar \equals \minus 2\left[ {\matrix{ {\lpar q_{\setnum{2}}^{\setnum{2}} \plus q_{\setnum{3}}^{\setnum{2}} \rpar } \tab { \minus \lpar q_{\setnum{1}} q_{\setnum{2}} \plus q_{\setnum{3}} q_{\setnum{4}} \rpar } \tab { \minus \lpar q_{\setnum{1}} q_{\setnum{3}} \minus q_{\setnum{2}} q_{\setnum{4}} \rpar } \cr { \minus \lpar q_{\setnum{1}} q_{\setnum{2}} \minus q_{\setnum{3}} q_{\setnum{4}} \rpar } \tab {\lpar q_{\setnum{1}}^{\setnum{2}} \plus q_{\setnum{3}}^{\setnum{2}} \rpar } \tab { \minus \lpar q_{\setnum{2}} q_{\setnum{3}} \plus q_{\setnum{1}} q_{\setnum{4}} \rpar } \cr { \minus \lpar q_{\setnum{1}} q_{\setnum{3}} \plus q_{\setnum{2}} q_{\setnum{4}} \rpar } \tab { \minus \lpar q_{\setnum{2}} q_{\setnum{3}} \minus q_{\setnum{1}} q_{\setnum{4}} \rpar } \tab {\lpar q_{\setnum{1}}^{\setnum{2}} \plus q_{\setnum{2}}^{\setnum{2}} \rpar } \cr} } \right]\comma](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqn21.gif?pub-status=live)
which is a homogeneous quadratic function of the components of the quaternion. Thus a discrete measurement of the error in the difference of the phase differentials due to changes in the attitude can be expressed as:
where v is an additive Gaussian random variable representing a white noise or delta-correlated stochastic process. For three-axis measurement of the attitude one would require three independent measurements which may be expressed as:
We observe that the process model given by Equations (6) and the measurement model given by Equation (23) are both nonlinear in the components of the quaternion. Generally the measurements cannot be cast as multiplicative quaternion error models. The estimation problem may now be stated: Given a sequence of noisy observations by Equation (21) we need to estimate the sequence of state vectors of the non-linear system driven by noise, defined by Equations (6). This problem is now dealt with by the using an UKF after briefly discussing the reasons for its choice in the next section.
4. NONLINEAR AND UNSCENTED KALMAN FILTERS
Most dynamic models employed for purposes of estimation or filtering of pseudo range errors or orbit ephemeris errors are generally not linear. To extend and overcome the limitations of linear models, a number of approaches such as the EKF have been proposed in the literature for nonlinear estimation using a variety of approaches. Unlike the Kalman filter, the EKF may diverge, if the consecutive linearizations are not a good approximation of the linear model over the entire uncertainty domain. Yet the EKF provides a simple and practical approach to dealing with essential non-linear dynamics. The main difficulty in applying the algorithm to problems related to the estimation of a spacecraft's orbit and attitude is in determining the proper Jacobian matrices. The UKF is a feasible alternative that has been proposed to overcome this difficulty, by Julier, Uhlmann and Durrant-Whyte (Reference Julier, Uhlmann and Durrant-Whyte2000), as an effective way of applying the Kalman filter to nonlinear systems. It is based on the intuitive concept that it is easier to approximate a probability distribution than it is to approximate an arbitrary nonlinear function or transformation, of a random variable.
The UKF gets its name from the unscented transformation, which is a method of calculating the mean and covariance of a random variable undergoing nonlinear transformation y=f(w). Although it is a derivative free approach, it does not really address the divergence problem. In essence the method constructs a set of sigma vectors and propagates them through the same non-linear function. The mean and covariance of the transformed vector are approximated as a weighted sum of the transformed sigma vectors and their covariance matrices.
Consider a random variable w with dimension L which is going through the nonlinear transformation, y=f(w). The initial conditions are that w has a mean 
and a covariance Pww. To calculate the statistics of y, a matrix χ of 2L+1 sigma vectors is formed. We have chosen to use the scaled unscented transformation (Julier Reference Julier2002) as this transformation gives one the added flexibility of scaling the sigma points to ensure that the covariance matrices are always positive definite.
Given a general discrete nonlinear dynamic system in the form,
where xk∊R n is the state vector, uk∊R r is the known input vector, yk∊R m is the output vector at time k. wk and vk are, respectively, the disturbance or process noise and sensor noise vectors, which are assumed to Gaussian white noise with zero mean. Furthermore Qk and Rk are assumed to be the covariance matrices of the process noise sequence, wk and the measurement noise sequence, vk respectively. The unscented transformations of the states are denoted as,
while the transformed covariance matrices and cross-covariance are respectively denoted as:
and
The UKF estimator can then be expressed in a compact form. The state time-update equation, the propagated covariance, the Kalman gain, the state estimate and the updated covariance are respectively given by:
Equations 28–32 are in the same form as the traditional Kalman filter and the extended Kalman filter. Thus higher order non-linear models capturing significant aspects of the dynamics may be employed to ensure that the Kalman filter algorithm can be implemented to effectively estimate the states in practice. For our purposes we adopt both the UKF approach to estimate the attitude quaternion and other states in the process model.
The UKF is based on approximating the probability distribution function rather than approximating a nonlinear function as in the case of EKF. The state distributions are approximated by a Gaussian probability density, which is represented by a set of deterministically chosen sample points. The nonlinear filtering using the Gaussian representation of the posterior probability density via a set of deterministically chosen sample points is the basis for the UKF. It is based on statistical linearization of the state dynamics rather than analytical linearization (as in the EKF). The statistical linearization is performed by employing linear regression using a set of regression (sample) points. The sigma points are chosen as the regression points. The mean and covariance at the sigma points then represent the true mean and covariance of the random variable with the particular Gaussian probability density. Thus when transformed to the nonlinear systems, they represent the true mean and covariance accurately only to the second order of the nonlinearity. Thus this can be a severe limitation of the UKF unless the nonlinearities can be limited to the first and second order in the process model.
In the predictor-corrector formulation of the UKF (PCUKF) the estimate and the covariance are corrected a second time. The post-estimate corrections are done by employing a pair of unscented transformations and noise-free model of the form:
In our application Equations 33 correspond to Equations 12 and 13. The transformed output mean and covariance now represent the final corrected estimate and its covariance. Although this process involves a second application of the unscented transformation, it facilitates the use of the standard UKF without the need for constraining the sigma vectors in any way. The complete filter may be expressed in terms of two nonlinear unscented transformations applied sequentially as:
The state estimate 
is given by the result of applying the first of the two unscented transformations in the sequence. Thus the filter is equivalent to two stage propagation accompanied by a single measurement update stage.
5. PCUKF BASED ATTITUDE ESTIMATION
The success of the application of the UKF depends largely on the approximation to the covariance which is estimated as a weighted linear sum of the covariance at the sigma points. When this approximation is such that the covariance is not positive definite, the UKF algorithm fails as the Cholesky decomposition is not possible. To ensure that this covariance is positive definitive it is essential to adjust the scaling parameter αus, if and when necessary. In the example below, αus was chosen to be a very small positive number in the case of both the transformations. Alternately one could start the first few steps by using an EKF and then switch to the UKF.
In the first instance it was attempted to estimate the attitude of the spacecraft by employing a single carrier phase sensor. As expected a single sensor was not able to produce a meaningful estimate of all the components of the attitude quaternion over a relatively large time frame (1440 minutes). This was due to the need for three vector measurements of non-zero, non-colinear attitude vectors for purposes of the attitude estimation from kinematic considerations.
The three receiver based carrier phase sensor was then simulated. The spacecraft, a navigation satellite in a typical GLONASS orbit, was assumed to performing sustained three axis rotations with a body axes angular velocity vector given by ![\bmomega \equals \left[ {\matrix{ {0{\cdot}1} \tab {0{\cdot}2} \tab {0{\cdot}3} \cr} } \right]{\rm \ rads\sol hr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022094933549-0623:S0373463309990075_eqnU5.gif?pub-status=live){kind=small}
. These are measured by three rate gyros and the measurements are assumed to be biased where the bias rate is driven by white noise. To establish the parameters for implementing the PCUKF the standard deviations were assumed to be the same order of magnitude as the rate gyros (0·002 arc-s/s) considered by Farrenkopf, Reference Farrenkopf1978 where the single axis attitude estimation problem was considered with attitude measurements provided by a star tracker (σ=20 arc-s) and an update rate of 0·1 per minute. The standard deviation of the carrier phase measurements were assumed to be between σ=50–70 arc-s. Independent simulations of the process and the measurements, including the rate-gyros were performed and standard deviations of the white noise sources were assumed to be 30 times larger than those assumed for the design calculations. Figure 1 illustrates the simulated and estimated attitude quaternion plotted to the same scale. Figure 2 shows the corresponding errors in the simulated quaternion and estimated quaternion over the same time frame.
Comparison of the simulated and the estimated attitude quaternion over a time frame of 1440 minutes.
Errors in the components of the simulated and the estimated attitude quaternion over a time frame of 1440 minutes.
Figure 3 illustrates the errors in the corresponding 3, 2, 1 sequence Euler angles, ψ, θ and φ. Attitude accuracies of the order of ±1°–±2° may be obtained by using antenna patches separated by distances of the order 0·2m. A tenfold decrease in the order of the accuracies (errors) is possible if the antenna patch distance is increased fivefold. In Figure 4 the growth of the simulated and estimated bias states, related to the gyro drift rate with the PCUKF are shown. In Figure 5 are shown the growth of the simulated and estimated bias states, related to the gyro drift rate without the second unscented transformation. It is clear that in the long term the bias states slowly tend to diverge from the corresponding simulated states.
Errors in the components of the simulated and the estimated attitude Euler angles φ, θ and ψ over a time frame of 1440 minutes
Growth and comparison of the simulated and estimated bias states in the gyros with the PCUKF.
Growth and comparison of the simulated and estimated bias states in the gyros without the second unscented transformation.
Comparing the PCUKF with the standard UKF, although there were practically no differences between the predicted mean and covariance of the estimate, it was observed that the bias states slowly began to diverge from the simulations in the case of the standard UKF. Thus the performance of the PCUKF was relatively more stable over longer time frame.
Figure 6 shows the simulated measurement error of a typical carrier phase sensor. The carrier frequency was assumed to be 1602 MHz. It was assumed that three pairs of patch antennas are located along mutually perpendicular three-axes and that each pair was separated by 18cms. Satellite navigation receiver based carrier phase measurements are subject to additional measurement errors, such as tropospheric and ionospheric errors, orbital ephemeris errors, multipath, frequent cycle slip occurrences, antenna phase centre variation and noise. In addition to measuring the arrival time of the satellite navigation signal using the code modulation, receivers that also measure the phase of the carrier frequency and the total phase change (both whole and partial cycles) since the initial measurement are currently available. This measure, referred to as the integrated carrier phase is statistically independent of the code measurement and is also about two orders of magnitude less noisy. Yet these error sources severely deteriorate attitude determination availability. Moreover these errors to the carrier phase cannot be modelled as multiplicative quaternion errors because of the nature of the probability density functions. To model these errors adequately the covariance of the errors in the measurements were assumed to be much higher than was assumed at the design stage. This fact is adequately reflected in Figure 6.
Simulated measured carrier phase error (in radians) of a typical carrier phase sensor.
Finally it must said that the filter was run over a much longer time frame (2882 minutes) and the performance of the filter did not deteriorate in spite of this long term operation. Thus the implementation of an attitude determination system over a relatively long time frame is successfully demonstrated.
6. CONCLUSIONS
In this paper a new predictor-corrector form of the unscented Kalman filter is presented and validated to maintain the normalization of the unscented mean quaternion estimate in the presence of additive disturbances. The method is a simple and straightforward extension of the UKF requiring no transformations of the states. We have also demonstrated the feasibility of a long term attitude determination system that is based on three vector measurements of non-zero, non-colinear vectors. Long term estimation of the attitude was facilitated by the use of the PCUKF. Comparisons with the results of several other studies, referenced in paper, clearly pointed to the superior performance of the PCUKF based attitude determination system, particularly in terms of the long term stability and operation of the filter.
