Regularization methods for the combination of heterogeneous observations using spherical radial basis functions

Various types of heterogeneous observations can be combined within a parameter estimation process using spherical radial basis functions (SRBF) for regional gravity field refinement. However, this process is in most cases ill-posed, and thus, 10 regularization is indispensable. We discuss two frequently used methods for choosing the regularization parameter which are the L-curve method and variance component estimation (VCE). Based on these two methods, we propose two new approaches for the regularization parameter determination, which combine the L-curve method and VCE. The first approach, denoted as ‘VCE + L-curve method’, starts with the calculation of the relative weights between the observation techniques by means of VCE. Based on these weights the L-curve method is applied to determine the regularization 15 parameter. In the second approach, called ‘L-curve method + VCE’, the L-curve method determines first the regularization parameter and it is set to be fixed during the calculation of the relative weights between the observation techniques from VCE. These methods are investigated based on two different estimation concepts for combining various observation techniques. All the methods are applied and compared in six study cases using four types of observations in Europe. The results show that the ‘VCE + L-curve method’ delivers the best results in all the six cases, no matter using SRBFs with smoothing or non-smoothing 20 features. The ‘L-curve method + VCE’ also gives rather good results, generally outperforming the cases just using the L-curve method or VCE. Therefore, we conclude that the newly proposed methods are decent and stable for regularization parameter determination when different data sets are combined and can be recommended regardless of the type of SRBFs used.


Introduction
Gravity field modeling is a major topic in geodesy, and it supports lots of applications including physical height system realization, orbit determination and solid earth geophysics.To model the gravity field, a proper approach needs to be set up to represent the input data as good as possible.The global gravity field is usually described by spherical harmonics (SH), due to the fact that they are fulfilling the Laplacian differential equation and are orthogonal basis functions on a sphere.However, in the spatial domain, a global coverage of data sets is not always fulfilled sufficiently, and it is well-known that SHs cannot Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.represent data of heterogeneous density and quality in a proper way (Schmidt et al., 2007).Regional gravity refinement is, thus, performed for combining different observation types such as airborne, shipborne or terrestrial measurements which are only available in specific regions.Different regional gravity modeling methods have been developed during the last decades, e.g., the statistical method of Least Squares Collocation (LSC) (see Krarup, 1970;Moritz, 1980;Pail et al., 2010) or the method of mascons (mass concentrations) (see Rowlands, 2005).The method based on spherical radial basis functions (SRBF) will be the focus of this work.
The fundamentals of SRBFs can be found amongst others in Holschneider et al. (2003) as well as Freeden and Michel (2004).
Due to the fact that SRBFs are isotropic and characterized by their localizing feature, they can be used appropriately for regional approaches to consider the heterogeneity of data sources; examples are given by Marchenko et al. (2003), Schmidt et al. (2007), Lieb et al. (2016).There are many factors in SRBF modeling that influence the accuracy of the regional gravity model, e.g., the shape, bandwidth, locations of the SRBFs and the extension of the data zone.Tenzer and Klees (2008) compared the performance of different types of SRBFs using terrestrial data, Bentel et al. (2013a, b) gave a comparison of nine different SRBFs in regional gravity field modeling based on simulated data.Bentel et al. (2013a) also studied the influence of point grids, and the results show that the differences between SRBFs are much more significant than the differences between different point grids.Another detailed investigation about the locations of SRBFs can be found in Eicker (2008), where methods for choosing a proper bandwidth were also introduced.Bentel (2013a) discussed the reasons for edge effects and Lieb (2017) provided a possible way to choose area margins in order to minimize edge effects.
After setting up all the factors, heterogeneous data sets are often combined within a parameter estimation process; two combination models (see Schmidt et al., 2015) are introduced and applied in this work.One model takes the relative weightings between the observation techniques into consideration while the other one relies on an equally weighted scenario.Regional gravity modeling is usually an ill-posed problem due to (1) the number of chosen basis functions, i.e. the SRBFs, (2) given data gaps, and (3) the downward continuation.Thus, regularization is in most cases inevitable in the parameter estimation process.Bouman (1998) discussed and compared different regularization methods, including Tikhonov regularization (Tikhonov and Arsenin, 1977), truncated singular value decomposition (TSVD, Xu, 1998), and iteration methods (Schuh, 1996).Choosing an appropriate regularization parameter is a crucial issue for a proper regularization.Generalized cross validation (GCV, Golub et al., 1979), L-curve criterion (Hansen, 1990;Hansen and Oleary, 1993) and variance component estimation (VCE, Koch, 1999;Koch and Kusche, 2002) are the three most commonly used methods for estimating the regularization parameter.
However, there are not many studies that compare the performance of each method, and the existing publications do not reach an agreement indicating which gives the best, i.e., the most realistic results.Kusche and Klees (2002) compared GCV and the L-curve method, and the results show that the L-curve criterion always yields over-smoothed solutions; the same results were indicated by Xu (1998).Naeimi (2013) showed that the L-curve method provides satisfactory results while VCE and GCV cannot regularize the regional solutions sufficiently.Bentel (2013) presented that the L-curve method leads to fairly good results for noise-free data but does not perform as good as VCE in the case of noisy observations.Naeimi et al. (2015) Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.investigated how the performance of the regularization method changes when different types of SRBFs are used.The L-curve method delivers the best results when a non-smoothing kernel (Shannon) is applied, while the opposite when smoothing kernels are used.Besides, Lonkhuyzen et al. (2001) showed that the knowledge of the variances of the observations is not a guarantee for obtaining good solutions.
Thus, the purpose of this study is to find out the best-performing method for regularization parameter determination.We will (1) compare the performance of the L-curve method and VCE based on the aforementioned two combination models, (2) propose two new methods which combine the L-curve method and VCE together, (namely 'VCE + L-curve method' and 'Lcurve method + VCE') and compare the results to the ones obtained using the L-curve method or VCE alone, and (3) test the stability of each method when different SRBFs are applied.This work is organized as follows: in Section 2, we present the fundamental concepts of SRBFs, different types of gravitational functionals and SRBFs are also represented briefly.Section 3 discusses the parameter estimation, the Gauss-Markov model as well as the two combination models.Section 4 is dedicated to the regularization method, the L-curve method, VCE and the two newly proposed combination methods.In Section 5, the study area, the simulated data used in the study are presented as well as the results.The performance of all five methods for choosing the regularization parameter is compared.Finally, a summary and conclusions will be given in Section 6.

