Application of the nonlinear optimisation in regional gravity field modelling using spherical radial base functions

The gravity field is a signature of the mass distribution and interior structure of the Earth, in addition to all its geodetic applications especially geoid determination and vertical datum unification. Determination of a regional gravity field model is an important subject and needs to be investigated and developed. Here, the spherical radial basis functions (SBFs) are applied in two scenarios for this purpose: interpolating the gravity anomalies and solving the fundamental equation of physical geodesy for geoid or disturbing potential determination, which has the possibility of being verified by the Global Navigation Satellite Systems (GNSS)/levelling data. Proper selections of the number of SBFs and optimal location of the applied SBFs are important factors to increase the accuracy of estimation. In this study, the gravity anomaly interpolation based on the SBFs is performed by Gauss-Newton optimisation with truncated singular value decomposition, and a Quasi-Newton method based on line search to solve the minimisation problems with a small number of iterations is developed. In order to solve the fundamental equation of physical geodesy by the SBFs, the truncated Newton optimisation is applied as the Hessian matrix of the objective function is not always positive definite. These two scenarios are applied on the terrestrial free-air gravity anomalies over the topographically rough area of Auvergne. The obtained accuracy for the interpolated gravity anomaly model is 1.7 mGal with the number of point-masses about 30% of the number of observations, and 1.5 mGal in the second scenario where the number of used kernels is also 30%. These accuracies are root mean square errors (RMSE) of the differences between predicted and observed gravity anomalies at check points. Moreover, utilising the optimal constructed model from the second scenario, the RMSE of 9 cm is achieved for the differences between the gravimetric height anomalies derived from the model and the geometric height anomalies from GNSS/levelling points.


