# Wavelet modelling of the gravity field by domain decomposition methods: an example over Japan

Isabelle Panet, Yuki Kuroishi, Matthias Holschneider

► To cite this version:

Isabelle Panet, Yuki Kuroishi, Matthias Holschneider. Wavelet modelling of the gravity field by domain decomposition methods: an example over Japan. Geophysical Journal International, 2011, 184, pp.203-219. 10.1111/j.1365-246X.2010.04840.x. insu-03606512

HAL Id: insu-03606512

https://insu.hal.science/insu-03606512

Submitted on 11 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wavelet modelling of the gravity field by domain decomposition methods: an example over Japan

Isabelle Panet,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Yuki Kuroishi<fnref n="3" /> and Matthias Holschneider<fnref n="4" />

Accepted 2010 October 5. Received 2010 October 1; in original form 2009 October 23

## SUMMARY

With the advent of satellite gravity, large gravity data sets of unprecedented quality at low and medium resolution become available. For local, high resolution field modelling, they need to be combined with the surface gravity data. Such models are then used for various applications, from the study of the Earth interior to the determination of oceanic currents. Here we show how to realize such a combination in a flexible way using spherical wavelets and applying a domain decomposition approach. This iterative method, based on the Schwarz algorithms, allows to split a large problem into smaller ones, and avoids the calculation of the entire normal system, which may be huge if high resolution is sought over wide areas. A subdomain is defined as the harmonic space spanned by a subset of the wavelet family. Based on the localization properties of the wavelets in space and frequency, we define hierarchical subdomains of wavelets at different scales. On each scale, blocks of subdomains are defined by using a tailored spatial splitting of the area. The data weighting and regularization are iteratively adjusted for the subdomains, which allows to handle heterogeneity in the data quality or the gravity variations. Different levels of approximations of the subdomains normals are also introduced, corresponding to building local averages of the data at different resolution levels.

We first provide the theoretical background on domain decomposition methods. Then, we validate the method with synthetic data, considering two kinds of noise: white noise and coloured noise. We then apply the method to data over Japan, where we combine a satellite-based geopotential model, EIGEN-GL04S, and a local gravity model from a combination of land and marine gravity data and an altimetry-derived marine gravity model. A hybrid spherical harmonics/wavelet model of the geoid is obtained at about 15 km resolution and a corrector grid for the surface model is derived.

**Key words:** Wavelet transform; Satellite geodesy; Geopotential theory.

## 1 INTRODUCTION

Knowing the Earth's geoid and gravity field at regional scales is essential for many applications. First, the conversion between GPS-derived and levelled heights becomes possible with the knowledge of the geoid. It allows us to monitor crustal movements over long periods of time beyond the advent of GPS, and thus to better understand crustal activities during the seismic cycle in active areas, such as subduction zones. Moreover, the geoid gives us the reference surface for ocean dynamics. An absolute geoid model with a sufficient accuracy can yield the absolute ocean currents, leading to a better understanding of the oceanic circulation and of the Western Boundary Currents (Gulf Stream, Kuroshio, for instance). Finally, a fine

knowledge of the gravity field gives insights into the structure of the lithosphere in areas of interest, such as rift basins, subduction zones, areas of intraplate volcanism and combined with other geophysical and geodetical observations, contributes to a better understanding of the geodynamic processes at stake.

The geoid is defined as the equipotential surface of the gravity field corresponding to the mean sea surface and an accurate gravity model, if available, can be used to determine the geoid accurately. In this paper, we focus on geoid and gravity field modelling over Japan. Because of the tectonic settings in and around Japan, located in a trench and island arc region where four tectonic plates converge, significant undulations of the gravity field and accordingly the geoid occur in a wide range of spatial scales. Consequently, Japan is one of

<footnote n="1">¹Institut Géographique National, Laboratoire LAREG, ENSG, 6/8 av. Blaise Pascal, Cité Descartes, Champs/Marne, 77455 Marne-la-Vallée Cedex 2, France.
E-mail: isabelle.panet@ensg.eu</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris (IPGP, Université Paris Diderot, CNRS), Bat. Lamarck, Case 7011, 35, rue Hélène Brion, 75205 Paris Cedex 13, France</footnote>

<footnote n="3">³Geospatial Information Authority of Japan, Space Geodesy Research Division, 1 Kitasato, Tsukuba, Ibaraki 305-0811, Japan</footnote>

<footnote n="4">⁴University of Potsdam, Department of Applied Mathematics, Am Neuen Palais 10, D14469, Potsdam, Germany</footnote>
the most challenging areas in the world for gravity modelling, which also makes it an excellent zone for testing new methodologies.

The emergence of dedicated-gravity satellite missions such as GRACE (launched in 2002) and GOCE (launched in 2009) highly improves the accuracy of the Earth's gravity field and geoid at long and medium wavelengths. These missions provide a global coverage of gravity field information with a high absolute accuracy in uniform fashion, to spatial resolutions of about 200 km for GRACE and about 110 km for GOCE.

Merging the data obtained by these missions with dense surface gravity data (land, marine and altimetric data) in an optimal way, in consideration of their respective error characteristics, is essential in order to derive the geoid and gravity models that will be used for geodetical and geophysical purposes. Such combination allows to highlight possible biases of the surface data at larger scales and thus to improve the quality of the surface data and consequently of a resulting geoid/gravity model. Different types of gravity field data sets, with various spatial distributions and spectral contents, can be combined by using a functional representation of the gravity field based on spherical wavelets, as shown by previous work (Freeden et al. 1998; Kusche et al. 1998; Panet et al. 2004, 2006; Schmidt et al. 2005, 2007; Klees & Wittwer 2007; Tenzer & Klees 2008). However, the ever increasing number of data on the one hand, and the large number of wavelets needed to derive high resolution models over wide areas on the other hand, stress the need for numerically efficient methods to compute a wavelet model of the gravity field. Moreover, the heterogeneous characteristics of both the gravity field, rugged in mountainous areas, smoother in the plains, and the data sets, affected by non-stationnary errors, require flexible modelling approaches. For instance, splines models with different parametrizations of the calculations for different areas have been developed by Eicker et al. (2006).

Here, we extend the Poisson wavelet modelling approach developed by Holschneider et al. (2003), Chambodut et al. (2005) and Panet et al. (2006), by introducing the theoretical framework of iterative domain decomposition approaches, that allow to solve a large problem by splitting it into smaller parts, and by applying these approaches to the case of regional wavelet modelling of the gravity field at higher spatial resolution than previously done. We explain how the use of subdomains, corresponding to the linear spans of subsets of wavelets at different scales and in different areas, allows to define data sets and regularization weights in a flexible way.

The structure of the paper is as follows. We first recall how to compute a wavelet model of the gravity field from various data sets. We then focus on the use of iterative domain decomposition methods. Then, we validate our approach with synthetic data. Finally, we apply the method to data over Japan.

# 2 WAVELET REPRESENTATION OF THE GRAVITY FIELD

## 2.1 Wavelet frames

Wavelets are functions localized in both space and frequency. They have been extensively described in the literature (for instance, Holschneider (1995) and Mallat (1999)). The wavelet theory has been extended to the spherical geometry for Earth sciences applications (Freeden et al. 1998; Holschneider et al. 2003). Here we use 3-D Poisson multipole wavelets, defined in spherical geometry. The wavelet $\psi^{a,m,\vec{e}}$ calculated at the point $\vec{x}$ on or outside the sphere of

<figure><img src="image_2.png" /><figcaption>Figure 1. Two Poisson multipole wavelets of order 3 on the sphere.</figcaption></figure>

radius $R$ is expressed as (Holschneider et al. 2003)

$$ \psi_{a,\vec{e}}^{m}(\vec{x}) = N_a \sum_{\ell} (a\ell)^m e^{-a\ell} \left( \frac{R}{\|\vec{x}\|} \right)^{\ell+1} Q_{\ell} \left( \frac{\vec{e}}{\|\vec{e}\|} \cdot \frac{\vec{x}}{\|\vec{x}\|} \right), $$

where $a$ is the scale, $\vec{e}$ is the centre of the multipole inside the Earth, $m$ is its order and $N_a$ is the normalization

$$ N_a = \left[ \int_{S(R)} \psi_{a,\vec{e}}^{m}(\vec{x}) \cdot \psi_{a,\vec{e}}^{m}(\vec{x}) \cdot ds(\vec{x}) \right]^{-\frac{1}{2}} $$

and $Q_{\ell} = (2\ell + 1) \cdot P_{\ell}$, with $P_{\ell}$ the Legendre polynomial of degree $\ell$. This expression can be rewritten in terms of multipoles located in $\lambda \cdot \frac{\vec{e}}{\|\vec{e}\|}$ inside the sphere, where $\lambda = e^{-a}$, according to (Holschneider et al. 2003):

$$ \psi_{a,\vec{e}}^{m} = (-1)^{m+1} a^m \left[ \frac{2}{a} (\lambda \partial_{\lambda})^{m+1} - (\lambda \partial_{\lambda})^m \right] \frac{1}{|\vec{x} - \lambda \frac{\vec{e}}{\|\vec{e}\|}|} $$

