2.  Data Analysis: Theory 

2.1  Introduction to GPS Measurements 

High-precision geodetic measurements with GPS are performed using the carrier beat phase, the output from a single phase-tracking channel of a GPS receiver.  It is the difference between the phase of the carrier wave implicit in the signal received from the satellite, and the phase of a local oscillator within the receiver.  The carrier beat phase can be measured with sufficient precision that the instrumental resolution is a millimeter or less in equivalent path length.   For the highest relative-positioning accuracies, carrier beat phase observations must be obtained simultaneously at each epoch from several stations (at least  two), for several satellites (at least two), and at both the L1 (1575.42 MHz) and L2 (1227.6 MHz) GPS frequencies. The dominant source of error in a phase measurement or series of measurements between a single satellite and ground station is the unpredictable behavior of the time and frequency standards ("clocks") serving as reference for the transmitter and receiver.  Even though the GPS satellites carry atomic frequency standards, the instability of these standards would still limit positioning to the several meter level were it not for the possibility of eliminating their effect through signal differencing.  

A second type of GPS measurement is the pseudo-range, obtained using the 300-m-wavelength CA ("coarse acquisition") code or 30-m-wavelength P ("protected") code transmitted by the satellites.  Pseudo-ranges provide the primary GPS observation for navigation but are not precise enough to be used alone in geodetic surveys.  However, they are useful for synchronizing receiver clocks, resolving ambiguities and repairing cycle slips in phase observations, and as an adjunct to phase observations in estimating satellite orbits. 

For a single satellite, differencing the phases (or pseudo-ranges) of signals received simultaneously at each of two ground stations eliminates the effect of bias or instabilities in the satellite clock.  This measurement is commonly called the between-stations-difference, or single-difference  observable.   If the stations are closely spaced, differencing between stations also reduces the effects of tropospheric and ionospheric refraction on the propagation of the radio signals  If the ground stations have hydrogen-maser oscillators (with stabilities approaching 1 part in 1015 over several hours), then single differences can, in principle, be useful, as they are for VLBI.  In practice, however, it is seldom cost effective to use hydrogen masers and single difference observations in GPS surveys.  Rather, we form a double difference by differencing the between-station differences also between satellites to cancel completely the effects of variations in the station clocks.  In this case the observations are just as accurate with low-cost crystal oscillators as with an atomic frequency standard (though the use of the latter may make editing a bit easier).  

Since the  phase biases of the satellite and receiver oscillators at the initial epoch are eliminated in doubly-differenced observations, the doubly-differenced range (in phase units) is the measured phase plus an integer number of cycles.  (One cycle has a wavelength of 19 cm at L1 and 24 cm at L2 for code-correlating receivers; half these values for squaring-type receiver channels.) If the measurement errors, arising from errors in the models for the orbits and propagation medium as well as receiver noise, are small compared to a cycle, there is the possibility of determining the integer values of the biases, thereby obtaining from the initially ambiguous doubly differenced phase an unambiguous measure of doubly differenced range.  Resolution of the phase ambiguities allows a more precise measure of the relative positions of the stations (see, e.g., Blewitt  [1989], Dong and Bock  [1989] ).

GAMIT incorporates difference-operator algorithms that map the carrier beat phases into singly and doubly differenced phases.  These algorithms extract the maximum relative positioning information from the phase data regardless of the number of data outages, and take into account the correlations that are introduced in the differencing process.  (See Bock et al.  [1986] and Schaffrin and Bock  [1988] for a detailed discussion of these algorithms.)

