Abstract
We solve a geodetic inverse problem for the determination of a distribution of point masses (characterized by their intensities and positions), such that the potential generated by them best approximates a given potential field. On the whole Earth, a potential function is usually expressed in terms of spherical harmonics which are basis functions with global support. The identification of the two potentials is done by solving a least-squares problem. When only a limited area of the Earth is studied, the estimation of the point-mass parameters by means of spherical harmonics is prone to error, since they are no longer orthogonal over a partial domain of the sphere. The point-mass determination problem on a limited region is treated by the construction of a local spherical harmonic basis that is orthogonal over the specified limited domain of the sphere. We propose an iterative algorithm for the numerical solution of the local point mass determination problem and give some results on the robustness of this reconstruction process. We also study the stability of this problem against added noise. Some numerical tests are presented and commented.
1 Introduction
One of the main problems in physical geodesy is the determination of the geoid: the equipotential surface for the gravitational potential which coincides with the mean ocean surface. During the past few decades there has been a growing interest in developing methods, with increasing levels of accuracy, for the determination of the high-resolution geoid (and/or regional geoid). Among the methods widely used by the geodetic community for the gravimetric determination of the geoid we mention the Stokes’ integration approach,[Citation1] the least-squares collocation approach,[Citation2, Citation3] and the point-mass inversion approach.[Citation4, Citation5]
The point-mass approach is based on the simple idea that, in principle, it is possible to find an equivalent distribution of discrete point masses placed below the Earth’s surface for which the gravitational potential field generated by them closely approximates the true gravity field in the region of interest. This approach can also be considered as a model reduction for which the representation by the inverse-distance Newtonian potential functions can be seen as simple alternative to the classical representation by spherical harmonics.
There are many methods based on the point-mass inversion approach which differ by the procedure used for the determination of the positions and intensities of the point masses. While some of these methods place the point masses on a fixed regular grid beneath the Earth’s surface and try to find their depths and intensities,[Citation6] other methods use freely positioned point masses and thus seek full positions and intensities of the point masses.[Citation4] We also note that the strategy used for the optimization of the positions and intensities differs from one method to another.
In this paper, we are interested in using the theory of localized functions to construct a distribution of buried point masses in such a way that the potential generated by them best fits, in the sense of least squares, a given potential gravitational field in a specified region of the Earth. In fact, in the case where the region of interest is the whole Earth, the orthonormal spherical harmonics basis is used to expand the given and synthesized potentials in order to formulate the problem of the determination of the positions and intensities of buried masses as an inverse problem. On the other hand, when only a limited region of the Earth is considered, the spherical harmonics basis is no longer orthogonal. Consequently, the spherical harmonic basis for the determination of the local point-mass problem is not appropriate and hence one has to find another basis which is localized in the region of interest. Once such a basis has been found, by expanding the given and synthesized potentials over this basis, the problem of the determination of the positions and intensities of buried masses can then be formulated as an inverse problem.
We note that the question of finding a local basis was first studied by Slepian, Pollak and Landau for an interval of the real line and for a rectangular region of the plane.[Citation7–Citation12] They discovered by serendipity a basis of functions (now called Slepian functions) whose energies are concentrated in the considered region. This procedure was very useful in several domains of applied mathematics and physics, notably geophysics,[Citation13, Citation14] cosmology [Citation15] and image processing.[Citation16, Citation17]
The problem of localization (or concentration) on a limited region of the sphere has been an active area of research for many years. We particularly mention the pioneering works of Simons and his coauthors [Citation18–Citation21] on spherical Slepian functions and their applications in geodesy. Other authors use different approaches to this problem. For example, Michel and his collaborators [Citation22–Citation24] use multiresolution splines to deal with the problem of localization. We also mention zooming-in approximation by locally supported wavelets.[Citation25]
The remaining of this paper is structured as follows. In Section 2, we start by gathering some background materials on spherical harmonics that will be used in the sequel. Then, we give the expansions of the given (observed) gravitational potential and the potential generated by the distribution of point masses (synthesized potential) in terms of spherical harmonics. In Section 3, we first formulate the problem of equivalent point-mass determination for the whole Earth as a minimization problem. Then, we use Slepian analysis to formulate the problem of the determination of equivalent point masses which generate a gravitational potential that closely approximates the given gravitational potential on a specified region. Some numerical results are presented and commented in Section 5. Finally, some conclusions and perspectives are given in Section 6.
2 The gravitational and synthesized potentials
2.1 Spherical harmonics
Let denote the unit sphere in and let denote a generic point on the sphere where is the co-latitude and is the longitude. On the space of square-integrable real functions on , we have the inner product defined by1 1 where denotes the surface element on the unit sphere . The associated norm is given by2 2 We recall that any square-integrable function on the unit sphere can be expanded in terms of spherical harmonics as [Citation26]3 3 whereis the spherical harmonic coefficient of associated with , the spherical harmonic of degree and order which is given by4 4 Here, is a normalized version of , the Legendre’s associated function of degree and order given by5 5 The spherical harmonics functions are orthonormal with respect to the inner product (Equation11 1 ), i.e. [Citation2]Finally, we recall that for two points and on the unit sphere, by the addition theorem of spherical harmonics we have the identity [Citation18]6 6 where is the spherical distance between and given by7 7
2.2 The gravitational potential in terms of spherical harmonics
The gravitational potential satisfies [Citation2]where is the Earth, is the mass density and is the universal gravitational constant. Furthermore, across the Earth’s surface we have continuity of both the potential and its normal derivative.
Since the gravitational potential is a harmonic function outside the Earth, for a point with spherical coordinates that lies outside the Earth we have the following expansion in spherical harmonics [Citation2]8 8 Here is the total mass of the Earth, is the equatorial radius of the Earth and is the geopotential coefficient of degree and order . The values of and depend on the choice of the geopotential model (e.g. EGM96, OSU91, etc.) at the considered zone. Note that the first term in the summation, i.e. for is which represents the potential generated by a homogeneous body of mass located at the origin.
Geopotential dedicated Earth-orbiting satellites, like CHAMP (CHAllenging Minisatellite Payload) and GRACE (Gravity Recovery and Climate Experiment), provide data that are incorporated to determine the spherical harmonics coefficients of the Earth gravitational potential. In the EGM96 model, the coefficients in this expansion are determined up to the degree , see e.g. [Citation27, Citation28]. The higher resolution model, EGM08, provides the coefficients up to the degree , see [Citation29].
2.3 The synthesized potential in terms of spherical harmonics
The inverse gravimetric problem of finding the mass density inside the Earth from information of the gravitational potential outside the Earth is not well posed in the sense of Hadamard. In fact, this problem does not have a unique solution as the associated Fredholm integral operator is of the first kind and has an infinite dimensional null space.[Citation22] This null space corresponds to the space of functions which are -orthogonal to the space, , of harmonic functions in . Furthermore, even when restricting this operator to , it still has a discontinuous inverse and thus stability is not guaranteed.
In the point-mass approach, the true mass density is approximated by a finite distribution of point masseswhere and are the mass and location of the point mass and is the Dirac distribution.
We here recall the following identifiability result [Citation30]:
Theorem 2.1
Let , be the solutions of the problemsIf on a subset of with non-empty interior then and (up to a permutation of the indices) and for .
It is interesting to note that this result does not hold for an infinite distribution of point masses. Lipschitz stability results for the problem of determining the intensities of the point-masses and their locations from measurements of the generated potential are given in [Citation31].
Similar to what we did for the potential generated with the mass density , we will expand the gravitational potential generated by the discrete mass density in spherical harmonics. Let be a distribution of points located at and with masses . Recall that the potential at a spacial point with coordinates generated by a point mass is given by9 9 where is the Euclidean distance between the two points and . If we denote by the spherical distance between and which is obtained by taking in (Equation77 7 ), then10 10 As the point masses are buried, we have , , and as a consequence, the inverse of the distance between and , admits a convergent Legendre series expansion. Indeed, if we set , then11 11 and therefore can be expanded into power series with respect to [Citation2]:12 12 Then, it follows from (Equation66 6 ) that13 13 The expression (Equation99 9 ) of the potential at the spacial point becomesThe potential at generated by the point masses is the sum of the potentials generated by the considered elementary buried masses14 14 On a selected sphere (i.e. for a fixed value of ), we seek the masses and positions of the buried points, so that the resulting generated potential best approximates the gravitational potential for the considered region. This involves solving an inverse problem which we formulate in the following section.
3 Formulation as inverse problems
3.1 Inverse problem on the whole Earth
The problem of determining the potential generated by a finite distribution of point masses of given positions and intensities can be solved by considering the expression of the elementary potential (Equation99 9 ), we gave in (Equation1414 14 ). Thus, the potential caused by a distribution of point masses at a spacial point with spherical coordinates , is the sum of the elementary potentials generated by all the point masses. Such a problem is called the forward problem.
The inverse problem is then to find the masses and positions of the point masses that generate a gravitational potential close to a given one, i.e. for which the difference between the given gravitational potential and the one generated by these masses is minimized. Eventually, some constraints on those unknowns need to be imposed to make the computations more stable. If we assume that the point masses co-latitudes and longitudes are distributed on a fixed regular grid on the unit sphere then the unknown parameters are the masses and radii . The approach we use here, is based on the expansion of the two potentials (i.e. the modelled and the predicted) in the same orthogonal basis. In the following sections, we describe the orthogonal basis to be used which depends on the region of interest: the whole Earth or on a limited region of the Earth.
To find the distribution of point masses which generates a potential field that best fits (in the sense of least squares) a given potential field over the whole Earth we introduce the real-valued function defined on byand we consider the following minimization problem15 15 That is we are interested in finding the masses and positions of the points so that the function attains its minimum.
Using the spherical harmonic expansions (Equation88 8 ) and (Equation1414 14 ) and the orthonormality of the spherical harmonics on , this minimization problem can be written as16 16 where is a vector with infinite numbers of rows and is a matrix function of with an infinite number of rows and columns. The entries of and arewhere . We shall assume that the positions of the point masses are such that the matrix has full column rank. Otherwise, we can arrange by adding/removing a point mass so that satisfies this hypothesis.
To solve the minimization problem (Equation1616 16 ), we use an alternating minimization algorithm which consists in an iterative procedure that solves alternatively a linear problem for the masses (which are assumed non-negative) and a non-linear problem for the positions . The main steps of this algorithm are as follows:
(i) | Given fixed positions for the masses, we solve17 17 This least-squares problem is linear but with non-negative constraints for the masses. This problem is solved, for example, by using the non-negative least squares algorithm (NNLS) of Lawson and Hanson [Citation32]. We denote by the solution of the problem (Equation1717 17 ). | ||||
(ii) | Given the masses , we solve (iteratively) for the problem18 18 The first-order optimality condition for this problem is the vanishing of the gradient of its objective function, i.e. Of course, to initialize the process, an initial position is set for all the point masses. |
3.2 Point-mass determination on a limited region of the Earth
The fact that, for a scalar product defined on a limited region of the Earth defined by19 19 the spherical harmonics basis is no longer orthogonal for the -inner product, makes the computations of the entries of the matrix and vector more complicated. A new basis of functions that are concentrated on the region of interest is now defined.
3.2.1 Slepian basis of localized functions
As stated earlier, a square-integrable function on can be represented asWe say that is band limited with band width if for all , i.e. it has the finite sum representation20 20 Following the pioneering work of Slepian, to find a basis of functions that are localized on a given region of the unit sphere, we seek functions such that the ratio21 21 between their energy over and their energy over is maximized. When the set has positive measure and is not equal to , we have for any .
This maximization problem is equivalent to an eigenvalue and eigenfunction problem. Indeed, if is band-limited with band width then (Equation2121 21 ) reduces to22 22 whereIf we denote by the matrix with coefficients and by the -spherical harmonic coefficient vector, then , and the maximization problem becomes the following algebraic eigenvalue problem23 23 Since the matrix is real, symmetric and positive definite, its eigenvalues , , are all positive. We shall assume that they are ordered such that . Furthermore, the symmetry of guarantees that the eigenvectors are real and orthonormal. Consequently,24 24 Every eigenvector , defines an associated spatially band-limited eigenfunction of the form (Equation2020 20 ). We remark that these functions are at the same time orthonormal with respect to the inner product (Equation11 1 ) and orthogonal with respect to the inner product (Equation1919 19 ), i.e.25 25 Now, if we multiply the eigenvalue Equation (Equation2323 23 ) by and sum over all and we deduce that the eigenfunction satisfies the homogeneous Fredholm integral equation of the second kindwhere is the integral operator associated with the kernelwhich is symmetric and depends only on the spherical distance, , between and and given in (Equation77 7 ).
3.2.2 Point-mass determination using the Slepian basis
The problem of point-mass determination on a limited region of the Earth is solved in an analogous manner to that we used for the one formulated on the whole Earth but using the Slepian basis functions instead of the spherical harmonics. The main steps are essentially:
(i) | Initialization We choose the geographic region to study and we compute its associated Slepian basis. Only the basis vectors associated with eigenvalues are kept. We also define the grid of the point-masses (see Section 3.2.1). | ||||||||||||||||||||||||||||||||||||||||
(ii) | Construction of the right-hand side We compute the coefficients of the Earth gravitational potential over the considered region relative to the spherical harmonics basis. We note that we used the geoid model EGM96 for which the potential decomposition on the spherical harmonics basis is truncated to degree . The spherical harmonic coefficients , , for this model are available from [Citation33]. The right-hand side of our problem is the vector of the spherical harmonics components of the gravitational potential (see Section 3.1). Then, we project this vector on the orthogonal Slepian basis. To start the iterative process, we fix the initial depth of all point masses to . | ||||||||||||||||||||||||||||||||||||||||
(iii) | Alternating minimization Set and repeat (a) and (b) until convergence:
|
4 Determination of a distribution of point masses on a grid
4.1 Minimum separation distance between point masses
To test the robustness of the results of the reconstruction process it is instructive to start by analysing the sensibility of the identification of two point masses (possibly with distinct masses) with respect to their separating distance and the ratio . In our computations, we take which is the radius of an observation sphere concentric with the Earth, and we consider two point masses and located on a regular grid buried at a constant depth . Let us denote by the potential generated by both point masses at a spacial point located on an observation sphere which has the same centre as the Earth and a radius . We restrict our analysis to the case where the point belongs to the planar section containing , and the geocenter. We then havewhere, as stated in (Equation1010 10 ), for .
As depends only on the spherical distance between and , we assume, without loss of generality, that the two point masses are on the same meridian, i.e. with fixed . Therefore, from (Equation1111 11 ) we obtainFurthermore, by shifting the latitude origin to the midpoint of the mass locations, and by denoting (see Figure ), we get andwhereFor , we consider the following two experiments. In the first one, the two point masses are assumed to have the same (fixed) longitude (i.e. they are positioned on the same meridian) and have different latitudes which vary in the interval . In the second one, the point masses are assumed to have the same (fixed) latitude (i.e. they are placed on the same parallel) and have different longitudes which vary in .
In Figure we display, for equal and distinct masses, the plots of the composite potential over the arc intersecting the plane containing the two point masses and the geocenter, with the observation sphere.
Figure shows that as the distance between the point masses decreases the curve of the potential generated by two point masses changes from a two-humped form to a one-humped form. This result holds true in the case of point masses with equal and different mass values. Furthermore, there is a critical distance below which the potential generated by the two point masses appears as it is the potential generated by a single point mass located at the midpoint and with a mass equal to the sum of the masses of the two points.
This phenomenon is observed in other fields such as optics and radar signals generated by very close point sources. The Rayleigh criterion [Citation34] is a widely accepted criterion for finding the minimum distance for which the point sources can be resolved from the signals. Consequently, for a chosen ratio , the minimal distance that permits us to distinguish point masses, can be determined to a prescribed precision using Rayleigh criterion.
In Figure , we plot the second derivative with respect to the spherical distance between the point masses. We remark that when the distance between the two point masses (in the case of equal and different masses) is less than , the second derivative is negative and the potential attains a local maximum at the midpoint. In this case, the identification process becomes unstable.
4.2 Stability of the reconstruction process
To check the stability of point-mass determination algorithm, we proceed as follows. First, we compute the vector of the problem (Equation1616 16 ) that corresponds to the potential generated by a given distribution of point masses. Second, we add a normally distributed noise vector to the vector and we determine, for several values of the noise level, the masses and positions that result from the minimization problem (Equation1616 16 ) with the unperturbed vector and the perturbed one , where is a vector of random normally distributed real values and is a scalar that represents the noise level times the magnitude of . Our experiments, which are performed on a regular grid with grid size , consisted of the following two situations: a single point mass and a pair of point masses (Figure ).
In the single point mass case, the vector contains the coefficients of the potential generated by . This experiment consists in recovering the point mass for data with and without noise. We remark that without noise we are able to recover the point mass accurately. As we add noise we observe that the recovered masses consists of a principal mass located at , and other masses located nearby and with much smaller mass than . We note that the sum of all identified masses is equal to .
In the two point masses case, the vector contains the coefficients of the potential generated by and . Without added noise we are able to accurately recover the masses. For a of noise level, we find two principal masses located at the same positions than our initial masses and some other masses nearby. However, when the noise level is , we are not able to recover the masses.
In Table , we show the proportion of retrieved principal masses and the initial masses in both cases (the value in the table means that we recover the initial point masses accurately). Values lesser than are explained by the presence of other masses around the principal ones. Table shows that when the noise level becomes high, our procedure fails to recover the initial masses.
5 Reconstruction results for a specified geographical region
The numerical experiments presented in this section are for a domain which corresponds to the spherical cap centred at the geographic point with coordinates , and with an epicentral distance of , see Figure .
The local basis associated to the domain is computed using the methodology described in Section 3.2.1 for a band-width . We have used Frederik J. Simons’ Matlab toolbox ‘Spatiospectral concentration on a sphere’.Footnote1
In the Figure , we present some eigenfunctions of the operator in decreasing order of significance. Eigenfunctions that are well concentrated within the domain are the ones associated with the eigenvalues close to 1. Whereas those that are weakly concentrated are those associated with small eigenvalues. Accordingly, the selected basis functions are chosen to be located inside the area under study and in an appropriate neighbourhood. It should be noted that the number of significant eigenfunctions increases as the bandlimit increases. Figure shows this behaviour for several bandlimit values.
We remind that our problem consists in finding point masses , defined by their masses, , and depths, , on a fixed buried grid . Here is the Earth’s mean radius. As stated in Section 3.2.2, we use a minimization algorithm in which we alternate (until convergence) between the two minimization problems (Equation1717 17 ) and (Equation1818 18 ). We start with a given , and for we seek solution of the non-negative minimization problem26 26 and solution of the non-linear minimization27 27 We note that (Equation2626 26 ) is a constrained minimization problem. To solve it we use the maximum entropy regularization method in proposed by Hansen [Citation35]. In this method, we need to solve the unconstrained problemwhere is the regularization parameter. The entropy terms ensure that the masses are positive. It has been shown that this regularization procedure is comparable to the Tikhonov regularization with non-negativity constraints.[Citation35] Selecting an appropriate regularization parameter is crucial and must be estimated from the given data. In our computation, we use the L-curve method [Citation36] for estimating the optimal value for the regularization parameter.
In the following experiments the depths of all points are initialized at . Our grid (Figure ) is chosen to be a regular grid formed by four concentric circles. On each circle there are 19 grid points. This choice satisfies our stability condition related to the minimal distance separating points for which the elementary potentials are distinguishable (see Section 4.1).
In Figure we present, at different values of the iteration number , the L-curve that yields the good choice of the regularization parameter at the iteration . As it can be seen, all these curves clearly exhibit L-shaped graphs plotting the solution norm against the residual norm.
To assess the accuracy of our reconstruction procedure we define the pointwise relative error as:where and are, respectively, the real and the point-mass gravitational potentials. Figure shows the relative error on the gravitational potential.
Our algorithm yields a point mass at each point of the regular grid, however only few have significant mass. We remark that the relative error values are less than on a small part of the domain, and less than on the rest of the domain.
6 Conclusion
In this paper, we have considered an important problem in geodesy, namely the inverse problem of point-mass determination relative to a limited region of the Earth. We have tackled this problem by using a localized Slepian basis. We have presented an iterative algorithm to numerically solve this problem. We have established that there is a minimum distance between the point masses that must be respected so that the identification process of the point masses works. We have computed this minimum distance of approach explicitly. We have presented several numerical examples and have studied the effects of added noise on the reconstruction process. The numerical tests of our point-mass determination algorithm on a specific geographical region show that the pointwise relative errors between the given potential and the one generated with the obtained point masses are acceptable.
Acknowledgments
The authors would like to thank the anonymous referees for their valuable comments and suggestions which helped them to improve an earlier version of this paper.
Notes
1 Available online from http://geoweb.princeton.edu/people/simons/software.html.
References
- Featherstone WE, Olliver JG. A method to validate gravimetric-geoid computation software based on Stokes’s integral formula. J. Geodesy. 1997;71:571–576.
- Hofmann-Wellenhof B, Moritz H. Physical geodesy. New York (NY): Springer; 2005.
- Tscherning CC. Collocation and least squares methods as a tool for handling gravity field dependent data obtained through space research techniques. J. Geodesy. 1978;52:199–212.
- Antunes C, Pail R, Catalão J. Point mass method applied to the regional gravimetric determination of the geoid. Stud. Geophys. Geod. 2003;47:495–509.
- Baur O, Sneeuw N. Assessing Greenland ice mass loss by means of point-mass modelling: a viable methodology. J. Geodesy. 2011;85:607–615.
- Lehmann R. Gravity field approximation using point masses in free depths. Bull. Int. Ass. Geodesy (IAG). 1995;1:129–140.
- Landau HJ, Pollak HO. Prolate spheroidal wave functions, Fourier analysis and uncertainty II. Bell Syst. Tech. J. 1961;40:65–84.
- Landau HJ, Pollak HO. Prolate spheroidal wave functions, Fourier analysis and uncertainty III. Bell Syst. Tech. J. 1961;41:1295–336.
- Slepian D. Prolate spheroidal wave functions, Fourier analysis and uncertainty IV. Bell Syst. Tech. J. 1964;43:3009–58.
- Slepian D. Prolate spheroidal wave functions, Fourier analysis and uncertainty V. Bell Syst. Tech. J. 1978;57:1371–430.
- Slepian D. Some comments on Fourier analysis, uncertainty and modeling. SIAM Rev. 1983;25:379–393.
- Slepian D, Pollak HO. Prolate spheroidal wave functions, Fourier analysis and uncertainty I. Bell Syst. Tech. J. 1961;40:43–64.
- Albertella A, Sansò F, Sneeuw N. Band-limited functions on a bounded spherical domain: the Slepian problem on the sphere. J. Geodesy. 1999;73:436–447.
- Lesur V. Introducing localized constraints in global geomagnetic field modelling. Earth Planets Space. 2006;58:477–483.
- Gòrski KM, Hivon E, Banday AJ, Wandelt BD, Hansen FK, Reinecke M, Bantelman M. HEALPix: a framework for high-resolution discretization and fast analysis of data distributed on the sphere. Astroph. J. 2005;622:759–771.
- Chung MK, Dalton KM, Evans AC, Davidson RJ. Tensor-based cortical surface morphometry via weighted spherical harmonic representation. IEEE Trans. Med. Imag. 2008;27:1143–1151.
- Maniar H. Mitra PP. The concentration problem for vector fields. Int. J. Bioelectromagn. 2004;7:142–145.
- Dahlen FA, Simons FJ. Spectral estimation on a sphere in geophysics and cosmology. Geophys. J. Int. 2008;3:774–807.
- Simons FJ. Slepian functions and their use in signal estimation and spectral analysis. In: Freeden W, Nashed MZ, Sonar T,editors. Handbook of geomathematics. Berlin Heidelberg: Springer-Verlag; 2010. Chapter 30; p. 891–923.
- Slobbe DC, Simons FJ, Klees R. The spherical Slepian basis as a means to obtain spectral consistency between mean sea level and the geoid. J. Geodesy. 2012;86:609–628.
- Simons FJ, Dahlen FA, Wieczorek MA. Spatiospectral concentration on a sphere. SIAM Rev. 2006;48:504–536.
- Michel V, Fokas AS. A unified approach to various techniques for the non-uniqueness of the inverse gravimetric problem and wavelet-based methods. Inverse Probl. 2008;24:1–25.
- Berkel P, Fischer D, Michel V. Spline multiresolution and numerical results for joint gravitation and normal-mode inversion with an outlook on sparse regularisation. Int. J. Geomath. 2011;1:167–204.
- Michel V, Wolf K. Numerical aspects of a spline-based multiresolution recovery of the harmonic mass density out of gravity functionals. Geophys. J. Int. 2008;173:1–16.
- Freeden W, Fehlinger T, Klug M, Mathar D, Wolf K. Classical globally reflected gravity field determination in modern locally oriented multiscale framework. J. Geodesy. 2009;83:1171–1191.
- Dahlen FA, Tromp J. Theoretical global seismology. Princeton (NJ): Princeton University Press; 1998.
- Rapp RH, Pavlis NK. The development and analysis of geopotential coefficient models to spherical harmonic degree. J. Geophys. Res. B. 1995;360:21885–21911.
- Tapley B, Ries J, Bettadpur S, Chambers D, Cheng M, Condi F, Gunter B,Kang Z, Nagel P, Pastor R, Pekker T, Pool S, Wang F. GGM02 - An improved earth gravity field model from GRACE. J. Geodesy. 2005;8:467–478.
- National Geospatial-Intelligence Agency. EGM2008 official webpage. Available from: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
- El Badia A, Ha-Duon T. An inverse source problem in potential analysis. Inverse Probl. 2000;16:651–663.
- Vessella S. Locations and strengths of point sources: stability estimates. Inverse Probl. 1992;8:911–917.
- Lawson CL, Hanson RJ. Solving least squares problems. Englewood Kliffs (NJ): Prentice-Hall; 1974.
- National Geospatial-Intelligence Agency. EGM96 official webpage. Available from: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
- Bochkov GN, Gavrilin AT, Gorokhov KV. Generalized Rayleigh criterion of two-point resolution. Radiophys. Quant. Elec. 2005;5:394–397.
- Hansen PC. Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems. DK-2800 Lyngby, Denmark: Technical University of Denmark; 2007.
- Hansen P. Analysis of discrete Ill-posed problems by means of the L-curve. SIAM Rev. 1992;34:561–580.