WO1999018677A1 - Gps augmentation system - Google Patents

Abstract

A system providing information that augments the navigational data transmitted by GPS satellites by providing fast and long-term corrections, a user differential range error, a grid ionospheric vertical error, and ionospheric corrections. The augmentation system includes one or more reference stations (2) having two or more reference receivers that independently receive signals from the GPS satellites and transmit the received information to at least one master station (8). Each master station has a set of processors, designated as correction and verification processors, each configured to receive the output signal from one of the reference receivers and calculate a set of augmentation data from that output signal. The augmentation system compares the augmentation data produced by the correction processor to the augmentation data produced by the verification processor to validate operation of augmentation system and transmits the validated augmentation data to users of the system.

Description

GPS AUGMENTATION SYSTEM

CROSS-REFERENCE TO RELATED APPLICATIONS

This application hereby claims the benefit under 35 U.S.C. § 119(e)(1) based on provisional U.S. Application Serial No. 60/028,731, filed on October 4, 1996, now abandoned.

BACKGROUND OF THE INVENTION Field of the Invention.

The present invention relates to a method of and apparatus for augmenting the Global

Positioning System.

Description of the Related Art

The NAVSTAR Global Positioning System (GPS) is a continuous, space-based navigation system that provides any suitably equipped user with highly accurate three-dimensional position, velocity and time information anywhere on or near the earth. GPS is basically a ranging system. To provide accurate ranging measurements, which are time-of-arrival measurements, very accurate timing is required in the satellites. To provide the required timing accuracy, the GPS satellites contain atomic clocks. The GPS satellites are in approximately 12 hour orbits (11 hours, 57 minutes and 57.27 seconds) at an altitude of approximately 11,000 nautical miles. The total number of satellites in the provide coverage at all locations on the earth as nearly to 100% of the time as possible. Each satellite transmits signals at two frequencies in L-Band, i.e., 1575.42 MHz (LI) and 1227.6 MHz (L2). The presence of the second frequency permits ionospheric refraction corrections by properly equipped users. These signals are modulated with synchronized, satellite-unique, pseudorandom noise (PRN) codes that provide instantaneous ranging capability. Those codes are modulated with satellite position, clock and other information, in order to provide users with that information which is also required for ranging. In the GPS navigation solution, three measurement capabilities (ranging, Doppler and integrated Doppler) exist. It is not uncommon for a specific user's equipment to use at least two of the three measurement capabilities to simultaneously solve for position, velocity, receiver clock offset and clock drift, in some cases using all satellites-in-view, which can number as high as twelve.

GPS provides two positioning services, the Precise Positioning Service (PPS) and the Standard Positioning Service (SPS). The PPS can be denied to unauthorized users, but the SPS is available free of charge to any user worldwide. Users that are crypto capable are authorized to use crypto keys to always have access to the PPS. These users are normally military users, including NATO and other friendly countries. These keys allow the authorized user to acquire and track the encrypted precise (P) code on both LI and L2 frequencies, and to correct for intentional degradation of the signal. Encryption of the precise code provides GPS with an Anti-Spoofing (A-S) capability. A- S is not meant to deny the P code to unauthorized users, but to prevent the spoofing of the precise code by an unfriendly force. But A-S does deny the P code to unauthorized users. Thus, A-S inhibits these users from correcting for ionospheric refraction, since the L2 signal only carries the P code, although there are "codeless cross-correlation" techniques that do allow a degraded form of this measurement. A- S does not prevent the use of the Coarse/ Acquisition (C/A) code, which is only carried on the LI signal. Intentional degradation, also known as Selective Availability (SA), is also included in the GPS signal to deny accuracy to an unfriendly force. Unfortunately, SA also denies accuracy to unauthorized users that are friendly, which is the entire civil community. The peace-time policy of the DoD is to provide an SPS accuracy of 100 meters, 2-dimensional horizontal accuracy 95% of the time. Either A-S or SA or both may be turned-on.

As shown in Fig. 1, GPS is comprised of three segments, i.e., space, user, and control segments. The monitoring and satellite control sites are dispersed around the world. The functions of these three segments are summarized in Table 1.

Table 1. GPS Segment Functions

Because of SA and A-S, the GPS system is not as accurate for some navigation applications (such as precision approach) as is desirable, and is only accurate to about 100 m. Accordingly, it would be desirable to augment the navigation data from the GPS system to provide improved accuracy, such as 10 m or less. In order to overcome the accuracy deficiency in GPS, as well as deficiencies in availability and integrity, the Federal Aviation Administration proposed a Wide Area Augmentation System (WAAS) to augment the GPS Standard Positioning System. Conceptually, GPS augmented with the WAAS would serve as a primary navigational system and provide information to all properly- configured aircraft in a particular region of interest (referred to as the service volume) to support all phases of flight, including en route through precision approach navigation and landing. In general, the WAAS was foreseen as providing two services: correcting and providing integrity information corresponding to the data received from GPS and Geostationary Earth Orbit (GEO) satellites (discussed below); and providing a ranging capability. The WAAS was foreseen as improving the integrity and availability of the GPS information, replacing and supplementing the current Category I Instrument Landing Systems (ILS), and providing other significant advancements. WAAS will not provide Category II/III precision approach capability.

With the WAAS, every runway that is suitable for instrument approaches would become a candidate to economically implement a Category I precision approach capability. Thus, aircraft within the service volume would have access to many more airports and runways, resulting in improved schedule reliability, reduced flight cancellations, and fewer diversions. Also, in high density terminal areas, there would be additional runways that could be used in instrument conditions and more secondary airports available to absorb the capacity demands.

As foreseen by the FAA, a GPS-based navigational system, augmented with WAAS, would also enhance safety by reducing cockpit workload and minimizing the possibility of mid-air collision and controlled flight into terrain. A GPS/WAAS system would also provide accurate vertical guidance and could be exploited by users who now have only the ability to perform nonprecision approaches. Current standard flight procedures could be redesigned to provide vertical guidance along a predetermined descent path to the minimum decision altitude based upon vertical guidance provided by GPS/WAAS.

WAAS would also provide the opportunity to optimize flight routing because the routings would no longer be dependent on the placement of ground-based navigational aides. The present airway system could be restructured to provide users with shorter routes and improved use of altitude and winds. Additionally, by exploiting the inherent flexibility in routing, alternate/parallel routes could be used to meet changing traffic situations and to improve recovery time after flow control restrictions, such as those caused by severe weather conditions. By increasing system capacity in high density areas, system delays would be reduced. Furthermore, improved navigational accuracy provided by the WAAS would offer the opportunity to incrementally reduce separation standards. Potential reductions include non-radar separations in en route airspace and terminal separations due to smaller obstacle clearance areas and protected airspace. These reduced separation standards would directly translate into increased system capacity and reduced delays with associated user cost savings. To date, no implementation of the WAAS concept has been developed.

Accordingly, it is an object of the present invention to provide an implementation of a Wide Area Augmentation System that achieves the objectives of the WAAS system proposed by the FAA.

SUMMARY OF THE INVENTION The current invention provides for a GPS augmentation system that is capable of providing a set of augmentation data, including long-term and fast corrections and ionospheric corrections, for a group of GPS satellites, wherein each GPS satellite broadcasts GPS signals that contain navigational data on both the LI and L2 frequency bands. The navigational data is preferably formulated to determine the range to a particular GPS satellite. Preferably, the augmentation system includes one or more reference stations, wherein each such reference station has two or more reference receivers for independently receiving the GPS signals from the group of GPS satellites. The reference stations preferably are capable of generating an output signal comprising the received navigational data. The presently preferred augmentation system further includes at least one master station that is in communication with the reference stations and receives the output signals emitted by the reference stations. Preferably, each master station has a set of processors, designated as correction and verification processors, wherein the correction processor is configured to receive the output signal from one of the reference receivers and calculates a set of augmentation data from that output signal, and the verification processor is configured to receive the output signal from a different one of the reference receivers and calculates the augmentation data using this output signal. The preferred augmentation system compares the augmentation data from the correction processor to the augmentation data from the verification processor to validate operation of augmentation system and transmits the validated augmentation data to users of the system.

In a preferred embodiment of the GPS augmentation system, each reference station has at least two primary reference receivers and at least one standby reference receiver, and wherein the master stations have the ability to detect a failure in one of the primary reference receivers and switch to receive the output signal data from the standby reference receiver in the event of such detected failure. In this embodiment, the correction and verification processors each compute standby clock and ionospheric biases for use when switching to receive the navigational data from the standby reference receiver. In another preferred embodiment, the GPS augmentation system is capable of detecting multipath error originating at a particular failed satellite. In this embodiment, the reference receivers comprise a narrow and a standard width correlator receiver, each capable of independently receiving the navigational data from the satellites, extracting multipath error in the navigational data, and computing and storing a pseudorange to each of the satellites. The system then computes the average value of the multipath error extracted from the navigational data prior to computing the pseudoranges over a period of time and compares the average values of the multipath error extracted to detect multipath error originating at a particular satellite.

Additionally, the present invention provides for a means to calculate a user differential range error that bounds the error of the long term corrections and fast corrections with a probability of 99.9% and a grid ionospheric vertical error that bounds the error of the transmitted ionospheric corrections with a probability of 99.9%.

In yet another embodiment, the GPS augmentation system includes at least one ground uplink station capable of receiving the augmentation data from the master stations and generating and transmitting uplinked signals that include the augmentation data to at least one GEO satellite, which then re-transmits the augmentation data. Preferably, the uplink station has signal control means for generating signals so that the GEO satellite transmissions mimic a conventional GPS signal.

BRIEF DESCRIPTION OF THE DRAWINGS To facilitate further discussion of the invention, the following drawings are provided in which:

Fig. 1 is a block diagram showing the major components of a prior art Global Positioning System;

Fig. 2 is a block diagram showing the major components of the GPS augmentation system in accordance with the invention;

Fig. 3 illustrates the preferred locations of the major hardware components of an augmented GPS navigational system for use in the continental United States;

Fig. 4 is a block diagram illustrating the major components of a reference station ("WRS") in accordance with the present invention;

Fig. 5 is a block diagram illustrating the major functional components of a reference station equipment ("WRE") in accordance with the present invention; Fig. 6 is a block diagram showing the carrier/ code tracking relationship in

WRE and GEO Uplink Subsystem ("GUS") receivers in accordance with the present invention;

Fig. 7 is a block diagram illustrating the major components of a master station ("WMS") in accordance with the present invention; Fig. 8 is a block diagram showing the major functional elements of the invention; Fig. 9 is a block diagram showing the processing and interface elements of the master stations in accordance with the present invention;

Fig. 10 is a block diagram showing the carrier phase pseudorange ambiguity resolution process in accordance with the present invention; Fig. 11 is a block diagram showing the carrier phase pseudorange ambiguity using a correction for ionospheric divergence process in accordance with the present invention;

Fig. 12 is diagram illustrating the process of calculating the ionospheric corrections in accordance with the present invention; Fig. 13 is a process flow diagram illustrating a Kalman filter as used in the invention;

Fig. 14 is an illustration showing an ionospheric grid for the continental United States;

Fig. 15 is an illustration showing the pierce point interpolation process in accordance with the invention;

Fig. 16 is a block diagram of the WAAS User Correction Process of the invention;

Fig. 17 is a block diagram of the WAAS Fast Correction Generation Process of the invention; Fig. 18 is a block diagram of the UDRE Computation Process of the invention; Fig. 19 is a block diagram of the Ionospheric GIVE Determination Process of the invention;

Fig. 20 is a process diagram for the processes of the Safety Monitor of the invention; Fig. 21 is a block diagram of the GEO and WNT clock steering process of the invention;

Fig. 22 is a block diagram showing the major components of the GEO Uplink Subsystem ("GUS") of the invention;

Fig. 23 is a block diagram showing the major components of the GUS Signal Generator Subsystem of the invention;

Fig. 24 is a block diagram showing the major components of the GUS RF Uplink of the invention.

Fig. 25 is a Kalman filter and control laws for the GUS signal control process of the invention; and Fig. 26 a block diagram of a Kalman filter for the GUS signal control process of the invention.

DETAILED DESCRIPTION OF THE INVENTION The present invention is a Wide Area Augmentation System ("WAAS") that augments the Department of Defense Global Positioning System (GPS) Standard Positioning Service. The system provides a signal to system users to support en route through Category I precision approach navigation and landing. Users of the system would include, for example, certified aircraft using the WAAS for any approved phase of flight, or any other appropriate users requiring improved navigational capabilities. The system provides two services: (1) correction and integrity data for GPS and Geostationary Earth Orbit (GEO) satellites (GEO satellites are not part of the GPS system and will be discussed in detail below) and (2) additional GEO ranging capability.

As shown in Fig. 2, the navigation data from GPS satellites 4 and GEO satellites 6 is received and processed at geographically dispersed WAAS Reference Stations 2 ( "WRSs" or "reference stations"). This data is forwarded to WAAS Master Stations 8 ("WMSs" or "master stations") which process the data to determine the integrity, differential corrections, and residual errors for each monitored satellite and for predetermined ionospheric grid points. This information is sent to GEO Uplink Subsystems (GUSs) 10 and uplinked along with GEO navigation messages to GEO satellites 6. GEO satellites 6 downlink this data to the users on the GPS LI frequency with a modulation similar to that used by GPS. Each of the above system components is connected via a Terrestrial Communications System 12.

As shown in Fig. 3, initially, the system preferably consists of 24 WRSs, 2 WMSs, and 4 GUSs, with each WRS communicating all of its data output to each WMS, and with each WMS communicating all of its data output to each GUS. Thus, if any system component were to become inoperative, system operation may be maintained due to the built-in system redundancy. The continental United States is covered by 3 GEO satellites, i.e., POR, which covers the western portion of the U.S., AOR-W, which covers the middle portion of the U.S., and AOR-E, which covers the eastern portion of the U.S. Each of the west coast GUSs is capable of communicating with either the POR or AOR-W satellites. Similarly, each of the east coast GUSs is capable of communicating with either the AOR-E or AOR-W satellites. Thus, for each GUS, there is at least one redundant GUS to communicate with each GEO satellite. As discussed above, reference stations 2 collect and transmit the GPS and GEO satellite information to master stations 8, which then integrate the data from the various WRSs, determine the integrity of the satellite ranging information, and calculate satellite and ionospheric correction data and GEO ephemeris and almanac data. The integrity and correction data generated by master stations 8 is then transmitted to each GEO satellite 6 using GUSs 10. In turn, this information is then broadcast by GEO satellites 6 using GPS-like signals and received by any participating aircraft 14 or other user.

Each master station 8 collects the data from reference stations 2 and performs various functions including: (1) determining ionospheric corrections; (2) satellite orbit determination; (3) determining satellite corrections; (4) determining satellite integrity; and (5) providing independent data verification and validation. Each master station 8 performs these calculations on independent processors within the master station using independent sets of data received from the reference stations, thereby allowing the master station to verify the results of its computations.

The GPS correction data calculated by master stations 8 is formatted into WAAS messages and sent to GUSs 10 to be uplinked, along with the GEO navigational messages, to GEO satellites 6. GEO satellites 6 in turn downlink the data using GPS-like signals on the GPS LI frequency for reception by participating users, thereby providing the users with the correction and integrity information. Since the GEO signal also appears as a GPS ranging signal, the users are provided with more ranging sources that GPS alone, thus improving continuity of function and availability.

Having described the general functionality of the invention, the reference stations, the master stations and the GUSs will now be described in further detail. WAAS REFERENCE STATIONS (WRSs

As discussed above, each WRS 2 performs the functions of data collection, reasonability checking, data processing, data recording, and data transferring. Referring to Fig. 4, each WRS consists of three WAAS Reference Equipment (WRE) units 16 which collect independent sets of data including GPS satellite data, GEO satellite data, and local tropospheric data, and transmits the data to each WMS 8. Independence of the data is ensured by gathering the data through independent hardware and is necessary to support the verification function performed by WMSs 8. The data is collected at a rate consistent with its expected level of variation; e.g., slowly changing weather condition data is collected less frequently than data from the GPS satellites. Prior to transmitting data to WMSs 8, each WRE 16 verifies the reasonability of its collected data. Failed data is marked as having failed a reasonability test and is forwarded to the WMSs.

Each WRE 16 in a WRS 2 is functionally identical and generates identical data threads, thereby providing sufficient redundancy to ensure the availability of at least two WREs 16 providing data to master stations 8 in the event that a WRE fails in a reference station, so that master stations 8 will continue to receive primary and secondary data sources, to allow for independent verification of the data. Reference stations 2 also include a computer control unit 18, having a keyboard 22 and monitor 20, and a remote status unit 24 capable of providing onsite audio and visual indicators for maintenance personnel.

To increase the availability of the data at each WMS 8, each WRS 2 transmits data through two independent nodes of the Terrestrial Communication System. Toward this end, each WRS 2 includes two communication routers 26, each of which receives, multiplexes, and transmits to each master station, all of the data from each WRE 16 in the WRS. As shown in Fig. 5, each WRE 16 consists of a WAAS receiver subsystem 28, an atomic clock 30, a GPS antenna subsystem 32, a WRE processor 34 interfaced to routers 26, and temperature/humidity and barometric pressure sensors 36 and 38. Each WAAS receiver subsystem 28 is connected to a standard GPS antenna (not shown) that includes a receiving antenna with associated collocated filters, amplifiers, and a voltage regulator. A preferred GPS antenna assembly is model 448078-0100 available from Wilcox Electric Incorporated. WAAS receiver subsystem 28 is preferably a 16-channel receiver that collects data from the GEO and GPS , including GPS Ll-C/A pseudorange data, GPS L1/L2 code differential data, GPS satellite navigational data, GEO signal-in-space, and maintenance and status data. In a preferred embodiment, each WRE processor 34 is preferably an IBM 7248-100 (43P), WAAS receiver subsystem 28 is preferably a NovAtel WAAS Receiver 044645-0001, and router 26 is preferably a Bay Networks ANH2. It is foreseen that other equipment may be substituted for the particular equipment described herein, provided that such substitute equipment performs the functions described herein. Atomic clock 30 provides a time source used to accurately compare the signals received from the GEO satellites with GPS time and, preferably, meets a stability requirement of at least one in 2 * 10¹³ parts over 24 hours. As discussed below, several atomic clocks are integrated together and employed to enable accurate satellite orbit determination. A preferred atomic clock is the model FTS 4040 A/RS from Frequency & Time Systems, Inc. Weather sensors 36 and 38 measure the air temperature, barometric pressure, and the relative humidity in the vicinity of the WRS. This data is forwarded to master stations 8 to be used in determining the effect of troposphere delay on the GPS and GEO satellite transmissions and to provide corrections for these effects.

Each WRE processor 34 receives data from the GPS antenna and performs the following functions: i) GPS/GEO data are collected and checked for reasonability, LI data at 1 Hz, and L1/L2 at 0.1 Hz; ii) GPS/GEO data are smoothed and averaged, LI data at 1 Hz, and L1 L2 data at 0.1 Hz; iii) GPS/GEO LI pseudorange differential corrections collected at 1 Hz; iv) pressure, humidity, and temperature readings are collected and checked for reasonability at 0.1 Hz; and wet and dry vertical tropospheric delays are calculated at 0.1 Hz.

Generally speaking, WAAS receiver 28 receives the GPS and GEO ranging signals and processes these signal as conventional GPS signals in the manner well-known in the art. However, as discussed below, WAAS receiver 28 performs several unique functions that lead to an overall improvement in the performance of the system. WAAS Receiver Functions

As shown in Fig. 6, the tracking of GPS signals in conventional GPS receivers consists of tracking the signal carrier (200) and signal PRN code (202) (envelope) simultaneously. Both the code and carrier have Doppler and Doppler rate components primarily due to motion of the GPS satellites. Code tracking in conventional receivers is aided by the measured carrier Doppler, which is inherently more accurate than code Doppler, scaled by 1/N to the code frequency (204). For example, for the LI C/A code, N is 1540. In this way, the bandwidth of the code tracking loop can be significantly narrowed, and thus provide much more accurate and more smooth code phase (pseudorange) measurements. In a typical GPS receiver, measurements are taken at a rate of l/20Hz or slower. These (smoothed) code phase measurements are correlated in time and closely time-spaced and therefore do not provide an independent measure of code status. The availability of independent measurements is critical for the WRS in order to timely sense pseudorange anomalies that may be different than those observed from the carrier phase. In the present invention, code phase measurement smoothing is performed later outside of the receiver at the WMS so that independent measurements are available at a 1 Hz rate. This is accomplished by controlling the code tracking bandwidth in the WRS Receiver via a command (as shown in Figure 6) to a value desired for obtaining independent measurements. Another significant feature preferably included within the receiver 28 is used to detect multipath error, a serious problem associated with receiving signals from GPS satellites. Multipath error arises when a direct signal from a satellite that impinges on a GPS receiver's antenna is distorted by reflections of the same signal. Such reflections can originate from many sources, including the ground or any objects near the receiver, such as a building. Any multipath distortion in the signal received at a reference station receiver from a GPS satellite may introduce errors in the pseudorange measurement between the receiver's antenna and the satellite and may decrease the accuracy of the corrections broadcast to the participating aircraft. Distortion from multipath error may be a significant problem for WAAS because the receiver antennas within the reference stations may be located on the rooftops of FAA air traffic control facilities in an environment containing many reflectors, e.g., the roof itself, additional antennas, and air conditioner units, etc.

