Abstract: Inverse problems which originate from potential field data may be solved using a combination of signal norm and data noise minimalisation. Since both the gravity and the magnetic potential fulfil the Laplace partial differential equation in outer space, they may be regarded as being elements of a reproducing kernel Hilbert Space (RKHS). This makes it possible to use the collocation method in order to construct a model of the outer potential field consistent with observational data of many different kinds.

The practical selection of the reproducing kernel may be done so that it represents the auto-covariance of the potential. This has the advantage, that the predictions and the error estimates may be given a statistical interpretation. Furthermore the square of the signal norm and the noise variance are balanced in a natural manner.

A similar model may be established if the source of the potential field may be modelled as a so-called quasi-harmonic function or as a simple density contrast, contingently with inversely correlated deeper located compensating masses. If the density contrast for example is the condensated ocean-bottom topography, empirical covariance functions may be determined from a-priori depth data. This may be used to select the appropriate reproducing kernel for the density.

1. Introduction.

Envisage the situation that we have a set of data associated with points in space

(generally at the Earths surface) which represent measurements of some functional applied on a real function. We will denote the function f. It could be the gravity or the magnetic potential, or the source producing one of the potentials. The data could be point or mean gravity or magnetic anomalies.

Inverse problems arise either when we want to determine the source from data related to the potential, or when we want to construct a representation of the potential valid at the Earth's surface from data collected in space- or aircrafts.

This type of problems is well-known in all branches of science, and solutions are known under various labels: collocation, regularization, optimal estimation, objective mapping etc.

The methods are normally classified as analytic or approximation methods or as statistical estimation methods. The aim of this paper is to show that the methods to a high degree are equivalent, and that we from the analysis of one method may understand better and hence modify another method.

We will only consider problems related to the interpretation of potential fields. They are nice, and give us the possibility for applying well developed mathematical structures.

2. Mathematical foundation.

Suppose the function f is an element of a finite or infinite dimensional linear vector space of functions from 3-dimensional space to the real line. Then the observations y_i are related to f and the data noise vector v through the equation

y and v are n-vectors, and A is an operator from the function space to n-dimensional real space. If the vector space is finite (m) dimensional, and A is a linear operator then A may be represented by a m x n matrix.

Example 1. Suppose the vector space is spanned by m indicator functions I_i, equal to one in a rectangular block (or sphere) and zero outside. Suppose the blocks do not overlap. Then a density model may be expressed as

where a_i is the density of the i'th block. If our observations are gravity disturbances g, then

where is the volume element in R³, P is the observation point, Q is the integration point, G the Newtonian gravity constant, z the vertical coordinate and P-Q is the distance between P and Q. (This is just Newtons inverse quare law applied).

In this case

where P_i is the point where the gravity disturbance is observed. We see that A_ij simply is the attraction in the i'th point of the j'th box or sphere.

Example 2. Suppose f is a periodic function on the interval from - to . Then f may be developed in the usual Fourier series, using the base functions cos(nx) and sin(nx) where n is a positive integer. The function may be used to represent the compensated density on a line at a given depth. In a plane one simply interchange x for the complex variable. Here the Fourier transform can be used to compute the gravity response in a fast manner, taking advantage of FFT-procedures if the data spacing fulfils some simple criteria, see e.g. Knudsen(1993).

Example 3. Suppose f is a harmonic function outside the Earth. The f may be approximated by the sum of an infinite series in solid spherical harmonics, f_ij, so that

where are the fully normalized Legendre functions, is the latitude, is the longitude, a is the Earth's mean radius and r the distance from zero. Then

In all the examples we may equip the spaces with an inner product. If this is the simple L₂ inner-product (integral over the product of the two functions) then the base functions I_i, cos(nx), sin(nx) and f_ij form orthonormal bases. The spaces become Hilbert spaces, if we in example 2 and 3 exclude all functions with infinite norm (the square sum of the a_ij is infinite).

The L₂-norm is not the only possible. We may use norms which in example 1 take into account the first and second order differences between the density values in between the blocks, and in example 2 and 3, we may integrate over sums of derivatives of the functions. Generally the consequence is that the base functions stay orthogonal, but they are scaled with a certain positive constant, b_i. The norm of the functions will then become the square sum of the fourier-coefficients a_i(j) divided by b_i. As a consequence, some more functions will be included or some functions will be excluded from the spaces related to example 2 or 3.

For the spaces we deal with in the examples we have or may construct a so-called reproducing kernel, K(P,Q). In this case we denote the space a reproducing kernel Hilbert space (RKHS).

