Application of Radial Basis Functions for Height Datum Unification

Abstract

Local gravity field modelling demands high-quality gravity data as well as an appropriate mathematical model. Particularly in coastal areas, there may be different types of gravity observations available, for instance, terrestrial, aerial, marine gravity, and satellite altimetry data. Thus, it is important to develop a proper tool to merge the different data types for local gravity field modelling and determination of the geoid. In this study, radial basis functions, as a commonly useful tool for gravity data integration, are employed to model the gravity potential field of the southern part of Iran using terrestrial gravity anomalies, gravity anomalies derived from re-tracked satellite altimetry, marine gravity anomalies, and gravity anomalies synthesized from an Earth gravity model. Reference GNSS/levelling (geometric) geoidal heights are used to evaluate the accuracy of the estimated local gravity field model. The gravimetric geoidal heights are in acceptable agreement with the geometric ones in terms of the standard deviation and the mean value which are 4.1 and 12 cm, respectively. Besides, the reference benchmark of the national first-order levelling network of Iran is located in the study area. The derived gravity model was used to compute the gravity potential difference at this point and then transformed into a height difference which results in the value of the shift of this benchmark with respect to the geoid. The estimated shift shows a good agreement with previously published studies.

1. Introduction

The availability of multi-resolution gravity data corroborates the idea of gravity data fusion to improve their spectral and spatial resolutions for gravity field modelling. Since the integrated dataset holds the inherent information of each contributing dataset, it is possible to enhance the accuracy of the resultant gravimetric model. Overland areas, terrestrial gravity contains the highest frequency of the gravity potential field. As most of the gravity field modelling techniques require (somewhat) regular data coverage, the data gaps are mostly filled with data synthesized from an Earth gravitational model (EGM). In coastal regions, besides terrestrial gravity on land, marine gravity data represent another valuable type of gravity data. However, their particular properties include a high noise level. In addition, marine gravity data may not have sufficient spatial coverage in some regions. In such cases, satellite radar altimetry observations are another valuable data which can be used to derive the medium-frequency gravity anomalies. Nevertheless, the main drawback of satellite radar altimetry data is their lower quality in areas close to coastlines. Indeed, the range between the satellite and water surface, which is one of the main observable quantities in satellite radar altimetry, is particularly erroneous in coastal regions. In such regions, the atmospheric propagation, rough surface topography, significant wave height and wind speed are the main sources of errors in gravity anomalies extracted from satellite radar altimetry [1]. Alternatively, by extracting gravity anomalies from re-tracked satellite radar altimetry, it is possible to reduce the errors in altimetry data. The “Re-tracking” technique is applied to refine the altimetry range; it refers to the analysis of radar pulses received in the satellite antenna [2].

The traditional method of gravity field modelling is to expand the external gravity potential into a spherical harmonic series. However, spherical harmonics require global data support, meaning that for currently available global gravity data, they cannot provide sufficient spatial localization; small regional data variations change the whole set of coefficients in the spherical harmonic expansion. Hence, basis functions with compact or quasi-compact supports are more suitable choices for local gravity field modelling considering various types of gravity data, their spatial distributions, spatially-limited areas of interest and required accuracy. Among available methods, such as the least-squares collocation (LSC) [3], least-squares modified Stokes (LSMS) method [4] and the fast Fourier transform approach [5], the radial basis functions (RBF) have been proven to be more suitable for the combined processing of different types of gravity data, e.g., [6,7,8,9,10,11,12,13,14,15,16]. In all these studies the RBF technique is used for gravity field modelling where in this study we investigate the application of this method for height datum unification. Height datum unification, using other methods, has been investigated in many studies, e.g., [17,18,19,20,21,22,23,24,25,26,27]. One way of height datum unification between regions is through estimating the offset at the levelling benchmark with respect to the global height datum, i.e., the global geoid.

