Thirty Years of Maintaining WGS 84 with GPS

2.1 WGS 84 Through the Decades

The WGS 84 TRF is realized through a globally distributed network of GNSS monitoring stations maintained by NGA. A list of these stations in a previous realization can be found in the 2014 NGA Standard 36, which describes the WGS 84 reference system (NGA, 2014). Updated positions for the current frame may be found on the website of NGA’s Office of Geomatics (NGA Office of Geomatics, 2024). Additionally, for completeness and convenience, the coordinates defining the last three WGS 84 TRF realizations are provided in Appendix A.

The first World Geodetic System (WGS), dubbed WGS 60, was created by the U.S. Department of Defense in the late 1950s as a “practical geodetic reference system that maintains consistency with the best scientific terrestrial reference system at the time but also retains some stability.” The need for a global system of navigation was driven largely by a desire to relate maps made in various local datums and to support intercontinental ballistic missiles (Slater & Malys, 1998). Note that these systems are distinguished by their defining parameters. Within each system, there have been multiple TRF realizations across the years, representing measured physical points in space consistent with the system.

This modern marvel was made possible in large part by the rise of satellite geodesy. WGS 60 and its successors WGS 66, WGS 72, and WGS 84, the currently active system, were realized through measurements of the Transit System (Slater & Malys, 1998). Described by Black et al. (1976) and Stansell (1971), this predecessor to GPS consisted of several satellites (five in 1976) that transmitted signals at 150 MHz and 400 MHz. The Doppler shift of these two frequencies was measured along with a precise orbit product, enabling users to identify their position anywhere in the world.

WGS 84 was originally developed by NGA’s predecessor agency, the Defense Mapping Agency (DMA), using coordinates for 1,591 locations derived from Doppler-shift-based Transit System measurements, surveyed up through December 1986. Additional surveys were conducted in 1987, including GPS-based measurements, to aid in the development of datum transformations (Defense Mapping Agency, 1992). Around the same time, DMA, which had been producing orbits for the Transit System using a network of tracking stations, began working with Applied Research Laboratories, the University of Texas (ARL:UT) to develop a ground tracking network for GPS. The first station was installed in December 1985, with four others to follow in 1986. These sites, along with others installed in the mid-1990s, remain active today as continuously operating GNSS monitoring stations that define the WGS 84 TRF (Renfro et al., 2012).

Following its initial release, the WGS 84 TRF has been updated multiple times. Because these updates do not fundamentally change the definition of the WGS 84 TRF, as described in NGA Standard 36 (NGA, 2014), they are referred to as different realizations. Table 1 provides a summary of these realizations, along with some details regarding when each realization was implemented. Other than the first WGS 84 TRF, the naming convention adopts the GPS week each realization was implemented. Because of operational constraints, most realizations of WGS 84 TRF were implemented at different times in NGA’s precise orbit product and the GPS control segment. Accuracies represent the per-component offset between a given realization of the WGS 84 TRF and its target ITRF, as calculated via a Helmert transformation.

In 1994, the Naval Surface Warfare Center-Dahlgren Division (NSWCDD), on behalf of the DMA, created the first GPS-only realization of the WGS 84 TRF. Data for this update came from five DMA/NGA GPS monitor stations along with observations provided from five U.S. Air Force GPS tracking stations. Station coordinates were estimated as part of a GPS orbit estimation process that included 24 IGS sites along with DMA and Air Force data. The ITRF91 positions of the IGS sites were used and eight IGS locations were held fixed so that the positions of the other sites could be estimated consistent with this reference frame. Finally, as the data span used was limited, according to Swift (1994), velocities were derived using the NNR-NUVEL1 model (Argus & Gordon, 1991).

This technique was used again for WGS 84 TRF realizations G873 (Malys et al., 1997; see also Cunningham & Curtis, 1996), G1150 (Merrigan et al., 2002), and G1762 (NGA, 2014). Note that between September 2012 and June 2015, nine NGA sites moved. The resulting set of coordinates, used in NGA’s precise orbit products for approximately five years, is designated WGS 84 (G1762’). Because the new coordinates remained aligned to ITRF2008 at the level of 2 cm and each site’s position was estimated individually, G1762’ is not considered a separate realization of the WGS 84 TRF (Malys et al., 2016).

For G1674, NGA employed a different strategy. In 1996, NGA began providing the IGS with 30-s receiver-independent exchange (RINEX) data from its monitoring stations. As a result, ITRF2008 included positions and velocities for the 11 NGA WGS 84 GPS monitoring stations. Rather than using tracking data from IGS sites, the ITRF positions for nine of the NGA GPS monitoring stations were held fixed. The positions for the remaining two NGA sites, 85405 (which had an unexplained 0.9-m offset) and 85413 (which had recently been impacted by an earthquake), and the six U.S. Air Force GPS tracking stations were then estimated as part of the GPS orbit determination process, as was done for previous realizations (Wong et al., 2012).