Because those wavelets can be identified as equivalent multipolar sources of gravity within the Earth, they are well suited to represent the potential field of the Earth's gravity. A wavelet family is constructed with two parameters: a scale parameter (defining the spatial extent), and a position parameter (defining the location). Fig. 1 shows two examples of such wavelets on the sphere having large and small scale parameters. The good localization of the functions makes them appropriate for local gravity field modelling.

The Earth's gravity potential can be arbitrarily well approximated by a linear combination of wavelets properly sampled in scale and position. We construct a family of multipole wavelets of order 3, following Chambodut et al. (2005) and Panet et al. (2006). This order offers a good compromise between localization in space and in frequency. A selection of scales should be made in order to provide a regular coverage of spectrum. This is obtained for a dyadic sequence of scales. Because the sphere is bounded, the wavelets at different scales are not truly dilated versions of each other. However, they asymptotically behave so when the scale becomes small. On a given scale, all wavelets are rotated copies of each other.

Regarding the position parameter, wavelets are sampled with their central locations at the vertices of a spherical mesh. The number of wavelets must increase as the scale decreases, because the number of degrees of freedom of the modelled spaces increases (see below). When the scale is divided by a factor of two, this dimension is
roughly multiplied by a factor of 4. Therefore, the spherical meshes are roughly four times finer as the scale decreases. The mesh on each scale should be chosen so that the wavelets have a coverage of the sphere as regular as possible. For local modelling application, the target area is limited and we may build simple meshes whose vertices are sampled regularly in longitude and latitude.

A family of wavelets thus sampled forms a frame. A frame is a complete set of functions which, contrary to a basis, may be redundant (i.e. overcomplete). The redundancy may be quantified in terms of the so-called frame bounds. An approximate procedure consists in dimension counting. The spectrum of the wavelets is peaked around its maximum, and most of its energy is contained in a limited spectral band around the maximum, comprised between degrees $l_1$ and $l_2$. Consequently, the wavelets can be approximated with band-limited functions by truncating their spectra to the spectral band between $l_1$ and $l_2$. These truncated wavelets can be represented with surface spherical harmonics with degrees between $l_1$ and $l_2$. The dimension $D$ of the linear span of the surface spherical harmonics with degrees between $l_1$ and $l_2$, denoted $W$, is finite. The frame (over-) completeness is estimated by comparing the number of truncated wavelets with the dimensions of the corresponding $W$ spaces (Holschneider et al. 2003). If the sampling of scales and positions is too dense, there are too many truncated wavelets as compared to the dimension of the band-limited space $W$ and this leads to overcompleteness. In contrast, if the sampling is too sparse, the wavelet set may be uncomplete and is no longer a frame. Panet et al. (2004) and Chambodut et al. (2005) discuss the spectral coverage of this family of wavelets and show that our selection (described in detail later in Section 4.1) is overcomplete at an estimated redundancy of 1.4 to a spatial resolution of 10 km. Note that, because of the overcompleteness of such representations, the spherical harmonics basis remains the most compact representation of the gravity field at a global scale. However, even if they are less compact, representations using local functions such as wavelets are of a great interest not only for regional modelling but also for global modelling, because they allow more easily to avoid local data errors to contaminate the estimation of the gravity field in other areas.

## 2.2 Consideration for application to local gravity field modelling

For the application to local gravity field modelling, we should make some special consideration. Different kinds of data, with different spatial coverage and error characteristics, are to be combined. Among them are gravity anomalies measured locally on the Earth surface and on board satellite measurements at the altitude of the satellite orbits. In the following, we presuppose that there are two data sets available: a set of local high-resolution gravity anomalies and a global geopotential model expanded in spherical harmonics up to degree 120, from dedicated-gravity satellite missions. From the latter we compute and use as a second set of data geopotential values at the ground level in the area considered. They contain much more reliable signals at longer wavelengths than local gravity anomaly data. In order to avoid the leakage inside the target area of possible edge effects related to the geopotential data modelling, we extend the data area of geopotential values by two degrees in each direction (north, south, east and west) as compared to the gravity anomaly data coverage. The width of the extension area is slightly larger than the resolution of the geopotential model. This empirical choice leads to satisfactory results in the target area.

Wavelength longer than the extent of the data coverage shall not be reliably recovered from local data and we use purely the global geopotential model at its lower frequency parts in the resulting model. Then, in the gravity field modelling, we take data residuals from the lower order parts of the geopotential model and combine them by the wavelet-based method, finally deriving a hybrid spherical harmonics/wavelet model. In the representation of the residual gravity field, therefore, we should include wavelets at scales no larger than a half of the computation area size. To better constrain the wavelet coefficients from local data, we limit wavelets to scales no larger than about a fourth of the area width.

## 2.3 Computation of wavelet coefficients

A least-squares inversion of the data sets is applied to the computation of the wavelet coefficients. Since the gravity potential is written as a linear combination of wavelets, the observation equations for each dataset, $p$, are derived from the functional relation between the data type and the gravity potential. In matrix notations, this reads $A_p \cdot x = b_p + e_p + s_p$, where $b_p$ is the measurement vector, $A_p$ the design matrix, $e_p$ the data errors, $s_p$ the projection error and $x$ the coefficients to be determined. The projection error is the component of the signal that cannot be modelled by a linear combination of the wavelets. It should be made as small as possible by a proper selection of the wavelets given the signal characteristics, or by a data pre-processing step. If the wavelet model resolution is smoother than the data resolution, neglecting $s_i$ leads to large residuals and to an aliasing of the higher frequencies of the data in the model. In what follows, we will suppose that the adequacy between the wavelets and the signal characteristics is good, and neglect this term. From the observation equations, we derive the normal equations: for each data set $p$, we have $N_p \cdot x = f_p$, where $N_p = A_p^t \cdot W_p \cdot A_p$ is the normal matrix, $W_p$ is the weight matrix based on the supposed measurement noise, and $f_p = A_p^t \cdot W_p \cdot b_p$ is the associated right-hand side. The normal equations are summed up for all data sets to form the global system and a regularization matrix $K$ is added with a parameter $\lambda$

$$ (N + \lambda K) \cdot x = f \quad (1) $$

with $N = \sum_p N_p$ and $f = \sum_p f_p$. The size of the normal equation system is $n_w \times n_w$, where $n_w$ is the total number of wavelets. We denote $G = N + \lambda K$ the regularized normal matrix.

The regularization is introduced for two reasons. First, the characteristics of the measurements, in terms of noise and spatial distribution, may lead to an ill-posed problem. In this case, a physical *a priori* information is needed to define the regularization matrix, which plays the role of filling the data gaps with *a priori* values and/or filtering out the data noise by assessing the expected amplitude of the physical signal. Usually, a condition expressing the spectral decrease of the field's energy is used (Panet et al. 2004; Chambodut et al. 2005). In the case of satellite data handling, the ill-posedness in the downward continuation may be regularized by the use of a surface data set, and the physical *a priori* is not always necessary. Second, if the wavelets frame is overcomplete, which is the case in our constructions with Poisson wavelets, there exists no unique solution $x$ to the problem, at least in the limit of infinitely many wavelets. For a finite number of wavelets, this non-uniqueness is mirrored by an ill-conditioning of the normal equations. In such cases, one has to add a purely numerical regularization.

In the Bayesian interpretation, the vector of wavelet coefficients $x$ is a random variable and the probability of a model $x$ follows a Gaussian distribution with covariance matrix $K^{-1}$. If the wavelet frame is overcomplete, $K$ is singular and the non-definite inverse
$K^{-1}$ should be replaced by a generalized inverse $K^{+}$ obtained by pseudo-inversion of $K$ (Härdle & Simar 2003). However, the pseudo-inversion involves the calculation of a singular values decomposition of $K$, which is time-consuming for large matrices. Consequently, from a practical point of view, a Tikhonov regularization of $K$ is used to impose a unique solution. This numerical regularization allows the optimal numerical calculation of the physical system by optimizing the trade-off between changing slightly the system to solve and improving its condition number. As the Kaula's rule (Kaula 1966) imposes an important spectral decrease of the gravity field's energy, a numerical Tikhonov regularization of the entire matrix $K$ may mask the physical regularization for the smaller scales components. Indeed, its effect is equivalent to a reduction of the *a priori* spectral decrease of the field's energy over the scales, by flattening the spectrum at the smallest scales. When we split the global problem into smaller, rather decorrelated subproblems, as done in a subdomain approach, we reduce the numerical ill-posedness due to the frame overcompleteness. In each subdomain, the regularization can be locally adjusted and if a physical regularization is needed, it is not hidden by the numerical one.

## 2.4 Interest of domain decomposition approaches

These considerations, together with the numerical challenge of handling the large normal systems that arise when a high-resolution model is sought over a wide area, lead us to develop domain decomposition approaches for the calculation of the wavelet models of the gravity field. The principle of such approaches, widely used for solving partial differential equations over large domains, is to restrict the global problem to smaller subspaces, hereafter called subdomains, solve locally the subproblems, and iterate until the local solutions diffuse to the whole space and the global solution is stable. This methods avoids the computation of the entire normal system, which may be time-consuming. Only small matrices need to be loaded in memory and inverted, which allows to solve large problems in an economic way. Another interest is that the restriction of the normal system to the different subdomains allows to reweight the data sets in a flexible way, taking into account their respective quality in the different wavebands and in the different areas. A subdomain problem may also need considerably less regularization, if no regularization at all, than the global problem. Thus, a subdomain-based parametrization of the regularization allows to extract a maximum of information from the data sets by avoiding a global physical regularization that may be locally inadapted or degraded by the numerical one.

