# On multiscale denoising of spherical functions: basic theory and numerical aspects.

Abstract. The basic concepts of selective multiscale reconstruction
of functions on the sphere from error-affected data is outlined for
scalar functions. The selective reconstruction mechanism is based on the
premise that multiscale approximation can be well-represented in terms
of only a relatively small number of expansion coefficients at various
resolution levels. A new proof, including non-bandlimited kernel
functions, of the pyramid scheme is presented to efficiently remove the
noise at different scales using a priori statistical information, i.e.
knowledge of the covariance function.

Key words. spherical wavelet theory, scalar multiscale approximation, pyramid scheme, spectral and multiscale variance-covariance model, hard and soft thresholding.

AMS subject classifications. 33C55, 42C40, 62-07, 65T60, 86A25.

1. Introduction. While standard Fourier methods in terms of spherical harmonics are very successful at picking out frequencies from a spherical signal, they are utterly incapable of dealing properly with data changing on small spatial scales. This fact has been well-known for years. To improve the applicability of Fourier analysis, various methods such as 'windowed Fourier transform' have been developed on the sphere to modify the usual Fourier procedure to allow analysis of the frequency content of a signal at each position (cf. [6, 8]). However, the amount of localization in space and in frequency is not completely satisfactory in the Fourier transform and its windowed variant. For example, geopotentials refer to a certain combination of frequencies, and the frequencies themselves are spatially changing. This space evolution of the frequencies is not reflected in the Fourier transform in terms of non-space-localizing spherical harmonics. Even the windowed Fourier transform contains information about frequencies over a certain area of positions instead of showing how the frequencies vary in space. With spherical wavelets, the amount of localization in space and in frequency is automatically adapted, in that only a narrow space-window is needed to examine high-frequency content, but a wide space-window is allowed when investigating low-frequency phenomena. The basic framework of this approach has been provided by the spherical wavelet theory developed by the Geomathematics Group at the University of Kaiserslautern during the last years (see http://www.mathematik.uni-kl.de/~wwwgeo/pub1.html)

When dealing with real geophysically relevant data it should be kept in mind that each measurement does not really give the value of the observable under consideration but that--at least to some extent--the data are contaminated with noise. That is, in order to successfully improve geomathematical field modeling, one main aspect is to extract the true portion of the observable from the actual signal. In consequence, a particular emphasis lies on the subject of denoising. This endeavor is precisely the goal in statistical function estimation. Here, the interest is to 'smooth' the noisy data in order to obtain an estimate of the underlying function. The Euclidean theory of wavelets provides signal processors with new, fast tools that are well-suited for denoising signals (for a survey the reader is e.g. referred to [13] and the references therein).

The objective of this article is to introduce multiscale signal-to-noise thresholding and to provide the wavelet oriented basis of denoising spherical data sets. First, we develop the corresponding theory of denoising spherical functions by non-rotation invariant wavelets (cf. [9] for the rotation-invariant case). With the basic introduction at hand, selective thresholding within a pyramid scheme is presented. The thresholding scheme is designed to distinguish between coefficients which contribute significantly to the signal, and those which are negligible. It should be noted that our approach is essentially influenced by the concept of sparse wavelet representations in Euclidean spaces (cf. [15, 2, 3]) and the stochastic model used in satellite geodesy (see e.g. [14]). Using a multiscale approach we are thus able to include detail information of small spatial extent while suppressing the noise in the approximation appropriately. A simple example of denoising geomagnetic field data will be given as an illustration.

2. Preliminaries. Let [R.sup.3] denote three-dimensional Euclidean space. For x, y [member of] [R.sup.3], x = [([x.sub.1], [x.sub.2], [x.sub.3]).sup.T], y = [([y.sub.1], [y.sub.2], [y.sub.3]).sup.T], the inner product is defined as usual by

x * y = [x.sup.T] y = [3.summation over (i=1)][x.sub.i][y.sub.i].

For all elements x [member of] [R.sup.3], x = [([x.sub.1], [x.sub.2], [x.sub.3]).sup.T], different from the origin, we have

x = r[xi], r = |x| = [square root of (x * x)] = [x.sup.2.sub.1] + [x.sup.2.sub.2] + [x.sup.2.sub.3],

where [xi] = [([[xi].sub.1], [[xi].sub.2], [[xi].sub.3]).sup.T] is the uniquely determined directional unit vector of x. The unit sphere in [R.sup.3] is denoted by [OMEGA]. If the vectors [[epsilon].sup.1], [[epsilon].sup.2], [[epsilon].sup.3] form the canonical orthonormal basis in [R.sup.3], the points [xi] [member of] [OMEGA] may be represented in polar coordinates by

[xi] = t[[epsilon].sup.3] + [square root of (1 - [t.sup.2])] (cos [psi][[epsilon].sup.1] + sin [psi][[epsilon].sup.2]), t = cos v, v [member of] [0, [pi]], [psi] [member of] [0, 2[pi]).

3. Spherical Harmonics. The spherical harmonics [Y.sub.n] of degree n are defined as the everywhere on [OMEGA] infinitely differentiable eigenfunctions of the Beltrami operator [[DELTA].sup.*] corresponding to the eigenvalues ([[DELTA].sup.*])^(n) = -n(n + 1), n = 0, 1, ..., where the Beltrami-operator is the angular part of the Laplace-operator [DELTA] in [R.sup.3]. As it is well-known, the functions [H.sub.n] : [R.sup.3] [right arrow] R defined by [H.sub.n](x) = [r.sup.n][Y.sub.n]([xi]), x = r[xi], [xi] [member of] [OMEGA], are homogeneous polynomials in rectangular coordinates which satisfy the Laplace-equation [[DELTA].sub.x][H.sub.n](x) = 0, x [member of] [R.sup.3]. Conversely, every homogeneous harmonic polynomial of degree n when restricted to [OMEGA] is a spherical harmonic of degree n. The Legendre polynomials [P.sub.n] : [-1, +1] [right arrow] [-1, +1] are the only everywhere on [-1, +1] infinitely differentiable eigenfunctions of the Legendre-operator (1 - [t.sup.2])[(d/dt).sup.2] - 2t(d/dt), which satisfy [P.sub.n](1) = 1. Apart from a multiplicative constant, the '[[epsilon].sup.3]-Legendre function' [P.sub.n]([[epsilon].sup.3]*) : [OMEGA] [right arrow] [-1, +1], [xo] [??] [P.sub.n]([[epsilon].sup.3] x [xi]), [xi] [member of] [OMEGA], is the only spherical harmonic of degree n which is invariant under orthogonal transformations leaving [[epsilon].sup.3] fixed. The linear space [Harm.sub.n] of all spherical harmonics of degree n is of dimension dim([Harm.sub.n]) = 2n + 1. Thus, there exist 2n + 1 linearly independent spherical harmonics [Y.sub.n,1], ..., [Y.sub.n,2n+1] in [Harm.sub.n]. Throughout this paper we assume this system to be orthonormal in the sense of the [L.sup.2]([OMEGA])-inner product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(d[omega] denotes the surface element). An outstanding result of the theory of spherical harmonics is the addition theorem

[2n+1.summation over (k=1)][Y.sub.n,k]([xi])[Y.sub.n,k]([eta]) = 2n+1/4[pi] [P.sub.n]([xi] x [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA].

The close connection between the orthogonal invariance and the addition theorem is established by the Funk-Hecke formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

H [member of] [L.sup.1][-1, +1], [xi], [zeta] [member of] [OMEGA], where the Legendre transform LT : H [right arrow] (LT)(H),H [member of] [L.sup.1][-1, 1], is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The sequence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is called the symbol of H. For more details about the theory of spherical harmonics the reader is referred, for example, to [12, 6].

We let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Of course,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that

dim([Harm.sub.0,...,m]) = [m.summation over (n=0)](2n + 1) = [(m + 1).sup.2].

As it is well known, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined by

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is the reproducing kernel of the space [Harm.sub.0,...,m] with respect to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Moreover it is worth mentioning that

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all [xi] [member of] [OMEGA] and all Y [member of] [Harm.sub.0,...,m].

In what follows we are mainly interested in results for the Hilbert space [L.sup.2]([OMEGA]) equipped with the inner product [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Any function of class [L.sup.2]([OMEGA]) of the form [H.sub.[xi]] : [OMEGA] [right arrow] R, [eta] [??] [H.sub.[xi]]([eta]) = H([xi] x [eta]), [eta] [member of] [OMEGA], is called a [xi]-zonal function on [omega]. Zonal functions are constant on the sets of all [eta] [member of] [OMEGA], with [xi] x [eta] = h, h [member of] [-1, +1]. The set of all [xi]-zonal functions is isomorphic to the set of functions H : [-1, +1] [right arrow] R. This gives rise to consider the space [L.sup.2][-1, +1] with norm defined correspondingly by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

as subspace of [L.sup.2]([OMEGA]).

The spherical Fourier transform H [??] (FT)(H), H [member of] [L.sup.2]([OMEGA]), is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

FT forms a mapping from [L.sup.2]([OMEGA]) onto the space [l.sup.2](N) of all sequences [{[W.sub.n,k]}.sub.(n,k)[member of]N] satisfying

[summation over ((n,k)[member of]N)][W.sup.2.sub.n,k] = [[infinity].summation over (n=0)][2n+1.summation over (k=1)][W.sup.2.sub.n,k] < [infinity],

where we have used the abbreviation

N = {(n, k)|n = 0, 1, ..., k = 1, ..., 2n + 1}.

The series

[summation over ((n,k)[member of]N)]F^(n, k)[Y.sub.n,k]

is called the spherical Fourier expansion of F (with Fourier coefficients F^(n, k), (n, k) [member of] N). For all F [member of] [L.sup.2]([OMEGA]) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

4. Convolutions. A kernel H : [OMEGA] x [OMEGA] [right arrow] R is called a square-summable product kernel if H is of the form

H([xi], [eta]) = [summation over ((n,k)[member of]N)]H^(n, k)[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])

such that

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In the case of rotational invariance of the kernel H, i.e. H^(n, k) = H^(n) for n = 0, 1, ..., k = 1, ..., 2n + 1, the last condition is equivalent to [l.sup.2](N)-summability (cf. [6]), i.e.

[[infinity].summation over (n=0)] 2n+1/4[pi] [(H^(n)).sup.2] < [infinity].

Assume that H is a square-summable product kernel and F [member of] [L.sup.2]([OMEGA]). Then the convolution of H against F is defined by

H * F = [[integral].sub.[OMEGA]] H (x, [eta]) F([eta]) d[omega]([eta]).

Two important properties of spherical convolutions should be mentioned: (i) If F [member of] [L.sup.2]([OMEGA]) and H is a square summable product kernel, then H * F is of class [L.sup.2]([OMEGA]). (ii) If [H.sub.1],[H.sub.2] are square-summable product kernels, then the convolution of [H.sub.1],[H.sub.2] defined by

([H.sub.1] * [H.sub.2])([xi] x [zeta]) = [[integral].sub.[OMEGA]][H.sub.1] ([xi] x [zeta])[H.sub.2]([xi] x [zeta]) d[omega]([eta])

is a square-summable product kernel with

([H.sub.1] * [H.sub.2])^(n, k) = [H^.sub.1](n, k) [H^.sub.2](n, k), (n, k) [member of] N.

We usually write [H.sup.(2)] = H * H to indicate the convolution of H with itself. [H.sup.(2)] is said to be the (second) iterated kernel of H. More general, [H.sup.(p)] = [H.sup.(p-1)] * H for p = 2, 3, ... and [H.sup.(1)] = H. Obviously,

([H.sup.(p)])^(n, k) = [(H^(n, k)).sup.p], (n, k) [member of] N, p [member of] N.

5. Multiscale Approximation. Next we consider a strictly monotonically decreasing sequence [{[[rho].sub.j]}.sub.j[member of]Z] of real numbers that satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(for example, [[rho].sub.j] = [2.sup.-j], j [member of] Z). The sequence [{[[rho].sub.j]}.sub.j[member of]Z] can be understood as a subdivision of the 'scale interval' (0,[infinity]) into a countable, strictly monotonically decreasing sequence. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be a family of square-summable product kernels satisfying the condition [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all j [member of] Z. Then, the family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of operators [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], defined by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is called a singular integral in [L.sup.2]([OMEGA]). [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is called kernel of the singular integral. If [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a kernel of a singular integral satisfying the conditions:

(i) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all (n, k) [member of] N,

(ii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all j [member of] Z and (n, k) [member of] N,

(iii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all (n, k) [member of] N,

then the corresponding singular integral [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called an approximate identity in [L.sup.2]([OMEGA]). It is known (see e.g. [6]) that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

provided that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is an approximate identity.

Our results immediately lead us to the following statement:

LEMMA 5.1. Assume that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a kernel constituting an approximate identity in [L.sup.2]([OMEGA]). Then the limit relation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

holds for all F [member of] [L.sup.2]([OMEGA]).

For J [member of] Z we set

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider a kernel [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] constituting an approximate identity in [L.sup.2]([OMEGA]). Assume that F is of class [L.sup.2]([OMEGA]). Then a simple calculation shows us that for all N [member of] N and J [member of] Z,

(5.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we have introduced the family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] by the spectral refinement condition

(5.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words,

(5.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

j [member of] Z, ([xi], [eta]) [member of] [OMEGA] x [OMEGA]. Hence, letting N tend to infinity we get the following multiscale reconstruction formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5.4)

for every J [member of] Z (in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]). Moreover, we find

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

hence,

(5.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Combining (5.4) and (5.5) we finally obtain the following result:

LEMMA 5.2. The multiscale representation of F [member of] [L.sup.2]([OMEGA])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

holds in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] provided that the so-called 'scaling function' [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] forms an approximate identity in [L.sup.2]([OMEGA]) and the 'wavelet' [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] satisfies the difference equation (5.3).

By construction, the wavelet theory leads to a partition of unity as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all (n, k) [member of] N. The class [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of all functions P [member of] [L.sup.2]([OMEGA]) of the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called the scale space of level j (with respect to the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], whereas the class [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of all functions Q [member of] [L.sup.2]([OMEGA]) of the representation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called the detail space of level j (with respect to the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. It is easily seen from (5.1) that

(5.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all j [member of] Z. But it should be remarked that the sum (5.6) generally is neither direct nor orthogonal (note that an orthogonal decomposition is given by the Shannon scaling function). The equation (5.6) can be interpreted in the following way: The set [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] contains a filtered ('smoothed') version of a function belonging to [L.sup.2]([OMEGA]). The lower the scale, the stronger the intensity of smoothing. By adding 'details' contained in the detail space [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the space [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is created, which consists of a filtered ('smoothed') version at resolution j + 1 (see [6] for more details of spherical theory, [5] for the application harmonic theory, and [11] for application in gravimetry).

Finally, it is worth mentioning that the scale spaces satisfy the following properties:

(i) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(ii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A collection of subspaces of [L.sup.2]([OMEGA]) satisfying (i) and (ii) is called a multiresolution analysis of [L.sup.2]([OMEGA]).

6. Examples. Singular integrals on the sphere are of basic interest in geomathematical applications. We essentially distinguish two types, namely bandlimited and non-bandlimited singular integrals.

6.1. Bandlimited Singular Integrals.

The Shannon Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with a strictly monotonically decreasing sequence of integers [{[[rho].sub.j]}.sub.j[member of]Z] satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(for example: [[rho].sub.j] = [2.sup.-j]).

The Smoothed Shannon Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [{[[rho].sub.j]}.sub.j[member of]Z] is defined as in the Shannon case and [{[[sigma].sub.j]}.sub.j[member of]Z] is a strictly monotonically decreasing sequence of integers satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [[tau].sub.j] is a strictly monotonically decreasing and continuous function of class C[[[sigma].sup.-1.sub.j], [[sigma].sup.-1.sub.j]], j [greater than or equal to] 0, such that

[[tau].sub.j]([[sigma].sup.-1.sub.j]) = 1, [[tau].sub.j] ([[sigma].sup.-1.sub.j]) = 0,

for example [[tau].sub.j](t) = 2 - [2.sup.-j]t with [[rho].sub.j] = [2.sup.-j-1] and [[sigma].sub.j] = [2.sup.-j].

6.2. Non-bandlimited Singular Integrals.

The Abel-Poisson Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Tikhonov Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where {[[sigma].sub.n,k]}(n,k)[member of]N is a sequence satisfying the following conditions:

(i) [[sigma].sub.n,k] [not equal to] 0 for all (n, k) [member of] N,

(ii) {[[sigma].sub.n,k]}(n,k)[member of]N is [l.sup.2](N)-summable, i.e.

[summation over ((n,k)[member of]N)] [[sigma].sup.2.sub.n,k] < [infinity].

It should be remarked that, even if we only consider bandlimited functions as data, one might want to use non-bandlimited kernel functions for the multiscale analysis. This is because of the characteristics of the kernel functions in the space domain. Non-bandlimited kernel functions show less oscillations in the space domain than bandlimited kernels, which is surely a desirable feature if data of high degree and order is to be analyzed.

7. Spectral Signal-to-Noise Response. Geoscientists mostly think of their measurements (after possible linearization) as a linear operator on an 'input signal' F producing an 'output signal' G

(7.1) [LAMBDA]F = G,

where [LAMBDA] is an operator mapping the space [L.sup.2]([OMEGA]) into itself such that

[LAMBDA][Y.sub.n,k] = [LAMBDA]^(n, k)[Y.sub.n,k], (n, k) [member of] N,

where the so-called symbol [{[LAMBDA]^(n, k)}].sub.(n,k)[member of]N] is the sequence of the real numbers [LAMBDA]^(n, k). Different linear operators [LAMBDA], of course, are characterized by different sequences [{[LAMBDA]^(n, k)}.sub.(n,k)[member of]N].

The 'amplitude spectrum' [{G^(n, k)}.sub.(n,k)[member of]N] of the response of [LAMBDA] is described in terms of the amplitude spectrum of functions (signals) by a simple multiplication by the 'transfer' [LAMBDA]^(n, k). For a large number of problems in geophysics and geodesy [LAMBDA] is a rotation-invariant operator, i.e. [LAMBDA]^(n, k) = [LAMBDA]^(n) for all (n, k) [member of] N.

7.1. Noise Model. Thus far only a (deterministic) function model has been discussed. If a comparison of the 'output function' with the actual value were done, discrepancies would be observed. A mathematical description of these discrepancies has to follow the laws of probability theory in a stochastic model (see e.g. [13]).

Usually the observations are not looked upon as a time series, but rather a function [??] on the sphere [OMEGA]('~' for stochastic). According to this approach we assume that

[??] = G + [??],

where [??] is the observation noise. Moreover, in our approach motivated by information in satellite technology, we suppose the covariance to be known

Cov [[??]([xi]), [??]([eta])] = E [[??]([xi]), [??](eta)] = K([xi], [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA],

where the following conditions are imposed on the symbol [{K^(n, k)}.sub.(n,k)[member of]N] of the kernel function ([xi], [eta]) [??] K([xi], [eta]) = [[summation].sup.[infinity].sub.n=0][[summation].sup.2n+1.sub.k=1] [Y.sub.n,k]([xi])[Y.sub.n,k]([eta])

(C1) K^(n, k) [greater than or equal to] 0 for all (n, k) [member of] N,

(C2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Condition (C2), indeed, implies in the case of rotational-invariance, i.e.

K^(n, k) = K^(n), n = 0, 1, ..., k = 1, ..., 2n + 1,

the [l.sup.(2)](N)-summability of the symbol [{K^(n, k)}.sub.(n,k)[member of]N], i.e.

[summation over ((n,k)[member of]N)] 1/4[pi] [(K^(n, k)).sup.2] = [[infinity].summation over (n=0)] 2n+1/4[pi] [(K^(n)).sup.2] < [infinity].

7.2. Degree Variances. Since any 'output function' (output signal) can be expanded into an orthogonal series of surface spherical harmonics

[??] = [??] = [summation over ((n,k)[member of]N)][LAMBDA]^(n, k)[??]^(n, k)[Y.sub.n,k] = [summation over ((n,k)[member of]N)] [??]^(n, k)[Y.sub.n,k]

in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], we get a spectral representation of the form

[??]^(n, k) = ([??])^(n, k) = [LAMBDA]^(n, k)[??]^(n, k), (n, k) [member of] N.

The signal degree and order variance of [??] = [??] is defined by

[Var.sub.n,k]([??]) = [(([??])^(n, k)).sup.2] = [[integral].sub.[OMEGA]][[integral].sub.[OMEGA]]([??])([xi])([??]) ([eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])d[omega]([eta])d[omega]([eta]).

Correspondingly, the signal degree variances of [??] = [??] are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

n = 0, 1, .... According to Parseval's identity we clearly have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Physical devices do not transmit spherical harmonics of arbitrarily high frequency without severe attenuation. The 'transfer' [LAMBDA]^(n, k) usually tends to zero with increasing n. It follows that the amplitude spectra of the responses (observations) to functions (signals) of finite [L.sup.2]([OMEGA])-energy are negligibly small beyond some finite frequency. Thus, both because of the frequency limiting nature of the used devices and because of the nature of the 'transmitted signals', the geoscientist is soon led to consider bandlimited functions. These are the functions [??] [member of] [L.sup.2]([OMEGA]), whose 'amplitude spectra' vanish for all n > N (N [member of] [N.sub.0], fixed). In other words, [Var.sub.n]([??]) = 0, for all n > N.

7.3. Degree Error Covariances. The error spectral theory is based on the degree and order error covariance

[Cov.sub.n,k](K) = [[integral].sub.[OMEGA]][[integral].sub.[OMEGA]]K ([xi], [eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])d[omega]([xi])d[omega] ([eta]), (n, k) [member of] N,

and the spectral degree error covariance

[Cov.sub.n](K) = [2n+1.summation over (k=1)] [[integral].sub.[OMEGA]] [[integral].sub.[OMEGA]]K([xi], [eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta]) d[omega]([xi])d[omega]([eta]), n [member of] [N.sub.0].

Obviously,

[Cov.sub.n,k](K) = K^(n, k).

In other words, the spectral degree and order error covariance is simply the orthogonal coefficient of the kernel K.

7.4. Examples of Spectral Error Covariances. To make the preceding considerations more concrete we would like to list two particularly important examples:

(1) Bandlimited white noise. Suppose that for some [n.sub.K] [member of] [N.sub.0]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [epsilon] is assumed to be N(0, [[sigma].sup.2])-distributed. The kernel is given by:

K([xi], [eta]) = [[sigma].sup.2]/[([n.sub.K] + 1).sup.2] [[n.sub.K].summation over (n=0)] 2n+1/4[pi] [P.sub.n]([xi] x [eta]).

Note that this sum, apart from a multiplicative constant, may be understood as a truncated Dirac [delta]-functional. It is known (see e.g. [10]) that for ([xi], [eta]) [member of] [OMEGA] x [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(2) Non-bandlimited colored noise. Assume that K : [OMEGA] x [OMEGA] [right arrow] R is given in such a way that K^(n, k) = K^(n) > 0 for an infinite number of pairs (n, k) [member of] N, the integral [[integral].sup.[delta].sub.-1] K(t)dt is sufficiently small (for some [delta] [member of] (1 - [epsilon], 1) for some [epsilon] > 0), and K([xi], [xi]) coincides with [[sigma].sup.2] for all [xi] [member of] [OMEGA].

Geophysically relevant examples are the following kernels:

(i) K([xi], [eta]) = [[sigma].sup.2]/exp(-c) exp(-c([xi] x [eta])), ([xi], [eta]) [member of] [OMEGA] x [OMEGA],

where c is to be understood as the inverse spherical correlation length (first degree Gauss-Markov model).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for some sufficiently large [J.sup.*] [member of] N (model of small correlation length). The family of locally supported singular integrals [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the case k = 0 this example is known as the Haar singular integral (more details about Haar wavelets can be found in [7]).

7.5. Spectral Estimation. Now we are in position to compare the signal spectrum with that of the noise.

Signal and noise spectrum 'intersect' at the so-called degree and order resolution set [N.sub.res] (with [N.sub.res] [subset] N). We distinguish the following cases:

(i) signal dominates noise

[Var.sub.n,k]([??]) [greater than or equal to] [Cov.sub.n,k](K), (n, k) [member of] [N.sub.res],

(ii) noise dominates signal

[Var.sub.n,k]([??]) < [Cov.sub.n,k](K), (n, k) [not member of] [N.sub.res].

Filtering is achieved by convolving a square-summable product kernel H with the 'symbol' [{H^(n, k)}.sub.(n,k)[member of]N] against [??]:

[??] = [[integral].sub.[OMEGA]]H(x, [eta])[??]([eta]) d[omega]([eta])

('^' denotes 'estimated'). In spectral language this reads

(7.2) [??](n, k) = H^(n, k)[??](n, k), (n, k) [member of] N.

Two important types of filtering are as follows:

(i) Spectral thresholding

(7.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [I.sub.A] denotes the indicator function of the set A. This approach represents a 'keep or kill' filtering, where the signal dominated coefficients are filtered by a square-summable product kernel, and the noise dominated coefficients are set to zero. This thresholding can be thought of as a non-linear operator on the set of coefficients, resulting in a set of estimated coefficients.

As a special filter we mention the ideal low-pass (Shannon) filter H of the form

(7.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In that case all 'frequencies' (n, k) [member of] [N.sub.res] are allowed to pass, whereas all other frequencies are completely eliminated.

(ii) Wiener-Kolmogorov filtering. Now we choose

(7.5) [??] = [summation over ((n,k)[member of]N)] H^(n)([??])^(n, k) [Y.sub.n,k]

with

(7.6) H^(n) = [Var.sub.n]([??])/[Var.sub.n]([??])+[Cov.sub.n](K), n [member of] [N.sub.0].

This filter produces an optimal weighting between signal and noise (provided that complete independence of signal and noise is assumed). Note the similarity to the rotational-invariant Tikhonov singular integral in (6.1).

8. Multiscale Signal-to-Noise Response. Consider a sequence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of square-summable product kernels constituting an approximate identity in [L.sup.2]([OMEGA]). Then we have verified that an 'output signal' [??] [member of] [L.sup.2]([OMEGA]) of an operator [LAMBDA] can be represented in multiscale approximation as follows

(8.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the equality is understood in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]-sense. The identity (8.1) is equivalent to the identity

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for every [J.sub.0] [member of] Z.

8.1. Scale and Position Variances. Denote by [L.sup.2](Z x [OMEGA]) the space of functions H : Z x [OMEGA] [right arrow] R satisfying

[[infinity].summation over (j=-[infinity])][(H(j, [eta])).sup.2] d[omega]([eta]) < [infinity].

[L.sup.2](Z x [OMEGA]) is a Hilbert space equipped with the inner product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

corresponding to the norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider a family of square-summable product kernels [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] constituting an approximate identity in [L.sup.2]([OMEGA]). From the multiscale formulation of an 'output function' [??] = [??] [member of] [L.sup.2]([OMEGA]) we immediately obtain (cf. [9])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The signal scale and space variance of [??] at position [eta] [member of] [OMEGA] and scale j [member of] Z is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The signal scale variance of [??] is defined by

[Var.sub.j]([??]) = [[integral].sub.[OMEGA]][Var.sub.j;[eta]]([??])d [omega]([eta]).

Obviously, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Expressed in the spectral language of spherical harmonics we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

With the convention Z = Z x [OMEGA] we are formally able to write

(8.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We mention that the Beppo-Levi Theorem justifies to interchange integration and summation. Note that all integrations are understood in the Lebesgue-sense.

8.2. Noise Model. Let K : ([xi], [eta]) [??] K([xi], [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA], satisfy the conditions (C1) and (C2) stated in Section 7.1. The error theory is based on the scale and space error covariance at [eta] [member of] [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The scale error covariance is defined by

[Cov.sub.j](K) = [[integral].sub.[OMEGA]][Cov.sub.j;[eta]](K)d[omega] ([eta]).

We obviously have in spectral language

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is clear from our stochastic model, i.e. from the special representation of the covariance as a product kernel, that the scale error covariance cannot be dependent on the position [eta] [member of] [OMEGA]. This is also indicated by the spectral formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Our error model is particularly useful for the proper handling of the satellite data in Earth's gravitational or magnetic potential determination (see [5] and the references therein).

8.3. Scale and Space Estimation. Signal and noise scale 'intersect' at the so-called scale and space resolution set [Z.sub.res] with [Z.sub.res] [subset] Z. We distinguish the following cases:

(i) signal dominates noise

[Var.sub.j;[eta]]([??]) [greater than or equal to] [Cov.sub.j;[eta]](K), (j; [eta]) [member of] [Z.sub.res].

(ii) noise dominates signal

[Var.sub.j;[eta]]([??]) < [Cov.sub.j;[eta]](K), (j; [eta]) [not member not] [Z.sub.res].

Via the multiscale reconstruction formula the (filtered) J-level approximation of the error-affected function [??] reads as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For J sufficiently large, [??] is well-represented by [([??]).sub.J]. In other words, all the higher-level coefficients are regarded as being negligible, i.e. [([??]).sub.J] [equivalent] [??].

9. Selective Multiscale Reconstruction. Similar to what is known in taking Fourier approximation, we are able to take multiscale approximation by replacing the (unknown) error-free function [LAMBDA]F of the representation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by (an estimate from) the error-affected function [??] such as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

J > [J.sub.0]. Computing the following coefficients at position [eta] [member of] [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

will, of course, require adequate methods of numerical integration on the sphere.

9.1. Numerical Integration on the Sphere. Many integration techniques are known from the literature (for a survey on approximate integration on the sphere see, for example, [6] and the references therein). In what follows we base integration on the approximate formulae associated to known weights [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and knots [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

('[equivalent]' always means that the error is assumed to be negligible). An example (cf. [9]) is equidistribution (i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]).

9.2. A Pyramid Scheme. Next we deal with some aspects of scientific computing. We are interested in a pyramid scheme for the (approximate) recursive computation of the integrals [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for j = [J.sub.0], ..., J - 1.

What we are going to realize is a tree algorithm (pyramid scheme) with the following ingredients: Starting from a sufficiently large J such that

(9.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we want to show that the coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (being, of course, dependent on the functiong_F under consideration) can be calculated such that the following statements hold true:

(i) The vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are obtainable by recursion from the values [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(ii) For j = [J.sub.0], ..., J

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For j = [J.sub.0], ..., J - 1

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Our considerations are divided into two parts, viz. the initial step concerning the scale level J and the pyramid step establishing the recursion relation:

The Initial Step. For a suitably large integer [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is sufficiently close to ([??])([eta]) for all [eta] [member of] [OMEGA]. Formally, the kernel [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] replaces the Dirac-functional [delta] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[delta]([xi], [eta]) = [[delta].sub.[epsilon]]([eta]) = [summation over ((n,k)[member of]N)][Y.sub.n,k]([xi])[Y.sub.n,k]([eta]ta)

and the series has to be understood in distributional sense. The formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are the reason why the coefficients for the initial step, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are assumed to be simply given in the form

(9.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Pyramid Step. The essential idea for the development of a pyramid scheme is the existence of kernel functions [[XI].sub.j] : [OMEGA] x [OMEGA] [right arrow] R such that

(9.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[[XI].sub.j] [equivalent] [[XI].sub.j+1] * [[XI].sub.j] (9.4)

for j = [J.sub.0], ..., J.

Note that for bandlimited scaling functions the kernels [[XI].sub.j], j = [J.sub.0], ..., J, may be chosen to be the reproducing kernels of the finite-dimensional scale spaces [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (cf. (3.1)), whereas in the non-bandlimited case [[XI].sub.j], j = [J.sub.0], ..., J, may be chosen such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Observing our approximate integration formulae we obtain in connection with relation (9.3)

(9.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now it follows by use of our approximate integration formulae and the assumption (9.4) that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, the coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be calculated recursively starting from the data [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for the initial level [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be deduced recursively from [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], etc. Moreover, the coefficients are independent of the special choice of the kernel (Observe that (9.5) is equivalent to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for n = 0, 1, ..., k = 1, ..., 2n + 1). This finally leads us to the formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] given by (9.2) and (9.5). In the bandlimited case (with [[XI].sub.j] chosen as indicated above) the sign "[equivalent]" can be replaced by "=" provided that spherical harmonic exact integration formulae of suitable degree are used (cf. [5]).

This recursion procedure leads us to the following decomposition scheme:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], ... allow the following reconstruction scheme of [??]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Once again it is worth mentioning that the coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] do not depend on the special choice of the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Moreover, the coefficients can be used to calculate the wavelet transforms [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for j = [J.sub.0], ..., J - 1 and all [eta] [member of] [OMEGA].

10. Scale Thresholding. Since the large 'true' coefficients are the ones that should be included in a selective reconstruction, in estimating an unknown function it is natural to include only coefficients larger than some specified threshold value.

In our context a 'larger' coefficient is taken to mean one that satisfies for j = [J.sub.0], ..., J and i = 1, ..., [N.sub.j]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Remark 10.1. In particular for "bandlimited white noise" of the form

K([eta], [xi]) = K([eta] x [xi]) = [[sigma].sup.2]/4[pi] [P.sub.0]([eta] x [xi]) = [[sigma].sup.2]/4[pi],

([eta], [xi]) [member of] [OMEGA] x [OMEGA] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (i.e. equidistributions), we find

[([k.sup.j.sub.i]).sup.2] = 2[square root of ([pi])]/[N.sub.j] [([[XI]^.sub.j](0, 1)).sup.2], j = [J.sub.0], ..., J, i = 1, ..., [N.sub.j].

For the given threshold values [k.sup.j.sub.i] such an estimator can be written in explicit form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, the large coefficients (relative to the threshold [k.sup.j.sub.i], i = 1, ..., [N.sub.j], j = [J.sub.0], ..., J - 1) are kept intact and the small coefficients are set to zero. Motivated by our former results the thresholding will be performed on [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [[??].sub.j;[eta]], j = [J.sub.0], ..., J - 1. The thresholding estimators of the true coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can thus be written in the form

(10.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the function [[delta].sup.hard.sub.[lambda]] is the hard thresholding function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The 'keep or kill' hard thresholding operation is not the only reasonable way of estimating the coefficients. Recognizing that each coefficient [[??].sub.j;[eta]] consists of both a signal portion and a noise portion, it might be desirable to attempt to isolate the signal contribution by removing the noisy part. This idea leads to the soft thresholding function (confer the considerations by [2, 3])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which can also be used in the identities (10.1) stated above. When soft thresholding is applied to a set of empirical coefficients, only coefficients greater than the threshold (in absolute value) are included, but their values are 'shrunk' toward zero by an amount equal to the threshold [lambda].

Summarizing all our results we finally obtain the following thresholding multiscale estimator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In doing so [([??]).sub.J] first is approximated by a thresholded [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], which represents the smooth components of the data. Then the coefficients at higher resolutions are thresholded, so that the noise is suppressed but the fine-scale details are included in the calculation.

11. Example. In order to illustrate the effectiveness of our multiscale denoising technique we present a simple example using synthetic geomagnetic data. It is clear that the method presented here can also be applied to non-synthetic data, but then it is disproportionately more difficult to compare the results with expected outcomes. For the purpose of our example we introduce geomagnetic coordinates X, Y and Z. X denotes the so-called north-, Y the east- and Z the downward-component. Using spherical polar coordinates and identifying 0 degree longitude with Greenwich and 0 degree latitude with the equator we end up with the following correspondence:

X [left and right arrow] [[epsilon].sup.t],

Y [left and right arrow] [[epsilon].sup.[psi]],

Z [left and right arrow] -[[epsilon].sup.r],

where [[epsilon].sup.t], [[epsilon].sup.[phi]] and [[epsilon].sup.r] are the usual unit vectors in spherical polar coordinates (for explicit representations see e.g. [6]). This means that X,Y and Z form a local triad with X always pointing towards the geographic northpole, Y pointing into the geographic east direction and Z always being directed towards the Earth's body.

From a bandlimited (up to degree and order 12, which is realistic for the geomagnetic main field) geomagnetic potential due to [1] we calculated the corresponding gradient field in geomagnetic coordinates, i.e. north ([[epsilon].sup.t]), east ([[epsilon].sup.[psi]]) and downward (-[[epsilon].sup.r]) components, which gave as noise-free data (note that we have used the low degree model just for illustrational purposes; it is clear that multiscale techniques are also valuable tools for handling phenomena of especially high degree). We then added some bandlimited white noise with variance [sigma] and bandwidth [n.sub.K] of approximately 0.9 and 60, respectively (see section 7.4). This resulted in noise of the order of magnitude 10.sup.0] [nT] in a field of the order of magnitude [10.sup.4] [nT]. It should be noted that, when looking at the pictures, the noise is not constant at the poles as one should expect it to be. This is due to our routine of adding the noise to the synthetic data. However, our results are not influenced by this, since during the process of decomposition and reconstruction each data point of the rectangular domain is weighted by integration weights due to [4]. These weights are constructed such that the poles do not contribute to the whole integration.

The noise signal then was decomposed and reconstructed using Shannon wavelets up to scale 4 (see construction principles in Section 6.1, where we have chosen [[rho].sub.j] = [2.sup.-j] and j = 1, 2, 3, 4). During the reconstruction process only those wavelet coefficients containing a predominant amount of the clear signal were used in accordance to our considerations in section 8.3. Fig. 11.1 shows the -[[epsilon].sup.r] component of the unnoised data, while Fig. 11.2 shows the absolute values of the added noise.

Figs. 11.3 and 11.4 show the denoised -[[epsilon].sup.r] component and the corresponding absolute error with respect to the unnoised data. Using our multiscale denoising technique the root-mean-square error of the noised data (w.r.t. the clear data), [([DELTA][[epsilon].sup.r.sub.noised]).sub.rms] = 1.13 [nT], has been reduced to [([DELTA][[epsilon].sup.r.sub.denoised]).sub.rms] = 0.35 [nT], which is an improvement of about 60 per cent.

Comparing Figs. 11.2 and 11.4 it can be seen how the rough structure of the noise has been smoothed out by the denoising process and how the peaks have been reduced throughout the whole data set. This example shows the functionality of our approach.

[FIGURE 11.1 OMITTED]

[FIGURE 11.2 OMITTED]

[FIGURE 11.3 OMITTED]

[FIGURE 11.4 OMITTED]

Acknowledgements. The support by German Research Foundation (DFG contract No. FR 761-10-1) is gratefully acknowledged.

REFERENCES

[1] C.J. CAIN, D.R. SCHMITZ, AND L. MUTH, Small-scale features in the earth's magnetic field observed by MAGSAT. Journal of Geophysical Research, Vol. 89, NO. B2 (1984), pp. 1070-1076.

[2] D.L. DONOHO AND I.M. JOHNSTONE, Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81 (1994), pp. 425-455.

[3] --, Adapting to unknown smoothness via wavelet shrinkage, J. Amer. Statistical Association, 90 (1995), pp. 1200-1224.

[4] J.R. DRISCOLL AND D.M. HEALY, Computing Fourier transforms and convolutions on the 2-sphere. Adv. Appl. Math., 15 (1994), pp. 202-250.

[5] W. FREEDEN, Multiscale Modelling of Spaceborne Geodata, B.G. Teubner, Stuttgart, Leipzig, 1999.

[6] W. FREEDEN, T. GERVENS, AND M. SCHREINER, Constructive Approximation on the Sphere (With Applications to Geomathematics), Oxford Science Publications, Clarendon, 1998.

[7] W. FREEDEN AND K. HESSE, On the multiscale solution of satellite problems by use of locally supported kernel functions corresponding to equidistributed data on spherical orbits. AGTM-Report 233, Laboratory of Technomathematics, Geomathematics Group, University of Kaiserslautern, 2000.

[8] W. FREEDEN, V. MICHEL, Constructive Approximation and Numerical Methods in Geodetic Research Today--an Attempt at a Categorization Based on an Uncertainty Principle, J. of Geodesy, 73 (1999), pp.452-465.

[9] W. FREEDEN, V. MICHEL, AND M. STENGER, Multiscale signal-to-noise thresholding, Acta Geod. Geoph. Hung., Vol 36(1) (2001), pp. 55-86.

[10] N.N. LEBEDEW, Spezielle Funktionen und ihre Anwendung, BI-Wissenschaftsverlag, Mannheim, Wien, Zurich, 1973.

[11] V. MICHEL, A Multiscale Method for the Gravimetry Problem--Theoretical and Numerical Aspects of Harmonic and Anharmonic Modelling. Doctoral Thesis, University of Kaiserslautern, Geomathematics Group, Shaker Verlag, Aachen, 1999.

[12] C. MULLER, Spherical Harmonics. Lecture Notes in Mathematics 17, Springer, Berlin, Heidelberg, 1966.

[13] R.T. OGDEN, Essential Wavelets for Statistical Applications and Data Analysis, Birkhauser, Boston, Basel, Berlin, 1997.

[14] R. RUMMEL, Spherical Spectral Properties of the Earth's Gravitational Potential and its First and Second Derivatives, Lecture Notes in Earth Sciences, 65 (1997), pp. 359-404.

[15] J.B. WEAVER, X. YANSUN, D.M. HEALY JR., AND L.D. CROMWELL, Filtering noise from images with wavelet transforms, Magnetic Resonance in Medicine, 24 (1991), pp. 288-295.

W. FREEDEN AND T. MAIER ([dagger])

* Received November 29, 2000. Accepted for publication May 31, 2001. Communicated by Sven Ehrich.

([dagger]) University of Kaiserslautern, Geomathematics Group, 67653 Kaiserslautern, P.O. Box 3049, Germany, email freeden@mathematik.uni-kl.de, tmaier@mathematik.uni-kl.de, wwwhttp://www.mathematik.uni-kl.de/~wwwgeo; Correspondence to W. Freeden

Key words. spherical wavelet theory, scalar multiscale approximation, pyramid scheme, spectral and multiscale variance-covariance model, hard and soft thresholding.

AMS subject classifications. 33C55, 42C40, 62-07, 65T60, 86A25.

1. Introduction. While standard Fourier methods in terms of spherical harmonics are very successful at picking out frequencies from a spherical signal, they are utterly incapable of dealing properly with data changing on small spatial scales. This fact has been well-known for years. To improve the applicability of Fourier analysis, various methods such as 'windowed Fourier transform' have been developed on the sphere to modify the usual Fourier procedure to allow analysis of the frequency content of a signal at each position (cf. [6, 8]). However, the amount of localization in space and in frequency is not completely satisfactory in the Fourier transform and its windowed variant. For example, geopotentials refer to a certain combination of frequencies, and the frequencies themselves are spatially changing. This space evolution of the frequencies is not reflected in the Fourier transform in terms of non-space-localizing spherical harmonics. Even the windowed Fourier transform contains information about frequencies over a certain area of positions instead of showing how the frequencies vary in space. With spherical wavelets, the amount of localization in space and in frequency is automatically adapted, in that only a narrow space-window is needed to examine high-frequency content, but a wide space-window is allowed when investigating low-frequency phenomena. The basic framework of this approach has been provided by the spherical wavelet theory developed by the Geomathematics Group at the University of Kaiserslautern during the last years (see http://www.mathematik.uni-kl.de/~wwwgeo/pub1.html)

When dealing with real geophysically relevant data it should be kept in mind that each measurement does not really give the value of the observable under consideration but that--at least to some extent--the data are contaminated with noise. That is, in order to successfully improve geomathematical field modeling, one main aspect is to extract the true portion of the observable from the actual signal. In consequence, a particular emphasis lies on the subject of denoising. This endeavor is precisely the goal in statistical function estimation. Here, the interest is to 'smooth' the noisy data in order to obtain an estimate of the underlying function. The Euclidean theory of wavelets provides signal processors with new, fast tools that are well-suited for denoising signals (for a survey the reader is e.g. referred to [13] and the references therein).

The objective of this article is to introduce multiscale signal-to-noise thresholding and to provide the wavelet oriented basis of denoising spherical data sets. First, we develop the corresponding theory of denoising spherical functions by non-rotation invariant wavelets (cf. [9] for the rotation-invariant case). With the basic introduction at hand, selective thresholding within a pyramid scheme is presented. The thresholding scheme is designed to distinguish between coefficients which contribute significantly to the signal, and those which are negligible. It should be noted that our approach is essentially influenced by the concept of sparse wavelet representations in Euclidean spaces (cf. [15, 2, 3]) and the stochastic model used in satellite geodesy (see e.g. [14]). Using a multiscale approach we are thus able to include detail information of small spatial extent while suppressing the noise in the approximation appropriately. A simple example of denoising geomagnetic field data will be given as an illustration.

2. Preliminaries. Let [R.sup.3] denote three-dimensional Euclidean space. For x, y [member of] [R.sup.3], x = [([x.sub.1], [x.sub.2], [x.sub.3]).sup.T], y = [([y.sub.1], [y.sub.2], [y.sub.3]).sup.T], the inner product is defined as usual by

x * y = [x.sup.T] y = [3.summation over (i=1)][x.sub.i][y.sub.i].

For all elements x [member of] [R.sup.3], x = [([x.sub.1], [x.sub.2], [x.sub.3]).sup.T], different from the origin, we have

x = r[xi], r = |x| = [square root of (x * x)] = [x.sup.2.sub.1] + [x.sup.2.sub.2] + [x.sup.2.sub.3],

where [xi] = [([[xi].sub.1], [[xi].sub.2], [[xi].sub.3]).sup.T] is the uniquely determined directional unit vector of x. The unit sphere in [R.sup.3] is denoted by [OMEGA]. If the vectors [[epsilon].sup.1], [[epsilon].sup.2], [[epsilon].sup.3] form the canonical orthonormal basis in [R.sup.3], the points [xi] [member of] [OMEGA] may be represented in polar coordinates by

[xi] = t[[epsilon].sup.3] + [square root of (1 - [t.sup.2])] (cos [psi][[epsilon].sup.1] + sin [psi][[epsilon].sup.2]), t = cos v, v [member of] [0, [pi]], [psi] [member of] [0, 2[pi]).

3. Spherical Harmonics. The spherical harmonics [Y.sub.n] of degree n are defined as the everywhere on [OMEGA] infinitely differentiable eigenfunctions of the Beltrami operator [[DELTA].sup.*] corresponding to the eigenvalues ([[DELTA].sup.*])^(n) = -n(n + 1), n = 0, 1, ..., where the Beltrami-operator is the angular part of the Laplace-operator [DELTA] in [R.sup.3]. As it is well-known, the functions [H.sub.n] : [R.sup.3] [right arrow] R defined by [H.sub.n](x) = [r.sup.n][Y.sub.n]([xi]), x = r[xi], [xi] [member of] [OMEGA], are homogeneous polynomials in rectangular coordinates which satisfy the Laplace-equation [[DELTA].sub.x][H.sub.n](x) = 0, x [member of] [R.sup.3]. Conversely, every homogeneous harmonic polynomial of degree n when restricted to [OMEGA] is a spherical harmonic of degree n. The Legendre polynomials [P.sub.n] : [-1, +1] [right arrow] [-1, +1] are the only everywhere on [-1, +1] infinitely differentiable eigenfunctions of the Legendre-operator (1 - [t.sup.2])[(d/dt).sup.2] - 2t(d/dt), which satisfy [P.sub.n](1) = 1. Apart from a multiplicative constant, the '[[epsilon].sup.3]-Legendre function' [P.sub.n]([[epsilon].sup.3]*) : [OMEGA] [right arrow] [-1, +1], [xo] [??] [P.sub.n]([[epsilon].sup.3] x [xi]), [xi] [member of] [OMEGA], is the only spherical harmonic of degree n which is invariant under orthogonal transformations leaving [[epsilon].sup.3] fixed. The linear space [Harm.sub.n] of all spherical harmonics of degree n is of dimension dim([Harm.sub.n]) = 2n + 1. Thus, there exist 2n + 1 linearly independent spherical harmonics [Y.sub.n,1], ..., [Y.sub.n,2n+1] in [Harm.sub.n]. Throughout this paper we assume this system to be orthonormal in the sense of the [L.sup.2]([OMEGA])-inner product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(d[omega] denotes the surface element). An outstanding result of the theory of spherical harmonics is the addition theorem

[2n+1.summation over (k=1)][Y.sub.n,k]([xi])[Y.sub.n,k]([eta]) = 2n+1/4[pi] [P.sub.n]([xi] x [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA].

The close connection between the orthogonal invariance and the addition theorem is established by the Funk-Hecke formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

H [member of] [L.sup.1][-1, +1], [xi], [zeta] [member of] [OMEGA], where the Legendre transform LT : H [right arrow] (LT)(H),H [member of] [L.sup.1][-1, 1], is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The sequence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is called the symbol of H. For more details about the theory of spherical harmonics the reader is referred, for example, to [12, 6].

We let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Of course,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

so that

dim([Harm.sub.0,...,m]) = [m.summation over (n=0)](2n + 1) = [(m + 1).sup.2].

As it is well known, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] defined by

(3.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is the reproducing kernel of the space [Harm.sub.0,...,m] with respect to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Moreover it is worth mentioning that

(3.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all [xi] [member of] [OMEGA] and all Y [member of] [Harm.sub.0,...,m].

In what follows we are mainly interested in results for the Hilbert space [L.sup.2]([OMEGA]) equipped with the inner product [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Any function of class [L.sup.2]([OMEGA]) of the form [H.sub.[xi]] : [OMEGA] [right arrow] R, [eta] [??] [H.sub.[xi]]([eta]) = H([xi] x [eta]), [eta] [member of] [OMEGA], is called a [xi]-zonal function on [omega]. Zonal functions are constant on the sets of all [eta] [member of] [OMEGA], with [xi] x [eta] = h, h [member of] [-1, +1]. The set of all [xi]-zonal functions is isomorphic to the set of functions H : [-1, +1] [right arrow] R. This gives rise to consider the space [L.sup.2][-1, +1] with norm defined correspondingly by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

as subspace of [L.sup.2]([OMEGA]).

The spherical Fourier transform H [??] (FT)(H), H [member of] [L.sup.2]([OMEGA]), is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

FT forms a mapping from [L.sup.2]([OMEGA]) onto the space [l.sup.2](N) of all sequences [{[W.sub.n,k]}.sub.(n,k)[member of]N] satisfying

[summation over ((n,k)[member of]N)][W.sup.2.sub.n,k] = [[infinity].summation over (n=0)][2n+1.summation over (k=1)][W.sup.2.sub.n,k] < [infinity],

where we have used the abbreviation

N = {(n, k)|n = 0, 1, ..., k = 1, ..., 2n + 1}.

The series

[summation over ((n,k)[member of]N)]F^(n, k)[Y.sub.n,k]

is called the spherical Fourier expansion of F (with Fourier coefficients F^(n, k), (n, k) [member of] N). For all F [member of] [L.sup.2]([OMEGA]) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

4. Convolutions. A kernel H : [OMEGA] x [OMEGA] [right arrow] R is called a square-summable product kernel if H is of the form

H([xi], [eta]) = [summation over ((n,k)[member of]N)]H^(n, k)[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])

such that

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In the case of rotational invariance of the kernel H, i.e. H^(n, k) = H^(n) for n = 0, 1, ..., k = 1, ..., 2n + 1, the last condition is equivalent to [l.sup.2](N)-summability (cf. [6]), i.e.

[[infinity].summation over (n=0)] 2n+1/4[pi] [(H^(n)).sup.2] < [infinity].

Assume that H is a square-summable product kernel and F [member of] [L.sup.2]([OMEGA]). Then the convolution of H against F is defined by

H * F = [[integral].sub.[OMEGA]] H (x, [eta]) F([eta]) d[omega]([eta]).

Two important properties of spherical convolutions should be mentioned: (i) If F [member of] [L.sup.2]([OMEGA]) and H is a square summable product kernel, then H * F is of class [L.sup.2]([OMEGA]). (ii) If [H.sub.1],[H.sub.2] are square-summable product kernels, then the convolution of [H.sub.1],[H.sub.2] defined by

([H.sub.1] * [H.sub.2])([xi] x [zeta]) = [[integral].sub.[OMEGA]][H.sub.1] ([xi] x [zeta])[H.sub.2]([xi] x [zeta]) d[omega]([eta])

is a square-summable product kernel with

([H.sub.1] * [H.sub.2])^(n, k) = [H^.sub.1](n, k) [H^.sub.2](n, k), (n, k) [member of] N.

We usually write [H.sup.(2)] = H * H to indicate the convolution of H with itself. [H.sup.(2)] is said to be the (second) iterated kernel of H. More general, [H.sup.(p)] = [H.sup.(p-1)] * H for p = 2, 3, ... and [H.sup.(1)] = H. Obviously,

([H.sup.(p)])^(n, k) = [(H^(n, k)).sup.p], (n, k) [member of] N, p [member of] N.

5. Multiscale Approximation. Next we consider a strictly monotonically decreasing sequence [{[[rho].sub.j]}.sub.j[member of]Z] of real numbers that satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(for example, [[rho].sub.j] = [2.sup.-j], j [member of] Z). The sequence [{[[rho].sub.j]}.sub.j[member of]Z] can be understood as a subdivision of the 'scale interval' (0,[infinity]) into a countable, strictly monotonically decreasing sequence. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be a family of square-summable product kernels satisfying the condition [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all j [member of] Z. Then, the family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of operators [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], defined by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is called a singular integral in [L.sup.2]([OMEGA]). [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is called kernel of the singular integral. If [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a kernel of a singular integral satisfying the conditions:

(i) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all (n, k) [member of] N,

(ii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all j [member of] Z and (n, k) [member of] N,

(iii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for all (n, k) [member of] N,

then the corresponding singular integral [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called an approximate identity in [L.sup.2]([OMEGA]). It is known (see e.g. [6]) that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

provided that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is an approximate identity.

Our results immediately lead us to the following statement:

LEMMA 5.1. Assume that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a kernel constituting an approximate identity in [L.sup.2]([OMEGA]). Then the limit relation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

holds for all F [member of] [L.sup.2]([OMEGA]).

For J [member of] Z we set

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider a kernel [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] constituting an approximate identity in [L.sup.2]([OMEGA]). Assume that F is of class [L.sup.2]([OMEGA]). Then a simple calculation shows us that for all N [member of] N and J [member of] Z,

(5.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we have introduced the family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] by the spectral refinement condition

(5.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words,

(5.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

j [member of] Z, ([xi], [eta]) [member of] [OMEGA] x [OMEGA]. Hence, letting N tend to infinity we get the following multiscale reconstruction formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5.4)

for every J [member of] Z (in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]). Moreover, we find

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

hence,

(5.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Combining (5.4) and (5.5) we finally obtain the following result:

LEMMA 5.2. The multiscale representation of F [member of] [L.sup.2]([OMEGA])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

holds in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] provided that the so-called 'scaling function' [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] forms an approximate identity in [L.sup.2]([OMEGA]) and the 'wavelet' [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] satisfies the difference equation (5.3).

By construction, the wavelet theory leads to a partition of unity as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all (n, k) [member of] N. The class [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of all functions P [member of] [L.sup.2]([OMEGA]) of the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called the scale space of level j (with respect to the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], whereas the class [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of all functions Q [member of] [L.sup.2]([OMEGA]) of the representation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is called the detail space of level j (with respect to the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. It is easily seen from (5.1) that

(5.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all j [member of] Z. But it should be remarked that the sum (5.6) generally is neither direct nor orthogonal (note that an orthogonal decomposition is given by the Shannon scaling function). The equation (5.6) can be interpreted in the following way: The set [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] contains a filtered ('smoothed') version of a function belonging to [L.sup.2]([OMEGA]). The lower the scale, the stronger the intensity of smoothing. By adding 'details' contained in the detail space [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the space [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is created, which consists of a filtered ('smoothed') version at resolution j + 1 (see [6] for more details of spherical theory, [5] for the application harmonic theory, and [11] for application in gravimetry).

Finally, it is worth mentioning that the scale spaces satisfy the following properties:

(i) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(ii) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A collection of subspaces of [L.sup.2]([OMEGA]) satisfying (i) and (ii) is called a multiresolution analysis of [L.sup.2]([OMEGA]).

6. Examples. Singular integrals on the sphere are of basic interest in geomathematical applications. We essentially distinguish two types, namely bandlimited and non-bandlimited singular integrals.

6.1. Bandlimited Singular Integrals.

The Shannon Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with a strictly monotonically decreasing sequence of integers [{[[rho].sub.j]}.sub.j[member of]Z] satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(for example: [[rho].sub.j] = [2.sup.-j]).

The Smoothed Shannon Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [{[[rho].sub.j]}.sub.j[member of]Z] is defined as in the Shannon case and [{[[sigma].sub.j]}.sub.j[member of]Z] is a strictly monotonically decreasing sequence of integers satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [[tau].sub.j] is a strictly monotonically decreasing and continuous function of class C[[[sigma].sup.-1.sub.j], [[sigma].sup.-1.sub.j]], j [greater than or equal to] 0, such that

[[tau].sub.j]([[sigma].sup.-1.sub.j]) = 1, [[tau].sub.j] ([[sigma].sup.-1.sub.j]) = 0,

for example [[tau].sub.j](t) = 2 - [2.sup.-j]t with [[rho].sub.j] = [2.sup.-j-1] and [[sigma].sub.j] = [2.sup.-j].

6.2. Non-bandlimited Singular Integrals.

The Abel-Poisson Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Tikhonov Singular Integral. The family [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where {[[sigma].sub.n,k]}(n,k)[member of]N is a sequence satisfying the following conditions:

(i) [[sigma].sub.n,k] [not equal to] 0 for all (n, k) [member of] N,

(ii) {[[sigma].sub.n,k]}(n,k)[member of]N is [l.sup.2](N)-summable, i.e.

[summation over ((n,k)[member of]N)] [[sigma].sup.2.sub.n,k] < [infinity].

It should be remarked that, even if we only consider bandlimited functions as data, one might want to use non-bandlimited kernel functions for the multiscale analysis. This is because of the characteristics of the kernel functions in the space domain. Non-bandlimited kernel functions show less oscillations in the space domain than bandlimited kernels, which is surely a desirable feature if data of high degree and order is to be analyzed.

7. Spectral Signal-to-Noise Response. Geoscientists mostly think of their measurements (after possible linearization) as a linear operator on an 'input signal' F producing an 'output signal' G

(7.1) [LAMBDA]F = G,

where [LAMBDA] is an operator mapping the space [L.sup.2]([OMEGA]) into itself such that

[LAMBDA][Y.sub.n,k] = [LAMBDA]^(n, k)[Y.sub.n,k], (n, k) [member of] N,

where the so-called symbol [{[LAMBDA]^(n, k)}].sub.(n,k)[member of]N] is the sequence of the real numbers [LAMBDA]^(n, k). Different linear operators [LAMBDA], of course, are characterized by different sequences [{[LAMBDA]^(n, k)}.sub.(n,k)[member of]N].

The 'amplitude spectrum' [{G^(n, k)}.sub.(n,k)[member of]N] of the response of [LAMBDA] is described in terms of the amplitude spectrum of functions (signals) by a simple multiplication by the 'transfer' [LAMBDA]^(n, k). For a large number of problems in geophysics and geodesy [LAMBDA] is a rotation-invariant operator, i.e. [LAMBDA]^(n, k) = [LAMBDA]^(n) for all (n, k) [member of] N.

7.1. Noise Model. Thus far only a (deterministic) function model has been discussed. If a comparison of the 'output function' with the actual value were done, discrepancies would be observed. A mathematical description of these discrepancies has to follow the laws of probability theory in a stochastic model (see e.g. [13]).

Usually the observations are not looked upon as a time series, but rather a function [??] on the sphere [OMEGA]('~' for stochastic). According to this approach we assume that

[??] = G + [??],

where [??] is the observation noise. Moreover, in our approach motivated by information in satellite technology, we suppose the covariance to be known

Cov [[??]([xi]), [??]([eta])] = E [[??]([xi]), [??](eta)] = K([xi], [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA],

where the following conditions are imposed on the symbol [{K^(n, k)}.sub.(n,k)[member of]N] of the kernel function ([xi], [eta]) [??] K([xi], [eta]) = [[summation].sup.[infinity].sub.n=0][[summation].sup.2n+1.sub.k=1] [Y.sub.n,k]([xi])[Y.sub.n,k]([eta])

(C1) K^(n, k) [greater than or equal to] 0 for all (n, k) [member of] N,

(C2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Condition (C2), indeed, implies in the case of rotational-invariance, i.e.

K^(n, k) = K^(n), n = 0, 1, ..., k = 1, ..., 2n + 1,

the [l.sup.(2)](N)-summability of the symbol [{K^(n, k)}.sub.(n,k)[member of]N], i.e.

[summation over ((n,k)[member of]N)] 1/4[pi] [(K^(n, k)).sup.2] = [[infinity].summation over (n=0)] 2n+1/4[pi] [(K^(n)).sup.2] < [infinity].

7.2. Degree Variances. Since any 'output function' (output signal) can be expanded into an orthogonal series of surface spherical harmonics

[??] = [??] = [summation over ((n,k)[member of]N)][LAMBDA]^(n, k)[??]^(n, k)[Y.sub.n,k] = [summation over ((n,k)[member of]N)] [??]^(n, k)[Y.sub.n,k]

in the sense of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], we get a spectral representation of the form

[??]^(n, k) = ([??])^(n, k) = [LAMBDA]^(n, k)[??]^(n, k), (n, k) [member of] N.

The signal degree and order variance of [??] = [??] is defined by

[Var.sub.n,k]([??]) = [(([??])^(n, k)).sup.2] = [[integral].sub.[OMEGA]][[integral].sub.[OMEGA]]([??])([xi])([??]) ([eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])d[omega]([eta])d[omega]([eta]).

Correspondingly, the signal degree variances of [??] = [??] are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

n = 0, 1, .... According to Parseval's identity we clearly have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Physical devices do not transmit spherical harmonics of arbitrarily high frequency without severe attenuation. The 'transfer' [LAMBDA]^(n, k) usually tends to zero with increasing n. It follows that the amplitude spectra of the responses (observations) to functions (signals) of finite [L.sup.2]([OMEGA])-energy are negligibly small beyond some finite frequency. Thus, both because of the frequency limiting nature of the used devices and because of the nature of the 'transmitted signals', the geoscientist is soon led to consider bandlimited functions. These are the functions [??] [member of] [L.sup.2]([OMEGA]), whose 'amplitude spectra' vanish for all n > N (N [member of] [N.sub.0], fixed). In other words, [Var.sub.n]([??]) = 0, for all n > N.

7.3. Degree Error Covariances. The error spectral theory is based on the degree and order error covariance

[Cov.sub.n,k](K) = [[integral].sub.[OMEGA]][[integral].sub.[OMEGA]]K ([xi], [eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta])d[omega]([xi])d[omega] ([eta]), (n, k) [member of] N,

and the spectral degree error covariance

[Cov.sub.n](K) = [2n+1.summation over (k=1)] [[integral].sub.[OMEGA]] [[integral].sub.[OMEGA]]K([xi], [eta])[Y.sub.n,k]([xi])[Y.sub.n,k]([eta]) d[omega]([xi])d[omega]([eta]), n [member of] [N.sub.0].

Obviously,

[Cov.sub.n,k](K) = K^(n, k).

In other words, the spectral degree and order error covariance is simply the orthogonal coefficient of the kernel K.

7.4. Examples of Spectral Error Covariances. To make the preceding considerations more concrete we would like to list two particularly important examples:

(1) Bandlimited white noise. Suppose that for some [n.sub.K] [member of] [N.sub.0]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [epsilon] is assumed to be N(0, [[sigma].sup.2])-distributed. The kernel is given by:

K([xi], [eta]) = [[sigma].sup.2]/[([n.sub.K] + 1).sup.2] [[n.sub.K].summation over (n=0)] 2n+1/4[pi] [P.sub.n]([xi] x [eta]).

Note that this sum, apart from a multiplicative constant, may be understood as a truncated Dirac [delta]-functional. It is known (see e.g. [10]) that for ([xi], [eta]) [member of] [OMEGA] x [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(2) Non-bandlimited colored noise. Assume that K : [OMEGA] x [OMEGA] [right arrow] R is given in such a way that K^(n, k) = K^(n) > 0 for an infinite number of pairs (n, k) [member of] N, the integral [[integral].sup.[delta].sub.-1] K(t)dt is sufficiently small (for some [delta] [member of] (1 - [epsilon], 1) for some [epsilon] > 0), and K([xi], [xi]) coincides with [[sigma].sup.2] for all [xi] [member of] [OMEGA].

Geophysically relevant examples are the following kernels:

(i) K([xi], [eta]) = [[sigma].sup.2]/exp(-c) exp(-c([xi] x [eta])), ([xi], [eta]) [member of] [OMEGA] x [OMEGA],

where c is to be understood as the inverse spherical correlation length (first degree Gauss-Markov model).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for some sufficiently large [J.sup.*] [member of] N (model of small correlation length). The family of locally supported singular integrals [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the case k = 0 this example is known as the Haar singular integral (more details about Haar wavelets can be found in [7]).

7.5. Spectral Estimation. Now we are in position to compare the signal spectrum with that of the noise.

Signal and noise spectrum 'intersect' at the so-called degree and order resolution set [N.sub.res] (with [N.sub.res] [subset] N). We distinguish the following cases:

(i) signal dominates noise

[Var.sub.n,k]([??]) [greater than or equal to] [Cov.sub.n,k](K), (n, k) [member of] [N.sub.res],

(ii) noise dominates signal

[Var.sub.n,k]([??]) < [Cov.sub.n,k](K), (n, k) [not member of] [N.sub.res].

Filtering is achieved by convolving a square-summable product kernel H with the 'symbol' [{H^(n, k)}.sub.(n,k)[member of]N] against [??]:

[??] = [[integral].sub.[OMEGA]]H(x, [eta])[??]([eta]) d[omega]([eta])

('^' denotes 'estimated'). In spectral language this reads

(7.2) [??](n, k) = H^(n, k)[??](n, k), (n, k) [member of] N.

Two important types of filtering are as follows:

(i) Spectral thresholding

(7.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [I.sub.A] denotes the indicator function of the set A. This approach represents a 'keep or kill' filtering, where the signal dominated coefficients are filtered by a square-summable product kernel, and the noise dominated coefficients are set to zero. This thresholding can be thought of as a non-linear operator on the set of coefficients, resulting in a set of estimated coefficients.

As a special filter we mention the ideal low-pass (Shannon) filter H of the form

(7.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In that case all 'frequencies' (n, k) [member of] [N.sub.res] are allowed to pass, whereas all other frequencies are completely eliminated.

(ii) Wiener-Kolmogorov filtering. Now we choose

(7.5) [??] = [summation over ((n,k)[member of]N)] H^(n)([??])^(n, k) [Y.sub.n,k]

with

(7.6) H^(n) = [Var.sub.n]([??])/[Var.sub.n]([??])+[Cov.sub.n](K), n [member of] [N.sub.0].

This filter produces an optimal weighting between signal and noise (provided that complete independence of signal and noise is assumed). Note the similarity to the rotational-invariant Tikhonov singular integral in (6.1).

8. Multiscale Signal-to-Noise Response. Consider a sequence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] of square-summable product kernels constituting an approximate identity in [L.sup.2]([OMEGA]). Then we have verified that an 'output signal' [??] [member of] [L.sup.2]([OMEGA]) of an operator [LAMBDA] can be represented in multiscale approximation as follows

(8.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the equality is understood in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]-sense. The identity (8.1) is equivalent to the identity

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for every [J.sub.0] [member of] Z.

8.1. Scale and Position Variances. Denote by [L.sup.2](Z x [OMEGA]) the space of functions H : Z x [OMEGA] [right arrow] R satisfying

[[infinity].summation over (j=-[infinity])][(H(j, [eta])).sup.2] d[omega]([eta]) < [infinity].

[L.sup.2](Z x [OMEGA]) is a Hilbert space equipped with the inner product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

corresponding to the norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider a family of square-summable product kernels [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] constituting an approximate identity in [L.sup.2]([OMEGA]). From the multiscale formulation of an 'output function' [??] = [??] [member of] [L.sup.2]([OMEGA]) we immediately obtain (cf. [9])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The signal scale and space variance of [??] at position [eta] [member of] [OMEGA] and scale j [member of] Z is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The signal scale variance of [??] is defined by

[Var.sub.j]([??]) = [[integral].sub.[OMEGA]][Var.sub.j;[eta]]([??])d [omega]([eta]).

Obviously, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Expressed in the spectral language of spherical harmonics we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

With the convention Z = Z x [OMEGA] we are formally able to write

(8.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We mention that the Beppo-Levi Theorem justifies to interchange integration and summation. Note that all integrations are understood in the Lebesgue-sense.

8.2. Noise Model. Let K : ([xi], [eta]) [??] K([xi], [eta]), ([xi], [eta]) [member of] [OMEGA] x [OMEGA], satisfy the conditions (C1) and (C2) stated in Section 7.1. The error theory is based on the scale and space error covariance at [eta] [member of] [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The scale error covariance is defined by

[Cov.sub.j](K) = [[integral].sub.[OMEGA]][Cov.sub.j;[eta]](K)d[omega] ([eta]).

We obviously have in spectral language

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is clear from our stochastic model, i.e. from the special representation of the covariance as a product kernel, that the scale error covariance cannot be dependent on the position [eta] [member of] [OMEGA]. This is also indicated by the spectral formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Our error model is particularly useful for the proper handling of the satellite data in Earth's gravitational or magnetic potential determination (see [5] and the references therein).

8.3. Scale and Space Estimation. Signal and noise scale 'intersect' at the so-called scale and space resolution set [Z.sub.res] with [Z.sub.res] [subset] Z. We distinguish the following cases:

(i) signal dominates noise

[Var.sub.j;[eta]]([??]) [greater than or equal to] [Cov.sub.j;[eta]](K), (j; [eta]) [member of] [Z.sub.res].

(ii) noise dominates signal

[Var.sub.j;[eta]]([??]) < [Cov.sub.j;[eta]](K), (j; [eta]) [not member not] [Z.sub.res].

Via the multiscale reconstruction formula the (filtered) J-level approximation of the error-affected function [??] reads as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For J sufficiently large, [??] is well-represented by [([??]).sub.J]. In other words, all the higher-level coefficients are regarded as being negligible, i.e. [([??]).sub.J] [equivalent] [??].

9. Selective Multiscale Reconstruction. Similar to what is known in taking Fourier approximation, we are able to take multiscale approximation by replacing the (unknown) error-free function [LAMBDA]F of the representation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

by (an estimate from) the error-affected function [??] such as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

J > [J.sub.0]. Computing the following coefficients at position [eta] [member of] [OMEGA]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

will, of course, require adequate methods of numerical integration on the sphere.

9.1. Numerical Integration on the Sphere. Many integration techniques are known from the literature (for a survey on approximate integration on the sphere see, for example, [6] and the references therein). In what follows we base integration on the approximate formulae associated to known weights [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and knots [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

('[equivalent]' always means that the error is assumed to be negligible). An example (cf. [9]) is equidistribution (i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]).

9.2. A Pyramid Scheme. Next we deal with some aspects of scientific computing. We are interested in a pyramid scheme for the (approximate) recursive computation of the integrals [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for j = [J.sub.0], ..., J - 1.

What we are going to realize is a tree algorithm (pyramid scheme) with the following ingredients: Starting from a sufficiently large J such that

(9.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we want to show that the coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (being, of course, dependent on the functiong_F under consideration) can be calculated such that the following statements hold true:

(i) The vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are obtainable by recursion from the values [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(ii) For j = [J.sub.0], ..., J

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For j = [J.sub.0], ..., J - 1

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Our considerations are divided into two parts, viz. the initial step concerning the scale level J and the pyramid step establishing the recursion relation:

The Initial Step. For a suitably large integer [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is sufficiently close to ([??])([eta]) for all [eta] [member of] [OMEGA]. Formally, the kernel [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] replaces the Dirac-functional [delta] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[delta]([xi], [eta]) = [[delta].sub.[epsilon]]([eta]) = [summation over ((n,k)[member of]N)][Y.sub.n,k]([xi])[Y.sub.n,k]([eta]ta)

and the series has to be understood in distributional sense. The formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are the reason why the coefficients for the initial step, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], are assumed to be simply given in the form

(9.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The Pyramid Step. The essential idea for the development of a pyramid scheme is the existence of kernel functions [[XI].sub.j] : [OMEGA] x [OMEGA] [right arrow] R such that

(9.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[[XI].sub.j] [equivalent] [[XI].sub.j+1] * [[XI].sub.j] (9.4)

for j = [J.sub.0], ..., J.

Note that for bandlimited scaling functions the kernels [[XI].sub.j], j = [J.sub.0], ..., J, may be chosen to be the reproducing kernels of the finite-dimensional scale spaces [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (cf. (3.1)), whereas in the non-bandlimited case [[XI].sub.j], j = [J.sub.0], ..., J, may be chosen such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Observing our approximate integration formulae we obtain in connection with relation (9.3)

(9.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now it follows by use of our approximate integration formulae and the assumption (9.4) that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, the coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be calculated recursively starting from the data [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for the initial level [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be deduced recursively from [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], etc. Moreover, the coefficients are independent of the special choice of the kernel (Observe that (9.5) is equivalent to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for n = 0, 1, ..., k = 1, ..., 2n + 1). This finally leads us to the formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] given by (9.2) and (9.5). In the bandlimited case (with [[XI].sub.j] chosen as indicated above) the sign "[equivalent]" can be replaced by "=" provided that spherical harmonic exact integration formulae of suitable degree are used (cf. [5]).

This recursion procedure leads us to the following decomposition scheme:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], ... allow the following reconstruction scheme of [??]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Once again it is worth mentioning that the coefficient vectors [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] do not depend on the special choice of the scaling function [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Moreover, the coefficients can be used to calculate the wavelet transforms [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for j = [J.sub.0], ..., J - 1 and all [eta] [member of] [OMEGA].

10. Scale Thresholding. Since the large 'true' coefficients are the ones that should be included in a selective reconstruction, in estimating an unknown function it is natural to include only coefficients larger than some specified threshold value.

In our context a 'larger' coefficient is taken to mean one that satisfies for j = [J.sub.0], ..., J and i = 1, ..., [N.sub.j]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Remark 10.1. In particular for "bandlimited white noise" of the form

K([eta], [xi]) = K([eta] x [xi]) = [[sigma].sup.2]/4[pi] [P.sub.0]([eta] x [xi]) = [[sigma].sup.2]/4[pi],

([eta], [xi]) [member of] [OMEGA] x [OMEGA] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (i.e. equidistributions), we find

[([k.sup.j.sub.i]).sup.2] = 2[square root of ([pi])]/[N.sub.j] [([[XI]^.sub.j](0, 1)).sup.2], j = [J.sub.0], ..., J, i = 1, ..., [N.sub.j].

For the given threshold values [k.sup.j.sub.i] such an estimator can be written in explicit form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In other words, the large coefficients (relative to the threshold [k.sup.j.sub.i], i = 1, ..., [N.sub.j], j = [J.sub.0], ..., J - 1) are kept intact and the small coefficients are set to zero. Motivated by our former results the thresholding will be performed on [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [[??].sub.j;[eta]], j = [J.sub.0], ..., J - 1. The thresholding estimators of the true coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can thus be written in the form

(10.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the function [[delta].sup.hard.sub.[lambda]] is the hard thresholding function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The 'keep or kill' hard thresholding operation is not the only reasonable way of estimating the coefficients. Recognizing that each coefficient [[??].sub.j;[eta]] consists of both a signal portion and a noise portion, it might be desirable to attempt to isolate the signal contribution by removing the noisy part. This idea leads to the soft thresholding function (confer the considerations by [2, 3])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which can also be used in the identities (10.1) stated above. When soft thresholding is applied to a set of empirical coefficients, only coefficients greater than the threshold (in absolute value) are included, but their values are 'shrunk' toward zero by an amount equal to the threshold [lambda].

Summarizing all our results we finally obtain the following thresholding multiscale estimator

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In doing so [([??]).sub.J] first is approximated by a thresholded [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], which represents the smooth components of the data. Then the coefficients at higher resolutions are thresholded, so that the noise is suppressed but the fine-scale details are included in the calculation.

11. Example. In order to illustrate the effectiveness of our multiscale denoising technique we present a simple example using synthetic geomagnetic data. It is clear that the method presented here can also be applied to non-synthetic data, but then it is disproportionately more difficult to compare the results with expected outcomes. For the purpose of our example we introduce geomagnetic coordinates X, Y and Z. X denotes the so-called north-, Y the east- and Z the downward-component. Using spherical polar coordinates and identifying 0 degree longitude with Greenwich and 0 degree latitude with the equator we end up with the following correspondence:

X [left and right arrow] [[epsilon].sup.t],

Y [left and right arrow] [[epsilon].sup.[psi]],

Z [left and right arrow] -[[epsilon].sup.r],

where [[epsilon].sup.t], [[epsilon].sup.[phi]] and [[epsilon].sup.r] are the usual unit vectors in spherical polar coordinates (for explicit representations see e.g. [6]). This means that X,Y and Z form a local triad with X always pointing towards the geographic northpole, Y pointing into the geographic east direction and Z always being directed towards the Earth's body.

From a bandlimited (up to degree and order 12, which is realistic for the geomagnetic main field) geomagnetic potential due to [1] we calculated the corresponding gradient field in geomagnetic coordinates, i.e. north ([[epsilon].sup.t]), east ([[epsilon].sup.[psi]]) and downward (-[[epsilon].sup.r]) components, which gave as noise-free data (note that we have used the low degree model just for illustrational purposes; it is clear that multiscale techniques are also valuable tools for handling phenomena of especially high degree). We then added some bandlimited white noise with variance [sigma] and bandwidth [n.sub.K] of approximately 0.9 and 60, respectively (see section 7.4). This resulted in noise of the order of magnitude 10.sup.0] [nT] in a field of the order of magnitude [10.sup.4] [nT]. It should be noted that, when looking at the pictures, the noise is not constant at the poles as one should expect it to be. This is due to our routine of adding the noise to the synthetic data. However, our results are not influenced by this, since during the process of decomposition and reconstruction each data point of the rectangular domain is weighted by integration weights due to [4]. These weights are constructed such that the poles do not contribute to the whole integration.

The noise signal then was decomposed and reconstructed using Shannon wavelets up to scale 4 (see construction principles in Section 6.1, where we have chosen [[rho].sub.j] = [2.sup.-j] and j = 1, 2, 3, 4). During the reconstruction process only those wavelet coefficients containing a predominant amount of the clear signal were used in accordance to our considerations in section 8.3. Fig. 11.1 shows the -[[epsilon].sup.r] component of the unnoised data, while Fig. 11.2 shows the absolute values of the added noise.

Figs. 11.3 and 11.4 show the denoised -[[epsilon].sup.r] component and the corresponding absolute error with respect to the unnoised data. Using our multiscale denoising technique the root-mean-square error of the noised data (w.r.t. the clear data), [([DELTA][[epsilon].sup.r.sub.noised]).sub.rms] = 1.13 [nT], has been reduced to [([DELTA][[epsilon].sup.r.sub.denoised]).sub.rms] = 0.35 [nT], which is an improvement of about 60 per cent.

Comparing Figs. 11.2 and 11.4 it can be seen how the rough structure of the noise has been smoothed out by the denoising process and how the peaks have been reduced throughout the whole data set. This example shows the functionality of our approach.

[FIGURE 11.1 OMITTED]

[FIGURE 11.2 OMITTED]

[FIGURE 11.3 OMITTED]

[FIGURE 11.4 OMITTED]

Acknowledgements. The support by German Research Foundation (DFG contract No. FR 761-10-1) is gratefully acknowledged.

REFERENCES

[1] C.J. CAIN, D.R. SCHMITZ, AND L. MUTH, Small-scale features in the earth's magnetic field observed by MAGSAT. Journal of Geophysical Research, Vol. 89, NO. B2 (1984), pp. 1070-1076.

[2] D.L. DONOHO AND I.M. JOHNSTONE, Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81 (1994), pp. 425-455.

[3] --, Adapting to unknown smoothness via wavelet shrinkage, J. Amer. Statistical Association, 90 (1995), pp. 1200-1224.

[4] J.R. DRISCOLL AND D.M. HEALY, Computing Fourier transforms and convolutions on the 2-sphere. Adv. Appl. Math., 15 (1994), pp. 202-250.

[5] W. FREEDEN, Multiscale Modelling of Spaceborne Geodata, B.G. Teubner, Stuttgart, Leipzig, 1999.

[6] W. FREEDEN, T. GERVENS, AND M. SCHREINER, Constructive Approximation on the Sphere (With Applications to Geomathematics), Oxford Science Publications, Clarendon, 1998.

[7] W. FREEDEN AND K. HESSE, On the multiscale solution of satellite problems by use of locally supported kernel functions corresponding to equidistributed data on spherical orbits. AGTM-Report 233, Laboratory of Technomathematics, Geomathematics Group, University of Kaiserslautern, 2000.

[8] W. FREEDEN, V. MICHEL, Constructive Approximation and Numerical Methods in Geodetic Research Today--an Attempt at a Categorization Based on an Uncertainty Principle, J. of Geodesy, 73 (1999), pp.452-465.

[9] W. FREEDEN, V. MICHEL, AND M. STENGER, Multiscale signal-to-noise thresholding, Acta Geod. Geoph. Hung., Vol 36(1) (2001), pp. 55-86.

[10] N.N. LEBEDEW, Spezielle Funktionen und ihre Anwendung, BI-Wissenschaftsverlag, Mannheim, Wien, Zurich, 1973.

[11] V. MICHEL, A Multiscale Method for the Gravimetry Problem--Theoretical and Numerical Aspects of Harmonic and Anharmonic Modelling. Doctoral Thesis, University of Kaiserslautern, Geomathematics Group, Shaker Verlag, Aachen, 1999.

[12] C. MULLER, Spherical Harmonics. Lecture Notes in Mathematics 17, Springer, Berlin, Heidelberg, 1966.

[13] R.T. OGDEN, Essential Wavelets for Statistical Applications and Data Analysis, Birkhauser, Boston, Basel, Berlin, 1997.

[14] R. RUMMEL, Spherical Spectral Properties of the Earth's Gravitational Potential and its First and Second Derivatives, Lecture Notes in Earth Sciences, 65 (1997), pp. 359-404.

[15] J.B. WEAVER, X. YANSUN, D.M. HEALY JR., AND L.D. CROMWELL, Filtering noise from images with wavelet transforms, Magnetic Resonance in Medicine, 24 (1991), pp. 288-295.

W. FREEDEN AND T. MAIER ([dagger])

* Received November 29, 2000. Accepted for publication May 31, 2001. Communicated by Sven Ehrich.

([dagger]) University of Kaiserslautern, Geomathematics Group, 67653 Kaiserslautern, P.O. Box 3049, Germany, email freeden@mathematik.uni-kl.de, tmaier@mathematik.uni-kl.de, wwwhttp://www.mathematik.uni-kl.de/~wwwgeo; Correspondence to W. Freeden

Printer friendly Cite/link Email Feedback | |

Author: | Freeden, W.; Maier, T. |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Date: | Jul 1, 2002 |

Words: | 7899 |

Previous Article: | An algorithm for nonharmonic signal analysis using Dirichlet series on convex polygons. |

Next Article: | A polynomial collocation method for Cauchy singular integral equations over the interval. |