Kurotamuno Peace Jackson, Lawrence Hart
Department of Surveying and Geomatics, Rivers State University, Port Harcourt, Nigeria.
DOI: 10.4236/ijg.2026.171004   PDF    HTML   XML   104 Downloads   389 Views   Citations

Abstract

Orthometric heights are obtained by Spirit leveling; however, this approach is not only slow but is also costly and time-consuming. Therefore, the uncertainties in the accuracy of geoid height models have hampered the use of GNSS heightening as an alternative method of determining Orthometric heights across Nigeria. This study was aimed at appraising the status of the geoid model in Nigeria based on gravity measurements, applying the remove-compute-restore (r-c-r) approach through the Fast Fourier Transform (FFT) tool. The GRAVSOFT tool was used in evaluating the residual height anomalies and terrain effects resulting from the residual terrain model (RTM). The residual height anomalies, (ζ_res) were computed by applying the modified Molodenskii Integral tampered 100% zero padding so as to overcome cyclic effects. The residual height anomalies were converted to geoid undulation by adding the ζ-to-N correction. The study revealed that, the difference between the height anomaly (ζ) and the geoid height (N) is about 24 cm and the geoid values across Nigeria range between 13.738 m and 26.052 m. There is a tilt of the geoid from the South-west to North-west part of the country. The gravimetric geoid is accurate to ±41 cm rmse when compared with the geometric geoid while the surface fitting revealed that the hybrid-geoid is accurate to ±3 cm rmse.

Keywords

Fast Fourier Transform, Gravimetric Geoid, Height Anomalies, Hybrid-Geoid, Residual Terrain Model

Share and Cite:

1. Introduction

The quest for an improved local geoid for Nigeria has been on for some time now. Geoid determination or modeling is one of the essential problems in geodesy. A key aim of physical geodesy is the determination of the geoid with accuracy at the centimeter (cm) level, matching the accuracy of Global Navigation Satellite System (GNSS) height determination. The combination of an accurate geoid model with GNSS coordinates plays a central part in achieving timely, high accuracy leveling results, mostly over long distances, unlike Spirit leveling, which is tedious, costly, and time-consuming. However, the up-to-date accuracy of global geoid modeling is a few decimeters to meters, which is not adequate for many scientific and engineering applications. Thus, there is a need to model the geoid locally.

The geoid heights, when compared to the ellipsoidal heights from GNSS, can result in orthometric heights accurate to those delivered by spiriting leveling. For practical purposes, most countries have established an officially fitted geoid model corresponding to the height system(s) in use. Such a national model can be used for height system unification, scientific purposes, and other regional tasks, but also as a national reference surface. One of the primary practical applications of geoid height (N) in surveying is to transform GNSS-derived ellipsoidal heights (h) to orthometric heights (H) as argued [1] using the simple algebraic relation as shown in Eq. 1:

$H=h-N$ (1)

However, there is no officially published geoid model in Nigeria for engineering and related applications. More so today, local modeling of the geoid, amongst other things, facilitates the connection between land datum and sea datum; this has been implemented in most developed countries with the terminologies such as Vdatum (Canada and US) and Vertical Offshore Reference Frame (VORF) in Great Britain, Ireland, and Northern Ireland [2]. It becomes very pertinent to emphasize the necessity for the determination of an acceptable geoid model for Nigeria, which will help and drive proper integration and effective use of the Global Navigation Satellite System (GNSS) over Nigeria for GNSS heightening. The geoid is evaluated using different methods such as geometric/leveling, astro-gravimetric, satellite, gravimetric, and a combination of any of the above methods. Each of these approaches has its peculiar strengths and weaknesses.

Thus, high-resolution regional/local models are still required for most practical purposes. There are many different methods proposed over the years for global, regional, and local geoid evaluation by gravimetric data, each with its own procedure, concepts, and philosophy. In addition, today all of these techniques integrate long-wavelength Global Geopotential Models (GGMs) with a combination of local terrestrial and airborne gravity data. However, they differ purely in the way they join these data sets.

There are different approaches to gravimetric geoid modeling such as Least Squares Collocation (LSC), Least Squares Modification of Stokes Integral with Additive Corrections, Remove-Compute-Restore (R-C-R) approach using Fast Fourier Transform (FFT) Technique, etc., each based on its philosophy, weaknesses, and strengths. In the gravimetric approach to geoid modeling, the geoid is modeled using the Stokes integral after gravity reduction as shown in Eq. 2 [3].

$N=\frac{R}{4\pi {\gamma }_0}{\displaystyle \underset{\sigma }{\iint }\Delta gS\left(\psi \right)\text{d}\sigma }$ (2)

where: $N$ : The Geoid height, ${\gamma }_0$ : Normal gravity, $\Delta g$ : Gravity anomaly, and$S\left(\psi \right)$ : Stokes kernel. The Stokes integral has two conditions that must be satisfied for the results to be true. Firstly, as the solution to the geodetic boundary-value problem (BVP), the Stokes integral requires the gravity anomalies to represent the boundary values on the geoid. That is, it must be referred to the geoid.