The majority of this information can be found in NGA Standard 36 (NGA, 2014), although the final date for implementation of G2296 in the GPS broadcast orbits was publicized via a notice advisory to NAVSTAR users (Navigation Center, U.S. Coast Guard, 2024). Note that the two most recent frames use a precise point positioning (PPP)-based alignment strategy (denoted as PPP in Table 1), which is described in detail in the following section. Figure 1 provides a map of the stations listed in Table A3 of Appendix A for the most recent realization of the WGS 84 TRF, along with their names in the ITRF2020.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Map of the WGS 84 GNSS monitoring stations provided in Table A3 of Appendix A Sites included in the ITRF2020 have their names marked on the map.

2.2 Alignment Through PPP

Starting with WGS 84 (G2139), released in NGA’s products on January 3, 2021 (NGA Office of Geomatics, 2021), NGA adopted a new method of aligning to the ITRF. Instead of solving for positions as part of an orbit estimation process, the IGS’s final precise ephemeris (PE) product (NASA, 2024) is used as a surrogate for the ITRF. These orbits represent a combination of solutions for nine analysis centers (Griffiths, 2018) and employ a timescale that remains close, but not directly tied, to GPS time (Konyk, 2021). Because the IGS orbit products are created in the IGS reference frame, itself containing hundreds of ITRF fiducial points, PPP combined with these ephemerides places a GNSS site in the ITRF. NGA adopted PPP because, with the software available at the time, it provided better alignment to the ITRF than the orbit-based network technique used for previous realignments and enabled direct velocity estimation.

Thus, to align the WGS 84 TRF to the ITRF, NGA processed raw receiver observations from WGS 84 GNSS monitoring stations using the ARL:UT PPP GPS/ RINEX ARL:UT Precise Position Estimator (GRAPE) software. This tool uses an extended Kalman filter to process dual-frequency GPS observables and provides centimeter-level positioning based on GPS PE data, satellite clock data, 30-s RINEX-format observation data, station and satellite antenna parameters, and geophysical models. This method produces, among other items, daily station positions (with associated uncertainties) in Earth-centered Earth-fixed (ECEF) and geodetic coordinates. For more details on the models used and processing strategy, see the work by Tolman (2008).

For WGS 84 (G2139) and WGS 84 (G2296), NGA employed GRAPE in conjunction with the IGS final orbit products and IGS-provided 30-s clock estimates to generate a daily position estimate for all WGS 84 GNSS monitoring stations from 2016 through an epoch close to the release date of the realization. From here, a station trajectory model was fit using a least-squares regression on the whole time series. This technique enables NGA to simultaneously estimate position and velocity terms for each WGS 84 TRF reference station in a computationally efficient manner.

Note that because the IGS final products are used to relate back to the ITRF, as opposed to orbits from one of the IGS reprocessing campaigns, before fitting a station trajectory, the appropriate 14-parameter transformation must be applied to each part of the data set to ensure a consistent frame. For example, between 2016 and 2021, the date range considered for WGS 84 (G2139), three different IGS reference frames were active: IGS08 (until January 28, 2017), IGS14 (from January 29, 2017 through May 16, 2020), and IGb14 (from May 17, 2020) (International GNSS Service, 2024).

To account for the different frames, the transformation parameters available from the ITRF’s website (IERS ITRS Center, 2023) were applied to the station position time series so that all points were in a consistent frame. It is worth noting that the translation portion of this transformation is on the order of 2 mm per component, making the correction small compared with the centimeter-level, daily uncertainty in the PPP process.

To verify the accuracy of this technique, NGA-produced GRAPE results were compared with positions in daily station-independent exchange (SINEX) products from the IGS. NGA currently provides the IGS with RINEX data from nine WGS 84 GNSS monitoring stations, allowing the IGS to solve for positions of the NGA stations in the IGS frame and providing direct relationships between the two frames. Figure 2 presents a visual representation of the agreement between the two solutions, displaying the magnitude of the distance between positions computed from January 1, 2016, through May 24, 2024. As can be seen, the median difference for all sites is approximately 0.5 cm, and the 75th percentile for all sites is within 1 cm.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

Total magnitude of the difference between the NGA-computed positions and their corresponding positions in the IGS SINEX files for the nine WGS 84 sites providing data to the IGS

The boxes represent the interquartile range (IQR), whiskers correspond to 1.5 times the IQR, and the violin plot highlights the spread in data.

To assess the long-term differences in position, a station trajectory model was fit using both the daily IGS SINEX positions and their NGA-computed counterparts. A comparison of the solution revealed differences in the least-squares position between the GRAPE and SINEX data at the 2024.0 epoch, demonstrating a mean value of less than 0.5 mm per component with a mean distance magnitude of 3.07 mm per station. These differences, evaluated at the 2024.0 epoch, are presented in Table 2.

The standard error, calculated from the covariance matrix output by the GRAPE PPP software, remains at low-to-mid tenths of a millimeter per component per station. This result indicates that the differences observed are statistically significant, with the values in Table 2 being greater than the uncertainty of the PPP process itself. However, the differences are practically insignificant from the perspective of the operational requirements of GPS and the purposes of safe navigation with GPS. For context, the fifth edition of the Global Positioning System Standard Positioning Service Performance Standard asserts an upper bound of 2.0-m accuracy 95% of the time (U.S. Department of Defense, 2020). Furthermore, this result demonstrates the suitability of GRAPE-derived data for station trajectory modeling, as the use of IGS- or NGA-calculated positions results in similar trajectory models.

2.3 Station Trajectory Modeling