In order to provide the maximum sensitivity to geometric parameters, the carrier phase must be tracked continuously throughout an observing session.  If there is an interruption of the signal, causing a loss of lock in the receiver, the phase will exhibit a discontinuity of an integer number of cycles.  This discontinuity may be only a few cycles ("cycle-slips") due to a low signal-to-noise ratio, or it may be thousands of cycles, as can occur when the satellite is obstructed at the receiver site.  Initial processing of phase data is often performed using time differences of doubly differenced phase ("triple differences", or "Doppler" observations) in order to obtain a preliminary estimate of station or orbital parameters in the presence of cycle slips.  The GAMIT software uses triple differences in editing but not in parameter estimation.  Rather, it allows estimation of extra free bias parameters whenever the automatic editor has flagged an epoch as a possible cycle slip or, in the program "quick" mode, whenever there is a gap in the data.  Various algorithms to detect and repair cycle slips are described by Blewitt [1990], and also in Chapter 6 of this document.

Although phase variations of the satellite and receiver oscillators effectively cancel in doubly differenced observations, errors in the time of the observations, as recorded by the receiver clocks, do not.  However, the pseudo-range measurements, together with reasonable a priori knowledge of the station coordinates and satellite position, can be used to determine the offset of the station clock to within a microsecond, adequate to keep errors in the doubly differenced phase observations below 1 mm. 

General background and an error analysis for the carrier-beat phase observable may be found in Chapter 5 of King et al.  [1985] and Chapter 2 of Feigl [1991].

2.2  Dual-Band Processing 

A major source of error in single-frequency GPS measurement is the variable delay introduced by the ionosphere.  For day-time observations near solar maximum this effect can exceed several parts per million of the baseline length. Fortunately, the ionospheric delay is dispersive and can be reduced to a millimeter or less (at mid-latitudes) by forming a particular linear combination (LC, sometimes called L3) of the L1 and L2 phase measurements:

	                      f LC   =  2.546 f L1   -  1.984 f L2
 
(See, e.g.,  Bender and Larden [1985], Bock et al. [1986], or Dong and Bock [1989])    Forming LC, however, magnifies the effect of other error sources.  On short baselines where the ionospheric errors cancel in between-station (single) differencing, it is preferable to treat L1 and L2 as two independent observables, rather than form the linear combination.  For baselines longer than a few kilometers, on the other hand, on which ionospheric errors are uncorrelated, it is preferable to form LC and completely eliminate the effects of the ionosphere.  In the general case, the optimal choice of dual-band observable must lie somewhere between these two extremes [Bock et al., 1986; Schaffrin and Bock, 1988].  That is, one must balance the amplification in noise introduced by forming LC against the benefits gained by eliminating ionospheric effects.  In practice, most networks of extent greater than a few kilometers are processed using the LC observable. 

In examining phase data for cycle slips, it is often useful to plot several combinations of the L1 and L2 residuals.  Single-cycle slips in L1 or L2 will appear as jumps of 2.546 or 1.984 cycles, respectively, in LC.  Single-cycle slips in both L1 and L2 (a more common occurrence) appear as jumps of 0.562 cycles in LC, which, though smaller, may be more evident than the jumps in L1 and L2 because the ionosphere has been eliminated.  If the L2 phase is tracked using codeless techniques, the carrier signal recorded by the receiver is at twice the L2 frequency, leading to half-cycle jumps when it is combined with full-wavelength data.  Hence, a jump of a "single" L2 cycle will appear as 0.892 in LC and simultaneous jumps in (undoubled) L1 and (doubled) L2 will appear as 1.654 cycles in LC.  Another  useful combination is the difference between L2 and L1 with both expressed in distance units:

		        f LG    =    f L2     -    0.779   f L1

sometimes called "LG" because the L2 phase is scaled by the "gear" ratio (f2/f1 = 60/77 = 1227.6/1575.42).  In the LG phase all geometrical and other non-dispersive delays (e.g., the troposphere) cancel, so that we have a direct measure of the ionospheric variations.   One-cycle slips in L1 and L2 are of course difficult to detect in the LG phase in the presence of much ionospheric noise since they are equivalent to only 0.221 LG cycles.