The approach we present in this paper is closely related to another iterative method named 'multigrid', which is an example of domain decomposition approach. Multigrid methods (Wesseling 1991) are based on the resolution of successive projections of the normal system on subdomains defined as coarse or fine grids. On the finer grids, the projected system is solved iteratively, whereas it can be solved exactly on the coarser grids. Multigrid techniques have been applied by Kusche (2001) to the case of gravity field modelling from simulated satellite gravity data. Successively solving the same problem approximated at different resolutions on the different grids, allows to considerably increase the convergence rate as compared to many iterative solvers, which may be slow because of long-wavelengths features in the modelled field. In contrast, the coarse grids allow to efficiently model these large scale components. As the hierarchical resolution of the grids clearly reminds the multiscale resolution of the wavelets, wavelets techniques have been used to construct multigrid iterations in the planar case by identifying the grids with

approximation spaces at different resolutions in the wavelet representation (Keller 2001). Here, we extend Keller's work on the plane to the spherical geometry and relate the iterative multilevel approach to the more general framework of domain decomposition algorithms. The presented approach works 'top/down/top', first computing the lower resolution components of the gravity field (top), progressively refining the spatial resolution of the model (down), and then re-estimating the lower resolution components (top). Data compression at different resolution levels is used in the process, in order to adjust the spectral content of the data to the calculated scale. This is in contrast with usual multiresolution approaches, such as developed by Schmidt et al. (2007), where the gravity components at progressively coarser resolutions are estimated and the calculation is not iterated. If the wavelets at different scales are orthogonal with respect to the quadratic Kaula/Energy form, iterations are not needed because the normal system is block-diagonal for proper data distributions. However, in the non-orthogonal case, the normal system is not perfectly block diagonal and iterations allow to progressively account for the correlations between wavelets at different scales. The order of the calculated scales, from the coarsest to the finest, with the corresponding data compression, allows to derive a good first guess of the solution before refining it.

# 3 DOMAIN DECOMPOSITION METHODS

## 3.1 Definition of the subdomains

Here we describe the iterative domain decomposition approach we apply to calculate a wavelet model of the gravity field over a wide area. The interested reader may refer to Chan & Mathew (1994) and Xu (1992) for a detailed presentation of domain decomposition algorithms. We recall here the main principles.

Let us introduce the space $H$ of the linear span of the wavelets, where $\sum$ is the mean Earth sphere. When we compute a wavelet model by the least-squares adjustment of data, we actually project the data on the image spaces of the design matrices $A_p$, which union is equal to $H$ for the ideal case of an infinitely dense sampling of the gravity potential at the data points. In the case of a discrete data distribution, these image spaces are samplings of $H$. Let us now partition this space $H$ into smaller subspaces $\{H_i; i = 1, \dots, q\}$, called subdomains, so that we have: $H = \sum_{i=1}^q H_i$. The subdomain $H_i$ is the space of all linear combinations of the wavelets in a subset $i$ of the total wavelet family, and the different subdomains may be overlapping or not. By taking advantage of the good localization of the wavelets in space and frequency we define our subdomains in the following way.

We first split $H$ into subspaces $\{H_{[a_i, 0]}; i = 1, \dots, n_{\text{scales}}\}$ corresponding to non-overlapping sets of wavelet scale indices $i$, where $H_{[a_i, 0]}$ is made by all wavelets on scale $a_i$; we refer to them as scale subdomains. If subdomain $H_{[a_i, 0]}$ is still composed of a large number of wavelets, we split it again into subdomains of smaller size $\{H_{[a_i, 0][j, 0]}: j = 1, \dots, n_{\text{blocks}(i)}\}$, again corresponding to non-overlapping sets of wavelets scale/position indices $i/j$. These subdomains are generated by subsets of wavelets at scale $a_i$, and verify

$$H_{[a_i, 0]} = \sum_{j=1}^{n_{\text{blocks}(i)}} H_{[a_i, 0][j, 0]}.$$

They are referred to as block subdomains, and the sum is direct for wavelets with a spectrum that is not band-limited. To be general, let us assume that all scale subdomains have a block subdomain
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Definition of the scale and blocks subdomains.</figcaption></figure>

splitting. We regard $n_{blocks(i)} = 1$ if there is only one block. The number of blocks may vary depending on the scale and, in practice, one partitions the space $H$ so that subspace $H_{[a_i, 0][j, 0]}$ is generated by a not too large number of wavelets (no more than a few thousands). Fig. 2 shows a typical example of scale and block subdomains definitions. In the following, we will call “non-overlapping spaces” the spaces generated by subsets of wavelets with non-overlapping sets of scale/position indices—the obtained spaces being in general not orthogonal with Poisson wavelets—and “overlapping spaces” the spaces generated by subsets of wavelets with overlapping sets of scale or position indices.

We now define an overlapping splitting of subspace $H_{[a_i, 0]}$ as made of a set of overlapping blocks $\{H_{[a_i, 0][j, \delta(i)]}, j = 1, \dots, n_{blocks(i)}\}$. Overlapping blocks $H_{[a_i, 0][j, \delta(i)]}$ are obtained by augmenting the non-overlapping blocks $H_{[a_i, 0][j, 0]}$ with wavelets located in overlap areas of width $\delta(i)$. The size of the overlap area $\delta(i)$ obviously depends on the scale $a_i$. However, for simplicity of notation, we will refer to it as $\delta$ here. It would also be possible to define a set of overlapping scale subdomains (although not done here), by extending the subdomains $H_{[a_i, 0]}$ with the wavelets on adjacent scales $a_{i-1}$ or $a_{i+1}$. Let us note $n_{[a_i, 0][j, \delta]}$ the number of wavelets that generate $H_{[a_i, 0][j, \delta]}$.

Let us now introduce the restriction and extension operators between $H$ and the subdomains. In the following, we will work in the space of wavelet coefficients $l^2(\Gamma)$, where $\Gamma$ is the coefficients index set and $l^2(\Gamma) = \{x; \|x^2\| = \sum_{n \in \Gamma} |x_n|^2 < +\infty\}$ and identify the subdomains $H_{[a_i, 0][j, \delta]}$ with subsets of wavelet coefficients. With each subdomain, we associate a rectangular matrix $R^t_{[a_i, 0][j, \delta]}$ of size $n_w \times n_{[a_i, 0][j, \delta]}$. This matrix is the extension by $n_w - n_{[a_i, 0][j, \delta]}$ counts of zeros of a vector $x$ of size $n_{[a_i, 0][j, \delta]}$ that belongs to $H_{[a_i, 0][j, \delta]}$. Its elements are thus 1 or 0. The transpose of this matrix is the restriction matrix to the subdomain $H_{[a_i, 0][j, \delta]}$. It acts on a vector of size $n_w$, and restricts it on $H_{[a_i, 0][j, \delta]}$ by holding only the entries that belong to $H_{[a_i, 0][j, \delta]}$. We also define restricted extension operators by

$$ \hat{R}^t_{[a_i, 0][j, \delta]} = R^t_{[a_i, 0][j, 0]} \cdot R_{[a_i, 0][j, 0]} \cdot R^t_{[a_i, 0][j, \delta]} $$

Matrix $\hat{R}^t_{[a_i, 0][j, \delta]}$ acts on a vector $x_{[a_i, 0][j, \delta]}$ of size $n_{[a_i, 0][j, \delta]}$ that belongs to $H_{[a_i, 0][j, \delta]}$, restricts it to the non-overlapping subdomain $H_{[a_i, 0][j, 0]}$ by eliminating all elements in the overlap area, and extends it by adding zeros to a vector $x$ of size $n_w$. Thus, in contrast with the simple extension operator $R^t_{[a_i, 0][j, \delta]}$, elements in the overlap area are replaced with zeroes in the case of the restricted extension operator. Finally, we apply weights when extending vectors from $H_{[a_i, 0][j, \delta]}$ to $H$, leading to a weighted, restricted extension operator:

$$ \tilde{R}^t_{[a_i, 0][j, \delta]} = w \cdot R^t_{[a_i, 0][j, 0]} \cdot R_{[a_i, 0][j, 0]} \cdot R^t_{[a_i, 0][j, \delta]} $$

where $w(p)$ is the inverse of the number of overlapping blocks to which the $p$th entry of $x_{[a_i, 0][j, \delta]}$ belongs (it corresponds to the redundancy of the computation due to overlap). Now, instead of solving the whole normal system directly, we will restrict it to the subdomains, derive partial subdomain solutions involving a part of the wavelet coefficient vector $x$, and progressively build the global solution $x$ by iterating the process. Such iterative algorithms are known as Schwarz algorithms (Xu 1992; Chan & Mathew 1994).

## 3.2 Principle of the Schwarz algorithms

Let us start with the non-regularized case, where the normal system to solve is: $N \cdot x = f$.

The restriction of the normal matrix $N$ to the subdomain $H_{[a_i, 0][j, \delta]}$ reads

$$ N_{[a_i, 0][j, \delta]} = R_{[a_i, 0][j, \delta]} \cdot N \cdot R^t_{[a_i, 0][j, \delta]} $$

