Practical Local Geoid Modelling of Benin City, Nigeria from Gravimetric Observations Using the Modified Stokes Integral

The geometric heights obtained from GNSS observations cannot be used for engineering works as they are not reduced to the geoid. This study presents practical local geoid modelling from gravimetric observations using the modified Stokes integral for engineering applications in Benin City. A total of 52 points were observed with GNSS receivers and a gravimeter to respectively obtain their positions and absolute gravity values. The theoretical gravity values of the points were computed on the Clarke 1880 ellipsoid to obtain their local gravity anomalies. The modified Stokes integral was applied to compute the geoid heights of the points. The combined topographic effect was applied to the computed geoid heights of the points to obtain their precise geoid heights. The mean of the precise geoid heights of the points was computed to obtain the local gravimetric geoid model of the study area. The determined geoid model was validated for its reliability as well as the accuracy using the RMSE index. It is recommended that the use of assumed, as well as handheld GPS receiver heights for engineering works should be totally abolished as this study has established the local geoid model of Benin City.


INTRODUCTION
The geoid is an equipotential (level) surface of the earth's gravity field which coincides with mean sea level (MSL) in the open oceans. As such, the geoid provides a meaningful reference frame for defining heights. The importance of accurately modelling the geoid has increased in recent years with the advent of satellite positioning systems such as the Global Navigation Satellite System (GNSS). GNSS provides height information relative to a best-fitting earth ellipsoid rather than the geoid (Seager et al., 1999;Yilmaz and Arslan, 2006). To convert ellipsoidal heights derived from GNSS to conventional (and meaningful) orthometric heights, the relationship between the geoid and the ellipsoid must be known (Kotsakis and Sideris, 1999;Yilmaz and Arslan, 2006). The fundamental relationship that ties ellipsoidal heights obtained from Global Navigation Satellite System (GNSS) measurements and heights with respect to a vertical datum established using spirit levelling and gravity data, is to the first approximation given by (Heiskanen and Moritz, 1967; Krynski and Lyszkowicz, 2006; as: Where h is the ellipsoidal height, H is the orthometric height, and N is the geoid undulation. Figure 1 shows the relationship between the orthometric, geoid and ellipsoidal heights.