In this article, a local gravity database is developed for regional gravity field modelling using the RBF technique in the southern coastal part of Iran which is used for the determination of the regional geoid model. The importance of determining the geoidal heights in this area is the location of the main benchmark of the first-order levelling network of Iran, called DNG1001. This point is connected to a tide gauge at the Persian Gulf and used to derive physical heights at the levelling network benchmarks. The Iranian levelling network is based on the geoid while the RBF technique is used to determine the quasigeoid. Therefore, all the gravimetric height anomalies must be converted to the geoidal heights. Recently [28] computed an RBF-based geoid model of Iran (IRG2016), where the method used in their study to parametrize the RBF parameters is different than the proposed methodology in this study. Besides terrestrial gravity data, the marine gravity anomalies are also used in the present study, which were not employed in IRG2016. Moreover, the root mean square error (RMS) of the differences between predicted geoidal height and observed ones at the GNSS/Leveling control points, even after a polynomial fitting, is 23 cm for the IRG2016 which is not reliable enough to be used for the purpose of height unification investigated in the present study.

This article is organized into 4 sections. In Section 2, the regional gravity field modelling using RBFs is reviewed. Section 3 focuses on the description of various gravity data used in this study and presents the results of the numerical study. Section 4 then summarizes the main findings of this article.

2. Theory

The spherical harmonic expansion of the disturbing gravity potential $T$ on or above the Earth’s surface is expressed as follows [29] (Sections 2–5):

$$T(r,\phi ,\lambda )={\displaystyle {\displaystyle \sum }_{n=2}^{\infty }}{\displaystyle {\displaystyle \sum }_{m=-n}^n}(\frac{R}{r}{)}^{n+1} {\overline{C}}_{nm} {\overline{Y}}_{nm}\left(\phi ,\lambda \right)$$

(1)

where ($r,\phi ,\lambda$) denotes the position of the computation point in the geocentric spherical coordinate system, ${\overline{C}}_{nm}$ are the fully-normalized Stokes coefficients and ${\overline{Y}}_{nm}$ are the fully-normalized surface spherical harmonics. Assuming the disturbing gravity potential is available on the reference sphere ${\Omega }_R$, the Fourier transform can be used to determine the Stokes coefficients ${\overline{C}}_{nm}$, i.e.,

$$\begin{array}{c}{\overline{C}}_{nm}={\displaystyle \underset{{\Omega }_R}{\int }}T\left(R,\phi {}^{\prime },\lambda {}^{\prime }\right) {\overline{Y}}_{nm}\left(\phi {}^{\prime },\lambda {}^{\prime }\right) \mathrm{d}\sigma \left(\phi {}^{\prime },\lambda {}^{\prime }\right)\end{array}$$

(2)

where $\mathrm{d}\sigma =\cos {\phi }^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }$ is the spherical surface element and $\left(R,\phi {}^{\prime },\lambda {}^{\prime }\right)\in {\Omega }_R$ is the position of the integration point on the geocentric sphere. Substituting Equation (2) into Equation (1), we get:

$$\begin{array}{c}T\left(r,\phi ,\lambda \right)={\displaystyle {\displaystyle \sum }_{n=2}^{\infty }}{\displaystyle {\displaystyle \sum }_{m=-n}^n}(\frac{R}{r}{)}^{n+1} {\displaystyle \underset{{\Omega }_R}{\int }}T\left(R,{\phi }^{\prime },\lambda {}^{\prime }\right) {\overline{Y}}_{nm}\left({\phi }^{\prime },\lambda {}^{\prime }\right) \mathrm{d}\sigma \left(\phi {}^{\prime },\lambda {}^{\prime }\right) {\overline{Y}}_{nm}\left(\phi ,\lambda \right)\\ ={\displaystyle \underset{{\Omega }_R}{\int }}T\left(R,{\phi }^{\prime },\lambda {}^{\prime }\right){\displaystyle {\displaystyle \sum }_{n=2}^{\infty }}(\frac{R}{r}{)}^{n+1} {\displaystyle {\displaystyle \sum }_{m=-n}^n}{\overline{Y}}_{nm}\left({\phi }^{\prime },\lambda {}^{\prime }\right) {\overline{Y}}_{nm}\left(\phi ,\lambda \right) \mathrm{d}\sigma \left(\phi {}^{\prime },\lambda {}^{\prime }\right)\end{array}$$

(3)

By considering the theorem of the spherical harmonics which relates spherical harmonics to the Legendre polynomials [29] (Sections 2–5), i.e.,

$${\displaystyle {\displaystyle \sum }_{m=-n}^n}{\overline{Y}}_{nm}({\phi }^{\prime },\lambda {}^{\prime }) {\overline{Y}}_{nm}(\phi ,\lambda )=\frac{2n+1}{4\pi R^2} {\overline{P}}_n(\cos \psi )$$