INTRODUCTION
The gravity field of the Earth represents the mass distribution and structure of the Earth's interior. Beside its geophysical applications, the gravity field plays a major role in geoid/height determination. Therefore, it is still an important subject to study and investigate for geodetic and geophysical purposes. Different methods have been developed and applied for regional gravity field modelling in different study areas, and they are in one way or another flavours of interpolation methods, except for the deterministic methods. The accuracy of each method is different from another and each has its own advantages and disadvantages.
Generally, there are two forms for interpolating the gravity data at the surface of the Earth. First, the classical interpolation that is based on spherical harmonics, which is used in global gravity field modelling (Pavlis et al., 2012). Second, interpolation of data in a small area on the Earth. In this case, spherical harmonics are not suitable interpolants because in a regional scale, the size of data area might be too small to resolve the low degrees. Hence, for the purpose of the regional gravity data fitting, local base functions like spherical radial basis functions (SBFs) provide better accuracy. The SBFs are known tools for gravity field modelling from satellite, airborne and terrestrial gravimetric data. There are several types of SBFs to regionally model the gravity field, such as the pointmass kernel, the radial multipoles, Poisson wavelets, and the Poisson kernel, which are all functions of the inverse Euclidean distance from the data points to the centres of these functions. The point-mass kernels have been used in this respect frequently for example by Barthelmes (1988), Heikkinen (1981), Lin et al. (2014), Vermeer (1995). Barthelmes and Dietrich (1991) showed that the number of required point-mass kernels can be reduced by optimising the location of their centres. The radial multipoles, which are the higher derivatives of the point-mass functions, were utilised for modelling the gravity field, see, e.g., Foroughi and Tenzer (2014), Marchenko (1998), Marchenko et al. (2001). Poisson wavelets, were introduced by Holschneider et al. (2003) and later used for local gravity field representing by, e.g., Chambodut et al. (2005), Klees et al. (2008), Klees and Wittwer (2007), Panet et al. (2006), Safari et al. (2014), Tenzer et al. (2012). In addition to the gravity field, Chambodut et al. (2005) applied the Poisson wavelets for modelling the magnetic field and proposed it as an efficient SBFs when data do not have global coverage and regular distribution. Spherical wavelets with different scales are also suitable tools to represent the gravity field and applied by, e.g., Schmidt et al. (2004Schmidt et al. ( , 2005. The centres of the SBFs play a significant role in the accuracy of the modelled gravity field and if they are determined optimally, all of the SBFs can be applied equally well with insignificant changes in accuracy . The disturbing potential of the Earth was represented as a linear combination of the SBFs by Klees et al. (2008) and a data-adaptive method was used to select the optimal number and locations of the SBFs. Tenzer et al. (2012) investigated the gravity field modelling using the Poisson wavelets of order 3 on different spherical equal-angular grids. For every choice of the grid, they found the optimal depth of SBFs inside the Bjerhammar sphere. They showed that applying topographical correction to the observed gravity data reduces the number of SBFs and improves the fit to the data. However, it implies a systematic bias in the adjusted gravity disturbance. Lin et al. (2014) presented a two-step method for regional gravity field modelling, applied an interval constrained Quasi-Newton minimisation technique to find the optimal locations of the centres of the point-mass kernels and after that, used the least-squares estimation for fitting the SBFs to the disturbing potential. Also, they investigated the effect of the number of point-masses and the direction of optimisation on the estimation. In spite of a high number of unknowns in this approach for the centres of the SBFs, the Genetic Algorithm (GA) was used by Mahbuby et al. (2017) first and then the conjugate gradient method was utilised to fit the SBFs to the synthetic data.
In the stochastic method of local gravity field modelling an ill-condition system of equations is created, which needs to be solved by regularisation. Orr (1995) investigated the regularisation forward selection of the SBFs' centres. Ophaug and Gerlach (2017) studied the above mentioned three parameterisation methods for regional geoid computation. They utilised the spline kernel as SBF and compared the three methods theoretically and numerically. Using synthetic noise-free data, all applied methods agreed on the millimetre level. Also, in several studies, various regional gravity field modelling methods have been used for the Auvergne test area in central France that some of them are mentioned. Yildiz (2012) interpolated the Gravity field and steady-state Ocean Circulation Explorer (GOCE) vertical gravity gradient data using collocation method. Janák et al. (2014) combined GOCE gravity gradiometry data with terrestrial gravity anomalies to solve the geodetic boundary-value problem numerically and obtained the standard deviation of 10 cm for their model. Abbak and Ustun (2015) developed a software based on Least-Squares Modification of Stokes (LSMS) formula to compute the geoid model. Lin et al. (2019) showed that free-positioned point-mass parameterisation for regional gravity field modeling performs better than fixed-positioned method in rough regions. Goyal et al. (2021) computed deterministic and stochastic quasigeoid models and compared the results with some previous studies.
Generally, determination of the optimal locations of the SBFs is a nonlinear leastsquares problem. The Gauss-Newton optimisation method is the simplest approach to solve a nonlinear least-squares problem. It can be viewed as a modification of Newton method with line search. For this purpose, the second-order term of the Hessian is excluded to obtain the descent direction i.e. the Hessian is approximated via a quadratic function in terms of Jacobean (Wright and Nocedal, 1999). Finding the descent direction in this method may lead to solving an ill-posed linear system of equations where a regularised Gauss-Newton method is needed (Qi-Nian, 2000). The convergence of this optimisation method in the Hilbert space was studied by Bauer et al. (2009). The above mentioned minimisation method is really efficient in terms of the processing time and memory requirements compared to the evolutionary algorithms like GA. Doicu et al. (2002) investigated the application of the iteratively regularised Gauss-Newton method to solve the nonlinear ill-posed inverse problems arising in atmospheric remote sensing.
One important issue in using the SBFs is the optimal position vector of the centre of functions, which is normally done by trial and error. However, in this study, the regularised Gauss-Newton optimisation technique is applied for this purpose, which is a nonlinear process. This approach makes this study different from those have been done so far. Regularised Gauss-Newton method is in the category of line search algorithms and convergence with low number of iterations is the most prominent feature of this approach. Stud. Geophys. Geod., 65 (2021) This method will be discussed in details. In order to show application of our method, two scenarios of interpolating the gravity anomalies with two types of SBFs are investigated and compared. In the first scenario, point-mass kernels are utilised for modelling of the gravity anomaly and in the second scenario, another kernel is defined which is a function of point-mass and its radial derivative. This kernel is applied for regional gravity field modelling. Furthermore, applying this kernel in the second scenario leads to estimation of the disturbing gravitational potential from gravity anomalies. In this study, first a closedloop test with synthetic data obtained from Earth Gravitational Model (EGM2008) (Pavlis et al., 2012) is done and thereafter, free-air terrestrial gravity anomalies in the Auvergne test area are used to construct the optimal regional gravity field model with the proposed method.

SPHERICAL RADIAL BASE FUNCTIONS AND POINT-MASS KERNELS
In this section, our mathematical tool, which is parameterisation (expansion) in terms of SBFs is presented. Also, the point -mass kernel and another kernel involving point-mass and its radial derivative are introduced and compared as two important base functions for gravity anomaly interpolation.

. 1 . S p h e r i c a l r a d i a l b a s e f u n c t i o n s
Consider a function of the coordinate vector x, expressed by the following linear combination of SBFs with the centre at y i : Here c i are the coefficients which should be estimated in an optimal way, m is the total number of the kernel functions, which are defined as follows Tenzer and Klees, 2008):

 
Tn i P  x y represents the Legendre polynomial of degree n, dot symbol    stands for inner product, R stands for the radius of the Bjerhammar sphere , In the case of known c i and the position of kernels, these can be used to determine the function   f x , which can be regarded as a forward problem. However, if c i are not known, in fact, each linear combination presenting   f x will be an observation equation having m unknown c i and known   f x values. Here, a system of equations can be organised and solved for c i over an area covered by the   f x values. The main issue is that the position of the kernels should be known for such a purpose.
Equation (2) is a general definition of a SBF, now consider that  n = 1, in this case, the kernel will change to the well-known point-mass kernel, having the following inverse relationship with distance: Here 2 i  x y means the distance between the point-masses and computation points. The frequency filtering property of this kernel depends on the location of the point-masses. If they are deep inside the reference sphere, 2 2 i y x will be small and by getting power n+1 it reduce the power of high frequency portion of the signal comparing to the case that the ration is larger. Therefore, deep point-masses can well approximate the smooth data or can be applied for determining the long wavelength structure of the data and vice versa. Figure 1 shows the plot of the point-mass function when the point-masses are respectively located at the depth of 25% and 15% of the radius of the reference sphere. As expected, the shallower point-masses can be used for recovering the higher frequencies of the data, as the kernel function will be sharper than the case where the point-masses are deep inside the sphere. 2 . 1 . 1 . G r a v i t y a n o m a l y i n t e r p o l a t i o n Now, consider that the gravity anomaly can be mathematically presented by the following linear combination of point-masses: where   , i  x y is the point-mass function given by Eq. (3).
According to the fundamental equation of Physical Geodesy, the relation between the gravity anomaly and the disturbing potential T is (Heiskanen and Moritz, 1967): If we assume the disturbing potential outside of the Bjerhammer sphere as a linear expansion of the point-masses, then considering Eq. (5), the gravity anomaly can also be expanded in the following way: where (7) Figure 2 shows the plot of the kernel (7), which is sharper comparing to the point-mass kernel with less bandwidth at the same penetration depth.
As we observe, there are two different scenarios of interpolating the gravity anomaly with different SBFs kernels  and , using Eq. (4) or Eq. (6). Expanding the gravity anomalies in terms of point-masses has the advantage of simplicity in Eq. (4) but it does  7)) at depth of 25% of the radius of the unit sphere. not produce any other functional of the gravity field. In the case of using Eq. (6), c i will also be used to provide the disturbing potential T. Indeed, when gravity anomaly data are written as in Eq. (6), the estimated c i can also be used in an expansion of point-mass functions to calculate the disturbing potential.