Secondly, there must be no masses present outside the geoid, making sure that 𝑔 on the geoid truly is the boundary surface, as the BVP’s used in physical geodesy are the solution to Laplace’s Equation, and Δ𝑉 ≠ 0 if masses exist above surface 𝑆 as opined by [3]. The gravity on the geoid and the gravity on the reference ellipsoid must be well-defined, such that there are no masses present above the geoid, to satisfy the requirements of the boundary value problems. Since gravity cannot be measured on the geoid, gravity observed on or above the Earth’s surface must be reduced to the geoid.

Even when gravity is measured at the geoid, all the mass on Earth above the geoid would have to be removed. This modifies the mass distribution of the Earth, making it difficult to define the actual geoid. The second condition of the Stokes integral; is overcome by applying the theory of Molodenskii. The terrain replaces the geoid as the boundary surface, and the result of the Stokes integral is the quasi-geoid instead of the geoid. While the quasi-geoid is approximately equivalent to the geoid in terms of shape, unlike the geoid, it has no physical meaning and is merely a mathematical surface of convenience [4]. In the remove-compute-restore (R-C-R) approach, as posited by [1], the long, medium, and short wavelength components are derived from global geopotential models, terrestrial gravity, and Digital Terrain Models (DTM). This approach adopts the procedure of computing the gravity anomalies and then the geoid model, taking into account the integration of the different wavelengths mentioned earlier.

The study is aimed at appraising the status of the geoid in Nigeria with the following objectives:

1) Compute the residual gravity anomaly ($\Delta g_{RES}$ ) from terrestrial gravity, including residual terrain models for the remove and restore effects,

2) Compute residual height anomalies from Molodenskii integral based on the FFT model and sum the resulting different effects to obtain height anomalies and also convert the height anomalies to classical geoid values so as to generate the regular geoid grid file based on FFT model,

3) Compare the gravimetric and geometric geoid and tailor the geoid file to the GNSS/leveling data using a seven-parameter model.

2. Data and Methods

2.1. Data Used for Research

The data used for this research are terrestrial gravity data, digital terrain model (DTM), and GNSS/leveling and were obtained from different sources. The terrestrial gravity dataset was obtained from Bureau Gravimetrique International [5] as shown in Table 1. These data were contributed by different organizations and individuals either for oil prospecting or for the development of a gravity geoid for the country, however, the data were domiciled at Bureau Gravimetrique International (BGI), and the fill-in gravity data were obtained from the International Centre for Global Earth Models [6]. The gravity field information received contains the absolute gravity value (g), free-air anomaly, bouguer anomaly and the height of the gravity stations. Statistically, more than ninety-five (95) percent of the gravity field information received contains this information.

Table 1. Statistics of gravity and anomalies data used for the research.

Types	Counts	Max. (mGal)	Min. (mGal)	Mean (mGal)	Std. Dev.	Remarks
Gravity	4,234	978203.44	977907.71	978057.46	36.16	BGI
Free-Air Anomaly	4,234	71.406	−70.960	−1.910	24.652	GRS80
Bouguer Anomaly	4,234	55.970	−85.260	−16.983	22.306	GRS80
Fill-in Dataset
Gravity	8,167	978286.21	977692.60	978073.84	128.22	EGM08
Free-Air Anomaly	8,167	72.60	−70.96	3.59	21.11	GRS80
Bouguer Anomaly	8,167	19.80	−105.05	−44.00	22.00	GRS80

Therefore, it was easier to evaluate the inner consistency between the information in the data given such as the height, free-air anomaly and the Bouguer anomaly, and to identify which normal gravity formula has been used (i.e. GRS67 or GRS80). Absence of internal consistency between constituents of a data record is a clear indication of erroneous gravity data. As it is most often not clear which constituent of a data record is erroneous, the internally non-consistent records could not be repaired and were therefore simply eliminated from the dataset. Nearly, five percent of data points were removed from the dataset, and it was discovered that the GRS67 ellipsoid was used in the computation, which was subsequently converted to align with the global ellipsoid GRS80.

The digital terrain model (DTM) sourced from ALOS Global Digital Surface Model “ALOS World 3D-30 m (AW3D30)” as reported by [7]. The SELECT subroutine was used to re-sample the DTM by low-pass filtering employing a moving average to create the detailed, coarse and reference grids. The data sets used for this study are shown pictorially in Figure 1.

Figure 1. Datasets used for this study.

2.2. Methodology

The remove-compute-restore (R-C-R) method was used in this study. According to [1], this method takes into consideration the geoid’s long, medium, and short wavelength components, which are obtained via Digital Terrain Models (DTM), terrestrial gravity, and global geopotential models. This method uses the previously described process of calculating the gravity anomalies and subsequently the geoid model, while accounting for the integration of the various wavelengths. Four main processes are involved in modeling the gravimetric geoid and fitting it to the geometric geoid model: i.e., removal, prediction, restoration, and tailoring/fitting.

2.2.1. Gravity Anomalies from Gravity Data

There are different methods proposed over the years for gravimetric geoid modeling. These methods depend on the approach of gravity reduction such as simple or complete Bouguer anomalies, Faye anomalies and free-air anomalies etc. The methods of the treatment of the topographic masses are distinguished amongst the various techniques of gravity anomalies. In this study the basic free-air anomalies were first employed in the preliminary stage. Nevertheless, the Bouguer anomalies were employed in the conversion of the height anomalies (ζ) to the classical geoid (N). The estimates Bouguer, $\Delta g_B$ , and simple Free air, $\Delta g_{FA}$ anomalies are given in Eqs. (3) and (5).