Regional gravity field modelling using SRBF
In general, a spherical basis function (,   ) related to a point   with position vector   on a sphere Ω  with radius R and an observation point  with position vector x can be expressed by , (1) (Schmidt et al., 2007), with  =  •  =  • [cos cos, cos sin, sin]  , where  is the spherical longitude,  is the spherical latitude,   =  •   ,   is the Legendre polynomial of degree n and   is a Legendre coefficient which specifies the shape of the SRBF.
With the spherical basis function (Eq.1), a harmonic signal () can be described as where K is the number of basis functions.The unknown coefficients   can be evaluated from the observations.As will be shown in the following subsection, using these coefficients, any functional of () can be described.

Gravitational functionals
Various functionals can be derived from the gravitational potential based on field transformations:

Disturbing potential
The disturbing potential T is defined as the difference between the total gravity potential W and the normal gravity potential U where the latter is the potential related to the level ellipsoid.The gravity potential W consists of two parts, the gravitational potential V and the centrifugal potential Z, i.e.

Gravitational potential difference
The satellite gravity field mission Gravity Recovery and Climate Experiment (GRACE, Tapley et al., 2004) 1.

Gravity disturbance
The gravity disturbance is generally used in airborne and terrestrial gravity field determination (Alberts, 2009).The gravity disturbance vector δ is expressed as the gradient of the disturbing potential In spherical approximation, the magnitude of the gravity disturbance can be written as

Gravity gradients
Equipped with a 3-axis gradiometer, the satellite mission Gravity Field and Steady-State Ocean Circulation Explorer (GOCE, Rummel et al., 2002) 1.