. 1 . . G a u s s -N e w t o n m e t h o d
After the determination of the initial regular grid of SBFs, the regularised Gauss-Newton method, which is in the category of the Quasi-Newton optimisation methods is applied. It implements a modification to find the Hessian and the gradient of the objective function in terms of Jacobian. All line search algorithms, e.g., Gauss-Newton method, include the following two steps: a) A descent direction k y for the objective function F is calculated in each iteration by the Newton method or a Quasi-Newton method. b) Thereafter, a suitable step size t is determined via an exact line search or a backtracking line search method until: Finding a descent direction alone does not ensure the reduction of the objective function. Thus, an appropriate step length t makes the line search algorithm reasonable. We use the backtracking line search with parameters, , as the following procedure (Boyd and Vandenberghe, 2004): In the classic Newton method to find a descent direction for the objective function   F y the following set of the linear equations should be solved: y stand for the Hessian and gradient of the objective function, respectively (Boyd and Vandenberghe, 2004). In least-squares problem, the objective function is: where   r y is the residual of the fitting at point y. The derivatives of   F y are (Wright and Nocedal, 1999): If the number of data is n, then   y J denotes an n  3m Jacobian matrix of first partial derivatives of the residuals with respect to first, second and third components of the kernel positions (i.e. y i1 , y i2 , y i3 ). The gradient g and Hessian H of the objective function can be expressed in terms of Jacobian by (Wright and Nocedal, 1999): Neglecting the second term of Eq. (14) we obtain: which is positive semidefinite but may be close to singular. Hence, a regularisation method can simply be applied for solving Eq. (10) and determining the updates to the initial vector y.