$\Delta g_{FA}=g+\delta y_A-{\gamma }_0+0.3086H$ (3)

where: $g$ : The observed terrestrial gravity., $\delta y_A$ : The atmospheric correction after [8], ${\gamma }_0$ : The normal gravity of a point on the ellipsoid (GRS80), $0.3086H$ : The simple Free-air gradient multiplied by the height of the gravity station above the geoid.

And the normal gravity (${\gamma }_0$ ) was evaluated numerically after Somigliana as:

${\gamma }_0={\gamma }_a\frac{1+k{\sin }^2\phi }{{\left(1-e^2{\sin }^2\phi \right)}^{1/2}}\text{wit}hk=\frac{b{\gamma }_b}{a{\gamma }_a}-1$ (4)

Where: ${\gamma }_0$ : The normal gravity on the ellipsoid at a point $\phi$ , ${\gamma }_a$ : The normal gravity at the pole, ${\gamma }_b$ : The normal gravity at the equator, $a,b$ : The semi-major and minor axes of the reference ellipsoid. The simple Bouguer anomalies were computed using Eq. (5):

$\Delta g_B=g-A_B+F-{\gamma }_0=g-0.1119H+F-{\gamma }_0$ (5)

where: $\Delta g_B$ : The simple Bouguer anomaly, $g$ : Observed gravity, $0.1119H$ : The Bouguer plate correction multiple by the height of the gravity station, $F$ : Normal vertical gravity gradient = 0.3086H, ${\gamma }_0$ : Normal gravity of a point on the ellipsoid (GRS80).

2.2.2. Evaluation of the Gravity Anomaly and Height Anomaly from Harmonic Coefficients

Global geopotential models (GGM) represent the solution of the boundary value problem for the global approach. In regional gravity field modelling, GGMs are used for computation of global, long-wavelength, contributions of different functionals such as the height anomaly $\zeta$ using [9]:

$\zeta \left(\phi ,\lambda ,r\right)=\frac{GM}{r\gamma }{\displaystyle \sum }_{n-2}^{n_{\max }}{\displaystyle \sum }_{m=0}^n\left({\overline{C}}_{nm}\cos m\lambda +{\overline{S}}_{nm}\cos m\lambda \right){\overline{P}}_{nm}\left(\sin \phi \right)$ (6)

where: ${\overline{C}}_{nm}$ and ${\overline{S}}_{nm}$ are the fully normalized geopotential coefficients of the anomalous potential, ${\overline{P}}_{nm}$ are the normalized Legendre functions, $n_{\max }$ is the maximum degree of expansion. $r=R+H$ . R is the mean radius of the Earth and H the height of the evaluation point. Gravity anomalies are computed using the equation:

$\Delta g_{GGM}\left(\phi ,\lambda ,r\right)=\frac{GM}{r^2}{\displaystyle \sum }_{n-2}^{n_{\max }}\left(n-1\right){\displaystyle \sum }_{m=0}^n\left({\overline{C}}_{nm}\cos m\lambda +{\overline{S}}_{nm}\cos m\lambda \right){\overline{P}}_{nm}\left(\sin \phi \right)$ (7)

The Geoid Undulation N is obtained by:

$N_{GGM}\left(\phi ,\lambda \right)=\frac{T-2\pi G\rho H^2}{\gamma }$ (8)

In a similar manner the height anomaly is obtained as:

${\zeta }_{GGM}=\frac{T-\text{Δ}{g}_{GGM}H}{\gamma }$ (9)

Where: H is the height of the evaluation points. In this research EGM 2008 was used as reference to a maximum degree of 2190.

2.2.3. Evaluation of “Remove” Terrain Effects

The RTM approach was chosen for the computation of the terrain effects because it yields an optimal result and the residual gravity anomaly and indirect effect are smaller and smoother than with other methods, according to analysis done by various scholars in geoid modeling [10]. Additionally, it demonstrates that RTM reduction is one of the most accurate methods of gravity reduction in gravimetric geoid computation. In this study, the Residual Terrain Model (RTM) was used to compute the terrain effects by prism integration. One of the most used mass reduction techniques, mostly for quasi-geoid determination, is the RTM. This approach eliminates the topography’s input and replaces it with an equivalent topography model.

As a result, masses fill up the gaps below the reference surface and the topographic masses above it are eliminated. The fine (detailed) resolution topography grid was averaged, and the average grid created by taking 9’ × 9’ moving averages over the coarse 3’ × 3’ mean height grid was then low pass filtered to create the reference surface. Geoid terrain impacts by prism integration were calculated using the FFT approach. In practice, the coarse topographic grid is utilized for the entire topography of the area being evaluated, while the detailed topographic grid is utilized up to a certain distance [11]. The resolution of the available detailed DTM determines the inner zone’s radius. A grid with a lesser resolution can be utilized for points outside of this inner region.

This approach is fully implemented in the GRAVSOFT software, particularly inside the TCFOUR program [11] and [13], which is frequently used in applications of gravity field modeling. Considering the coarse/detailed grid system speeds up computation. Using this program, a bicubic spline interpolation approach is used to further densify the topographic data in a small inner zone surrounding the computation point. This allows for the integration of the inner zone’s frequently significant effects [11]. The RTM terrain was directly evaluated by the TCFOUR program by Fourier (FFT) method.

