1. INTRODUCTION
The integration between the Global Positioning System (GPS) and inertial navigation systems (INS), long holding the interest of navigation researchers, has recently drawn more attention due to increasing concerns over GPS signal vulnerability. INS systems operate on dead reckoning principles and are naturally complementary to GPS systems, which are based on positioning. Moreover, inertial sensors are impervious to either jamming or interference, both of which can affect GPS receivers. The two systems fundamentally differ in that when the GPS signal is available (interference, jamming or signal blockage is not present), position output is stable and reliable. In contrast, INS position outputs tend to drift over time due to the integration of imperfect acceleration and angular rate measurements. Nevertheless, the two systems can be combined in such a way that INS errors are calibrated by GPS when satellite signals are available. Any subsequent, brief GPS satellite outages are then bridged by the relatively accurate INS position output.
Traditionally, in tightly coupled GPS/INS integration, GPS is utilized as an aiding sensor for INS calibration. The difference between GPS and INS information realized in the range domain is formed as the input to a complementary Kalman filter (KF) for navigation error estimation. In general, there are two mechanizations: one uses the difference between computed INS range measurements and measured GPS range measurements to calibrate INS errors (Ko, Reference Ko2000) (Scherzinger, Reference Scherzinger2000) (Marty and Pagnucco, Reference Marty and Pagnucco1992) (Moafipoor and Brzezinska, Reference Moafipoor, Brzezinska and Toth2004), the other uses differences in carrier-phase measurements in time to calibrate the INS (mainly the velocity states) (Moafipoor and Brzezinska, Reference Moafipoor, Brzezinska and Toth2004) (Farrel, Reference Farrell2002). However, in these approaches, the temporal and spatial decorrelation in GPS measurements is not strictly modelled.
This work adopts a tightly couple GPS/INS integration system which includes rigorous temporal and spatial decorrelation models for GPS measurement errors: INS is embedded in a GPS navigation filter to streamline the filter time propagation. By fusing all GPS and INS estimates through state dynamics in one centralized extended Kalman filter (EKF) (Gao and Krakiwsky, Reference Gao, Krakiwsky, Abousalem and McLellan1993), GPS positioning and INS calibration are accomplished simultaneously in a single measurement update, which has exactly the same observation matrix for the measurement update as a GPS-only system.
Prior research in (Johnson and Lewantowicz, Reference Johnson and Lewantowicz1990) and (Cunningham and Lewantowicz, Reference Cunningham and Lewantowicz1988) has been conducted in which GPS and INS interact through state dynamics. The GPS component was designed primarily as a discrete-time positioning system, leading to limited performance in estimating the system dynamic states (velocity and acceleration). Furthermore, GPS measurement errors were not rigorously accounted for. In (Moafipoor and Brzezinska, Reference Moafipoor, Brzezinska and Toth2004) and (Farrel, Reference Farrell2002), time-differential carrier phase was introduced to achieve more accurate INS calibration, in particular for the velocity error. This approach has some practical issues in that the absence of any single carrier phase measurement (for example due to a temporary loss of phase lock) leads to the loss of two time-differenced carrier phase measurements. However, more importantly, it did not address the time correlations of tropospheric, ionospheric and carrier-phase multipath errors. Instead it was assumed these errors were to be cancelled or uncorrelated in the time-differenced carrier-phase measurements. Unfortunately, this assumption cannot be held true for applications demanding high integrity, such as precision approaches, and its use may affect the integrity of the navigation solution, even under fault free conditions. Reference (Jekeli, Reference Jekeli2001) has suggested models for atmospheric delays in GPS/INS integrated systems, and illustrates the effect of using white noise and first order Gauss Markov processes to model ionospheric delays. Reference (Scherzinger, Reference Scherzinger2000) implements a popular tightly coupled scheme for GPS/INS navigation system with a focus on solving integer ambiguities without modelling any atmospheric decorrelation using state augmentation. The tightly coupled GPS/INS navigation system in (Ko, Reference Ko2000) includes the code multipath in the system state vector to explore the integrated system performance with various measurement update schemes (the carrier-smoothed codes, the code and carrier update and the carrier-phase riding) for Local Area Augmentation System (LAAS). However, the spatial decorrelation in GPS measurements between reference station and user receivers was still not dealt with rigorously.
To address these issues, this paper first presents the dynamic integration of GPS and INS that uses both pseudorange and carrier phase measurements and rigorously accounts for temporal and spatial decorrelation of GPS measurement errors (such as tropospheric, ionospheric and multipath errors) through state augmentation. The realization for the unified GPS/INS dynamic system is designed to ensure Kalman filter integrity under fault-free error conditions. Further, time-differenced measurements are not needed in the integrated filter, so the absence of a given single GPS measurement does not preclude the use of time-adjacent measurements. INS calibration and GPS positioning are simultaneously achieved using normal pseudorange and carrier-phase measurements (which capture information of relative motion over time) and the knowledge of accumulated correlations between GPS and INS states (which are implicitly accounted for in the propagation of the integrated system dynamics). In addition, cycle ambiguity states are explicitly present in the state realization, which can enable seamless interaction with external integer fixing algorithms and quick recovery of floating cycle ambiguity estimates for satellites that have been newly acquired or re-acquired.
In addietion to the fault-free integrity and implementation issues addressed above, this paper also evaluates the ability of simple fault detection functions, based on the filter innovations, to help provide integrity under certain fault conditions. For example, it is shown that the EKF innovations vector, a natural element of the integrated estimator, can provide high-integrity detection of abrupt carrier-phase cycle slips. These events can occur frequently when GPS signals are subject to interference or jamming. Furthermore, a seldom-used detection algorithm, based on the integral of the EKF innovations, is evaluated against other existing GPS fault detection methods. It is shown that the algorithm outperforms existing methods in the detection of slowly developing ranging errors, such as those caused by ionospheric and tropospheric anomalies.
In this paper, a differential GPS (DGPS) architecture is adopted as a framework for the development. The application specifically targeted is the Local Area Augmentation System (LAAS), which is a high integrity DGPS navigation system being developed for aircraft precision approach and landing (RTCA, 2004). All GPS-related states are derived for relative (reference station to aircraft) DGPS error models. Extension to other widely used GPS implementations, such as the Wide Area Augmentation System (WAAS) or stand-alone GPS, can be achieved by modifying the GPS error models and coordinate systems appropriately.
2. Individual GPS and INS system error models
2.1. Single Frequency DGPS System with Code and Carrier Measurements
In order to integrate GPS and INS system dynamics, a DGPS system model suitable for time propagation is needed. Currently, GPS signals on the L1 frequency (1575·42 MHz) are the only ones that are protected by Aeronautical Radio Navigation Services (ARNS) and readily accessible to the civil user. Therefore, only L1 code-and-carrier DGPS measurements are used in the development.
A newly developed DGPS error model with the appropriate time-propagating structure is adopted from (Pervan and Chan, Reference Pervan, Chan, Gebre-Egziabher, Pullen, Enge and Colby2003) and is reviewed briefly here:
where δiΔρk is the single difference L1 code measurement error between user and reference station for satellite i at time k. δΔxk is the error in the relative (user minus reference) position vector, and iek is the satellite line-of-sight (los) vector. For the later integration of these measurements with an inertial navigation system, GPS positions and los vectors are most conveniently expressed in the north-east-down (NED) local-level coordinates at the location of the reference station. Other terms including δΔb u,k, δiΔT k, δiΔI k,δiΔm ρ,k and iΔερ,k are single difference receiver clock error, tropospheric and ionospheric spatial decorrelation errors (between user and reference), single difference code multipath error and single difference receiver code-tracking noise, respectively. Among these error sources, the receiver clock error and the tropospheric and ionospheric decorrelation errors exist in the single difference carrier-phase measurement iΔφ kλ (metre) as well, where λ is the wavelength of L1 carrier. The magnitudes of carrier multipath iΔm φ,k and receiver carrier-tracking noise iΔεφ,k are much smaller compared with their code counterparts. In addition, ionospheric decorrelation has an effect on the carrier that is opposite in sign relative to the code. Finally, an additional integer cycle ambiguity iΔN appears in the carrier-phase measurement.
The relative ionospheric and tropospheric delays can be approximated using the LAAS decorrelation models in (McGraw and Murphy, Reference McGraw, Murphy, Brenner, Pullen and Van Dierendonck2000).
Tropospheric decorrelation model:
where δΔn is the residual refractivity uncertainty, h 0 is the troposphere scale height, h(k) is the user height at time k, E(i, k) is the elevation angle in radians for satellite i at time k. The first-order approximated tropospheric decorrelation error modelled by (2) is written in simplified form as the product of decorrelation factor iΔT dc, k and refractivity uncertaintyδΔn.
Ionospheric decorrelation model:
where iΔV ig is the vertical ionospheric gradient (VIG), typically expressed in mm/km, x air,h(k) is the horizontal distance between the user and reference station at time k. R E is the earth radius, and h I is the ionospheric shell height (350km). Similar to the simplified notation for tropospheric decorrelation, the ionospheric decorrelation error can also be approximated to first order by the product of the decorrelation factor iΔI dc,k and the vertical ionospheric gradient error δiΔV ig.
Combining (1)–(3), a measurement error model describing code and carrier measurement errors for one satellite can be formed:
![\hskip-12\eqalign{ \left[ {\matrix{ \delta\hskip3 {^{i} \rmDelta \rho _{k} } \cr \delta \hskip3 {^{i} \rmDelta \phi _{k} \lambda } \cr} } \right] \equals \tab \left[ {\matrix{ \minus {^{i} \varepsilon _{k} } \tab 1 \tab {^{i} \rmDelta T_{dc\comma k} } \tab {^{i} \rmDelta I_{dc\comma k} } \tab 0 \tab 1 \tab 0 \cr \minus {^{i} \varepsilon _{k} } \tab 1 \tab {^{i} \rmDelta T_{dc\comma k} } \tab \minus {^{i} \rmDelta I_{dc\comma k} } \tab \lambda \tab 0 \tab 1 \cr} } \right] \cr \tab \times \left[ {\matrix{ {\delta \rmDelta x_{k} } \tab {\delta \rmDelta b_{u\comma k} } \tab {\delta \rmDelta n} \tab {\delta ^{i} \rmDelta V_{ig} } \tab {\delta ^{i} \rmDelta N} \tab {\delta ^{i} \rmDelta m_{\rho \comma k} } \tab {\delta ^{i} \rmDelta m_{\phi \comma k} } \cr} } \right]^{T} \cr \tab \plus \left[ {\matrix{ {^{i} \rmDelta \varepsilon _{\rho \comma k} } \cr {^{i} \rmDelta \varepsilon _{\phi \comma k} } \cr} } \right] \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn4.gif?pub-status=live)
If there are n satellites in view, the complete set of estimated DGPS state errors is:
The complete measurement error model corresponding to the errors in the estimated DGPS state SG in (5) can be expressed as:
where δZk, Hk and εG,k are defined as:
![\delta {\bf Z}_{k} \equiv \left[ {\matrix{ {\delta {\bf Z}_{\rmDelta \rho \comma k} } \cr {\delta {\bf Z}_{\rmDelta \phi \comma k} } \cr} } \right]\comma {\bf H}_{k} \equiv \left[ {\matrix{ {{\bf H}_{\rmDelta \rho \comma k} } \cr {{\bf H}_{\rmDelta \phi \comma k} } \cr} } \right]\comma \bmepsiv _{G\comma k} \left[ {\matrix{ {{\bmepsiv}_{\rmDelta \rho \comma k} } \cr {{\bmepsiv} _{\rmDelta \phi \comma k} } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn7.gif?pub-status=live){kind=small}
and the measurement error vectors δZΔρ,k,δZΔφ,k and the corresponding observation matrices HΔρ,k, HΔφ,k in the (7) above are:
![{\bf H}_{\rmDelta \rho \comma k} \equiv \left[ {\matrix{ { \minus ^{\setnum{1}} {\bf e}_{k} } \tab 1 \tab {^{\setnum{1}} \rmDelta T_{dc\comma k} } \tab {^{\setnum{1}} \rmDelta I_{dc\comma k} } \tab {} \tab {{\bf 0}_{u\lpar n \minus \setnum{1}\rpar } } \tab {} \tab {} \tab \cr \vdots \tab \vdots \tab \vdots \tab \tab \ddots \tab \tab {{\bf 0}_{n} } \tab {{\bf I}_{n} } \tab {{\bf 0}_{n} } \cr {\minus ^{n} {\bf e}_{k} } \tab 1 \tab {^{n} \rmDelta T_{dc\comma k} } \tab {{\bf 0}_{L{\rm \lpar }n \minus\setnum{1}{\rm \rpar }} } \tab \tab {^{n} \rmDelta I_{dc\comma k} } \tab \tab \tab \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn9.gif?pub-status=live)
![{\bf H}_{\rmDelta \phi \comma k} \equiv \left[ {\matrix{ { \minus ^{\setnum{1}} {\bf e}_{k} } \tab 1 \tab {^{\setnum{1}} \rmDelta T_{dc\comma k} } \tab { \minus ^{\setnum{1}} \rmDelta I_{dc\comma k} } \tab \tab {{\bf 0}_{u\lpar n \minus \setnum{1}\rpar } } \tab \tab \tab \cr \vdots \tab \vdots \tab \vdots \tab \tab \ddots \tab \tab {\lambda {\bf 0}_{n} } \tab {{\bf 0}_{n} } \tab {{\bf 0}_{n} } \cr { \minus ^{n} {\bf e}_{k} } \tab 1 \tab {^{n} \rmDelta T_{dc\comma k} } \tab {{\bf 0}_{L{\rm \lpar }n \minus \setnum{1}{\rm \rpar }} } \tab \tab { \minus ^{n} \rmDelta I_{dc\comma k} } \tab \tab \tab \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn10.gif?pub-status=live)
where 0U(n-1) and 0L(n-1) are the upper and lower triangular part of a size n-1 square zero matrix, and In is an n×n identity matrix (to save some space, a square matrix is noted by its dimension only once. i.e. using In instead of In ×n ). The code and carrier receiver tracking noise components are stacked in the εΔρ,k and εΔφ,kn×1 noise vectors.
Widely used models for DGPS error dynamics assume that the local refractivity and vertical ionospheric gradients errors are constant during an aircraft approach, and that the code and carrier multipath are adequately described by first order Gauss-Markov Random Processes (GMRP).
The system error dynamic model can be written in continuous time as:
The dynamic matrix FG and the process noise vector WG are:
![{\bf F}_{G} \equiv \left[ {\matrix{ {{\bf 0}_{\lpar \setnum{5} \plus \setnum{2}n\rpar } } \tab {{\bf 0}_{\lpar \setnum{5} \plus \setnum{2}n\rpar } { _{\times} { _n}} } \tab {{\bf 0}_{\lpar \setnum{5} \plus \setnum{2}n\rpar } _{ {_\times} {_n}} } \cr {{\bf 0}_{n \times \lpar \setnum{5} \plus \setnum{2}n\rpar } } \tab {\minus {\bf I}_{n} \sol \tau _{\rmDelta \rho } } \tab {{\bf 0}_{n} } \cr {{\bf 0}_{n \times \lpar \setnum{5} \plus \setnum{2}n\rpar } } \tab {{\bf 0}_{n} } \tab {\minus I_{n} \sol \tau _{\rmDelta \phi } } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn12.gif?pub-status=live)
![{\bf W}_{G} \equiv \left[ {\lpar {{\bf C}_{n}^{n_{r} } {{\bf v}_{e } \rpar ^{\rm T} \ f_{cik}\ {\bf 0}_{\setnum{1} \times \lpar \setnum{1}\plus \setnum{2}n\rpar }\ {\bf n}_{\rmDelta \rho } ^{\rm T}\ {\bf n}_{\rmDelta \phi } ^{\rm T} } \right]^{\rm T}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn13.gif?pub-status=live){kind=small}
where τΔρ and τΔφ are the time constants and nΔρ and nΔφ are the driving white noise vectors in the GMRPs for single difference code and carrier multipath models respectively. ve is the user ground velocity, which is not available to the GPS-only user, and will be the basis of integration with INS in the later discussion of the subject along with the coordinate transformation matrix 
. fclk is the instantaneous receiver clock frequency offset.
In order to use a discrete-time extended Kalman filter to perform state propagation and estimation, a continuous-to-discrete conversion of the DGPS system error dynamics (11) can be done in a straightforward manner, and the discrete-time DGPS system error time propagation is expressed in simplified notation as:
Throughout this error analysis, the state estimation errors at time k, δSG,k, are related to the true states, SG,k, by:
For Kalman filter implementation, the error propagation in time and the covariance of state estimation error can be formed as:
where 
, 
and 
are defined as:
and 
is the state estimate prior to the measurement at epoch k and 
is the state estimate after the measurement at epoch k-1.
After time propagation from epoch k-1 to k, the measurements at epoch k are available for EKF measurement update of the DGPS state vector estimate and the corresponding state covariance. The state covariance after the EKF measurement update is:
![{\bf \hat{P}}_{S_{G} \comma k} \equals \lpar {\bf I} \minus {\bf K}_{k} {\bf H}_{k} \rpar { \bar{\bf P}}_{S_{G} \comma k} \lpar {\bf I} \minus {\bf K}_{k} {\bf H}_{k} \rpar ^{\rm T} \plus {\bf K}_{k} {\bmSigma}_{\varepsilon _{G} } {\bf K}_{k} ^{\rm T}]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn18.gif?pub-status=live){kind=small}
where Kk is the Kalman filter gain and 
is the covariance matrix of the errors on available measurements. They are defined as:
![\eqalign{ \tab {\bf K}_{k} \equiv { \bar{\bf P}}_{S_{G} \comma k} {\bf H}_{k} ^{\rm T} \lpar {\bmSigma}_{\varepsilon _{G} }\minus {\bf H}_{k} { \bar{\bf P}}_{S_{G} \comma k} {\bf H}_{k} ^{\rm T} \rpar ^{ \minus \setnum{1}} \cr \tab {\bmSigma}_{{\bmepsiv}_{G}} \equiv {\rm E}\left[ {\bmepsiv}_{G\comma k \minus \setnum{1}} {\bmepsiv}_{G\comma k \minus \setnum{1}} ^{\rm T} } \right] \equals {\rm E}\left[ {\bmepsiv _{G\comma k} {\bmepsiv}_{G\comma k} ^{\rm T} } \right] \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn19.gif?pub-status=live)
2.2. INS Navigation
The general inertial navigation equation expressed in a local-level navigation frame, n (in north-east-down), is well known and briefly described in this section without detailed derivations (Titterton and Weston, Reference Titterton and Weston2004) (Jekeli, Reference Jekeli2001):
where 
is the time derivative of the ground velocity (ve) with respect to the navigation frame n (denoted by the left superscript) and expressed in navigation frame n coordinates (denoted by right superscript within parentheses), f(n) is the specific force expressed in navigation frame n coordinates. Similarly, wie(n), wen(n) and g1(n) are all expressed in navigation frame n coordinates, where wie(n) is the earth's rotation rate; wen(n) is the angular velocity of the navigation frame relative to the earth (also known as transport rate); g1(n) is the local gravity vector.
Euler angle propagation is adopted as the basic mechanism for vehicle attitude propagation, which is sufficient for analysis purposes:
![{ \dot { \bmPsi}} \equals {\bf F}_{EU} {\bf w}_{nb}^{\lpar b\rpar } \equals {\bf F}_{EU} \left[ {{\bf w}_{ib}^{\lpar b\rpar } \minus {\bf C}_{n}^{b} \left( {{\bf w}_{ie}^{\lpar b\rpar } {\bf w}_{en}^{\lpar b\rpar } } \right)} \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn21.gif?pub-status=live){kind=small}
where Ψ is the Euler angle vector consisting of roll (Φ), pitch (θ), and yaw (Ψ) angles. FEU is the transformation matrix to translate the instantaneous body-to-navigation-frame angular velocity,wnb(b), into the Euler angle rotation rate 
(Titterton and Weston, Reference Titterton and Weston2004) (Jekeli, Reference Jekeli2001). wib(b) is the gyro measured angular velocity of the body relative to an inertial frame and Cn b is the transformation matrix from navigation frame to body frame based on the three Euler angle rotations.
The linearized inertial error dynamics equations can also be found in most inertial textbooks (Titterton and Weston, Reference Titterton and Weston2004) (Jekeli, Reference Jekeli2001); therefore they are shown here without further explanation:
![\delta { \dot { \bmPsi}} \equals \minus \left[ { {{\bf w} _{in}^{\lpar n\rpar } \times } \right]\delta {\bmpsi} \minus {\bf C}_{b}^{n} \delta {\bf w}_{ib}^{\lpar b\rpar } \plus \delta {\bf w}_{in}^{\lpar n\rpar }](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn23.gif?pub-status=live){kind=small}
where win(n) is the navigation frame angular rate relative to an inertial frame. It can be decomposed as the sum of wie(n) and wen (n) as is done in (21). The error in the earth rotation rate is usually so small for general inertial navigation applications that it can be ignored. Furthermore, the error on the transport rate wen (n) can be expressed in terms of user ground velocity error as:
![\delta {\bf w}_{e n}^{\lpar n\rpar } \equals \left[ {\matrix{ 0 \tab {\displaystyle{1 \over {\lpar R_{P} \plus h\rpar }}} \tab 0 \cr {\displaystyle{{ \minus 1} \over {\lpar R_{M} \plus h\rpar }}} \tab 0 \tab 0 \cr 0 \tab {\displaystyle{{ \minus \tan \lpar l_{{\rm at}} \rpar } \over {\lpar R_{M} \plus h\rpar }}} \tab 0 \cr} } \right]\left[ {\matrix{ {\delta v_{N} } \cr {\delta v_{E} } \cr {\delta v_{D} } \cr} } \right] \equiv {\bf F}_{V\setnum{2}T} \delta {\bf V}_{e }^{\lpar n\rpar }](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn24.gif?pub-status=live)
where R P is the prime radius of curvature of the earth, R M is the meridian radius of curvature, and h is the local ellipsoidal height.
Substituting (24) into (23), the equations for attitude error dynamics can be re-written as:
![\delta { \dot { \bmPsi}} \equals \minus \left[ {{\bf w}_{in}^{\lpar n\rpar } \times } \right]\delta {\bmPsi} \minus {\bf C}_{b}^{n} \bi w_{ib}^{\lpar b\rpar } \plus {\bf F}_{V\setnum{2}T} \delta {\bf v}_{e }^{\lpar n\rpar }](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn25.gif?pub-status=live){kind=small}
Inertial sensor errors generally propagate into navigation system errors through δf(b) in (22) and δwib(b) in (23), which include IMU measurement bias error δba/g, misalignment error δmis,a/g, scale factor error δsf,a/g and measurement noise εa/g (subscript a represents accelerometer while subscript g represents gyroscope). To be consistent with the sign of state error definition in (15), IMU sensor errors are defined as:
At this point, instead of pursuing optimal system performance for one specific type of IMU error model, a generalized IMU error model, using first order GMRP models to bound the inertial sensor bias variations and process noise to handle other residual errors, is adopted in this work. Sensitivity analysis has shown that for a navigation grade IMU, modelling scale factor and misalignment errors as white process noise inputs rather than initial unknown biases in state, which is the actual physical mechanism, shows almost no difference in the system performance. The benefit of state reduction is the reason for adoption of this method for handling scale factors and misalignment in this work. Equation (27) expresses the dynamics of inertial bias state errors:
where τg and τa are the time constants of GMRP models for gyroscope and accelerometer biases respectively. ng and na are the corresponding driving white noise vectors.
Another important factor affecting INS performance is the accuracy of local gravity information. A popular model (Titterton and Weston, Reference Titterton and Weston2004) (Gebre-Egziahber, Reference Gebre-Egziahber2001) constituting of the deflections of vertical (DOV) and vertical gravity anomaly is used in this work and the error analysis is shown as:
![\delta {\bf g}_{l}^{\lpar n\rpar } \equals \left[ \openup3{\matrix { 0 \cr\vskip-2 0 \cr {{\displaystyle{{ \minus 2g_{\setnum{0}} } \over {R_{\setnum{0}} }}}} \cr} } \right]\delta h \plus \left[ {\matrix{ { \minus g_{\setnum{0}} } \tab 0 \tab 0 \cr 0 \tab { \minus g_{\setnum{0}} } \tab 0 \cr 0 \tab 0 \tab { \minus 1} \cr} } \right]\left[ {\matrix{ {\delta \xi } \cr {\delta \eta } \cr {\delta \rmDelta g_{\setnum{0}} } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn28.gif?pub-status=live)
where δh is the vertical position error, δξ and δη are the DOV errors in north-south and east-west directions respectively. g 0 and R 0 are the nominal values of local gravity and the earth's radius. δΔg 0 is the error in the local gravity magnitude. Simplified notations for the vectors and matrices in (28) are as adopted as follows:
![{\bf F}_{h\setnum{2}V} \equiv \left[ {\matrix{ 0 \tab 0 \tab { \minus 2g_{\setnum{0}}\sol R_{\setnum{0}} } \cr} } \right]^{\rm T} \comma {\bf F}_{gm\setnum{2}V} \equiv \left[ {\matrix{ { \minus g_{\setnum{0}} } \tab 0 \tab 0 \cr 0 \tab { \minus g_{\setnum{0}} } \tab 0 \cr 0 \tab 0 \tab \minus1 \cr} } \right]{\rm \ and\ }\delta {\bf g}_{gm} \equiv \left[ {\matrix{ {\delta \xi } \tab {\delta \eta } \tab {\delta \rmDelta g_{\setnum{0}} } \cr} } \right]^{\rm T}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqnU10.gif?pub-status=live)
Equation (28) can now be re-written as:
Spatial variations of the gravity anomalies are handled in a similar way to inertial bias variations. A first order GMRP is used to model the variations of the DOVs and vertical gravity anomaly. The dynamics of modelled gravity errors can be illustrated as:
where τgm is the time constant of GMRP for the gravity error model, which is obtained by dividing the correlation distance of gravity anomaly by vehicle velocity. ngm is the vector of driving white noise.
Collecting (22), (25), (27) and (30), the dynamics of INS errors can be formed as:
![\left[ {\matrix{ {{}^{n}\delta { \dot{\bf v}}_{e}^{\lpar n\rpar } } \cr {\delta {\dot\bmPsi} } \cr {\delta {\dot{\bf b}}_{g} } \cr {\delta {\dot{\bf b}}_{a} } \cr {\delta { \dot{\bf g}}_{gm} } \cr} } \right] \equals \left[ {\openup1\matrix{ {{\bf 0}_{\setnum{3}} } \tab {{\bf f}^{\lpar n\rpar } \times } \tab {{\bf 0}_{\setnum{3}} } \tab { \minus {\bf C}_{b}^{n} } \tab {\bf F} _{gm\setnum{2}V} \cr {{\bf F}_{V\setnum{2}T} } \tab { \minus {\bf w}_{in}^{\lpar n\rpar } \times } \tab {\bf C}_{b}^{n} \tab {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \cr {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab { \minus {\displaystyle{{{\bf I}_{\rm \setnum{3}} } \over {\tau _{g} }}}} \tab {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \cr {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab { \minus {\displaystyle{{I_{\rm \setnum{3}} } \over {\tau _{a} }}}} \tab {{\bf 0}_{\setnum{3}} } \cr {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab { \minus {\displaystyle{{{\bf I}_{\setnum{3}} } \over {\tau _{gm} }}}} \cr} } \right]\left[ {\matrix{ {\delta {\bf v}_{e}^{\lpar n\rpar } } \cr {\delta {\bmpsi}} \cr {\delta {\bf b}_{g} } \cr {\delta {\bf b}_{a} } \cr {\delta {\bf g}_{gm} } \cr} } \right] \plus W_{I}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn31.gif?pub-status=live)
where
![W_{I} \equals \left[ {\matrix{ {{\bf C}_{b}^{n} \lpar \minus {\bmdelta} {\bf s}_{f\comma a} \minus {\bmdelta } {\bf m}_{is\comma a} \minus {\bmepsiv}_{a} \rpar \plus {\bf F}_{{\rm h}\setnum{2}V} \delta h} \cr { \minus {\bf C}_{b}^{n} \lpar \minus {\bmdelta}{\bf s}_{f\comma g} \minus {\bmdelta} {\bf m}_{is\comma g} \minus \rmvarepsilon _{g} \rpar } \cr {{\bf n}_{g} } \cr {{\bf n}_{a} } \cr {{\bf n}_{gm} } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqnU11.gif?pub-status=live)
Using simplified notation, equation (31) can be expressed as:
3. GPS/INS integration in the dynamic model
The time derivative of GPS position with respect to the ground is the user ground velocity, which is exactly the information derived from INS. Therefore, the error dynamics of the GPS positions is derived:
where the n r frame is the local-level navigation frame at the location of the reference station used to generate the unit line-of-sight vectors in equation (1).
The GPS system dynamics is connected to INS through (33) with a coordinate transformation 
:
The coordinate transformation is necessary for a mathematically rigorous derivation, yet it can be ignored in actual implementation if the distance between the user and reference station is small. The direct connections of GPS and INS through system dynamics in (34) and also through gravity compensation error due to height in (29) result in strong correlations across GPS and INS states. This bonding among Kalman filter states is the foundation that simultaneously enables excellent INS calibration and GPS positioning at the same time through the EKF measurement update process. Fusing the dynamics of the integrated system in the continuous-time domain is accomplished by combining (11) and (32), and moving ve from WG and Fh2Vδh from WI to the state dynamics as specified by (32) and (29) respectively. The resulting centralized EKF state error dynamics are described by:
![\left[ {\matrix{ {\delta { \dot{\bf S}}_{I} } \cr {\delta { \dot{\bf S}}_{G} } \cr} } \right] \equals \left[ {\matrix{ {{\bf F}_{I} } \tab {{\bf F}_{G\setnum{2}I} } \cr {{\bf F}_{I\setnum{2}G} } \tab {{\bf F}_{G} } \cr} } \right]\left[ {\matrix{ {\delta {\rm S} _{I} } \cr {\delta {\rm S} _{G} } \cr} } \right] \plus \left[ {\matrix{ {{\bf W}_{I\setnum{1}} } \cr {{\bf W}_{G\setnum{1}} } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn35.gif?pub-status=live){kind=small}
where FI2G, FG2IWI1and WG1 are defined as:
![{\bf F}_{G\setnum{2}I} \equiv \left[ {\matrix{{ {\bf 0}_{\setnum{3}} } \tab \cdots \tab \cdots \cr {{\bf 0}_{\setnum{3}} } \tab \cdots \tab \cdots \cr {{\bf 0}_{\setnum{3}} } \tab \cdots \tab \cdots \cr {{\bf 0}_{\setnum{3}} } \tab \cdots \tab \cdots \cr {{\bf F}_{X\setnum{2}V} {\bf C}_{n}^{n_{r} ^{{\rm T}} } } \tab {{\bf 0}_{\setnum{3}} } \tab \cdots \cr} } \right]_{\setnum{15} \times \lpar \setnum{5} \plus \setnum{4}n\rpar } {\bf F} _{X\setnum{2}V} \equiv \left[ {\matrix{ {{\bf 0}_{\setnum{1} \times \setnum{3}} } \cr {{\bf 0}_{\setnum{1} \times \setnum{3}} } \cr { \minus {\bf F}_{h\setnum{2}V}^{\rm T} } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqnU13.gif?pub-status=live)
![\hskip-87{\bf W}_{I\setnum{1}} \equiv \left[ {\matrix{ {C_{b}^{n} \lpar \minus\delta s_{f\comma a} \minus \delta m_{is\comma a} \minus \varepsilon _{a} }\rpar \cr { \minus C_{b}^{n} \lpar \minus \delta s_{f\comma g} \minus \delta m_{is\comma g} \minus \varepsilon _{g} } \rpar \cr {n_{g} } \cr {n_{a} } \cr {n_{g m}} \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn36.gif?pub-status=live)
![{\bf F}_{I\setnum{2}G} \equiv \left[ {\matrix{ {{\bf C}_{n}^{n_{r} } } \tab {{\bf 0}_{\setnum{3}} } \tab \cdots \cr 0 \tab {{\bf 0}_{\setnum{1} \times \setnum{3}} } \tab \cdots \cr {{\bf 0}_{\setnum{3}} } \tab {{\bf 0}_{\setnum{3}} } \tab \cdots \cr \vdots \tab \vdots \tab \ddots \cr} } \right]_{\lpar \setnum{5} \plus \setnum{4}n\rpar \times \setnum{15}} {\bf W}_{G\setnum{1}} \equiv \left[ {\matrix{ {{\bf 0}_{\setnum{3} \times \setnum{1}} } \cr {f_{clk} } \cr {{\bf 0}_{\lpar \setnum{1} \plus \setnum{2}n\rpar \times \setnum{1}} } \cr {{\bf n}_{\rmDelta \rho } } \cr {{\bf n}_{\rmDelta \phi } } \cr} } \right]](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqnU14.gif?pub-status=live)
and n is the number of satellites in view.
The continuous-time system error dynamics (35) of the centralized integrated GPS/INS EKF can be converted into discrete-time system error propagation equations for a given IMU sampling interval Δt. Generally, the IMU sampling interval Δt is much smaller than the time span between GPS measurement time t k at epoch k and measurement time t k-1 at epoch k-1 (Δt⩽t k−t k-1):
![\hskip-9\eqalign{ \displaylines\tab { \left[ {\matrix{ {\delta {\bar{\bf S}}_{I\comma t_{{k \minus \setnum{1}}} \plus \rmDelta t} } \cr {\delta {\bar{\bf S}}_{G\comma t_{{k \minus \setnum{1}}} \plus \rmDelta t} } \cr} } \right] \equals \left[ {\matrix{ {{\bmPhi}_{I\comma k \minus \setnum{1}} } \tab {\bmPhi}_{G\setnum{2}I\comma k \minus \setnum{1}} \cr {\bmPhi _{I\setnum{2}G\comma k \minus \setnum{1}} } \tab {\bmPhi}_{G\comma k \minus \setnum{1}} } \cr} } \right]\left[ {\matrix {\delta {{\bf \hat{S}}_{I\comma k \minus \setnum{1}} } \cr {\delta {\bf \hat{S}}_{G\comma k \minus \setnum{1}} } \cr} } \right] \plus \left[ {\matrix{ {{\bf W}_{I\setnum{1}\comma k \minus \setnum{1}} } \cr {{\bf W}_{G\setnum{1}\comma k \minus \setnum{1}} } \cr} } \right]} \cr \tab\hskip10pc\downarrow \cr\tab \lsqb {\rm propagation\ from\ }t_{k \minus \setnum{1}} {\rm \ to\ }t_{k} {\rm \ with\ sampling\ interval\ \rmDelta }t{\rm \rsqb } \cr \tab\hskip10pc\downarrow \cr \tab\left[ {\matrix{ {\delta {\bar{\bf S}}_{I\comma k} } \cr {\delta {\bar{\bf S}}_{G\comma k} } \cr} } \right] \equals \left[ {\matrix{ {{\bmPhi}_{I\comma t_{k} \minus \rmDelta t} } \tab {{\bmPhi}_{G\setnum{2}I\comma t_{k} \minus \rmDelta t} } \cr {\bmPhi}_{I\setnum{2}G\comma t_{k} \minus \rmDelta t} \tab {{\bmPhi}_{G\comma t_{k} \minus \rmDelta t} }} } \right]\left[ {\matrix{ {\delta {\bar{\bf S}}_{I\comma t_{k} \minus \rmDelta t} } \cr {\delta {\bar{\bf S}}_{G\comma t_{k} \minus \rmDelta t} } \cr} } \right] \plus \left[ {\matrix{ { {\bf W}_{I\setnum{1}\comma t_{k} \minus \rmDelta t} } \cr { {\bf W}_{G\setnum{1}\comma t_{k} \minus \rmDelta t} } \cr} } \right] \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn37.gif?pub-status=live)
When the measurements at epoch k are available, the EKF measurement update can be performed to calibrate INS error states and estimate GPS positions by a simple measurement model, equivalent to that used if only the DGPS position is estimated:
![\delta {\bf Z}_{k} \equals \lsqb {\bf 0}_{\setnum{2}n \times \setnum{15}} \quad {\bf H}_{k} \rsqb \left[ {\matrix{ {\delta {\bar{\bf S}}_{I\comma k} } \cr {\delta { \bar{\bf S}}_{G\comma k} } \cr} } \right] \plus {\bmepsiv} _{G\comma k}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn38.gif?pub-status=live){kind=small}
where δZk, Hk and εG,k are defined in (7)–(10).
Using the system error dynamics and measurement error models from (37) and (38), the time propagation and the measurement update of the integrated system covariance can be performed as described for the DGPS EKF process in (16)–(19) with the corresponding transition matrix in (37), and observation matrix and measurement error in (38).
INS velocity and attitude calibration performance was quantified using covariance analysis for a nominal LAAS straight-in approach at O'Hare international airport in Chicago. The values used for nominal navigation grade IMU and gravity error models are displayed in Table 1. The worst GPS satellite geometry, defined by the largest position dilution of precision (PDOP) observed using the DO-229D constellation (RTCA, 2006) for 24 hours, was selected for the covariance analysis. Figure 1 shows the performance of INS calibration with a nominal GPS model displayed in Table 2, in which the aircraft launched the precision approach 37 km away from the runway end, and the INS was calibrated all the way through landing. The results show that less than 3 mm/sec standard deviation of calibrated velocity errors and less than 0·01 degree standard deviation of calibrated roll and pitch angle errors can be achieved within a 20 second calibration period after system initialization (about 1·5 km travelling distance in Figure 1). The standard deviation of the azimuth angle error can be calibrated down to less than 0·1 degree after a longer period of approximately 4 min (about 18·5 km travelling distance). This velocity estimation performance is comparable to that obtained using time-differenced carrier-phase measurements in (Moafipoor and Brzezinska, Reference Moafipoor, Brzezinska and Toth2004) and (Farrel, Reference Farrell2002) (which show standard deviation in the range of 2–5 mm/sec). However in the case of this filter realization, fault-free integrity of the actual state estimate is ensured because temporal and spatial correlations of GPS measurement errors have been properly accounted for. Velocity estimation is implicitly achieved through the system dynamics. Further, as noted earlier, the potential for unnecessary measurement loss due to time-differencing is eliminated. Not surprisingly, these results show superior performance relative to more conventional complementary KF implementations, for example (Marty and Pagnucco, Reference Marty and Pagnucco1992) and (Soltz and Donna, Reference Soltz, Donna and Greenspan1988), which use time-differenced pseudorange to obtain sub-metre level velocity estimation accuracy.
The nominal velocity (left) and attitude (right) calibration performance of the integrated navigation system with error model parameters in Table 1.
Parameter values of nominal navigation grade IMU and gravity error models.
(hr=hour, ARW=angle random walk, μg=micro g, nmi=nautical mile, std=standard deviation, E-W=east – west direction, N-S=north – south direction, d τ=τgm·ν h where νh is horizontal ground speed.)
Parameter values of the nominal GPS model.
4. system validation
The mathematical models and the EKF implementation derived in the previous section for the high-integrity tightly-coupled GPS/INS integration were validated using archived flight data provided by the FAA William J. Hughes Technical Center. The data consists of GPS and IMU measurements recorded over eight approaches during a LAAS flight test at Memphis International Airport on Sep. 28th 2006.
Raw GPS measurements in the archived flight data were recorded with a 5 Hz output rate by two L1/L2 capable NovAtel OEM-4 receivers, one on an N47 aircraft and the other at a ground reference station on the airport. To reconstruct truth trajectories for system performance evaluation, a carrier-phase double-difference (DD) wide-lane position history was generated for the aircraft. In this regard, ionosphere-free carrier-phase DD wide-lane integers were estimated and fixed by filtering wide-lane carrier plus narrow-lane code (Heo and Pervan, Reference Heo and Pervan2006). Given the fixed DD wide-lane integers, carrier-phase DD wide-lane positioning was executed using least-squares estimation, resulting in a position error at or below the 10-cm level. Figure 2 shows all eight trajectories using the carrier-phase DD wide-lane DGPS positioning. The LGF antenna is the reference antenna for all DGPS data.
True trajectories of eight approaches in flight data reconstructed by carrier-phase DD wide-lane DGPS positioning.
To ensure that the DD wide-lane carrier integers were fixed correctly, and that the reconstructed ‘truth’ trajectory is correct, carrier-phase DD wide-lane measurement residuals were computed. These residuals were generated by removing the estimated geometric ranges, which are the projections of the computed ‘true’ positions onto the line-of-sight vectors, and the fixed wide-lane integers from corrected carrier-phase DD wide-lane measurements. A number of smaller corrections to the differential carrier phase measurements were also applied as specified in (Lawrence, Reference Lawrence1996). As long as the removed position projections and integers are correct, the main components in the residuals should be DD ionospheric spatial decorrelation and DD wide-lane carrier-phase multipath. The DD ionospheric residual error under normal ionospheric activity is expected to be the dominant error source, with the magnitude of the standard deviation around 10 cm or less at the beginning of each approach, when the aircraft is far from the reference station; the ionospheric error is expected to reduce gradually from here while the carrier-phase multipath error remains at a similar magnitude all the way to the end of the runway. The standard deviation of DD wide-lane carrier-phase multipath error generally has a magnitude around 2 to 3 cm or less. Figure 3 shows the resulting residuals computed for one example approach. The behaviour of the residuals is consistent with expectations, and similar behaviour is observed for the other approaches. This is evidence that the reconstructed position trajectories represent a basis for truth within approximately 10 cm. Again, the purpose of these truth trajectories is to help validate the integrated navigation system.
Carrier-phase DD wide-lane measurement residual check for approach one.
Although the aircraft was also equipped with an onboard navigation grade IMU (Honeywell HG 2001GD), the archived raw data was of insufficient precision to be modelled by the typical navigation specifications in Table 1. To overcome this problem, stored yaw, pitch and roll angle output, computed directly by the aircraft inertial system, were used as a truth reference to re-evaluate the statistical gyro error model parameters for the raw angular rate measurements. The deviations between the stored reference attitude and the pure INS-only propagated Euler angle attitude (free-coasting attitude generated from the raw angular rate measurements) are presented in Figure 4. To account for the attitude propagation errors exhibited in Figure 4, the gyroscope error model parameters are selected so that the free-coasting attitude covariance propagation provides a reasonable match to the actual free-coasting attitude errors. These matched model parameters are listed in Table 3.
Deviations between the reference attitudes and the pure INS-only propagated attitudes for three Euler angles: roll δφ, pitch δθ and yaw δψ.
Parameter values of IMU and gravity error models.
The archived acceleration measurements were of insufficient quality to be treated in the same manner as in the gyro model re-evaluation above. However, in this case artificial acceleration data was generated from the reference ‘truth’ position and attitude histories. The reference position history was interpolated (re-sampled) with high frequency in the Latitude, Longitude, and Altitude (LLA) domain with simulated gravity anomalies added. The reference attitudes were then used to convert the simulated accelerations from the aircraft NED navigation frame to the aircraft body frame. Finally, simulated sensor errors were added into the body frame accelerometer measurements. The simulated sensor error models and gravity anomaly error models used are defined in Table 3. Following this process, a “partially-simulated” version of IMU flight data was produced. The free-coasting INS position errors and the corresponding coasting position standard deviations, generated by propagating the partially simulated IMU data, are shown in Figure 5 for one example approach. From this point on, whenever the IMU flight data is mentioned, the acceleration data contained therein are all obtained in this manner.
Deviations between the true positions and the pure INS-only propagated positions from approach one.
The performance of the GPS/INS integrated system was evaluated for a specific LAAS application, in which the INS is calibrated by GPS measurements when the aircraft enters LAAS service volume until reaching the final approach fix (FAF) point. After that, the calibrated INS becomes the sole navigation system leading the aircraft to the decision height (DH), which is specified to be 50 feet in this work, to eliminate the concern of vulnerability of the GPS signals to radio frequency interference (RFI) during the most critical precision landing stage.
Figure 6 shows the lateral and vertical INS coasting position errors at DH for all eight approaches. The predicted standard deviations of position errors stay reasonably close to the actual errors. The results support the performance analysis of the tightly-integrated system described in the previous section, and validate the covariance analysis by showing that the predicted errors are sensibly near the actual errors.
Final coasting position errors at DH in all approaches.
In addition, the results for position, velocity and attitude errors during an entire approach (without coasting) shown in Figures 7 and 8 further validate the covariance performance analysis of the GPS/INS integration architecture.
Position errors of GPS/INS integrated system during approach 2.
Velocity (left) and attitude (right) errors of GPS/INS hybrid system during approach 2.
5. Innovations-based detection
Given that precision approach is a potential application for the GPS/INS integrated navigation system, navigation integrity must be ensured. Considering that INS calibration performance is highly dependent on the quality of GPS carrier-phase measurements, one or more effective detection mechanisms to protect the integrity of carrier-phase measurements in real time are necessary. A carrier-phase cycle slip is the most common fault affecting carrier-phase GPS users, and will cause INS calibration errors if it happens undetected.
There are different ways to detect carrier-phase cycle slips, which include well-known Receiver Autonomous Integrity Monitoring (RAIM), carrier-phase Relative RAIM (RRAIM) (Heo and Pervan, Reference Heo, Pervan, Pullen, Gautier, Enge and Gebre-Egziabher2004), and direct comparison of carrier-phase measurements with other information, such as INS (Colombo and Bhapkar, Reference Colombo, Bhapkar and Evans1999) (Lee and O'Laughlin, Reference Lee and O'Laughlin1999). In general, additional processing is needed in order to detect cycle-slip errors. However, one advantage of the GPS/INS integration architecture is a built-in cycle-slip detection mechanism using the KF innovation vector (Groves, Reference Groves2008). The innovation vector is an intermediate product of the KF process, and the squared magnitude of its normalized carrier-phase elements provides a simple test statistic for the detection of abrupt changes in the carrier-phase measurements.
To investigate the general detection performance of built-in cycle slip detection, a half-cycle slip occurring on one satellite measurement in the middle of an approach was simulated. A half-cycle magnitude was chosen to simulate a difficult-to-detect scenario (because it is the smallest possible magnitude for the slip). The LAAS precision landing scenario described above is considered here again. In Figure 9, the integrated navigation system was initialized about 37 km away from the runway end (corresponding to the far right of each plot), which is the nominal LAAS service radius. The half-cycle slip was introduced at a distance of about 20 km away from the runway end, and the system outputs were affected immediately. The estimation errors caused by the cycle slip persist over time and did not disappear even when the aircraft reached FAF, where INS coasting begins. In this simulated cycle-slip case, with the smallest slip magnitude (1/2 cycle), the unsettled errors in position and velocity states could greatly affect the accuracy of INS propagation during the final approach and landing. This clearly demonstrates the necessity of system integrity monitoring on GPS measurements during INS calibration.
Simulated integrated system responses to a half-cycle slip occurring on one satellite carrier-phase measurement.
The test statistic of the built-in carrier-phase innovation detector, T s,inno(k), can be expressed mathematically as:
![\eqalign{ \tab{ T}_{s\comma inno} \lpar k\rpar \equiv \left\Vert {{\bf r}_{\rmphi \_inno\comma k} } \right\Vert_{norm}^{\setnum{2}} \sim \chi ^{\setnum{2}} \lpar n_{k} \rpar \cr \tab\left\Vert {{\bf r}_{\rmphi \_inno\comma k} } \right\Vert_{norm}^{\setnum{2}} \equals \left\Vert {{\bf r}^{\rm T}_{\rmphi \_inno\comma k{\rm \ } } {{\bf P}_{\rmphi \_inno\comma k}\hskip-3 ^{{{\rm \ } \minus \setnum{1}{\rm \ }}} } {\bf r}_{\rmphi \_inno\comma k} } \right\Vert^{\setnum{2}} \cr \tab\left[ {\matrix{ {{\bf r}_{\rho \_inno\comma k} } \cr {{\bi r}_{\rmphi \_inno\comma k} } \cr} } \right] \equals {\bf Z}_{k} \minus {\bf H}_{sys\comma k} {\bar{\bf S}}_{sys\comma k} \cr \tab \left[ {\matrix{ {{\bf P}_{\rho \_inno\comma k} } \tab {{\bf P}_{\rmphi \rho \_inno\comma k} } \cr {{\bf P}_{\rmphi \rho \_inno\comma k} } \tab {{\bf P}_{\rmphi \_inno\comma k} } \cr} } \right] \equals {\bf H}_{sys\comma k} { \bar{\bf P}}_{sys\comma k} {\bf H}_{sys\comma k} {\bmSigma} _{\varepsilon _{\rm G}} \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqn39.gif?pub-status=live)
where Zk, Hsys,k, 
and 
are defined according to (38):
![\eqalign{\tab {\bf H}_{sys\comma k} \equiv \lsqb {\bf 0}_{\setnum{2}n \times \setnum{15}} \quad {\bf H}_{k} \rsqb \comma \quad {\bar{\bf S}}_{sys\comma k} \equiv\left[ {\matrix{ {{\bar{\bf S}}_{I\comma k} } \cr { {\bar{\bf S}}_{G\comma k} } \cr} } \right] \cr \tab { \bar{\bf P}}_{sys\comma k} \equiv E\lsqb \delta {\bar{\bf S}}_{sys\comma k} \delta {\bar{\bf S}}_{sys\comma k} ^{T} \rsqb \cr}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022041624836-0158:S0373463311000038_eqnU17.gif?pub-status=live)
The test statistic T s,inno(k) has a normalized chi-squared probability density function with n k degrees of freedom (DOF), which is the number of available carrier-phase measurements at epoch k (T s,inno(k)~χ2(n k)). The detection threshold D th,inn o(k) is conservatively set to provide a false alarm probability of 2×10−6/30, which is based on the most stringent LAAS category III continuity risk requirement (2×10−6 for 30 sec) (RTCA, 2004). The thresholds for other detection algorithms, used for comparison purposes below, are defined according to the same requirement. The probability of false alarm (FA) under normal conditions (NC), P(FA|NC), is related to the detection threshold and allowable continuity risk as follows:
where the integral is the incomplete gamma function.
One example of the time responses to five consecutive half cycle slips on one SV from various existing error detection algorithms, which include code-phase RAIM, carrier-phase RRAIM, carrier-phase RRAIM with INS augmentation, are shown in the top graph of Figure 10 to be compared with the response of the built-in carrier-phase innovation detector, which is shown at the bottom of Figure 10. Every half cycle slip error (stepped up from 0 to 0·48 m in five steps) was detected with high confidence by the built-in detector and the corresponding carrier integer was reset accordingly. The detection performance of the built-in detector is obviously superior to others.
Simulated typical time responses of various error detection algorithms to half-cycle slip error on one satellite carrier-phase measurement.
Before other types of anomalies are considered, another detection algorithm is evaluated as well. It will be shown shortly that the algorithm has advantages in detecting slowly varying errors inherent in GPS measurements. The detection algorithm is based on the integral of the carrier-phase KF innovation – in other words, the summation over time of T s,inno(k). It will be called the innovation integration detection algorithm.
The test statistic of the innovation integration detection algorithm, T s,integral(k), also has a chi-squared probability density function since each individual element of the summation is fundamentally a chi-squared distributed random variable. The DOF of T s,integral(k) is the sum of the available carrier-phase measurements within the integration period. The mathematical expression for T s,integral(k) is:
The threshold is plotted in Figure 11 along with a Monte Carlo simulation of 1000 approaches under normal error conditions.
Monte Carlo simulation for the innovation integration detection algorithm with 1000 approaches under normal error conditions.
A time-of-fault-onset integrity analysis was performed for the built-in innovation and the innovation integration detectors for the half-cycle slip scenario considered above. In this analysis the error occurrence time (previously corresponding to 20 km distance) is varied. The probability of mis-detection, P md, of the innovation and the innovation integration detections are evaluated for all possible times of occurrence from system initialization to the final approach fix (FAF) point in 1 second increments. The probability of mis-detection is the product of the minimum probability of no-alarm (P NA) for each detection algorithm (from approach initiation to FAF) and the probability that final coasting positions at DH exceed LAAS integrity requirements, which are defined by the vertical alert limit (VAL) and lateral alert limit (LAL). The probability of mis-detection for an error occurring at epoch k, P md(err k), is computed as:
where P NA(k+n) is the probability of no-alarm at k+n epoch, δx vεr(DH) is the final-coasting vertical position error at DH and δx lat(DH) is the final-coasting lateral position error at DH.
In general, monitor detection performance is evaluated by comparing the probability of mis-detection with the allowable integrity risk allocated for the type of system failure targeted by the monitor. In the case at hand, cycle-slips fall under the category of LAAS H2 integrity risk (RTCA, 2004). The integrity risk allocated for this analysis for cycle-slip error is 75% of 1×10−9/10, which is derived from the LAAS category III signal-in-space integrity requirement, 1×10−9, by using an allocation method similar to category I in (RTCA, 2004) (75% for faulty conditions), assuming 10 types of satellite faults evenly divided. The VAL is 10 metres and LAL is 17 metres for LAAS category III.
Figure 12 shows that the built-in innovation detection algorithm outperforms the innovation integration detection algorithm on cycle slip detection performance. Note that the mis-detection probability for the innovation detection algorithm is below the requirement at the initial possible occurrence of 37 km, and then rapidly decreases to negligible values for shorter occurrence distances. In contrast, the mis-detection probability for the innovation integration detection algorithm to abrupt carrier measurement errors generally increases with time (decreasing distance to runway) because of the increasing number of chi-squared variables unaffected by the fault in the test statistic.
Simulated typical time responses of various error detection algorithms to half-cycle slip error on one satellite carrier-phase measurement.
6. innovation integration algorithm analysis
While the built-in innovation detection mechanism is clearly superior in the detection of abrupt changes, the innovation integration detection algorithm has the potential to reliably detect slowly growing errors (Groves, Reference Groves2008). This can be particularly important for applications such as LAAS, where atmospheric error can exhibit such behaviour.
Currently LAAS is a single frequency DGPS navigation system. Nominal ionospheric and tropospheric decorrelations have been accounted for in the aircraft positioning process. Abnormal ionospheric and tropospheric decorrelation errors have been observed and may cause hazardous misleading information (HMI) to LAAS users, which is dangerous to category I users and catastrophic to category III users (Lee and Luo, Reference Lee, Luo and Pullen2006) (Huang and van Graas, Reference Huang and van Graas2006). As currently envisioned, LAAS category I system integrity is ensured by implementing a number of monitors at the LAAS ground facility (LGF) (US DoT FAA, 2002); category III is still under development.
A copious amount of research has targeted these abnormal decorrelation errors to protect the integrity of LAAS category I users (Lee and Luo, Reference Lee, Luo and Pullen2006) (Huang and van Graas, Reference Huang and van Graas2006) (Gratton and Pervan, Reference Gratton and Pervan2006). However, the introduction of an inertial coasting phase to maintain robustness to RFI (i.e., the application investigated in this paper) requires additional monitoring and analysis to ensure the integrity of the INS calibration process through such anomalous atmospheric events. More specifically, some abnormal errors may cause potentially large coasting position errors due to erroneous INS calibration, while they might pose no threat to GPS-only users. The innovation integration provides additional detection capability which no other current detection algorithm can match in terms of detecting the threats that impact inertial coasting.
The simulated impacts on the hybrid navigation system from ionospheric storms will focus on the slow-moving storm fronts which are hardest to detect at the LGF (Luo and Pullen, Reference Luo, Pullen, Ene, Qiu, Walter and Enge2004). Based on the ionospheric threat space derived from the observations obtained from WAAS and the National Geodetic Survey Continuously Operation Reference Stations (NGS-CORS, or just CORS) to date (Luo and Pullen, Reference Heo, Pervan, Pullen, Gautier, Enge and Gebre-Egziabher2004), the threat space parameters and the corresponding equivalent static fronts seen by an approaching aircraft (i.e., the errors caused by the equivalent static fronts are the same as the slow moving fronts to an approaching aircraft) are displayed in Table 4.
Ionospheric storm threat space parameters.
Single-difference (SD) ranging errors due to weather storms must be considered as well, based on the observed tropospheric errors across short baselines (Huang and van Graas, Reference Huang and van Graas2006). For example, reference (Huang and van Graas, Reference Huang and van Graas2006) provides evidence for short-baseline transient differential errors in double-difference (DD) slant ranging of up to 40 cm during periods of severe troposphere activity.
A very large number of SD ranging error cases resulting from abnormal ionospheric and tropospheric decorrelations on both code and carrier-phase measurements were simulated by comprehensively sampling the ionospheric threat space and modelling tropospheric differential error transits as triangle and trapezoidal pulses with varying amplitudes and slopes. The pulse model is consistent with error transits actually observed in (Huang and van Graas, Reference Huang and van Graas2006). The total number of simulated ionospheric threats is 1876, with the assumption of only one satellite in view being affected by an ionosphere storm, and number of tropospheric threats is 18 055 with the assumption that up to half of the satellites in view can be affected by a troposphere storm. In this analysis the worst-case nominal (PDOP) satellite geometry is again used.
Figure 13 shows hard-to-detect examples of simulated ionospheric and tropospheric threats described above. Each line represents one threat which will cause the indicated SD zenith error as a function of the distance to the runway. The simulated ranging errors are then scaled by applying the appropriate obliquity factor. All the simulated threats are assumed not to affect the LGF measurements (i.e., they are not detectable by LGF). The location of LGF is conservatively chosen to be 5 km downwind from the runway end to avoid taking advantage of very short baseline distances (where nominal DGPS decorrelation errors are generally negligible, so that the system performance is improved). This is also the reason why some SD errors still exist when the aircraft reaches the runway end.
Examples of simulated ionospheric (left) and tropospheric (right) threats that are difficult to detect.
An ionospheric or tropospheric threat causes HMI for the hybrid navigation system if the evaluated detection performance satisfies the definition below:
where P event(e rrk) is the prior probability of the simulated error case happening and P md(e rrk) is the probability of mis-detection as defined in (42). Note that in the previous detection performance analysis for cycle slip events, P event(errk) was very conservatively assumed to be equal to one.
When a threat does not cause HMI given a specific detection algorithm, we can say the system integrity is protected from this threat by the algorithm. The total number of the protected threats for each evaluated detection algorithm is used as a metric to compare the relative detection performance of the algorithms. Because the appropriate values for prior probabilities for ionospheric and tropospheric storms are still an open question, prior probability is parameterized in this analysis. Three prior probabilities values used for ionospheric threats are 10−3, 10−4, 10−5 and 10−2, 10−3, 10−4 for tropospheric threats.
Five detection algorithms were evaluated against all the simulated ionospheric and tropospheric threats: 1) innovation detection, 2) innovation integration detection, 3) smoothed code RAIM detection, 4) carrier-phase relative RAIM detection; 5) carrier-phase relative RAIM detection with INS augmentation.
Figure 14 (left) shows the number of the protected ionospheric threats (y-axis) which can be detected exclusively by one detection algorithm (algorithm numeric identifier on x-axis). The second detection algorithm, the innovation integration detection algorithm, clearly demonstrates superior detection capability.
(Left) Relative detection performance to detect simulated ionospheric threats. (Right) Relative detection performance to detect simulated Type I tropospheric threats.
Figures 14 (right) and 15 confirm the capability of the innovation integration detection algorithm to capture tropospheric threats that the rest cannot, even though multiple satellites can be affected by troposphere storms. In these figures, three bars are shown for each detection algorithm, each bar corresponding to a different number of affected satellites. These results are presented in three groups according to the different types of simulated SD tropospheric errors. They are symmetric triangle shape (Type I), asymmetric triangle shape (Type II) and trapezoidal shape (Type III, including symmetric and asymmetric). As seen in the figures, the innovation integration detection algorithm has excellent detection performance for all three groups of simulated tropospheric threats.
Relative detection performance to detect simulated Type II (left) and Type III (right) tropospheric threats.
7. CONCLUSION
A tightly-coupled GPS/INS integrated navigation architecture was described, implemented and tested. The unified system, based on the direct fusion of GPS and INS through Kalman filter state dynamics, has rigorously accounted for the spatial decorrelations in DGPS measurements. By fusing all GPS and INS estimates through state dynamics in one centralized extended Kalman filter, the integrity of GPS positioning and INS calibration is ensured under fault-free error conditions.
Validation of the integrated navigation system was done by post-processing flight data from the FAA Technical Center to verify the accuracy of the covariance analyses performed.
Finally, a GPS/INS fault detection algorithm, based on the integral of Kalman filter innovations, was analyzed and evaluated against other existing GPS fault detection algorithms. It was shown that the innovation integration detection algorithm outperforms other RAIM-based methods in detection of slowly developing ranging errors, such as those caused by ionospheric and tropospheric anomalies.
ACKNOWLEDGMENTS
The authors would like to thank Professor Demoz Gebre-Egziabher for his helpful advice on the modeling of inertial navigation systems. We would also like to express our appreciation to the Federal Aviation Administration (FAA) for supporting this research and the FAA William J. Hughes Technical Center for providing the flight data.