If precise (P-code) pseudorange is available for both GPS frequencies, then a "wide-lane" (WL) combination of L1, L2, P1, and P2 can be formed which is free of both ionospheric and geometric effects and is simply the difference in the integer ambiguities for L1 and L2:

		WL = n2 - n1  = f L2  - f L1   +   (P1 + P2 ) (f1 - f2)/(f1 + f2) 

The WL observable can be used  to fix cycle slips in one-way data [Blewitt, 1990] but should be combined with LG and doubly differenced LC to rule out slips of an equal number of cycles at L1 and L2. 

Using wide-lane observations averaged over a satellite pass is the most efficient approach to resolving the wide-lane (n2 - n1 )  phase ambiguities (see, e.g. Blewitt, 1989; Dong and Bock, 1989) though clock jumps and inter-channel receiver biases can sometimes produce spurious results.  To provide a check the algorithm used by GAMIT employs both the pseudoranges and the "phase wide-lane" with ionospheric constraints to resolve the wide-lane ambiguities [Bender and Larden, 1985; Dong and Bock, 1989; Feigl et al., 1993]. 

2.3  Modeling the Satellite Motion

A first requirement of any GPS geodetic experiment is an accurate model of the satellites' motion.  The (3-dimensional) accuracy of the estimated baseline, as a fraction of its length, is roughly equal to the fractional accuracy of the orbital ephemerides used in the analysis.  The accuracy of the Broadcast Ephemerides computed regularly by the Department of Defense using pseudorange measurements from 5 stations is typically 5-10 parts in 107 (10-20 m), well within the design specifications for the GPS system but not accurate enough for the study of crustal deformation.  By using phase measurements from a global network of over 50 stations, however, the International GPS Service for Geodynamics (IGS) [Beutler et al., 1994a], is able to determine the satellites' motion with an accuracy of 5–10 parts in 109 (10–20 cm).  For GPS surveys prior to about 1991, the global tracking network was much smaller but can still be used to achieve accurate results for regional surveys.  If we include in our analysis observations from widely separated stations whose coordinates are well known (from VLBI, SLR, or global GPS measurements), the fractional accuracy of the baselines formed by these stations is transferred through the orbits to the baselines of a regional network.  For example, a 10 mm uncertainty in the relative position of sites 2500 km apart introduces an (approximate) uncertainty of only 1 mm in the components of a 250 km baseline.  This scheme can be used successfully even with regional fiducial sites, transferring, for example, the relative accuracy of 250-500 km  baselines to a network less than 100 km in extent.  

The motion of a satellite can be described, in general, by a set of six initial conditions (Cartesian position and velocity, or osculating Keplerian elements, for example) and a model for the forces acting on the satellite over the span of its trajectory.  To model accurately the motion, we require knowledge of the acceleration induced by gravitational attraction of the sun, moon, and higher order terms in the Earth's gravity field, and some means to account for the action of non-gravitational forces due, for example, to solar radiation pressure and gas emission by the spacecraft's batteries and attitude-control system.  For GPS satellites non-gravitational forces are the most difficult to model and have been the source of considerable research over the past 10 years (see Colombo [1986] Lichten and Bertiger [1989], Beutler et al. [1994b] for more discussion). 

In principle, a trajectory can be generated either by analytical expressions or by numerical integration of the equations of motion; in practice, numerical integration is almost always used, for both accuracy and convenience.  The position of the satellite as a function of time is then read from a table (ephemeris) generated by the numerical integration.  The equations of motion and numerical integrator used by GAMIT were adapted from the Planetary Ephemeris Program, developed originally at MIT's Lincoln Laboratory.  A detailed description of the equations and algorithms may be found in  Ash [1972].