If the scalar (inner) product is denoted (.,.) then

The function is reproduced by the kernel. The kernel will be equal to the product sum of the base functions evaluated in the points P, Q, respectively,

where we have used the orthonormality of the base functions.

Example 4. The spaces in the examples 1-3 are or may become RKHS. Using the L₂-norm the indicator functions divided by the square-root of their volume will form an orthonormal system. The sin(nx), cos(nx) functions are orthonormal and so are the solid spherical harmonics, when the L₂-norm integration area is the surface of the mean Earth Sphere.

The reproducing kernels becomes

for example 1, where v_i is the volume of the i'th block,

for example 2, where we see that the kernel is only dependent of x-y. However, the kernel is not a correct reproducing kernel, because the function for x = y becomes infinite. This is a consequence of the fact that a space with a L₂-norm may contain functions with spikes. In order to become a RKHS, a stronger norm must be introduced, for example also integrating over the first and second derivatives. This will have as a consequence that each cos-term is multiplied by a factor proportional to a quantity which decreases towards zero faster than i^-1.

For example 3 we have,

where P_i are the usual Legendre polynomials, r' is the distance of Q from 0 and is the spherical distance between P and Q. The expression may be further simplified using that the inverse-distance may be expanded in a series in Legendre functions, see Moritz (1980). The important fact here is, that the kernel now only depends on r, r' and . Here again we have the problem that the kernel becomes infinite for a=r=r', i.e. at the boundary. Again a norm which takes into account the derivatives up to a certain order will solve this problem.

In a reproducing kernel Hilbert space the problem of determining an approximation to f from n observations is easily solved if the observations are related to f in a linear manner, i.e. through linear functionals, L_i(f) = y_i. The "optimal" solution is the projection spanned on the n-dimensional sub-space spanned by the so-called representers of the linear functionals, L_i(K(P,Q)) = K(L_i,Q). The projection is the intersection between the subspace and the subspace which consists of functions which agree exactly with the observations, L_i(g)=y_i. (See Figure 1).

f {L_if=y_i}

{L_if=0}

0

{a_iK(L_i,P)}

Figure 1. The Collocation estimates the "point" with shortest distance from zero agreeing with the observations, and is simultaneously the orthogonal projection of f on the subspace spanned by the basefunctions K(L_i,P).

Then

with

If the data contain noise, then the elements _ij of the variance-covariance matrix of the noise-vector is added to K(L_i,L_j). The solution will then both minimalize the norm of f and the noise variance. If the noise is zero, the solution will agree exactly with the observations. This is the reason for the name collocation. Upper limits for the approximation error may be calculated if the norm of f is known,

This is not very useful, since the upper limit generally is very pessimistic (see Tscherning, 1986).

The collocation method may be used to solve inverse problems related to situations where the data are distributed at altitude, e.g. collected in space- or aircrafts. The problem we are left with is then which norm ? The next section should give an answer to this problem.

Example 5. For the internal density distribution of the Earth, it is difficult to construct an appropriate Hilbert space. - There are too many possible functions. If the start point are all usual polynomials xⁱy^jz^k (i.e. no discontinuities !) then we may construct possible Hilbert spaces, see Tscherning (1993). They will have base functions equal to the surface spherical harmonics (cf. example 3) multiplied by Jacobi polynomials of kind (0,2). (The Jacobi polynomials are orthogonal polynomials here of a variable equal to r/a, see also Ballani et al., 1993).

3. The statistical model.

The reproducing kernel Hilbert space has an equivalent interpretation as a stochastic process. The functionals L_i are stochastic variables, and the reproducing kernel K(P,Q) is the covariance function. The solution to the collocation problem is also the optimal statistical solution. However, in order to obtain a solution which is not optimal for all functions in a given RKHS, but for the specific function we want to model, then the covariance function (or the equivalent reproducing kernel) must be derived from this function. Here the idea is to find a function which in a least squares sense minimalizes the errors for the specific function we are dealing with.

On the other hand, statistics normally deals with repetitions, and covariances are estimated from repetitions. Here we have a problem, since there is only one Earth available. But we may regard the Earth as a new outcome of a random generator, if it is turned randomly around zero. This corresponds to the introduction of a hypothesis of stationarity for a time-series, where the covariance is supposed to be only dependent of the time difference (compare example 4). This is the basis for the concept of the so-called empirical covariance function COV(P,Q), which expresses the covariance of two values f(P) and f(Q). For a function expanded in solid spherical harmonics, the covariance function related to the function f becomes

where P_i are the Legendre-polynomials, the spherical distance between P and Q and