While a TRF is physically realized through a series of measured points on the ground, each realization is published by a set of station trajectories that represent nonlinear equations seeking to account for plate tectonic motion, antenna movement, hardware changes, post-seismic deformation due to earthquakes, and other unmodeled geophysical phenomena. Equation (1) presents the form of this trajectory model used by NGA in aligning WGS 84 to the ITRF. Note that this model includes a simplified post-seismic deformation model compared with that used in the ITRF2020 (Altamimi et al., 2023) and described by Brown (2014). In both references, additional terms are applied to more accurately account for the nonlinear motion following an earthquake. NGA only included a single logarithmic term because it was found to adequately model the two WGS 84 stations impacted by substantial earthquakes:

$\begin{array}{cc}x\left(t\right)& =x_0+\dot{x}t+\sum _i\left[A_i\cos \left({\omega }_{i}t\right)+B_i\sin \left({\omega }_{i}t\right)\right]\\ & +\sum _q{C}_q\Theta \left(t-t_q\right)\log \left(1+\frac{t - t_q}{{\tau }_q}\right)+\sum _j{D}_j\Theta \left(t-t_j\right)\end{array}$ Equation (1)

In Equation (1), plate tectonic motion is accounted for by x₀, which represents the Cartesian position component (and generalizes to y or z) at the epoch of the fit, the station’s velocity $\dot{x}$, and the time t (in decimal years) from a given epoch. To handle empirically observed seasonal behavior, periodic terms of amplitude A_i and B_i are included, with ω_i representing one of the two seasonal frequencies with an annual (365.25 days) or semi-annual (182.625 days) period. Earthquakes are accounted for by using a logarithmic term, with t_q being the time of an earthquake, τ_q denoting the relaxation parameter, C_q being the magnitude of the post-seismic deformation, and Θ indicating a step function that begins at the time of the earthquake. Finally, a general step function (Θ) is included to account for any number of jumps in the time series of magnitude D_j occurring at t_j, which could include both earthquakes and station antenna movements.

Ever since WGS 84 (G730), the first realization of the WGS 84 TRF to use GPS, all station positions for the GNSS monitoring sites that define a WGS 84 TRF realization have been reported with a position and velocity. For G1762 and before, velocities were either inherited from nearby ITRF sites or calculated from a global model such as NNR-NUVEL-1 (Argus & Gordon, 1991). This approach was necessary because the few-week span used in the position estimation process was insufficient to independently estimate velocities. Starting with WGS 84 (G2139), however, the PPP-based technique, combined with years of GPS observations, allowed NGA to estimate both position and velocity for each of the WGS 84 GNSS monitoring stations. These positions are given in Table A2 in Appendix A.

The PPP time series used to align WGS 84 (G2139) and (G2296) to the ITRF both began on January 1, 2016. This date was chosen because it represents the point at which NGA sites began flowing data to the IGS following a lengthy hiatus (Maggert, 2015b). During this time span, NGA’s stations in South Australia (85402 in Table A2), Ecuador (QUI3 in Table A2), and New Zealand (MRL1 in Table A2) experienced significant events that required specific processing steps.

The antenna for 85402 moved in mid-2019. For G2139, because this change occurred close to the end of the time span, the velocity was estimated via a linear model fit to data before the move. The new position was computed from an average of 14 days of GRAPE-calculated positions after the move. In G2296, a station discontinuity was estimated following Equation (1) along with the rest of the least-squares process, enabling the full data set to be used for velocity estimation.

NGA’s Ecuador and New Zealand stations both experienced significant earthquakes in 2016, resulting in position discontinuities and nonlinear motion due to the post-seismic deformation. For the Ecuador site, the primary event was a magnitude-7.8 earthquake on April 16, 2016 (U.S. Geological Survey, 2016a), which moved the station by 4.1 cm. The New Zealand event on November 11, 2016 (U.S. Geological Survey, 2016b), also a magnitude-7.8 event, affected the station more severely, initially displacing the antenna by 38 cm.

For G2139, these events were handled by visually inspecting when the station motion became linear after the earthquake and only using data after that point to estimate the station’s position and velocity. Because both displacements occurred early in the time series used for G2139, sufficient data were available to enable this approximation. G2296 improved on this method by estimating the post-seismic deformation parameters, using the linear fit from the model as the final output. This result is visualized for both sites in Figure 3, which shows the (normalized) positions, the full G2296 model, and the linearized G2139 and G2296 positions reported in Tables A1 and A2. Both models represent the data well, but the more advanced modeling in G2296 provides an estimate that appears to better represent the long-term behavior of the station.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

Total distance $\left(\Delta \overrightarrow{r}\right)$ from each site’s initial position in 2016 for the daily PPP-derived position, full G2296 model, and G2139 and G2296 linearized position and velocity terms for NGA’s Ecuador station (QUI3) and New Zealand station (MRL1)

The 14-parameter transformation provided in Table 3 was applied to the G2139 model to bring all data into a consistent frame. The gray shaded region represents the time range used to estimate the G2139 velocity values.

2.4 Updating Satellite Antenna Phase Offsets