. . I n t e r p o l a t i o n o f t h e r e s i d u a l g r a v i t y a n o m a l i e s w i t h S B F s
Now, let us explain how this method can be applied for interpolating gravity anomalies. Here, we define residual gravity anomalies which are the difference between the observed anomalies and those generated from the EGM2008 (Pavlis et al., 2012) global gravity model to degree and order 360. This means that the long wavelength portions of the gravity anomalies are reduced by the EGM2008 model. Our goal is to perform a remove-compute-restore (RCR) technique for improvement of our interpolation procedure (Abbak et al., 2012;Featherstone et al., 2004;Sjöberg, 2005). The long wavelengths are removed, and the interpolation is done on residual anomalies and the long wavelength will be restored to the interpolated anomalies.
where OBS g  and EGM g  are the observed gravity anomalies and those computed from EGM2008, respectively.
In the first scenario, the interpolation of the residual gravity anomalies by point-mass functions, which is a nonlinear least-squares problem, turns into an inequality constrained minimisation problem as follows: The inequality constraint means that the initial values of y i are properly selected within the Bjerhammar sphere with the radius R. This means that even after updating the initial values they should remain inside the Bjerhammar sphere.
This problem is solved in two phases. Phase 1 is a pre-analysis process which we look for a proper initial value for the optimisation. Phase 2 involves two steps, i.e. optimisation and solution of the system of linear equations, to find the expansion coefficients in Eq. (4).
In the second scenario, for regional gravity field modelling the disturbing potential is expressed as a linear combination of the point-mass kernels and consequently the gravity anomaly is shown by expansion of Eq. (6). Thus, the interpolation is done by solving the following optimisation problem: where d 1 and d 2 are the upper and lower limits of i y , this means that the norm of updated initial positions should be between these limits. Also, This problem is an all-direction optimisation problem, which is sensitive to the chosen upper and lower depth limits (d 1 and d 2 ) (Lin et al., 2014). Considering the logarithmic barrier function and the central path parameter in the objective function, the minimisation problem will change to: where t  stands for the parameter of the central path. The central path tends towards the optimal solution of the problem as the parameter t  grows (Boyd and Vandenberghe, 2004). This problem is also solved in two steps, that is assigning appropriate initial values and improving these values during the optimisation process, and finding the corresponding coefficients.