(4)

where $\psi$ is the spherical distance between the two points, Equation (3) may be written as follows:

$$\begin{array}{c}T(r,\phi ,\lambda ) ={\displaystyle \underset{{\Omega }_R}{\int }}T\left(R,{\phi }^{\prime },\lambda {}^{\prime }\right){\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}\frac{2n+1}{4\pi R^2} {(\frac{R}{r})}^{n+1} {\overline{P}}_n(\cos \psi ) \mathrm{d}\sigma (\phi {}^{\prime },\lambda {}^{\prime })\end{array}$$

(5)

Introducing the following kernel function:

$$\Phi \left(r,R,\psi \right)={\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}\frac{2n+1}{4\pi R^2} (\frac{R}{r}{)}^{n+1} {\overline{P}}_n(\cos \psi )$$

(6)

We may write:

$$\begin{array}{c}T\left(r,\phi ,\lambda \right)={\displaystyle \underset{{\Omega }_R}{\int }}T\left(R,{\phi }^{\prime },\lambda {}^{\prime }\right) \Phi \left(r,R,\psi \right) \mathrm{d}\sigma \left(\phi {}^{\prime },\lambda {}^{\prime }\right)={\left(\Phi \ast T\right)}_{{\Omega }_R}\left(r,\phi ,\lambda \right)\end{array}$$

(7)

Equation (7) is based on the convolution theorem in spherical analysis on the sphere with the kernel function $\Phi \left(r,R,\psi \right)$. Assuming that the functions $T$ and $\Phi$ belong to the same Hilbert space, their convolution can be represented by a linear combination of the kernels defined in Equation (6), see [30] (Chapter 3). Hence, we have:

$$T\left(r,\phi ,\lambda \right)={\displaystyle {\displaystyle \sum }_{i=1}^M}{\alpha }_i {\Phi }_i\left(r,r_i,{\psi }_i\right)$$

(8)

where ${\alpha }_i$ are the scaling coefficients and $M$ is the number of the kernel functions used in the approximation. The kernel function $\Phi$, which depends only on the spherical distance $\psi$ between the two points, is called the radial basis function (RBF). The triad $\left(r,\phi ,\lambda \right)$ is the 3-D position of the computation point located outside the reference sphere of radius R, and the triad $(r_i=R,{\phi }_i,{\lambda }_i)$ denotes the 3-D position of RBFs located on the reference sphere. The location of RBFs is called the RBF centre. In general, a reproducing RBF kernel function is defined by:

$$\Phi \left(r,R,\psi \right)={\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}\frac{2n+1}{4\pi R^2} (\frac{R}{r}{)}^{n+1} k_n {\overline{P}}_n\left(\cos \psi \right)$$

(9)

In Equation (9), $k_n$ are called the Legendre coefficients which specify the kernel type and influence its behaviour in the frequency domain.

In our example, however, gravity disturbances are not observable and gravity anomalies are used instead as input data which are modelled by RBFs by employing the fundamental equation of physical geodesy. Assuming the location of RBFs is known (this issue will be discussed later in detail), a system of linear equations may be constructed to find the scale coefficients of Equation (8):

$$\begin{gathered} \left[\begin{array}{c}{L}_1\\ L_2\\ ⋮\\ L_{N-1}\\ L_N\end{array}\right]=\left[\begin{array}{c}\begin{array}{ccc}{\varphi }_{1,1}& \cdots & {\varphi }_{1,M}\end{array}\\ \begin{array}{c}\begin{array}{ccc}{\varphi }_{2,1}& \cdots & {\varphi }_{2,M}\end{array}\\ \begin{array}{ccc}⋮& \ddots & ⋮\end{array}\\ \begin{array}{ccc}{\varphi }_{N-1,1}& \cdots & {\varphi }_{N-1,M}\end{array}\end{array}\\ \begin{array}{ccc}{\varphi }_{N,1}& \cdots & {\varphi }_{N,M}\end{array}\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ \begin{array}{c}{\alpha }_2\\ ⋮\\ {\alpha }_{N-1}\end{array}\\ {\alpha }_N\end{array}\right], \\ L=A\alpha , \end{gathered}$$