It is actually a block of $N$, comprising the entries related to the wavelets which linear span is the subspace $H_{[a_i, 0][j, \delta]}$.

Starting from an initial guess of the solution, we build a sequence of estimate $x^{k+1} = x^k + M \cdot (f - N \cdot x^k)$. In the parallel (additive) version of the algorithm, matrix $M$ is

$$ M = \sum_{i=1}^{n_{\text{scales}}} \sum_{j=1}^{n_{\text{blocks}(i)}} \tilde{R}^t_{[a_i, 0][j, \delta]} \cdot N_{[a_i, 0][j, \delta]}^{-1} \cdot R_{[a_i, 0][j, \delta]} $$
In other words, we first restrict the right-hand side to $H_{[a_i, 0][j, \delta]}$ and update it from previous estimates of solution

$$f_{[a_i, 0][j, \delta]}^k = R_{[a_i, 0][j, \delta]} \cdot (f - R_{[a_i, 0][j, \delta]} \cdot N \cdot x^k).$$

Note that the update matrix $R_{[a_i, 0][j, \delta]} \cdot N$ is actually a band of $N$. Then, we solve locally the problem

$$N_{[a_i, 0][j, \delta]} \cdot x_{[a_i, 0][j, \delta]}^{k+1} = f_{[a_i, 0][j, \delta]}^k.$$

The last step is to extend the obtained coefficient vector $x_{[a_i, 0][j, \delta]}^{k+1}$ to the full size $n_w$

$$\mathrm{Ext}[x_{[a_i, 0][j, \delta]}^{k+1}] = \tilde{R}_{[a_i, 0][j, \delta]}^t \cdot x_{[a_i, 0][j, \delta]}^{k+1},$$

where $\mathrm{Ext} [\cdot]$ denotes the extension of a vector to a vector of size $n_w$. At the end of iteration $k+1$, we add all extended subdomain solutions $\mathrm{Ext}[x_{[a_i, 0][j, \delta]}^{k+1}]$ to derive the global vector $x^{k+1}$, and we iterate the process.

In this version of the algorithm, all subdomains are computed simultaneously, the global solution $x^{k+1}$ being updated only after all $(k+1)$th iterations of local solution end. This parallel (also called additive) version of the algorithm converges slower than the sequential version (also called multiplicative), in which the global solution $x^{k+1}$ is updated immediately after a subdomain solution $x_{[a_i, 0][j, \delta]}^{k+1}$ becomes available, and the $(k+1)$th iteration of the next subdomain solution is computed only after the global solution has been updated.

We mix parallel iterations with sequential ones for finding a hybrid algorithm with good convergence and parallelization properties, that is illustrated on Fig. 3. We thus carry out sequential iterations over the scale subdomains and on a given scale, we carry out parallel iterations over the blocks. This is the reason why we do not use overlapping scale subdomains, but make use of overlapping blocks. The overlap speeds up the convergence in the parallel iterations, which converge more slowly, whereas the sequential iterations over scales, which converge faster, can be carried out without overlap.

Such method will converge quickly if the subdomains are not too much correlated, reflecting the sparse structure of the normal matrix. In the case of wavelets that are well localized in both space and frequency, matrix $N$ is sparse and the definition of the subdomains reflects its structure. Therefore, it is quite appropriate for a hybrid algorithm to be used.

### 3.3 Approximations of the subdomain normals

As we progressively refine our solution during the iteration process, we do not need to use exact subdomain normal matrices at the beginning of the process, but apply adequate approximations to the subdomains normals. For that, we follow the approach by Minchev et al. (2009). We construct a 3-D, spherical tiling of the space on and outside of the Earth's sphere and approximate locally the values of the wavelets within each elementary tile, also called a 'cell' in the following. The 3-D tiling is based on a spherical mesh and its radial extension. The 3-D space is divided into radial layers of triangular prism cells sectioned by this mesh, so that the heights of the cells is similar to their width. The spherical mesh is derived by subdividing the facets of a regular icosahedron with respect to the sphere. For more details about this construction, the reader is referred to Chambodut et al. (2005). We obtain a regular tiling with no singularity at the poles. The deviation of mesh lengths is between -25 and +10 per cent of the mean length, with most of lengths within 10 per cent deviation. Depending on the level of subdivision, one obtains deferent level of fineness of the mesh.

There are different ways to approximate the wavelets locally within the cells. Here we use the average value of the wavelet inside the cell, but more generally, one may also consider 3-D Taylor series development at varying orders, leading to different precisions of approximation. The use of order-0 development or

<figure><img src="image_4.png" /><figcaption>Figure 3. Iteration schemes over the subdomains.</figcaption></figure>
averaging actually projects the normal system on a 3-D-grid and is equivalent to data filtering in the cell. Then, the precision of approximation must be adjusted depending on the wavelet scale to be computed. On each scale, we associate a mesh at a resolution similar to the wavelet scale.

As we carry out the iterations from the coarsest wavelet scales to the finest ones, compressing the data at the resolution of each calculated scale allows to derive a good first guess of the solution. Indeed, this process leads, at each step, to a bandpass filtering of the data in the spectral band roughly corresponding to the spectral coverage of the calculated wavelet scale, thus allowing to avoid the leakage of higher or lower resolution signal that may arise due to the spectral overlap of the wavelets at consecutive scales. Then, as we progress in iteration, we gradually refine the approximations of the subdomain normals at all scales. Namely, as the subdomain solutions at finer and finer scales become available, the solutions at the coarser scales have to be updated on a precision level compatible with the desired precision of the overall solution.

## 3.4 Data reweighting per subdomain

The quality of the wavelet model depends on the choice of an appropriate weighting of the different data sets, introduced in the matrices $W_p$. In the Bayesian interpretation, the data set $p$ provides a realization of a gaussian random variable with covariance $W_p^{-1}$ and mean $b_p$. If the observations $b_p$ are exact, then the covariance tends to zero. If they are degraded by a white noise, with rms $(\alpha_p)^{-0.5}$, the weight matrix $W_p$ is diagonal and equal to $\alpha_p I$. However, in addition to purely random errors, real data sets are often affected by unmodelled systematic biases and correlated errors. For instance, altimetry-derived gravity anomalies are known to show local biases close to the coastlines. Around Japan, the presence of such correlated errors has been stressed by Kuroishi & Keller (2005). Moreover, the amplitude of the errors is not perfectly known in general, and a calibration process is required. To deal with these errors, it has been suggested to add parameters to the model (Schaffrin & Iz 2001; van Loon & Kusche 2005), and to use fully populated covariance matrices reflecting the correlated structure of the errors. A simple approach to approximately compensate for the non-random character of the unmodelled errors and recalibrate the errors is to downweight the data sets (Lerch et al. 1991).

To handle properly the data sets errors, it is useful to search for a basis B in which the covariance matrix $W_p$ is diagonal. Such discrete basis is made of sampled functions at the data points, with as many functions as data points, and is orthonormal. The noise coefficients covariances in the space domain (Dirac basis), in $W_p$, are related to the noise coefficients covariances $D_p$ in the basis B as follows:

$$W_p = F_p^t \cdot D_p \cdot F_p,$$

where $F_p$ and $D_p$ are square matrices of size the number of data $n_p$ in the data set $p$. $F_p$ is an orthogonal matrix containing the basis B functions sampled at the data points, and verifying $F_p^t \cdot F_p = I$. The inverse of $W_p$ thus verifies

$$W_p^{-1} = F_p^t \cdot D_p^{-1} \cdot F_p.$$

If the error corresponds to a stationary process, then it is known that $D_p$ is diagonal when B is the discrete Fourier basis. If the error is locally stationary, that is to say, without any abrupt transition inside the different intervals of stationnarity, then $D_p$ will be close to a diagonal matrix if the basis B is made of local Fourier bases over these intervals of stationnarity (Mallat 1999). More generally, the diagonalization of $W_p$ may be obtained through its principal components analysis. In the case of a well chosen fully discrete orthogonal

wavelet basis B, associated to a discrete orthogonal wavelet transform of the noise (DWT), $D_p$ will be close to block-diagonal if the error is actually made of a sum of not too correlated components at different scales and in different areas. Note that a rather regular data sampling is necessary in order to build such basis. One may wonder how different such DWT basis B would be from the Poisson multipole wavelets frame used in our modelling. To account for the high-frequency (resp. low frequency) noise, it would include smaller scales (resp. larger scale) wavelets than the Poisson multipole frame, which is at contrary oversampled at the data points. If the scales and blocks subdomains introduced in Section 3.1 are such that the errors are almost stationary in each subdomain, and not too correlated between different subdomains, one can reasonably assume that $D_p$ should be almost block diagonal with blocks corresponding to these subdomains.

Inserting $W_p^{-1}$ in the normal system then leads to computing the discrete correlations between the Poisson multipole wavelets and the orthogonal wavelets used to model the error, and it is a reasonable approximation to neglect the products between wavelets at different scales and in different blocks. Consequently, the subdomains normal matrices for the data set $p$, $N_{p, \{a_i, 0\}[j, \delta]}$, verify

$$N_{p, \{a_i, 0\}[j, \delta]} = A_{p[a_i, 0][j, \delta]}^t \cdot F_{p[a_i, 0][j, \delta]}^t \cdot D_{p[a_i, 0][j, \delta]}^{-1} \cdot F_{p[a_i, 0][j, \delta]} \cdot A_{p[a_i, 0][j, \delta]}$$