While the WGS 84 TRF is realized by the position–velocity model described above, few users of WGS 84 are surveying close enough to these sites to use the publicly available GPS observation data along with a differential positioning technique to reference their data to a WGS 84 GNSS monitoring station. Instead, NGA recommends that users with high-accuracy positioning needs use the NGA PE product (NGA Office of Geomatics, 2020) along with a PPP-like positioning algorithm to arrive at WGS 84 positions for GPS-derived data. However, this process requires the application of consistent satellite antenna phase offsets (APOs), and the values used by NGA have changed over the years.

Historically, NGA used propriety manufacturer-measured APOs provided by the U.S. Air Force (and later U.S. Space Force) in its orbit estimation process. Starting with WGS 84 (G2139), NGA adopted the IGS14 APO values (Rebischung, 2016) for all but the GPS Block III satellites (NGA Office of Geomatics, 2021). This decision was motived by several factors. First, the IGS on-orbit estimates of the GPS Block II satellites represented an improvement in accuracy as compared with the measurements available from the manufacturer. Second, at the time, NGA had access to the now-public Lockheed Martin dual-frequency measured values for GPS Block III satellites (Lockheed Martin, 2021), which NGA could use to account for differences among individual satellites. When comparing NGA and IGS orbits using a method like that employed by the IGS Analysis Center Coordinator (Zajdel et al., 2023), this choice to use a different APO calibration method for GPS Block III satellites led to a several-nanosecond clock offset for these satellites, creating a potential issue for high accuracy if not taken into account.

In WGS 84 (G2296), NGA adopted the IGS20 APO values (Villiger, 2022), specifically those with SINEX code IGS20_2283, for all GPS satellites. This time around, the IGS-calculated Block III values were adopted in part because Lockheed Martin made public the values previously used by NGA. Perhaps more importantly, however, as successive ITRF realizations change by less than 1 cm at the overall frame level, minor differences in scale become noticeable. To more closely align with the ITRF scale and because the IGS-derived APO values are designed to match this value (International GNSS Service, 2023), NGA and the U.S. Space Force jointly adopted the IGS20 APO values in both the NGA PE and the GPS operational control segment (NAVCEN, 2024). As of this writing, the most convenient publicly accessible source of historic NGA APO values is an article by Malys et al. (2020), which includes the APO values used by the U.S. Space Force and NGA up through March 28, 2021.

Because the majority of the APO occurs in the radial direction, the primary impact of their adoption on the PE relates to the estimated satellite clock state, owing to the nature of GPS measurements: a measured distance is the transit time of the signal multiplied by the speed of light. Thus, a 1-m error in distance or a 3-ns error in clock offset leads to the same measurement collected by the user. Consequently, satellite clock offset estimates usually absorb modeling differences along the radial direction. As a result, the APO changes are effectively aliased into clock biases.

As part of the WGS 84 (G2296) update process, NGA produced a PE by using the ITRF2020-aligned positions provided in Table A3 of Appendix A and incorporating IGS20 APO values. This prospective product was then compared against NGA’s official product to empirically determine the impact of transitioning to the use of the IGS20 APO values. The results of this experiment are shown visually in Figure 4, which highlights that after a median bias is removed to account for differences in the underlying timescales, the residual clock differences are close to the theoretical value (gray line in the plot) calculated by the difference between APO values. As expected, adopting the IGS20 values for GPS Block III satellites resulted in a change of approximately 2.85 ns in the clock estimates.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Clock differences between NGA PE products for GPS Week 2286 created using WGS 84 (G2139) positions and APO values and WGS 84 (G2296) positions and APO values

The top plot contains only GPS Block III satellites, whereas the bottom plot contains all other satellites. Note the different scales on the y-axis for both plots. In the plots, the box bounds the 25th and 75th percentiles, the black line in the box represents the median, and the whiskers correspond to the 5th and 95th percentiles. The gray line represents the expected clock offset as calculated from the differences in APO values. This figure was reproduced with permission from work by the NGA Office of Geomatics (2024).

3.1 Transformation Parameters

When comparing successive realizations of a reference frame, NGA uses the Helmert transformation, as it provides the necessary parameters to convert coordinates from one frame to another while preserving the underlying inner geometry and simultaneously estimating translation, rotation, and scale factors (Petit & Luzum, 2010). With successive iterations of the ITRF, the IERS has also created an online tool providing 14-parameter Helmert transformations between any two realizations of the ITRF (IERS ITRS Center, 2023).

Historically, NGA has provided a 7-parameter Helmert transformation along with each new frame update to quantify how much the frame has changed and to provide users with a means for updating positions from one realization to the next. While these parameters are largely scattered across the references highlighted in the previous section, Kelly and Dennis (2021) recently collected this information in a single document, independently computing 14-parameter transformations to convert between successive WGS 84 TRF realizations from WGS 84 (G730) through WGS 84 (G1762), published in Table 2 of that work.

As noted by Kelly and Dennis (2021), at the time of writing, there was not enough public information regarding WGS 84 (G2139) for the authors to calculate a transformation to this frame. Additionally, WGS 84 (G2296) was released after the paper’s publication. To help fill this gap in documentation, Tables A1, A2, and A3 of Appendix A in this work present positions of the WGS 84 GNSS monitoring stations that define the G1762’, G2139, and G2296 realizations of the WGS 84 TRF.

