# Encyclopedia of Geodesy

Living Edition
| Editors: Erik Grafarend

# Disturbing Potential from Gravity Anomalies: From Globally Reflected Stokes Boundary Value Problem to Locally Oriented Multiscale Modeling

• Matthias Augustin
• Christian Blick
• Sarah Eberle
• Willi Freeden
Living reference work entry
DOI: https://doi.org/10.1007/978-3-319-02370-0_124-1

## Keywords

Gravity Anomaly Local Support Cubature Formula Gravity Disturbance Spherical Approximation
These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

## Definition

Stokes problem, Fourier series expansion in terms of outer harmonics, classical global solution by convolving gravity anomalies against the Stokes kernel , regularization of the Stokes kernel, local multiscale approximation.

## Introduction

The traditional approach of physical geodesy (cf., e.g., Heiskanen and Moritz, 1967; Moritz, 2015) starts from the assumption that scalar gravity intensity is available over the whole Earth’s surface. The gravitational part of the gravity potential can then be regarded as a harmonic function outside the Earth’s surface. A classical approach to gravity field modeling was conceived by G.G. Stokes (1849). He proposed reducing the given gravity accelerations from the Earth’s surface to the geoid (see, e.g., Listing, 1878), where the geoid is a level surface, e.g., its potential value is constant. The difference between the reduced gravity disturbing potential, i.e., the difference between the actual and the reference potential, can be obtained from a (third) boundary value problem of potential theory. M.S. Molodensky (Molodensky et al., 1960) proposed to improve Stokes’ solution by “reducing” the gravity anomalies, given on the Earth’s surface, to a “normal level surface” (telluroid). In both cases, the calculation via the associated integral formulas is usually performed in spherical approximation although concepts of ellipsoidal realization are available (see Grafarend et al., 2015 and the references therein). In fact, H. Moritz (Hofmann-Wellenhof and Moritz, 2006) mentioned that the reference surface is never a sphere in any geometrical sense but always an ellipsoid. As the flattening of the Earth is very small, the ellipsoidal formulas can be expanded into power series in terms of the flattening so that terms containing higher orders can be neglected. “In this way one obtains formulas that are rigorously valid for the sphere, but approximately valid for the actual reference ellipsoid as well” (see Hofmann-Wellenhof and Moritz, 2006).

For practical evaluation, the Stokes convolution integral between Stokes kernel and gravity anomalies must be replaced by approximate cubature formulas using certain integration weights and knots. The approximate integration formulas are the essential problem in the framework of globally determining the disturbing potential and, subsequently, the geoidal height following Bruns’ concept (see Bruns, 1878). In fact, we are confronted with the following dilemma: On the one hand, Weyl’s law of equidistribution (cf. Weyl, 1916, Cui and Freeden, 1997) tells us that numerical integration and equidistribution of the nodal points are mathematically equivalent. This law holds true for any reference surface, i.e., telluroid, ellipsoid, as well as sphere. In order to get better and better accuracy in approximate integration procedures, we thus need dense, globally over the whole reference surface equidistributed datasets. On the other hand, even nowadays, observations in sufficient data width and quality are only available for certain parts of the Earth’s surface, and there are large areas, particularly at sea, where no suitable data are given at all. In fact, terrestrial gravity data coverage now and in the foreseeable future is far from being satisfactory and totally inadequate for the purpose of high-precision geoid determination. As a consequence, Stokes’ type integral formula and its improvements based on Molodensky’s idea cannot be applied on a global basis neither in an ellipsoidal nor in a spherical framework. We have to observe the specific heterogeneous data situation. A mathematical way out is an adequate multiscale method providing a “zooming in” approximation in adaptation to the data distribution and density.

In this contribution our particular goal is a local high-resolution gravitational model reflecting the available data obligations as far as possible. Since the flattening in a local approach is negligibly small, a calculation in spherical approximation is canonical. For simplicity, we restrict ourselves to error-free data. A multiscale signal-to-noise ratio method handling noisy data is proposed, e.g., in (Freeden and Maier, 2002).