(10)

where $N$ is the number of observations, $M$ is the number of radial basis functions, ${\Phi }_{N,M}$ are the kernel values used to construct the coefficient matrix, $L$ is the vector of gravity anomalies and $\alpha$ is the vector of unknown scale coefficients to be solved using this system of linear equations.

Different types of kernel functions as RBFs may be used for regional gravity field modelling, see, e.g., [7,9]. It was shown by [31] that almost the same accuracy of gravity field modelling can be achieved for different types of RBFs when the depth of the basis functions is chosen properly. Here the Poisson wavelet kernel is used. This type of the RBF kernel was introduced by [32] and used first in numerical experiments by [33]. This kernel has been proven as a suitable type of the kernel function for gravity field modelling [33]. The Poisson wavelet kernel reads:

$$\begin{gathered} \mathsf{\Phi }\left(\mathit{r},{\mathit{r}}^{\prime }\right)=\frac{1}{4\pi R^2}\left(2{\chi }_{m+1}+{\chi }_m\right), \\ {\chi }_m={\left(\left|\mathit{r}\right|\frac{\partial }{\partial \left|\mathit{r}\right|}\right)}^m\frac{1}{\left|\mathit{r}-{\mathit{r}}^{\prime }\right|} \end{gathered}$$

(11)

where $m$ is the order of the kernel function (here the third order is used). The location of base functions and their bandwidths (Equation (9)) must be defined accurately to properly model the gravity field. Once these critical parameters are defined, the scaling coefficients (Equation (8)) can be estimated using an optimization technique (e.g., Equation (10)). A regular grid, containing the RBF centres, may be defined based on the chosen number of RBFs. The grid is usually slightly larger than the target area to avoid edge effects in the solution. The number of RBFs is chosen according to the RMS of discrepancies between predicted gravity anomalies (using the RBF model) and observed ones at control points. The bandwidth of the Poisson kernel is controlled by the depth of RBFs below the Bjerhammar sphere, see, [7], and needs to be determined, too. The general cross-validation approach is usually used to find the bandwidth of RBFs, although it was numerically shown by [27] that minimizing the RMS value of the discrepancies at the control points could be used to find the optimal bandwidth of RBFs with the sufficient accuracy.

The solution of the system of linear equations of Equation (10) may be found by using the least-squares (LS). The system of linear equations, which is singular in most practical cases, may be caused by using a larger number of RBFs than needed, data gaps, and wrong choice of data bandwidths. One part of the singularity can physically be removed by choosing an optimal number of RBFs as well as by the proper choice of their depths using the RMS criteria described above. The remaining singularity can be removed using a regularization (smoothing) parameter ($\lambda$) in solving Equation (10), see [34]:

$$\left(A+\lambda I\right)\alpha =L.$$

(12)

3. Numerical Study

The selected study area for investigating the performance of RBFs for the unification of height datums spreads over the southern part of Iran, around the Qeshm Island, in the rectangle limited as follows: ${53}^{°}<\lambda <{58}^{°}$ and ${25}^{°}<\phi <{29}^{°}$, see Figure 1.

The terrestrial gravity data provided by the Bureau Gravimetrique International (BGI) and by the National Cartographic Center (NCC) of Iran is the dominant source for providing the high-frequency component of the gravity field. The target area consists of 5427 points with inhomogeneous spatial distribution, see the red and green points in Figure 2. Due to difficulties and high expenses related to collecting terrestrial gravity, typically they are not well distributed over the test area. Terrestrial gravity anomalies, especially when different datasets are merged together, often include many blunders and duplicate points; thus, removing these points is crucial for their use in high-frequency gravity field modelling. In our test case, 56 duplicate points and 102 outliers were eliminated from the database.

The land areas with no terrestrial gravity information are filled by gravity anomalies synthesized from EGM. In our case, EGM2008 [35] up to the full degree/order (2160) is used to fill the gap areas. Generally, Iran is a challenging area for using EGMs as there were not many gravity data provided to global gravity field modellers. [36] compared the terrestrial gravity data of Iran with values synthesized from different EGMs and the EGM2008 model was chosen as the best performing model in Iran with an RMS of discrepancies of 12.9 mGal. A very close value of RMS was also achieved for the EIGEN-6C4 model [37]. The free-air gravity anomalies in the gap areas, mostly in the eastern part of the study area (6752 points) are synthesized on $1^{\prime }\times 1^{\prime }$ using the isGrafLab software (in the tide-free system) and masked out the sea area [38].

