Home Articles Applications of Kalman Filter in WAAS

Applications of Kalman Filter in WAAS

M. R. Sivaraman

M. R. Sivaraman

Rajiv Ranjan Sahay

Rajiv Ranjan Sahay

Satcom & IT Applications
Area,Space Applications Centre
Ahmedabad -380015, India.
[email protected]

The Kalman filter being a multi-input, multi-output, recursive digital filter that produces estimates of states of a system, which are optimal in the mean square sense, finds application in two ways in GPS based navigation systems like WAAS. The first application is estimating the ionospheric TEC from the GPS observables and the other is in the Time Synchronization of the atomic clocks placed in the RIMS.

The Kalman filter is a multi-input, multi-output, recursive digital filter that can optimally estimate, in real time, the states of the system based on its noisy outputs. The estimates are statistically optimal in the sense that they minimize the mean square estimation error. This filter finds application in GPS based navigation systems, (like Wide Area Augmentation System for Category I precision approach of aircrafts), in two ways, namely, estimation of the ionospheric Total Electron Content (TEC) and Time Synchronization of the Range and Integrity Monitoring Stations (RIMS). A brief approach of using the Kalman filter for estimation of TEC and synchronization of the atomic clocks at the RIMS is given below. Firstly, we look at the problem of estimation of ionospheric TEC in section 2 and 2.1, then we focus on the time synchronization of the various RIMS in section 3.1.

Estimation of Ionospheric TEC
The ionospheric Total Electron Content (TEC) can be estimated from the GPS observables, namely pseudorange and carrier-phase measurements. However, in the estimation of the ionospheric TEC form the GPS observables several instrumental systematic effects such as the biases in the GPS satellites and the GPS receivers on the ground must be modeled. The Kalman filtering technique can be adopted for estimating the instrumental biases as well as the TEC at each GPS station using dual GPS data.

There will be a bias for each of the two GPS frequencies ( f1 =1575.42 and f2 =1227.60 MHz) and their difference, the differential instrumental bias, will produce systematic errors in the estimates of the ionospheric delays. For accurate estimates of the TEC these differential instrumental biases have to be removed.

The GPS carrier phases Likj and pseudoranges Pikj , in range units, can be modeled using the following formulae [Sardon et al]:
L1kj = rij + c(dti- dtj) + ditrop j – diion 1j – l1bi1j
L2kj = rij + c(dti- dtj) + ditrop j – diion 2j – l1bi2j         (1)
Pi1j = rij + c(dti- dtj) + ditrop j + diion 1j + dq1j + dq1i
Pi2j = rij + c(dti- dtj) + ditrop j + diion 2j + dq2j + dq2i
Here, the subscripts k=1,2 refer to frequency fk.

rji the distance from the receiver to the satellite

c the speed of light

dtj, dti receiver and satellite clock offsets, respectively

ditrop j tropospheric delay

lk is the carrier wavelength at frequency fk

bikj initialization constant at frequency fk , lumping together the carrier phase ambiguity and the satellite and receiver instrumental phase delay biases.

dqki , dqkj satellite and receiver instrumental group delay biases at frequency fk , respectively.

diion kjionospheric delay at frequency fk which can be modeled as [Sovers and Fanselow,1987 ]

Figure 1. Coordinates Y and c of the point P in the chosen reference system

dion = af I         (2)

where I is the integrated electron density along the propagation path.

af = 40.3/ f2 , f is the frequency of the signal

We can group the frequency-dependent and the non-dispersive terms separately in equations (1) . Grouping the non-dispersive terms together we get:

pij = rij + c(dti-dtj) + ditrop j         (3)

Now eqn(1) can be written as:

Likj = pij – akIij-lkbikj
Pikj = pij + akIij+qikj         (4)

where ak °afk and qikj = dqkj+dqik includes both the receiver and satellite instrumental group delay biases for each pair station j and satellite i.

However, in this system, the number of unknowns is greater than the number of equations. So we distribute the effect of instrumental delays among other terms defining other variables as:

The term Iji is called as the biased ionospheric term and it contains information about the ionospheric delay Iji and the receiver and satellite differential group delay biases kj, ki .

The cycle-slip detection and estimation of the modified initialization constant terms bkji is made using the P code pseudorange measurements as discussed in Blewitt [1990]. Thus, after the removal of ambiguities, we will only consider the more precise carrier phase observables, with the ambiguities removed.

Now, the next step is the separation of the contribution _of the ionospheric delay Iji and that of the differential instrumental biases kj and ki in Iji .