In modeling the phase and pseudorange observations, we must take into account not only the motion center of mass of the satellite but also meter-level offsets between the center of mass and the phase-center of the transmitting antenna.  These offset are negligible for regional networks but can introduce centimeter-level errors for baselines approaching an Earth radius.  Also important are temporary phase excursions of several decimeters lasting up to a half-hour during the maneuvers the satellites execute to keep their solar panels facing the Sun when the orbital plane is nearly aligned with the Earth-Sun direction.  For the satellites in each orbital plane, this alignment occurs for several weeks twice a year, the so-called "eclipse season".  Yoaz Bar-Sever and colleagues at JPL have spent considerable effort developing models of the satellites' orientation, even to point of making the behavior more predictable by getting DoD to apply a small bias about the yaw axis—a change that was implemented gradually between June, 1994, and November, 1995.  See Bar-Sever [1996] for a complete discussion.

2.4  Estimation by Least Squares  

GAMIT incorporates a weighted least squares algorithm to estimate the relative positions of a set of stations by fitting to doubly differenced phase observations.  If tracking data are available from a global network of stations, then the least-squares solution can be rendered non-singular by constraining only the longitude of one station.  In the early stages of processing, however, one often has available data from only a regional or continental network.  In this case it is necessary to fix or constrain at least two horizontal and one vertical site coordinates, and it is sometimes useful to apply a priori constraints to the coordinates of several stations and/or the initial conditions of the satellite orbits (see, e.g., Lindqwister et al. [1990]; Shimada and Bock  [1992]).  The approach usually employed is to tightly constrain the coordinates of at least three "fiducial" stations.  These stations should form a geometric figure which extends as widely as possible in both the north-south and east-west directions.  In the initial analysis, it is also necessary to estimate an independent set of L1 and L2 (or L2-L1)  phase biases.  Finally, a zenith delay parameter for each site can be estimated to account for tropospheric errors.  

Since the functional (mathematical) model relating the observations and parameters is non-linear, the least-squares fit for each session may need to be iterated until convergence, i.e., until the corrections to the estimated station coordinates and other parameters are negligible.  Data can be viewed and edited interactively without iteration since CVIEW uses the pre-fit residuals, partial derivatives, and parameter adjustments to compute and display "predicted" post-fit residuals (see Section 6.5).  Automatic editing is more effective near convergence, however, since AUTCLN uses pre-fit residuals.

In the first stage of GPS processing, the data are usually aggregated by "sessions", each from a few to 24 hours long.  This grouping greatly facilitates bookkeeping and editing.  Once the data are clean, however, it is important to combine many sessions and campaigns in order to get the best estimate of station positions and velocities.  This is most easily accomplished by saving the adjustments and covariance matrices generated by the GAMIT solutions and "stacking" these together in program GLOBK, which uses a Kalman filter to estimate all of the "global" parameters—station coordinates and velocities, orbital initial conditions and non-gravitational force parameters, and Earth rotation parameters [Herring, 1995].  Since GLOBK does not estimate phase ambiguities, it is important that they be resolved in the GAMIT solution.  In cases where the observed network is weak, this may mean performing a multi-session solution in GAMIT or tightly constraining the coordinates in GAMIT using estimates from a previous multi-session (or even multi-year) GLOBK analysis.  Whatever the constraints used in resolving the phase ambiguities, GAMIT always performs an additional solution in which very loose constraints are applied to all parameters; it is this solution, unbiased by the GAMIT constraints (unless they have led to incorrectly resolved ambiguities), which is the primary input to GLOBK (see Herring et al., [1991]; Feigl et al., [1993]).

2.5  References  

Ash, M. E., Determination of Earth Satellite Orbits, Tech. Note 1972-5, 258 pp., Massachusetts Institute of Technology, Lincoln Laboratory, Lexington, 1972.

Bar-Sever, Y., A new model for GPS yaw attitude, Journal of Geodesy, 70, 714–723, 1996. 

Bender, P. L. and D. R. Larden, GPS carrier phase ambiguity resolution over long baselines, in Goad C. C. (ed), Proceedings of the First International Symposium on Precise Positioning with the Global Positioning System, Volume 1, National Geodetic Survey , Rockville, Maryland, 357-361, 1985.