Gravitational functionals in terms of SRBFs
In this study, all observations are simulated in the sense of disturbing gravity field quantities, i.e. disturbing potential differences ∆ = (  ) − (  ) for GRACE, the first order radial derivatives   for terrestrial and airborne observations as well as the second order radial derivatives   for GOCE.For each type of observation, the adapted basis functions are listed in Table 1.
Basis functions adapted to other functionals of the disturbing potential which are not used here are listed in Koop (1993), Lieb et al. (2016) and Lieb (2017).

Types of spherical radial basis functions
Since it is not possible to reach perfect localization in both the spectral and spatial domain due to the uncertainty principle (Freeden, 1998;Ozawa, 2003), we want to find an appropriate compromise between these two domains.Different types of SRBFs can be found amongst others in Schmidt (2007), Bentel (2013a, b).Three types of SRBFs are studied in this work, including functions with smoothing features (Blackman and Cubic Polynomial) and without smoothing features (Shannon).
The Shannon function has the simplest representation; its Legendre coefficients are given by The Blackman function is derived from the Blackman window; its Legendre coefficients are given by where In case of the Cubic Polynomial (CuP) function, the Legendre coefficients are given by a cubic polynomial, namely is a certain degree to which the SRBFs are expanded, representing the cut-off degree in the frequency domain.These three functions and their referring Legendre coefficients for   = 255 are plotted in Figure 1, it visualizes the characteristics in the spatial and the spectral domain correspondingly.
In the spatial domain, the Shannon function has the sharpest peak but the strongest oscillations; the Blackman function has less oscillations than the Shannon and the CuP function has the weakest ones.However, in the spectral domain, the Shannon function gets the strongest localization due to its exact band limitation without losing any spectral information; the Blackman function has a smoothing decay at the higher frequencies of the function; the CuP function has an even stronger smoothing decay and thus, extracts less spectral information compared to Shannon and Blackman.Therefore, in this study, we use the Shannon function for estimating the coefficients   within the analysis step to reduce the loss of signal content, and then use the Blackman function for the synthesis step to reduce erroneous systematic effects due to oscillations.The same experiments will be applied using the CuP function as well to test if different SRBFs will affect the performance of the regularization method.

Parameter Estimation
To determine the unknown coefficients   from Eq. ( 2), the method of parameter estimation is used in this study.This process allows different types of observations to be combined considering their individual strength and favorable features (Schmidt et al., 2015).

Gauss-Markov Model
For one single observation, i.e. a functional of the disturbing potential T, the observation equation reads (,   ) represents the adapted SRBFs as listed in Table 1.Collecting the observations ( 1 ), ( 2 ), …, (  ) in the  × 1 observation vector f, the Gauss-Markov model The least-squares adjustment can be applied to the model (Eq.17) as long as the design matrix is of full column rank (Schmidt et al., 2015).Then the solution reads  ̂= (  ) −1    (18) Due to the aforementioned three reasons, the normal equation matrix  =    is ill-posed or even singular.For handling this problem, we introduce an additional linear model as prior information.  is the  × 1 expectation vector of the coefficient vector ,   is the corresponding error vector and (  ) is the  ×  covariance matrix of the prior information with   2 the unknown variance factor and   the given positive definite weight matrix.
Combining the two models (Eq.17) and (Eq.20) yields the extended linear model Now the least-squares adjustment can be applied and leads to the normal equations The variance factors  2 and   2 can either be given as prior information or estimated within a VCE, then the solution reads  ̂= (   +   ) −1 (   +      ) (23) wherein  =  2   2 ⁄ is the regularization parameter, see Koch and Kusche (2002) and Schmidt et al. (2007).

