A combined regional Geopotential Model using optimized global Gravity Field Solutions

To develop a gravimetric geoid, a Global Geopotential Model (GGM) is required to minimise the truncation error arising from using the Stokes integral with a limited number of gravity data points. The choice of a best-fitting GGM determines the accuracy of a gravimetric geoid solution. Selecting a suitable GGM is a rigorous process, requiring both internal and external evaluation of all GGMs available at the International Center for Globa Earth Models (ICGEM). Moreover, GGMs perform differently depending on the wavelength, and it is difficult to obtain a GGM that performs best across the full harmonic spectrum. In this study, a combined GGM is developed from a selection of the most recent and high-resolution GGMs covering Peninsular Malaysia. The selected models are first synthesized harmonically to obtain geoid undulations at collocated GNSS-levelled points, and free air anomalies at randomly sampled points across the study area. These quantities are compared with the observed geoid undulations and point gravity anomalies interpolated from a grid of free air anomalies. The best performing GGMs are then used to produce a combined GGM, by selecting the spherical harmonic coefficients with the best characteristics for every degree. The signal and error spectra of the new GGM are compared with the selected geopotential models. The combined GGM produced a higher cumulative signal to noise ratio (SNR) of 4402.669 compared to all the selected GGMs, with XGM2016 and Eigen-6C following suit with SNR of 4139.561 and 4092.462, respectively. Besides, the new combined GGM performed better across the whole harmonic spectrum than all selected GGMs. The use of combined GGMs in geoid modelling, instead of a single GGM may be more desirable because they can improve the quality of results.


Introduction
The Stokes integral formula is widely used to create a gravimetric geoid model, by computing geoidal undulations using gravity measurements at places on the earth's surface This formula requires gravity data covering the whole earth, yet, discrete gravity data is available only within a spherical cap in practice. A Combination of a global geopotential model (GGM, plural GGMs) with gravity data can IOP Publishing doi: 10.1088/1755-1315/1051/1/012001 2 reduce the truncation error caused by this weakness. For optimal estimation of a gravimetric geoid, it is crucial to choose a GGM capable of recovering gravity field functionals (e.g., gravity anomalies, geoid undulations, and vertical deflections) that are comparable to those generated from GNSS-leveling and terrestrial gravity data.
To date, high-quality global gravity field models have been created from advanced satellite data or derived values with a spatial resolution of 9km, corresponding to a maximum degree of 2190 [1]. The International Centre for Global Earth Models (ICGEM) is a service that provides the scientific community with a high-tech repository of static and temporal global gravity field models of the earth. About 180 GGMs are now available on the ICGEM's website (http://icgem.gfzpotsdam.de/ICGEM/-ICGEM.html).
Several studies have been conducted to determine the applicability of GGMs in Malaysia in general, and Peninsular Malaysia, in particular. In [2], a comparison was made between OSU91A and EGM96 GGMs. It was concluded that the EGM96 model could recover the short frequency signals of both the geoid undulations and the gravity anomalies better. In [3], the EIGEN-6C4 model was used to compute a gravimetric geoid model over Peninsular Malaysia because it provided the best fit amongst the combined models. In the same study, the GO_CONS_GCF _2_SPW_ R4 was considered best amongst the satellite only GGMs.
With the number of geopotential models growing every year, it will be difficult for the user to select the most optimum model for their regional modelling without testing each model. Moreover, different geopotential models perform differently in different wavelengths [4] and hence, the best overall geopotential model selected for an area may be weaker in some wavelengths, and vice-versa. Instead of using a single geopotential model, it would be desirable to get the best out of all available geopotential models.
This study aims to develop a combined geopotential model by extracting the best spectral information from a selection of the most recent and high-resolution GGMs from the International Center for Global Earth models (ICGEM) for future geoid modelling in Peninsular Malaysia. The models are first filtered using GNSS geoid undulations and observed gravity anomalies, to obtain the best performing model(s) from which to extract the spectral information. The new model's internal characteristics are tested together with other geopotential models using spectral analysis. Results show that combined geopotential models perform better than the individual geopotential models, and may improve the quality of regional geoid modelling.

Data used
The most recent and highest resolution GGMs, GNSS-levelled points, and point free-air anomalies were employed in this investigation. The following subsections provide an overview of the data set.

