GNSS Solutions • May/June 2013
Calculating Time OffsetsHow do you deal with timing differences between GNSSes?GNSS Solutions” is a regular column featuring questions and answers
about technical aspects of GNSS. Readers are invited to send their
questions to the columnist, Dr. Mark Petovello, Department of Geomatics
Engineering, University of Calgary, who will find experts to answer
them. mark.petovello@ucalgary.ca
Share via: Slashdot Technorati Twitter Facebook Q: How do you deal with timing differences between GNSSs? A: All GNSSs inherently depend on precise timekeeping to measure the satellite/receiver time of flight of signal propagation with sufficient accuracy to compute ranges/distances for multilateration calculations. Each GNSS ground segment therefore dedicates considerable effort to maintaining a highly stable atomic time scale as well as the corresponding offset to global standards such as UTC (Coordinated Universal Time). Nevertheless, even after accounting for gross differences between the different GNSS time scales — for example, GPS is approximately synchronized to UTC modulo one second whereas GLONASS is synchronized with UTC less three hours — differences still exist among the time scales at the 10s or 100s of nanosecond level. When scaled by the speed of light, these seemingly small time differences manifest as very large ranging errors that degrade positioning accuracy. The various GNSSs already plan to broadcast the offset between their system time and the other system times. Unfortunately, especially for the newly launched systems, this information is either not available or is not being transmitted in the navigation message. (GLONASS does broadcast their offset relative to GPS, but currently not on all satellites. On May 3, the European Space Agency announced that the four Galileo inorbit validation satellites had begun broadcasting the European GNSS system’s time offset from GPS.) Consequently, until such information is readily available on all signals of all satellites in all systems, users are left to handle the time differences on their own. This article addresses how they can do this. Before dealing with the offsets between the GNSS time scales, let’s briefly review the pseudorange measurement equations. For the purpose of this article, we can write the pseudorange measurement, P, to the ith satellite for a single system (denoted with subscript “sys”), as Equation 1 where ρ is the geometric distance between the receiver and the satellite, b_{sys} is the receiver clock bias for the subscripted system, and ε is the combined effect of all measurement errors. Equation (1) implicitly assumes that each satellite’s clock is corrected to its system’s time scale using data broadcast in the navigation message. A more specific definition for the receiver clock bias is: the difference between the receiver’s estimate of time and the time maintained by the particular GNSS to which the pseudorange is measured. By extension, because every GNSS has a slightly different time scale, each system will have different receiver clock biases even if all of the measurements are made by a single receiver. This article looks at how to estimate the different clock biases (or related parameters). Assuming measurements from two different GNSSs (denoted with subscripts 1 and 2) are available, equation (1) can be used to write Equation 2 However, we can express the clock bias for the second system in terms of the bias of the first system as follows: Equation 3 where Δb_{12} is the time offset (or intersystem offset or, simply, offset) between the two systems. In this scenario, b_{1} can be interpreted as the “reference” system bias and the offset to each other system is computed relative to this. Note that because each GNSS maintains a highly stable time scale, the time offset between any two systems changes vary slowly, typically at a few nanoseconds per day — or on the order of 10 fs/s (i.e., femtoseconds per second; femto = 10^{15}). Obviously, equation (3) can be substituted into equation (2). The question is whether this is worth doing. The answer to this question is addressed next in the context of leastsquares and Kalman filtering.
LeastSquares Estimation If equation (3) is used, measurements from the reference system are only related to a single clock bias, as in the singlesystem case. However, measurements to nonreference system satellites are related to the clock bias for the reference system and the corresponding intersystem offset. This increases the complexity of the software. Therefore I believe that using equation (2) directly is slightly easier to implement, because a measurement to a satellite in a particular system is only related to the clock bias for that particular system.
Kalman Filtering Equation 4 This can then be converted to the discretetime equivalent which is given by Equation 5 Equation 6 For details about the meaning/ interpretation of the above terms and equations and how they are obtained from the continuoustime equations, the reader is referred to one of the many textbooks available on Kalman filtering. Here, the focus is on how to derive models to handle the different GNSS time scales. With this in mind, a very common clock model used for singlesystem GNSS is Equation 7 where d is the receiver clock drift and, w_{b} and w_{d} are the process noise values for the clock bias and clock drift states, respectively. The spectral densities of the latter two terms are defined based on the quality of the receiver’s oscillator. For simplicity, we assume herein that the two process noise terms are uncorrelated, which is not true in general. However, accounting for the correlation between the two processes is well understood and can be easily added. Starting Simple. As a starting point for the multisystem case we will, for the time being, use the substitution in equation (3) such that the state vector consists of the receiver clock bias for system one (b_{1}), the intersystem offset between the two systems (Δb_{12}), and a single clock drift. Two things are worth noting at this point. First, for simplicity we are excluding any position or velocity states from the equations. From the perspective of the system model (i.e., equation (4)), the position and velocity states are unrelated with the clock states and thus can be added or removed as necessary. Of course, any practical system will include such states when processing real data, but these will only serve to complicate the discussion that follows here. Second, by only including a single clock drift state, the model implicitly assumes that both GNSSs have the exact same frequency, which, as mentioned in the opening paragraphs, is generally not the case. We will revisit this assumption again later. Mathematically, the system model for the Kalman filter under consideration can be written as Equation 8 Equation 9 where q_{b} is the spectral density of the clock bias of the reference system, q_{d} is the spectral density of the clock drift, and Δt is the time difference between epochs k and k‒1. In other words, we treat the intersystem offset as a random constant. Adding extra systems to this model is therefore straightforward because one need only add an extra constant to the state vector (and corresponding rows and columns of zero to the process noise matrix). One Clock Bias Per System. As with the leastsquares case, a bit more housekeeping is associated with this model because we need to first define a reference system and then keep track of what intersystem offsets are being estimated. Consequently, we prefer to estimate a “separate” clock bias for each system. To do this, we transform the state vector in equation (8) as follows: Equation 11 Equation 12 where the second state, b_{2}, is obtained directly from equation (3). From equation (11), we can write the following two equations: Equation 13 Equation 14 Equation 15 which is of the same form as equation (4) and can thus be used to derive an equivalent discretetime Kalman filter. Doing so for the state vector in equation (12) gives Equation 16 Equation 17 Equation 18 The peculiar thing about equation (18) is that the process noise for clock bias states is perfectly correlated. In other words, the matrix is rankdeficient. Filter designers would normally shudder at this thought since perfect correlation amongst state estimates — that is, perfect correlation within P — normally leads to numerical instabilities. In this case, however, we are only dealing with the uncertainty added in the Kalman filter’s prediction step, not the uncertainty of the states themselves (i.e., P is still full rank because the measurements from each system still only update the clock bias for that system). In fact, the perfect correlation between the clock biases makes sense because this model was derived from the (albeit inaccurate) assumption that all system time scales have the same frequency, and thus predictions of their respective clock biases are equally uncertain. Equations (16) to (18) are generally easier to program than equations (8) to (10) because — as with the leastsquares case — a measurement from a particular system only relates to the clock bias for that system. Back to Reality. We now return to the inaccurate assumption that all GNSS time scales have the same frequency. Fortunately, since the GNSS time scales drift with respect to each other at a very low rate — again, on the order of nanoseconds per day — the above models would likely be sufficient for relatively short processing periods. For example, if an application processes 24 hours of data at a time, the accumulated drift between systems would be on the order of a few meters (i.e., a few nanoseconds scaled by the speed of light). If this is less than the magnitude of other errors — and the application can tolerate this level of error — then using the foregoing models probably will not adversely affect results. By extension, shorter processing runs will be even less affected. However, what if an application requires a more rigorous model? In that case, the best approach would be to estimate a separate clock drift for each system. For the twosystem case discussed earlier, the system model might look as follows: Equation 19 where the subscript on the clock drift states denotes the system to which it applies. Also, note that each clock drift state is driven by the sum of a common process, w_{d}, and a systemspecific process, w_{Δdi}_{ }(where i denotes the system). The common term models the receiver’s oscillator while the systemspecific term models the drift of each system’s time scale (and is thus likely very small, by comparison). The drawback of this model is that the size of the state vector increases by two for every new system used. An alternative but less rigorous approach would be to update equation (8) to model the intersystem offset (Δb_{12}), not as a random constant, but as a randomwalk process driven by a very small spectral density. Doing this and applying the transformation in equation (12) will yield the same transition matrix as in equation (17), but the process noise matrix in equation (18) would have to be increased by the following matrix Equation 20 where q_{Δb} is the spectral density used to drive the intersystem bias. The net effect of matrix (20) would be to ensure that the process noise for the different clock biases are not fully correlated, thus allowing them to change more independently of each other.
Summary
As these offsets are included in the various navigation messages, the models presented here should — theoretically, at least — degenerate to the singlesystem case. Of course, this assumes that the accuracy of the broadcast offsets is sufficiently high. If not, the above models may still be needed to estimate the residual offsets. Copyright © 2017 Gibbons Media & Research LLC, all rights reserved. 