Recently, NGA wrote a program to estimate a 14-parameter transformation that provides the rate of change for each parameter between successive realizations of the WGS 84 TRF. A mathematical derivation of the algorithm used is provided in Appendix A. As a self-consistency check, NGA applied this process to reproduce the ITRF2014-ITRF2020 transformation parameters applied to the PPP data used in calculating the WGS 84 (G2296) model. With only the WGS 84 GNSS monitor stations as input, nearly all of the parameters were identically zero. The three nonzero values were translation in the z direction (–0.1 mm), the scale parameter (0.01 parts per billion [ppb]), and the scale rate (–0.01 ppb/year). These differences are small relative to the ITRF2014-ITRF2020 transformation parameters, however, and are within acceptable error for the comparatively low number of stations involved.

With this approach, NGA provides, for the first time, a 14-parameter transformation between the last three WGS 84 TRF realizations as well as corresponding uncertainties. These transformations are given in Tables 3 and 4, with the T terms corresponding to translation, R indicating rotation in milliarcseconds (mas), and Scale given in parts per billion. The uncertainties given are formal; that is, they are the square root of the diagonal values from the covariance matrix of the ordinary least-squares estimation process employed.

Because the WGS 84 TRF is defined with 17 stations and the differences between ITRF2020 and ITRF2014 are very small, some of the uncertainties are of the same order of magnitude as the values computed. Additionally, for the scale and rate terms in Tables 3 and 4, the cross-correlation components of the covariance matrix are larger than the diagonal terms. This result suggests that the individual components may not be completely resolved; however, it is worth noting that the uncertainties in the rate terms are comparable to those provided with the ITRF2020 (Altamimi et al., 2023). Despite this, and as will be shown in the next section, the new frame represents a detectable improvement in alignment to the ITRF and the transformations provided are expected to remain within the WGS 84 TRF’s operational accuracy requirements (U.S. Department of Defense, 2023).

Note that MRL1 (New Zealand) was excluded from both transformation calculations because of the substantial nonlinear motion occurring following the 2016 earthquake referenced above. Similarly, the parameters in Table 4 were calculated without NGA station 85415 (Western Australia) because this station did not exist while G1762 was active. Finally, this transformation was performed with the WGS 84 (G1762’) coordinates. Because these coordinates remain aligned to ITRF2008 and each station’s new position was computed individually, the parameters in Table 4 are valid for the entire duration for which WGS 84 (G1762) remained in effect in NGA’s PE product.

3.2 Historic Agreement

NGA is responsible for maintaining WGS 84 to support safe, dependable navigation. According to MIL-PRF-32722 (U.S. Department of Defense, 2023), signed on January 5, 2023, NGA is required to evaluate the need for a new WGS 84 TRF realignment should one of four conditions be met:

1. The IERS releases a new realization of the ITRF

2. The root-sum-square (RSS) of a daily 7-parameter transformation between the IGS and NGA precise orbits exceeds 2 cm over a six-month period computed at a mean Earth radius

3. The RSS of the east, north, or up components of PPP positions calculated with NGA precise orbits differs by more than 2 cm over a six-month period relative to their ITRF positions

4. The root-mean-square (RMS) of the NGA sites in the ITRF differs by 2 cm over a six-month period in the east, north, or up component between the NGA-calculated position using IGS orbit and clock products and their corresponding IGS SINEX coordinate

The first item is outside of NGA’s control and has historically triggered the process for a new WGS 84 TRF realization. The last item is primarily a self-consistency check and addresses whether NGA’s process is able to replicate the same positions as the IGS, which is effectively demonstrated in Figure 2 of this paper. The remainder of this section focuses on items two and three, both of which represent meaningful comparisons between the ITRF and the WGS 84 TRF. These items will be explored by directly calculating the distance between reference frames for sites shared between the ITRF and the WGS 84 TRF, by comparing the IGS and NGA PE, and by leveraging NGA’s PE along with raw IGS observation data to calculate WGS 84 positions for ITRF sites.

For consistency between plots, two different metrics will be used to measure agreement. The first is a “Helmert distance” (ℋ), which is described by Equation (4) in the work by Malys et al. (2020). This equation is reproduced in Equation (2) of this work, with r_e being the mean Earth radius as defined in NGA Standard 36 (NGA, 2014) at 6,378,137 m, and the T, R, and s parameters being the Helmert transformation parameters, such as those provided in Tables 3 and 4:

$\mathscr{H}=\sqrt{T_x^2+T_y^2+T_z^2+s^2{r}_e^2+\frac{2r_e^2}{3}\left(R_x^2+R_y^2+R_z^2\right)}$ Equation (2)

Equation (2) presents a clever way to translate a 7-parameter (or 14-parameter) transformation into the average impact at the surface of the Earth. This is accomplished by integrating over 4π steradians under the assumption that the two frames are close enough that one can apply a first-order small-angle approximation, leading to the closed-form solution given in Equation (2). This result is not formally equivalent to the RSS described in item 2, as it weights the rotation component by a factor of 2/3 inside the square root term, but it is used here because it provides a more physically intuitive interpretation of this transformation.