Fig: 1: Relationship between Orthometric, Geoid and Ellipsoidal Heights
Source:  Geoid Modelling is carried out by various methods such as the gravimetric, geometric, Astro-geodetic, transformation and gravimetric-geometric methods. The gravimetric method can be carried out by the application of the well-known Stokes-integral, equation (2) and the use of accurately determined absolute gravity data (Heiskanen and Moritz, 1967, Eteje, 2015 and. Where N is geoid undulation, g  is gravity anomaly,    S is Stokes function,  is normal gravity on the reference ellipsoid and R is mean radius of the earth. The application of the Stokes integral for local geoid modelling requires the solving/integration of the modified Stokes integral given in equation (2). Here, the geoid height of each point in a study area is computed and corrected for the combined topographic effect to obtain a precise geoid height of the point.

Integration of Stokes's Formula
According to , using the modified Stokes integral given in equation (2), the geoid heights of points can be computed if their gravity anomalies and geographic coordinates are known. Featherstone and Olliver (1997) gave the integration of equation (2), as well as the Stokes integral as the mean radius of the earth, r = R of the points.

Surface Spherical Radius Computation
The surface spherical radius, o  is computed as (Shrivastava et al., 2015) ) cos( cos cos sin sin cos

Theoretical Gravity Computation
To obtain the local gravity anomalies of points in a study area, the normal, as well as the latitude gravity, is computed on a specified ellipsoid. That is, the ellipsoid adopted for geodetic computation in the area or region of study. Where,

Gravity Anomaly Computation
The gravity anomaly, , g  is the difference between the observed gravity value (g) reduced to the geoid, and a normal, or theoretical, computed gravity value ( o  ) at the mean earth ellipsoid, where, the actual gravity potential on the geoid equal the normal gravity potential at the ellipsoid, at the projection of the same terrain point on the geoid and the ellipsoid respectively, that is (Dawod, 1998 Considering the nature of the topography of the earth surface, which is irregular in shape, there are two basic types of gravity anomalies (free air and Bouguer anomalies). In this study, it was only the free air correction that was applied.

Free Air Correction
This is the first step for reducing topography effects. It simply corrects for the change in the elevation of the gravity meter, considering only air (hence a free-air) being between the meter and selected datum. According to Aziz et al. (2010), this correction is added to the observed gravity because the increased radial distance of the station from the centre of the Earth results in a lower observed gravity value than if the station were at the local datum. The formula to calculate the magnitude of the reduction in practice is given by Eteje Where, H = Station orthometric height in metres g = Mean value of gravity (980500 mGal) r = Mean radius of the Earth

Mean Radius of the Earth Computation
The mean radius of the earth, r = R was computed using: Where M is the radius of curvature along the meridian section and N is the radius of curvature in prime vertical. The formula for computation of the radius of curvature in prime vertical, N is given as (Ono, 2009) (9) while that for computation of the radius of curvature in meridian section, M is given as (Kotsakis, 2008

Combined Topographic Effect Computation
To obtain a precise geoid height of a point, the combined topographic effect is calculated and applied to the computed geoid height of the point. The formula for the computation of the combined topographic effect, Where G is the earth gravitational constant,  is density, R is the mean radius of the earth and H is the orthometric height of observation point which can be obtained from the DTM of the area.

The Geoid Model
The final gravimetric geoid model is the mean of the geoid heights computed with equations (3) and (11)

II. METHODOLOGY
The adopted methodology was divided into different stages of data acquisition, data processing, and results presentation and analysis. Figure 3 shows the adopted methodology flow chart.

Data Acquisition
A total of 52 points were used in the study. The points included two primary control stations (XSU 92 and XSU 100 were respectively located in Edo College and School of Nursing premises). The other 50 points were selected along the major roads of the City (See Figure 4). Spirit levelling was carried out on 3 of the 50 points for validation purpose. GNSS observations were carried out using CHC 900 dual-frequency GNSS receivers to obtain the coordinates and ellipsoidal heights of the points. The observations were carried out relative to control station XSU 92 using the static method (See Figures 5 and 6). Selected Points (RR01) at Ring Road The selected points were observed with a gravimeter (SCINTREX CG-5 Autograv Gravimeter) to obtain their absolute gravity values. The observations were carried by an expert, a Geophysicist from Mountain Top University, Ibafo, Ogun State. The gravity observations of the points were carried out in seven different loops relative to a point whose absolute gravity value was known and located within the Benin City Airport premises (See Figures 7 and 8).   Figure 9 shows the validation points' levelling loops.

Data Processing
The GNSS observations were respectively downloaded and processed with HcLoader and Compass Postprocessing software to obtain the positions and the ellipsoidal heights of the points. The coordinates and the ellipsoidal heights of the points were processed in Minna datum. The gravity observations of the points were processed by the expert who carried out the observation to obtain their absolute gravity values. All the necessary corrections such as drift correction, etc were applied during the processing. The theoretical gravity values of the points were computed on the local (Minna) datum ellipsoid (Clarke 1880 ellipsoid) using the latitude coordinates of the points, as well as equation (5). The gravity anomalies of the points were computed by finding the differences between the absolute gravity values of the points and their corresponding theoretical gravity values, as well as using equation (6). The computation of the free air correction requires the application of the orthometric heights of the points. And these were obtained by interpolation using the orthometric heights and the absolute gravity values of the two primary control stations (XSU 100 and XSU 96). The orthometric heights of the points were interpolated as there was no Digital Terrain Model (DTM) of the study area. The free air correction was applied to the computed gravity anomalies of the points using equation (7). The free air and the Bouguer gravity anomalies of the points were computed but the free air gravity anomalies were used in the study. This is because the geoid heights of the two primary control stations obtained from their known orthometric and ellipsoidal heights approximated the geoid heights of the stations computed using the free air gravity anomalies, as well as equation (3). The gravimetric geoid heights of the points were computed with the geographic coordinates, free air gravity anomalies and the theoretical gravity of the points using equation (3). The computation of the gravimetric geoid heights of the points required the application of the spherical radius and the mean radius of the earth. The spherical radius and the mean radius of the earth were respectively computed using equations (4) and (8). Also, the computation of the mean radius of the earth required the computation of the radius of curvature in prime vertical and in meridian section using equations (9) and (10) respectively. The computed gravimetric geoid heights of the points using equation (3) were co-geoid heights. To obtain precise gravimetric geoid heights of the points, the combined topographic effect has to be computed and applied to the co-geoid heights. The combined topographic effect was computed using equation (11). The final local gravimetric geoid model of the study area was obtained by finding the mean of the precise geoid heights of the points using equation (12). The spirit levelling of the validation points was reduced using the height of instrument method, as well as collimation method. The RMSE, as well as the accuracy of the model, was computed by finding the square root of the mean of the differences between the model and the observed orthometric heights of the validation points and the two control stations using equations (13).

III. RESULTS PRESENTATION AND ANALYSIS 3.1 Analysis of the GNSS Observation Results
The DGPS observations were carried out to obtain the coordinates and ellipsoidal heights of the selected points. The DGPS observations were processed using Compass post-processing software. From the processing of the DGPS observations results, it was seen that the processed observations passed both the Network Adjustment Test and the X-Square (Chi-square) Test. This implied that the normal matrix generated was a regular one and inverted accordingly for the calculation of residuals.  /dx.doi.org/10.22161/ijaems.512.1  ISSN: 2454-1311 Table 1 presents the closure errors/accuracy of the two loops of the validation points levelling. The levelling of the validation points was done to obtain the orthometric heights of the validation points. The levelling was carried in two loops. The first loop started from XSU100 to TP1 and closed back on XSU100 while the second loop started from TP1 through TP2 to TP3 and closed back on TP1.

Analysis of the Validation Points Levelling
From Table 1 it is seen that the closure error of the first loop is -0.006m while the second loop closure error is 0.009m which were within millimetres standard. The results were accepted as the closing errors of the two loops. The high accuracy of the levelling was as a result of the fairly flat topography of the study area, the observer's know-how and the equipment used.  Table 2 and Figures 10a and 10b respectively present the coordinates of the selected points, their corresponding computed local gravimetric geoid heights and the determined geoid model, and their surface and contour plots. The gravimetric geoid heights of the points were computed using the integration of modified Stokes' integral. The gravity anomalies of the points used for the computation of the local gravimetric geoid heights of the points were computed on the Clarke 1880 ellipsoid which is the ellipsoid adopted for local geodetic computation in Nigeri a. From Table 2, it is seen that the determined geoid model (mean of geoid heights) is 2.066m. This implied that to convert ellipsoidal heights of points to orthometric heights in Benin City, 2.066m will be subtracted from the ellipsoidal heights of the points. Also, from Figures 10a and 10b, it is seen that at the centre, there are depressions with small cliffs closed to them while at some distances from the centre there are small depressions which implies that the geoid heights of the study area vary with no constant value.  done to present the consistency, as well as the reliability of the determined geoid model. It can be seen from Table  3 that the computed RMSE, as well as the accuracy of the determined gravimetric geoid model, is 0.412m. This implies that ellipsoidal heights can be converted to orthometric height with an accuracy of 41cm using the determined geoid model. It can also be seen from Figure  11 that the model and the observed orthometric heights of the validation points/stations are identical in shape which also implies the high reliability, as well as consistency of the determined gravimetric geoid model.

IV. CONCLUSIONS AND RECOMMENDATIONS
1. The local geoid model of Benin City has been determined to be 2.066m using the gravimetric method of integration of modified Stokes formula. 2. The study has shown that ellipsoidal heights can be converted to orthometric heights with an accuracy of 0.412m using the determined geoid model. 3. It is recommended that the determined geoid model should be applied whenever ellipsoidal heights are to be converted to practical, as well as orthometric heights in Benin City. 4. It is also recommended that the use of assumed, as well as handheld GPS receiver heights should be totally abolished as this study has established the local geoid model of Benin City.