It is well-known that position errors due to multipath distortion may be evaded using a number of techniques. For example, the proposed invention utilizes a receiver consisting of an array of narrow correlators configured so as to provide a multipath estimation delay lock loop (MEDLL). This combination of narrow correlators and MEDLL provides a substantial reduction in the errors prOoduced by multipath effects, particularly as compared to the use of standard width correlator receivers. However, if certain types of satellite failure occur, multipath error may also originate at the satellite itself. If receiver 28 and the user's receiver were to have the same width correlators, multipath errors originating at the satellite would be eliminated because the reference receiver and the receiver located on the aircraft would each extract the same amount of error. However, in the preferred WAAS configuration, the user's receiver may be configured using any width correlator, thereby complicating the ability to distinguish multipath distortion originating with the satellite, e.g., due to a malfunctioning satellite, from distortion due to reflections near the receiver antenna. In order to overcome this problem, receiver 28 in reference stations 2 has both narrow and standard width correlators. Because both correlators receive signals from the same antenna, multipath error originating near the antenna will cause an incorrect pseudorange to be calculated by the two correlators. At any point in time, the narrow correlator will remove more multipath error than the standard-width correlator. Over time, e.g., a few hours, as the GPS satellites move in relation to the receiver, the nearby multipath effects will average out to approximately zero. However, the multipath effects originating from the satellite will remain constant and, thus, if the average value of multipath error removed by the two correlators is not zero, then the multipath error must be originating with the satellite. Thus, by using a receiver with both narrow and standard width correlators, the problem of multipath originating at the satellite may be detected despite the presence of multipath from reflectors located near the receiver antenna. By storing and comparing the range results from these two correlators over time, the system can detect multipath distortion originating from a failed satellite and then notify the participating aircraft not to use or rely upon the data originating with the failed satellite. WAAS MASTER STATIONS (WMSs Each WMS 8 performs the functions of correction processing, satellite orbit determination, integrity determination, verification, validation, and WAAS message generation. Once per second, each WMS 8 collects the data received from all WRSs 2 and processes this data to perform these functions. This processing is performed on all available WRS data and results in the transmission to all GUSs of a formatted 250-bit WAAS message once per second. The timing of WMS processing is scheduled to allow broadcast of the resulting WAAS message from the GEO satellite coincident with the desired 1 -second Coarse/Acquisition (C/A) code epoch. The WAAS validates the Signal-in-Space by checking the downlinked messages to ensure that they are identical to those transmitted and by comparing navigation position solutions from WAAS/GPS with the surveyed WRS locations. Each WMS includes an Operation and Maintenance console from which control over the WAAS can be exercised via a computer-human interface at the WMS. To avoid conflicts, only one WMS within WAAS can be designated as the controlling WMS (active Operation and Maintenance console) at any one time. As shown in Fig. 7, each master station 8 includes three similar processors: a corrections processor 40, a verification processor 42, and an operation and maintenance processor 44, the latter including a display, a keyboard, a printer, and a modem. All three processors 40, 42, and 44 are connected to the terrestrial communication system using a Backbone Concentrator Node (BCN) router and an Ethernet hub. An additional processor 50 is provided for a safety monitor. Each workstation 40, 42, and 44 is preferably based on and IBM RISC System/6000 Model 39 or other appropriate multi-user, multitasking, numeric-intensive processor.

The correction and verification workstations 40 and 42 each contain identical software (the correction and verification program), but operate independently on independent data. As shown in Fig. 8, and described in detail below, the correction and verification program determines ionospheric corrections, satellites orbits, satellite corrections, satellite integrity and provides a record and playback feature for archiving data received and processed by the WMS. Additionally, a safety monitor, also discussed below, resides on a separate processor in the master station. Each master station 8 also includes a GPS antenna 46 and a GPS clock 48 connected to correction and verification workstations 40 and 42. WMS Error Correction

Each WMS 8 determines correction or augmentation data to be broadcast to the user via uplinks at the GUSs. These corrections consist of Ionospheric Grid Point (IGP) Delays and Grid Ionospheric Vertical Error (GIVE) estimates, and long term corrections, fast corrections and associated User Differential Range Error (UDRE) for each GPS satellite. Also, Clock Steering commands are generated to keep WAAS network time within 50 nanoseconds of GPS time and to ensure that the GEO clock is within the range of the type 9 GEO navigation message. Fig. 9 is a block diagram showing operation of the WMSs. Each element of Fig. 9 will now be discussed in further detail. Process Input Data GOO)

To ensure the integrity of the correction generation process, two sets of independent input data from each reference station 2 (from 2 of the 3 WREs at the reference station) are processed by master station 8, one to determine correction data, and one to determine verification data to validate the accuracy of the correction data. Also, since it is necessary to switch immediately to a standby thread of data in case of a failure of a main thread, i.e., a source of the correction or verification data, a third source of backup data is provided at each WRS and the receiver clock biases and ionospheric L1/L2 bias are computed in real-time.

A Process Input Data Task 300 in both the correction processor and the verification processor selects and monitors all three input data threads. These two tasks communicate with each other to ensure that independent data from all reference stations is maintained between the two processors. Also, both processes compute the WRE clock and ionosphere biases to ensure that a switchover can be made to the standby thread immediately in case of a failure in either the current correction thread or verification thread of data.

The standby WRE clock bias and frequency bias are computed using the LI C/A code pseudorange and carrier phase measurements from both the standby thread and the thread being used in the correction process. Real-time receiver clock biases of the online units are computed by the Real- Time Orbit Determination Model. The LI clock bias (Sa_{fo i>jmLi j C}) for the standby WRE thread is

computed as

where:

and where: ^ε _L\j.c i^s h bias residual computed from the received pseudorange measurements from

the primary and standby threads each second with a C/No > than the bias cutoff limit. PRLIJJ i^s the measured LI C/A code pseudorange.

&^l _fo,bi_asLi _,Λ ^ιs th estimated receiver clock bias for the primary thread

A denotes a primary thread

C denotes the standby thread data.

N number of measurements and estimates made from both the Primary and Standby

WRE threads in an estimation interval of 30 minutes. This computation is preferably done as a running 30 minute average. The variance is computed as:

A similar process is performed using the carrier phase on LI to determine the standby WRE frequency bias. The WRE L2-L1 bias (x_L2._{L biasJ} )«<»*, for the standby WRE thread is computed as

(^XL2-Llbύs,i )^B°*S = -(! /

tt² - ^Λl, - X_L2._UJUt,j ~ u-LlJΛuj)* ~(Σ_M ^T L2-L ) W ³) where:

(T_L2__LUJ)A measured slant delay difference for the primary thread,

(^X _L2-_L\ i_asj ^A computed L2-L1 bias for the primary thread,

(^X _L2-_{L\ b}i_as )^λ computed L2-L1 bias for the satellite.

^τ _L2-_{L j})^B measured slant delay difference for the standby WRE thread,

N number of measurements and estimates made from both the Primary and Standby WRE threads.

If desired, other types of standby bias filters could be implemented to replace the averaging process. This would eliminate the necessity of keeping older measurements in the data. For example, the filter could be a low pass filter or a Kalman filter .

Preprocessing/Multipath Removal (302)

The raw pseudorange measurements from the GPS satellites must be corrected for several different types of errors before they can be used in the orbit determination and correction determination processes. These errors include relativity, troposphere, ionosphere, and multipath. Data Preprocessing/Multipath process 302 is included in correction and verification processors 40 and 42 to determine corrections for these errors.

Troposphere correction

A preferred method of determining the troposphere correction uses the Berman/Chao tropospheric model, which is discussed in detail in i). Chao, C.C. "A New Method to Predict Wet Zenith Range

Correction From Surface Measurements." JPL Technical Report 32-1526; ii) Chao, C.C. "The

Tropospheric Calibration Model for Mariner Mars 1971." JPL Technical Report 32-1587; and iii) Studenny, Jon. "Summary Of Known Tropospheric Delay Models," Canadian Marconi Company, the contents of each of which are incorporated herein by reference. This technique uses local measurements of temperature, pressure and humidity to derive a slant delay correction that is computed using the following technique: The range error in meters is computed as:

ΔR = AME) + _ZW(E)F_W(E) 4)

where Δ^ andΔ_ZH, are the dry and wet components, respectively, of the zenith value in meters, and

F_D (E) and F_w (E) are the dry and wet elevation sensitivities, respectively. These components are computed as Δ_ZD = 2.276 x l0^_3 5)

,1.23 ,146

Δz_™w = 470- + 1705α-^ 6)

F_D(E) =

. _, π 0. n0m014_3 7) tanZY 0.0445 I

^F^^{E) ~} . __. 00..0000003355 ⁸⁾ sm£ H tan£ + 0.017

7.617-™_ e₀ = 35.65 x RHx 10 ^τ 9)

In these equations, P is the barometric pressure in millibars, T is the temperature in

Kelvin, E is the elevation angle of the satellite, a is the temperature lapse rate (K/km), e₀ is the surface

partial pressure, and RH is the relative humidity (0 < RH < 1). Each of these values is determined by appropriate sensors on each WRS, which are transmitted to the WMSs along with the pseudorange data from the GPS satellites. Ionosphere correction

The GPS ionosphere correction is determined from the L2-L1 pseudorange difference measured at the WRSs. The GEO ionosphere correction is formed from the estimated grid delays determined in the Kalman filter. The process used to determine the GEO LI pseudorange ionospheric correction is the same as that specified in the Minimum Operational Performance Standard for Global Positioning System/Wide Area Augmentation System Airborne Equipment, RTCA Document No. RTCA/DO-229, January 16, 1996) , the content of which is incorporated herein by reference, for the airborne user ionospheric correction application.

The relativity, earth rotation corrections as well as the propagation of ephemeris information for the orbit determination and correction processing are performed as defined in the GPS Interface Control Document (ICD GPS - 200 Global Positioning System Standard Positioning Service Signal Specification, dated June 2, 1995). Multipath corrections

Although most of the pseudorange corrections are straight forward and common to many GPS applications, the technique of multipath deteirnination is unique. In WAAS, the WRSs are located at Air Traffic Control Centers with environments not very suitable for obtaining accurate GPS and GEO pseudorange measurements because of multipath, even when using receivers with multipath elimination technologies. Therefore, the WMSs may have to process pseudorange measurements with undesirable multipath errors. A common method of reducing the effects of multipath is to use carrier phase pseudorange measurements instead of code phase pseudorange measurements because carrier phase measurements have much smaller multipath errors (less than 5 centimeters). However, since the carrier phase measurements are ambiguous, code phase pseudorange measurements must be used to resolve the ambiguities. This is accomplished by smoothing the difference between the code phase and carrier phase measurements and adding the smoothed difference, which is the ambiguity, to the raw carrier phase measurements. This is illustrated in Fig. 10. The problem with this method is that the code and carrier pseudoranges diverge from each other over time because of the changing ionospheric delay, by a rate of twice the rate-of -change of pseudorange delay through the ionosphere. Thus, the smoothing time-constant is limited by the time-constant of the ionospheric delay. This time-constant is usually not long enough to always provide good multipath reduction for a stationary receiver antenna. Except for real-time integrity monitoring of pseudoranges, the WMSs always use ambiguity-resolved carrier phase pseudoranges for computing wide area differential corrections. However, to increase the time-constant of the ambiguity smoothing, the carrier phase measurements are first corrected for the ionosphere code/ carrier divergence to provide an altered LI carrier phase pseudorange that matches a multipath-free, noise free LI code phase pseudorange. This is accomplished by utilizing L2 carrier phase measurements in addition to the LI carrier phase measurements. This improved process is shown in Fig. 11 (γ is the ratio of

the squared LI and L2 frequencies, or 1.646944444444.)

The difference between LI code and carrier pseudoranges (sometimes called code/ carrier divergence), at time tk, is as follows (in meters): PR, _k -cφ,.,//_. = η_τ -cη_φ //, + μ_τ -cμ_Φ / , +2δPR_{I k} +cb, f 10)

where φι,jt is the carrier phase measurement (in cycles), subscript τ indicates code (in meters), subscript φ indicates carrier phase, η indicates noise, μ indicates

is the LI frequency, c is the speed of light, &_φι is the LI carrier phase ambiguity bias (in cycles), and 2δPRn ,* is the

ionospheric divergence between code and carrier at LI (in meters). The ionospheric divergence is cleanly measured with an ambiguity bias as (including carrier phase noise and multipath terms)

2c (φ.., Φ,,

2δPR, φ/

⁽ι-^γ) /.

2cb, 11)

= 2δPR, ₊- Φ2 ^2cfyι ^2criφi,t _| 2cμ_{φ2 t}

+- ^2cn Φ2,* 2cμ

Λ(ι-y) /,(ι-γ) /,(ι-γ) Λ(ι-γ) /,(ι-γ) Λ(ι-γ) whenever both the LI and L2 carrier phase measurements (φ and φ_2Jv) are available. Then,

Equation 10 is corrected to a pseudorange residual as

Δ Λ_.w = ««_...* -cφ //, -2δPR_φ/

^I_Φ , ^2cη u ^cτlφ2.* , i.* 2cμ _U ²cμ ₂,„

= ^τlτi,* 7— + 7; — ΓTT_I — s- + μ_τ — ^+

/, .(i-^γ) Λ(ι-^y) ^τ /, /,(ι-^y) Λ(ι-γ)

2c6 ₂ 2c6_φl c6 _l

⁺Λ(I-Y) /,(I-Y) ⁺ /,

where the noise and multipath terms and the three bias terms are lumped together into one ambiguity bias term plus single multipath and noise terms that will be dominated by that of the code phase pseudorange measurements.

The lumped ambiguity bias term is unknown. It is deteirnined by smoothing the one-second residuals of Equation 12 to arrive at

δP^~R_bi2.k_*y = (l - K_{bl2 lc} )δP^~Rbi2.k + K_{b k}APR_{τi i} 13)

where

δ R.n.o = Δ R._lι0 14) and the filter gains are

K_b{₂,k ⁼ TΓTT ' ^< limitl2

! 15)

~ l2, o — , . ' * — limitl2

^Λliraitl2 ⁺ * where /cim_ti2 is on the order of 600 (10 minute time constant) or more.

The LI pseudorange is then based upon ambiguity-resolved carrier phase measurements where

PR = {^ /Λ) ^{+ 2δPR}w.* +δ R 2.^„ 16)

The measurement variances of the carrier phase pseudorange measurements are also determined recursively. That is, following Equation 13, the measurement variance at time tk is

σ_b ² _{u l} = (l - K_{bi2 k} ) σ_b ² _{l2 k}

17)

where

^σ.12,0 = ^στl.O I⁸) and where the input noise variance is dominated by the code pseudorange noise variance and the code multipath error variance. That is,

The input noise variance is that provided by the WRS receiver along with the LI code pseudorange measurements. The input multipath error variance is an a priori estimate of that variance, determined from historic data, that is a function of elevation angle Ek. The multipath error is not independent from sample-to-sample, so it would not be smoothed to the degree that the noise error will be. Thus, its variance is increased to prevent an artificial reduction in the smoothing process by incorporating a multiplier (also a function of the elevation angle) that accounts for a variation in the correlation time constant of the error. This multiplier is shown in the brackets of Equation 19, which is valid for multipath error time-constants [l/Kmiφ,k(Ek)] that are small with respect to the smoothing time constant.

Thus, the present WAAS WMS Multipath Reduction Pre-Processing mitigates the problem of having WRS antenna environments not very suitable for obtaining accurate GPS and GEO pseudorange measurements because of multipath, even when using receivers with multipath elimination technologies. Since L2 carrier phase measurements are also available from the same WRS antennas, long term smoothing of carrier phase ambiguities can be accomplished to reduce these severe multipath effects. Ionosphere Estimation (306)

The ionosphere delay as an error source in GPS ranging is second only in magnitude to Selective Availability (SA), which is a technique used by the Department of Defense to reduce the ranging accuracy of the GPS system. The ionosphere delay is inversely proportional to the square of the frequency being transmitted and can thus be measured using dual frequency L1 L2 GPS receivers. Referring to Fig. 12, in order to meet CAT I landing requirements the system of the invention broadcasts ionosphere corrections at specified ionospheric grid points (IGPs) 400 over the WAAS coverage area. These ionosphere corrections are computed from data gathered from L1 L2 receiving reference stations 2 located about the coverage region.

The ionosphere estimation process uses an adaptive Kalman filter to estimate the error in the Bent model vertical delays and their change between estimation epochs at the selected IGPs. These error model estimates are then used to correct the computed Bent model vertical delay at the grid point to form an estimate of the Ionospheric delay at each Ionospheric Grid Point (IGP). The model for this Kalman filter uses the interpolation process provided in the aforementioned Minimum Operation Performance Standards.. While it is foreseen that a least-squares algorithm could have been used in addition to or in lieu of the Kalman filter, because some IGP delays may not always be observable, the state vector of the least-squares process would have to change dynamically. This unobservability can be handled much easier in a Kalman Filter, as can sequentially available measurements and adaptivity. Kalman Filter Implementation

The Kalman Filter, as described in R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering. Second Edition, J. Wiley & Sons, Inc., 1983, the content of which is incorporated herein by reference, is implemented for the ionospheric correction process as shown in Fig. 13. It begins with an initial estimate of the state vector and a corresponding initial error covariance matrix. When a set of measurements is available at time t£, the inner loop of Fig. 13 is exercised until all measurements at that time are processed. Then the outer loop is exercised to update to the time of the next set of measurements, at which time the inner loop is again exercised. Delay Estimation Measurements

The ionospheric delay estimation process uses dual frequency GPS semi-codeless P code pseudorange and carrier phase measurements collected at the WRS locations. Semi-codeless P code is described in A. J. Van Dierendonck, "Innovation: Understanding GPS Receiver Technology ~ A Tutorial," GPS World. Volume 6, Number 1, January 1995, pp. 34 - 44, the content of which is incorporated herein by reference. Because these semi-codeless measurements can be noisy with possible drop-outs due to low signal-to-noise conditions or scintillation, C/A code code/carrier divergence measurements are also used as a backup to the L1-L2 carrier phase measurements. These much-more- robust divergence measurements are the difference between the C/A code pseudorange measured from code correlation and pseudorange measured from the carrier phase. Because the ionosphere is a dispersive media that acts like a wave guide, the code is delayed (group delay) while the carrier phase advances. This is because, for such a media, the group and phase velocity of the signal are related as V =c² 20)

This difference measures twice the ionospheric delay, but with a carrier cycle ambiguity, making it somewhat useless as a delay measurement. However, it can be used to measure twice the rate- of-change of the ionosphere. Since the C/A code is used for this measurement, the signal-to-noise ratio is much higher than that of the semi-codeless measurements and much more robust in a scintillating environment.

Measurement Preprocessing

Prior to applying the measurements to the estimation filter, one-second L1-L2 pseudorange measurements are smoothed against the L1-L2 carrier to support the filter update rate of once per 30 seconds. This is sufficient for tracking the rate-of-change of the group delay during ionospheric storms. In the present system, IGP delays are preferably broadcast every 2 to 5 minutes so that there is no reason to update the filter at a rate higher than 30 seconds. C/A carrier/code divergence measurements are preferably collected every second, and then smoothed once per 30 seconds. Estimation Model

The estimation model consists of a state vector comprising Bent Model IGP vertical delay errors and their associated rate of change. The state vector is preferably that shown in Table 2, Equation 1.1, where τ/ is the Bent Model vertical delay error at IGP i at time tfc, {τAt)_ιk is the change in that vertical delay error (at IGP i between time updates) and N is the

number of pre-established IGPs. The measurement vector is preferably either that of Equation 1.2, if L1-L2 carrier phase measurements are available, or that of Equation 1.3 if the C/A code divergence measurements are being used, where

is the slant L1-L2 delay through IPPy (as shown in Fig. 12, an IPP 402 is an ionospheric Pierce Point, which is the point at which a signal from a WRS to a GEO pierces the ionospheric grid ) at time _V is the maximum number of EPPs in the 30 second epoch (that is, the total number of GPS satellites x the total number of WREs smoothed measurement received in the 30 second epoch), Fj is the obliquity factor for IPP