with $F_{p[a_i, 0][j, \delta]}$ (resp. $D_{p[a_i, 0][j, \delta]}^{-1}$) the restriction of $F_p$ (resp. $D_p^{-1}$) to the subdomain $H_{[a_i, 0][j, \delta]}$. Approximating $D_{p[a_i, 0][j, \delta]}$ with a constant diagonal matrix, and with the product $F_{p[a_i, 0][j, \delta]}^t \cdot F_{p[a_i, 0][j, \delta]} \sim \text{Id}_{[a_i, 0][j, \delta]}$, we obtain the following expression for the subdomain normals

$$N_{p[a_i, 0][j, \delta]} = \sigma_{p[a_i, 0][j, \delta]} \cdot A_{p[a_i, 0][j, \delta]}^t \cdot A_{p[a_i, 0][j, \delta]}$$

where $\sigma_{p[a_i, 0][j, \delta]}$ is the approximate variance of the error for the data set $p$ and the subdomain $H_{[a_i, 0][j, \delta]}$.

Estimates of $\sigma_{p[a_i, 0][j, \delta]}$ for the different subdomains and data sets are obtained by variance components analysis of the DWT of the residuals of the data sets. If we note $v_p$ the residuals to the data set $p$ of size $n_p$, we have (Koch 1986; Kusche 2003):

$$v_p^t \cdot W_p^{-1} \cdot v_p = n_p - \text{Trace}[N^{-1} \cdot N_p].$$

This may also be written: $(F_p \cdot v_p)^t \cdot D_p^{-1}$. $(F_p \cdot v_p) = n_p - \text{Trace}[N^{-1} \cdot N_p]$, where $F_p \cdot v_p$ contains the DWT of the residuals. If we select the components of $F_p \cdot v_p$ corresponding to the scale and blocks subdomains, denoted $v_p[a_i, 0][j, \delta]$, by application of a matrix $F_p[a_i, 0][j, \delta]$ of size $n_p \times n_p$ which entries are the entries of $F_p$ for scales and positions roughly corresponding to the domain $H_{[a_i, 0][j, \delta]}$, and 0 otherwise, we have

$$\sigma_{p[a_i, 0][j, \delta]} \approx v_p^t[a_i, 0][j, \delta] \cdot v_p[a_i, 0][j, \delta] / (n_p[a_i, 0][j, \delta] - \text{Trace}[N_{[a_i, 0][j, \delta]}^{-1} \cdot N_{p[a_i, 0][j, \delta]}])$$

where $n_p[a_i, 0][j, \delta] = \text{Trace}[(F_p[a_i, 0][j, \delta])^t \cdot F_p[a_i, 0][j, \delta]]$, and $N_{p[a_i, 0][j, \delta]}$ is the subdomain normal matrix for data set $p$.

Such weighting naturally leads to assess a larger influence to the satellite gravity derived information at the largest wavelet scales, and to the surface gravity data for the smaller scales wavelets. With this choice of parameters, the iterative algorithm will lead to a combined wavelet model that cannot perfectly adjust partially inconsistent data sets. The discrepancies between the wavelet model and the data sets reflect how tightly the model is adjusted to these data sets for the different subdomains. They can be used to define corrector models to the data sets, accounting for their systematic errors. Then, the method is applied again to the corrected data sets to derive an improved wavelet model.
## 3.5 The regularization

The algorithm may be directly applied to the regularized normal system. However, in the case of handling a large number of blocks and wavelet coefficients, the condition number of the normal matrix should be small enough in order for the solution to converge rapidly. This requires to increase the weight of the regularization matrix especially for high-resolution surface data sets, more than what one expects. Consequently, we choose to apply an iterated regularization approach (Engl 1987). In such an approach, the regularized normals $G \cdot x = f$ are solved for the subdomains, but the right hand side is updated by using the non-regularized normals $N$. This is equivalent to a progressive removal of the initial regularization as the iteration continues, the regularization being finally controlled by the number of iterations. This approach is applied for the iterations over the scales subdomains, whereas we applied a full update of the right hand side, using the regularized blocks normals, in the blocks iterations.

The number of iterations over the scales is chosen as follows. We applied a simple approach, using fixed iteration cycles over the scales and stopping the iterations when the calculations have converged to a wavelet model that fits the data within the noise level. Iterating too much destabilizes the solution and degrades the wavelet model quality. Indeed, as we progress in the iterations, the amount of signal to fit decreases, whereas the noise is still there, leading to a decrease of the signal-to-noise ratio. This may destabilize the computations, and it becomes necessary, either to increase the regularization parameter, or to stop the iterations. Although not

done here, it would also be possible to define the stopping criterium more precisely, for instance by using adaptative iteration cycles, where the order of scales is not fixed in advance but depends on intermediate computational results (Wesseling 1991), and applying through the iterations the L-curve or generalized cross-validation methods on the subdomains.

## 3.6 Discussion of the parameters choice

The performance of the method depends on an appropriate choice of the following parameters: subdomains definition, size of the overlap areas, weighting of the extension/restriction operators, approximation of the subdomains normals, condition number, iteration scheme over the scales and stopping criterium for the iterations. Sections 4 and 5 provide a typical parametrization as an example for practical applications. Here we provide some elements to explain these choices.

With the notations of Section 3.2, the error $\epsilon^n$ on the solution $x$ verifies after $n$ iterations

$$ \|\epsilon^n\| \le \|(Id - M \cdot N)^n\| \cdot \|\epsilon^0\|, $$

where $\epsilon^0$ is the initial error. Thus, a reduction of the initial error by a factor $10^{-d}$ is obtained for $n \approx \frac{d}{\log_{10}[\rho(Id - M \cdot N)]}$, where $\rho(Id - M \cdot N)$ is the spectral radius of the matrix. It predicts the worst-case error reduction over the iterations and must be smaller than 1. Higher convergence rates are obtained for small positive values of $\rho$, corresponding to a case where $M$ is a good approximate inverse of $N$, and the convergence rate decreases when $\rho$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Geographic distribution of the synthetic data. Left-hand panel: potential data. Right-hand panel: gravity anomaly data.</figcaption></figure>
approaches 1. Consequently, in order to ensure good convergence, and as mentioned above, the definition of the subdomains must be done in order to follow as closely as possible the structure of the normal system, and the local approximations should not degrade too much the local normals.

Whereas the scale subdomains appear naturally, the *definition of the spatial blocks* may be done in different ways. In this paper, we show a simple parametrization where the blocks are limited by meridians and parallels, but a more flexible spatial splitting can be used. We noticed here that during the iterations, the residuals are located all around the blocks edges, and they are smaller at the centre of the blocks. As a consequence of this heterogeneous distribution, a definition of the blocks that respects the structure of the gravity signal, should improve the convergence rate by avoiding the blocks edges to cross large gravity anomalies. In contrast, we also considered fully entangled spatial domains. We observed in this case a more homogeneous distribution of the iteration residuals, but a worse convergence rate as compared to the geographical blocks.

*Overlap between blocks* introduces a coupling between the local

solves and allows a faster diffusion of the local solutions between subdomains: the larger the overlap, the faster the convergence (Frommer & Szyld 2000). To be efficient, at a given wavelet scale, this coupling is done on an area of width a few times the wavelet scale. In the examples we show, we extend the non-overlapping blocks with an area of width $\delta$ equal to four times the wavelet scale $a$, leading to a total overlap between blocks equal to $\delta a$. This coupling also introduces a redundancy $q$ in the determination of the wavelet coefficients for wavelets centred in the overlap areas, and it has been observed heuristically that the abovementioned choice of weight $w = q^{-1}$ in the definition of the weighted restricted extension operator $\tilde{R}$ is particularly efficient (Cai & Sarkis 1999; Frommer & Szyld 2000). We thus follow this choice in all our tests.

Another parameter of the method is the *level of approximation of the subdomains normals*, that actually allows to filter the data at the spatial resolution of the used spatial mesh. If the level of approximation is too crude, then the iterations will not converge because $M$ provides a too bad approximation of the inverse of $N$. This is the reason why the size of the 3-D mesh used in approximating

<table><caption>**Table 1.** Scales of wavelets in wavelet frames used in local modelling.</caption><thead><tr><th>Level</th><th>Scale (dimensionless)</th><th>Spatial scale (km)</th><th>Number of wavelets</th><th>Area covered</th></tr></thead><tbody><tr><td>6</td><td>0.046875</td><td>300</td><td>380</td><td>25/49°N, 129/153°E</td></tr><tr><td>7</td><td>0.023438</td><td>150</td><td>1406</td><td>25/49°N, 129/153°E</td></tr><tr><td>8</td><td>0.011719</td><td>75</td><td>2401</td><td>29/45°N, 133/149°E</td></tr><tr><td>9</td><td>0.005859</td><td>38</td><td>9604</td><td>29/45°N, 133/149°E</td></tr><tr><td>10</td><td>0.0029297</td><td>20</td><td>38220</td><td>29/45°N, 133/149°E</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Geographic distribution of residuals of the synthetic data. Left-hand panel: potential residuals to degree 120. Right-hand panel: gravity anomalies residuals at 15 km resolution.</figcaption></figure>
the subdomain normals must be sufficiently refined as the iterations progress. In the examples presented in this paper, the last iterations are carried out on the exact subdomains normals. However, when the number of data is orders of magnitude larger than the number of wavelet coefficients, replacing the exact normals with approximated ones, such as $ρ(Id - M·N)$ remains smaller than 1, leads to an approximation of the exact solution at a level of precision that can be kept under control and is fixed by a trade-off between the number of iterations and the data compression rate (Minchev *et al.* 2009).