Combination models
To combine different types of heterogeneous data sets for regional gravity field modeling, combination models (CMs) need to be set up (see e.g.Schmidt et al., 2015).
In general, let   with  = 1, … ,  be the observation vector of the p th observation technique, such as   = [  ( 1 ),   ( 2 ), … ,   (   )]  ,   and   are the corresponding error vector and the design matrix.Note that for different techniques, the data are observed as different gravitational functionals and thus, the adapted SRBFs as discussed in the Sect.
2.1 should be applied accordingly, and For the combination of the P observation techniques, including the additional linear model for the prior information (Eq.20), an extended Gauss-Markov model can be formulated (Lieb, 2017) [ where   is the   ×   positive definite weight matrix of the  ℎ observation technique.
CM 1: We assume that for each technique p the variance factor   2 is the same, i.e.  1 2 =  2 2 = ⋯ =   2 =  2 .Hence, with Thus, this combination model is transferred into the single observation model (Eq.21) and the estimation of the coefficient vector  ̂ can be obtained from the Eq. ( 23) if the regularization parameter  is known.
CM 2: Since different data sets have different spatial resolution and spectral characteristics, the assumption made in CM 1 is not always the most accurate case.Therefore, they can be combined in a way which takes the individual variance component of each observation technique into account.
Applying the least-squares method to Eq. ( 25), the extended normal equations read Solving Eq. ( 27) with given values for the variance factors, we obtain +     ) , such that the regularization parameter  =  1 2   2 ⁄ (Koch and Kusche, 2002), and the factors ⁄ express the relative weightings of the observation vector   with respect to  1 .

Choice of the regularization parameter
A critical question of regularization is the selection of an appropriate regularization parameter  (Kusche and Klees, 2002).In the following, the L-curve method and the VCE will be explained in more detail.

L-curve method
The L-curve is a graphical procedure for regularization (Hansen, 1990;Hansen and OLeaary, 1993;Bouman, 1998;Hansen, 2000).Plotting the norm of the regularized solution ‖ ̂ −   ‖ against the norm of the residuals ‖ ̂‖ = ‖ ̂ − ‖ by changing the numerical value for the regularization parameter  shows a typical 'L-curve' behavior, i.e. it looks like the capital letter 'L' (see Fig. 3).The corner point in this 'L-shaped' curve means a compromise of the minimization of the solution norm and the residual norm, and thus can be interpreted as the 'best fit' point that corresponds to the desired regularization parameter.

VCE
Variance component estimation is a useful method when several data sets need to be combined in a parameter estimation procedure (Koch, 1999;Naeimi, 2015).The variance components are estimated by an iterative process, starting from initial values for   2 ,   2 and ending in the convergence point.The estimations read where   and   are the partial redundancies computed following Koch and Kusche (2002), and the residual vectors  ̂ and  ̂ are given by

Combination of VCE and the L-curve method
Two ways of combining VCE and the L-Curve method are discussed and applied in this study, namely 'VCE + L-curve method' and 'L-curve method + VCE'.

'VCE + L-curve method'
Figure 2 illustrates the procedure of the 'VCE + L-curve method'.In the first step, the VCE is applied according to the combination model CM 2. This step gives the regularization parameter  VCE and the relative weighting factors ω  .The weighting factors ω  are then used in the L-curve method to regenerate a new regularization parameter  (Fig. 3).In this case, the coefficient vector  ̂=  ̂ is calculated for a group of changing regularization parameters  using Eq. ( 30).
Thus, the final solution is computed using Eq. ( 30) with the weights ω  and the new regularization parameter  from the Lcurve criterion.

'L-curve method + VCE'
Figure 4 illustrates the procedure of the 'L-curve method + VCE'.In opposite to the 'VCE + L-curve method', in the 'L-curve method + VCE' the L-curve method is applied according to the combination model CM 1 first.A regularization parameter  L−curve is obtained in the first step, and it is used for defining the value of   2 in the variance component estimation.
After each iteration within the VCE, the value of   2 is set to  1 2 / L−curve again, with the new value of  1 2 obtained in this iteration.In this case, the regularization parameter  calculated from the L-curve method will be kept, but the relative weighting factors ω  are recomputed in each iteration step.The final solution is computed using Eq. ( 30) with the relative weights ω  and the regularization parameter  L−curve .
To summarize, the purpose of these two proposed methods is to bring the L-curve method and VCE together, and test if the regularization results will be improved.'VCE + L-curve method' fixes the relative weights of each observation technique first and tries to find a 'best fit' regularization parameter; while 'L-curve + VCE' fixes the regularization parameter first and then tries to find the relative weights for each observation technique.