Following [12] the topographic effect on gravity of the RTM reduction is computed using Eq. (10):

$\text{Δ}{g}_{\text{RTM}}=G\rho {\displaystyle \underset{}{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{}{\int }}{\displaystyle \underset{z=h_{ref}\left(x,y\right)}{\overset{z=h\left(x,y\right)}{\int }}\frac{z-h_p}{{\left[{\left(x_Q-x_P\right)}^2+{\left(y_Q-y_P\right)}^2+{\left(z_Q-h_P\right)}^2\right]}^{3/2}}\text{d}{x}_Q\text{d}{y}_Q\text{d}{z}_Q}}}$ (10)

2.2.4. Evaluation of “Restore” Terrain Effects

The restore step follows after residual quasi-geoid ${\zeta }_{\text{Δ}g}$ has been obtained in the compute step. The GGM and RTM contributions on the height anomalies ${\zeta }_{GGM}$ and ${\zeta }_{RTM}$ are computed using the same input parameters that were used in the remove step for computation of $\text{Δ}{g}_{GGM}$ and $\Delta g_{RTM}$ contributions on gravity anomalies. The height anomaly, ${\zeta }_{GGM}$ from the GGM and height anomaly, ${\zeta }_{RTM}$ from the DTM are then added to the computed residual quasi-geoid ${\zeta }_{res}$ . The RTM terrain effects on the quasi-geoid is evaluated using Eq. (10) but evaluated in Eq. (11) by FFT, here reproduced in short form:

${\zeta }_{restore}={\zeta }_{RTM}=\frac{G\rho }{\gamma }{\displaystyle \underset{h_{ref}}{\overset{h}{\int }}\frac{1}{r}\text{d}x\text{d}y\text{d}z}$ (11)

where $r={\left[{\left(x_P-x\right)}^2+{\left(y_P-y\right)}^2+{\left(H_P-z\right)}^2\right]}^{1/2}$

2.2.5. Evaluation of Reference Fields from the Topography

The evaluation of the terrain effects, especially by the RTM approach, required coarse, detailed and a reference grid. The purpose of the reference grid is to provide optimal smoothing of the geoid effects. Thus, a reference grid height around 100km resolution was computed using the TCGRID subroutine. The TCGRID does two operations: it first makes an average grid and then filters the average grid with a moving average operator. The topography reference field was evaluated using TCGRID module in the GRAVSOFT program.

2.2.6. Evaluation of Reduced Gravity Values

The compute step, which determines the residual gravity field by subtracting the global and local gravity effects, comes after the remove step in the gravimetric geoid evaluation using the R-C-R approach. The computation stage uses residual gravity field data as input. The data are often gravity anomalies. Since there are multiple topographic reduction techniques available, there are various variations of the R-C-R methodology. Gridding can be done either before or after the removal step, and the input data can be either Faye or RTM anomalies.

The reduced gravity values were obtained using the following steps:

1) Start from the observations on the Earth’s surface $\Delta g_{FA}$ .

2) Compute long-wavelength contribution to the gravity $\Delta g_{GGM}$ .

3) Compute RTM terrain effects on gravity $\Delta g_{terr.effects}=\Delta g_{RTM}$ .

4) Obtain residual gravity anomalies $\Delta g_{RES}=\Delta g_{FA}-\Delta g_{GGM}-\Delta g_{RTM}$ .

5) Interpolate (grid) residual gravity anomalies $\Delta g_{RES}$ .

2.2.7. Evaluation of the Modified Stokes/Molodenskii Integral

In Stokes integration, each gravity anomaly Δg from the grid within the integration radius is multiplied with the Stokes function$$ S\left(\psi \right)$$. Stoke’s function can be computed spectrally using Legendre polynomial series or analytically by closed expressions. The spectral Stokes’ function was employed in this research as shown as given in Eq. (12):

$\text{Spectral}:S\left(\psi \right)={\displaystyle \sum }_{n=2}^{\infty }\frac{2n+1}{n-1}{P}_n\cos \left(\psi \right)$ (12)

However, the continuous integration of gravity anomalies across the entire sphere should be used to evaluate Stokes’ integral; yet, these values are not available for integration over the entire sphere, particularly in regional gravity field modeling. Because of this, the integration can only be done inside the spherically defined area (${\sigma }_0$ ) and not over the entire sphere (σ). Since local gravity data are known to contain errors in long wavelengths, information is thus lost in the remote zone where data are not integrated ($\sigma -{\sigma }_0$ ), resulting in long wavelength truncation errors. Long wavelength inaccuracies from gravity measurements therefore spread into the geoid under evaluation.

Molodenskii demonstrated that the modified Stokes’ kernel function, which integrates the long wavelength contribution in the form of low-degree coefficients (up to degree L) from satellite-derived GGMs and terrestrial gravity data, may minimize the truncation error of the remote zone. The Stokes kernel function is altered in this study in accordance with [14]. The low degree Legendre polynomials are eliminated from the unaltered Stokes’ function by the Wong and Gore modification. The height anomaly ($\zeta$ ), is evaluated as given in Eq. (13):

${\zeta }_{res}=\frac{R}{4\pi \gamma }{\displaystyle \underset{\sigma }{\iint }\left(\text{Δ}{g}_{res}+G_1\right)S\left(\psi \right)\text{d}\sigma }$ (13)