Global Geopotential models
A total of 30 GGMs were chosen for the study, all of which either have a maximum degree not less than 360 or were produced within the last ten years. The GGMs span a variety of models with various input information (e.g., terrestrial gravity, satellite tracking and altimetry data), as well as degree and order diversity. The present study includes models that had been examined and recommended in earlier Malaysian studies. In Table 1, the features of the GGMs, which were retrieved from the ICGEM website, are listed.

Gridded free air anomalies
A 30m by 30m resolution grid of surface free air anomalies was used. This grid was computed from terrestrial and air-borne data supplied by the Department of Survey and Mapping Malaysia (DSMM), and processed as explained in [5]. Table 2 shows the statistical indicators of the gravity data that was used in the investigation.

GNSS-Levelling data
The selected geopotential models were validated using 173 GNSS-levelled sites acquired from the Department of Surveying and Mapping Malaysia (DSMM). Table 3 contains the height statistical indicators, whereas Figure 1 depicts the position of the data points.  Figure 1. Map of free air anomalies used in the study. Red triangles indicate the positions of GNSS-levelled benchmarks.

Synthesis of Gravity functionals
A Global geopotential model is a collection of dimensionless, fully normalized spherical harmonic coefficients ̅ and ̅ as well as their errors ̅ and ̅ , that may be used to simulate the earth's gravity field. These coefficients were calculated using a combination of terrestrial and satellite gravity measurements or a study of satellite observations [6]. Geoid undulations, height anomalies, gravity anomalies, gravity disturbances, and vertical deflections, as well as other gravity field functionals implied by the related GGM, can all be calculated using the coefficients. The disturbing potential, is commonly represented by the expansion [7][8][9][10]: - where is the earth's gravitational constant, is the semi-major axis of the normal reference ellipsoid, is the distance of a point from the origin, ̅ and ̅ are the fully normalized spherical harmonic coefficients of the disturbing potential, are the fully normalized associated Legendre functions for degree and order , and ( , ) are spherical polar coordinates of point . The coefficients ̅ ( = 0) are referred to an ellipsoid of a given flattening. Geoid undulations may be obtained by combining equation (1) with the Brun's equation to obtain [8][9][10][11][12]: - where is the product of the universal gravitational constant and the earth's mass, is a scaling parameter related to a specific GGM, is the geopotential model's maximum degree, and r, , are the geocentric coordinates of the geoid-reduced computation point.

Zero-Degree Terms
With equations (2) and (3), the terms 0 and ∆ 0 are components of the zero-degree terms for the GGM geoid undulations and gravity anomalies, respectively, in relation to the reference normal ellipsoid. They allow the geoid undulation and gravity anomalies synthesized from GGMs to be linked to a specific equipotential surface with 0 and values, by accounting for the changes in masses and potential between the geopotential model utilized and the reference ellipsoid. They may be computed from the formulae [9,13,14]: - where the normal ellipsoid is represented by the parameters and 0 [15] to create the Somigliana-Pizzeti normal gravity field.

Selection of Permanent Tide
Global Geopotential models are provided either in the zero-tide or tide-free, or in both versions [16]. The mean-tide potential cannot be employed in gravity field modelling because it comprises the permanent tide-generating potential, which is generated by masses exterior to the Earth. For an explanation of the pros and cons of using the different permanent tides, the reader is referred to e.g. [16][17][18]. Concerning spherical harmonic coefficients, only the 20 coefficient is affected by the permanent tide. Depending on the desired permanent tide to be used, the 20 coefficients of the selected GGMs may be transformed into the new system using the relation [19,20]: - where − 20 and − 20 are the spherical harmonic coefficients in the tide-free and zero-tide systems, respectively, and 2 = 0.3 is the adopted second-degree love number.
In Peninsular Malaysia, surveying and mapping activities are referred to the Peninsular Geodetic Vertical Datum (PMGVD), which is based on 10 years of tidal observation data [21,22]. The datum was transferred from 0 at the Port Kelang tide gauge to 3.624m at a memorial monument using a combination of precise levelling and gravity survey. Due to lack of luni-solar correction to precise leveling in Peninsular Malaysia, orthometric heights are automatically reported in the mean tide system [18]. In order to have all quantities in the same permanent tide, these heights may be converted into the tide free system and vice-versa using the formulae [22,23]: where − and − are the tide-free and mean-tide levelling heights, respectively and is the geodetic latitude of the GNSS-levelled benchmark.