are the so-called degree-variance.

Note the similarity between the covariance function and the reproducing kernels in example 4. We realize, that by selecting a reproducing kernel equal to or closely approximating a covariance function, we have implicitly selected an inner product ! However, one problem in the expression for the covariance function is that it seems like we have to know the function completely (all the a_ij), and not just a finite set of observations. A practical solution to this problem will be described below.

For Fourier series a similar expression gives us the auto-covariance function. It is similar to the kernel in the L₂-norm, but each term is multiplied with the elements of the power-spectrum.

The estimation of a covariance function is done by forming mean values of product sums of observations having spherical distance within a given (small) interval. This is exactly like a covariance is estimated for a normally distributed quantity. The problem is then that we need infinitely many estimates in order to determine the infinitely many degree-variances or constituents in the power-spectrum.

The solution is to adopt a hypothesis about the decrease of the quantities as the degree i goes to infinity. An important factor is here that the degree-variances must sum to a finite number, equal to the signal variance. This means that they must go faster to zero than 1/i. This may be repaired, by supposing that they go exponentially to zero like qⁱ, where q < 1. Many different "laws" have been proposed, see Moritz (1980). One of the most popular combines a polynomial decrease in i with an exponential factor. The clever choice of a model for the degree-variances, may lead to very simple expressions for the covariance function (Knudsen, 1993, eq. 31)

where this is the covariance function between two gravity anomalies, g, at the surface of the Earth. is the density, h₀ is the depth to a density layer producing the gravity anomalies, s is the distance between P and Q, and c and D are determined empirically, so that the covariance function fit the empirically determined discrete values, see Figure 3.

Finally it is worth noticing, that the models of the degree-variances may allow us to find out which mathematically defined norm we have minimalized. Typically norms which minimalize a weighted combination of first and second derivatives (curvature like for splines) are used implicitly when the degree-variances decrease like the degree raised to a power < -1. An exponential decrease corresponds to strong smoothness conditions. For harmonic functions this means that the set of harmonicity is larger than the set on which the integration is performed when calculating the inner product. (This is what has lead to the concept of the Bjerhammar-sphere, a sphere inside the Earth, as bounding the set of harmonicity, see Moritz (1980)).

"Pure" mathematical methods may, however, also be given a statistical interpretation. The use of the indicator functions in example 1 with the L₂-norm corresponds to having used statistically independent blocks. This is quite unlikely in reality, and alternatives have therefore also been proposed, see Tscherning (1991).

If the pure mathematically defined inner products are used, the error estimate (eq.(13)) expres

ses the upper bound for the error. There is no natural way to balance the noise variance minimalization with respect to the norm of the function. The weight ratio between the square of the two norms is a factor which must be determined iteratively. Here the use of the normal equation matrix

enables the calculation of error estimates (variance of true value minus estimated)

In the diagonal of the normal equation matrix we have the sum of the signal and the noise variance. We can say that we here have obtained a natural balance between the signal and the noise.

4. Statistical models for Density distributions.

For the density distribution of the Earth, the problem of determining the empirical covariance function is nearly unsolvable. It is as difficult as the inverse problem itself. Attempts by Strykowski (1994) of estimating empirically the statistical properties of the density distribution show that data with many different probability distributions may occur in the same geological region.

We may a-priori adopt a "law", relating the exterior field and the density distribution. Examples are the so-called quasi-harmonic functions (Tscherning & Suenkel, 1981). Here a one-to-one relationship is established between the degree-variances of the exterior gravity potential and the density distribution developed in quasi-harmonic base functions. These functions are the interior solid spherical harmonics multiplied by a non-zero function of the radial distance r, only.

The use of these functions were not successful (see Hein et al. 1991), because they concentrated the density near to the earths surface. A way out of this is to use other classes of base functions, which use the Jacobi functions with r as an argument, cf. section 2.

A more straightforward situation occurs, if we suppose that the sources are topographic variations about a density interface, like base-rock or ocean depth, see Knudsen (1993). In the latter case we may even be able to estimate the auto-covariance of the depth variations, the gravity anomalies and the cross-covariance function between the two quantities. These functions may then be represented analytically by simple expressions,

where this is the auto-covariance function for height variations of a density interface and

gives the cross-covariance of the height variations and the gravity anomalies.

The inverse problem of finding topographic variations from gravity anomalies may be solved using collocation with the functions given in eq. (16),(19) and (20) equal to K(L_i,L_j) as discussed in section 2, see Tscherning et al. (1994).

5. Conclusion.