where: $G_1=\frac{R^2}{2\pi }{\displaystyle \iint \frac{\left(H-H_P\right)}{l_0^3}\Delta g\text{d}\sigma }$

The quasi-geoid is derived by summing the effects from the different gravity field as shown in Eq. 14:

$\zeta ={\zeta }_{res}+{\zeta }_{RTM}+{\zeta }_{GGM}$ (14)

The quasi-geoid is converted to the classical geoid using the relation in Eq. 15:

$N=\zeta -\frac{\text{Δ}{g}_{B}H}{\gamma }$ (15)

where $H_P$ : The evaluation point and H being the reference point.

3. Results and Discussion

The statistical results of the residual gravity anomaly, residual terrain effects are given on Table 2.

Table 2. Statistics of the RTM terrain effects (Remove effects).

Gravity Field/Unit: mGal	Max.	Min.	Mean	Std. Dev.	Remarks
$\Delta g_{FA}$	72.599	−70.960	3.592	21.111	1
$\Delta g_{GGM}$	58.292	−61.591	6.963	17.326	2
$\Delta g_{FA}-\Delta g_{GGM}$	71.245	−92.356	−2.190	14.719	3
$\Delta g_{RTM}$	140.445	−154.540	−0.079	9.971	4
$\Delta g_{FA}-\Delta g_{GGM}-\Delta g_{RTM}$	131.766	−92.318	−0.851	14.243	5
RTM ${\zeta }_{RTM}$	1.999	−0.485	0.033	0.173	7
RTM Geoid Effects ${\zeta }_{RTM}$ (Niger Republic)	1.473	−0.306	−0.014	0.069	8

The RTM impacts should have a mean value and standard deviation that are nearly zero, according to the statistics in Table 2. The RTM effects on the geoid in this study are smooth because, as Table 2 illustrates, the mean value is close to zero. The acquired results are in line with the expected RTM anomalous property, which has a small indirect influence and the smallest power at short wavelengths. Furthermore, using the same methodology, [15] findings from the neighboring Niger Republic produced results that were similar to those shown in Table 2.

Additionally, it is preferable that the remaining gravity field quantity following the remove step be smooth and free of outliers or systematic biases in the data before gridding. The residual gravity anomalies are smooth, according to the data in Table 2, and they are consistent with readings from neighboring countries that used the same methodology.

The results of the residual gravity anomaly $\Delta g_{RES}$ was gridded to a 3′ x 3′ grid using GEOGRID subroutine and this was input into the SPFOUR (GRAVSOFT) for the evaluation of the residual height anomaly ${\zeta }_{RES}$ after Molodenskii modified after Wong-Gore kernel modification. The statistic of the output is shown on Table 3.

Table 3. Procedure and statistics of the reduced height anomaly.

Gravity Field/Unit:(m)	Max.	Min.	Mean	Std. Dev.
FFT Contribution ζres	2.060	-2.195	-0.004	0.446

In this study, the geoid, as defined by Eqs. 14 and 15, was initially assessed as a quasi-geoid and then transformed into a geoid by applying the quasi-geoid to the geoid correction (ζ-to-N). Mountainous regions, including the Plateau and Taraba States in North Central and East of the research area, have the highest values of the (ζ-to-N) in absolute terms. These modifications have maximum range values greater than 24 cm. This modification is minor and occasionally disregarded in low-lying and coastal regions. Although it hasn’t been accomplished yet, geodesists worldwide are working to create sub-centimeter geoid. The total accuracy of the geoid was evaluated for this study using GNSS/leveling, whose dependability is unknown. An estimate of the projected error budget of the geoid might have been obtained by performing an error propagation analysis to assess the undulations’ error if the GNSS/leveling error had been known.

But given this clear fact, we assumed that the GNSS/leveling data to be errorless. As a result, Table 4 displays the comparison and overall level accuracy of the geoid for this study, which is 41 cm (decimeters level). The Nigerian Primary Geodetic Network’s GNSS/leveling data, which mostly comes from the orthometric heights from the geodetic leveling network, is the main source of the significant disparity between the geometric and gravimetric geoids. Furthermore, unmodeled errors from the DEM utilized for this study, terrestrial gravity data, and the interpolation of the reduced gravity anomalies used in the Stokes or Molodenskii integral may potentially be responsible for the decimeter level precision attained.

Table 4. External assessment of gravimetric geoid with GNSS/levelling.