j as described below, x_Carjk i the carrier phase change since the last time update, Δτdivjk is the C/A code code/carrier divergence change since the last time update, t_„5e«t 'A: is the Bent model delay for the pierce point, and ΔxBentjk i the Bent model delay rate of change, all at time tk. In general, not all of these 2M (2M because in each epoch both a delay error and rate of change error are formed) measurements are available at each time step. Those that are not available are preferably not processed, i.e., nothing is done. Since the measurement errors are independent, they can be processed one at a time as scalars. Because the system processes the data as scalars, throughput of the system is improved.

Table 2- State Vector and Measurement Vector Equations

Computing Slant Ionospheric Delay

The vertical delay at a grid point or pierce point can be converted to a slant line-of-sight delay by multiplying that vertical delay by the obliquity factor ^, as

(^λ/> - Φff ) = ^FPP * >PP{^X ' *>PP) ²¹) where τ_vpp is the interpolated vertical delay at the user-to-satellite IPP derived as described herein, and

¹ R. cosE^{λϊ F}» ~ 1 - 22)

R_e +h,

The Kalman Filter is partitioned into two processes ~a time update and a measurement update. The time update is processed after the completion of all measurement processing at time tj_r. Time Update The state vector x# in the Kalman Filter includes all of the estimated IGP Bent model error delays and their changes, including those that are not currently observable, but will be at some time. The delay change estimates of the non-observable states are allowed to "coast" when the IGP is no longer observable, and declared as not-monitored by the WMS in the WAAS messages. The time update of the process is given by the "OUTER LOOP" path of Fig. 13, where # is the state vector made up of N IGP Bent model vertical delay errors and N IGP Bent model vertical error delay changes. The last N elements of the vector b are made up of a priori second differences of vertical delay at the IGPs over the time update interval. That is, the last Ν elements of the vector are made up of

"i*N.k ^{= τ} BaιtJ.k+\ ~ ^^τBauj.k ^{+ τ}Bmtj,k-l *^••>) where βent,i_,k i the vertical delay at IGP i at time t as predicted by the Bent Model. The Bent Model is described in further detail in S. K. Llewellyn and R. B. Bent, "Documentation and Description of the Bent Ionospheric Model," AFCRL-TR-73-0657, AD 772733, 1973, the content of which is incorporated herein by reference. The first N elements of b£ are always zero. This vector varies with solar activity (sun-spot number, solar flux), IGP location, time of year, and time of day. The state transition matrix is

Φ* = 0 24)

where I/ is an NxN identity matrix for N IGPs. The associated covariance update is also given in the "PROPOGATE" block of Fig. 13. P# is the estimation error covariance matrix prior to the time update and P_t'_tl is the same after the time update.

The broadcast IGP delays are computed as the sum of the first half of the Kalman filter state vector (the Bent model delay error estimates for each IGP) and the predicted Bent model delay at time tø for each IGP. A Priori Delay Model

As discussed above, the Bent Model is used as the basis of the ionosphere estimation process. This process takes advantage of the knowledge of the dynamics of the ionosphere. The b# is mainly used to keep the changes in ionospheric delay states close between updates to minimize the update uncertainty, and when not observed for better re-initialization after they again become observable. Measurement Update

A measurement is modeled as a linear combination of the four (4) IGPs surrounding the IPP of the measurement. It includes measurement noise and known WRE and satellite biases. All are divided by the appropriate obliquity factor to convert them to vertical delays. Equation 1.4 of Table 2 represents the L1-L2 pseudorange measurements, Equation 1.5 represents the L1-L2 carrier phase measurements and Equation 1.6 represents the C/A code divergence measurements. The FF-'s are the weighting functions described in the Ionospheric Pierce Point Delay Interpolation section below and i(n) denotes the IGP number related to IGP n surrounding the IPP. wu- i k is the L1-L2 measurement noise, bias - 2 is the sum of the known WRE and satellite L1-L2 biases, w_car , is the carrier phase measurement noise and wn_vjk is the divergence measurement noise. These measurement models map the prediction of the IGP delays into a predicted IPP vertical delay for the purpose of defining measurement residuals when differenced from the actual measurement. The Δτ_i_V jk/2 and W_di_V /jfc/2 quantities can be replaced with change in L1-L2 carrier phase and its noise, if available.

The measurement vectors hjk of the "MEASUREMENT PREPROCESSING" block of Fig. 13 are made up of the coefficients of the IGP states given in Equations 1.4 through 1.6 of Table 2. These map the vertical IGP delay states into the slant range EPP delay measurements. The equation for each measurement in a cell surrounded by IGPs is identical to the interpolation model used for ionospheric pierce point delay interpolation described below.

The associated measurement uncertainties are represented in

which is made-up of diagonal elements representing the variances of the individual smoothed L1-L2 measurements (first M) and the individual divergence or carrier phase measurements (second M). Both are converted to variances of vertical delay measurements. The first M variances will be dominated by the noise of the semi-codeless measurements. The second M variances will have an additive multipath component, if divergence measurements are used, weighted by the inverse of the tan² of elevation angle, as well as a small measurement noise component. If the L1-L2 carrier phase measurements are available, the delay change measurement variances will be reduced significantly. The noise variances are estimated using C N_o measurements accompanying the data from the WREs, which will effectively de-weight noisy measurements coupled with low elevation angle satellites. This will have the effect of influencing IGP delays that are already established by very little, while also providing IGP delays that have not been established, but with a larger uncertainty.

One of the features of the Kalman Filter processing independent measurements is that the measurements can be processed sequentially. Thus, if a measurement is missing because of a drop-out, it is simply skipped in the processing. Measurements are deweighted in the filter gain computation by the expected variance of the residual jk evaluated as shown in the "EDIT & COMPUTE KALMAN GAIN" block of Fig. 13, where rjk is they^'th diagonal element of Rjfc. The gain matrix for thisy'th measurement is also shown in that block.

The updated covariance and state vector are shown in the "UPDATE STATE" and "UPDATE COVARIANCE" blocks of the inner loop of Fig. 13, where

-

for L1-L2 measurement zn-L2jh ^an<^

for carrier phase change and divergence change measurements, respectively. The equations in the "UPDATE STATE zjk" block of Fig. 13 can be rewritten as

x_k = x'_k +k_Jk(_Zjk -h_jkx_k' )

where the terms within the brackets are known as "measurement residuals" or "innovations." Their associated variances are the α ^'s. These quantities are saved and used in the "ADAPT PROCESS

NOISE" process. The L1-L2 versions are archived for use in WRE and satellite L1-L2 bias error estimation processes. The covariance and state vector are updated as shown in the inner loop "COPY" block between each measurement update, which is the same as a zero-second time update. The inner loop computations are repeated until all measurements collected at time tk are processed. Filter Initialization The initial IGP delay and its uncertainty are determined from the Bent model. That is, the initial states i are

^Xι,0 - ^ B ttJ.O 29)

* 30)

and the initial error covariance matrix is

»_'_f-___Y__'x' ^T 31)

⁰ liooj °^x°

Adaptive 0 Determination

The Kalman filter implementation, as shown in the Ionospheric Delay Estimation Kalman Filter (Fig.13), amounts to using "postcorrection" measurement residuals for Q adaptation. Also, ε is used for the residual ratio. The adaptation process is performed in parallel with the Kalman filter using measurements collected during the period after the state and error covariance are updated using the previous period's measurements. After some filtering, the residual ratios are collected into cells to be combined for adaptation of the rows and columns of the Q related to the IGPs surrounding the cells, as shown in the Example Grid and Grid Cells for CONUS (continental United States) of Fig. 14.

The square of the pseudorange residual ratios at time t , for n_m measurements in Cell w, is as follows:

where

or

ε ._t _ = ^■(_____________ . ,_. . "m 34) ^h_.W ^+,

Similar equations apply to the n_m phase change measurements (carrier phase change or

divergence change), providing 2n_m squared-residual ratios spRjk and ε^R jk- The non-weighted

pseudorange residuals are also used in the GIVE computations discussed below.

In order that the adaptation of Qt does not respond to noisy conditions, the 2n_m squared-

residual ratios are then smoothed with a time constant ar k that varies with the variance of the

measurement. That is, the larger the variance, the longer the time constant. The filter is a first-order filter of the form

'PR.jk 'PR.jk 35)

and

where

^ PR.jk ~ ^E PR.jk 37)

and

^ PR k ~ ^EάPR k 38)

whenever a given measurement is available the first time at time t£. The scale factor a is a predetermined

constant based on filter performance analysis, and Δt is d e time between measurements. If rjk, the variance of the measurement, is constant, the filtered variance is the variance of the squared-residual-ratio divided by 2arjk-

The 2n_m smoothed squared-residual ratios at time tk are then averaged over the n_m samples to provide one residual-ratio number for the Cell m as

1 "" ^E-Pk ~ -ι ^Z PR.jk 39) ⁿm y=l for the pseudorange measurements, and

1 ""

^E'&PR.mk ⁼ 2-_t ^SWR,jk 40)

for the phase change measurements.

These residual ratio numbers are then limited as l/E_lim < E_M,mλ < E_lira 41)

and l/E^ ≤ E^ ≤ Eπ. 42)

where Eι,_m is tentatively set to 4/3. This limiting prevents any abrupt changes in the Qk- Determination of the Cir and Dfr The Qk are part of the adaptive process that adjusts the statistics of the filter to respond to ionospheric model deficiencies and measured anomalies, where

where the Q,,_k are defined as follows:

Q22* = »kQben,kVk 44)

Q_ι = C_tQ_J24C_t 46)

where

Qbe ,,k = (z^ t>2kt>^T2k + Q_m 47)

where b₂£ is the second half of the b vector with elements defined in Equation 23 and/? is a percentage error expected in the Bent Model (-25%). Ck and D# are diagonal matrices used to adapt the Qk described below and Q_mιn is a minimum diagonal matrix that prevents zero elements when the Bent Model has no dynamics.

At this point, the adaptation scales the new Q (using the change measurement average squared-residual ratio as an example) based upon what appeared to be post-correction measurement error variance mismatches on the last iteration, given as

"_PΛ,m*'^m* 22.*-lΔ"m*^αΔ/'Λ,-* _ p «, ⁿΛ/^,Λ,n_"22,*-l ^DΔPΛ,mt where the i_eJ>R_tmk interpolates the IGPs surrounding Cell m to the centroid pierce point and the Δ _mk is a diagonal matrix with entries of Is except for those elements associated with the 4 IGPs. This is so that only the elements of Qbent.k associated with those 4 IGPs are changed. The solution for the ΔD,~£ i ^not deterministic since there are 4 unknowns and only one known (Ε._ό R_>mk)- Even if one built up a matrix of these scalars covering the entire grid, there would always be more IGPs than cells. Thus, the result would always be not deterministic. Thus, just as with the GIVEs below, some variables must be constrained so as to be able to determine an answer. Accordingly, the four non-unity elements of a diagonal matrix AD_mk are set to

δd_mkj = ^E^_mk 49)

so that D_n, = ΔD_miD_m,„_, 50)

where

*>_mo = 51)

This approach applies equal weighting to each of the 4 IGPs. The approach is applied to all M cells, so that

If any cell has no pierce points (or squared-residual ratios), then set D km ⁼ IN- This same technique is applied to the pseudorange squared-residual ratios, where

^hP*.m*^ΔC-*D*Q22.A-l^D*^ΔCm*^h *._m* = E PR.mk 53)

" PR,mk **Q 22,*-l "k " PR.mk

The four non-unity elements of a diagonal matrix AC_mk are set to

^δCmS = PR.mk 54)

so that

*-"m* ^{— Δ}^-m* _{m t}_, 55)

where

After applying to M cells, then

If any cell has no pierce points (or squared-residual ratios), then set km = IN These new Ck and D# matrices are then applied to the new Qbent.k at time tk as described in Equations 43 through 47 above. Pierce Point Location Determination

The following process describes the weighting and pierce point functions that are used throughout the ionosphere estimation process. Considering the satellite and user locations, the user must first determine the location of the Ionospheric Pierce Point (IPP) of the signal path from satellite to the receiver. A technique for determining this location is provided in J. L. Junkins, G. W. Miller and J. R. Jancaitis, "A Weighting Function Approach to Modeling of Irregular Surfaces," Journal of Geophysical Research. Volume 78, No. 110, April 1973, the content of which is incorporated herein by reference. Ionospheric Pierce Point Delay Interpolation

Although the data broadcast to the user is as vertical IGP delays, these points do not generally correspond with the user's computed IPP locations. Thus, it is necessary for the user to interpolate from the broadcast IGP delays to that at the user's computed IPP locations as shown in Figure 15. Given four nodes of a cell of the IGP grid described above that surround the user's EPP to a satellite, the user can interpolate from the nodes to the user's pierce point using the following technique. A weighting function approach for modeling irregular surfaces provides a simple procedure for approximating an irregular surface from regularly spaced ionospheric grid point vertical delay data, as discussed in the aforementioned "A Weighting Function Approach to Modeling of Irregular Surface.." The mathematical formulation for interpolated vertical IPP delay ι_Vpp §pp, pp) as a function of IPP latitude φ-, and longitude λpp is

(ΦpYpp) = Σ> pp>: )^τw ⁵⁸)

where the general equation for the weighting function is

W(x,y) = x²y²(9 - 6x -6y + 4xy) 59) and τ_v; are the broad cast grid point vertical delay values at four corners of the IGP grid, as shown in Fig. 15. In particular, x_vpp is the output value at desired pierce point pp, whose geographical coordinates are φp , λp_p,

W,(x,y) = W(x,y) 60)

W_l(x,y) = W(\-x,y) 61)

W,{x,y) = W{l-x,\-y) 62)

W_t{x ,y) = W(x,\ i-y) 63)

* ^■„ ^{= λ}PP ~ y 64)

Δλ Δλ»

^XPP - 66) λ₂ -λ, longitude grid interval

Ionosphere Bias Estimation (304)

L1-L2 biases can exist in both the WRS receivers and GPS satellites. Since the biases are essentially constant (except for configuration changes), they will appear as biases in the system of the invention. Accordingly, the residuals from the Ionosphere Estimation Kalman filter are used in a separate Kalman filter to estimate these biases in a manner similar to the techniques described in 1) B. D. Wilson and A. J. Mannucci, "Instrument Biases in Ionospheric Measurements Derived from GPS Data," Proceedings of ION GPS-93, Salt Lake City, UT, September 22 - 24, 1993, pp. 1343 - 1351 ; and 2) B. Wilson and A. Mannucci, "Extracting Ionospheric Measurements from GPS in the Presence of Anti-Spoofing," Proceedings of ION GPS-94, Salt Lake City, UT, September 21 - 23, 1994, pp. 1599 - 1608, the contents of each of which are incorporated herein by reference. While in these references the biases were included in the ionospheric delay state vector, the large difference between the time constants of these biases and the ionosphere allows the state vector to be partitioned into two state vectors, and thus, two Kalman filters.

The bias estimation process is based on the ionosphere estimation residuals in the following way. The ionosphere estimation residuals have the form of: ε^ SV^ + WRE^ + η

for a non-gold standard receiver.

The residual for a measurement made from the gold standard receiver (described below) is

^_} = SV_buuJ + η

where ε_tJ is the residual, SVbias.i is the satellite bias, WRE bias J is the reference station receiver

bias, and h is Gaussian noise.

The Bias Estimation Kalman filter removes the noise term. The use of the gold standard receiver allows the ambiguity in the first equation to be resolved since the second equation will only have a residual associated with a satellite. The filter based on the residuals with only satellites can then resolve the bias ambiguities in all measurements by referencing all biases to the gold standard receiver. Accordingly, in the system of the invention, a single Gold

Standard receiver is calibrated periodically and is used as the reference in the Ionosphere Bias

Estimation Filter.

Satellite Orbit Determination (308) Since the system of the invention uses GEOs as ranging sources, it is necessary to compute ephemeris and almanac information that can be uplinked to the system users so that the user can determine the satellites' positions. Also due to message bandwidth limitations, the broadcast corrections consist of slow and fast corrections. The received data must be separated to allow computations of these two forms of corrections. Finally, a method was developed to check the integrity of the broadcast GPS ephemeris information by comparison to the computed orbits from the satellite orbit determination process.

The Satellite Orbit Determination process is used to compute the GPS and GEO satellite navigation information. This data consists of estimates of the GEO and GPS position and velocity vectors and satellite clock offset and drift states. The system also provides reference station clock bias, and drift clock states as well as tropospheric residual estimates. The implementation of the Satellite Orbit Determination Process is a near real-time extended Kalman filter. The term near-real time is used since the update rate of the Kalman filter is preferably 5 minutes using data collected and smoothed from the previous five-minute period. The preferred method of orbit determination is to the RTOD™ developed by Logicon Ultra Systems. A batch type orbit determination process could also be used but this would not allow real-time detection of anomalous satellite behavior. This failure to detect an anomaly would constitute a form of HMI (discussed below).

The orbit determination data is used as the basis for the long term corrections, as the source of the GEO ephemeris and almanac information, and for the .GEO and WAAS Network Time ("WNT") clock steering processes. Long Term Corrections (310)

WAAS must compute long term corrections (slowly varying errors in ephemeris and clocks), GEO almanacs and ephemeris. The long term corrections are formed by propagating ahead in time (preferably 15 minutes from the present time) the GPS and GEO broadcast navigation data and the computed GPS and GEO navigation data. The computed GPS and GEO navigation data values are derived from the satellite orbit determination process. These values are then differenced to form the WAAS long term ephemeris and clock corrections. The computed GEO navigation data is also used as the basis for the GEO ephemeris data and the GEO almanac data. Also, the GEO data is propagated for up to 12 hours and sent to the GUS for use in the autonomous ranging mode (discussed below). Fast Corrections (314)

Ultimately, the purpose of WAAS in a precision approach environment is to provide the WAAS users with accurate wide area differential conections consisting of longer term ionospheric model corrections, longer term satellite ephemeris and clock corrections and shorter term satellite clock corrections. The WAAS user collects the long term corrections from one GEO satellite and applies these corrections to pseudoranges along the line-of-sight (LOS) to the measured satellite. (See Fig. 6)

The WAAS user cannot mix corrections from multiple GEO satellites for two reasons - their long term corrections are not synchronized and may differ, and short term corrections maybe computed from a slightly different time scale. These long term corrections are fixed (or slowly propagated) for their period of applicability. The user then applies the most currently valid fast corrections (the user cannot use fast corrections that are older than those for which there exists integrity information). These fast corrections must be collected from the same GEO. The user expects that these fast corrections match the long term corrections that were computed from collections from the same GEO. That is, when the fast corrections are combined with the long term corrections, the combined correction provides the accuracy specified in the integrity information that the user also receives from mat same GEO. The integrity information is only valid for corrections received from that GEO. Thus, the WAAS WMSs must generate accurate and timely fast corrections that match the longer term corrections and the integrity information that accompanies the fast corrections.

Because the phasing of the data provided by each GEO does not necessarily match that of other GEOs, or the fact that the already broadcast longer term corrections could have been generated at another WMS, each WMS computes fast corrections for each GEO independently based upon the long term corrections being broadcast by the applicable GEO. The WMS does this by applying the same process applied by the user except that the fast corrections are generated instead of applied. This is shown in Fig. 17. Fast Correction Generation

The fast corrections, which absorb the errors of all other corrections previously applied, are computed for each GEO footprint. The fast correction for satellite i is defined as

I n, c sv = ψ ^{+ δα} _/>(' " L . ∑^ΔPR,, 68) ⁿi *=1 where the first term is the "currently broadcast" (CB) long term clock correction, ak is defined below, «/ is the number of WRSs (with valid measurements) in view of the GEO from which the currently broadcast corrections came, and

APR_ik = PR_ik{t)-R_{ik gps}{t)-l_ik .

+cAt_maji

where, for WRS k, PRik is the multipath eliminated carrier phase pseudorange described above, Rik,gps is the range from SV i to WRS k computed from the known WRS location and the decoded GPS ephemeris data (with the same Issue of Data (IOD) being broadcast from the GEO), l £ is the line-of-sight unit vector from WRS k to SV i, which is dotted into die long term satellite position correction vector that is currently broadcast, Atsv_,i_,gps ^ΪS the clock correction from the decoded GPS navigation data (with the same IOD), ARβR_jk is the Earth's rotation correction, cAτt_r0p_ti is the tropospheric delay correction computed from meteorological sensor data, cδAτ_V(_rθp k is the vertical tropospheric delay residual at WRS k estimated in the orbit determination (OD) estimator, AϊψBS.k i the WRS k clock offset estimated in the OD estimator and cAxi_ono is the LI ionospheric delay computed from the currently broadcast Ionospheric Grid Point (IGP) data for die pierce point between WRS k and SV i. Weighted Error Computation The system of the invention de-weights measurements based upon both prior and measured variances of the measurement errors. For example, the Ionospheric Grid Point delay estimation technique and the orbit determination technique are both Kalman filters that always de-weight measurements based upon their error variances. Likewise, the fast correction determination uses a weighted average over all WRSs (in view of the applicable GEO) based upon the measurement error variances.

Using weighting, the statistical error samples become normalized, where APR!_k = a_ikAPR_ik ; k= \, ..., ni 70)