. . 1 . S o l v i n g m i n i m i s a t i o n p r o b l e m b y t h e G a u s s -N e w t o n m e t h o d
In the previous section, we wrote two minimisation problems (17) and (20). Two important issues in such an optimisation problem are the position of the point-masses and the kernels . Here, first, a truncated Gauss-Newton method is applied to determine the optimal positions of kernels and second, regularisation method is utilised to find the unknown coefficients of the expansion c i . Let us write the following general form for the objective function (20): where The gradient and Hessian of   i F y are: In the theory of Newton and Quasi-Newton optimisation methods, convergence requires the step direction to be a descent direction which will happen once the Hessian is positive definite. In practice, the pure Newton method may fail to converge to a minimum in regions where the objective function is not convex. The second derivative matrix in this problem is not necessarily positive semidefinite, and it is possible that the direction, which results from the solution of Eq. (10) does not indicate a reduction in the objective function. So the truncated Gauss-Newton method for this respect is applied. In this approach, we solve Eq. (10) via the conjugate gradient (CG) method iteratively and terminate the implementation if the negative curvature is encountered. The following algorithm shows how to solve an optimisation problem with truncated Gauss-Newton method (Wright and Nocedal, 1999): o Compute a search direction p k (inexact Gauss-Newton step) by applying the CG method to o The Newton step p k is defined as the final CG iteration;  Set In this section, we present a general scheme for applying our method. This scheme has two phases: Step 1: 1. Considering the initial location of the centres of SBFs on a regular spherical grid inside the Bjerhammar sphere. 2. Estimating the coefficients via truncated singular value decomposition (TSVD) method of regularisation. 3. Comparing different initial grids and finding the proper initial location of the kernels.
Step 2: 1. Given the central path parameters  > 1, t  and m, which is the total number of the inequality constraints, i.e. twice the number of kernels, backtracking parameters , , t, initial values y 0 form phase 1 and a stop criterion , i.e., a tolerance for function value decrease in two consecutive iterations. 2. Calculation of the regularised Gauss-Newton or the truncated Gauss-Newton descent directions (in the first and second scenarios, respectively); i.e. The interpolation matrices in Eqs (4) and (6) are defined as: , for i = 1, …, n, and j = 1, …, m . (27) Utilising the decomposed form of  A or  A and the least-squares solution and TSVD method for regularisation, the scaling coefficients c i will be estimated (Hansen, 1990(Hansen, , 1994. Hence, using the Gauss-Newton minimisation algorithm and a TSVD regularisation method, the optimised locations of the centres of the point-masses or the kernels , and the coefficients of Eqs (4) and (6), will be calculated.
Thereafter, the height anomaly will be estimated at GNSS/levelling points in the study area. The geometric height anomaly  geometric is defined as the difference between the geodetic height h, and the normal height H normal at GNSS/levelling points: Adding long wave effects of the disturbing potential from EGM2008, to the residual disturbing potential computed via the superposition of the point-mass kernels in the second scenario, the total disturbing potential at each evaluation point x, is obtained. The gravimetric height anomaly is defined as the ratio of the total disturbing potential to the gravitational acceleration of the reference ellipsoid  (Moritz, 1980): Usually there is a systematic difference between the gravimetric and geometric height anomalies due to the systematic deformations of the national height network and tide system and also long-wavelength errors in the gravimetric height anomalies .

NUMERICAL RESULTS
In this section, a closed-loop test based on synthetic data is presented to ensure the proper performance of the method. To do so, both stated scenarios in Section 2.2 are applied with regular spherical grid data provided by EGM2008 with the resolution of 0.02  0.02. The Auvergne area which is limited between the longitudes 2.53.5E and the latitudes 4546N is used for testing our approach. The short wavelengths of the freeair gravity anomalies in the study area are obtained by subtracting the computed values up to degree and order 2190 from the long wavelengths up to degree and order 360. Moreover, 100 points are excluded from the provided synthetic dataset as test data to evaluate the performance of the method.
In order to implement the first scenario, namely the interpolation of the residual freeair gravity anomalies in the study area, a regular spherical distribution of point-mass kernels at an initial depth of 15 km below the Bjerhammar sphere is considered. The number of utilised SBFs is 10% of the synthetic data. The criterion for considering the initial depth (d) is to reduce an objective function f o that considers the simultaneous root mean square error (RMSE) of reconstruction at test data and norm of the obtained coefficients   c d from solving the inverse problem at each depth.
In Fig. 3, we can see the values of the objective function at different depths, and the initial depth of 15 km with a regular spherical grid is the best introduced initial value.
Choosing the proper initial depth not only makes the accuracy of the initial reconstructed model suitable even without optimising the kernel location, but also leads to find the optimal location of the kernels with fewer iterations. As can be deduced from Fig. 4, the objective function in Eq. (17) is minimised in 7 iterations.
Finally, after the optimisation process (17), the optimal depth and proper horizontal positions of point-masses are achieved; see Fig. 5. The noteworthy point in this figure is that since the synthetic data distribution is regular, the optimal horizontal positions of the point-masses are close to the introduced initial ones, but their depths vary significantly around the initial depth. In fact, a major contribution to improve the accuracy of the optimal model compared to the model made with the initial location of the point-masses is related to the optimal depths. Also, the closed-loop test based on the introduced synthetic data is done to investigate the efficiency of the stated method for the second scenario. Here, the kernels  are used to parameterise the gravity anomalies and point-masses are applied to model the height anomalies at the study area with the same coefficients. The purpose of the second scenario is to compare the performance of the point-mass kernels and the kernels  in regional gravity field modelling. Moreover, the second scenario, i.e., parameterisation of the residual gravity anomalies based on the kernels , has advantages over interpolation with point-masses. The advantages of the second scenario is that, first, the utilised base  30)); the value of 15 km with a regular spherical grid is considered as the best initial depth.
Stud. Geophys. Geod., 65 (2021) functions has a better ability to model the details, as is evident by the comparison of Fig. 1a and Fig. 2. Because at the same depth, the kernel  is sharper than point-mass and of the less bandwidth. Second, the disturbing potential and gravimetric height anomaly can be modelled with the same coefficients. Thus, the expansion of the residual gravity anomalies in terms of the kernels  with regular spherical distribution at various constant depths is investigated. The number of utilised base functions is the same as the first scenario, i.e., 10% of the synthetic dataset. The objective function (30) is plotted against different depths and the proper initial depth is archived. Figure 6 illustrates the values of the objective function for different depths. As is evident, 20 km is the best initial depth and a depth of 10 to 30 km has led to appropriate values of the objective function (30) and therefore the values (d 1 and d 2 ) in Eq. (18) are considered 10 and 30 km, respectively.
The constrained optimisation problem (18) is solved with the proposed method in Section 2.2.1 with only 13 iterations and decreasing objective function (20) is illustrated in Fig. 7. The stop criterion is  = 0.01, for both problems (17) and (18).
The optimal depths of the kernels  and their proper horizontal positions are shown in Fig. 8. For synthetic dataset, the radius of the Bjerhammar sphere (R) of 6367.4 km is considered. As can be deduced from the figure, the radius of each optimised kernel  varies from R  d 2 to R  d 1 . Also, the initial regular horizontal distribution has changed to some extent.
In the end, the final optimised model is compared with the primary model with initial SBFs' location in terms of accuracy. As can be seen from Table 1, the model made with the initial location of the kernels, i.e., the regular spherical distribution and the appropriate initial depth obtained from the analysis of Figs 3 and 6, has good accuracy. Moreover, with a small number of iterations, the position vectors of the kernels can be optimised to obtain the optimal model with very high accuracy. The expression "misfit" is used for the root mean square of the differences between the reconstructed values by the model at the observation points and the observed values. The geometric height anomalies and the gravimetric height anomalies, introduced in Section 2.3.1, are calculated at 16 scattered GNSS/levelling points and shown in Fig. 9a. Also, the differences between the gravimetric height anomalies, calculated by the regional gravity field modelling method used in the second scenario, and the geometric height anomalies are shown in Fig. 9b. The GNSS/levelling points in the study area can be seen as the black stars in Fig. 10a. Root mean square error of the differences is 6 cm for the closed-loop test based on the synthetic data. Also, the mean and the standard deviation of the differences between gravimetric and geometric height anomalies are 5 and 2 cm, respectively.
After examining the performance of the method for the synthetic dataset, the terrestrial free-air gravity anomaly data over the study area is utilised to represent the regional gravity field model. The long wavelengths of the gravity data are removed by the global gravity model up to degree and order 360. In fact the higher frequencies of the gravitational signal (residual gravity anomalies) are interpolated. The total number of these residuals is 3699. In order to validate our gravity models over the study area, their accuracies are tested against 100 check points, which are not included in our interpolation process. The check points are selected to be almost scattered throughout the entire study area; see Fig. 10b. The minimum, maximum and mean value of the geodetic heights (above the reference ellipsoid WGS84) are 272.64 and 1714.92 and 710.52 m, respectively, showing somehow that the area is relatively rough. Figure 10a, shows the map of geodetic heights over the study area and Fig. 10b is the map of the distribution of the gravity data, the locations of the check points and GNSS/levelling data. The minimum, mean, maximum, and standard deviation of the residual gravity anomalies are, respectively, 47.65, 8.43, 66.86, and 23.93 mGal.  To interpolate the observed terrestrial residual gravity anomalies with point-masses, the regularised Gauss-Newton method with TSVD is applied with 1089 point-masses, which are about 30% of the total number of data. The synthetic data has a regular distribution unlike the real one. In addition, the real data contains higher gravitational frequencies, which need to be modelled accurately. Therefore, a larger number of base functions is required for gravity field modelling using the real data. SBFs are highly dependent on spatial distribution of the data. Any change in the geographical location of data, leads to a new set of SBFs. When the number and values of the base functions change, due to a new data distribution, then a new search interval should be considered for the base function depths. This is the reason that the estimated distribution and depths of the SBFs from the synthetic data differ from the real one.The advantage of our approach is that an accurate gravity anomaly model can be obtained with almost any initial values and few point-masses. However, in this method, the gravity field functional such as the height anomaly is not available. Here, the radius of the Bjerhammar sphere is considered 6367.4 km and the initial grid of point-masses with depths from 60 to 200 km with 10-km interval are tested as the initial depths for a regular grid. The regular grids with their depths are considered as initial values for our nonlinear optimisation, which is also trial and error method for finding the depth and resolution of grids when using the SBFs such as point-mass kernels. These initial values are applied to the nonlinear optimisation process to determine the optimal position and depth of the point-masses. This means that they will no longer be on a regular grid and a constant depth. Table 2 summarises the misfit and RMSE between the modelled gravity anomalies and our 100 check points before and after optimisation of the position of point-masses. Initial grids with a depth of 110 to 150 km lead to primary models with an RMSE of less than 5 mGal at check points. The initial grid at a depth of 120 km, or radius of 6250 km led to the most accurate primary model. As expected, optimisation of the depth and horizontal position of the point-masses with this initial values, also leads to the most accurate model with a misfit of 1.6 mGal and an RMSE of 1.7 mGal at check points. The optimal parameters are obtained in only 6 iterations. Optimisation with low number of iterations is the advantage of the applied method. Significant improvements are seen in the misfits and RMSE after optimisation with all introduced initial regular grids at constant depths. Primary with initial point-mass kernel parameters 0.3 0.4 Optimised with optimal point-mass kernel parameters 0.1 0.2 Primary with initial parameters of the kernels  0.2 0.2 Optimised with optimal parameters of the kernels  0.1 0.1 Fig. 9. a) Gravimetric and geometric height anomalies, and b) the differences between the gravimetric and geometric height anomalies at 16 GNSS/levelling points. The obtained optimal model is considered for presenting the residual gravity anomaly, see Figs 11 and 12. As can be seen in Fig. 11a, the optimal radius of the 1089 SBFs varies from 6224.4 to 6267.4 km. Also, the horizontal variations of the kernels' positions due to the optimisation process with respect to the initial regular grid is demonstrated in Fig. 11b.
The estimated residual gravity anomaly model with optimal parameters, ranging from 50 to 70 mGal is illustrated in Fig. 12b. In addition, its differences with respect to the observed gravity anomalies (see Fig. 12a) are shown in Fig. 12c. As can be observed from this figure, most of the significant differences are around the longitude of 3E and in almost all latitudes in the study area. Comparing Fig. 12a with Fig. 10, it can be concluded that the observed gravity anomalies are highly correlated with geodetic heights.
3 . 2 . L o c a l d i s t u r b i n g p o t e n t i a l a n d h e i g h t a n o m a l y m o d e l l i n g Here, we use the truncated Gauss-Newton method with TSVD to find the optimal position vectors of the kernels  and the coefficients of the expansion. The prominent advantage of this method is the ability to calculate the disturbing potential in addition to the gravity anomaly in the study area. To solve this problem, the limits (d 1 and d 2 ) that are used to construct the constraints in Eq. (18) are important, otherwise, the disturbing potential cannot be modelled. To diagnose the proper limits (d 1 and d 2 ), we have to calculate RMSE of parameterisation at check points by each primary grid of SBFs at a specific depth. For various numbers of the kernels , the proposed constrained optimisation is accomplished, and the optimal distribution of the kernels and their corresponding coefficients are calculated. As can be seen in Table 3, the initial regular grid, which is at a depth of 150 km with 1089 kernels, i.e., about 30% of the data, has the least misfit and RMSE at the check points. Therefore, it is selected as the initial value and the SBFs locations are optimised in the constrained optimisation process. The results show a misfit of 1.4 mGal and an RMSE of 1.5 mGal at the check points. Limits (d 1 and d 2 ) used in the optimisation are set according to Table 3. The initial depth is increased in steps of 10 km to obtain the appropriate depth. Since the initial depth from 130 to 170 km has resulted in an RMSE of less than 4.5 mGal at the check points, this interval is considered as the search interval for the optimal depths of the kernels  and subsequently, due to optimisation in the Cartesian coordinate system, the horizontal distribution of SBFs will also be optimised. Figure 13a shows the optimal kernels depths in the search interval from 130 to 170 km, which is equivalent to a radius from 6237.4 to 6197.4 km. Figure 13b shows the optimal horizontal distribution of the kernels  in the study area. Another point from the comparison of Tables 2 and 3 is that if the same number of SBFs is used, the proper depth of the kernels Fig. 11. a) Optimal radius of point-masses, and b) horizontal distribution of point masses after optimization.
 is more than the proper depth obtained for the point-masses. Furthermore, the gravimetric height anomalies for various numbers of the kernels  are calculated and RMSE of the differences between the obtained height anomalies from the estimated disturbing potentials and geometric height anomalies at 16 GNSS/levelling points in the area are listed in Table 4. Fig. 12. a) The observed gravity anomalies, b) the gravity anomaly model based on 1089 optimal point-mass kernels, and c) the difference between the model and observed gravity anomalies. Table 4 shows the appropriate initial depth and limits (d 1 and d 2 ) for optimisation problem (18) for different number of the kernels . Several important points can be deduced from this table. First, increasing the number of kernels to a threshold can improve the accuracy of the primary and optimal models, but utilising more kernels has no significant effect on RMSE of the gravity anomaly at check points or RMSE of the differences between gravimetric and geometric height anomalies. Also, applying more kernels leads to increased computational complexity and requires more memory. The mentioned threshold in this study is considered to be 1089 SBFs or equivalent to 30% of data. According to Table 4, using more kernels has not any significant effect in reducing the RMSE at the check points and the bias between geometric and gravimetric height anomalies at GNSS/levelling points. This bias is 9 cm from the presented optimal model in the study area, which is 3 cm more than the value estimated by the synthetic data. Also, the mean and the standard deviation of the differences between gravimetric and geometric height anomalies are 8 and 2 cm, respectively. Due to the presence of higher frequencies in the real data and irregularities in data distribution, it is expected that the modelling accuracy will be less than the case using the synthetic data. Furthermore, the system of equations, which is solved for estimating the coefficients from the real data is more illconditioned than that from the synthetic data. Furthermore, as shown in Table 4, the appropriate initial depth and the value of the limits (d 1 and d 2 ), decrease with increasing the number of the kernels .