Degree variances and error degree variances.
Using spherical harmonic coefficients, the signal degree variances of the anomalous potential may be calculated from [24,25]: - and the error degree variances by [24,25]:- where quantities are as defined previously.
Equations (8) and (9) are applicable to all functionals of the gravity field provided the appropriate eigenvalue, is inserted in the equations. The factor corresponding to the various functionals of the disturbing potential are given in Table 4. The signal degree variance signifies the amount of signal power implied by all the coefficients within a specific degree and is commonly referred to as the power spectrum. Conversely, the error degree variance is an expression of how much signal power error of a given anomalous quantity exists for all the coefficients of a specific degree. The variation of power spectra with the degree may therefore be used to describe the rate of decay of the anomalous signal as the degree increases.

3.4.2.
Root mean square error. The root mean square (RMS) by degree of gravity functionals may be computed from [4,24,25]: - while the overall RMS may be obtained from [24,25]: -

(11)
In this study, the root mean squares of the geoidal undulations and gravity anomalies were estimated by the wavelength of the selected geopotential models. Four wavelength types were selected as defined in [4]: -1.
Long wavelengths: gravity field information ranging from degrees n=2 to n=10, equivalent to a linear half-wavelength of 2000km and greater. 2.
Intermediate wavelength: for gravity field information ranging from degrees n=11 to n=100 equivalent to a linear half-wavelength of 200 to 2000km. 3.
Short wavelength: gravity field information ranging from degrees n=101 to n=1000 equivalent to a linear half-wavelength of 20 to 200km. 4.
Very short wavelength: gravity field information ranging from degrees 1001 to ∞ equivalent to a linear half-wavelength less than 20km. The total RMS of the point gravity field functional may be expressed in terms of its wavelength components from [4]: -= ⌈ 2,10 2 + 11,100 2 + 101,1000 2 + >1000 2 ⌉ 1 2

Signal to Noise ratio.
The signal-to-noise (SNR) ratio may be computed both cumulatively or by degree. The SNR by degree may be obtained from [26]: - where variables are as previously defined.

Comparison with GNSS-levelled heights
3.5.1. Computation of bias. Observed Geoid heights, may be computed from spirit-levelled orthometric heights, and GNSS-measured ellipsoidal heights, ℎ using the famous equation [27]: - These heights are independent of the gravimetrically derived geoid heights, making them appropriate for evaluating geoid models, including GGMs. Equation (2) was used to determine the model-derived geoid undulations implied by the GGMs. The differences or bias for each GGM was computed for statistical analysis using the following equation [27]: - where is the geoid height obtained from the GGM. The basic statistical indicators, such as the mean, minimum, maximum and standard deviation, were then obtained.

Detection of outliers.
To assess the quality of the GNSS-levelled points, the observed geoid undulations may be compared with undulations implied by a high resolution GGM as in equation (15). Standardized residuals are then computed using: - where is the point residual and is the standard deviation.
Observations with a standard residual >3 were considered as outliers in this study.

Fitting of parametric models.
To minimize the effect of systematic errors and datum inconsistencies [28,29] between the GNSS-levelling derived geoid undulations and those derived from GGMs, various parametric models were fitted into the observations using the equation [27]:- where denotes a vector of the unknown parameters, and is a design matrix corresponding to the unknown coefficients of a pre-selected parametric model. In this study, three, four, five and seven parametric models were used, respectively as follow [30]: -Three parameters: Four parameters: Five parameters: Seven parameters: where and are the flattening and first eccentricity, respectively of the reference ellipsoid and is given by √1 − 2 2 . The unweighted least-squares adjustment was used since there was no information on the accuracy of both the observed and the synthesized geoid undulations. The solution to the least-squares adjustment was of the form [31,32]: - where is the vector of residuals, the standard deviation, = − , the degrees of freedom, the number of evaluation points and the number of parameters in terms of the parametric model.

Comparison with gravity anomalies.
To carry out the assessment using gravity anomalies, 1000 points were randomly selected within the study area, and free-air gravity anomalies interpolated at their positions, using the MATLAB in-built function interp2. At the same points, gravity anomalies were obtained from synthesis of the selected GGMs using equation (3), with the resulting residual anomalies computed without terrain effects using the equation [33]: -

Combined GGM
The combined GGM was obtained by comparing the signal to noise ratio of the geoid undulation and gravity anomaly spectra by degree, for all the selected GGMs. The coefficients ̅ and ̅ , (and their errors), whose GGM provided the largest signal to noise ratio for a particular degree were selected for the combined GGM. The mass and the radius of the combined model was calculated from a weighted mean of the values given for the selected GGMs.