Data description
The data used in this study are obtained from the ICCT (Inter-Commission Committee on Theory) Joint Study Group 0.3, part of the IAG (International Association of Geodesy) programme running from 2011 to 2015.The observation data are simulated from the Earth gravitational model EGM2008 (Pavlis et al., 2012) and are provided along with simulated observation noise.
The standard deviation of the white noise is set to 8 • 10 −4 m 2 s 2 ⁄ for GRACE, 10 mE for GOCE, 10 μGal for the terrestrial data and 1 mGal for the airborne data.The study area chosen here is 'Europe' where the validation data are also provided on geographic grid points in terms of disturbing potential values T with different grid resolutions (30'×30' and 5'×5') and different spherical harmonic resolutions (maximum degree 250 and 2190).
Figure 5 illustrates the available observation data as well as the validation data.Three types of observations are included: 1. Satellite data: provided along the real satellite orbits of GRACE (green tracks in Fig. 5) and GOCE (red tracks).GRACE data span a one month period, and GOCE data cover a full repeat cycle of 61 days.

Model configuration
A Remove-Compute-Restore approach is applied in this study, i.e., from each type of observations, the background model EGM96 (Lemoine et al., 1998) up to spherical harmonic degree 60 is removed and restored after computation.The background model serves additionally as prior information, and in this case, the expectation vector   can be assumed to be the zero vector (Lieb, 2017).We assume that the coefficients have the same accuracies and are uncorrelated, thus,   = , where  denotes the identity matrix.Further, we set   =  by assuming the measurement errors to be uncorrelated and the observations to have the same accuracy.
In the analysis step we use the Shannon function for estimating the vector  ̂ of the unknown coefficients   related to the grid points   within the area Ω  of computation (see Fig. 6) from the measurements available within the area Ω  of observations.In the following synthesis step the Blackman function is used for calculating the output gravity functional within the area Ω  of investigation.It has to be mentioned that the points   within the area Ω  of computation are defined by a Reuter grid.
Margins  have to be defined between the three areas to minimize edge effects in the computation process (Lieb et al., 2016).
In this study, we conducted the experiments using different margin sizes, and the one which gives the smallest RMS error is finally chosen.
The aforementioned five methods for choosing the regularization parameter (Table 2) are applied to six groups of data sets (Table 3), respectively.
The computed disturbing potential   is compared with the corresponding validation data   and assessed following three criteria: 1. Root mean square error (RMS) of the computed disturbing potential   with respect to the validation data   over the 2. Correlation coefficient between the estimated coefficients   collected in the vector  ̂ and the validation data   .
The reason that this correlation can be used as a criterion is that the estimated coefficients   reflect the energy of the gravity field at their locations.The energy   at location   can be expressed by , (34) (Lieb, 2017).For   → ∞ and   = 1, the relation (Eq.34) equals approximately   =   2 .
The same criterion is used as a quality measure by Bentel et al. (2013a) and Naeimi et al. (2015).
3. Correlation coefficient between the recovered gravity field   and the validation data   .

Results
For the sake of brevity, only the results of two study cases (A and F) are detailed here.However, results obtained from all study cases, including the RMS error and the correlations between the estimated coefficients   and the validation data   of each method are summarized in the Tables 6 and 7, respectively.To test the stability of our new methods, the same experiments are also applied using the CuP function for analysis and synthesis as a comparison scenario.The results are listed in the Tables 8 and 9.