The second metric for measuring agreement is $\Delta \overrightarrow{r}$, representing the RSS of the difference between the Cartesian components of a station or satellite in two frames. This term is equivalent to a geometric sum of the east, north, and up components specified in MIL-PRF-32722 and will be larger than the difference for any individual component. As a result, this metric represents a good upper bound for items 3 and 4.

The analysis in this section focuses on the performance of the last three realizations of the WGS 84 TRF, motivated in part by a desire to quantify the accuracy of the PPP-based technique adopted by NGA in 2021. Table 5 provides a list of transition dates for both the WGS 84 and IGS reference frames, along with their target ITRF realizations during the time period studied, with the IGS dates taken from the IGS website (International GNSS Service, 2024).

Note that an updated version of the ITRF, dubbed ITRF2020-u2023, was released in 2024. This update involved adding stations and updating discontinuities and earthquake models and remains aligned to ITRF2020 with transformation parameters of zero (IERS ITRS Center, 2024). Because it was released after the submission of this paper and because the update should be formally equivalent to the initial ITRF2020 release at the frame level, the authors have not incorporated any analysis of this update in this work.

3.2.1 Direct Comparison

Since at least 1993 (Neilan, 1993), NGA (then DMA) has provided 30-s RINEX observations from its GPS monitoring stations to the IGS, with five NGA sites being included in epoch 1992 experiments (Noll, 1993). In 1996, the DMA Bahrain GPS station (BAHR) began sending daily data to the IGS network (Chase, 1996). In 2008, this data delivery was expanded to include a total of 11 NGA sites (Wiley, 2008). Because of antenna movements during the WGS 84 (G1762) era and IGS’s naming conventions, the IGS sites corresponding to NGA locations were decommissioned (Maggert, 2015a), and new IGS sites representing the new NGA antenna locations were created (Maggert, 2015b).

The result of this longstanding partnership is that ITRF2008 and ITRF2020 have included WGS 84 GNSS monitoring stations in the models published by the IERS. Evaluating these station trajectory models at both past and future times provides a means to directly and unambiguously compare the reference frames using the positions that define each realization, along with a way to estimate how closely they will remain aligned. Such a comparison is shown in Figure 5, which presents both ℋ and the median value of $\Delta \overrightarrow{r}$ for the NGA and IGS frames relative to their respective target ITRF realization listed in Table 5.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

A direct comparison of the WGS 84 (left) and IGS (right) reference frames relative to their target ITRF

Note that because of the station movements necessitating the WGS 84 (G1762’) designation described in Section 2.1, both the G1762’ and G2139 coordinates provided in Tables A1 and A2 were compared with their ITRF2020 values transformed into the target frame listed in Table 5 using the 14-parameter transformations provided by the ITRF. The bottom plot displays the median difference between all stations used in this comparison. The vertical axis units are centimeters for the WGS 84 comparison and millimeters for the IGS comparison. Dashed lines represent values past 2025.

Because antenna movements performed around 2013 necessitated changing the IGS site ID, these new locations were not included in the ITRF2014 realization and were no longer consistent with their ITRF2008 values. Thus, for the WGS 84 (G1762’) and (G2139) frames, the ITRF2020 position of these antennae were converted to ITRF2014- and ITRF2008-equivalent values using the transformation parameters provided by the ITRF. Additionally, NGA stations QUI3 and MRL1 were not included in the G1762’ comparison because of the 2016 earthquakes mentioned previously. Note also that the seasonal ITRF terms were not applied, as neither the WGS 84 nor IGS reference frames (Villiger, 2022) include these terms in published station trajectory models.

Figure 5 clearly shows that the IGS frames remain within millimeters of their ITRF target values in the overall frame sense (ℋ) and remain well under 1 cm for all realizations in terms of the daily median station position differences. This result may be expected, given the fact that these IGS sites define the GNSS portion of the ITRF. There is a noticeable increase in $\Delta \overrightarrow{r}$ for the IGb14 frame, but this increase may be due to the fact that the IGb14 realization included updated earthquake models and velocity estimates, whereas the ITRF2014 was not updated after its publication. Regardless, this plot provides evidence that the IGS realization of the ITRF is a very good surrogate, thereby supporting NGA’s decision to use the IGS precise orbit products to calculate ITRF positions for its WGS 84 GNSS monitoring stations.

Because only nine sites directly relate the two frames together, the Helmert transformation between the WGS 84 and ITRF is significantly more sensitive to modeling differences at individual stations. In comparison, there are hundreds of connections between the IGS reference frame and the ITRF. For example, in this calculation, 195 sites are in both the ITRF2020 and IGS20 reference frames.

Regardless, both WGS 84 (G2139) and WGS 84 (G2296) have remained within 2 cm of the ITRF, as shown by the values of ℋ in Figure 5, and have remained within a few centimeters overall of the ITRF at most sites. Based on the change in site position when the system was updated from G2139 to G2296, it is likely that much of the improvement in station-level accuracy to the ITRF in G2296 is due to improved estimation on the NGA’s part of nonlinear station motion at MTV2 and MRL1. Additionally, these two frames represent an improvement in both alignment and station position accuracy from WGS 84 (G1762’).