where the weights are defined, for measurement PRik error standard deviation σ_η/ for satellite i, as

which have the property

£«_ = «, 72) t-l

The statistics of the new weighted samples are

^σ«_ = ^aWpR_lk ⁷³) and μ_PR;k = a_ikμ_PRΛ 74)

In general, in the long run, the mean of the errors is zero.

The purpose of die WMS Fast Correction Generation Process is to minimize die WAAS

User's total differential correction error, not to minimize the error in the current WAAS error estimation process. Fast Corrections are computed on a "per GEO" basis after applying currently broadcast longer term corrections because that is how the WAAS users apply the corrections. Furthermore, die corrections are also more accurate because they are generally applied to a smaller region and don't have to absorbed large spatial decorrelation errors associated with an entire WAAS service volume. The smaller region is limited to the union of me WAAS service volume and the applicable GEO footprint. UDRE Determination (316)

The WAAS System Specification and RTCA WAAS MOPS Signal Specification require the broadcast of User Differential Range Error (UDRE) parameters that bound the differential pseudorange error of the combined fast corrections and long term satellite corrections with a probability of 0.999. Furthermore, for integrity purposes, tiiese bounds must reflect (upward) changes within a time-to-alarm of 5.2 seconds. In mat short period of time, it is impossible to obtain enough independent error samples to determine such a bound (with a probability of 0.999) without being extremely conservative - so much so mat precision approach availability and continuity would be significantly compromised because of the larger conservative UDREs. On the other hand, increasing die number of independent error samples will increase me amount of time it takes to determine the UDRE, which may not give an indication of current differential correction error.

There are a number of approaches to solving die problem of not having enough independent error samples. One approach is based purely on non-parametric statistics using no prior information. The criticism of this approach is diat it takes too many samples and is too slow and far too conservative. This is typical of any non-parametric (distribution free) approach. The next less conservative approach is to assume an error probability distribution form (such as normal), compute sample means and variances and determine error bounds. Wi i tiiis approach, die only prior information used is that of the assumed distribution form. This approach also tends to be quite conservative unless a large number of samples are available (but not as large as required for me non-parametric approach). Another approach is to base die computation of UDREs using post-correction error samples as well as a priori (prior) covariances. The problem wi i this approach is how to combine these two sets of data so diat bod die bound probability and the time-to-alarm are met. The system of the invention uses and alternative approach known as Bayesian Statistics, Bayesian Inference or Bayesian Data Analysis. In the Bayesian approach, not only is the error distribution form assumed, but mere is some knowledge of the enor distribution statistical parameters, such as die mean and variance. However, even ese "known" parameters have a certain amount of uncertainty assigned to d em (as probability densities). That is where post-correction (a posteriori) error samples come in to play. These a posteriori (posterior) error samples are used to update or defend the prior assumptions. In general, die Bayesian approach is less conservative tiian die above mentioned approaches.

The Bayesian solution solves the UDRE determination problem very well in that there is definitely some knowledge of the prior statistics as estimation error covariance matrices. These covariance matrices not only reflect me geometry and measurement availability, but are themselves somewhat dependent on die posterior data in that die estimation processes are adaptive. That is, mey adapt to measurement residuals, but widi some time constant delay. Furthermore, prior statistics also include the effects of measurement error variances.

In addition to using prior statistics, post-correction error samples are weighted inversely widi die known, measured error variances so diat die UDREs are not penalized witii measurement errors that have notiiing to do widi the signal-in-space errors, except as tiiey affect die error estimation process. This is consistent witii error estimation processes that always de-weight measurements based on tiieir error variances, eitiier using weighted averages or via Kalman Filters. This is especially important when using low elevation satellites for correction parameter estimation. Definition of Error Bound

The basic definition of an Error Bound B is independent of the statistical approach used in estimating it. Only d e value of die bound will have a different level of conservatism attached to it depending upon d e approach. Usually (except in the distribution free case), 5 = |μ| + r^*S 75)

where μ is the sample mean and S is the square root of the sample variance over n enor samples and I* is

die 99.9% probability factor determined from the underlying probability distributions, where, for n error samples *;, witii equal weighting,

μ = ~∑x_i 76)

and

S^ -^∑lΛ - μ)² 77)

The probability factor F is a function of the number of error samples, assumed sample probability distribution, and, in the case of Bayesian Statistics, prior distributions, which are usually a function of probability densities of a prior mean μo, a prior standard deviation σ_o and other parameters describing their uncertainty. If there were an infinite number of error samples, Equation 75 would be approximately |μ| + 3.29σ assuming diat the error samples were drawn from a normal distribution.

In die distribution free case, die bound is a function of ordered statistics, where the error samples are ranked and d e bound is a function of rank. For a high probability bound witii a limited sample set, such as the WAAS case, this usually means that the bound is die worst case error sample. For a 99.9% bound, tiiousands of enor samples are required in die distribution free case.

An example of die next less conservative approach, which assumes an error probability distribution form, is shown in Conker, R. S., El-Arini, M. B., Albertson, T. W., Klobuchar, J. A. and Doherty, P. H., "Development of Real-Time Algorithms to Estimate the Ionospheric Error Bounds for WAAS," Proceedings of ION GPS-95, Palm Springs, CA, September 12 - 15, 1995, the content of which is incorporated herein by reference, and applied to a bound for die Grid Point Ionospheric Vertical Eπors (GIVEs). In that case, a normal distribution was assumed with an unknown mean and unknown standard deviation. In that paper, the authors based the probability factor on the techniques presented in Hahn, G. J., and Meeker, W. Q., Statistical Intervals: A Guide for Practitioners. John Wiley & Sons, Inc., New York, 1991; and Odeh, R. E. and Owen, D. B., Tables for Normal Tolerance Limits. Sampling Plans. and Screening. Marcel Dekker, Inc., New York, 1980. For this case, the probability factor T* takes on a value of 5.43 for a error sample size of 30 for a 99.9% bound. This value is quite large considering that, for a normal distribution, that μ ± 3.29σ bounds 99.9% for an infinite error sample size. Unfortunately, die probability factor reduces very slowly for larger finite error sample sizes, even if they are available. To approximate the value of I as a function of n, use the following function of the statistics of the sample mean and variance for 99.9% bound (assuming samples drawn from a normal distribution):

T^' « -J m² + (³-²⁹ _>/(« - / xo._w(,o- ._B-i )^{2 7g)} where t _.999b.n- is the 0.9996th quantile point of the t distribution (the distribution of die sample mean)

with n-1 degrees of freedom and X_{o.99 (},_o-__),„-i is the 0.999(10-6) quantile point of the chi-square

distribution (the distribution of the inverse of the sample variance) with n-\ degrees of freedom. The digits b and 10 - ό (l < b < 9) provide a tradeoff for apportioning the total 0.999^th quantile to the sample mean and standard deviation, and can be chosen to minimize T*. The first term in the brackets represents the uncertainty in the sample mean, which is usually dominated by the second term, but not always. The value under the radical of d e second term represents the uncertainty in die sample variance. If the uncertainty in the sample variance is much higher than that of the sample mean, a much higher quantile can be assigned to the bound on the sample mean than is assigned to the sample variance, resulting in an joint quantile equal to, for example, 0.9999 x 0.9991 « 0.999, since the sample mean is independent of the sample variance. For n = 30, Equation 78 yields 5.437, not far from die value used in the aforementioned "Development of Real-Time Algorithms to Estimate the Ionospheric Error Bounds for WAAS." For smaller sample sizes, d e results of Equation 78 become quite large. For example, using n = 6, representing a UDRE at a worst case location over a 6-second interval, Equation 78 yields 17.02. This is unacceptable. Using n = 100 yields 4.217, suggesting that the value does not diminish quickly as die number of samples exceeds 30. Thus, it is desirable to develop a more optimistic bound determination, given that prior information about the statistics of die errors is available. Bavesian Statistics

Bayesian Statistics provide a method for accounting for prior information. Bayes' Theorem is described in Box, G. E. P. and Tiao, G. C, Bavesian Inference in Statistical Analysis. Wiley Classics Library Edition, John Wiley and Sons, Inc., New York, 1992, the content of which is incorporated herein by reference, and generally provides as follows: Suppose that x =

..., x_n) is a vector of n observations whose probability distribution p(x\θ) depends upon die values of k parameters Θ = (θi, ..., θβ . Suppose also that Θ itself has a probability distribution p(Θ). Then,

Given the observed data vector x, the conditional distribution of Θ is

Also,

p{x)

= c-* = J/>(x|Θ)/>(Θ)_ 81)

where E[ (Θ)] is the mathematical expectation of/(Θ) with respect top(Θ). Thus, Equation 80 becomes

Equation 78, or equivalently, Equation 82, is usually referred to as Bayes Theorem. In this expression, p(&), which describes what is known about Θ without the knowledge of the data, is call die prior distribution of Θ. Correspondingly, ρ(@\x), which describes what is known about Θ given the knowledge of the data, is call the posterior distribution of Θ given x. The quantity c is merely a "normalizing" constant to ensure that die posterior distribution p(Θ\x) integrates to 1.

Ifp(Θ|x) can be determined for some prior distribution of Θ, then certainly a bound B can be found such that

The key here is that tiiere is a knowledge of die prior error mean μo (normally zero for die WAAS estimation processes) and a standard deviation of the enor, σ₀, each with some uncertainty, from the error estimation processes and measurement error estimates. Because of the nature of these processes (long time constants and the properties of measurement noise), we can reasonably assume that the errors are normally distributed, at least in die long term. Probability densities such as given in the above equations have been derived for reasonably defined prior densities for μ and σ², given μo and σo, as shown for example in Gelman, A., Carlin, J. B., Stem, H. S., and Rubin, D. B., Bavesian Data Analysis. Chapman & Hall, London, 1995, die content of which is incorporated herein by reference. However, Equation 78 is likely take on a different form, where

_? = β(μ,S,σ^μ₀,H,0.999) 84)

plus any other parameter used to define uncertainties in the prior densities. For example, infened from using a prior normal distribution of the mean wid a mean of μo and a variance of σo/κ₀ and a prior scaled-

inverse-chi-square distribution of the variance with vo degrees of freedom and a scale of σo, yielding the bound

κ„ _{+ T'} I ^vo^σo ^■ + , (" " S² , κ₀«(μ - μ₀)

85)

« + κ₀ -μ₀ + - n + κ_n v₀ + « v₀ + « (« + κ₀)(v₀ + «)

which is an estimate of the function B = Iμl + 3.29σ , assuming a normal distribution, where now

T ^W Λ/ ^0.999*,_K„ / 'V ^K» ) ⁺ (3-^{2 V}. / XO.999(!0-4),_V. j 86)

v_o and κo defines the "spread" or uncertainty in die prior variance and mean, respectively, — the larger they are, the less spread there is. Note that the quantity under the radical is a weighting of the prior variance and the posterior error sample variance, plus an enor sample mean-squared term representing the uncertainty in the prior mean. For larger values of vo and κo, which is reasonable for WAAS application, this method has promising performance for me small error sample sets diat are necessary for a short time-to-alarm. The new estimate of the variance σ_π ² is die quantity under die radical witii degrees of freedom v₀ + n. The new

estimate of the mean is μ_n, whose absolute value is the first term of Equation 85, has the variance

°7(^K _o +")-

The parameter κo is used to define die uncertainty in the prior mean to be a percentage of that of the prior error distribution itself. The parameter v₀ is used to define die uncertainty in the prior

variance, where

E(θ) = - al , for vo > 2 87)

, -²

and

For example, using values v_o = 12 and κo = 4 give a standard deviation of the prior mean of

0.5σ₀ and a standard deviation of the prior σ² of 0.6σ² yields die bound

Assuming there are 30 enor samples witii a prior mean of zero to compare with the results of the aforementioned paper entitled "Development of Real-Time Algorithms to Estimate the Ionospheric Enor Bounds for WAAS" (but using quantiles of 0.9999 and 0.9991 for die mean and standard deviation contributions), results in

B = |0.882μ| + 4.94^0.286σ² + 0.690S² + 0.084μ² 90)

which is much more optimistic, provided that the prior variance is less than the posterior enor sample variance and error sample mean-squared.

Assume there are six error samples widi a prior mean of zero to emulate a fast response for UDRE at e worst case location, results in

B = |0.6μ| + 6.65^0.667σ² + 0.278S² + 0.133μ² 91)

which relies more heavily on the prior variance. However, even if the error sample mean and variance were zero, the bound would be 5.43σo, which is still quite conservative. If the prior variance were small, however, d e reliance of the error sample mean and variance would be quite reasonable.

The values of vo and κo are configuration parameters that are adjusted based upon experience.

UDRE Computations UDRE is based upon the statistics of die enors in the fast conections described above. The associated total conection enor for usery, who is not necessarily close to WRS k, including the ionospheric model enor, is then δPR_y = cδAt_{S¥Λ ast} + δR_{e≠ IJ} + cδΔτ,_ono„ + δPR_{latency ι} 1 I ^ ' r

= -∑^fl η,_«*,„ -(δ^_A,_t + cδτ_{Mp jfc} - cδΔt_{TO5 ιl})_oo - cδΔτ,^

"ι *=1

H^δ _0β ^+cδΔτ-._y ^+δ < latency,! 92)

\ cδΔτ_{(0Λ0 V} ∑α,*cδΔτ tono.ik + δPR l,atency,!

', *-l

representing four basic enor sources — measurement enors, OD enors, ionospheric delay enors and latency due to Selective Availability (SA), including die spatial deconelation enors due to die ephemeris and the ionosphere. This UDRE computation process is illustrated in Fig. 18. Equation 92 includes the effect of the ionospheric model enor as seen by die user.

Equation 92 provides the basis for defining the prior statistics for the Bayesian Estimator based upon known measurement enor variances, OD covariances, the IGP vertical delay estimation filter covariances and suspected latency enors. However, because UDRE is defined as die bound at die worst case user location, over a short 6-second Fast Conection update interval, the second and third terms of Equation 92 represent bias-like enors that are not independent from sample to sample. In fact, except for some anomalies, they would not change from one 6-second interval to the next. Thus, at each user location, the enor should be treated as having a mean value with an uncertainty.

Including the total ionospheric model enor, die prior pseudorange variance for satellite i at a

user location^ is

2 , 2 ^σ R,-J = ~ Σ ^flU^ση.~ ⁺ ^latency,! + ^CT eph,tropo, clock, i tono j 93)

*-l where σ^ik is the standard deviation of the measurement enor. The last two terms represent bias

uncertainties. They are not directly measurable, except for users that would be located at the WRSs. Thus, to that degree, they represent spatial deconelation enors. Based upon the estimation filters implemented in the WMS, the prior variances for these last two terms are described here. Prior Ephemeris. Troposphere and WRS Clock Enor Variances

The prior variance of die total user enor due to the ephemeris enor combined with die effects of the ephemeris enor, me residual troposphere and WRS clock on the Fast Conection enor is

f ² eph,tropo, clock, ij

^{= i};^cxA ⁺ — ∑Σ^a,k^aιl^h,k X,_u 94) ⁿ, k=l 1=1

where

h_rt = [i; 1/sinE,, -l] 95)

where 1 (1 jfc) is the line-of-sight (LOS) unit vector between user/ (WRS k) and satellite , and

δX_y

96)

is a 5x1 vector of satellite position estimation enors and WRS vertical troposphere and clock estimation enors. C_x is the 5x5 covariance matrix of the position enor of satellite i and WRSsy and k

tropospheric delay and clock enors from the OD filter, where

where C_x is the 3x3 satellite i position covariance matrix, C_xw is the 3x2 cross-covariance matrix for

satellite and WRS k and C_w is the 2x2 cross-covariance matrix for WRS k and WRS /. If & = /, this

latter matrix becomes the 2x2 covariance matrix for WRS k. Note that the variance of Equation 94 models the spatial deconelation enor between the user and the weighted-average over all WRSs. This weighted-average enor is absorbed in the satellite clock conection. Prior Ionospheric Delay Conection Enors

Since UDRE enors are conelated to the ionospheric model enors, these model enors must be considered here. This is because the weighted average of all WRS ionospheric pierce point enors using the broadcast IGP delays is also absorbed in the fast clock conections for satellite i. Other effects, such as ionospheric spatial deconelation is assigned to the GIVEs. Thus, given the effect of the weighted average, the remaining pierce point enor for usery to satellite i is

where Fy is the obliquity factor for the pierce point between location./ and satellite i, the W^ s are the interpolation weighting factors defined in the Minimum Operational Performance Standards for the Global Positioning System/Wide Area Augmentation System Airborne Equipment, Appendix A, RTCA Document No. RTCA DO-229, January 16, 1996, the content of which is incorporated herein by reference, cδiQp_n ij is the enor in the IGP n vertical delay estimate for those IGPs that sunound the

user's pierce point ij and cδτgPn',ιk is the enor in the IGP « 'vertical delay estimate for those IGPs that sunound the WRS's pierce point ik. The conesponding prior variance of the ionospheric conection effect on the user's Fast Conection, including the effect of the "weighted-average" ionosphere on the weighted solution for the satellite clock conection, is then as follows:

4

= ∑Σ^W _nW_lC_{CPnUj +}-∑∑a_dca_ιmF_lkF_m∑∑W,W_rC_CP,_l, ,ιkm 99) n-1 _I i=l m=l n'=l /'=!

^" - Σ ^a*^F> Σ Σ ^W (C_CPnn;,jk + < „ )

'l t=l B=l n'-l the cδτQPn s are the vertical delay conection enors at the Grid Points that sunound pierce point ij, «_ is the number of WRSs in the GEO footprint observing satellite i, CQp_n[ H is the enor covariance between Grid Points n and / for those Grid Points sunounding pierce point ij (if n = /, it is the variance), and CQPnj. km i the enor covariance between Grid Points n' and /' for those Grid Points sunounding both pierce points ik and im. The first of the three terms is the prior ionospheric enor variance for user pierce point ij, and it will be assigned as the prior variance for GIVE evaluation. It is not part of UDRE. The second term is the prior variance of the weighted-average slant pierce point enor over the « WRSs. It is part of the prior variance of UDRE. The third term makes up the cross-covariance between the user's pierce point enor and the weighted-average of all WRS pierce point enors. It could be assigned to either the GIVE or to UDRE or, conservatively, assigned to neither. Since the difference between different Grid Point enors are certainly positively conelated, this third term is negative. Therefore, if it is neglected, the end result is conservative. Certainly, for WRS pierce points far from the user's pierce point, the term is negligible, since the cross-covariances are small. Assign one-half of it to both GIVE and UDRE. Cross-covariances diminish with distance, thus, the fact that it is subtracted never causes the resulting variances to be negative.

Assigning die sum of die second term and ti ird term to UDRE, yields

1 "' "' 4 4

° IOΠOJUDREJJ ⁼ T _ -( ^{a a}m^Fιk^ιm 2__ _ _. "n^,"^l^GP ^,l t ⁿι *-l m-I n'-l /'-I . «„..

Treatment of Bias Enors

At the WRSs, die enors represented in Equations 94 and 100 are biases that can be measured wid noisy measurements. That is, the enor at WRS k is simply

ε«.,_ = μ_t + η_ 101) where, once measured, has a standard deviation of σ^ik- Furthermore, the results of the previous 6-second interval are used to provide a prior mean and variance to provide a Bayesian estimate of d e bound. However, it is appropriate here to keep intervals independent of each otiier to detect short term anomalies. It is also appropriate to first estimate die bound for each WRS location based upon Equation 94 and tiien add on the effects of spatial deconelation from that location. The next step would be to then pick die worst over all WRSs, including die spatial deconelation. Since the latency enor is the same at all locations, it could be added in at any time. Prior Mean Enors

In order to detect short term anomalies, here, it is assumed that tiiere is no prior mean information. To fit tiiis into the structure of the Bayesian approach, that is the same as setting κo = 0. Prior Variance Enors

Since the sum of the ephemeris, tropospheric delay, WRS clock and ionospheric delay enors are treated as a bias over die 6-second interval, the only prior variance left is that of the measurement enor. Since this variance is based upon measured values, its uncertainty is quite small. Thus, it is appropriate to use a large vo. Spatial Deconelation Enors

The variance of d e spatial deconelation enor between WRS k and Usery due to ephemeris enor can only be estimated in die prior statistics domain since ere are no measurements at the Usery location. This variance is

je£l_k (lJδX, - βδX_.)² = max[(l, -1,)^TC_X,(1, -1,,) 102)

where the parameters are as defined above for Equation 94. Ωk is a set of user locations in a region assigned to WRS k. The spatial deconelation enors due to the ionosphere are assigned to the GIVEs. Latency Enor Variance Since there is no latency in determining prior UDRE variance, the latency term need not be applied until after UDRE is determined. However, since it is the same for all WRSs and user locations, it is applied to the bound equation presented later. The worst case latency is

α_?.(Δt)²

ZPRla,e cy.t = ^ - 103)