Beutler, G., I. I. Mueller, and R. E. Neilan, The International GPS Service for Geodynamics:  development and start of official service on January 1, 1994, Bulletin Geodesique, 68, 39–70, 1994a. 

Beutler, G., E. Brockmann, W. Gurtner, U. Hugentobler, L. Mervart, and M. Rothacher, Extended orbit modeling techniques at the CODE Processing Center of the International GPS Service for Geodynamics (IGS):  theory and initial results, Manuscripta Geodaetica, 19, 367–386, 1994b 

Blewitt, G., Carrier phase ambiguity resolution for the Global Positioning System applied to geodetic baselines up to 2000 km, Journal of Geophysical Research, 94, 1187-1203, 1989.

Blewitt, G., An automatic editing algorithm for GPS data, Geophysical Research Letters, 1990.

Bock Y., S. A. Gourevitch, C. C. Counselman III, R. W. King, and R. I. Abbot,  Interferometric analysis of GPS phase observation,  Manuscripta Geodaetica, 11, 282-288, 1986. 

Colombo, O. L, Ephemeris errors of GPS satellites, Bulletin Geodesique, 60, 64–84, 1986.

Dong, D.-N., and Y. Bock, GPS network analysis with phase ambiguity resolution applied to crustal deformation studies in California, Journal of Geophysical Research, 94, 3949-3966, 1989.

Fiegl, K. L, R. W. King, T. A. Herring, M. Rotchacher, A scheme for reducing the effect of selective availability on precise geodetic measurements from the Global Positioning System, Geophysical Research Letters, 18, 1289–1292, 1991.

Feigl, K. L., Geodetic Measurement of Tectonic Deformation in Central California, Ph. D. thesis, Massachusetts Institute of Technology, 223 pp., 1991.

Feigl, K. L, D. C. Agnew, Y. Bock, D.-N. Dong, A. Donnellan, B. H. Hager, T. A. Herring, D. D. Jackson, R. W. King, S. K. Larsen, K. M. Larson, M. H. Murray, Z.-K. Shen, Measurement of the velocity field in central and southern California,Journal of Geophysical Research, 98, 21667–21712, 1993.

Herring, T. A., D.-N. Dong, and R. W. King, Submilliarcsecond determination of pole position by GPS measurements, Geophysical Research Letters, 18, 1893-1896, 1991.

Herring, T. A., GLOBK:  Global Kalman filter VLBI and GPS analysis program Version 3.2  Internal Memorandum, Massachusetts Institute of Technology, Cambridge, 1995. 

King, R. W., J. Collins, E. M. Masters, C. Rizos, A. Stolz, Surveying with GPS, Monograph No. 9, School of Surveying, University of New South Wales, Sydney, 1985; reprinted by Ferd. Dummlers Verlag, Bonn, 1987.

Lichten, S. M., and W. J. Bertiger, Demonstration of sub-meter GPS orbit determination and 1.5 parts in 108 three-dimensional baseline accuracy, Bulletin Geodesique, 63, 167-189, 1989.

Lindqwister, U. J., S. M. Lichten, and G. Blewitt, Precise regional baseline estimation using a priori orbital information, Geophysical Research  Letters, 17, 219-222, 1990.

McCarthy, D. D., IERS Standards (1992), IERS Technical Note 13, Observatoire de Paris, July, 1992. 

Murray, M. H., Global Positioning System Measurement of Crustal Deformation in Central California, Ph. D. thesis, Massachusetts Institute of Technology, 223 pp., 1991.

Schaffrin, B., and Y. Bock, A unified scheme for processing GPS phase observations, Bulletin Geodesique, 62, 142–160, 1988. 

Shimada, S., and Y. Bock, Crustal deformation measurements in central Japan determined by a GPS fixed-point network, Journal of Geophysical Research, 97, 11437–12455, 1992.