ID		S/N		Lat. dd		Long. dd		h (m)		H (m)		NGNSS (m)		NGRAV (m)	NGNSS-NGRAV (m)
A16		1.		10.1249		12.3759		786.6040		768.7680		17.8360		18.0830	0.2470
A39		2.		11.2889		10.4174		495.4990		474.7460		20.7530		21.3500	0.5970
C16		3.		6.1371		9.0267		628.4300		607.4730		20.9570		21.1740	0.2170
CFL56		4.		11.8532		13.1164		357.4050		340.6800		16.7250		17.1610	0.4360
CFH66		5.		6.1731		6.7501		57.3680		37.2140		20.1540		19.7030	−0.4510
CFA33A		6.		6.6269		3.3231		70.4120		47.0930		23.3190		23.7260	0.4070
H4	7.		7.4627		8.6032		362.9030		341.7780		21.1250		20.5410			−0.5840
L10	8.		7.2044		3.3449		197.8920		172.9870		24.9050		25.6210			0.7160
R43	9.		12.0212		5.8946		477.2400		456.4400		20.8000		20.9280			0.1280
U81	10.		6.7870		6.4857		215.8380		193.7200		22.1180		23.0190			0.9010
U78	11.		6.7764		7.2634		424.4120		402.0030		22.4090		23.2000			0.7910
ZVS3003	12.		4.8480		7.0478		35.3270		16.6140		18.7130		18.3650			−0.3480
L8	13.		7.9562		3.6267		464.8940		439.4110		25.4830		25.3410			−0.1420
L18	14.		8.5956		4.3915		458.4220		434.3450		24.0770		23.4510			−0.6260
D13	15.		10.2756		4.3907		290.6730		268.1180		22.5550		22.1010			−0.4540
K01	16.		12.0143		8.5664		504.2310		482.9880		21.2430		20.7170			−0.5260
C014	17.		6.2028		8.6243		148.3030		127.2320		21.0710		21.2820			0.2110
N123A	18.		11.4412		6.9920		769.7730		747.9120		21.8610		21.5390			−0.3220
N032	19.		9.1060		7.2016		708.8040		686.0680		22.7360		22.0860			−0.6500
A001	20.		9.1888		12.4961		217.9380		200.5570		17.3810		18.3160			0.9350
C036	21.		7.9988		10.9930		543.6250		524.4710		19.1540		19.3050			0.1510
A21	22.		10.4553		11.6280		758.9630		740.3310		18.6320		18.7310			0.0990
AP1	23.		4.8695		6.9779		33.7200		14.8081		18.9119		18.8199			−0.0920
HS8	24.		4.7651		7.0166		26.0280		6.9860		19.0420		19.0350			−0.0070
PT 4 Emma	25.		4.7984		7.0056		30.6930		11.6906		19.0024		18.9524			−0.0500
PT 8 Emma	26.		4.8338		7.0070		26.7890		7.8509		18.9381		18.8381			−0.1000
PT 5 Emma	27.		4.8069		7.0094		29.3740		10.3801		18.9939		18.9339			−0.0600
PT 9 Emma	28.		4.8366		7.0153		29.1410		10.1660		18.9750		18.8860			−0.0890
PT 2 Abdul	29.		4.8443		7.0395		32.6400		13.6539		18.5861		18.5161			−0.0700
PT 3 Abdul	30.		4.8408		7.0313		26.7500		7.7697		18.9803		18.9093			−0.0710
AP4	31.		4.8683		6.9899		35.8490		16.9261		18.9229		18.8349			−0.0880
PP5	32.		4.8703		7.1089		38.8020		19.7522		19.0498		19.0028			−0.0470
PT 4 Abdul	33.		4.8372		7.0229		32.8420		13.8392		19.0028		18.9558			−0.0470
PT 3 Emma	34.		4.7902		7.0023		25.1950		6.2283		18.9667		18.8907			−0.0760
PT 6 Emma	35.		4.8155		7.0098		34.5140		15.4366		19.0774		19.1044			0.0270
													RMS			41 cm

Additionally, the data showed that the quasi-geoid had a tilt and/or slope, supporting the conclusions of earlier research conducted in the region by [16] and [17]. As seen in Figure 2(a) and Figure 2(b), this tilt and/or slope extends from the north central and southwestern regions to the research area’s edges. The degree of proximity between the two surfaces was revealed by the fitness between the gravimetric and geometric geoid. The central and southwest regions of the research area had the highest values of the undulations, and the geoid values throughout Nigeria show that the geoid is above the ellipsoid.

(a)

(b)

Figure 2. (a): Quasi-geoid over Nigeria; (b): Geoid over Nigeria.

The degree of agreement between the geometric and gravimetric geoids, that is, (N^GNSS-N^GRAV) is shown pictorially Figure 3. This clearly shows whether the geometric geoid is below and/or above the gravimetric geoid.

Figure 3. The residuals between the geometric and gravimetric geoid over Nigeria.

However, the recent GNSS/leveling dataset acquired from the office of the Surveyor’s General Rivers State was used exclusively to evaluate the geoid file, revealing a centimeter level of accuracy as shown in Table 5 to demonstrate that the modelled geoid in this study and the level of accuracy achieved for the gravimetric geoid across Nigeria is at centimeter level.

Table 5. External assessment of gravimetric geoid with GNSS/levelling (Rivers state).

ID	S/No	Lat. dd	Long. dd	h (m)	H (m)	NGNSS (m)	NGRAV (m)	NGNSS-NGRAV (m)
AP1	1.	4.8695	6.9779	33.7200	14.8081	18.9119	18.8199	−0.0920
HS8	2.	4.7651	7.0166	26.0280	6.9860	19.0420	19.0350	−0.0070
PT 4 Emma	3.	4.7984	7.0056	30.6930	11.6906	19.0024	18.9524	−0.0500
PT 8 Emma	4.	4.8338	7.0070	26.7890	7.8509	18.9381	18.8381	−0.1000
PT 5 Emma	5.	4.8069	7.0094	29.3740	10.3801	18.9939	18.9339	−0.0600
PT 9 Emma	6.	4.8366	7.0153	29.1410	10.1660	18.9750	18.8860	−0.0890
PT 2 Abdul	7.	4.8443	7.0395	32.6400	13.6539	18.5861	18.5161	−0.0700
PT 3 Abdul	8.	4.8408	7.0313	26.7500	7.7697	18.9803	18.9093	−0.0710
AP4	9.	4.8683	6.9899	35.8490	16.9261	18.9229	18.8349	−0.0880
PP5	10.	4.8703	7.1089	38.8020	19.7522	19.0498	19.0028	−0.0470
PT 4 Abdul	11.	4.8372	7.0229	32.8420	13.8392	19.0028	18.9558	−0.0470
PT 3 Emma	12.	4.7902	7.0023	25.1950	6.2283	18.9667	18.8907	−0.0760
PT 6 Emma	13.	4.8155	7.0098	34.5140	15.4366	19.0774	19.1044	0.0270
							RMS	7 cm

