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:

(1)

where δⁱΔρ_k is the single difference L1 code measurement error between user and reference station for satellite i at time k. δΔx_k is the error in the relative (user minus reference) position vector, and ⁱe_k 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, δⁱΔT _k, δⁱΔI _k,δⁱΔm _ρ,k and ⁱΔε_ρ,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 ⁱΔφ _kλ (metre) as well, where λ is the wavelength of L1 carrier. The magnitudes of carrier multipath ⁱΔm _φ,k and receiver carrier-tracking noise ⁱΔε_φ,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 ⁱΔ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:

(2)

where δΔn is the residual refractivity uncertainty, h ₀ 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 ⁱΔT _{dc, k} and refractivity uncertaintyδΔn.

Ionospheric decorrelation model:

(3)

where ⁱΔ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 _dc,k and the vertical ionospheric gradient error δⁱΔV _ig.

Combining (1)–(3), a measurement error model describing code and carrier measurement errors for one satellite can be formed:

(4)

If there are n satellites in view, the complete set of estimated DGPS state errors is:

(5)

The complete measurement error model corresponding to the errors in the estimated DGPS state S_G in (5) can be expressed as:

(6)

where δZ_k, H_k and ε_G,k are defined as:

(7)

and the measurement error vectors δZ_Δρ,k,δZ_Δφ,k and the corresponding observation matrices H_Δρ,k, H_Δφ,k in the (7) above are:

(8)

(9)

(10)

where 0_U(n-1) and 0_L(n-1) are the upper and lower triangular part of a size n-1 square zero matrix, and I_n is an n×n identity matrix (to save some space, a square matrix is noted by its dimension only once. i.e. using I_n instead of I_{n ×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:

(11)

The dynamic matrix F_G and the process noise vector W_G are:

(12)

(13)

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. v_e 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 . f_clk 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:

(14)

Throughout this error analysis, the state estimation errors at time k, δS_G,k, are related to the true states, S_G,k, by:

(15)

For Kalman filter implementation, the error propagation in time and the covariance of state estimation error can be formed as:

(16)

where , and are defined as:

(17)

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:

(18)

where K_k is the Kalman filter gain and is the covariance matrix of the errors on available measurements. They are defined as:

(19)

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):

(20)

where is the time derivative of the ground velocity (v_e) 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⁽ⁿ⁾ is the specific force expressed in navigation frame n coordinates. Similarly, w_ie⁽ⁿ⁾, w_en⁽ⁿ⁾ and g₁⁽ⁿ⁾ are all expressed in navigation frame n coordinates, where w_ie⁽ⁿ⁾ is the earth's rotation rate; w_en⁽ⁿ⁾ is the angular velocity of the navigation frame relative to the earth (also known as transport rate); g₁⁽ⁿ⁾ is the local gravity vector.

Euler angle propagation is adopted as the basic mechanism for vehicle attitude propagation, which is sufficient for analysis purposes:

(21)

where Ψ is the Euler angle vector consisting of roll (Φ), pitch (θ), and yaw (Ψ) angles. F_EU is the transformation matrix to translate the instantaneous body-to-navigation-frame angular velocity,w_nb^(b), into the Euler angle rotation rate (Titterton and Weston, Reference Titterton and Weston2004) (Jekeli, Reference Jekeli2001). w_ib^(b) is the gyro measured angular velocity of the body relative to an inertial frame and C_n ^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:

(22)

(23)

where w_in⁽ⁿ⁾ is the navigation frame angular rate relative to an inertial frame. It can be decomposed as the sum of w_ie⁽ⁿ⁾ and w_en ⁽ⁿ⁾ 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 w_en ⁽ⁿ⁾ can be expressed in terms of user ground velocity error as:

(24)

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:

(25)

Inertial sensor errors generally propagate into navigation system errors through δf^(b) in (22) and δw_ib^(b) in (23), which include IMU measurement bias error δb_a/g, misalignment error δm_is,a/g, scale factor error δs_f,_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:

(26)

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:

(27)

where τ_g and τ_a are the time constants of GMRP models for gyroscope and accelerometer biases respectively. n_g and n_a 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:

(28)

where δh is the vertical position error, δξ and δη are the DOV errors in north-south and east-west directions respectively. g ₀ and R ₀ are the nominal values of local gravity and the earth's radius. δΔg ₀ is the error in the local gravity magnitude. Simplified notations for the vectors and matrices in (28) are as adopted as follows:

Equation (28) can now be re-written as:

(29)

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:

(30)

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. n_gm is the vector of driving white noise.

Collecting (22), (25), (27) and (30), the dynamics of INS errors can be formed as:

(31)

where

Using simplified notation, equation (31) can be expressed as:

(32)

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:

(33)

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 :

(34)

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 v_e from W_G and F_h2Vδh from W_I to the state dynamics as specified by (32) and (29) respectively. The resulting centralized EKF state error dynamics are described by:

(35)

where F_I2G, F_G2IW_I1and W_G1 are defined as:

(36)

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):

(37)

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:

(38)

where δZ_k, H_k 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.

Figure 1. The nominal velocity (left) and attitude (right) calibration performance of the integrated navigation system with error model parameters in Table 1.

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.)

Table 2. 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.

Figure 2. 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.

Figure 3. 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.

Figure 4. Deviations between the reference attitudes and the pure INS-only propagated attitudes for three Euler angles: roll δφ, pitch δθ and yaw δψ.

Table 3. 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.

Figure 5. 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.

Figure 6. 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.

Figure 7. Position errors of GPS/INS integrated system during approach 2.

Figure 8. 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.

Figure 9. 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:

(39)

where Z_k, H_sys,k, and are defined according to (38):

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)~χ²(n _k)). The detection threshold D _{th,inn o}(k) is conservatively set to provide a false alarm probability of 2×10⁻⁶/30, which is based on the most stringent LAAS category III continuity risk requirement (2×10⁻⁶ 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:

(40)

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.

Figure 10. 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:

(41)

The threshold is plotted in Figure 11 along with a Monte Carlo simulation of 1000 approaches under normal error conditions.

Figure 11. 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:

(42)

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 H₂ integrity risk (RTCA, 2004). The integrity risk allocated for this analysis for cycle-slip error is 75% of 1×10⁻⁹/10, which is derived from the LAAS category III signal-in-space integrity requirement, 1×10⁻⁹, 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.

Figure 12. 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.

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.

Figure 13. 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:

(43)

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(e_rrk) 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⁻³, 10⁻⁴, 10⁻⁵ and 10⁻², 10⁻³, 10⁻⁴ 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.

Figure 14. (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.

Figure 15. 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.