We have seen that we for the solution of inverse problems have to our disposal strong functional analytic and strong statistical methods, which to a certain extend are

equivalent. From statistical information (the covariance function) we may select an appropriate reproducing kernel. If we have selected a reproducing kernel, we may give this a statistical interpretation as a covariance function, and thereby judge whether the corresponding inner product gives us a solution with a smoothness similar to what we expect. Also the much favoured method of singular value decomposition is closely related to collocation. Its statistical interpretation (uncorrelated parameters) shows how cautious one must be in order not to introduce unrealistic geophysical constraints, see Sanso et al (1986, Appendix).

Analytic types of constraints (quasi-harmonicity, single-layer density) which lead to unique solutions of the inverse problems also have inherent problems. The solution is probably not to introduce new mathematical tools, but to take into account more data sources (like seismic data) when dealing with inverse potential field problems.

6. Exercises.

A: Use the program grcol (Knudsen, 1993) for determining depths to the ocean bottom in the Mediterranean, see Tscherning et al. (1994).

B: Use the reproducing kernel Hilbert space spanned by 3 spheres each with constant, but different density, with radius 500 m located at depth 1 km, on a line with 1 km between the centers. There are two observations of gravity disturbances at height 0, located in the middle between the spheres, having observed gravity of 10 and -8 mgal, respectively. What is the density in g/cm³ for the spheres.

C: Use the geocol10 program from the GRAVSOFT package (Tscherning et al. 1994) for quasi-harmonic inversion of a spherical harmonic expansion (gpm2).

Acknowledgement: This paper originated as lecture notes prepared for the Interdisciplinary Inversion Summer School 1, Mønsted Chalk Mines, Denmark, Aug. 14-19, 1994.

References:


Ballani, L., J.Engel and E.Grafarend: Global base functions for the
mass density in the interior of a massive body (Earth), Manuscripta
Geodaetica, Vol. 18, pp. 99-114, 1993.


Hein,G., F.Sanso', G.Strykowsky and C.C.Tscherning: On the Choice
of Norm and Base Functions for the Solution of the Inverse Gravimetric
Problem. Ricerche di Geodesia Topografia Fotogrammetria, Vol. 5,
pp. 121-138, CLUP, Milano, 1989.

Knudsen, P.: Integrated inversion of gravity data. Final Report
Norsk Hydro R&D Project, Kort & Matrikelstyrelsen, Geodetic Division
Technical report No. 7, Koebenhavn, 1993.

Moritz, H.: Advanced Physical Geodesy. H.Wichmann Verlag, Karlsruhe,
1980.

Sanso',F., R.Barzaghi and C.C.Tscherning: Choice of Norm for the Density
Distribution of the Earth. Geoph.Jour.Royal Astr. Soc., Vol. 87, pp.
123-141, 1986.

Strykowski, G.: Statistical analysis of the log data. Proc.
Interdisciplinary Inversion Workshop 3, Odense 1994, Ed. P. Sibani,
pp. 63-84, Odense, 1994.


Tscherning, C.C.: Functional Methods for Gravity Field Approximation.
Lecture Notes in Earth Sciences, Vol. 7, pp. 3-47, Mathematical and
Numerical Techniques in Physical Geodesy, Springer-Verlag, 1986.

Tscherning, C.C.: Density-Gravity Covariance Functions Produced
by Overlapping Rectangular Blocks of Constant Density. Geophysical
Journal International, Vol. 105, pp. 771-776, 1991.

Tscherning, C.C.: Isotropic reproducing kernels for the inner of a
sphere and their use as density covariance functions. Proc.
Interdisciplinary Inversion Workshop 2, Copenhagen, May 19, 1993,
pp. 71-76, University of Copenhagen, 1993.

Tscherning, C.C., P.Knudsen and R.Forsberg: Description of the
GRAVSOFT package. Geophysical Institute, University of Copenhagen,
Technical Report, 1991, 2. Ed. 1992, 3. Ed. 1993, 4. ed, 1994.

Tscherning, C.C., R.Forsberg & P.Knudsen: First experiments with
improvement of depth information using gravity anomalies in the
Mediterranean Sea. GEOMED Report no. 4, Ed: Arabelos & Tziavos,
pp. 133-148, Thesssaloniki, 1994.

Tscherning, C.C. and H. Suenkel: A Method for the Construction of
Spheroidal Mass Distributions consistent with the harmonic Part of the
Earth's Gravity Potential. Manuscripta Geodaetica, Vol. 6, pp.
131-156, 1981.
Last change 1997-12-01 by cct.