Our considerations are based on the work (Freeden et al., 1998; Freeden and Schreiner, 2006; Freeden and Wolf, 2009; Freeden and Gerhards, 2012; Freeden, 2015). The illustrations are essentially taken from (Freeden and Wolf, 2009) and the PhD thesis (Wolf, 2009).

## Stokes Wavelets

We begin our work with the recapitulation of the global Stokes’ approach in spherical approximation. Let ΩR be the sphere with radius R around the origin and the gravity anomaly Δg ∈ C(0) R ) with
$${\displaystyle \underset{\Omega_R}{\int}\Delta g\left(\boldsymbol{x}\right) ds\left(\boldsymbol{x}\right)}=0$$
(1)
and
$$\begin{array}{ll}{\displaystyle \underset{\Omega_R}{\int}\Delta g\left(\boldsymbol{x}\right)\left({\boldsymbol{\upvarepsilon}}^k\cdot \boldsymbol{x}\right)} ds\left(\boldsymbol{x}\right)=0,\hfill & k=1,2,3\hfill \end{array}$$
(2)
be given. Here, ds is the surface element and ε (1), ε (2), ε (3) are the canonical cartesian unit vectors in ℝ3. Then, the disturbing potential $$T:\overline{\Omega_R^{\mathrm{ext}}}\to \mathrm{\mathbb{R}}$$ is the unique solution of the exterior Stokes boundary-value problem (see also Freeden, 1978; Freeden and Wolf, 2009; Wolf, 2009):
1. (i)

T is continuously differentiable in $$\overline{\Omega_R^{\mathrm{ext}}}$$ and twice continuously differentiable in $${\Omega}_R^{\mathrm{ext}}$$ i.e., $$T\in {\mathrm{C}}^{(1)}\left(\overline{\Omega_R^{\mathrm{ext}}}\right)\cap {\mathrm{C}}^{(2)}\left({\Omega}_R^{\mathrm{ext}}\right),$$

2. (ii)

T is harmonic in Ω R ext , i.e., Δ x T = 0 in Ω R ext

3. (iii)

T is regular at infinity,

4. (iv)
$${\displaystyle {\int}_{\Omega_R}T\left(\boldsymbol{y}\right){H}_{-n-1,k}^R\left(\boldsymbol{y}\right)} ds\left(\boldsymbol{y}\right)=0,n=0,1, k=1,\dots, 2n+1,$$

5. (v)
$$\begin{array}{ll}-\frac{\boldsymbol{x}}{\left|\boldsymbol{x}\right|}\cdot {\nabla}_{\boldsymbol{x}}T\left(\boldsymbol{x}\right)-\frac{2}{\left|\boldsymbol{x}\right|}T\left(\boldsymbol{x}\right)-\frac{2}{\left|\boldsymbol{x}\right|}T\left(\boldsymbol{x}\right)=\Delta g\left(\boldsymbol{x}\right),\hfill & \boldsymbol{x}\in {\Omega}_R\hfill \end{array}.$$

$${\Omega}_R^{\mathrm{ext}}$$ denotes the exterior of the sphere $${\Omega}_R$$.