Results of the Cross Validation of the GNSS/Levelling Data

The disparity between the derived and observed orthometric heights nationwide was fitted or tailored using a seven-parameter empirical model as shown in Eqs. 16 to 18:

gravimetric geoid can be expressed as:

$l_i=h_i-H_i-N_i=A_{i}x-v_i$ (16)

$Ax=\left(N_{GPS/leveling}-N_{grav}\right)+v$ (17)

where A is the design matrix, x is the vector of unknown parameters and v is the vector of the unknown random errors coming from GPS and levelling observations and geoid itself. Even though the choice of the appropriate parametric model, AX depends on the distribution, density and the quality of the data used, in this thesis, we used a simplified 7-parameter version of the usual 7-parmeter similarity datum shift transform model to derive the corrector surface model after [3]. It has also been used in previous studies by [18]. The 7-parameter model can be expressed mathematically as:

$\begin{array}{c}{A}_{i}x=\cos {\phi }_i\cos {\lambda }_i{x}_1+\cos {\phi }_i\sin {\lambda }_i{x}_2+\sin {\phi }_i{x}_3+\frac{\cos {\phi }_i\sin {\phi }_i\cos {\lambda }_i}{W_i}{x}_4\\ +\frac{\cos {\phi }_i\sin {\phi }_i\sin {\lambda }_i}{W_i}{x}_5+\frac{{\sin }^2{\phi }_i}{W_i}{x}_6+x_7\end{array}$ (18)

where ${\phi }_i,{\lambda }_i$ , the geodetic coordinates of the GPS/leveling points, e is the eccentricity of the reference ellipsoid and $W_i=\sqrt{1-e^2{\sin }^2{\phi }_1}$ . The vector of the unknown parameters are $x_1,x_2,x_3,\cdots$ and $x_7$ .

A cross-validation of the GNSS/leveling data subset was performed. This is to guarantee that, following the interpolation and transformation process, the customized or hybrid geoid may produce orthometric heights that are dependable to a specific degree of precision. In order to determine the associated orthometric heights, the corresponding geoid values were interpolated from the hybrid geoid file and then deducted from the ellipsoidal heights of these points. As indicated in Table 6, the hybrid geoid (post-fit) is accurate to 3 cm root mean square (r.m.s.) when the derived or transform orthometric height is compared with the observed orthometric height.

Table 6. Results of the cross-validation.

Station’s Attribute	Lat. (dd)	Long (dd)	h(m) Observed	H (m) Observed	H* (m) Derived	(H-H*) (m)	Source
C21	6.9687	9.2468	344.6080	324.6000	324.5950	0.005	OSGOF
A24	10.6038	11.3397	624.6253	605.6290	605.6293	−0.003	OSGOF
D29	11.3880	5.2173	526.7924	505.5490	505.5474	0.002	OSGOF
N25	8.9988	8.0874	591.9590	570.1200	570.1160	0.004	OSGOF
CBL10	5.5393	5.7379	25.8125	4.7930	4.8005	−0.008	OSGOF
L03	7.4173	3.5209	290.3270	264.2840	264.2980	−0.014	OSGOF
PT 7 Emma	4.8239	7.006	33.3790	14.3716	14.4470	−0.075	SG’s Rivers
Uniport	4.8937	6.9144	29.7120	10.8670	10.8654	0.002	SG’s Rivers
PP9	4.8883	7.1445	33.5700	14.4602	14.5040	−0.044	SG’s Rivers
						RMS	3 cm

4. Conclusions

Distortions or residuals are always present in any normal geoid derivation, fitting, and transformation process from 3D GNSS-derived ellipsoidal height to gravity-related 1D orthometric height. The transformation process through interpolation, systematic mistakes caused by the various contributors to the geoid, and inconsistent datums could all be to blame for this. For example, it was found that the gravity and gravity anomaly datasets acquired from BGI were calculated using the GRS67 reference ellipsoid, but the global ellipsoid, such as WGS84 or GRS80, is now used for gravimetric geoid computation. Long wavelength inaccuracy would have resulted from ignoring this. Over Nigeria, the effect is roughly 0.87 mGal. Thus, conversion was carried out from GRS67 to GRS80 prior to computation.

