Recent Investigations Toward Achieving a One Centimeter Geoid
Daniel R. Roman and Dru A. Smith
NOAA National Geodetic Survey, 1315 East-West Highway, Silver Spring MD 20910 USA
Abstract. New research is discussed that relates to improving techniques of hybrid geoid generation and the use of Helmert condensation of surface free-air gravity anomalies. New approaches are examined for modeling signal remaining between a regional gravimetric geoid, G99SSS, and a database of GPS-derived ellipsoid heights at 6169 spirit-leveled benchmarks. One approach reduced the remaining unmodeled signal by one centimeter or about 20% over the approach used for GEOID99. Helmert condensation also offers a solution to improving regional geoids. This approach models the gravitational effects of the terrain masses using fast Fourier transforms for outer-zone (>4) contributions and spherical, space-domain calculations for the inner-zone. The effects of the 3D masses are removed from the surface anomalies, which are then downward continued to the geoid where the effects of the 2D masses (condensed onto the geoid) are restored. This approach removes the approximation of Faye gravity anomalies and permits downward continuation in the smoother gravity field where the terrain has been removed. Preliminary investigations have begun on this approach that demonstrate the viability of modeling the 3D and 2D terrain masses as compared to other approaches.
Keywords. Geoid, GPS, leveling, Helmert gravity anomalies.
1 G99SSS/GEOID99 processing stream
To address these recent investigations, a review of the process used to develop the gravimetric geoid, G99SSS, and the hybrid geoid, GEOID99, is necessary. Understanding the steps involved in their development will give insight into the areas that have been subsequently studied for improvement.
The steps taken to develop G99SSS (Smith and Roman, 2000a) were nearly the same as those used for the earlier G96SSS model (Smith and Milbert, 1999) where Faye anomalies were used as an approximation of Helmert anomalies. Terrain corrections (TC's) were better determined specifically for the Pacific Northwest region by using the NGSDEM99 one arcsecond model (Smith and Roman, 2000b). More generally, TC's were improved through use of latitudinal bands to minimize errors associated with convergence of the meridian. Finally, the database containing the GPS-derived ellipsoid heights at spirit-leveled bench marks (GPSBM's) had nearly doubled to 6169 points in 1999. The distribution of these points through the United States was much better with regions that had previously been unrepresented now having at least minimal coverage.
Geoid heights determined from these GPSBM's represent independent assessments from the G99SSS model. The residual values between the GPSBM's and G99SSS geoid values provide the basis for defining a conversion surface that is applied to G99SSS to create GEOID99. This surface is generated using homogeneous least squares collocation (LSC) where single values are used for the correlation length, signal variance, and error or data variance. A plot of the empirical covariance function derived from the residuals is given in Figure 1. Plus symbols (+) denote 50 km bins of empirical data while the curve is determined from best fit homogeneous values of an 18.2 cm signal with a 400-km correlation length. The mismatch occurs due to two apparent signals arising from sources in the GPS, the leveling, or the G99SSS geoid model (note the double hump in the empirical data).
Determining a better model for these residuals thus became one primary area of focus for research and is discussed in the next section. Modeling both the apparent signals is necessary to generate better future hybrid geoids as well as to better understand errors in the GPSBM's or the underlying gravimetric geoid model. Several approaches have thus far been examined including the use of minimum curvature gridding of the residuals with a low pass filter applied and the use of iterative-LSC to remove bands of the residual signal present in Figure 1 in different steps.
While better modeling of these residuals will result in better near-term results from hybrid geoids, the other major area of research has been to directly implement Helmert anomalies instead of using Faye anomalies as an approximation, which is discussed in section 3. The focus has been to divide the terrain's gravitational effects into inner- and outer-zones. The inner-zone, having the greatest effect, will be modeled in the data domain. The outer-zone will be modeled using FFT's on a global basis. Removal of these 3D gravitational effects from surface free-air gravity anomalies, downward continuation of the residual anomalies to the geoid, and restoration of the 2D gravitational effects will accomplish the Helmert condensation onto the geoid. Several areas related to this research are ongoing and will be discussed in terms of expected results and limitations.
2 Refinements of hybrid methodology
The first attempt to analyze the GPSBM-G99SSS residuals involved simply gridding the 6169 points at one arcminute using the minimum curvature algorithm of Smith and Wessel (1990). An aspect ratio of 0.7547 (cos 41) was selected to reduce the effects of the convergence of the meridian, and a tension of 0.75 with a convergence limit of 0.1 mm was also selected. The resulting grid was filtered in 0.25 increments starting from 0.0(unfiltered) and going through 4.0 (444-km full wavelength). The resulting grids were then compared to the original 6169 GPSBM's to determine how much power was lost, and a curve representing the tradeoff between power-loss and filtering was generated (top panel in Figure 2). The point of inflection was located in the second horizontal derivative (bottom panel of Figure 2) to determine that about 300 km was the filtering level required to remove the least signal at the longest filter. Wavelengths below that are assumed to be purely noise and wavelength above it are assumed to be purely signal. The comparison of this filtered grid with the GPSBM data resulted in a 5.1 cm RMS error, indicating that this performed worse than in GEOID99.
An iterative-LSC approach was next examined. This approach uses homogeneous LSC values for each iteration but uses different values in each iteration. The first iteration used a correlation length of 550 km with 15.0 and 5.2 cm used for the signal and data sigma values, respectively, and is shown in Figure 3. This signal sigma represents the point where the lower "hump" in Figure 1 would intercept the y-axis, representing the power associated with the broader features in the residuals. The a priori data sigma value was selected to match the a posteriori value generated by the covariance matrix and represents the "noise" of both the uncorrelated observation errors and the correlated signal occurring at wavelengths significantly shorter than the correlation length.
When this grid of correlated components was added to the removed trend and bias, the resulting intermediate hybrid geoid was again compared to the GPSBM's. A grid reflecting transformations between the NAD83 and ITRF96(1997.0) reference datums was also incorporated. The residuals that resulted from comparing the intermediate geoid with GPSBM's thus reflected the above "noise," and these new residuals were modeled to capture the shorter wavelength components. In Figure 4, the second covariance function again uses homogeneous values of a 33-km correlation length, 3.0 cm signal sigma and 2.3 cm data sigma. There was no trend or bias in these residuals, because these had already been addressed by the first conversion surface. The resulting grid of correlated values represented a second conversion surface that was applied to the intermediate geoid to generate a final hybrid geoid. This was compared again to the GPSBM's and resulted in a 3.3 cm RMS remaining signal with 2.4 cm RMS being correlated at a 12-km correlation length. All of these values are improved when compared to the results of the single pass LSC approach taken for GEOID99, which were a 4.6 cm RMS remaining signal with 2.6 cm RMS being correlated at a 23-km correlation length.
The signal represented by the second conversion surface is displayed in Figure 5. Narrow features with significant lateral extent (up to thousands of km's) can be seen to criss-cross the United States. Hence, the second conversion surface doesn't just match random errors but models a correlated signal that may represent systematic errors in the NAVD88 leveling. The signal shown in Figure 5 has a strong spatial correlation to the distribution of NAVD88 adjustment nodes, and this possible relationship is being investigated further.
The iterative-LSC approach will be examined further to see if improvements in gridding will better capture signal in the second iteration and for evidence of a relationship with the NAVD88 leveling adjustment. Other approaches that will be examined include the use of variable data sigma values for the a priori error estimates in the covariance function. As little information is currently available on local geoid and leveling quality, the quality of GPS information will be used as the primary means to weight these error estimates.
3 Development of a direct Helmert approach
In addition to examining the generation of the hybrid geoid, reduction of surface free-air gravity anomalies by Helmert condensation is examined. This involves the removal of 3D mass effects, downward continuation, and restoration of effects of 2D mass effects on the geoid. The approach taken in G96SSS and G99SSS was to approximate this procedure using Faye anomalies. Preliminary investigations support the hypothesis that the order that the downward continuation occurs in the procedure is irrelevant, as long as anomalies are accounted for in the proper space. In Figure 6, three approaches are given that illustrate this point.
In the left panel, the 3D masses are removed, the 2D masses are restored on the geoid, and the surface gravity observation is reduced through Helmert space to the geoid. In the middle panel, the 3D masses are removed, the surface gravity observation is downward continued in "no-terrain" space, and the 2D masses are restored at the geoid. In the right panel, the surface observations are downward continued in real space, the 3D masses are removed, and the 2D masses are restored at the geoid.
The last approach is complicated by the necessity of downward continuation through the 3D masses and was not considered further. The first approach is made difficult because the high frequency component represented by the terrain has been restored to the gravity values in the form of the condensed 2D masses. Downward continuation of this high frequency signal would require great care in selection of filters. The middle approach was selected because downward continuation of the smoother gravity signal (with the 3D masses removed and before the 2D masses are restored) should be more readily accomplished.
Both the 3D and 2D masses will be broken into two zones for determining their total gravitational attraction. The outer zone will be calculated outside a cap radius globally using FFT's, while the inner zone will be determined in the spatial domain. A nominal range of 4 degrees has been selected to separate these zones based on other investigations (Smith, Robertson, and Milbert, 2000). A balance must be struck between having an inner zone that is overly large, requiring excessive computational time, and minimizing the effects of the Gibbs phenomenon at the cap radius truncation. The former would require that the cap radius be small so as to minimize the number of space domain calculations required. The latter would require that cap radius be as large as possible, because the gravitational effects of the masses drop off by 1/r2. Additionally, equivalent space and frequency domain tapers must be devised to ensure that two data sets may be added without duplicating signal.
Figure 7 shows preliminary results for the outer zone computation using FFT's for the gravitational effects of the masses outside of a 4 cap radius from any given point. As an example of the effect of the cap, note the 30 milligal signals both East and West of the Andes mountains (blue-green) in South America and the lower 20 milligal signal along the Andes (purple-blue). The terrain effect of the Andes is masked when you are located on them but present in the immediately adjacent regions.
The space domain calculations will be in spherical coordinates to better reflect physical relationships between the masses and the observation points. This area will be investigated later this year using the Maui HPCC facility.
4 Other future developments
In addition to these major research areas, other efforts are under way to improve future geoid products. A joint North American geoid project has been initiated and may involve the United States, Canada, Mexico, Denmark, and several countries in Central America and the Carribean. This project is expected to add further data sets from the involved countries (both archived and recent observations) into the generation and refinement of future geoids.
Within the United States, the improvements represented by NGSDEM99 are expected to continue such that a one arcsecond digital elevation model will be available for all regions inside the conterminous United States. Three arcsecond models are already available on a global basis and will provide height information for all other regions where such one arcsecond data are unavailable. Additional improvements of United States data are expected in the next few years as a national adjustment will provide an improved ellipsoidal height datum. This improvement should significantly reduce the uncorrelated errors associated with residuals deriving from gravimetric geoid and GPSBM differences. These reductions would directly improve the generation of hybrid geoids because the remaining correlated signal would be represented more clearly in the residuals and then be more effectively modeled.
Impending satellite gravity missions will also provide several iterations of global gravity models with increasingly higher spatial resolutions. The generation of G96SSS and G99SSS involved the remove-and-restore technique using EGM96. EGM96 is also a global, spherical-harmonic model and complete to degree and order 360. However, it was constructed with satellite-only mission data through degree 72 with higher harmonics determined from terrestrial data. Future gravity missions will provide satellite-only models complete through degree 120 and (later) 250 with anticipated 1-cm commission errors at 333 km and 160 km wavelengths, respectively. Use of these independent data sets as reference fields will enable regional solutions for the North American geoid that have increased validity.
5 Conclusions
Recent investigations for modeling the differences between geoid values generated from gravimetric data and those determined using GPS-derived ellipsoid heights on spirit-leveled benchmarks (GPSBM's) provide an improved solution over older modeling methods. These improvements can be implemented now and have reduced the RMS differences from 4.6 to 3.3 cm when the final hybrid geoid model is compared to the GPSBM's. This reduction in RMS difference is likely due to use of shorter correlation lengths to model and remove more signal in these final residuals. While the amplitude of the correlated signal present in these residuals only reduces by 0.2 cm, the correlation length is cut from 23 to 12 km, indicating that more of the remaining power in the residuals is now uncorrelated and related to observation errors.
However, it is desirable to have an overall geoid accuracy of one centimeter. A 1-cm geoid model would be necessary for obtaining GPS-derived orthometric heights at the 2-cm and 5-cm levels anticipated as a follow up to the National Geodetic Survey's national height modernization initiative (http://www.ngs.noaa.gov/PUBS_LIB/NGS-58.html). Hence, other solutions must be investigated.
Directly implementing Helmert condensation by determining the gravitational effects of 3D masses, downward continuation, and restoration of the gravitational effect of the 2D masses is under study. The determination of the masses' gravitational effect has been separated into inner- and outer-zones. The outer-zone is determined outside of a cap radius of 4globally using FFT's, while the inner-zone is determined in the space domain to preserve the relative positions and physics of the discrete mass elements. This method would make fewer simplifying assumptions and would more rigorously determine the free-air gravity anomalies on the geoid for calculation of the geoid.
Additionally, expected improvements in reference global models, increased spatial resolution of digital elevation models, and more extensive gravity anomaly observations will all improve future geoid models. As they are used in current procedures, densely sampled terrestrial gravity observations increase the spatial resolution of the regional geoid and correct underlying errors in the reference model. The anticipated commission errors of impending satellite-only models are expected only to be one centimeter for as low as 160 km spatial resolution. This may eliminate the need of generating a hybrid geoid and offer the possibility of analysis of the NAVD88 leveling network.
Acknowledgments
The National Imagery and Mapping Agency (NIMA) provided a major portion of the NGS land gravity data, and was instrumental in the creation of various 3" and 30" digital elevation data grids in use today. NIMA was also a partner in a joint project with Goddard Space Flight Center, that developed the EGM96 global geopotential model.
References
Smith DA and DG Milbert, The GEOID96 high-resolution geoid height model for the United States, J. Geod., 73: 219-236.
Smith DA, DS Robertson, and DG Milbert, 2000, Gravitational attraction of local crustal masses in spherical coordinates, J Geod (in press).
Smith DA and DR Roman, 2000a, GEOID99 and G99SSS: One arc-minute geoid models for the United States, J. Geod. (in review).
Smith DA and DR Roman, 2000b, A new high resolution DEM for the Northwest United States, Surv Land Inf Syst, (in press).
Smith WHF, and P Wessel, 1990, Gridding with continuous curvature splines in tension, Geophysics, 55(3): 293-305.