T is determined by Stokes integral formula
$$T\left(R\boldsymbol{\xi} \right)=\frac{1}{4\pi R}{\displaystyle \underset{\Omega_R}{\int } St\left(R\boldsymbol{\xi}, R\boldsymbol{\eta} \right)\Delta g\left(R\boldsymbol{\eta} \right)ds\left(R\boldsymbol{\eta} \right)}=\frac{R}{4\pi }{\displaystyle \underset{\Omega}{\int } St\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)\Delta g\left(R\boldsymbol{\eta} \right)\;ds\left(\boldsymbol{\eta} \right)},$$
(3)
with
$$St\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)=1-5\boldsymbol{\xi} \cdot \boldsymbol{\eta} -6{\left(S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^{-1}+S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}+\frac{1}{{\left(S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}\right)$$
(4)
is the Stokes kernel , ξ, η ∈ Ω = Ω1, where we have used the abbreviation
$$\begin{array}{ll}S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)=\frac{\sqrt{2}}{\sqrt{1-\boldsymbol{\xi} \cdot \boldsymbol{\eta}}},\hfill & 1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ne 0\hfill \end{array}.$$
(5)
To regularize the improper integral Eq. 3, we replace the zonal kernel S by the space-regularized zonal kernel (see, e.g., Freeden and Schreiner, 2006; Freeden and Wolf, 2009; Wolf, 2009)
$${S}^{\rho }(t)=\left\{\begin{array}{ll}\frac{R}{\rho}\left(3-\frac{2{R}^2}{\rho^2}\left(1-t\right)\right),\hfill & 0<1-t\le \frac{\rho^2}{2{R}^2},\hfill \\ {}\frac{\sqrt{2}}{\sqrt{1-t}},\hfill & \frac{\rho^2}{2{R}^2}<1-t\le 2.\hfill \end{array}\right.$$
(6)
Clearly, the function (depicted in Figure 1) S ρ is continuously differentiable on the interval [−1, 1], and we have (see Freeden and Wolf, 2009; Wolf, 2009)
$$\left({S}^{\rho}\right)^{\prime }(t)=\left\{\begin{array}{ll}\frac{2{R}^3}{\rho^3},\hfill & 0\le 1-t\le \frac{\rho^2}{2{R}^2},\hfill \\ {}\frac{1}{\sqrt{2}{\left(1-t\right)}^{\frac{3}{2}}},\hfill & \frac{\rho^2}{2{R}^2}<1-t\le 2.\hfill \end{array}\right.$$
(7)
Furthermore, the functions S and S ρ are monotonically increasing on the interval [−1, 1), such that S(t) ≥ S ρ (t) ≥ S(−1) = S ρ (−1) = 1 holds true on the interval [−1, 1). Considering the difference between the kernel S and its linearly regularized version S ρ , we find
$$S(t)-{S}^{\rho }(t)=\left\{\begin{array}{ll}\frac{\sqrt{2}}{\sqrt{1-t}}-\frac{R}{\rho}\left(3-\frac{2{R}^2}{\rho^2}\left(1-t\right)\right),\hfill & 0<1-t\le \frac{\rho^2}{2{R}^2},\hfill \\ {}0,\hfill & \frac{\rho^2}{2{R}^2}<1-t\le 2.\hfill \end{array}\right.$$
(8)
It can be shown (Freeden and Schreiner, 2006) that the following lemma holds:

### Lemma 1

For F ∈ C(0)(Ω) and S ρ defined by Eq. 6, we have
$$\underset{\rho \to 0+}{ \lim}\underset{\boldsymbol{\xi} \in \varOmega }{ \sup}\left|{\displaystyle \underset{\varOmega }{\int }S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)F\left(\boldsymbol{\eta} \right) ds\left(\boldsymbol{\eta} \right)}-{\displaystyle \underset{\varOmega }{\int }{S}^{\rho}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)F\left(\boldsymbol{\eta} \right) ds\left(\boldsymbol{\eta} \right)}\right|=0.$$
(9)
To obtain another useful convergence result, we observe that for all t ∈ [−1, 1) with $$1-t\le \frac{\rho^2}{2{R}^2}$$
$$\begin{array}{ll} \ln \left(\frac{1}{S(t)}+\frac{1}{{\left(S(t)\right)}^2}\right)- \ln \left(\frac{1}{S^{\rho }(t)}+\frac{1}{{\left({S}^{\rho }(t)\right)}^2}\right)=& \ln \left(1+S(t)\right)\hfill \\ {}& - \ln \left(1+{S}^{\rho }(t)\right)\hfill \\ {}& -2\Big( \ln \left(S(t)- \ln \left({S}^{\rho }(t)\right)\right)\hfill \end{array}$$
(10)
and, thus,
$$\left| \ln \left(\frac{1}{S(t)}+\frac{1}{{\left(S(t)\right)}^2}\right)- \ln \left(\frac{1}{S^{\rho }(t)}+\frac{1}{{\left({S}^{\rho }(t)\right)}^2}\right)\right|=O\left(\left|S(t)-{S}^{\rho }(t)\right|\right).$$
(11)
This leads to the following result:

### Lemma 2

Let S be the singular kernel given by $$S(t)=\frac{\sqrt{2}}{\sqrt{1-t}}$$ and let S ρ , ρ ∈ (0, 2R], be the corresponding (Taylor) linearized regularized kernel defined by Eq. 6. Then
$$\underset{\rho \to 0+}{ \lim }{\displaystyle \underset{-1}{\overset{1}{\int }}\left| \ln \left(1+S(t)\right)- \ln \left(1+{S}^{\rho }(t)\right)\right|dt=0,}$$
(12)
$$\underset{\rho \to 0+}{ \lim }{\displaystyle \underset{-1}{\overset{1}{\int }}\left| \ln \left(\frac{1}{S(t)}+\frac{1}{{\left(S(t)\right)}^2}\right)- \ln \left(\frac{1}{S^{\rho }(t)}+\frac{1}{{\left({S}^{\rho }(t)\right)}^2}\right)\right|dt=0,}$$
(13)
$$\underset{\rho \to 0+}{ \lim }{\displaystyle \underset{-1}{\overset{1}{\int }}\left({\left(S(t)\right)}^2-{\left({S}^{\rho }(t)\right)}^2\right)}\sqrt{1-{t}^2}dt=0.$$
(14)
The regularization given in Eq. 6 leads us to the following regularized global representation of the disturbing potential corresponding to gravity anomalies as boundary data (see Freeden and Wolf, 2009):
$${T}^{\rho}\left(R\boldsymbol{\xi} \right)=\frac{R}{4\pi }{\displaystyle \underset{\Omega}{\int }S{t}^{\rho}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)\Delta g\left(R\boldsymbol{\eta} \right) ds\left(\boldsymbol{\eta} \right)}$$
(15)
$$\begin{array}{l}S{t}^{\rho}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)= 1-5\boldsymbol{\xi} \cdot \boldsymbol{\eta} -6{\left(S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^{-1}\hfill \\ {}+ {S}^{\rho}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{{\left({S}^{\rho}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}+\frac{1}{S^{\rho}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}\right)\hfill \\ {}=1-5\boldsymbol{\xi} \cdot \boldsymbol{\eta} -6\frac{\sqrt{1-\boldsymbol{\xi} \cdot \boldsymbol{\eta}}}{\sqrt{2}}\hfill \\ {}\begin{array}{l}\hfill \\ {}+\left\{\begin{array}{ll}\frac{R}{\rho}\left(3-\frac{2{R}^2}{\rho^2}\left(1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)\hfill & \hfill \\ {}-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(1+\frac{R}{\rho}\left(3-\frac{2{R}^2}{\rho^2}\left(1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)\right)\hfill & \hfill \\ {}+6\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{R}{\rho}\left(3-\frac{2{R}^2}{\rho^2}\left(1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)\right),\hfill & 0\le 1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \le \frac{\rho^2}{2{R}^2},\hfill \\ {}\frac{\sqrt{2}}{\sqrt{1-\boldsymbol{\xi} \cdot \boldsymbol{\eta}}}-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1-\boldsymbol{\xi} \cdot \boldsymbol{\eta}}{2}+\frac{\sqrt{1-\boldsymbol{\xi} \cdot \boldsymbol{\eta}}}{\sqrt{2}}\right),\hfill & \frac{\rho^2}{2{R}^2}<1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \le 2,\hfill \end{array}\right.\hfill \end{array}\end{array}$$
(16)
for ξ, η ∈ Ω and ρ ∈ (0, 2R]. Here, we have made use of Eq. 10.

With Lemma 1 and Lemma 2, we obtain

### Theorem 3

Suppose that T is the solution of the Stokes boundary-value problem. Let T ρ , ρ ∈ (0, 2R], represent its regularization as in Eq. 15 . Then
$$\underset{\rho \to 0+}{ \lim}\underset{\boldsymbol{\xi} \in \Omega}{ \sup}\left|T\left(R\boldsymbol{\xi} \right)-{T}^{\rho}\left(R\boldsymbol{\xi} \right)\right|=0.$$
(17)
The linear space regularization technique enables us to formulate multiscale solutions for the disturbing potential from gravity anomalies. For numerical application, we have to go over to scale-discretized approximations of the solution to the boundary-value problem. For that purpose, we choose a monotonically decreasing sequence $${\left\{{\rho}_j\right\}}_{j\in {\mathrm{\mathbb{N}}}_0}$$, such that
$$\begin{array}{ll}\underset{j\to \infty }{ \lim }{\rho}_j=0,\hfill & {\rho}_0=2R.\hfill \end{array}$$
(18)
A particularly important example, which we use in our numerical implementations below, is the dyadic sequence with
$$\begin{array}{lll}{\rho}_j={2}^{1-j}R,\hfill & j\in \mathrm{\mathbb{N}},\hfill & {\rho}_0=2R.\hfill \end{array}$$
(19)
It is easily seen that $$2{\rho}_{j+1}={\rho}_j$$,$$j\in {\mathrm{\mathbb{N}}}_0$$, is the relation between two consecutive elements of the sequence. In correspondence to the sequence $${\left\{{\rho}_j\right\}}_{j\in {\mathrm{\mathbb{N}}}_0,}$$ a sequence $${\left\{S{t}^{\rho j}\right\}}_{j\in {\mathrm{\mathbb{N}}}_0}$$ of discrete versions of the regularized Stokes kernels Eq. 16, so-called Stokes scaling functions, is available. Figure 2 shows a graphical illustration of the regularized Stokes kernels for different scales j.
The regularized Stokes wavelets, forming the sequence $${\left\{WS{t}^{\rho_j}\right\}}_{j\in {\mathrm{\mathbb{N}}}_0}$$, are understood to be the difference of two consecutive regularized Stokes scaling functions, respectively,
$$\begin{array}{ll}WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)=S{t}^{\rho_{j+1}}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)-S{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right),\hfill & j\in {\mathrm{\mathbb{N}}}_0.\hfill \end{array}$$
(20)
These wavelets possess the numerically nice property of a local support. More specifically, the function $$\boldsymbol{\eta} \mapsto WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)$$, η ∈ Ω, vanishes everywhere outside the spherical cap $${\Gamma}_{\rho_j^2/2{R}^2}\left(\boldsymbol{\xi} \right)$$.
Explicitly written out, we have
$$\begin{array}{l}WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)=\hfill \\ {}\left\{\begin{array}{ll}{S}^{\rho_{j+1}}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{S^{\rho_{j+1}}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}+\frac{1}{{\left({S}^{\rho_{j+1}}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}\right)\hfill & \hfill \\ {}-{S}^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)+3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{S^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}+\frac{1}{{\left({S}^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}\right),\hfill & 0\le 1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \le \frac{\rho_{j+1}^2}{2{R}^2},\hfill \\ {}S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)-3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}+\frac{1}{{\left(S\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}\right)\hfill & \hfill \\ {}-{S}^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)+3\boldsymbol{\xi} \cdot \boldsymbol{\eta} \ln \left(\frac{1}{S^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)}+\frac{1}{{\left({S}^{\rho_j}\left(\boldsymbol{\xi} \cdot \boldsymbol{\eta} \right)\right)}^2}\right),\hfill & \frac{\rho_{j+1}^2}{2{R}^2}<1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \le \frac{\rho_j^2}{2{R}^2},\hfill \\ {}0,\hfill & \frac{\rho_j^2}{2{R}^2}<1-\boldsymbol{\xi} \cdot \boldsymbol{\eta} \le 2.\hfill \end{array}\right.\end{array}$$
(21)
Let J ∈ ℕ0 be an arbitrary scale. Suppose that $$S{t}^{\rho_J}$$ is the regularized Stokes scaling function at scale J. Furthermore, let $$WS{t}^{\rho_j}$$, j = 0,…,J, be the regularized Stokes wavelets as given by Eq. 21. Then an easy manipulation shows that
$$S{t}^{\rho_J}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)=S{t}^{\rho_{J_0}}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)+{\displaystyle \sum_{j={J}_0}^{J-1}WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)}.$$
(22)
The local support of the Stokes wavelets within the framework of Eq. 22 should be studied in more detail: Following the sequence given by Eq. 19, we start with a globally supported scaling kernel $$S{t}^{\rho_0}=S{t}^{2R}$$. Then we add more and more wavelet kernels $$WS{t}^{\rho_j}$$, j = 0,…,J, to achieve the required scaling kernel $$S{t}^{\rho_J}$$. It is of particular importance that $$\boldsymbol{\eta} \mapsto WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right)$$, ξ ∈ Ω fixed, are ξ-zonal functions and possess spherical caps as local supports. Clearly, the support of the wavelets $$WS{t}^{\rho_j}$$ becomes more and more localized for increasing scales j. In conclusion, a calculation of an integral representation for the disturbing potential T starts with a global trend approximation using the scaling kernel at scale j = 0 (of course, this requires data on the whole sphere, but the data can be rather sparsely distributed since they only serve as a trend approximation). Step by step, we are able to refine this approximation by use of wavelets. The increasing spatial localization of the wavelets successively allows a better spatial resolution of the disturbing potential T. Additionally, the local supports of the wavelets have a computational advantage since the integration is reduced from the entire sphere to smaller and smaller spherical caps. Consequently, the presented numerical technique becomes capable of handling heterogeneously distributed data sets in adaptation to their mutual spacing.
All in all, keeping the space-localizing properties of the regularized Stokes scaling and wavelet functions in mind, we are able to establish an approximation of the solution of the disturbing potential T from gravity anomalies Δg in the form of a zooming-in multiscale method . A low-pass filtered version of the disturbing potential T at the scale j in an integral representation over the unit sphere Ω is given by
$$\begin{array}{ll}{T}^{\rho_j}\left(R\boldsymbol{\xi} \right)=\frac{R}{4\pi }{\displaystyle \underset{\Omega}{\int}\Delta g\left(R\boldsymbol{\eta} \right) S{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right) ds\left(\boldsymbol{\eta} \right)},\hfill & \boldsymbol{\xi} \in \Omega, \hfill \end{array}$$
(23)
while the j-scale band-pass filtered version of T leads to the integral representation
$$\begin{array}{ll}W{T}^{\rho_j}\left(R\boldsymbol{\xi} \right)=\frac{R}{4\pi }{\displaystyle \underset{\Gamma_{\rho_j^2/2{R}^2}\left(\boldsymbol{\xi} \right)}{\int}\Delta g\left(R\boldsymbol{\eta} \right) WS{t}^{\rho_j}\left(\boldsymbol{\xi}, \boldsymbol{\eta} \right) ds\left(\boldsymbol{\eta} \right)},\hfill & \boldsymbol{\xi} \in \Omega .\hfill \end{array}$$
(24)