Estimation of Ionospheric Delay and Instrumental Differential Group Delay Biases
The ionospheric TEC at any observation direction can be modeled as a function F describing the vertical TEC on the subionospheric point of observation, mapped to the line of sight by means of a slant function S:

Iji(t) = S(eji) F(dji, t)         (9)

where t is the observation time, eij is the elevation of the satellite i as seen from the receiver j , dij is the geocentric direction of the subionospheric point.

The slant function S, given by Sovers and Fanselow [1987] can be used here.

where e is the satellite elevation, h1 and h2 are the lower and the upper height of the ionosphere and RE is the mean radius of the earth.

Figure. 2: Local ionosphere around the zenith of the station j. P is subionospheric point of satellite i as seen from station j, and Z is subionospheric point of the zenith of station j.

We consider the geocentric reference system as shown in Figure 1 with axis z towards the solar center. The biased ionospheric term can be modeled as (refer to Figure 2):

Iji (t) = S(eji) [Aj(t)+Bj(t)dYji+Cj(t)dcji]+kj+ki         (11)

where dYji=Yji-Yj, dcji=cij-ci and Aj(t), Bj(t), Cj(t) represent the coefficients of the local linear approximation to the global function F. Aj(t) is the vertical TEC for the zenith direction of station j.

The parameters for which we should solve are the differential instrumental biases and Aj(t), Bj(t), Cj(t) at each station j and for each time t.

Eqn. (11) does not allow an unambiguous solution for kj and ki, so one way out is to choose an arbitrary receiver kr. Therefore, the estimated satellite and receiver are kri=ki+kr and kj,r=kj-kr .

The parameters A, B, C, kj,r, kri are estimated using a Kalman filter. The state equations that are used are:

Aj(t2) = Aj(t1) + Bj(t1)[ Yj(t2) – Yj(t1)] + Cj(t1)[cj(t2) – cj(t1)]
Bj(t2) = Bj(t1)
Cj(t2) = Cj(t1)
kj,r(t2) = kj,r(t1)
kri(t2) = kri(t1)

Time synchronization Implementation
In US WAAS, GPS Common View Time Transfer Technique is used to estimate the difference between the Reference station clocks and the Master clock located at the MCC. The basic equation used to estimate the clock offset is the following pseudorange residual equation:

{(Dk,m = Drk.1k,m + Dbm-DBk+ vk,m)k=1}m=1,M           (1)

In eqn. (1), as stated above, Dr k is the vector that connects the true location of the kth satellite, to its location according to the broadcast navigation message. 1k,m denotes the unit vector from the kth satellite to the mth reference station. Additionally, Dbm and DBk are the offsets in the Reference station and the satellite clocks with respect to the GPS time. The measurement noise is given by vk,m. K is the no. of satellites in view of the mth Reference station and M is the total no. of Reference stations.

As seen from eqn. (1), considering only one Reference station, i.e. m=1, we can say that there are 5 unknowns. They are, the three components of Drk namely, the along track error, the cross track error and the radial track error in the ephemeris of the kth satellite, and Dbm as well as DBk . Since, in the Common View Time Transfer Technique, a single satellite is tracked by all M Reference stations, in subsequent equations i.e. m=2, ….., M, the number of unknowns is incremented by one each time, namely, the corresponding Dbm . Thus we see that the total number of unknowns is 5+M-1 and the number of equations are only M. Moreover, the observations of the clock offsets between the k th satellite and the mth Reference station are embedded in noise vk,m . The clock offset estimation problem can be formulated in terms of the well known system equation and the noisy observation equation of the Kalman filter. Hence, this problem can be easily solved by using the Kalman filtering technique.

In this paper, we have briefly outlined the applications of the Kalman filtering technique in Wide Area Augmentation Systems (WAAS). This technique has the potential of efficiently solving the problems of estimation of ionospheric Total Electron Content as well as the Time Synchronization of the atomic clocks at the RIMS.


  • Blewitt G. An automatic editing algorithm for GPS data, Geophys. Res. Lett., 17(3),199-202,1990.
  • Herring, T. A., C. R. Gwinn and I. I. Shapiro, Geodesy by radiointerferemetry: The application of Kalman filtering to VLBI data, J. Geophys. Res., 95(B8),12561-12581,1990.
  • Sovers, O. J. and J. L. Fanselow, Observation model and parameter partials for the JPL VLBI parameter estimation software MASTERFIT-1987, Jet Propulsion Lab. Publ. 83- 89, Rev. 3,1-60,1987.