Results and Discussion
External validation was performed on all 30 GGMs utilizing 173 GNSS-leveled locations and observed free-air anomalies. On the GNSS-leveled points, an outlier detection mechanism was used by comparing the observed geoid undulations with those obtained from EGM2008, one of the highest resolution GGMs. Using the approaches of section 3.5.2, one outlier was identified and removed.

Computation of observed gravity field functionals
From the known orthometric and ellipsoidal heights of the GNSS-levelled stations, observed geoid heights were estimated using equation (14). 1000 plots were then randomly selected within the study area and free air anomalies were interpolated from the free air anomaly grid. The statistics of the observed GNSS undulations and interpolated free air anomalies are shown in Table 5.

Synthesis of gravity field functionals
Equations (2) and (3) were used to synthesize geoid undulations and free air anomalies from the geopotential models at the GNSS-levelled points and at the 1000 random points, respectively, using MATLAB functions developed by the authors. The earth's geocentric gravitational constant ( ), and other constants, were derived from the GGMs.

Computation of Zero-Degree Terms
The zero-degree terms for the geoid undulation and free air anomalies were first computed for each GGM using equations (4) and (5), and the results were added to the synthesized gravity field functionals. A value of =3.986004415e14 was used in all GGMs. The reference ellipsoid provided the reference gravitational constant, 0 , the mean Earth radius, , the normal potential, 0 , and the mean normal gravity, γ.

Comparison using GNSS-levelled points
The differences between the geometric and the GGM-based geoid undulations were computed using equation (15), and Table 6 shows the statistical indicators of the results. In the table, the GGMs are ranked in ascending order of the standard deviations obtained. As shown by the mean values in Table 6, there is evidence of some bias between the geopotential of the Malaysian vertical datum's zero-height surface and the equipotential surface specified by the IERS conventional value Wo = 62636856.00 m2/s2, which was used in the development of most of the selected GGMs. Long and medium wavelength inaccuracies in the spherical harmonic coefficients are most likely to blame for these discrepancies [36,37]. Least-squares parametric fitting was used to model the systematic errors using the models mentioned in section 3.5. The residuals were then combined with the observations (h-H), and the results were compared to the GNSS-level undulations. As indicated in Table 9, this improved the standard deviations. Among the selected GGMs, the XGM_2019_2159 and XGM_2019e produced the lowest standard deviation of 0.071m each of the difference between GNSS based and geopotential-based undulations. Generally, most of the higher resolution GGMs (>360) performed well with standard deviations < 0.1m. Among the Satellite only GGMs, GO_CONS_GCF_2_TIM_R6, GOCO06s and GO_CONS_GCF_2_DIR_R6 performed best with standard deviations of 0.134. 0.135 and 0.137m respectively. After removal of the systematic errors, the standard deviation improved in all models, with the higher resolution GGMs still performing better than lower resolution GGMs. When the parametric models were compared, the seven-parametric model produced lower standard deviations for all GGMs, hence it was used to rank the GGMs' performance in terms of undulations.
. Overall, the XGM_2019e, XGM_2019e_2159 and SGG-UGGM1 produced the best fit with standard deviations of 0.0615, 0.0615 and 0.0617 m, respectively, with all the high resolution GGMs obtaining < 9cm standard deviations. For the satellite only GGMs, GO_CONS_GCF_2_DIR_R6,

Comparison using free air anomalies
The free-air gravity anomalies synthesized using the geopotential models were compared to the observed free-air anomalies at the 1000 random points. The residual anomalies were produced after subtracting the synthetic anomalies from the observed anomalies, and the statistical results are reported in Table 8. As can be observed from Table 8, the higher resolution GGMs outperformed the lower ones. EGM2008, SGG-UGM-1 and EIGEN-6C produced the best fit in terms of gravity anomalies in Peninsular Malaysia because of their low standard deviations of 8

Computation of combined GGM
From the previous evaluation, ten of the best ranked GGMs in each category were selected for the purpose of comparing their spectral behavior, in terms of the signal to noise ratio of the disturbing potential. The selected GGMs are depicted in Table 9. 4.6.1. Selection of permanent tide. The tide-free system was used in computation of the new geopotential model, since this is the system used for the ITRF coordinate system to evaluate geopotential models. All geopotential models were therefore converted into the tide-free system by changing the 20 coefficient as in equation (6). The results of the conversion for the selected models are depicted in Table  10.