### Theorem 4

Let $${T}^{\rho_{J_0}}$$ be the regularized version of the disturbing potential at some arbitrary initial scale J 0 as given in Eq. 23 , and let $$W{T}^{\rho_j}$$, j = J 0 , J 0 + 1,…, be given by Eq. 24 . Then, the following reconstruction formula holds true:
$$\underset{J\to \infty }{ \lim}\underset{\boldsymbol{\xi} \in \Omega}{ \sup}\left|T\left(R\boldsymbol{\xi} \right)-\left({T}^{\rho_{J_0}}\left(R\boldsymbol{\xi} \right)+{\displaystyle \sum_{j={J}_0}^{J-1}W{T}^{\rho_j}\left(R\boldsymbol{\xi} \right)}\right)\right|=0.$$
The multiscale procedure (wavelet reconstruction) as developed here can be illustrated by the following scheme:
$$\begin{array}{ll}W{T}^{\rho_{J_0}}\hfill & W{T}^{\rho {J}_0+1}\hfill \\ {} \searrow \hfill & \searrow \hfill \\ {}{T}^{\rho {J}_0}\to +\to \hfill & {T}^{\rho {J}_0+1}\to +\to {T}^{\rho {J}_0+2}\dots .\hfill \end{array}$$
(25)
Consequently, a tree algorithm based on regularization in the space domain has been realized for determining the disturbing potential T from locally available data sets of gravity anomalies Δg. An example is shown in Figure 3. The fully discretized multiscale approximations have the following representations
$$\begin{array}{ll}{T}^{\rho_J}\left(R\boldsymbol{\xi} \right)\simeq \frac{R}{4\pi }{\displaystyle \sum_{k=1}^{N_J}{w}_k^{N_J}\Delta g\left(R{\boldsymbol{\eta}}_k^{N_J}\right)} S{t}^{\rho_J}\left(\boldsymbol{\xi}, {\boldsymbol{\eta}}_k^{N_J}\right),\hfill & \boldsymbol{\xi} \in \Omega \hfill \end{array},$$
(26)
$$\begin{array}{ll}W{T}^{\rho_j}\left(R\boldsymbol{\xi} \right)\simeq \frac{R}{4\pi }{\displaystyle \sum_{k=1}^{N_j}{w}_k^{N_j}\Delta g\left(R{\boldsymbol{\eta}}_k^{N_j}\right)} WS{t}^{\rho_j}\left(\boldsymbol{\xi}, {\boldsymbol{\eta}}_k^{N_j}\right),\hfill & \boldsymbol{\xi} \in \Omega \hfill \end{array},$$
(27)
where $${\boldsymbol{\eta}}_k^{N_j}$$ are the integration knots and $${w}_k^{N_j}$$ the integration weights. Whereas the sum in Eq. 26 has to be extended over the whole sphere Ω, the summation in Eq. 27 has to be computed only for the local supports of the wavelets (note that the symbol $$\simeq$$ means that the error between the right-and the left-hand side can be neglected).
Figures 4, 5, and 6 show that the method presented here solves a dilemma in geodesy: Common global solution methods need ever denser, globally equidistributed data sets over the whole sphere Ω R to obtain a better approximation quality (according to Weyl’s Law of Equidistribution). However, the reality is quite different. On the one hand, we have large gaps in data sets, particularly at sea. On the other hand, there are some regions where the accuracy and density of available data sets is quite remarkable. The solution offered by our wavelet method is to start with a coarse, global approximation, e.g., of the disturbing potential using a scaling function of scale J 0 and add local refinement in the form of band-pass filtered versions using Stokes’ wavelets. This can be realized as these wavelets only have a compact support. This procedure allows the incorporation of heterogeneous data sets in a way that locally improves the approximation of the disturbing potential despite nonequidistributed data sets.