Finally, utilising 1089 kernels , with the proper bound limits R  d 2 , R  d 1 for the radius of the functions and an initial regular grid at the depth of 150 km inside the Bjerhammar sphere, the optimal gravity anomaly model is constructed in only 7 iterations. The difference between the presented model and the observed gravity anomaly is shown  Fig. 14. Here again, it should be mentioned that there is a high correlation between the observed geodetic heights and gravity anomalies. The scatterplot of them is shown in Fig. 15. The correlation coefficient between them is 0.91 and no significant correlation exists between the estimated residuals and the heights. In the presented optimal models, the high values of residuals are mostly seen around the longitude of 3E.

CONCLUSIONS
Two types of constrained optimisation problems were reviewed for local gravity anomaly modelling and compared in terms of accuracy and efficiency. The aim of presenting the first problem was only to interpolate the residual gravity anomalies by simple point-mass functions and not to consider other gravity field's functionals. This problem by choosing the strictly feasible and appropriate initial value is converted to an unconstrained problem which results in an accurate model of the local gravity anomalies Fig. 13. The same as in Fig. 11, but for SBFs of type .
with point-mass kernels. The aim of solving the second problem was to model the gravity anomaly along with the disturbing potential. The residual gravity anomalies are expanded in terms of the kernels . Also, the obtained coefficients can be used for modelling the disturbing potential with point-mass kernels. In fact, this problem is a constrained optimisation problem and sensitive to the chosen bounds of the radius of the kernels  as the constraints.
In order to compare the performance of these two types of SBFs, the same number of both types of functions are used to solve the simulation problem with synthetic data and modelling with terrestrial gravimetric data. The number of each type of functions is considered to be 10% of the dataset for interpolation of synthetic data and 30% of the dataset for interpolation of the terrestrial gravimetric data. The results of both studies reveal that, first, the kernels  performed somewhat better than point-masses in terms of modelling accuracy. Second, using the computed coefficients in the second scenario, the disturbing potential and consequently gravimetric height anomaly can be calculated with acceptable accuracy. Hence, the proposed method with any type of SBFs leads to an optimal model which improves the accuracy of the primary models with regular distribution of kernels at a constant depth, significantly. Table 4. Root mean square error (RMSE) of the gravity anomalies and height anomalies for different methods and different numbers of the kernels  with their corresponding depth limits d 1 and d 2 .
Usually, there is a considerable bias between geometric and gravimetric height anomalies obtained from the interpolation. No evident improvements can be found in the RMSE of the differences between geometric and gravimetric height anomalies via the optimisation process. However, optimisation always implies better accuracy in gravity anomaly modelling. This is because the height anomalies are not sensitive to the short wavelengths. The bias between the geometric and gravimetric height anomalies at 16 Fig. 14. The same as in Fig. 12, but for SBFs of type .
GNSS/levelling points in the study area is obtained 6 cm for synthetic dataset and 9 cm for terrestrial gravimetric dataset with the final optimal model. Increasing the number of the kernels  to a threshold provides a significant improvement on the accuracy of the gravity anomaly models and gravimetric height anomalies constructed by both primary grids and optimal kernels but using more functions than this threshold is not worthwhile in terms of accuracy and memory, in particular, it does not have a significant effect on improving the accuracy. Finally, another important result is that the proper initial grid depth inside the Bjerhammar sphere and the appropriate bound limits utilised to construct the constraints in the constrained optimisation problem are reduced by increasing the number of the kernels .
Open Access: This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/license/by/4.0), which permits use, duplication, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.