2.2. Modeling the Motions of the Satellites and Stations

A first requirement of any GNSS 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 navigation (“broadcast”) orbit computed by the GNSS control centers using pseudorange measurements from < 15 tracking stations is typically 1–5 parts in \(10^7\) (2–10 m), well within the navigation design specifications for the system but not accurate enough for the study of crustal deformation. By using phase measurements from a global network of over 100 stations, however, the International GNSS Service (IGS) [Beulter et al., 1994], is able to determine the satellites’ motion with an accuracy of 1 part in \(10^9\) (2 cm; 5–20 cm in earlier years; see http://acc.igs.org). For GPS surveys prior to 1994, the global tracking network was much smaller but can still be used to achieve accurate results for regional surveys. If we estimate orbital parameters and include in our analysis observations from widely separated stations whose coordinates are well known, 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 uncertainty of only ~ 1 mm in the components of a 250 km baseline. This scheme can be used successfully even with regional fiducial stations, transferring, for example, the relative accuracy of 250–500 km baselines to a network less than 100 km in extent, a helpful approach with surveys conducted prior to the availability of precise global orbits.

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 the gravitational attraction of the Sun, Moon, higher order terms in the Earth’s static gravity field, solid-Earth and ocean tides, and some means to account for the action of non-gravitational forces due to direct solar radiation pressure, radiation reflected from the Earth, radiation emitted by the spacecraft’s radio transmission, and gas emission by the spacecraft’s batteries and attitude-control system. For all GNSS satellites non-gravitational forces are the most difficult to model and have been the source of considerable research over the past 20 years (see Colombo [1986], Beutler et al. [1994], and Ziebart et al. [2002] 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. In GAMIT the integration is performed by program arc using equations given by Ash [1972].

Besides the orbital motion of a satellite, we must take into account meter-level offsets between its center of mass and the phase-center of the transmitting antenna, including temporary 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 the US Defense Department to apply a small bias about the yaw axis to GPS satellites—a change that was implemented gradually between June 1994 and November 1995. See Bar-Sever [1996] and Kouba [2009] for a discussion of GPS yaw, and Montenbruck et al. [2015] for other GNSS.

The position of the ground station in the Earth-centered inertial system defined by the satellites’ orbits is affected by a number of geophysical phenomena. These include the Earth’s rotation, precession and nutation of the spin axis in inertial space, motion of the spin axis with respect to the crust (“wobble”), luni-solar solid-body tides, and loading of the crust by ocean tides and the atmosphere. All of these effects are incorporated into the model of the phase and pseudorange observations computed in program model.

In modeling the phase and pseudorange observations, we must also take into account changes in the apparent distance due to variations in the phase centers of the transmitting and receiving antennas. With matched ground antennas in a regional network, these effects nearly cancel, but for longer baselines and/or different antenna types they can amount to several centimeters in estimated heights. Models for the phase center offsets (PCOs) and variations (PCVs) with elevation and azimuth for most commonly used ground antennas have been determined by electro-mechanical measurements [http://gnpcvdb.geopp.de; http://www.ngs.noaa.gov/ANTCAL], and for the satellite antennas by analysis of global tracking data [Schmid et al., 2007].

Also part of the phase and pseudorange model is the propagation delay caused by the neutral part of the Earth’s atmosphere. This effect is represented by a time-dependent “zenith delay”, a “mapping function” that extends the delay to other elevation angles, and a simple function for north-south and east-west gradients. The zenith hydrostatic (“dry”) delay (ZHD) can be well represented by surface pressure, available either from an empirical function of latitude, longitude, and season (the GPT2 model of Lagler et al. [2013]) or at six-hour intervals from a numerical weather model as computed by TU Vienna (VMF1) [Boehm et al., 2006] and distributed in yearly grid files by MIT for GAMIT users. The hydrostatic and “wet” mapping functions are also part of both the GPT2 and VMF1 models. The component of the zenith delay due to water vapor (zenith “wet” delay, ZWD) and the local gradients cannot be accurately determined even from a numerical weather model and must be estimated from the GPS data. The details of these models are described in Chapter 7 of the GAMIT Reference Manual. As noted above the first-order effect of the ionosphere on the signal delay is made nearly negligible by combining L1 and L2 measurements. However, during periods of high ionospheric activity, 2nd and 3rd order effects at the millimeter to centimeter level are possible and can be modeled in GAMIT using the approach described by Petrie et al. [2010].