Study case A
GRACE observations and terrestrial observations with a 30' resolution are combined.The maximum degree in the expansion in terms of SRBF is set to   = 300 for the Shannon function (Eq.12) in the analysis step as well as  1 = 250 and   = 350 for the Blackman function (Eq.13) in the synthesis step.The margin  between the different areas (Fig. 6) is chosen to 4 o .Five solutions are performed according to the five parameter choice methods 1) to 5) listed in Table 2.For each solution, the estimated coefficients   , the calculated disturbing potential   as well as its difference to the validation data are plotted in Fig. 7.The results for the three criteria measures from above are listed in Table 4.
The correlations between the modelled gravity field   and the validation data   for all parameter choice methods are rather satisfying.However, the CM 1 based methods give slightly smaller correlations than the others.The last three methods are comparable and provide better results than the CM 1 based methods with respect to both RMS values and correlations.The lowest RMS error is obtained from the 'VCE + L-curve method' which is 3.61 m 2 s 2 ⁄ .This method also delivers the highest correlation between the estimated coefficient   and the validation data."L-curve method + VCE' gives the second-best RMS value which is 3.64 m 2 s 2 ⁄ , followed by 'VCE based on CM 2'.The same rank applies to the correlation between the estimated coefficient   and the validation data.The largest RMS value and smallest correlation are both obtained from parameter choice methods based on CM 1.Among the methods based on CM 1, the L-curve method provides smaller RMS value compare to VCE, but smaller correlation factor as well.
It is worth clarifying that the solution obtained from the 'L-curve method + VCE' is not unique.Due to the fact that the regularization parameter  L−curve is fixed during VCE, the results change when it refers to different observation techniques.
Here, two solutions are obtained by setting GRACE data and the terrestrial data as the reference observation technique, respectively.In this study case, the better result (smaller RMS value and larger correlation) is the one when choosing the terrestrial data as the reference.Generally, a better solution is obtained when the terrestrial data are chosen as the reference, and the results for 'L-curve method + VCE' listed in this paper are always the better ones.

Study case F
In case F, the maximum degree in the expansion in terms of SRBF is set to   = 1050 for the Shannon function (Eq.12) in the analysis step as well as  1 = 900 and   = 1100 for the Blackman function (Eq.13) in the synthesis step.The margin  between the different areas (Fig. 6) is chosen to 2 o .For each method, the estimated coefficients   , the calculated disturbing potential   as well as its difference to the validation data are plotted in Fig. 8.The results for the three criteria measures from above are listed in Table 5.
Compared to the study case A, the results in the study case F shows a general improvement, in terms of all three criteria.The correlations between the estimated coefficients and the validation data are promising, which can be observed in Fig. 8.It shows that the estimated coefficients reflect the energy of the recovered gravity signal rather well.'VCE + L-curve method' still provides the smallest RMS error 0.78 m 2 s 2 ⁄ as well as the highest correlations.'L-curve method + VCE' gives RMS value 0.80 m 2 s 2 ⁄ , followed by 'VCE based on CM 2' with 0.82 m 2 s 2 ⁄ .The same rank applies to the correlation between the estimated coefficient   and the validation data.
The largest RMS values and the smallest correlations are still both obtained from parameter choice methods based on CM 1.
The differences between CM 1 and CM 2 based methods are considerable.It could be caused by the large variation between the spectral resolution of GRACE, GOCE, terrestrial and airborne data.Therefore, giving each observation technique a relative weight in the combination might help to provide better results.Among the two methods based on CM 1, the L-curve method provides worse RMS value in comparison to VCE but, a much more significant correlation between   and   .