The order of the iterations over the wavelet scales, can be fixed in advance, as done in this study, or adapted as we progress in the calculation. Multiscale wavelet representations are very close to multilevel approaches such as multigrids, where Schwarz iterations are applied between the different grid subdomains. We can identify the sum of subdomains for all the wavelet scales coarser than $a_l$, namely $\sum_{i=0}^{l} H_{[a_i, 0]}$, with the grid of resolution $a_l$. The iteration schemes between grids of varying resolutions are numerous, and they have been widely studied and compared in the multigrid literature (see for instance Wesseling 1991). Consequently, to design the iteration cycles over the wavelet scales subdomains, we followed the classical iterations schemes between grids. They are represented on Fig. 3 that summarizes and illustrates the applied algorithm.

More generally, there is a trade-off between (i) the condition number of the system $N$, (ii) how closely the restrictions to the local subdomains approximate $N$ and (iii) the number of iterations to be carried out. An improved condition number allows to use smaller blocks and/or to reduce the size of the overlap areas while keeping fast convergency. On the other hand, if the condition number is too bad, the iterations may never converge. A classical criterium for stopping the iterations is to estimate the reduction factor, which is equal, at the $k$th iterate for a given scale $a_i$, to the ratio $\frac{\|x_{[a_i, 0]}^k - x_{[a_i, 0]}^{k-1}\|}{\|x_{[a_i, 0]}^{k-1} - x_{[a_i, 0]}^{k-2}\|}$, and stop iterating when it gets smaller than a threshold. We apply this criterium to find when to stop iterations between the spatial blocks.

# 4 VALIDATION WITH SYNTHETIC DATA

## 4.1 White noise case

We first validate this approach with synthetic data over Japan. Two sets of synthetic data are prepared for that purpose.

The first one is a set of gravity potential values as given by EGM2008 (Pavlis *et al.* 2008), up to complete degree 120 (~170 km resolution): 5448 synthetic data on the ground level regularly distributed in the area between latitudes 27°N to 47°N and longitudes 131°E to 151°E. White noise (rms = 1 m² s⁻²) is added. This noise level is a little lower than that of the satellite-only geopotential models as given by their error spectra up to degree/order 120.

The second data set is a 3 × 3 min grid of 103 041 free-air gravity anomalies on the ground level, also computed from EGM2008, over the area between latitudes 29°N to 45°N, and longitudes 133°E to 149°E. The resolution of this grid (5.5 km) corresponds to typical resolutions of surface gravity data and, accordingly, we can validate the approach under realistic situations. However, in order to simplify our test, the grid gravity anomalies are smoothed by truncating EGM2008 to degree 1000, and by applying damping starting at degree 730 with a half amplitude at degree 835. This corresponds to a spatial resolution of about 24 km. Finally, white noises are added to the data set. At the resolution of the computed wavelet

model (about 15 km), the amplitude of the noise is 0.35 mGal. This value is chosen within the range of terrestrial gravity data precision. As mentioned above, the low-frequency components of EGM2008 have been removed from both data sets. Fig. 4 shows the geographical distribution of the synthetic data.

The wavelet frame used is composed of five levels of scales as given in Table 1. Scales of the wavelets are about 300, 150, 75, 38 and 20 km, corresponding to the depths of the equivalent multipolar sources below a sphere of radius 6370 km, the average Earth radius. Because of the ellipticity, the depths with respect to the Earth's ellipsoid varies by ±3 km in the studied area, modifying slightly the wavelet scales on the ground level. Although the finest scale of the wavelets is 20 km, the wavelets can capture signal down to about 15 km resolution, because of the smooth decay of their spectrum. Finally, let us underline that, in the preparation of the synthetic potential data, we truncate EGM2008 at degree 120. Accordingly, we should also truncate the spherical harmonics expansion of the wavelets at the same degree for consistency. On the other hand, the anomaly data are treated as real observations independent of the potential data (though we apply a damping at the highest degrees to reduce the aliasing effects in the analysis). Therefore, we do not need to truncate the harmonics expansion of the wavelets in the construction of the gravity anomalies observation equations.

On levels 8, 9 and 10, the area is divided into four, 16 and 36 blocks, respectively, and the computation of the wavelet coefficients

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Histograms of residuals between the wavelet model and the synthetic data. Top panel: potential residuals at degree 120. Bottom panel: gravity anomalies residuals.</figcaption></figure>
are made in each block. Regarding the iteration scheme, we first test one FMG-cycle. At each step of the cycle the data are reduced in cells at a resolution corresponding to the finest wavelet scale that has been already computed. Next, we apply about 40 V-cycle without any data reduction. On the finest scale, the total number of cells exceeds slightly the number of data themselves and, consequently, the reduction of the data in the cells becomes meaningless. On each scale, the number of block iterations ranges from a few tens to a few hundreds when there is more than one block. The number of iterations for the blocks is controlled by the stopping criterium given in Section 3.6, and we stopped iterating over the scales when the misfit between the wavelet model and the synthetic data reached the synthetic noise level. The weights to the data are given homogeneously in space according to the noise level. In the spectral domain, larger weights are assigned to the geopotential data at the largest scales and to the free-air data at the smaller scales.

The residuals to the respective data sets are shown geographically in Fig. 5 and their histograms are given in Fig. 6. The rms of residuals as in Fig. 6 are well restored to the applied noise levels and indicate a good performance in iteration: 1.02 m² s⁻² with an average of 0.03 m² s⁻² for the potential data, and 0.36 mGal with an average of 0.01 mGal for the gravity anomaly data. The spatial pattern of the residuals seems to be that of white noise in both the potential and gravity anomalies data.

Only slightly noticeable is the feature with a larger amplitude around east Hokkaido, where EGM2008 (and the actual gravity field) contains strong signals at high frequencies. During the iteration process we learn that the residuals there show a small, but systematic trend of some 0.05 mGal in amplitude at earlier stages

of iteration, but which has been gradually disappeared as iteration goes.

## 4.2 Coloured noise case

As a second validation of the method with synthetic data, we assume the existence of realistic systematic errors as coloured noise in the anomaly data, with a localized distribution in space. Such types of errors may exist near the coasts in altimetry-derived gravity anomalies or in ship-borne gravity surveys. Then, we additionally append such errors as shown in Fig. 7 to the free-air gravity anomaly data, but we do not modify the potential data, both obtained in the preceding section. In particular, we do not use the full variance/covariance matrix information for the EGM2008 model. This choice was made in order for the test to follow the settings of the application presented in the next section, for which the variance/covariance matrix of the geopotential model is not available. Finally, the same analysis as that of the previous section is carried out to these data.

After completion of the one FMG-cycle and two V-cycle iterations, the residuals to the respective data are shown geographically in Fig. 8. Same as for the white noise test, larger weights are assigned to the geopotential data at the largest scales and to the free-air data at the smaller scales in the spectral domain. In addition, to account for the spatial variability of the systematic errors, the 300 km scale subdomain has been split into four blocks of equal size (two blocks in latitude and two blocks in longitude) and data weights at this scale are computed for each block. Taking as a reference the northeastern block, the gravity anomalies are downweighted by a factor of 1.5

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Coloured noise added to the gravity anomaly data (left-hand panel). Spherical harmonics expansion of the noise up to degree/order 120 (right-hand panel)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Geographical distribution of residuals to the synthetic data. Left-hand panel: potential residuals up to degree 120. Right-hand panel: gravity anomalies at 15 km resolution.</figcaption></figure>

in the northwestern block, and by a factor of 3 in the two southern blocks, where the data noise is the largest.

The result of the computation shows notorious patterns correspondingly to the location of the coloured noise in the residuals to the gravity anomaly data, whereas the residuals to the potential data are rather uncorrelated. Small artefacts are observed in the vicinity of the main patterns in the map of residuals to the gravity anomaly data. This is due to the fact that the potential data allow to highlight only a part of the total error: its spherical harmonics components up to degree/order 120, leading to Gibbs oscillations (see Fig. 6). Part of these artefacts may also arise from the approximations involved in the estimation of constant noise variances per subdomain.

Applying a low-pass filter to the map of residuals to the gravity anomalies in order to keep only the dominant patterns, we build a corrector model for the anomaly data, represented on Fig. 9. By using the corrected anomaly data thus obtained, we start again the computations and derive an improved wavelet model. As a corrector has been applied, in order to gain computational time, we do not re-estimate different weights for the spatial blocks, but assign spatially uniform weights for the different scales subdomains. Fig. 10 displays the geographical distribution of the residuals to the two data sets. We observe no clear pattern in the residuals to the potential data, but a pattern of white noise. The rms of the residuals, $1.03 \text{ m}^2 \text{ s}^{-2}$, is compatible to the level of the assumed white noise. The residuals to the gravity anomaly data become also closer to the white noise structure, with rms 0.39 mGal and slight correlated patterns left. Even if the noise recovery is not perfect, we conclude that our approach leads to a fairly good model given the synthetic

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Corrector model to the gravity anomaly data set.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Geographical distribution of residuals to the synthetic data after applying downweights to the anomaly data. Left-hand panel: potential residuals up to degree 120. Right-hand panel: gravity anomalies at 15 km resolution.</figcaption></figure>