where a$A i taken to be 19 mm/sec² , as defined in Minimum Operational Performance Standard for Global Positioning System/Wide Area Augmentation System Airborne Equipment. RTCA Document No. RTCA/DO-229, January 16, 1996, the content of which is incorporated herein by reference, and Δt is the time between the fast conection determination and the end-time of applicability of the fast conection. Prior UDRE Variance In this case, the prior UDRE variance is simply the measurement noise variance. That is,

σ___Λ£Λ = σ² _,rt 104)

Posterior Determination of UDRE

Given the prior UDRE variance with an unknown mean, die posterior update of die mean and variance is determined using actual measured data from die WRSs. For each 1 -second measurement from all WRSs k in the GEO footprint observing satellite i, enor samples for each WRS are simply the individual quantity of Equation 69 of d e fast conection description minus the weighted average over all WRSs in die footprint used to compute the fast conection, both evaluated at the same time using d e user time projections specified in die WAAS Signal Specification , which is provided in die aforementioned Minimum Operational Performance Standard for Global Positioning System/Wide Area Augmentation System Airborne Equipment. That is,

_,„(') = APR_ik{t)- _ai,APR_u{t) 105)

Determination of UDRE Bound The UDRE Bound for 6 enor samples, a prior κo = 0, a prior v_o of 20 and equal quantiles of 0.9995 for the mean and standard deviation is,

B_i = maxf|μ_| + ^3.5σ² _UDRE^ + 838S² _λ + (3.29)²σ²,_,,_t + δPR² _atencyΛ ) 106)

The value of T* = 6.6 is carried into the radical to be combined witii the prior and enor sample variances. Since the enor samples are at the worst case user location and because of the small sample set, the sample mean and the prior variance will play a significant role in the determination of _? . The spatial deconelation effect is expected to be small. The enor sample variance is likely to be small because the only enor sample variation would be due to measurement noise.

Equation 106 provides a quite pessimistic UDRE bound because of the limited number of enor samples and d e fact tiiat measurement noise from a WRS tiiat views the satellite at a low elevation angle, which is not a signal-in-space enor, dominates tiie bound. Other Approaches for Determination of UDRE Bound

Two other approaches diat are not as conservative and rely less on die WRSs with low elevation angles to die satellite will now be described. The first is one that does rely on previous sample means so that the mean enor is not completely unknown. The second takes advantage of weighting over the WRSs. Using Previous Sample Means

In this case, an exponentially decaying average of past sample means is used, where

Δt, Λ update μ_ = i - μ +^-μ* 107)

T_c j where T_c is the time constant of the decaying average, At_Upd_at_e is the time between updates of the sample mean, and where

μ_ = μ_ 08) is the new prior mean after using die previous one in determining the bound. The conesponding variance of tiiis exponentially averaged sample mean is

Initially, starting with a new WRS,

μ_ = μ,* HO) and

⁾

Widi Ko = 1, n = At_Upd_ate ^~ 6 and v₀ = 20, now T* - 6.16, with quantiles of 0.99965 and 0.99935 respectively for the sample mean and variance. Using Equation 85 with the prior mean and variance defined from the last exponential average update, but adding in the spatial deconelation and latency effects,

This bound not only has the slight advantage that T* is smaller, but more of an advantage because die prior variance is significantly less. Note tiiat there is still heavy reliance on die posterior mean, which will still respond to anomalies. The disadvantage is still d e fact diat a WRS widi a low elevation angle to d e satellite still has dominance, but to a much lesser degree. Taking Advantage of Weighting Over the WRSs

Weighting over WRSs, as defined above in the equations of the fast conection description, cannot be used in die previous two approaches for determining die bound because the worst case location could be easily de-weighted to not have an effect on d e bound. However, if the bound for each WRS is first found using the first approach, but not including the spatial deconelation and latency effects, and then weighting the bounds defined as

B* = |μ_| + P^ _k + 8 8S², ; k = 1,...«,- 113)

for bounds widi no prior information, or

B_Λ /29.18σ²; + 73S², +1.25(μ_li - μ- )² 114)

using exponentially decaying average of prior sample means such as is used in Equation 112. These individual bounds are then averaged over WRSs to arrive at the weighted bound

where σ² is the variance of the residual weighting enor determined from prior covariance matrices to

account for enors in Equations 94 and 99 (at the WRSs) not included in d e weighted bound. That is,

^σ .. + ^σ _β,ω)Λ_?.ι* ) | ^{i 6})

Note that if there is equal weighting (all aik = 1), Equation 116 equates to zero. Actually, certain terms of this Equation could be negative for over- weighted WRSs. However, that will only counteract over-weighted bounds. The key here is diat all the bounds include anomalies in the satellite and the use of Equation 115 brings back in die effect of die ephemeris and ionosphere at d e worst case WRS location. Note that die spatial deconelation between die WRSs and the user locations is still accounted for in Equation 115.

As die foregoing demonstrates, die WAAS WMS UDRE Determination Process of the invention uses Bayesian statistical methods to optimally blend a priori and a posteriori information of enors in the wide area differential conections to arrive at UDREs that are not overly-pessimistic because of an inadequate number of measured enor samples. Various variations to the UDRE Determination Process as well as the availability of "tuning" parameters allows for optimizing the technique as desired. GIVE determination (318)

Along with the Ionospheric Grid Point delay, an enor bound must be provided that bounds the conection by 99.9 %. This enor bound is known as the Grid Ionospheric Vertical Enor ("GIVE"). Due to the sampling rate of the ionosphere estimation process, it is very difficult to obtain enough samples to compute a 99.9% enor bound in a timely fashion. The process for determining GIVEs is shown in Fig. 19. Grid Points & Grid Cells Fig. 14 shows a sample grid for CONUS showing 38 cells sunounded by grid points. Each of a user's or WRS's pierce points lies witiiin one of these cells. The vertical delay at one of these pierce points is interpolated from die four sunounding grid point vertical delays. Thus, it follows tiiat die Grid Point Ionospheric Vertical Delay Enors (GIVEs) should be based upon prior and posterior statistics from within the cell tiiat is sunounded by die grid points. The so-called ionospheric delay truth is measured using semi-codeless L2 pseudoranges differenced with LI pseudoranges, conected for estimated differenced-measurement biases. The L2 pseudoranges can be quite noisy at low elevation angles. Because of die nature of these measurements, it is appropriate to de-weight the measurements. Since die evaluation of GIVE is in the vertical delay domain, it is also appropriate to also convert the measurements to vertical delay measurements with die obliquity factor Fy as well.

Relationship Between GIVE and UTVE

The GIVEs are defined to provide a bound on the User Ionospheric Vertical Enor (UIVE) that could occur anywhere in a cell sunounded by four GIVEs when interpolated from the grid point vertical delays. UTVE is itself interpolated from die GIVEs. For convenience and simplicity, define, for cell m,

UIVE_m = ∑W_n -GIVE_n,m 117) π=l where the GIVEs sunound die cell and die W_n are as defined in Paragraph A.4.4.7.2 of the aforementioned Minimum Operational Performance Standard for Global Positioning Svstem/Wide Area Augmentation System Airborne Equipment. Common GIVEs can be associated with more than one cell, so they are not usually defined exclusively for a given cell. Note, also, that there are many more grid points than there are cells. Thus, enors observed in a given cell must be distributed to more than one GIVE. The metiiod of doing this may vary as long as Equation 117 provides die 99.9% bound of the enors within the cell. The deterniination of the GIVEs from die UIVEs is die last step in Fig. 19. Essentially equation 118 is solved for GIVE as the last step in determining die GIVEs for each IGP.

< [Hani mGlVEo 118)

where:

UIVEiMii) user ionospheric vertical enors that could occur anywhere in M cells GIVE win grid ionospheric vertical enors which provides a bound on the-UTVEs

[/_ ]_(**> matrix of weighting factors(PT,'_ )for the IGPs of the sunounding cells containing

the user location for the centroid pierce point position for each cell. M number of cells

N number of grid points A description of the complete processing shown in Fig. 19 is provided below.

Compute Centroid Pierce Point Position A weighted centroid pierce point is determined based on die measured pierce point in d e cell. For this associated pierce point, the UWE_m is determined at tiiis point for the cell. The method of determining the centroid pierce point is:

^XPP.m = ~ Σ ^X _PP,J 1 9) ⁿm ιx =l and

where (xpp j^ppjj) is the pierce point between WKEj and satellite /and the by are the weights defined in die next section. After evaluating Equations 119 and 120, Xpp_tm ^aΑ^ypp,m must be limited to the interval

[0, 1] before using them to define die Wfc This procedure is repeated for each cell. The procedure works its way around until all cells are covered that have pierce points in die cells. In cases where no centroid pierce point can be computed, die GIVEs are set to not monitored if adjacent cells cannot provide die necessary information to compute a GIVE value..

Compute Weighted Measurement Enors

At times, d e measurement enor can shadow the actual estimation enor. In the computations performed in die WAAS WMS to estimate die grid point delays, methods are used diat always de-weight measurements based upon botii prior and measured variances of die measurement enors.

Thus, it is not preferable to also weight die post-conection enors with the inverse of the measurement enor variances in deterrnining die true post-conection enors.

Using weighting, the n_m statistical enor samples in Cell m become normalized, where x_u' = b_lΛ ; ixj= \, ..., n_m 120)

where the weights are defined, for measurement ij enor standard deviation σ_ηy in Cell m at die pierce point for satellite i and WREj

which have the property

The measurement enor statistics of the new weighted samples are all equally

σ % 124)

Σ k=\ and

= 125)

where

where die subscript e indicates the estimation enor component diat comes from die grid point estimation processes. In general, over d e long haul, both of the mean enors should be zero.

The new posterior enor samples are now normalized. However, in this case, the estimation enor variances for the prior statistics are not weighted since tiiey are extrapolated from the grid point estimation enor covariance matrix that is already in d e vertical domain. Compute Prior Variance for UIVE for Cell m

The prior variances for the user ionospheric vertical enors are defined at die WRE/satellite pierce points diat occur in every observable cell m. For a given cell m, die maximum prior variance of tiiose pierce points is given as

127)

where the estimation enor covariance parameters CGP_ΠI._IJ is the enor covariance between grid points n and / for those grid points sunounding pierce point ij ( if /.=/, it is the variance). This information is

derived from the covariance matrix (P) from the Ionosphere Estimation Kalman filter. Compute Posterior Variance and Sample Mean for UTVE for Cell m

The posterior sample mean and variance for the user vertical ionospheric enor in observable Cell m are computed from n_m weighted measurements at the WRE/satellite pierce points that occur in that

cell per 10-second interval. For a given Cell m, these weighted measurement enors are given as

tono.ij ij ^₁-,_2,,, - ∑ ^δτ Gfn,!j ixj = \, ..., n_m 128)

^( -y) for vertical grid point estimates c τQp_n>ij sunounding the cell and conected dual frequency

pseudorange differences PRL\-L₂,ij collected since the last update, where

y = /..//_₂ = (77/60)² 129)

Equation 128 applies to determining the sample mean and sample variance for n_m measurements for the

UIVE for Cell m, which then would be used for determining GINEs to be broadcast using measurements

collected over the update interval. The cδτQp_n jaxe those estimated just prior to generating the data for

the broadcast. For a maximum 5-minute update interval, there can be a significant number of samples within a cell. Compute UINE Upper Bound for Cell m

Using Equation 90 as an example, the UIVE Bound for 30 samples in Cell m would be

UIVE_m = |0.882μ + 4.94j0.286_CT ² _VIono,m + 0.690S² _m + 0.084μ² 130)

Since there are more samples here than for the UDRE at the worst case location, the sample mean and die sample variance will play a more significant role in die determination of ULVE_m than will die prior variance of the worst pierce point in that cell. The GIVEs sunounding d e cell are then determined using the procedure described above. Compute Spatial Deconelation Enor

The GIVEs determined above are based upon measurements and prior variances at pierce points in a given cell. Since these measurements do not necessarily occur at or near a grid point, spatial deconelation effects must be added to the GIVEs. Rather than basing these effects upon the actual measurements, it is more appropriate to base them on the vertical delay estimates at those pierce points sunounding the grid point n given the processed measurements. Defining a spatial deconelation coefficient as

where «_pn « is the number of pierce points (within the update interval) in the cells adjacent to the grid point (1, 2, 3 or 4, depending upon where the grid point is), and djk is the distance between they^'&th pair of pierce points. Then defining the spatial deconelation enor at grid point n as

sd_GPn = Vsd_GPn x ^■ Σdkn 132)

"pp.n *=1 which is the spatial deconelation coefficient of Equation 131 times the average distance between the pierce points sunounding the grid point and the grid point. The value of Equation 132 is then added to the GIVE for grid point n before it is quantized into die GIVE indicator (GIVEI) of Table A-9 of the aforementioned Minimum Operational Performance Standard for Global Positioning Svstem/Wide Area Augmentation System Airborne Equipment. Determining UTVE Alarms

The posterior statistics are also updated more often to determine if there are any changes to a UTVE derived from cunently broadcast data that would wanant an alarm condition. In this case, " )CB 133)

are weighted UTVE enors determined for the cunently broadcast data using the same measurements used in Equation 128 but with different weights b j. The weights for determining the future GIVEs are based

upon measurement variances of all measurements during an update interval. Here, the weights are based upon measurement variances of measurements collected during the cuπent broadcast interval. The magnitude of the weighted enors are compared to

UIVE_m = ∑W„ -(GIVE_n,m)_cB 134)

where GIVE_nιtn are the cunently broadcast GIVEs at the grid points sunounding the applicable Cell m. If the magnitude of the weighted enor exceeds the UIVE_m twice in a row, the GIVEI for one of the grid points is increased to the next level, starting with the grid point closest to applicable pierce point, then the next closest, etc., until the UIVE_m is no longer exceeded. This constitutes an alarm condition for the applicable grid points, for which the GIVEI can only be lowered with a new GIVE determination of the type described above. Timing of Fast Conection. UDRE. and GIVE Computations The integrity time to alarm limit is 5.2 seconds for the WAAS system. The integrity time is defined as the instance that the integrity parameter (for example, satellite UDREs with associated fast conections or GIVEs with associated IGP delays) becomes out of tolerance to the time the user receives the broadcast message warning of the hazard. Also, since the GEO uplink message is timed to the WAAS Network Time (WNT) 1 second epoch, the uplink messages must by synchronized to the 1 second WNT epoch. The master station must ensure that both the processes to determine the conections and the timing of when the messages are sent to the GUS are timed properly. The WMS solution to the timing problems described above is to incorporate a GPS receiver at the master station as a timing source. This source is used to define the 1 second C/A code epoch at the master station. By using an interrupt from a GPS receiver at the WMS, the WMS is synchronized with the GUSs and ultimately GPS time to within 50 nanoseconds. This synchronization allows the WMS software (including the independent verification functions) to meet the WMS software allocated 0.5 second time to alarm budget which is the critical timing budget parameter in the WAAS system. This is accomplished by pre-processing most of the estimation and conection data. In this case pre-processing means that the data is processed as it arrives. The critical pseudorange conection UDRE, and GIVE processing described above is performed based on synchronization on the 1 second interrupt using die buffered results of the pre-processing functions. This data is then forwarded to the safety monitor for the independent verification checks and generation of the final conection message. . This allows the data to be computed, verified ,and validated in the safety monitor in die 1 second epoch just preceding the uplink epoch thus ensuring the that the time to alarm is met as well as the proper timing to the GUS for transmission in the next uplink epoch. The fact that the safety monitor is completely data driven based on the 1 second interrupt in the conection and verification processors also simplifies the safety monitor design. Clock Steering (320)

The Conection and Verification process provides clock steering for both the WAAS Network Time (WNT) and the GEO clocks. This process is shown in Fig.21. By specification, the WNT must be within 50 nanoseconds of GPS time where WNT is the ensemble average of the WRE receiver clocks. The apparent clock offset and drift of the GEO is located in the navigation message of the WAAS. Under certain conditions, this value could become excessively large, due to a slow accretion of drift over a long time. If this clock offset becomes excessively large, it will not fit within the message structure of the navigation message, thereby denying the users the ability to use the GEOs as a ranging source. To avoid this problem, the WMS sends a rate at which to conect the GEO clocks to the GUSs. The WNT Clock Steering process ensures that WNT time is within 50 nanoseconds of GPS time by calculating a time steering conection command. In general, the ensemble average of all GPS satellite clocks, conected with their respective navigation message data, is defined as GPS time. Thus, the ensemble average of all WAAS GPS satellite clock conections is the difference between WNT and GPS time. Ensemble averaging the long term clock conections over all GPS satellites is sufficient to determine that time difference. That is, the average time difference over a period of a day is:

1 ^NGPS 1 ^Nr, WNT_CPS = — ∑ — ∑δ ₀,( 135)

* * C_PS ι=l -" " TV m=l where NJ¹ is the number of slow conections broadcast over a period of a day for satellite i, and NβPS is

the number of GPS satellites.

Because both the GPS satellites and WNT are based on stable clocks, a daily update is sufficient. A time steering conection can then be computed as a WNT frequency shift for the next day equal to the quantity of Equation 135 divided by 86,400 seconds. This frequency shift is applied to the time solution of each WRS in the orbit determination process, which will automatically steer all clock conections toward GPS time as well as WNT (GEOs' time base), thus ensuring that WNT time is within 50 nanoseconds of GPS time.

The GEO clock steering process calculates a GEO clock steering command. To ensure that the GEO clock offset with respect to WNT is within the range of the type 9 GEO navigation message, the master station sends a rate (GEO clock steering command) at which to conect the GEO clock to the GUS.

GEO satellite clock steering is performed once each 24 hours unless required more frequently to remain within the range of the type 9 GEO satellite Navigation Message clock conection parameters.

The equation to determine the GEO clock steering command parameters is shown below. GEO satellite clock states (offset and drift) are available from the orbit determination filter every 5 minutes. The orbit determination process provides these as clock conections for the GEO satellite navigation type 9 message. The clock offset is propagated forward in time using the broadcast clock offset and drift conections. The clock offset must be maintained within the range of the parameters to be broadcast in the GEO type 9 navigation message. These ranges are: clock offset: 9.536743E-07 seconds clock drift: 1.16415E-10 seconds per second. If the propagated clock offset and drift exceed the limits specified by the following interval, a clock steering rate command for GEO satellite i is issued to the GUS:

(9.89996E - 07 seconds - <S_Iim,,_,C£0) > Δt„_o, > -(9.89996E - 07 seconds - Jt_Iim,G£0) 136)

where _<5t_limC£0 configuration limit to provide a margin for uncertainty and an adequate lead time to

minimize the magnitude of the required steering rate. The configuration limit is initially set a nominal value of 4.89996E-07 seconds.

The value of 9.89996E-07 seconds represents the range of the GEO satellite navigation clock offset and drift propagated forward in time 300 seconds beyond the broadcast time of applicability. The propagation time conesponds to the 300 second OD update for the GEO satellite clock offset and drift information. The equation for propagating the GEO satellite clock offset forward in time is:

Δt_ao . = δa_Cf0J(t_0J) + δ _Cfl(t_QJ)*At_i 137)

where δa_G/0J(t_0J) represents the broadcast clock offset for GEO satellite i with time of applicability tQ

and -«_c/1(t_o ) represents the broadcast clock drift for GEO satellite / with time of applicability tg . The propagation interval is calculated by:

t; ⁼ 'o+300secj ^— *0j 13o) where tø is the time of computation by orbit determination filter.

The clock drift rate steering conection is calculated by:

δ d„o_rc_tk, . = — ^ ₃₀₀ 2ι_^. seconds p ver second 139) h,

where: Δt . is the quantity from equation 137 calculated at the previous 5-minute SOD update epoch.

The C&V processor computes and sends the GEO satellite clock steering command to the controlling GUS. Upon receipt of the response to the steering command from the GUS, the C&V then sends a command to the orbit determination filter to update the GEO clock states. The orbit determination filter is notified of the steering so that it remains synchronized with the new ranging data that will be received from the reference stations after the controlling GUS starts to steer the GEO clock. Safety Monitor (50)