## Conclusion

So far, as pointed out in Hofmann-Wellenhof and Moritz (2006), much more gravity anomalies than gravity disturbances are available and are being processed. In the future, we may expect a change in the practice of physical geodesy because of GPS. In this respect, the approach as presented here for the Stokes problem can be formulated for the Neumann problem as well (see, e.g., Wolf, 2009; Freeden and Gerhards, 2012), i.e., we are able to make the transition from the globally reflected Neumann problem to locally oriented multiscale modeling in an analogous way.

## References and Reading

1. Bruns, E. H., 1878. Die Figur der Erde. Publikationen des Königlichen Preussischen Geod ä tischen Instituts. Berlin, P. Stankiewicz Buchdruckerei.Google Scholar
2. Cui, J. and Freeden, W., 1997. Equidistribution on the sphere. SIAM 18, 595–609.Google Scholar
3. Freeden, W., 1978. An application of a summation formula to numerical computation of integrals over the sphere. Bulletin Géodésique, 52, 165–175.
4. Freeden, W., 2015. Geomathematics: Its role, its aim, and its potential. In Freeden, W., Nashed, Z., and Sonar, T. (eds.), Handbook of Geomathematics. 2nd edn pp. 3–78. Heidelberg: Springer.Google Scholar
5. Freeden, W., and Gerhards, C., 2012. Geomathematically Oriented Potential Theory. Boca Raton: Chapman & Hall/CRC.
6. Freeden, W., and Maier, T., 2002. On multiscale denoising of spherical functions: basic theory and numerical aspects. Electronic Transactions on Numerical Analysis (ETNA), 14, 56–78.Google Scholar
7. Freeden, W., and Schreiner, M., 2006. Local multiscale modelling of geoid undulations from deflections of the vertical. Journal of Geodesy, 79, 641–651.
8. Freeden, W., and Wolf, K., 2009. Klassische Erdschwerefeldbestimmung aus der Sicht moderner Geomathematik. Mathematische Semesterberichte, 56, 53–77.
9. Freeden, W., Gervens, T., and Schreiner, M., 1998. Constructive Approximation on the Sphere (With Applications to Geomathematics). Oxford: Oxford Science Publications/Clarendon Press.Google Scholar
10. Grafarend, E. W., Klapp, M., and Martinec, Z., 2015. Spacetime modelling of the Earth’s gravity field by ellipsoidal harmonics. In Freeden, W., Nashed, Z., and Sonar, T. (eds.), Handbook of Geomathematics, 2nd edn, pp. 381–496, Heidelberg: Springer.Google Scholar
11. Heiskanen, W. A., and Moritz, H., 1967. Physical Geodesy. San Francisco: W.H. Freeman.Google Scholar
12. Hofmann-Wellenhof, B., and Moritz, H., 2006. Physical Geodesy, 2nd edn. Wien/New York: Springer.Google Scholar
13. Listing, J. B., 1878. Neue geometrische und dynamische Constanten des Erdkörpers. Nachrichten von der Königlichen Gesellschaft der Wissenschaften und der Georg-Augusts-Universität zu Göttingen, pp. 749–815.Google Scholar
14. Molodensky, M. S., Eremeev, V. F., and Yurkina, M. I., 1960. Methods for study of the external gravitational field and figure of the Earth. Trudy TsNIIGAiK, Geodezizdat, Moscow, p. 131 (English translat.: Israel Program for Scientific Translation, Jerusalem, 1962).Google Scholar
15. Moritz, H., 2015. Classical physical geodesy. In Freeden, W., Nashed, Z., and Sonar, T. (eds.), Handbook of Geomathematics. 2nd edn, pp. 253–290, Heidelberg: Springer.Google Scholar
16. Stokes, G. G., 1849. On the variation of gravity on the surface of the Earth. Transactions of the Cambridge Philosophical Society, 8, 672–695.Google Scholar
17. Weyl, H., 1916. Über die Gleichverteilung von Zahlen mod Eins. Mathematische Annalen, 77, 313–352.
18. Wolf, K., 2009. Multiscale modeling of classical boundary value problems in physical geodesy by locally supported wavelets. PhD thesis, University of Kaiserslautern, Geomathematics Group.Google Scholar