data characteristics, and that the method developed is valid for local gravity field modelling from different sets of data with different error characteristics.

Note that increase in weights at intermediate scales of the potential data would not improve the resulting model. At that scales the anomaly data contain reliable signals, but would be assigned to smaller weights, leading to worsen the anomaly data fitted to the pattern, if any, of white noise in the potential data.

Finally, it is important to recognize the fact that the pattern of the gravity anomaly residuals does not perfectly reproduce the coloured noise by the applied wavelet method with constant weights per sub-domain at one hand, and to keep in mind the difficulty in determining proper relative weights to respective data on the other hand. In addition, we should remember the fact that it is not possible to fully rely on the potential data at the large wavelet scales because they also include their own imperfection.

## 5 APPLICATION OVER JAPAN

We finally apply the method to real data over Japan: a high-resolution local gravity model from Kuroishi & Keller (2005) and a spherical-harmonics model of the geopotential, EIGEN-GL04S, complete to degree 150 (Biancale *et al.* 2005).

The geopotential model is developed only from GRACE and LAGEOS measurements. Its cumulative error at degree 120, in terms of rms, is estimated to about 0.8 m² s⁻², corresponding to

about 8 cm in geoid height error. We calculate 5448 potential values from the model up to degree 120, regularly spaced on the ground level in the same area as that of the synthetic validation.

The local gravity model is based on a combination of land and marine gravity data and satellite altimetry derived gravity anomalies from KMS2002 (Andersen & Knudsen 1998). We decimate it on a grid of 3 × 3 min and take 103 041 Faye anomalies on the ground level. The geographical distribution of Faye anomaly is shown in Fig. 11. Highest frequency undulations below 10 km of wavelength have been damped by a moving-average filter before applying the wavelet analysis. We subtract low-degree components of EIGEN-GL04S from both data sets and apply to those data sets the method with the same parameter setting as those of the validation tests. The iteration scheme here is one FMG-cycle, followed by two V-cycles, and itis repeated three times on progressively corrected data sets (see below). For each iteration scheme, we do not iterate the variance components estimation of weights, because, in contrast to the synthetic tests, iterations produce weights associated to an overfitting of the potential data and underfitting of the surface gravity data. Replacing the real surface gravity with the synthetic gravity anomalies used in Section 4.2, for which we are sure that there is a good adequacy between the data and the wavelet model, leads to the same conclusion. A possible explanation might be related to the fact that, when the noise belongs to the same three dimensional wavelet space as the signal, it is more difficult to distinguish it than when it is orthogonal to the signal. The first variance components computation is considered to roughly provide
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Local gravity anomaly model for Japan from Kuroishi and Keller (2005).</figcaption></figure>

<table><caption>**Table 2.** Relative weighting of the data sets used in gravity field modelling over Japan (first computation, on non-corrected data sets).</caption><thead><tr><th>Level</th><th>EIGEN-GL04S data weights<br>(m<sup>2</sup> s<sup>-2</sup>)</th><th>Surface gravity weights<br>(mGal)</th></tr></thead><tbody><tr><td>6</td><td>2</td><td>30</td></tr><tr><td>7</td><td>1</td><td>3</td></tr><tr><td>8</td><td>0.7</td><td>2.6</td></tr><tr><td>9</td><td>0 (no data used)</td><td>1.3</td></tr><tr><td>10</td><td>0 (no data used)</td><td>1.5</td></tr></tbody></table>

a reasonable weighting scheme. At the very end of the computation, we stopped the last V-cycle before it ends, to avoid a destabilization of the calculation.

After a first run, the residuals to the respective data are shown geographically in Fig. 12. The corresponding data weighting is given in Table 2. The downweighting of the surface gravity is important at the largest scales, mostly constrained by the EIGEN-GL04S data. The rms of the residuals are 0.9 m² s⁻² with a bias of −0.1 m² s⁻² for the potential data, and 1.20 mGal with a bias of 0.70 mGal at 15 km resolution for the anomaly data. These rms are reasonable in consideration of the data precision. Systematic residuals at large scales are remarkable particularly south of the Japanese islands and obvious along the coastal areas of the main island, Honshu as well.

These features at the ocean in the residuals are likely to reflect the systematic errors in the anomaly data controlled by the altimetric model. Kuroishi (2009) shows similar results by comparing GGM02C/EGM96 with the local gravity model and develops a

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Geographical distribution of residuals. Left-hand panel: potential residuals to degree 120. Right-hand panel: gravity anomalies at 15 km resolution.</figcaption></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Correctors to the local gravity anomaly model. Left-hand panel: corrector estimated after a first series of iteration, Right-hand panel: the final corrector used.</figcaption></figure>

<table><caption>**Table 3** Relative weighting of the data sets used in gravity field modelling over Japan (second computation, with corrected gravity anomaly grid using the corrector from Fig. 12, left-hand panel).</caption><thead><tr><th>Level</th><th>EIGEN-GL04S data weights (m<sup>2</sup> s<sup>-2</sup>)</th><th>Surface gravity weights (mGal)</th></tr></thead><tbody><tr><td>6</td><td>2</td><td>10</td></tr><tr><td>7</td><td>1</td><td>3</td></tr><tr><td>8</td><td>0.5</td><td>2.3</td></tr><tr><td>9</td><td>0 (no data used)</td><td>1.3</td></tr><tr><td>10</td><td>0 (no data used)</td><td>1.5</td></tr></tbody></table>

highly improved gravimetric geoid model for Japan, JGEOID2008, after removal of such errors from the local gravity model. This demonstrates how the uniformity of accuracy of the GRACE-derived static gravity model contributes much to the detection of areas of degraded quality in the local gravity data.

Based on the discussion, we try to correct the anomaly data for an improved combination. First, we exercise a low-pass filter to the anomaly residuals at the resolution of EIGEN-GL04S. The left panel in Fig. 13 shows the model corrector obtained, which is subtracted from the anomaly data. Then we apply again the developed method to the corrected data sets. The corresponding weights are described in Table 3. This second computation allows to refine the corrector model for the surface gravity, leading to a new corrector model represented on the right panel of Fig. 13. A last computation is carried out on the corrected data sets, with relative weighting of 1.

The residuals to the geopotential data and to the corrected anomaly data are represented geographically in Fig. 14. The rms of residuals are 0.80 m$^2$ s$^{-2}$ for the potential data, and 0.50 mGal at 15 km resolution for the corrected anomaly data. We find that no significant bias remains in both residuals and these rms are consistent with the estimated levels of data noise.

In the plot of the residuals to the corrected anomaly data, on the right-hand panel in Fig. 14, features only on quite small scales are dominant. This indicates that the resolution of the combined wavelet model is a little coarser than that of the anomaly data. In addition, we observe some edge effects, especially in the northern and southern boundaries in the case of the anomaly data. The same tendency is also visible in the case of the potential data. This shows that the inversion slightly lacks stability in these areas.

Finally, we note that constructing the corrector model by simple low-pass filtering of the residuals of the surface data with respect to the combined wavelet model is a very simple approach. As the errors in the residuals patterns are not only coloured in spectrum, but also localized in space, one may consider for instance using wavelet analyses to construct the corrector model from the residuals maps, in a similar fashion as done by Kuroishi & Keller (2005). In the future, we will thus work on improvements of the construction of the corrector models. Then, another interesting element is to take into account the relative weights that finally come out of the last computation. Indeed, at this stage, all or almost all longer wavelength errors have been removed, thanks to progressively refined corrector models, from the surface gravity data. Consequently, the weights obtained should provide an estimate of the signal quality in the
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Geographic distribution of residuals in the final combination. Left-hand panel: potential residuals to degree 120. Right-hand panel: residuals of corrected gravity anomalies at 15 km resolution.</figcaption></figure>

corrected data. This will be particularly interesting when combin-
ing the future GOCE satellite gravity data with surface gravity data
sets—all the more than in the case of a combination of real measure-
ments, such weights can be better assessed than for a combination
with a harmonic geopotential model, as explained above.

6 CONCLUSION

We have developed an iterative method to combine various kinds
of gravity data into a wavelet model of the geopotential, taking into
account the error characteristics of the different data sets in different
wavebands and locations. This method was validated with synthetic
data and then applied to real data over Japan: local high-resolution
gravity anomaly data and a GRACE-derived global model, EIGEN-
GL04S. We obtained a hybrid spherical harmonics/wavelet model
of the geopotential over Japan at about 15 km resolution and the
residuals to the respective data underlined biases on medium scales
between the two data sets, which suspected origin is errors in the
anomaly data. We then corrected the anomaly data by subtracting
the evidenced biases and repeated the method again to the corrected
data sets, resulting in an improved hybrid model of the gravity field
over Japan. The method may be improved in the future, for instance
by using a more flexible definition of the subdomains, depending
on the variability of the gravity signal and of the noise patterns.
Another improvement is needed to reach very high spatial resolu-
tions. The wavelet functions at a given scale, as they appear on the
topographic level, are indeed smoother than on the mean sphere due