The software contained in the present implementation of WAAS is a mixture of software specifically developed for the system and commercial off-the-shelf/non-developmental software. The commercial software includes large operating systems, network software, and router software as well as embedded firmware used in the various receivers, clocks and signal generators. Typically, such commercial software is not developed using the rigorous standards used for safety critical avionics software and, therefore, does not take into consideration the effects of the severity of the hazard associated with the software's failure. The FAA has established a series of hazard levels according to severity of the effect of software failures and has established criteria for categorizing the impact of a particular failure against these levels. The four levels of failure, in descending order of severity, are catastrophic, hazardous, major, and minor. The FAA established that in a WAAS system, any software must be certified to the level commensurate with its failure as prescribed in RTCA DO-178B and must be developed also as prescribed in DO-178B. DO-178B defines levels of software certification conesponding to the hazard levels.

The two highest level failure scenarios for the WAAS as anticipated by the FAA are any failure or combination of failures that lead to a loss of continuity of the WAAS signal-in- space or that would result in the transmission of hazardously misleading information. Hazardously misleading information occurs when conections enor is greater than the outer tunnel boundary and when the UDRE/GIVE functions do not provide the conect upper bound to the enor. This condition can potentially send an aircraft beyond the outer tunnel boundary and is considered a hazardous condition for both en route/non-precision approach and precision approach.

The continuity of a function consists of continuity of navigation (providing the required accuracy) and continuity of integrity. Continuity of integrity can be provided by broadcasting integrity information from the ground or by the availability of Fault Detection and Exclusion (FDE) in the airborne system. A loss of navigation accuracy or a loss of the integrity function results in a loss of continuity of function. A loss of either function precludes the use of the system as a primary means of navigation.

The prefened solution for the present WAAS system for maintaining the continuity of function when a loss of broadcast integrity information occurs is to provide a continuous ranging signal to the aircraft so that FDE is available to the aircraft to provide the integrity function. The method of assuring that a continuous ranging signal is provided to the participating aircraft is described below with respect to the GEO uplink subsystem.

The prefened solution for assuring that hazardously misleading information is not transmitted is to provide a safety monitor 50 which verifies the conectness of the information. The safety monitor 50 contains safety-critical software which is certified to the appropriate DO- 178B level.

The capabilities of the prefened WAAS system are implemented such that critical integrity (i.e., probability of transmitting hazardously misleading information) and continuity requirements are met through the use of appropriately certified software. The remaining capabilities that are not safety critical may be accomplished through the use of commercial software wherever possible.

The WAAS system of the present invention meets both the integrity and continuity requirements. Integrity is maintained through the use of safety monitor 50 software that guarantees that hazardously misleading information is not broadcast to users. Continuity of navigation for en-route/non-precision approach is met by guaranteeing either broadcast integrity information or continuous broadcast of GEO ranging and associated navigation messages in the event of a system failure.

As shown in Fig. 20, WAAS safety monitor 50 performs three main functions. First, safety monitor 50 independently verifies the conection data (501) calculated by the conection and verification processors 40 and 42. Second, safety monitor 50 determines if the pseudorange enor, after applying conections, is bounded by the appropriate combination of UDRE and GIVE (502). Third, safety monitor 50 performs message generation (503) Each WRE forwards to safety monitor 50 the local differential conections and LI/, which the GPS receiver at the WRE at wraps in a 24-bit Cyclic Redundancy Check (CRC). This allows safety monitor 50 to authenticate that the information has not been altered by the commercial operating system executing on conection and verification processors 40 and 42 or by the communications software on the network. The information includes a time-tag that enables a WAAS Message Processor at each GUS to determine if the information it is receiving is timely. A CRC is appended to the data by safety monitor 50 before its transmission to the WAAS Message Processor. This allows the broadcast information to be validated before transmission to the satellite. This is to ensure that the data has not been altered as it is passed through the verification processor 40, and was transmitted through the routers and over the network before finally being received by the WAAS Message Processor. Safety monitor 50 also performs automatic integrity checks to ensure that the WAAS is meeting time-to-alarm and message broadcast rates.

Fourth, safety monitor 50 provides the capability to verify conections from two independent processing threads (the conection and verification processors), to check orbit determination calculations, to generate and verify broadcast messages, and to validate all broadcast conections against truth data.

The validation of UDRE and GIVE is performed using enors in position based on the differential conections computed by the WRS receivers. These differential conections are computed from the pseudorange and geometric range computations made in the receiver after all WAAS broadcast conections have been applied. By applying the cunently broadcast WAAS conections, the computed position enor is equivalent to the total enor of the WAAS broadcast conections. The total enor of the prefened WAAS is bounded by the combination of UDRE and GIVE. Safety monitor 50 performs a process to validate that the computed UDRE and GIVE values (before broadcast) do indeed bound the WAAS system enor computed from the differential conections.

The first step in the validation process is to compute the user position enor based on the user's weighted least squares solution, a process described in the aforementioned Minimum Operational Performance Standard for Global Positioning System/Wide Area Augmentation System Airborne Equipment. The range residuals or (differential conections) are computed in the WRE receiver as:

= ^ro s,, - ^rpred,i(^Xs_Ur_Vey) 140)

where r₀ ~ ,^■ is conected for i) ionosphere using the broadcast ionosphere conections based on

the WAAS IGP model, ii) troposphere using local meteorological measurements, iii) SV clock using broadcast conections and iv) other normal pseudorange conections (Earth's rotation, etc.); and ^rpred_,i i^s computed using surveyed position and WAAS conected satellite position (using

broadcast data).

The position and clock enor is computed using a weighted least squares process:

141)

where the row for a satellite of the G matrix is defined as:

g, = [cos 4 cos E_t sin A, cos £, sin _T, l] 142) and the G matrix includes both the position and clock terms as defined below.

Δr = G[Δx^τ cA*J^T = [G_to l][Δx^τ cΔt_wra]^T

143) and the weighting matrix W weights n SV measurements is defined as:

and the 99.9% weighting terms preferably are:

δr₀ ² ₉₉₉, = UDRE + [3.29c ^(E,)]² + f{F²UIVE²) 145)

The 99.9% covariance matrix is formed as:

C, = (G^TWG)^"' 146)

Safety monitor 50 can validate the position Vertical (VUNE) and Horizontal

(HUNE) navigation enor using the data from the 99.9% covariance matrix as follows:

where the delta position are predefined system limits.

The next step in the UDRE and GIVE validation process is to convert the total position enors back into the range domain. This is accomplished by removing the clock enors from the position solution and rotating the conesponding position enor into the range domain as shown below: Δr' = G^Δx = G^A.Δr 148)

This step removes the clock and other common enors and includes both UDRE and UTVE enors in the conesponding residual.

The UTVE residuals, which are used for validation of the computed GIVE are computed next. The process uses the same gains (G) used to determine position enor. The ionospheric conection model range residuals are computed in the receiver as:

fe, = ^~ P^RL ~ ^Ctesv,) - "model, 149)

These residuals are then used to compute position enors in the same manner as defined for the differential conections. These position enors are then also rotated into the range domain in a similar fashion thus removing common bias components: ^Δr _f = G_AlAx_uιyE = G_Δ. A.Δr_vn,_£ 150)

The computed GIVEs are then compared to the range residuals as shown. If the inequality does not hold, the associated GIVEs for the UTVE are increased until the equation is true:

Ar U'IVE, l ≤ F. - UIVE. 151)

where the UIVEj is computed from the computed GIVE values for the IGPs sunounding the

satellite's pierce point.

The ionospheric residuals are then subtracted from the total range enor residuals to form the residuals for UDRE comparison: Δr " ,„ = Δr ' - Δr ^ = G_Δl A_Δl (Δr - Δr^_£ ) 152)

These are then compared to the computed UDRE values. The computed UDRE values must be less then the UDRE range residuals computed from the receiver differential conections. If this inequality (Equation 153) does not hold true, the computed UDRE values are increased until the following inequality holds true:

Any SV measurement with significant enor will add enor to the clock solution. Therefore, if a large single enor is sensed, that satellite is removed from the solution and the process of validation is re-computed without the satellite. That particular satellite is then flagged as "don't use."

In order to accomplish the UDRE and GIVE validation function, the WAAS receiver data is fed directly into safety monitor 50. The integrity of this data is maintained, both at WRSs 2 and across the communication links. The latter is achieved through enor control algorithms such as a CRC implemented at the transmit and receive ends of the link. This software must be certified to RTCA DO-178B Level B.

Safety monitor 50 also ensures that the messages broadcast to the satellites themselves are free of enors. This is achieved by performing message generation within safety monitor software. The pseudorange check is also used to provide an independent check of "don't use" messages, which helps to ensure continuity of navigation.

From an integrity standpoint, the WRSs data input to the conections processing is also input to safety monitor 50 and certified enor control is incorporated throughout the communication network. In order to accomplish this, the Local Differential Conections (LDC) input data for the process comes from a verified source, such as a GPS receiver. The GPS receiver generates a message with a 24-bit CRC that is passed through the network with additional enor control appended by the communications process. The additional enor control is stripped away before it is passed to safety monitor 50 but the 24-bit CRC is decoded in safety monitor 50. Therefore, the integrity of the data is guaranteed as it passes through the communication paths of the software with a lower level of criticality. The receiver itself must also be certified to RTCA DO-178B Level B since it is the source of raw data.

Safety monitor 50 performs additional functions such as a time-to-alarm and monitoring check. The check is initiated by sending pseudorange data to the conection and verification processes for a satellite that does not exist in the pseudorandom noise (PRN) mask. Both processors generate a fast conection message for the satellite which is compared in safety monitor 50 and a fast conection message is then generated which is sent in the Type-2 or Type- 24 message. Safety monitor 50 then monitors the system timing from the start of the broadcast of the message to the time of receipt of the message back at safety monitor 50. This automatic integrity test also tests the conection processing functions in the conection and verification processors.

Safety monitor 50 validates the broadcast messages in two ways. First, it stores the broadcast messages and verifies the proper transmission by checking the sent messages against the messages received at the WREs. Second, all messages have certain broadcast rates that must be maintained, i.e., fast conections every six seconds. These message rates are verified for all WAAS broadcast messages by safety monitor 50.

Safety monitor 50 must reside in either the GEO uplink system or the WMSs because safety monitor 50 uses data from both reference stations 2 and WMSs 8. However, communication bandwidth considerations make it preferable that safety momtor 50 be located in WMSs 8. The addition of safety monitor 50 at each WMS validates the outputs from the conection and verification processors 40 and 42 and satisfies the hazardously misleading information safety requirement. The safety momtor function is implemented on a small platform, for example, an Intel 80486, with a DO-178B Level B operating system or a safety critical kernel. Keeping the size small and the functionality limited are key factors to assuring the ability to limit the types and quantity of failure modes. The incorporation of safety monitor 50 certified to

RTCA DO-178B Level B, consistent with the identified hazard condition category for integrity, enables a lower certification level in the complex conections processing software. This software, including the real-time orbit determination package RTOD™, is certified to RTCA DO-178B Level C. GROUND UPLINK STATIONS (GUSs)

Each GUS 10 performs the functions of broadcast and ranging. As shown in Fig. 22, each GUS 10 receives a 250-bit formatted WAAS message once each second from each master station 8 in the WAAS system. To improve the availability of WAAS messages, each GUS 10 is connected to two WAAS backbone nodes on the TCS, which is the network over which data is transfened between the WRSs, the WMSs, and the GUSs, and preferably consists of the FAA LINCS network. GUS 10 selects one WMS as its message source and encodes the received message using a 1/2 rate forward enor conecting convolutional code. The resulting 500-symbol message is modulated on a GPS-type signal and uplinked to GEO satellite 6. GEO satellites 6 are preferably the INMARSAT satellites. Each GEO satellite 6 is served by two GUSs: one operating as the primary uplink and the other operating as a hot standby. The two GUSs serving a GEO satellite are operationally independent and located at geographically diverse Ground-Earth Stations ("GES") separated by a minimum of 300 miles. A GES is a facility consisting of one or more GUSs and provides shelter, power, and operations and maintenance services for the GUS. The GEO satellite "bent-pipe" transponder shifts the frequency of the signal and broadcasts it to the users. Transition between primary and backup GUSs is initiated, when necessary, to maintain the availability of the WAAS Signal in Space.

Each GUS 10 consists of a signal generator subsystem 250 and an RF uplink 252. As shown in Fig. 23, signal generator subsystem 250 includes an atomic clock 254 (preferably model FTS 4040 A, available from Frequency & Time Systems), a GUS Processor 256 (with its associated keyboard and monitor), GUS communication routers 258, a WAAS Receiver Subsystem 260 (preferably NovAtel WAAS Receiver 044645-0001), a GPS antenna subsystem 262, and a WAAS generator 264 (preferably model 7201, from Stanford Telecom). The signal generator subsystem selects a WMS to receive the WAAS message, encodes (FEC) the message, modulates the message on an uplink signal using C/A type code, and controls the timing and frequency of the message transmission for use as a ranging signaL RF uplink 252 is responsible for converting the IF signal received from signal generator subsystem 250 and broadcasting it to GEO satellite 6. RF uplink 252 also receives signals from GEO satellite 6 for use by signal generator subsystem 250 and GUS Processor 256. Additionally, RF uplink 252 includes a frequency distributor amplifier that provides a frequency standard back to signal generator subsystem 250.

Each GUS 10 includes a signal control process to control die uplink signal to die GEO satellite 6. The GUS signal control process controls the uplink in such a way as to make the signal relayed through the GEO satellite 6 mimic the signal emanated from GPS satellites. As discussed below, the process for die Backup GUS is slightly different from the process used for d e Primary GUS. Fundamentally two parameters of the signal transmitted to die GEOs 10 are controlled. The first is die rate of change in die pseudorandom noise (PRN) code rate, which is measured by the WAAS receivers as the pseudorange. Optimally, the controlled uplink code rate should equal the measured downlink pseudorange rate with a conection for ionosphere and steering of the clock by a master station. The second parameter controlled is the carrier frequency such that it remains in a 1540 ratio with the code rate with a conection for the ionosphere.

The GUS signal control process, as shown in Fig.25 uses a Kalman Filter to estimate the enors in the cuπent signal. These estimates are used by classic first and second order control laws to compute the parameters that die signal generator uses to time its signal. The process uses L1 L2 measurements from the WAAS Receiver 260 of die GUS to determine the downlink C/A code phase rate enor and the incoherence between the downlink C/A code phase and die downlink carrier phase. The rationale for using a Kalman filter to generate estimates of die signal enors is that the Kalman filter is a statistically optimal estimator that can run in real-time and does not require storage of past data as batch least squares solutions typically require. Thus, for the real-time operation of GUS 10, it is an excellent technique. Classic control laws are used instead of merely applying the outputs of the Kalman filter to die signal generator in order to produce a smooth transition in removing the enors that die Kalman estimates. This architecture is unique and inventive in terms of solving this problem.

As is known in die art, a Kalman filter is a template that requires the matrices to be rationally filled in order to work for an application. The various matrices and measurement calculations are original and designed in order to generate enor estimates that are easy to convert to signal generator controls.

The inputs to the GUS signal control process not only include the code and carrier phase measurements of WAAS receiver 260, but certain inputs from master station 10 and the primary GUS 10 in the case of the backup GUS 10. During primary closed-loop operations, frequency conections from the master station must also be accepted in order to control the GUS/GEO clock offset and drift. The inputs to this processing are described in die following table, which also includes inputs for die Backup GUS, designated widi a *:

Input Source Description Units

PR received, LI (PR_Ll ) GUS Receiver Received C/A code pseudorange at the LI meters frequency PR received, C_d0« (via L2)GUS Receiver Received C/A code pseudorange at the C meters downlink frequency (Cd_0Wn via the L2 receiver channel) φ received at the LI GUS Receiver Received C/A carrier phase at the LI LI cycles frequency (φ_u ) frequency

O_Li GUS Receiver Standard Deviation of the Received LI C/A meters code pseudorange

CTCdown GUS Receiver Standard Deviation of the Received L2 meters

(C_do ι_>)C/A code pseudorange

PR LI wide band GUS Receiver Wide band LI pseudorange correction to meters correction generate independent 1 second measurements

"Uplink Pseudorange Primary GUS The 4 signal generator control values used by meters, m/sec phase and phase rates, the

Uplink Carrier phase primary GUS. and phase rates

"Primary Uplink and Primary GUS Primary Uplink and Downlink smoothed meters

Downlink smoothed ionospheric delay ionospheric delay

*ECEF x, y, z and WMS (Message Earth-Centered-Earth Fixed Position, velocity meters, m/sec, derivatives of GEO and to Type 9 and acceleration of GEO and time of m/sec² and parameters) applicability seconds of week

*ECEF x,y,z coordinates WMS Earth-Centered-Earth Fixed Position of both meters of both GUSs GUS antennas

GEO frequency correction WMS GEO/ GUS frequency correction meters/ sec (cAfcus)

Table 3

The GUS signal control process (in the primary GUS) uses the inputs described in Table 3 to estimate the code phase, code phase rate, carrier phase and carrier phase rate for the signal generator. Measurements are formed for an input to a four (4) state Kalman Filter to estimate these four quantities. After this estimation process, a control process, described later, is used to derive commands to the signal generator for die uplink signal. A block diagram of a Kalman Filter is shown in Fig. 26. The cycle is performed once every second after the measurements are received from me receiver and in die case of the backup, the primary GUS 10 and master station 8.

The states of this Kalman Filter state vector x at time t are as follows:

1) δPR - uplink/downlink pseudorange difference, conected for the ionosphere and GUS time

offset conection in meters =x\_tk',

2) uplink/downlink pseudorange rate difference, conected for die ionosphere and GUS dt time drift conection in meters/second = Λ_ ;

< 3) PR - smoothed difference between received pseudorange and carrier phase in

meters, conected for die ionosphere = x^k, and

4) - difference between received pseudorange rate and carrier frequency in meters

per second, conected for die ionosphere = x^k.

State 2 is the estimate of the Doppler enor to be applied to the control process for the command to die signal generator to equate uplink and downlink Doppler. The measurement for tiiis state is conected for GUS time rate conection and die measured downlink and uplink ionosphere delay rate-of- change. State 4 is the downlink code/carrier coherency enor, conected for downlink ionosphere delay rate- of-change. States 1 and 3 are in the conected pseudorange and carrier phase measurement domain, scaled appropriately. In effect, the Kalman Filter is deriving smootiied derivatives.

The measurements of measurement vector z at time tk for the filter is formed as follows:

^Z\.k ^{= E}δPR_k ^{= £} PR_k ~ ^£ _PR₀ ~ ^C^GVS.k

— — 154)

= PRn ._k - ^PR _U,k - &^t™°-C.4 +cΔtiono,Ll.„ - - „. ~cAt_{Gυs k}

in meters, where PJL _Λ = SF_C 155)

is the integral up to time tk of commanded at time t_n after scaling and

truncation required for die Signal Generator, cΔtGusΛ ⁼ cunent time conection for die GUS witii respect to WAAS Network Time ("WNT") in meters, which is an integration of GEO/GUS frequency conections received from master station 8, as incorporated in the Signal Generator commands (to cancel that accumulated in the Signal Generator as it is not an enor). This is to prevent the Kalman filter from tracking the WMS frequency conection as an enor. Equation 155 can be replaced with a measured output of the Signal Generator that supplies the same information, scaled appropriately. The ionospheric delay conection for LI is

At_Mk = -i-(/^>Λ_L - PRc^s) meters 156)

where y = /__. '/c . = °-lS831242 and c = 2.99792458^χl0⁸ meters/second. cΔtiono i : is filtered using a

low-pass digital filter with a time constant Tr (on the order of 60 seconds) of die form

CΔtiono.Ll,* = cΔtkmo,u,*-ι + — cAt_{ioao Ll k} 157)

initialized with Equation 156 using the first measurements. The C-band uplink conection is: cΔtiono,c_ = = 0.05956cΔti_Ono.Li meters 158)

Measurement 2 is

^z _2k = ^εs_Φk = ^ε _k ~ ^ε ₀ = ^PR .k ~ ,k ^~ 2cΔtio„o,_u - ε_φt meters 159)

at time t#. The initial formed measurements (εpR and ε_φ) are subtracted from all subsequent measurements

because of pseudorange biases and carrier cycle ambiguities. z₂£ is reset any time the LI phase

measurements from the receiver are reset because of a sensed cycle slip (Lock Time is reinitialized), in which case die first measurement after the reset is discarded. PRu is calculated in botii measurements as

the sum of the pseudorange measurement and d e wideband conection received from the receiver.