Additionally, it is anticipated that any gravity reduction procedure for gravimetric geoid modeling will result in residual gravity anomalies that are smooth, have negligible indirect effects, and have geophysical significance. It should be noted that the defects present in the Nigerian vertical network have not been eliminated by the fitting and transformation procedure outlined in this work (Table 6). However, it has only harmonized the corresponding points on the two surfaces (N^GNSS and N^GRAV), producing results that are compatible throughout the research region and repeatable and consistent. Using Spherical Harmonic Coefficients (EGM08) as a reference field, the standard remove-restore method was used to generate the geoid. Spherical Fast Fourier Transform (FFT) and Residual Terrain Modelling (RTM) by prism integration were applied. The FFT’s use is based on the idea that it is quick and effective, especially when performing calculations across wide areas where the traditional integration of the Stokes/Molodenskii integral in the space domain takes a long time.

After a thorough evaluation as a quasi-geoid, the gravimetric geoid model over Nigeria was transformed into a classical geoid by applying the 𝑁−𝜁 adjustment. The gravimetric geoid’s strength and details are reflected in the geoid model at short range, while the leveling networks’ trends are reflected at long range. It is crucial to remember that the GNSS/leveling data throughout Nigeria, not the geoid model, is what makes it difficult to utilize GNSS to determine suitable gravity-related (Orthometric) height. It should be noted that systematic mistakes and land uplifts are likely to cause tilts in the final geoid model. Therefore, the geoid surface from fitting is a corrective surface rather than an equipotential surface. A graphic user interface (GUI) using the interpolation program Height Transformation Model, which enables height users to interpolate geoid values and then convert ellipsoidal heights to orthometric heights or vice versa, provides access to the geoid surface and the underlying gravimetric geoid file.

Acknowledgements

The authors wish to thank the Bureau Gravimetrique International (BGI) for the release of terrestrial gravity for this research and acknowledge the producers of the GRAVSOFT software that was used in the computation.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sansò, F. and Sideris, M.G. (2013) Geoid Determination: Theory and Methods. Springer.[CrossRef]
[2]	Adams, R., IIiffe, J., Ziebart, M., Turner, J. and Oliveira, J.F. (2006) Joining up Land and Sea: UKHO/UCL Vertical Offshore Reference Frame. Hydro International, 10, 7-9.
[3]	Hofmann-Wellenhof, B. and Moritz, H. (2006) Physical Geodesy. Springer.[CrossRef]
[4]	Sjoberg, L.E. and Bagherbandi, M. (2017) Gravity Inversion and Integration: Theory and Applications in Geodesy and Geophysics. Springer.[CrossRef]
[5]	Bureau Gravimetrique International (2015) Land and Marine Gravity Data.
[6]	Ince, E.S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., et al. (2019) ICGEM—15 Years of Successful Collection and Distribution of Global Gravitational Models, Associated Services, and Future Plans. Earth System Science Data, 11, 647-674.[CrossRef]
[7]	Alaska Satellite Facility (2020) ALOS PALSAR—Digital Elevation Model (12.5 m Spatial Resolution. https://asf.alaska.edu/data-sets/derived-data-sets/alos-palsar-rtc/alos-palsar-radiometric-terrain-correction
[8]	Wichiencharoen, C. (1982) The Indirect Effects on the Computation of Geoid Undulations. Report No. 336, Dept. of Geodetic Science and Surveying, The Ohio State University, Columbus.
[9]	Barthmeles, F. (2013) Definition of Functionals of the Geopotential and Their Calculation from Spherical Harmonic Models: Theory and Formulas Used by the Calculation Service of the International Centre for Global Earth Models (ICGEM). Deutsches GeoForschungZentrum GFZ: Scientific Technical Report STR09/02. [Google Scholar] [CrossRef]
[10]	Stelios, P.M. (2010) Gravity, Geoid and Earth Observation. Springer.
[11]	Forsberg, R. (2010) Terrain Effects in Geoid Computations. Lectures Note, International School for the Determination and Use of the Geoid, Saint Petersburg, 28 June-2 July 2010, 159-181.
[12]	Tscherning, C.C. (2012) Geoid Determination by 3D Least-Squares Collocation. In: Sansò, F. and Sideris, M., Eds., Geoid Determination, Springer, 311-336.[CrossRef]
[13]	Forsberg, R. (1984) A Study of Terrain Reductions, Density Anomalies and Geo-physical Inversion Methods in Gravity Field Modelling. Department of Geodetic Science and Surveying, The Ohio State University, No. 355.
[14]	Wong, L. and Gore, R. (1969) Accuracy of Geoid Heights from Modified Stokes Kernels. Geophysical Journal International, 18, 81-91.[CrossRef]
[15]	Yahaya, S.I. and Azzab, D.E. (2018) High-Resolution Residual Terrain Model and Terrain Corrections for Gravity Field Modelling and Geoid Computation in Niger Republic. Geodesy and cartography, 44, 89-99.[CrossRef]
[16]	Ezeigbo, C.U. (1991) Estimation Models of Geoid in Nigeria Using Doppler Satellite Observations. In: Rapp, R.H. and Sansò, F., Eds., International Association of Geodesy Symposia, Springer New York, 230-240.[CrossRef]
[17]	Ezeanaka, E.E. (1987) A Study of the Local Geoid Using Integral Equivalents of Stokes and Vening Meinesz Equivalents. Master’s Thesis, University of Nigeria.
[18]	Fotopoulos, G. (2003) An Analysis on the Optimal Combination of Geoid, Ortho-metric and Ellipsoidal Height Data. PhD Thesis, University of Calgary.