Results of all study cases
Table 6 lists the RMS values of each method obtained from all the study cases.Table 7 lists the correlations between the estimated coefficients   and the validation data   of each method obtained from all the study cases.The best results (smallest RMS value and largest correlation) in each study case are bold-typed and the second bests are italicized.
Based on these results, the following conclusions can be drawn: 1. From the five parameter choice methods considered here, 'VCE + L-curve method' performs the best, and always gives the smallest RMS error as well as the largest correlation.
2. 'L-curve method + VCE' and 'VCE based on CM 2' also show a good and stable performance, especially when the spectral resolution of the combined data sets differs.'L-curve method + VCE' generally outperforms 'VCE based on CM 2' slightly.
3. The results in terms of RMS value and correlation are consistent, i.e., in most cases, the regularization parameter choice method which gives the smaller RMS error also delivers the larger correlation.
4. Generally, results provided by CM 1 based methods are not as good as the others.The larger the spectral resolution between each observation technique is, the worse these methods perform (e.g.case E and F).However, for case D where the combined data sets have similar resolution, 'the L-curve method based on CM 1' perform even better than 'L-curve method + VCE' and 'VCE based on CM 2'. that the larger the spectral resolution differs between each type of observation technique, the better 'VCE based on CM 2' and 'L-curve method + VCE' perform compared to the L-curve method or VCE based on CM 1; which is reasonable.
We also carried out the same investigation using the CuP function for comparison scenario to test the dependency of our new methods against the type of SRBFs used.In this scenario, the performance of the L-curve method is slightly lower than the one obtained using the Shannon function.This behavior is consistent with the results from the literature.Our newly proposed 'VCE + L-curve method' still provides the best results in all the six cases and the 'L-curve method + VCE' is also performing well.Thus, we are able to conclude that our new methods provide good regularization results for different observation combinations and are stable regardless of the type of SRBF used.From our investigation, we conclude that 'VCE + L-curve method' is the best choice among those five methods for the determination of the regularization parameter.
In future, a primary concern is to apply the newly devised methods using more types of SRBFs, so that the performance of different SRBFs can be compared while making sure that the differences in results are not coming from the regularization method.In addition, further investigations are planned for using real observations instead of simulated data sets.
Competing interests.The authors declare that they have no competing interests.

(
deterministic part) (stochastic part) can be set up.In the deterministic part,  = [( 1 ), ( 2 ), … , (  )]  is the  × 1 vector of the observation errors and  = [(,   )] is the  ×  design matrix containing the corresponding basis functions.In the stochastic part, () is the  ×  covariance matrix of the observation vector f with  2 the unknown variance factor and P the given positive definite weight matrix.Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.
Finally, two new proposed methods are presented as a combination of VCE and the L-curve method.Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.

2.
Terrestrial data: provided on a regular grid with two different resolutions, one over an area of 20 o × 30 o (latitude × longitude) with a grid spacing of 30' (blue shaded area) and the other one over an inner area of 6 o × 10 o with a grid spacing of 5' (yellow shaded area).3. Airborne data: provided on two different flight tracks, one over the Adriatic Sea (magenta shaded area) and the other one over Corsica connecting Southern Europe with Northern Africa (cyan shaded area).Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.

5 .
The results also indicate that when the Shannon function is used for analysis and the same combination model CM 1 is used, the L-curve method generally outperforms VCE (case A, B, C, D, E regarding the RMS value and case B, C, D, E, F regarding the correlation).Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 1 : 5 Figure 2 :
Figure 1: Different SRBFs in the spatial domain (top, ordinate values are normed to 1) and the spectral domain (bottom) for   = .

Figure 8 : 5 Solid
Figure 8: The estimated coefficients   (left column), the recovered disturbing potential (mid column) and the differences w.r.t the validation data (right column) for study case F. The results are obtained using regularization methods 'the L-curve method based on CM 1' (first row), 'VCE based on CM 1' (second row), 'VCE based on CM 2' (third row), 'VCE + L-curve method' (forth row) and 'L-curve method + VCE' (fifth row); see also Table 5 for numerical results.5 observed the tensor ∆ of the gravity gradients   with ,  ∈ {, , }, i.e. all second-order derivatives

Table 1 : Adapted basis functions for each type of observation observation adapted basis function B
∞ =0

Table 6 : RMS value results of each method for different study cases
Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.

Table 9 : Correlations between the estimated coefficients and the validation data of each method using the CuP function (as a comparison scenario w.r.t Table 7)
Solid Earth Discuss., https://doi.org/10.5194/se-2019-60Manuscript under review for journal Solid Earth Discussion started: 5 April 2019 c Author(s) 2019.CC BY 4.0 License.