Optimization of harmonic coefficients.
Using the spherical harmonic coefficients, the signal to noise ratio of the disturbing potential was computed per degree for all GGMs selected. The disturbing potential was used since all other gravity field functionals are its derivatives, and therefore related to them. For every degree, up to a maximum of 2500, the coefficients and their error estimates producing the larger signal to noise ratio were selected for the combined GGM. The mass of the earth and reference radius were adopted from a weighted mean of the values used for the selected GGMs. Values of 3.9860044150e14 and 6378136.3874 were obtained for the earth's gravitational constant and its radius, respectively. Table 11 shows the GGMs whose coefficients were selected for the respective degree and order. From the results of Table 11, the combined model was derived from 6 GGMs, with Eigen-6C containing most of the very short wavelength coefficients, according to the classification of section 3.4. In Table 12, a sample of the normalized spherical harmonic coefficients of the model, which was called PM_com_1c in this study, is shown. Figure 2 and Figure 3 depict the geoid undulations and free air anomalies, respectively, implied by the new model.

Validation of New Geopotential model
The new geopotential model was validated by comparing its spectral information with the other GGMs used in this study. Using the procedures of section 3.4, the signal and error spectra of different gravity functionals were compared up to a maximum degree of 2500. Table 13 shows the results of the SNR of the signal and error spectra by wavelength, respectively, in terms of both gravity anomaly and geoid undulation.  From the analysis of Table 13, the new combined geopotential model produced the largest signal to noise ratio in the disturbing potential, gravity anomaly and geoid undulation which, expectedly, had the same values across the whole spectrum under investigation. The same can be said of Figure 4 and Figure  5 where the signal to noise ratios are plotted by harmonic degree for the gravity anomaly and geoid undulation, respectively.

4.7.2.
Comparison of GGMs using error spectra. The cumulative errors of the gravity anomaly and geoid undulation signal are depicted in Figure 6 and Figure 7, respectively, by harmonic degree according to equation (12). From the above figures, the new GGM exhibits a lower cumulative error spectra in terms of both the gravity anomalies and geoid undulations.

Conclusion
This study developed a combined geopotential model from a selection of the best performing, most recent and highest resolution GGMs available at the ICGEM. The models were selected based on their fit to a set of collocated GNSS-spirit levelled points and a grid of gravity anomalies covering the Peninsular Malaysia.
From the geometrical comparison of selected GGMs, it was revealed that the XGM_2019_2159 and XGM_2019e produced the lowest standard deviation of 0.071m each of the difference between GNSS based and geopotential-based undulations. Most of the higher resolution GGMs (>360), however, performed well with standard deviations < 0.1m. Among the Satellite only GGMs, GO_CONS_GCF_2_TIM_R6, GOCO06s and GO_CONS_GCF_2_DIR_R6 performed best with standard errors of 0.134. 0.135 and 0.137m respectively. After removal of the systematic errors, the standard deviation improved in all models, with the higher resolution GGMs still performing better than lower resolution GGMs. Overall, the XGM_2019e, XGM_2019e_2159 and SGG_UGM_1 produced the best fit with standard errors of 0.0614, 0.0615 and 0.0619m, respectively. For the satellite only GGMs, GO_CONS_GCF_2_DIR_R6, GO_CONS_GCF_2_TIM_R6 and GO_CONS_GCF_2_SPW_R4 were the best performers with standard errors of 0.1014, 0.1019 and 0.1022 m respectively. In terms of gravity anomalies, EGM2008, SGG-UGM-1 and EIGEN-6C produced the best fit because of their low standard errors of 8.149, 8.187 and 8.349 mGal, respectively.
The new GGM, which was named PM_com_1c, was computed by selecting from the best performing GGMs, the harmonic coefficients producing the largest signal to noise ratio of the disturbing potential. After comparison of PM_com_1c with the other GGMs, it was revealed that it produced larger cumulative signal to noise ratio in the disturbing potential, gravity anomaly and geoid undulation, than all the individual GGMs, in all parts of the harmonic spectrum. This study recommends the use of PM_com_1c or other regionally combined GGM for gravimetric geoid modelling in Peninsular Malaysia. It is expected that, when such regional geopotential models are used to provide the long wavelength components of the gravity field, the quality of regional geoid models will improve.