The altimetry-based gravity data provide the medium-frequency information of the gravity field. Indeed, emitting a series of radar chirp signals, a satellite altimeter senses ranges reflected from the sea surface or inland lakes which consequently allow for the recovery of the sea surface slopes required for gravity data extraction. Despite the fact that difficult conditions in coastal areas affect the altimetry observations, the quality of satellite radar altimetry data can be enhanced by applying a post-processing algorithm known as waveform re-tracking [39]. In marine parts of the study area, the gravity anomalies were extracted from the global $1^{\prime }\times 1^{\prime }$ marine gravity anomaly data provided by [40]. The blue points in Figure 2 show the spatial distribution of the 20,951 marine gravity points. All available gravity data and their sources are summarized in

Table 1. The summary of the input gravity data.

Free Air Gravity Data	Source	Number of Points
Terrestrial	BGI-NCC	5427
Marine (Re-tracking)	Re-Tracked satellite altimetry	20,951
Global model	EGM2008 up to d/o 2160	6752

.

Figure 3 shows the variation of the merged gravity anomalies over the study area. Please note that an interpolation method is used to plot the image; however, the gravity data mentioned above are used at the scattered points for RBF parametrization. The statistics of the free-air gravity anomalies and the reference field used are summarized in

Table 2. The statistics of the available free air gravity anomalies and the reference field in the study area.

	Min (mGal)	Max (mGal)	Mean (mGal)	STD (mGal)
Free-air gravity anomalies	−211.24	348.57	−13.14	65.50
Reference field	−145.12	174.77	−12.17	57.40

.

The long wavelength part of the gravity field was removed from the merged gravity anomalies by synthesizing the free-air gravity anomalies at the scattered points using EGM2008 up to the degree and order 360, see Figure 4.

To find the optimal number of RBFs, as well as their depth below the Bjerhammar sphere, 5% of the input gravity data (1600 points), were excluded from the gravity database. These points were chosen randomly. Different numbers of control points were chosen, and trial and error showed that excluding 5% of the input gravity data as control points would give a reliable output which was also confirmed by [12]. The control points at this stage are chosen to find the optimal number of RBFs and will be added in the final analysis. As mentioned above, the RBFs in this study are located on the regular grid. In the first step, the depth of RBFs is set to a fixed value of 10 km below the Bjerhammar sphere and a varying number of RBFs is tested to construct grids with the different spatial resolution and predict gravity at the control points. The optimal number of RBFs is estimated based on the minimum RMS of the discrepancies between the predicted free-air gravity anomalies at the control points and observed values. As it can be seen from Figure 5a, adding more than 6500 RBFs does not improve the quality of the gravity field solution, in contrary, more RBFs can worsen the solution as larger data noise is modelled. To determine the depth of RBFs below the Bjerhammar sphere, 6400 RBFs located on the regular coordinate grid and the same RMS criteria were employed to find the optimal depth of the RBF centres below the Bjerhammar sphere. Figure 5b shows the variation of the RMS value with respect to different depths. The optimal depth according to this figure is estimated as $13 \mathrm{km}$.

Once the optimal number of RBFs and their depth below the Bjerhammar sphere are known, the system of linear Equation (12) can be solved. The smoothing parameter in this equation was estimated using the L-Curve method and the regularization tools of Matlab [41]. This value is set to $1.3487\times {10}^{-5}$, see Figure 6.

By determining the scaling coefficients, the 3-D configuration and parameters of RBFs can be used to represent the regional gravity field of the study area; see Equation (8). The differences between predicted and observed gravity anomalies are shown in Figure 7 and their statistics are summarized in

Table 3. The statistics of the discrepancies between the predicted and observed gravity anomalies.

	Min [mGal]	Max [mGal]	Mean [mGal]	STD [mGal]
$\Delta g^{obs}-\Delta g^{RBF}$	−23.90	21.13	0.00	1.41

.

The disturbing gravity potential of the test area can be converted to height anomalies using the Bruns formula [29] (Equations (2)–(144)):