Finally, because these parameters are modeled via station trajectories, it is possible to propagate them well into the future and estimate when the WGS 84 TRF is expected to violate the 2-cm threshold outlined by MIL-PRF-32722. For WGS 84 (G2139), ℋ would pass this value in early 2025. The updated frame, G2296, is expected pass this threshold sometime in mid-2034. While this approach does not achieve the level of stability desired for the ITRF of 0.1 mm/year (Altamimi et al., 2023), it remains robust and accurate for safe real-time navigation, where uncertainties for real-time users are expected to be on the order of meters (U.S. Department of Defense, 2020).

3.2.2 Surrogate Comparison

While a direct comparison of the station trajectory models that define a reference frame realization, as performed in Section 3.2.1, is the most unambiguous method for assessing accuracy, most users will access the respective reference frames via precise orbit products. These products represent the most convenient means by which the ITRF, IGS, and WGS 84 reference frames are distributed from their fiducial sites to any point on the Earth.

NGA, then DMA, began producing precise GPS orbits in January 1989, using the OMNIS software developed by the NSWCDD. NSWCDD has been supporting Department of Defense efforts to produce high-accuracy orbits for navigation satellites since the 1960s, including systems such as Transit, NTS-1 and NTS-2 demonstration satellites, and GPS. In 2012, OMNIS was succeeded by EPOCHA, which has been in use by NGA through the present (Swift, 2018).

As part of a regular quality control process since 1995, NGA has computed a daily Helmert transformation and RMS of $\Delta \overrightarrow{r}$ between the NGA and IGS final GPS orbit products, available on NASA’s Crustal Dynamics Data Information System repository spanning several decades (NASA, 2024). Results from this comparison are presented in Figure 6, which displays both ℋ and the RMS of $\Delta \overrightarrow{r}$ spanning nearly 30 years of NGA and IGS production. Because software upgrades also coincide with improvements in modeling, these updates are displayed as vertical lines on the plot. Finally, colors display all GPS-derived WGS 84 TRF realizations to date.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

Comparison of the NGA and IGS final orbit products

Colors represent the active NGA reference frame, and vertical lines mark NGA orbit software updates, with O corresponding to OMNIS and E corresponding to EPOCHA. Note that both plots use a log scale.

As both products represent a best estimate of each satellite’s position, it is not possible to identify which is more accurate from this data set. It is also worth noting that this approach in assessing performance differs from the user range error (URE) method, used in reports such as a performance evaluation of the GPS constellation (Renfro et al., 2021) and derived by Renfro et al. (2024). Rather than evaluating the raw difference between products, the URE is designed to quantify the impact to the user by estimating the portion of the orbit error aligned with the line of sight to a user. As a result, a difference of 10 cm in the orbit estimate of two PE products does not translate to a 10-cm difference in URE. For example, the RMS of $\Delta \overrightarrow{r}$ between the NGA and IGS orbit products for GPS Week 2320 is 10.1 cm, whereas the orbit URE is 3.4 cm.

Figure 6 clearly shows that over the past 30 years, agreement between the NGA and IGS orbits has improved by nearly an order of magnitude. In the G730 era, the median frame difference (ℋ) was 14.0 cm, and the median orbit difference was 92.0 cm. By comparison, the most recent frame, G2296, exhibits a median frame difference of 1.4 cm with a median orbit difference of 10.4 cm.

From the perspective of reference frame alignment, there has not been substantial improvement in the agreement between the NGA and IGS orbit products since G1150. Although there is an increase in ℋ around 2010, the two frames have generally remained within 2 cm of each other, with G1762 being the first to have a median below this threshold at 1.60 cm. Additionally, Figure 6 suggests that changes in NGA’s software have not impacted the underlying reference frame. One exception is the switch from OMNIS to EPOCHA on February 26, 2012, which coincided with improved alignment.

Software changes have had a greater impact on the orbit differences, however, as shown by the RMS of $\Delta \overrightarrow{r}$ in Figure 6. Indeed, the switch from OMNIS to EPOCHA improved agreement between NGA and IGS orbits by roughly 5 cm. Through successive iterations of EPOCHA, agreement continued to improve until reaching approximately 8 cm in EPOCHA V6, V7, and V8. Curiously, at the beginning of 2023, there is an increase in the RMS of $\Delta \overrightarrow{r}$ between the NGA and IGS products. This increase does not correspond with a new version of NGA software or the adoption of IGS20 in the IGS products; it is suspected that the cause is an undocumented change in the quality control process used by NGA.

At the moment, the source of persistent modeling differences between the IGS and NGA orbit products is not well characterized. While it would be convenient to compare the models used by NGA (NGA Office of Geomatics, 2020) with those of the IGS, direct comparison is difficult because the IGS orbit product represents a combination of multiple analysis centers (Griffiths, 2018). As reference frame accuracies continue to improve, identifying and evaluating these differences may be an area of future research.

One of the most obvious modeling differences between the NGA and IGS is the use of different APO values. As previously discussed, prior to WGS 84 (G2139), NGA’s PE used proprietary APO values provided by the U.S. Air Force (and, later, the U.S. Space Force) and its contractors. With G2139, NGA adopted the IGS14 APO values for all GPS satellites except for the GPS Block III satellites; however, it is important to note that this approach was not implemented in NGA’s products until March 28, 2021, 84 days after the new reference frame coordinates were implemented in EPOCHA (NGA Office of Geomatics, 2021). Finally, in G2296, NGA adopted the IGS20 APO values for all GPS satellites (NGA Office of Geomatics, 2024). These changes in modeling did not result in substantially different values for $\Delta \overrightarrow{r}$ in Figure 6, supporting the expectation that APO values, a key factor in achieving optimal alignment to the ITRF scale, primarily impact a satellite’s clock estimate.

