2.3. Parameter Estimation

GAMIT (program solve) incorporates a weighted least squares algorithm to estimate the relative positions of a set of stations, orbital and Earth-rotation parameters, zenith delays, and phase ambiguities by fitting to doubly differenced phase observations. Since the functional (mathematical) model relating the observations and parameters is non-linear, GAMIT produces two solutions, the first to obtain coordinates within a few decimeters, and the second to obtain the final estimates. (See the discussion in Section 2.2 of the GAMIT Reference Manual.)

In current practice, the GAMIT solution is not usually used directly to obtain the final estimates of station positions from a survey. Rather, we use GAMIT to produce estimates and an associated covariance matrix (“quasi-observations”) of station positions and (optionally) orbital and Earth-rotation parameters which are then input to GLOBK or other similar programs to combine the data with those from other networks and times to estimate positions and velocities [Dong et al., 1998, Feigl et al., 1993]. GLOBK uses a Kalman filter (equivalent to sequential least squares if there are no stochastic parameters in the solution) which operates on covariance matrices rather than normal equations and hence requires that you specify a non-infinite a priori constraint for each parameter estimated (see, e.g., Herring et al. [1990]). In order not to bias the combination, GAMIT generates the solution used by GLOBK with loose constraints on the parameters. Since phase ambiguities must be resolved (if possible) in the phase processing, however, GAMIT also generates several intermediate solutions with user-defined constraints before loosening the constraints for its final solution. These steps are described in detail in Section 3.4 of the GAMIT Reference Manual.

In parameter estimation based on least-squares, the conventional measure of goodness-of-fit is the \(\chi^2\) (chi-square) statistic, defined for uncorrelated data as the sum of the squares of each observation residual (post-fit observed minus computed observation, “\(\text{o} - \text{c}\)”) divided by its assigned uncertainty. In a GPS analysis parameter correlations arise (even if not represented in the original data weights) so the computation of \(\chi^2\) in GAMIT or GLOBK involves a complex matrix operation (see Dong et al. [1998]), but the idea is the same. The value of \(\chi^2\) is usually normalized by dividing by the “degrees of freedom” (\(df\)), the number of observations minus the number of parameters estimated, so that the ideal value for properly weighted, independent random observations is 1.0. Conversely, with whatever a priori weight is assigned to the observations, multiplying the estimated parameter uncertainties by the square root of \(\chi^2/df\) yields the “formal errors” (“formal standard deviations”) of the parameter. With white noise, these formal uncertainties will be realistic and will be inversely proportional to the square root of the number of observations (i.e., will depend strongly on the sampling interval).

The uncertainties in a GNSS analysis, however, cannot be treated with white noise statistics because errors with temporal correlations dominate both the phase observations and estimates of station coordinates (the quasi-observations input to GLOBK). In the phase residuals (from cview or the sky plots produced by sh_gamit), the visible noise—from multipathing and tropospheric fluctuations—is typically correlated over spans of 15–30 minutes. This implies that only samples taken at these intervals are independent, and, to a first approximation, we would get realistic uncertainties by multiplying the formal errors by the square root of the ratio of this interval to the sampling interval used—e.g., for the 2 minute sampling commonly used in solve, we would increase the uncertainties by a factor of 3–5. There also errors with longer correlation times that do not show up in the residuals but are absorbed into the parameter adjustments. Assessing the magnitude of these errors requires us to use noise visible in the residuals (phase or coordinates) to infer the character of the noise at lower frequencies. There are a number of excellent studies of the character of the errors discussed in the presentation “Error_Analysis.ppt” from our most recent workshops on the Documentation page of the GAMIT/GLOBK web site (e.g., Mao et al. [1999], Dixon et al. [2000], and Williams [2003]). Once you have adopted a particular weighting of the data, it is often possible to use external knowledge of the expected behavior of the coordinates or velocities to validate the uncertainties; see McClusky et al. [2000], Davis et al. [2003], and McCaffrey et al. [2007].

In GAMIT/GLOBK there are several ways you can control the uncertainties you obtain for coordinates and velocities, and it is important for you to keep clearly in mind how each of these operates. The uncertainties generated by solve and passed to GLOBK in the h-file are determined by the a priori error assigned to the phase observations and by the sampling interval—solve does not rescale by the square root of \(\chi^2/df\) (“postfit nrms” in the q-file). In the initial (“preliminary”) solution, we normally assign an uncertainty of 10 mm to each one-way L1 phase. By Equation (1), the assigned uncertainty in an LC phase becomes 32 mm. The mean rms of one-way LC residuals is typically ~ 6–9 mm, so the nrms computed by solve is 0.2–0.3. In the second (“final”) solution, we normally reweight the observations using a constant and elevation-dependent term computed in data editing by program autcln from the actual (one-way LC) phase residuals. In order to keep the overall weighting approximately the same as with the 10 mm constant error, the values computed by autcln (ATELV table in file autcln.post.sum) are multiplied by an arbitrary factor of 1.7 (in script sh_sigelv) before being input to solve (via the N-file). We chose to use inflated values of the a priori phase error and not rescale by the nrms in order to generate coordinate uncertainties that (in the presence of correlated noise) are approximately realistic with 2-minute sampling. An equally valid approach would be to rescale by the nrms (i.e. make \(\chi^2/df = 1.0\)) and compensate later for the unrealistically low coordinate uncertainties. With whatever weighting you use in solve, you can increase the coordinate uncertainties used by GLOBK by rescaling all covariances on the h-file or by adding white noise or random-walk noise to the variances of individual stations. We discuss in Section 5.3 why the latter approach is usually preferred.

To obtain meaningful estimates of crustal motion, it is necessary to define a reference frame by imposing constraints on the solution. These come in two common flavors: With finite constraints, we assign realistic a priori uncertainties to the coordinates and velocities of one or more sites. This is the only option available in GAMIT and it is also available in GLOBK. A minimum (non-redundant) constraint with this approach would involve restricting translation by fixing or tightly constraining the three coordinates of one site, and restricting rotation by constraining Earth orientation. Constraining additional sites provides redundancy, but can distort the network and remove your ability to detect errors in those sites except through an increase in \(\chi^2\). The second type of constraint, available through the glorg program of GLOBK, is generalized. With this approach, we choose as large a set as possible of sites with good a priori values and minimize their adjustments while estimating an overall translation, rotation, and scale (Helmert parameters) of the network. Since all of the frame-defining sites are free to move, outliers can be readily detected and removed. Moreover, with generalized constraints there can be no internal distortion of the network: all realizations of the reference frame with a given data set will differ only by a translation and rotation. See Dong et al. [1998] for a mathematical description of each of these approaches.