$$\zeta \left(r,\phi ,\lambda \right)=\frac{T\left(r,\phi ,\lambda \right)}{\gamma }$$

(13)

where $\zeta$ is the height anomaly derived by RBF parameterization of the disturbing gravity potential $T$. The height datum in Iran is the geoid, i.e., the height anomalies must be converted to the geoidal heights. This can be done using either the classical formula and/or with more explicit formulas of the geoid-to-quasigeoid separation (cf., [42,43,44]:

$$N\left(R,\phi ,\lambda \right)-\zeta \left(r,\phi ,\lambda \right)\cong \frac{\Delta g^B}{\overline{\gamma }} H$$

(14)

where $H$ is the orthometric height, $\Delta g^B$ is the simple Bouguer anomaly, $N$ is the geoidal height. The gravimetric geoidal heights computed by the RBF parameters are shown in Figure 8. There are 63 GNSS/levelling points available in the study area, see Figure 8. To assess independently the accuracy of the gravimetric geoid model derived in this study, the predicted geoidal heights at the GNSS/levelling points were compared with the geometric geoid.

Table 4. The statistics of the geoidal height test at the GNSS/levelling points.

	Min (m)	Max (m)	Mean (m)	STD (m)
$N^{GNSS/Lev}-N^{RBF}$	−0.11	−0.18	−0.120	0.04

summarizes the statistics of these differences.

As mentioned in the introduction, the purpose of this study is to evaluate the offset of the main benchmark of the Iranian levelling network with respect to the geoid. Particularly, some studies investigated the unification of the Iranian height datum, e.g., Ref. [23] computed the gravity potential value at the above-mentioned levelling benchmark and reported the difference of 0.9 ${\mathrm{m}}^2/{\mathrm{s}}^2$ with respect to the value of $W_0$ which translates to the $9 \mathrm{cm}$ shift between the DNG1001 benchmark and the global geoid. [27] computed the same value using the RBF method and estimated a shift of $21 \mathrm{cm}$. This value is way off the value reported in [23]. One reason for this difference may be using the erroneous satellite altimetry data in the processing done by [27] for determining the geoid. In this article, however, beside using the same terrestrial gravity data as in [27], gravity anomalies derived from re-tracked satellite altimetry are used instead of only satellite altimetry derived anomalies, a new set of terrestrial gravity data is added, and the RBF processing scheme is advanced by using scattered gravity observations (rather than using gridded gravity anomalies). Using these parameters the offset of the local levelling network of Iran is estimated and

Table 5. The estimated offset of the geoid at the DNG1001 point for the levelling network of Iran.

Reference Point	Location	${\mathit{N}}^{\mathit{r}\mathit{e}\mathit{f}}\left(\mathit{m}\right)$	Offset (m) [23]	Offset (m) [27]	Offset (m) This Study
DNG1001	$\lambda =56.1831$ $\phi =27.1658$	−28.961	0.090	0.210	0.115

compares the estimated offset of the DNG1001 benchmark with respect to the geoid computed by three different studies.

4. Summary and Conclusions

We investigated the applicability and reliability of merging multi-resolution gravity data in precise gravimetric geoid modelling in coastal regions. To accomplish this purpose, we generated a new gravity database consisting of terrestrial gravity anomalies and gravity anomalies derived from re-tracked satellite radar altimetry data. Furthermore, the data gaps were filled with EGM2008-based (in the tide-free system) gravity anomalies were computed by the GrafLab software. After producing the integrated gravity database in the form of an equiangular grid using cubic spline interpolation, the remove-compute-restore technique was applied to prepare the required data for regional gravity field recovery. By choosing the Poisson wavelet of the order 3 as the base function, the parameterization of the gravity field model was done with RBFs. In order to find the optimum values of the RBF parameters, minimizing the RMS value of the predicted gravity at excluded control points were utilized.

The shift of the Iran levelling network with respect to the geoid model is computed in this study and compared with previously published studies. The result of this study (11 cm shift) confirms the finding of [23] (9 cm shift) which are both quite far from the result of [27] (21 cm shift). The reason of this difference may be in using the satellite altimetry data-only in the coastal area whereas, in the present study, the satellite altimetry data are refined using re-tracking to reduce their errors in the coastal region of the test area. Besides, a new dataset of terrestrial gravity data is added and the RBF parametrization of the gravity field is advanced.