3.2.3 Transferred Comparison

The final method presented for assessing alignment between the WGS 84 TRF and the ITRF is to determine the difference between the position of ITRF fiducial sites calculated with NGA’s PE and their modeled position in the target ITRF frame. This method may be the most effective way to verify the utility and accuracy of any orbit product in its capacity to distribute a particular reference frame, as this simulates what a high-accuracy user would experience in using said product.

The analysis presented here begins with the adoption of WGS 84 (G1762) on October 16, 2013, and ends on June 29, 2024, capturing the first several months of G2296. WGS 84 positions were calculated for all available IGS “core” network sites using the NGA PE and GRAPE software described previously. The IGS core network represents a subset of all available IGS sites and was chosen for this purpose because it “was designed to have the best possible global distribution and temporal stability.” The IGS core network is provided as a series of primary stations with a set of substitute sites (Villiger, 2022). For this analysis, positions for all primary stations and their substitutes were computed if daily observation data were available.

One limitation of this study is the fact that the analysis software used was only able to process RINEX 2.11 data. The number of sites processed under each reference frame realization is shown in Table 6. Well over 100 sites are included for all but the G2296 comparison, as most IGS sites have transitioned to newer versions of RINEX. Regardless, the 66 available sites provided sufficient data to estimate a daily Helmert transformation between frames and to compute daily position differences.

View this table:

• View inline
• View popup

TABLE 6

Statistics of the PPP Comparison

The median columns represent the median of the daily Helmert and position offsets plotted in Figure 7. The rightmost column represents the number of IGS core sites successfully processed by NGA using GRAPE for the WGS 84 reference frames and found in IGS daily SINEX products for the IGS reference frames.

In addition to calculating station positions with NGA orbits and GRAPE, the SINEX positions of the IGS core sites were compared with their respective ITRF values in order to identify the relative accuracy of the NGA and IGS position solutions. For both the NGA and IGS comparisons, positions that differed by more than 1 m or seven median absolute deviations from the ITRF position were removed when daily Helmert transformations were calculated. Additionally, to separate position and frame offsets, this daily 7-parameter Helmert transformation was applied to all station positions. Finally, to prevent outliers from obscuring overall trends, the median $\Delta \overrightarrow{r}$ was computed as a means of identifying how much, in general, one would expect an individual station’s position to differ from its ITRF value. Figure 7 summarizes the results of this computation, and Table 6 provides the median values for ℋ and $\Delta \overrightarrow{r}$.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

Frame and position differences between the named WGS 84 TRF and IGS frame realizations and their target ITRF realization, as listed in Table 5

The vertical lines and gray text boxes represent updates to NGA’s EPOCHA software.

As shown in Figure 7, the past three WGS 84 TRF realizations have remained within 3 cm of their corresponding ITRF target frame. This difference was particularly large for G1762 during two periods of time. The first instance ended in late 2015 and was driven by a large x-translation component of the Helmert transformation. Given that this deviation was corrected with the implementation of EPOCHA V2.4 in operations, it is possible that the new software took advantage of updated station positions following antenna movements associated with the WGS 84 (G1762’) coordinates discussed previously.

More curious is the period of time between when EPOCHA V5 became active in production and when G2139 was implemented. During this period, the Helmert distance between the WGS 84 TRF and ITRF began to drift and self-correct. The cause remains unclear; however, it is worth noting that during this time, the z-rotation component of the daily Helmert transformation follows a similar trajectory. This trend appears to primarily impact the frame realization, however, as the increase in position differences during the same period appears to be largely linear.

While there has not been sufficient time to evaluate the long-term performance of WGS 84 (G2296), the behavior of WGS 84 (G2139), the first realization completed using PPP-derived station trajectories, suggests that the new technique can achieve sufficient accuracy and stability levels to satisfy a 2-cm distance threshold in ℋ and $\Delta \overrightarrow{r}$. The frame and position differences remained smaller than those in G1762, and crucially, the median of both metrics remained below the threshold. Additionally, while there appears to be some level of increasing divergence in ℋ, the position differences appear to remain stationary, with a median of 1.79 cm for G2139 and 1.48 cm for G2296. If G2139 is representative of this new technique, then it appears that the WGS 84 TRF and ITRF can be successfully aligned to within 1 cm with a drift of 0.89 ± 0.08 mm/y, estimated via linear regression.

Based on the same comparison metrics, the IGS positions are substantially closer to their target ITRF values. For all IGS frames, ℋ remains below 1 cm, and the median of $\Delta \overrightarrow{r}$ remains below 0.5 cm. Interestingly, there is a periodic component that can be seen in the comparison between the ITRF2014 and IGS14 and IGb14 frames that does not appear to be present when compared with the ITRF2020, which includes periodic models. Overall, this finding demonstrates that the IGS positions convey an accurate representation of the ITRF and supports NGA’s decision to use IGS products as a surrogate for the ITRF when realigning the WGS 84 TRF.