|
Progress Towards A Centimeter Level Orthometric Heights Using a Single GPS Receiver 2. Data Used In The Computations 2.1. Gravity data Gravity data for Dubai were surveyed by the German contractor "Hansa Luftbild" in the period 2000-2001. The gravity net was referenced to three absolute gravity stations established in Dubai. The heights of the gravity points were measured with fast static GPS. As part of the gravity processing, gravity values and ellipsoidal heights of the gravity points were converted into conventional free-air and Bouguer gravity anomalies, using the EGM96 geoid model (with the new geoid model the gravity anomalies will change a fraction of a mgal, but with the GPS levelling fit applied this will have no consequence for the geoid). The Bouguer anomalies are smooth in mainland Dubai, but a very large gradient goes through Hatta, cf. Fig. 1.  Fig.1a. Bouguer anomalies of Dubai main area (5 mgal contour interval)  Fig 1b. Bouguer anomalies of Hatta (10 mgal contour innterval) The large horizontal gradient of the gravity anomalies through Hatta makes the gravimetric geoid computation from gravity in this region very difficult - the data area is just too small. Gravity data from neighbouring areas (other emirates of U.A.E. and Oman) were not available for this project. All gravity was checked for outliers, and the marine gravity data compared to the satellite altimetry, to check for possible datum errors. The comparison statistics is shown in Table 1. The statistics show that the region has a quite rough gravity field, and that the marine data fits well.
2.2. Digital elevation models The digital elevation models available consisted of a 30" x 30" DEM from NGDC, USA, covering the complete region, and a number of detailed heights provided by Dubai Municipality. The detailed heights were averaged into a 100 m x 100 m basic DEM, and the 30" DEM used to fill in the missing regions. 2.3. GPS and leveling data A set of approx. 3750 leveled benchmarks with GPS ellipsoidal heights were made available for the computation. The GPS data were tied into the ITRF base network of Dubai, and the levelling referred to a fundamental tide gauge at Port Rashid, Dubai. Most levelling was third order, with some points levelled by trigonometric methods, which was not for this Geoid computation purpose. Many GPS points were repeated RTK measurements (with a 5 cm acceptance limit), while other points in build-up areas were actually determined using classical techniques from nearby GPS points. At points with GPS and leveling a GPS geoid value was derived by: NGPS = hGPS - Hlevelling It should be pointed out that these geoid heights, opposed to the geoid heights determined from global models and gravity data, refer to a local vertical datum, due to the sea-surface topography at the reference tide gauge. In connection with the gravity observations, a levelling line was observed around the perimeter of the Dubai main area, and GPS observations (for gravity station heights) were done in connection with this. The eastern and southern part of the perimeter levelling line GPS was done using rapid static techniques. However, baselines were relatively short, and it appears that the accuracy was good enough also for geoid use (3-5 cm for most points). The perimeter GPS geoid data have therefore also been used for constraining the final geoid, after an iterative editing of outliers, outlined in detail later 3. Computation of The Gravimetric Geoid The Dubai precise geoid has been computed in two steps:
The method of least-squares collocation has also been used for an alternative preliminary geoid computation, primarily to assess the overall accuracy of the geoid. Collocation has the advantage that both steps can be done at once, with a semi-automatic rejection of outliers. The disadvantage is the computational expense, since a very large system of linear equations has to be solved. The collocation error estimates showed that the geoid in the center of mainland Dubai is accurate to 2 cm r.m.s., while in Hatta expected errors are at the 10 cm level. All computations outlined in the sequel have been done by programs of the GRAVSOFT package, © KMS and University of Copenhagen. 3.1 Gravimetric geoid computation: "remove" step In the gravimetric geoid computation, the geoid signal N is constructed by subdividing it into three parts N = N1 + N2 + N3 Where first part comes from a spherical harmonic expansion complete to degree and order 360 (EGM96, cf. Lemoine et al, 1996), the second part from the topography, and the third part from the contributions of "residual" gravity (i.e., gravity anomalies minus the global field contribution and gravimetric terrain effects). The EGM96 geopotential model values are computed in a grid (GRAVSOFT program "harmexp"), and observed free-air anomaly gravity data reduced using this linear interpolation in this grid using ("geoip"). Terrain effects have been removed in a consistent residual terrain model (RTM) data reduction using prisms ("tc" program), taking into account the topographic irregularities relative to a mean height surface, for details see Forsberg (1984). Table 2 shows the effects of including the 100 m DEM data, compared to the 1 km NGDC DEM. The difference is only significant in Hatta. 
Two different resolutions of the RTM mean height surface was tested: 50 km and 100 km (computed by "tcgrid"). In addition a complete topographic-isostatic reduction (land only) was also tested, to see if the terrain- and EGM96-reduced data in Hatta could be sufficiently detrended or bias reduced. This was not possible, and it appears that the small Hatta enclave sits by chance on a major geological gravity anomaly, and that the topographic-isostatic gravity field do not provide a good fit to anomalies. The statistics of the terrain reductions are shown in Table 3, for the DM data sets as well as the complete data set, including marine gravimetry and altimetry. An atmospheric correction has been applied to the gravimetry data, in accordance with theoretical requirements for solving the geodetic boundary values problem. 
The 30' RTM reduction was selected for the subsequent geoid processing, as it fits EGM96 best concerning data bias in Hatta, and should (in theory) be consistent with a perfect EGM. 3.2. Gravity to geoid conversion by FFT Residual gravity anomalies are subsequently converted into residual height anomalies by spherical FFT (Strang van Hees, 1990). The method in principle evaluates Molodensky's integral 
where g1c is the first term of the Molodensky series for the terrain-reduced field, and z the quasi-geoid. The integral is in practice identical to Stokes' integral, as the g1c-term is very small, and may for most practical purposes be neglected. For the computations here we have treated z and the classical geoid N as identical; the difference is less than 2 cm in Hatta, and since the leveling data have not been adjusted in geopotential numbers, the theoretical type of geoid is not an issue, especially considering that GPS-leveling is used to constrain the geoid in the end. In the used FFT method, the Molodensky/Stokes' integral is written as a spherical convolution in latitude and longitude (f,?) for a given reference parallel fref, and by utilization of a number of bands a virtually exact convolution expression may be obtained by a suitable linear combination of the bands. For each band the convolution expressions are evaluated by 
where Sref is a (modified) Stokes' kernel function, and * and F the two-dimensional convolution and Fourier transform, respectively, for details see Forsberg and Sideris (1993). In the actual implementation of the method, the data are girded by least-squares collocation ("geogrid"), and a 100% zero padding is used to limit the periodicity errors of FFT ("spfour"). A Wong-Gore kernel modification has been used for spherical harmonic degrees less than 60 to limit the long-wavelength errors. The collocation gridding was done using a correlation length of 20 km assuming a free-air anomaly noise of 0.5 mgal. Fig. 2 shows the empirical covariance models determined from the reduced data. The gravimetric geoid has been computed in the region 22.5-27.5 N, 52.5 - 58.5 W at a 1' resolution in latitude and longitude. Use of a more detailed grid spacing do not improve the geoid since the residual geoid signal at very short wavelengths is essentially nil.  Fig. 2. Covariance function of reduced gravity data (dashed line isostatic). 3.3. Restoring of topography and EGM96 for the final gravimetric geoid The topographic RTM restore signal (?2) is evaluated by FFT methods, using the first-order mass-layer approximation to the RTM geoid effect, cf. Forsberg, 1985 ("tcfour"). Table 4 shows the statistics for the restore steps. Table 5 shows a comparison to DM GPS-leveling, using a slightly edited data set where obvious outliers greater than 0.5 m were removed. It is seen that the gravimetric geoid indeed makes a large improvement over EGM96 (the mean value is not of significance, as the GPS geoid heights refers to a local datum).
|