Except for the time offset conection, die absolute pseudorange or carrier phase differences are not important, only the rate-of-change. However, die control loop will have die effect of driving tiiese formed measurements towards zero. In fact, die changes in State 3 over periods of time are important performance parameters (code/carrier phase coherency).

The filter model includes matrices describing the estimation process for the states and measurements described above as well as die state transition matrix, die initial process noise matrix, the initial enor covariance matrix, the measurement matrix and d e measurement noise matrix. These matrices are defined later. The state transition matrix is:

1 Δt 0 0

0 1 0 0

Φ = 160)

0 0 1 Δt

0 0 0 1

where Δt is the time update interval between uncoπelated 1 Hz measurements or 1 second. The initial process noise matrix is

where aR_tk is the cunent line-of-sight acceleration to the GEO and ho, h._\, and Λ_-2 are Allan Variance parameters describing die stability of die GEO local oscillator, where

Λ₀ = 2.0x10 ,^'-22 162) h_Λ = 7.2x10 ^■23 163)

During orbit maneuvers, a ι.k is set to die appropriate value for an orbit maneuver, unless the maneuver is embedded in die GEO x, y, z and derivatives data as will be die usual situation. The initial enor covariance matrix is

'11.* 0 '12,* 0

0 δv\ 0

P' = 165)

'12,* 0 '22,* 0

0 0 0 c² χ l0^"12

where the r k are the elements of the R matrix described below, δv is the line-of-sight velocity uncertainty of the GEO, and die uncertainty of the carrier/code coherency is die frequency offset of the GEO's local oscillator, tentatively taken to be 1 part in 10⁶. There is no uncertainty in states 1 and 3 because their solutions are immaterial. The filter is simply estimating tiieir derivatives, which do have uncertainties. Also, the velocity of the GEO will be fairly well known, so that the imtial uncertainty in state 2 will also be quite small. The measurement matrix maps the states into the measurements and is

166)

representing that states 1 and 3 are measured wid d e formed (but noisy) conelated measurements, where the measurement noise covariance is represented as

where, neglecting the small carrier phase measurement noise,

^σ>, 168)

σ! = ^{■ •}{Wι-r)-ιf ^■ PR, + σ PR 169)

W-r c, _»> .

and

A2^ff, σ„, = fl[ι- r_/(ι-r)lr_/(ι-r)-ψ - σ 170) τ, (ι-r)² PR,

where σ². = 10σ² L,lm„eα„-.,*, meters 171)

and

^■ PR, <*«,J = 10σ ^{■ 2}C_Λmmeas,k m "^«e"t~er^«s 172)

where the <Jmeas,k ^are reported by the receiver along with d e measurements at time tk- The pseudorange measurements are provided as two numbers each (a nanow loop bandwidth measurement plus a wideband conection) which are added together. The resulting bandwidtii is 10 times the nanow loop bandwidti used by the receiver to compute the standard deviation of the measurement. Thus, the wideband conected measurement has 10 times the variance.

The Kalman Filter equations for updating die measurements are defined as follows, starting with an initial state vector estimate of

and its initial enor covariance matrix of Equation 165, die two formed measurements are applied one at a time as follows: The Kalman Filter gain vector at time t for measurement zjk is calculated as

A_t = HP_t'H^T +R, 174)

where Ak is the predicted covariance of die measurement residuals (or innovations) for measurement vector zk, H is the measurement matrix, and P₄ is the predicted enor covariance. The gain matrix is then

K_{i =} P;H^TA-' 175)

The state vector and enor covariance are then updated by incorporating measurement k as i₄ = (I - K₄H)X₄ + κ₄z₄ = χ'_k +κ₄(z₄ -Hi,) 176)

P₄ = (I-K₄H)P₄' 177)

where the term z_k - Hx_k' is the measurement residual vector (or innovation vector).

After processing all measurements at time tk, the state vector and enor covariance matrix are projected to die time of the next 1 -second measurements as x;₊₁ = Φx₄ 178)

P _+I = ΦP,Φ^T + Q, 179)

where Qk is adapted from die measurement residual vector as described below. Filter divergence may be caused by model deficiencies or by anomalous behavior such as unknown orbit maneuvers, bad measurements or satellite oscillator anomalies. Bad measurements may be

detected by comparing d e absolute value of measurement residuals against a ti reshold of 4 a_jk , where

ajk is they^'tii diagonal element of A£. If die tiireshold is exceeded, the measurement is typically rejected. If the threshold is exceeded several times in a row, it is a symptom of a sudden anomalous condition or filter divergence. In this case, the state associated with die measurement should be reinitialized. That is, if die measurement

exceeds the tiireshold number of times within a given period, measurement z k and state x_\k are reset to 0 and rows and columns 1 and 2 of P£ are reset to that of P_o. If the measurement residual for z₂ exceeds the tiireshold amount of times, measurement z₂ and state Λ:₃£ are reset to 0 and rows and columns 3 and 4 of P are reset to that of Po.

Slowly changing anomalies may be accommodated by adapting die process noise to capture the anomaly. This adaptation is a function of die smoothed measurement residual vector and its relationship to their predicted covariance. That is, smooth the residual to avoid responding to noise measurements, as

where die time constant TC_mi_n < TCk⁼ arjjk < TC_max is controlled by the pseudorange measurement variance times a scalar a. The minimum time constant is Δt, so tiiat a > Atlr_max, where r_max is die maximum expected pseudorange measurement noise variance. Then, the residual ratio is calculated:

e_k = RESkA; RESk 181) On average, the value of the residual ratio should be close to 1. This is accomplished by increasing or decreasing d e elements of Qk- The residual ratio is then bounded as β < e_k < l/β 182)

i.e., a bounded version of Equation 181 with 0 < β < 1 thus preventing adaptation from occurring too fast (the result of Equation 181 could ti eoretically be 0). Qk may then adapted as

Q* = e*QiY* 183) where the acceleration terms in the previous Q are scaled to die new acceleration to form Q'. The values of a and β are system configuration parameters.

The outputs of the Kalman Filter are only estimates of the signal control enors. These enors are then passed through a control process to provide commands to die Signal Generator with d e addition of open loop timing commands from die master station 8. The commands described here are tiiose output from the "Control Laws" Block of Figure 25. The code rate (in meters/second) command is determined via a second order loop. The second order loop is used to account for any range acceleration due to GEO motion which in a first order loop would be a constant bias error:

f_{co e},_k = f_code. -₁ + /_«*._-. / + 2ζω _N<code At ^• X- + / ^ 184)

f_code._k = f_codc_k-y +

^• X- 185)

where a>N_,code ^ιs the corner frequency of the second order code frequency code control loop in radians/second and _W_MS is the GEO clock conection rate from die WMS. (ύM.code is related to die time constant of die control loop as

TC_code =

186)

It is sufficient to use a first order loop aided with die code frequency for the control of the carrier frequency. The carrier frequency command if carrier_, k in Hz defined at he C_up frequency) _£V cr J C_v ^ω N, carrier^ „ learner,* = carrier, k-\ + 4 187)

C

_ f_Cψ f carrier, k ~ J code,k ⁺ ^1 carrier, k \ ύθ)

C where &N,carrier is the comer frequency of the first order carrier frequency code control loop in

radians/second. ®Ncarrier ^ιs related to the time constant of the control loop as follows:

TC_carrier = l/ω_N>carrier 189)

Equation 187 is the accumulated coherency enor which should converge to and follow the satellite's oscillator offset defined at die uplink frequency.

As previously discussed, in the WAAS architecture, there are 2 geographically diverse GUSs 10 (primary and backup) serving each GEO satellite 6. This ensures that the GEO ranging signal will not be lost, even in the case of a natural disaster. However, if backup GUS 10 does not have the same timing as primary GUS 10 at transition within approximately 100 nsec, the users' receivers will lose lock on the signal and will take a number of minutes to reacquire the signal from the new primary GUS. Thus, it is important for the two GUSs 10 to be very closely synchronized in time. Otiier designs for a GUS include a detailed hardware synchronization and periodic resynchronization to maintain the GUS clocks within die required time. The prefened process, on die otiier hand, calculates the difference in GUS clocks by using knowledge of the primary GUS controls and die range of die satellite to botii the primary and backup GUSs, as well as die normal pseudorange measurements. This modifies die measurement matrix equation from die primary GUS process:

where R is the line of sight range vector from the GEO to the GUS calculated as die magnitude of die difference in the ECEF coordinates of d e GEO satellite and die GUS. The ECEF coordinates of die GEO must be calculated every second using the following:

a(At)² s = s₀ + v₀Δt + - 191)

where s, v and a are the ECEF vectors in the GEO navigation message and Δt is the time since the time of applicability in the GEO navigation message.

The control law for the code frequency is the same, with equationl92 replacing equation 185.

fcode.k = fcode.k-X + A«fe,*-.^Δ' + ²ζ^ύ} _N,code^At x, +ω N .code 192)

Once the backup GUS switches to the primary mode, a small transient will occur. The master station will determine fast conections and a new GEO Navigation Message (containing the residual clock enor) for the GEO's signal until it (the master station) can establish conections for die new primary GUS.

The backup GUS process is designed to remove the clock difference between the primary and backup GUS as well as generally follow the code rate that would be conect if the backup GUS were to begin operation as a primary GUS. This process includes the novel function of having the backup and primary GUS synchronized to each other to generate the same pseudorange at a user location (whether the primary GUS is operating conectly or not). In conventional receivers, the LI and L2 receivers, which are coherent and which have proportional dopplers, are used to conect for ionospheric delays. In otiier words, die L2 code frequency and code Doppler are basically the same as tiiose received on LI (1.023 MHz), as tiiey differ only by d e rate-of-change of the ionospheric delay. In die system of the invention, receivers are used in GUSs 10 to track the broadcast signal from GEO Satellites 6. As discussed above, die GEO signal consists of two frequencies that are tracked , i.e., the normal LI frequency and a lower powered C-Band frequency (3630.42 MHz) by the GUSs. Refening to Fig. 23 and 24, GUS WAAS receiver 260 translates the C-Band signal to the GPS L2 frequency to allow tracking in die normally available L2 channel of die receiver. However, die C-band signal from die GEO is not coherent. As a result, when the C-band signal is translated to an L2 signal, die L2 signal is not coherent. Thus, the carrier-to-code Doppler ratio of the converted C-band signal is much different tiian on eitiier the LI or L2 signals in the GPS receiver. This may prevent the L2 receiver from tracking die code of die converted C-band signal. That is, the N in Fig. 6, which is 1200 for a C/A code on the L2 signal, is inconect. For the C-Band signal, it should be approximately 3548.8. If the standard ratio for L2 tracking is used, i.e., 1200, the code Doppler aiding will be wrong by a factor of 2.957, which would result in severe code tracking difficulties.

Referring to Fig. 6, in order to overcome ti ese problems, receiver 260 includes means for adjusting the bandwidtii of the receiver (206). As a result, code phase measurements may be performed at a greater bandwidtii than the 1/20 Hz of conventional receivers. Code phase measurement smoothing is performed outside of die receiver at the WMS. By controlling the code tracking bandwidtii in GUS WAAS Receiver 260 via a command to a value desired for obtaining independent measurements, preferably 1 Hz, each code phase measurement is independent of die others, which overcomes die shortcomings of nanow bandwidtii receivers. Receiver 260 also includes means for prograrriming die carrier-to-code Doppler ratio (208) value N . In tiiis way, even tiiough a receiver channel is configured to track an L2 frequency, the carrier-to- code frequency ratio can be programmed to track the incoherent code Doppler scaled from the C-Band Doppler (or any other frequency such as Ku-Band). In a prefened embodiment, the carrier-to-code Doppler ratio is set at 3548.8 to enable the receiver to track the L2 signal which is translated from die C- Band Doppler . Thus, the GUS receiver of the invention permits independent 1 -second code state monitoring that is not be available in conventional receivers and permits tracking by the receiver of die GEO satellite C-Band code that is translated to L2. AUTO RANGING Another significant function preferably performed by GUSs 10 is to provide continuity of en route/nonprecision approach, which can be supported by ensuring that a continuous GEO ranging signal is provided to all users. This reversionary mode configuration requires that a GEO ranging signal and GEO navigation message be provided for at least an hour after the loss of WAAS conections data. GPS signals augmented by WAAS GEO ranging signals and associated GEO navigation messages are sufficient to provide the accuracy required for en route/nonprecision approach in the complete absence of WAAS conections. Furthermore, given the availability of the GEO ranging source, continuity of fault detection by the aircraft is supported with better than 99.986% availability throughout the service volume (assuming seven GEO satellites, 2.5 years of GEO restoration time, and no barometric aiding).

The two main items required for a continuous ranging signal are the satellite ephemeris and a stable pseudorange. The ephemeris gives the location of the satellite. The range/pseudorange is provided via GUS 10. This allows the GEOs 6 to be used as a ranging source in the UDRE position computation process. For GUS 10 to be able to provide an autonomous ranging signal, it must be able to have the ephemeris data for the satellite that is valid throughout the time period of the autonomous ranging mode. The master station 8 propagates the orbit of the satellite every five minutes for the next two hours, forming two hours worth of five minute ephemeris data for the GEO 6. Each master station 8 sends this data to GUS 10 every five minutes where it is stored by GUS 10 for use in the event that a master station 8 fails.

Two possible embodiments to solve the continuity problem include a stand-alone GUS configuration (reversionary mode) and a networked GUS configuration. A stand-alone GUS, without any external communication, can provide a continuous GEO ranging signal and an hour or more of GEO navigation message. This means that in the event of a complete loss of conections, the system can provide more than an hour of accurate en route/nonprecision approach navigation, during which time system maintenance can be performed.

Alternatively, a networked GUS configuration requires that every GUS, plus selected reference stations, be interconnected through a communication network. This alternative configuration could provide an accurate GEO navigation message indefinitely and could provide supplemental fault detection and exclusion from the ground. These additional capabilities, although useful, are not necessarily required. A networked GUS configuration necessitates that continuity of the communication link be guaranteed and that the user's receiver be modified to accommodate supplemental ground-based fault detection and exclusion capability. The additional complexity of the networked GUS configuration is not justified and, therefore, the stand-alone GUS is the prefened configuration.

If the GUS 10 becomes unable to communicate with any master station 8 in the system, it enters autonomous ranging mode. It ceases validation and monitors the satellite to see if the standby GUS 10 is online and became the primary GUS 10. If no broadcast is detected, the primary GUS 10 that has lost communications reverts to on-line and provides ranging information. The GUS 10 continues the primary ranging GUS control process and broadcasts the stored ephemeris information, null messages (used since one message must be sent every second), and fast conections that mark all satellites as not monitored with the scheduling of message types as specified in the WAAS specification. During a loss of system SIS, GEO 6 will provide the required guidance performance.

However, in order to provide the full safety potential of the reversionary ranging mode with respect to simultaneous common mode software failures, all critical software/firmware in GUS 10, including the commercial software, is preferably developed or certified to DO-178B Level B. The WAAS Message Processor 266 provides monitoring and controlling of the GEO ranging signal as well as preparing broadcast messages for WAAS generator 264. WAAS Message Processor 266 serves to ensure that signal monitor and control functions and GEO navigation message broadcast operate continuously. Finally, the WAAS Receiver 260 to GEO 6 channel is preferably monitored to prevent two GUSs from broadcasting simultaneously in the reversionary mode. SERVICE VOLUME MODEL The present system preferably includes an off-line systems engineering analysis tool, defined as a Service Volume Model ("SVM"), to evaluate and analyze the total system performance of the system. Preferably, the SVM provides analysis of position accuracy, availability, continuity, and integrity. The SVM may also be used to validate the various processes within the system. Certain prefened functions of the SVM pertaining to evaluation of accuracy and continuity of accuracy are described as follows.

The Instantaneous Availability Level, at a particular location, of a WAAS, or any other Global Navigation Satellite System, is defined as the probability that the satellite system meets the navigation requirement at that instant of time at that location. For the present system, where the GPS constellation of a number (MGPS) of satellites is augmented with a number (M_GEO) of Geostationary satellites, the IAL may be expressed by the following equation:

A,( l,t) ≡

wherein Pj^GPS and P_j ^GEO are, respectively, the Markov probabilities of (Mcps-i) usable GPS satellites and (M_GEO-J) GEO satellites. The availability of the navigation function is a function not only of the quantity of healthy GPS/GEO satellites in view at a particular user's location, but also the geometry relative to the user. These factors, in combination over time, are major contributors in determining the continuity of accuracy. The summation is done over all possible combination of GPS and GEO failures. Note that N',j is the total combinations having i GPS failures and j GEO failures:

^'NΛ and the combination operator denotes the total cases of k SV's taken from Ν service

volumes as:

The Boolean function bool(.) is used the count the number of cases where the requirement is met:

^bθθl(

The navigation requirement for the present system includes both an accuracy requirement and a continuity of accuracy requirement. The accuracy requirement, A, can be either a horizontal accuracy requirement or a vertical accuracy requirement or both:

A = (HUNE, \d _JR_ΩM_{X jr}S_v,ιj.. ≤ H max )' (](VUNE, \σ,ιj .. ≤ V max )⁷ f)(GEO y^v!^s ≥ GEO mn .r )j) 197) '

wherein HUNE and VUNE are, respectively, the horizontal and vertical user navigation enors,

and GEO-" is the number of visible Geostationary satellites. Hmax and N^ are defined as the

horizontal and vertical accuracy (one-sigma) thresholds, and GEOmin is defined as the minimum visible GEO satellites.

Continuity of accuracy is defined as the probability that the accuracy requirement is met throughout the flight duration, provided that the accuracy requirement was met at the start of the flight.

The problem heretofore with attempting to adequately define the navigational requirement for WAAS is that only the accuracy requirement, A, has been used in the generic IAL in equation (192), i.e.,

198)

even though the navigation requirement should include both the accuracy requirement and the continuity requirement. To date, there have been some attempts to approximate the continuity requirement by simply requiring 5 or 6 satellites to be visible. This approximation is invalid because it depends on the number of satellites rather than the accuracy requirement. Thus, there has been no mathematical formulation of the navigation requirement that incorporates continuity of accuracy.

A new concept, called Instantaneous Continuity Level (ICL), C,_j>n(f), is

introduced at a time instant t for a particular case n of i GPS and j GEO failures. The navigation requirement R(ij,n) is deemed met if the ICL exceeds the continuity requirement, C_mj_n.

A_j (M,t) - 199)

For a particular case n of i GPS and j GEO failures, Cj_j,_n(t), is defined as:

^MGPS ^{~ l M}GEO ~ ^j _GPc _GEn C. (t) ≡ ∑ ∑ Q^U Q^{U U}K.. 200)

U,n _{χ = Q 0} x v ιj,xy

wherein Q_x ^Gps is the probability of having exactly x additional GPS failures per hour (given the probability of a single GPS failure per hour of P_GPS ) and Q_y ^GE0 is the probability of having exactly y additional GEO failures per hour (given the probability of a single GEO failure per hour of PGEO), i.e.,

-yGEO M - Λ

GEO p y- GEO ^■j - y

^'y,j y GEOV i^l- ^ypGGEEOθ) 203) Note that:

wherein N^Tj_j,xy is the total combination of cases where there are x GPS failures and y GEO failures after there were already i GPS and j GEO failures.

Because most losses of continuity of accuracy will occur as a result of a loss of one or more ranging sources during a flight interval, it is these situations and their probabilities that must be exhaustively modeled. The K,j,_Xy term in equation (205) is simply the conditional probability of meeting the accuracy requirement, A , given that there are i+x GPS andj+y GEO failures. This conditional probability is weighted by the compound probabilities of having i+x GPS and j+y GEO failures.

When the probability of single satellite failure, P_GPS and P_GEO. is set to 0 or when the continuity requirement is removed, i.e. Cmi_n = 0, the availability of continuity of accuracy in equation 199 is identical to the availability of accuracy only in equation 198. For a given flight path beginning at time T and lasting for duration D, the

Continuity Condition at the beginning time T is computed as follows. First, the Instantaneous Availability Indicator, β, (R,l,t), at each time-space step is computed. Then, the continuity condition is evaluated over the entire time-space: the continuity condition is met only if the IAI at every time-space point along the flight path is 1. A simple way to determine continuity condition over a flight duration is to ignore the space dimension, i.e., continuity is considered only in the time dimension while the location is stationary. This is acceptable for short-duration analyses, such as for precision approach, in which the location does not change significantly during the flight. In such cases the continuity condition is determined as follows.

The continuity condition is met over the flight duration D = M_D x AT if and only if the Instantaneous Availability Level at each point on the flight duration meets the Minimum Instantaneous Availability Level. In other words, the Continuity Condition Indicator (CCI), χCR,t,D) , at location 1 over the flight duration (t,t+D) is defined as:

A process to determine continuity over a flight duration D taking into account both time and space changes may also be defined. Since, for statistical purposes, one has to account for all possible flight paths (or trajectories), it is not unreasonable to consider all locations within the "reach" of the aircraft during the flight. For example, for an aircraft flying at 600 mph, the geographical area of concern will grow gradually from the original location at time t to some location 600 miles away after 1 hour. First, the Instantaneous Availability Levels are computed over a circle of a 600-mile radius centered at the location of interest with a fine grid (2 degrees by 2 degrees ) and over the time interval of interest (plus the flight duration to cover the flight duration at the end of the interval) in steps of, for example, 2.5 minutes. If one is interested only in surface (i.e., lat/long) rather than volume (i.e., lat/long and altitude), then a three-dimensional anay for an Instantaneous Availability Indicator (IAI), β, [lat] [Ion] [time] is

calculated. This tiiree-dimensional anay pertains to a "cone" in time and space wherein the base of the cone grows from 0 to the 600-mile coverage and the height of the cylinder is 1 hour. To meet the continuity condition, all IAI's on the surface of the cone must be unity. Thus, if s aircraft speed (e.g. 600 mph )

AT - calculation time steps r ( t + AT ) - radius of coverage at time t + k Δr = sΔ

L(t+kDT) = number of calculating locations at t+kDT on the circle of radius r, where L( t) = 1

then the Continuity Condition Indicator may be expressed by:

χ(M,l,t,D) + kAT) 207)

A trajectory is a special case of time-space dimension because the aircraft location at every calculation time step is known. Therefore, the Instantaneous Availability Indicators, computed at every time-space point on the trajectory in time step of AT , must be unity. The Continuity Condition Indicator may be expressed by:

χ( l,t,D) = bool = M_n 208)

The computation of IAL based on equation 193 involves all combinations of GPS and GEO failures using All satellites in their respective constellations. Since failures of hidden satellites (i.e. satellites that are not visible) do not affect the navigation performance, an equivalent but computationally efficient expression may be defined as:

wherein V_GPS and V_GEO are, respectively, the number of visible GPS and GEO satellites at time t. In the above equation, Wj^GPS and Wj^GE0 are calculated from the Markov probabilities using the hidden GPS satellites H_Gps = _Gps - V_Gps and hidden GEO satellites H_GEO = M_GEo - V_GEo-

where

^VGPS _f pc ^VGEO _GFO _r w ™ = ∑ W ^t-^U = 1 212) ι = 0 ' 7 = 0 ^J

The counting is based on combinations of visible satellites instead of the total number of satellites:

(V. fγ \

N^V = GPS GEO 213) i J

Considerable time saving is realized since, on the average for GPS, only about one-third of the satellites are visible. Considering the case of 24 GPS and 3 GEO satellites, the brute-force evaluation of equation 193 for each time-space point will take 2²⁷= 134,217,728 combinations of satellite failures. For an average of, say, 10 visible satellites, the evaluation of equation 194 will take only 2¹⁰ = 1,024 combinations. All things being equal, this represents a computational saving of 131 ,071 : 1.

The computation of ICL based on equation 199 involves all combinations of GPS and GEO failures using all satellites in their respective constellations. Since failures of hidden satellites (i.e., satellites that are not visible) do not affect the navigation performance, an equivalent but computationally efficient expression may be defined as:

N^V v,^χy κ ∑ bool[A. y^v, y ιj,xy,m (0. 215)

N' m = \ y, wherein V_Gps and e number of visible GPS and GEO satellites at time t. w GPS V_GE_O represent th and W_j ^GE0 may then be calculated from the compound probabilities using the hidden satellites

HGPS = M_Gps - GPS and H_GEO ⁼ MGEO - VGEO-

and

The counting is based on combinations of visible satellites instead of total number of satellites:

Considerable time savings may be realized since, on the average for GPS, only about one-third of the satellites are visible. Considering the case of 24 GPS and 3 GEO satellites, the brute- force evaluation of equation 193 for each time-space point will take 3²⁷ = 7.6 trillion of satellite failures. For an average of, for example, 10 visible satellites, the evaluation of equation 194 will take only 3¹⁰ = 59,049 combinations. All things being equal, this represents a computational saving of 129,140,163:1.

Thus, the concept of instantaneous continuity level incorporates both accuracy requirement and continuity of accuracy requirement and may be used to calculate of instantaneous availability level for the WAAS. Additionally, the instantaneous availability level using only the accuracy requirement is a special case of the instantaneous availability level using continuity of accuracy (when the continuity requirement Cm is not required). Furthermore, the equivalent but computationally efficient formulation of instantaneous availability level using continuity of accuracy dramatically reduces the computation time.

While the description set forth above sufficiently discloses and describes the present invention, the following documents attached hereto as Appendix A, the contents of each of which are incorporated herein by reference for all purposes, provide further detailed descriptions of system and processes of the invention: 1) Interface Control Document for the Navstar GPS Space Segment / Navigation User Interfaces, Drawing No. ICD-GPS-200;

2) Global Positioning System Standard Positioning Service Signal Specification, Nov. 5, 1993; 3) Interface Requirement Specification for the Wide Area Augmentation System (WPN

305677-0047), Feb. 2, 1996;

4) Wide Area Augmentation System Driving Requirements Presentation Material;

5) Software Requirements Specification for the Conection and Verification CSCI for the Wide Area Augmentation System (WPN 305677-0047), Jan. 31, 1996; 6) Software Requirements Specification for the Data Collection Processing CSCI for the

Wide Area Augmentation System (WPN 305677-0047), Jan. 30, 1996; and

7) Software Requirements Specification for the GUS Processing CSCI for the Wide Area Augmentation System (WPN 305677-0047), Jan., 1996.

More generally, although the present invention has been described in detail with reference to certain prefened embodiments thereof, variations exist that are within the scope of the invention and, therefore, the invention should not be limited to the description of the prefened embodiments described herein.

Ill

Claims

CLAIMS We claim:

1. A GPS augmentation system which provides augmentation data for a plurality of GPS satellites, wherein each GPS satellite broadcasts GPS signals comprising navigational data on LI and L2 frequencies, the augmentation system comprising:

(i) at least one reference station, each reference station comprising at least two reference receivers for independently receiving the GPS signals from the plurality of GPS satellites and for generating an output signal comprising the navigational data, wherein the navigational data is sufficient to determine the range from the GPS satellite to the reference receivers;

(ii) at least one master station in communication with the reference stations for receiving the output signals from the reference stations, each master station comprising conection and verification processors, and with respect to each reference station, the conection processor receiving the output signal from one of the reference receivers and calculating the augmentation data, the verification processor receiving the output signal from a different one of the reference receivers and calculating the augmentation data;

(iii) means for comparing the augmentation data from the conection processor to the augmentation data from the verification processor to validate operation of augmentation system; and (iv) means for transmitting the validated augmentation data.

2. The GPS augmentation system of claim 1 wherein the augmentation data comprises long term conections and fast conections for each of the plurality of satellites.

3. The GPS augmentation system of claim 1 wherein each reference station verifies the reasonability of the received navigational data prior to generating the output signal.

5 4. The GPS augmentation system of claim 1 wherein each reference station comprises at least two primary reference receivers and at least one standby reference receiver, and wherein the master stations further comprise means for detecting a failure of one of the primary reference receivers and for switching to receive the output signal data from the standby reference receiver in the event of such failure. 0

5. The GPS augmentation system of claim 4 wherein the conection and verification processors each compute standby clock and ionospheric biases for use when switching to receive the navigational data from the standby reference receiver in the event of a failure of one of the primary reference receivers. 5

6. The GPS augmentation system of claim 5 wherein each reference receiver generates a pseudorange measurement based on a C/A code on the LI frequency, the conection and verification processors each compute a pseudorange conection for each satellite, and the standby clock bias is computed as follows:

where: ^╬╡L\,j,C ⁼ "^RL\,ij.A ~ " fo,biasLl,j,A ~ "^RL\.ij,C > and where:

A denotes data from one of the primary reference receivers;

C denotes data from the standby reference receiver; ╬╡ _Ll Γûá _c is a bias residual computed from the pseudorange measurements from the standby

receiver each second;

P^RLl_,ij i^s the pseudorange measured by the reference receivers; δa_{fo bUuL JtΛ} is an estimated receiver clock bias for one of the primary receivers;

i represents the satellite; j represents the location of the reference station; and

N is the number of measurements and estimates made from the data from both the primary and standby reference receivers in an estimation interval.

7. The GPS augmentation system of claim 5 wherein the conection and verification processors each compute an Ll and an L2 bias for each satellite, and the standby ionospheric bias is computed as follows:

(*_2-_l.fcw )^to* = ^"(I ' ^ t ╬ú*., ^U2 - ^Li┬½ ~ ^XL2-LU,ias,j ~ ^XL2-UJ,iasj )^{A "}(Γêæ₄., ^T L2-Ll.ij )┬╗]

where:

(τ_L2__{Ll J}J)A is a measured slant delay difference for the data from one of die primary reference

receivers; (^r_₂-╬╣_╬╣ i_j )^{B ╬╣s a} measured slant delay difference for the data from the standby reference

receiver; (^χ _L2-_Lu_>huj)^A *^{s a} computed L2/L1 bias for the primary reference receiver, representing a

hardware propagation delay between the Ll and L2 frequencies; (^X _L2-_L\,bi_{as j})^Λ i^{s a} computed L2/L1 bias for the satellite, representing a hardware propagation

delay between the Ll and L2 frequencies; and N is the number of measurements and estimates made from the data from both the primary and standby reference receivers in an estimation interval.

8. The GPS augmentation system of claim 1 wherein the GPS signals comprise a signal carrier and a signal PRΝ code, the reference receivers further comprise means for measuring and smoothing the phase of the PRΝ codes, and each of the master stations further comprise means for unsmoothing the PRΝ code phase measurements to provide an independent measure of the status of the signal PRΝ code.

9. The GPS augmentation system of claim 1 wherein each reference receiver further comprises means for detecting multipath enor originating at a particular failed satellite.

10. The GPS augmentation system of claim 9, wherein the reference stations calculate a pseudorange from each reference receiver to each GPS satellite, the failed satellite causes multipath enor resulting in an inconect pseudorange to be computed for that particular satellite, and wherein the means for detecting multipath enor originating at a particular failed satellite comprises:

(i) a nanow and a standard width conelator receiver, each receiver independently receiving the navigational data, extracting multipath enor from the navigational data, and computing a pseudorange to each of the satellites;

(ii) means for storing the pseudoranges computed by the nanow width conelator receiver and the pseudoranges computed by the standard width conelator receiver over a period of time;

(iii) means for computing the average value of the multipath enor extracted from the navigational data prior to computing the pseudoranges; and (iv) means for comparing the average values of the multipath enor extracted over a period of time to detect multipath enor originating at the particular satellite.

11. The GPS augmentation system of claim 10 further comprising means for broadcasting the detection of multipath enor originating at a particular satellite.

12. The GPS augmentation system of claim 2 wherein the reference stations measure a pseudorange from each reference receiver to each GPS satellite, and wherein the master stations further comprise: means for determining GPS satellite orbits, the determination having an associated covariance, and means for determining ionospheric grid point delays, the determination having an associated covariance, and wherein the long term conections and the fast conections each have an estimation enor, the system further comprising: means for applying the fast and long term conections to the pseudorange measurements; means for combining the estimation enors of the fast and long term conections; means for calculating a user differential range enor representing an enor bound on the combination of the estimation enors of the fast and long term conections; and means for calculating estimates of the calculated user differential range enor using conection enor residual samples representing the difference between the range to the satellite and the reference station pseudorange measurements after applying the fast and long term conections, and the a priori covariance data associated with the means for determining GPS satellite orbits and with the means for determining ionospheric grid point delays so that the user differential range enor bounds the enor of the long term conections and the fast conections with a probability of 99.9%.

13. The GPS augmentation system of claim 12 wherein the means for calculating a user differential range enor comprises: means to compute and store the conection enor residual samples; and Bayesian statistical means to optimally blend the stored conection enor residual samples with the a priori covariance data associated with the means for determining satellite orbits and with the means for determining ionospheric grid point delays.

14. The GPS augmentation system of claim 1 wherein the master stations further comprise means for determining ionospheric grid point delays at each of a plurality of ionospheric grid points, the determination having an associated covariance, and wherein the augmentation data further comprises ionospheric conections and a grid ionospheric vertical enor that bounds the ionospheric grid point delays by 99.9%, the system further comprising: (i) means for computing the grid ionospheric vertical enor; and (ii) means for calculating estimates of the computed grid ionospheric vertical enor using ionospheric enor residual samples representing the difference between the ionospheric grid point delays and the transmitted ionospheric conections for each particular ionospheric grid point, and a priori covariance data associated with the determination of ionospheric grid point delays so that the grid ionospheric vertical enor bounds the enor of the ionospheric conections with a probability of 99.9%.

15. The GPS augmentation system of claim 14 wherein the means for computing the grid ionospheric vertical enor comprises: means to compute and store the ionospheric enor residual samples; and Bayesian statistical means to optimally blend the stored ionospheric enor residual samples with the a priori covariance data associated with the means for determining ionospheric grid point delays.

16. The GPS augmentation system of claim 1 further comprising safety monitor means for independently validating the augmentation data to ensure that hazardously misleading augmentation data is not transmitted by die system.

17. The GPS augmentation system of claim 16 wherein the safety momtor means further comprises means for comparing the augmentation data computed by the conection processor to the augmentation data computed by the verification processor.

18. The GPS augmentation system of claim 16 wherein the augmentation data further comprises a user differential range enor, a grid ionospheric vertical enor, and long term and fast conections, the safety monitor means further comprising means for verifying that the user differential range enor and grid ionospheric vertical enor bound the enor in the transmitted long term and fast conections.

19. The GPS augmentation system of claim 1 further comprising at least one ground uplink station, the ground uplink station receiving the augmentation data from the master stations and generating and transmitting uplinked signals comprising the augmentation data to at least one GEO satellite, wherein the GEO satellite is capable of downlinking the augmentation data.

20. The GPS augmentation system of claim 19 wherein the ground uplink station comprises signal control means for generating signals so that the GEO satellite transmissions mimic a conventional GPS signal, and the uplinked and downlinked signals have a PRN code and a carrier frequency, the signal control means comprising: (i) means for controlling the rate of change of the uplinked PRN code so that it is essentially equivalent to the rate of change of the downlinked PRN code, the downlinked PRN code being adjusted for ionospheric conections; and

(ii) means for controlling the carrier frequency so that it remains in a 1540 ratio with the downlinked PRN code.

21. The GPS augmentation system of claim 19 wherein the ground uplink stations comprise a primary and a backup ground uplink station, and wherein the backup ground uplink station is synchronized in time to the primary ground uplink station to ensure that the backup ground uplink station is capable of transmitting a signal having a PRN code essentially equivalent in timing as the PRN code transmitted by the primary ground uplink station.

22. The GPS augmentation system of claim 21 wherein the GEO satellites broadcast GEO signals comprising an Ll frequency signal and a C-Band frequency signal and wherein each ground uplink station further comprises:

(i) a GUS receiver configured to receive the GEO signals;

(ii) means for translating the C-Band frequency signal of the received GEO signal to a GPS L2 frequency to allow tracking in an L2 channel of the GUS receiver; and

(iii) means for adjusting the code tracking bandwidth of the GUS receiver to track the GEO signals.

23. The GPS augmentation system of claim 22 wherein the GUS receiver further comprises means for programming a carrier-to-code Doppler ratio value to track C-band signals translated to an L2 frequency band with an L2 receiver.

24. The GPS augmentation system of claim 21 wherein the ground uplink stations further comprise means for generating and transmitting a GEO ranging signal and a GEO navigational message.

25. A method for providing augmentation data for a plurality of GPS satellites, wherein each GPS satellite broadcasts GPS signals comprising navigational data, the method comprising the steps of:

(i) providing at least one reference station, each reference station comprising at least two reference receivers;

(ii) independently receiving the GPS signals from the plurality of GPS satellites by each of the reference receivers;

(iii) generating an output signal comprising the navigational data, wherein the navigational data is sufficient to determine the range from the GPS satellite to the reference receivers;

(iv) providing at least one master station in communication with the reference stations, each master station comprising conection and verification processors, and with respect to each reference station, the conection processor receiving the output signal from one of the reference receivers and calculating the augmentation data, the verification processor receiving the output signal from a different one of the reference receivers and calculating the augmentation data; (v) comparing the augmentation data from the conection processor to the augmentation data from the verification processor to validate operation of augmentation system; and (vi) transmitting the validated augmentation data.

26. The method of claim 25 wherein the navigation data comprises pseudoranges, long term conections, and fast conections for each of the plurality of satellites, the long term conections and fast conections each having an estimation enor, the method further comprising the steps of:

(i) determining GPS satellite orbits, the determination having an associated covariance;

(ii) determining ionospheric grid point delays at each of a plurality of ionospheric grid points, the determination having an associated covariance;

(iii) determining the long term and fast conections; (iv) applying the fast and long term conections to the pseudoranges; (v) combining the estimation enors of the fast and long term conections;

(vi) computing a user differential range enor representing an enor bound on the combination of the estimation enors of the fast and long term conections;

(vii) computing conection enor residual samples by computing the difference between the true range and the pseudorange measurements after applying the fast and long term conections; and

(viii) calculating estimates of the computed user differential range enor using the conection enor residual samples and a priori covariance data associated with the means for determining satellite orbits and ionospheric grid point delays, so that the user differential range enor bounds the enor of the long term conections and the fast conections with a probability of 99.9%.

27. The method of claim 26 wherein the step of computing a user differential range enor comprises:

(i) computing and storing the conection enor residual samples; and (ii) optimally blending, using Bayesian statistics, the stored conection enor residual samples with the a priori covariance data associated with the means for determining satellite orbits and with the means for determining ionospheric grid point delays.

28. The method of claim 25 wherein the augmentation data further comprises ionospheric conections and a grid ionospheric vertical enor that bounds the ionospheric grid point delays by 99.9%, the method further comprising the steps of:

(i) determining ionospheric grid point delays at each of a plurality of ionospheric grid points, the determination having an associated covariance;

(ii) computing the grid ionospheric vertical enor; and

(iii) calculating estimates of the computed grid ionospheric vertical enor using ionospheric enor residual samples representing the difference between the ionospheric grid point delays and the transmitted ionospheric conections for each particular ionospheric grid point, and a priori covariance data associated with the determination of ionospheric grid point delays so that the grid ionospheric vertical enor bounds the enor of the ionospheric conections with a probability of 99.9%.

29. The method of claim 28 wherein the step of computing the grid ionospheric vertical enor comprises:

(i) computing and storing the ionospheric enor residual samples; and (ii) optimally blending, using Bayesian statistics, the stored ionospheric enor residual samples with the a priori covariance data associated with the determination of ionospheric grid

point delays.

30. A method for improving the accuracy of the navigational data received by the reference stations of the GPS augmentation system of claim 1, the method comprising the steps of:

(i) each of the reference stations measuring, from each of the GPS satellites, an Ll code phase pseudorange, an Ll carrier phase pseudorange, and an L2 carrier phase pseudorange, the pseudoranges representing the range from the GPS satellite to the reference receiver; (ii) conecting the Ll and L2 carrier phase pseudoranges for ionosphere code and carrier divergence, providing an altered Ll carrier phase pseudorange that matches a multipath- free, noise-free Ll code phase pseudorange;

(iii) calculating the difference between the Ll code phase pseudorange and the altered Ll carrier phase pseudorange; (iv) calculating an ambiguity by smoothing the difference between the Ll code phase pseudorange and the altered Ll carrier phase pseudorange; and

(v) adding the ambiguity to the altered Ll carrier phase pseudorange to generate an unambiguous Ll pseudorange.

31. The method of claim 28 wherein the step of calculating the ambiguity includes the step of incorporating a multiplier that accounts for a variation in the conelation time constant of the enor.

32. A method of minimizing an enor in the synchronization between a ranging signal and a GPS signal while maintaining the coherence of the ranging signal, the ranging signal being generated and re-transmitted by a ground station and transmitted through a geosynchronous satellite, the method comprising the steps of: the ground station receiving the pseudorange data originally generated by the ground station and transmitted through the geosynchronous satellite; the ground station calculating using a Kalman filter estimates of the enor between the synchronization of the received ranging data and the GPS signal; using the enor estimates, the ground station generating control parameters for controlling a signal generator in the ground station to modify the ranging signal so as to minimize the enor in the synchronization between the ranging signal and the GPS signal; and the ground station transmitting the modified pseudorange signal to the geosynchronous satellite for re-transmission.