to the upward continuation between the sphere and the topographic
surface. If the difference is small at low resolution, is becomes
quite large when reaching topographic scales and the design of the
wavelet frame must be refined to model gravity anomalies given on
the topography at resolutions below 15 km.

The method has broad perspectives of applications. First, it should be used to combine directly the satellite gravity products from the GRACE and GOCE missions with surface and airborne gravity data. This should allow to improve the regional gravity models derived, and contribute to the validation of the satellite gravity data using ground measurements. Moreover, because we can handle any functional of the gravity potential, the method should be very useful to also combine locally or globally the GRACE and GOCE measurements themselves, together with possible future data from possible future satellite gravity missions. Finally, the combined gravity models thus obtained can be easily analysed using wavelets for geodynamic purpose, as shown by Panet et al. (2006). Performing their joint analysis together with other geophysical data such as magnetic and seismology data, will allow to improve our understanding of the inner structure of our planet and of the geodynamic processes at stake, from local to regional and global scales.

ACKNOWLEDGMENTS

We thank the CNES space agency for financial support of this work
through the TOSCA committee. Part of this work was done during
the CNES post-doctoral fellowship of Isabelle Panet. The comments
by Jürgen Kusche and anonymous reviewers were of great value
for improving our manuscript. All maps were plotted using the GMT software (Wessel & Smith 1995). This is IPGP contribution number 3095.

## REFERENCES

Andersen, O.B. & Knudsen, P., 1998. Global marine gravity field from the ERS-1 and Geosat geodetic mission altimetry, *J. geophys. Res.*, **103**, 8129–8137.

Biancale, R., Lemoine, J-M., Balmino, G., Loyer, S., Bruisma, S., Perosanz, F., Marty, J-C. & Gegout, P., 2005. Three years of decadal geoid variations from GRACE and LAGEOS data, CD-Rom, CNES/GRGS product.

Cai, X-C. & Sarkis, M., 1999. A restricted additive Schwarz preconditioner for general sparse linear systems, *SIAM J. Scient. Comput.*, **21**, 792–797.

Chambodut, A., Panet, I., Mandea, M., Diament, M., Jamet, O. & Holschneider, M., 2005. Wavelet frames: an alternative to the spherical harmonics representation of potential fields, *Geophys. J. Int.*, **168**, 875–899.

Chan, T. & Mathew, T., 1994. Domain decomposition algorithms, *Acta Numer.*, 61–143.

Eicker, A., Mayer-Gürr, T. & Ilk, K-H., 2006. A global CHAMP gravity field by merging regional refinement patches, in *Proceedings of the Joint CHAMP GRACE Science Team Meeting, Potsdam*, 2004, electronic publication http://www-app2.gfz-potsdam.de/pb1/JCG/jcg_index.html.

Engl, H.W., 1987. On the choice of the regularization parameter for iterated Tikhonov regularization of ill-posed problems, *J. Approx. Theory*, **49**, 55–63.

Freeden, W., Gervens, T. & M., Schreiner, M., 1998. *Constructive Approximation on the Sphere (With Applications to Geomathematics)*, Oxford Science Publication, Clarendon Press, Oxford.

Frommer, A. & Szyld, D., 2000. An algebraic convergence theory for restricted additive Schwarz methods using weighted max norms, Research Report BUGHW-SC 00/3, Bergische Universität GH Wuppertal.

Härdle, W. & Simar, L., 2003. *Applied Multivariate Statistical Analysis*, Springer-Verlag, Berlin Heidelberg, New York.

Holschneider, M., 1995. *Wavelets: An Analysis Tool*, Oxford Sciences Publications, Oxford.

Holschneider, M., Chambodut, A. & Mandea, M., 2003. From global to regional analysis of the magnetic field on the sphere using wavelet frames, *Phys. Earth planet. Inter.*, **135**, 107–124.

Kaula, W.M., 1966. *Theory of Satellite Geodesy*, Blaisdell, Waltham.

Keller, W., 2001. A wavelet approach for the construction of multi-grid solvers for large linear systems, in *Proceedings of the IAG Sc. Ass., IAG Symposium Vistas for Geodesy in the New Millenium*, Vol. 125, pp. 265–270, ed. Sanso, Springer, Berlin.

Klees, R. & Wittwer, T., 2007. Local gravity field modelling with multipole wavelets, in *Proceedings of the IAG Sc. Ass., IAG Symposium Dynamic Planet*, Vol. 130, pp. 303–308, ed. Tregoning, Rizos, Springer, Berlin.

Koch, K-R., 1986. Maximum likelihood estimate of variance components, *Bull. Geod.*, **60**, 329–338.

Kuroishi, Y., 2009. Improved geoid model determination for Japan from GRACE and a regional gravity field model, *Earth, Planets, Space*, **61**(7), 807–813.

Kuroishi, Y. & Keller, W., 2005. Wavelet improvement of gravity field-geoid modeling for Japan, *J. geophys. Res.*, **110**, B03402, doi:10.1029/2004JB003371.

Kusche, J., Ilk, K.H., Rudolph, S. & Thalhammer, M., 1998. Application of spherical wavelets for regional gravity field recovery—a comparative study, in *Proceedings of the IAG Sc. Ass., IAG Symposium Geodesy on the Move*, Vol. 119, pp. 213–218, eds Forsberg, R., Feissel, M. & Dietrich, R., Springer, Berlin.

Kusche, J., 2001. Implementation of multigrid solvers for satellite gravity anomaly recovery, *J. Geodyn.*, **74**, 773–782.

Kusche, J., 2003. A Monte-Carlo technique for weight estimation in satellite geodesy, *J. Geodyn.*, **76**(11–12), 641–652.

Lerch, F.H., Marsh, J.G., Klosko, S.M., Patel, G.B., Chin, D.S., Pavlis, E.C. & Wagner, C.A. 1991. An improved error assessment for the GEM-T1 gravitational model, *J. geophys. Res.*, **96**, 20 023–20 040.

Mallat, S., 1999. *A Wavelet Tour of Signal Processing*, 2nd edn, Academic Press, San Diego.

Minchev, B., Chambodut, A., Holschneider, M., Panet, I., Scholl, E., Mandea, M. & Ramillien, G., 2009. Local multipolar expansions for potential fields modeling, *Earth, Planets, Space*, **61**(10), 1127–1141.

Panet, I., Jamet, O., Diament, M. & Chambodut, A., 2004. Modelling the Earth's gravity field using wavelet frames, in *Proceedings in the IAG Sc. Ass., IAG Symposium GGSM 2004*, Vol. 129, pp. 48–53, Porto, Portugal, 30 Aug.-3 Sept. 2004, eds Jekeli, C., Bastos, L. & Fernandes, J., Springer, Berlin.

Panet, I., Chambodut, A., Diament, M., Holschneider, M. & Jamet, O., 2006. New insights on intraplate volcanism in French Polynesia from wavelet analysis of GRACE, CHAMP and sea-surface data, *J. geophys. Res.*, **111**, B09403, doi:10.1029/2005JB004141.

Pavlis, N., Holmes, S.A., Kenyon, S.C. & Factor, J.K., 2008. An Earth gravitational model to degree 2160: EGM2008, *Presented at the 2008 General Assembly of the European Geosciences Union*, Vienna, April 13–18 2008.

Schaffrin, B. & Iz, H.B., 2001. Integrating heterogeneous data sets with partial inconsistencies, in *Proceedings of the IAG Sc. Ass., IAG Symposium Gravity, Geoid and Geodynamics 2000*, Vol. 123, pp. 49–54, Banff, Alberta, Canada, 31 July-4 Aug. 2000, ed. Sideris, M., Springer.

Schmidt, M., Kusche, J., Shum, C.K., Han, S.-C., Fabert, O. & van Loon, J. 2005. Multiresolution representation of regional gravity data, in *Proceedings of the IAG Sc. Ass., IAG Symposium GGSM 2004*, Porto, Portugal, 30 Aug.-3 Sept. 2004, pp. 167–172, eds Jekeli, C., Bastos, L. & Fernandes, J., Springer, **129**.

Schmidt, M., Fengler, M., Mayer-Gürr, T., Eicker, A., Kusche, J., Sanchez, L. & Han, S.-C., 2007. Regional gravity modeling in terms of spherical bases functions, *J. Geodyn.*, **81**(1), 17–38.

Tenzer, R. & Klees, R., 2008. The choice of the spherical radial basis functions in local gravity field modeling, *Stud. Geophys. Geod.*, **52**(3), 287–304.

van Loon, J.P. & Kusche, J., 2005. Towards an optimal combination of satellite data and prior information, in *Proceedings IAG Sc. Ass., IAG Symposia*, Vol. 130, Cairns, Australia, 22–26 August 2005, pp. 345–353, eds Tregoning, P. & Rizos, C., Springer.

Wesseling, P., 1991. *An Introduction to Multigrid Methods*, ed., Wiley, Philadelphia, 294pp.

Wessel, P. & Smith, W.H.F., 1995. New version of the Generic Mapping Tool released, EOS, Trans. AGU, **76**, 329

Xu, J., 1992. Iterative methods by space decomposition and subspace correction, *SIAM Rev.*, **34**(4), 581–613.