# Local multi-polar expansions in potential field modeling

B. Minchev, A. Chambodut, M. Holschneider, I. Panet, E. Schöll, M. Mandea, G. Ramillien

► To cite this version:

B. Minchev, A. Chambodut, M. Holschneider, I. Panet, E. Schöll, et al., Local multi-polar expansions in potential field modeling. Earth Planets and Space, 2009, 61, pp.1127-1141. 10.1186/BF03352965. insu-03604890

HAL Id: insu-03604890

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

Submitted on 10 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
Local multi-polar expansions in potential field modeling

B. Minchev<fnref n="1" />*, A. Chambodut<fnref n="2" />, M. Holschneider<fnref n="1" />, I. Panet<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, E. Schöll<fnref n="6" />, M. Mandea<fnref n="7" />†, and G. Ramillien<fnref n="8" />‡

(Received September 21, 2008; Revised March 7, 2009; Accepted June 28, 2009; Online published November 30, 2009)

The satellite era brings new challenges in the development and the implementation of potential field models. Major aspects are, therefore, the exploitation of existing space- and ground-based gravity and magnetic data for the long-term. Moreover, a continuous and near real-time global monitoring of the Earth system, allows for a consistent integration and assimilation of these data into complex models of the Earth’s gravity and magnetic fields, which have to consider the constantly increasing amount of available data. In this paper we propose how to speed up the computation of the normal equation in potential field modeling by using local multi-polar approximations of the modeling functions. The basic idea is to take advantage of the rather smooth behavior of the internal fields at the satellite altitude and to replace the full available gravity or magnetic data by a collection of local moments. We also investigate what are the optimal values for the free parameters of our method. Results from numerical experiments with spherical harmonic models based on both scalar gravity potential and magnetic vector data are presented and discussed. The new developed method clearly shows that very large datasets can be used in potential field modeling in a fast and more economic manner.

**Key words:** Potential fields (gravity, geomagnetism), inverse problem, spherical harmonics, satellite data, size reduction.

# 1. Introduction

Our knowledge of the Earth’s system is far from being complete, and therefore gaining deeper insight into global processes and their interactions is one of the most urgent challenges in geo-sciences. The development of accurate and reliable global models must go hand in hand with the continuous global monitoring of the gravity and magnetic field, and their changes occurring at regional to global spatial scales, and over timescales ranging from extremely short duration events to long duration ones. The long time series of geodetic, geodynamic, and magnetic parameters obtained from such an observing system provide crucial information about the Earth system, i.e. about the shape and deformation of the Earth, its internal structure and material composition, the behavior of the hydrosphere, cryosphere and atmosphere, the Earth’s variable rotation, the static and time-variable gravity and magnetic fields. Moreover, nowadays, a dense global Earth’s observing system, aiming to

high temporal and spatial resolution and accuracy, has to integrate and combine the required large variety of complementary ground- and space-based observational data. The new satellite epoch for gravity and magnetic observation has brought us a new dimension in the modeling process—the amount of data to be considered.

Indeed, the launch of the Danish satellite Ørsted in February 1999 has initiated a new era for geopotential field research, known as *International decade of geopotential fields*. Ørsted satellite was followed by other gravity and/or magnetic missions, as CHAMP, SAC-C and GRACE, launched in July 2000, November 2000 and March 2002, respectively. These satellites have provided a large amount of gravity and/or magnetic data.

In a first step of processing, selection criteria are applied to remove measurements with insufficient accuracy. In gravity field modeling, data are flagged according to their quality, which depends on the phase breaks of the measurement system, the filling of data gaps, the quality of applied corrections (e.g. clock corrections), etc. Only the most valid data are selected, with a high enough signal-to-noise ratio.

Considering the magnetic field modeling, identical procedures are used. Several selection criteria (geomagnetic activity indices, local time—see Chambodut et al. (2003) for details) are usually applied to the initial datasets.

Although ignoring observations with insufficient accuracy can change the number of measurements used in the

<footnote n="1">¹Department of Applied and Industrial Mathematics, University of Potsdam, Potsdam, Germany</footnote>

<footnote n="2">²Institut de Physique du Globe de Strasbourg (UMR 7516 CNRS, Université de Strasbourg/EOST), Strasbourg, France</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, CNRS, Paris, France</footnote>

<footnote n="4">⁴Space Geodesy Research Division, Geographical Survey Institute, Tsukuba, Ibaraki, Japan</footnote>

<footnote n="5">⁵Institut Géographique National—Laboratoire LAREG—Marne-la-Vallée, France</footnote>

<footnote n="6">⁶Institute of Theoretical Physics, Berlin University of Technology, Germany</footnote>

<footnote n="7">⁷Helmholtz-Zentrum Potsdam, Deutsches GeoForschungsZentrum—GFZ, Potsdam, Germany</footnote>

<footnote n="8">⁸Laboratoire d'Etudes en Géophysique et Océanographie Spatiales, Observatoire 048 Midi-Pyrénées, Toulouse, France</footnote>

<footnote>*Now at Global Life, Allianz SE, 70178 Stuttgart, Germany.
†Now at Université Paris Diderot—Institut de Physique du Globe de Paris Géophysique spatiale et planétaire, Paris, France.
‡Now at Géodésie Spatiale, CNRS DTP UMR 5562, OMP, Toulouse, France.</footnote>
model, after this selection step, the overall size of the datasets is still very large. In addition, not only the satellite data, but also surface data should be considered. High resolution field modeling, at global or regional scales, indeed requires to combine the ‘medium’ resolution satellite measurements with very dense surface datasets comprising ground-based, altimetric, airborne and marine measurements. This increases considerably the size of the datasets from which the potential fields models are computed. For such high resolution modeling, and for global modeling as well, localized basis functions like spherical wavelets or radial basis functions may be used (see for instance Holschneider et al., 2003; Chambodut et al., 2005; Panet et al., 2006; Lesur et al., 2007) and the availability of fast algorithms is crucial.

In this context, modeling the potential fields of the Earth becomes a challenging task with higher and higher demands of computing memory and time. The number of data to be modeled is huge and in order to be able to fully benefit from the wealth of information contained in those measurement flows, we should aim at treating them as they are, as much as we can.

To cope with the large number of satellite measurements, and in particular with the future 100 millions and 300 millions of measurements, respectively from the forthcoming 20 months-GOCE (gravity) and 4 years-Swarm (magnetic) missions, methods have been developed to avoid the storage of large normal matrices involving a large number of spherical harmonics. Those methods are based on the iterative approach (see for instance Kusche, 2000; Schuh, 2000; Keller, 2001; Nakajima, 2001; Reubelt et al., 2004). They make use of a matrix that is representative of the normal system, but easier to compute. Such matrix allows to compute an approximate solution of the normal system, that is iteratively refined. It can also be used as a pre-conditioner of the normal system, allowing to speed-up convergence rates of iterative solvers (conjugate gradient, ...). For instance, a block-diagonal approximation of the normal system may be obtained when processing gravity data in a time-wise or in a space-wise manner. In the time-wise approach, the measurements are treated as time-series along the orbit and the gravity potential is represented using inclination functions related to the spherical harmonics. In this case, the normal matrix exhibits a block-diagonal structure for a circular, repeat orbit with constant inclination and no data gaps (Colombo, 1984). The approximation of the real orbit with a circular one thus leads to an approximative block-diagonal solver that may be used within an iterative method. This allows important computational savings (see for instance the paper by Klees et al. (2000)). However, by nature such methods do not lead to a reduction of the size of the underlying problem. Moreover, they cannot be applied to localized basis functions.

In the space-wise approach, the measurement positions are directly treated as points in the 3D space, and the modeling functions (spherical harmonics, wavelets, ...) are expressed as functions of positions. Then, a widely used approach to reduce the size of the datasets to a numerically reasonable size is to project them onto a grid. This preprocessing may be applied to both satellite and surface data.

The measurements are often down-sampled, and then interpolated on a regular spherical grid at average satellite or surface altitudes. Note that such a regular gridding can also be used to build a pre-conditioner of the normal system in an iterative resolution framework (Schuh, 2000), since it leads to a block-diagonal approximation of the normals. The data gridding is implemented locally, cell by cell, from the original data, using least-squares collocation and a-priori covariance functions (Moritz, 1980; Langel and Hinze, 1998). Both, area-mean or point values can be obtained. A spherical harmonics model then can be computed from this gridded dataset using quadrature formula for the spherical harmonics on the grid or a least-squares adjustment or collocation (Colombo, 1981; Sanso and Tscherning, 2003). Migliaccio et al. (2004, 2006), applied this method to recover the gravity field from synthetic GOCE data. Another example is the computation of the EGM96 model, complete up to degree and order 360, and constrained with satellite and surface measurements. Lemoine et al. (1998), used 30' mean surface gravity anomalies instead of the original surface measurements: more than 30 millions of terrestrial point gravity anomalies over the continents, in addition to the dense satellite altimetric coverage of the oceans! Recently, under the World Digital Magnetic Anomaly Map Project, an international effort has been made to integrate all available near-surface and satellite magnetic anomaly data into a global map database. In Maus et al. (2007), the authors use a least-square collocation method to combine different surveys of near-surface data and then merged the obtained grid with marine and aeromagnetic line-leveled tracks. For both gravity and magnetic applications, the main drawback of the gridding approach, is that the local interpolations and reduction to a sphere modify the stochastic properties of the data in a complex way and lead to a loss of information. Consequently, it is not clear whether the subsequent possible inversion will lead to the optimal estimator in the least-squares sense. This is the reason why, in our approach, we stick closely to the original data. We propose to approximate the modeling functions used in the representation of the potential field, but not the observation data itself. Base on this approximations we construct an iterative method of stationary type, which under curtain conditions (see Section 3), always converges to the solution of the underlying inverse problem. Thus, there is no risk of changing the stochastic properties of the original data set. Similar technique was also used by Ditmar et al. (2003), where the authors have developed an interpolation methodology implemented within an iterative solver, which allows to approximate the normal system in such a way that the interpolation error is kept under control. The exact values of the spherical harmonics are computed on the vertices of a 3D spherical mesh (this step can be very fast for a convenient mesh geometry), and the values of the spherical harmonics at the measurement points are approximated by interpolation using the closest vertices. Another example is Han (2004), who make use of Taylor series approximations of the Legendre functions in latitudinal bands to speed-up the computation of the normal system in a conjugate gradient approach.

In this paper, we propose a new method for modeling
<figure><img src="image_2.png" /><figcaption>Fig. 1. Space discretization based on the icosahedron. Projection of the icosahedron vertexes onto a sphere (**left**); Shell discretization at satellite altitude (**right**).</figcaption></figure>

large amounts of data, which takes advantage of the rather smooth behavior of the potential fields. Very often, the measurements strongly oversample the smooth modeling functions, which locally behave as low degree polynomials. We design our method based on two main requirements: (i) the data should be pre-processed in a compact and reasonable way; (ii) the inversion scheme should allow fast and local computations.

In this context, we develop local multipole approximations of the modeling functions at satellite altitude. We first design a 3D tiling of the measurement space. Next, we replaced each model function by a power series of local multipolar expansion, around a defined point inside each volume of the discretized space. In this way, we obtain a piecewise decomposition of the potential field at the altitude of the measurements. The coefficients of each local expansion are used to precondition the underlying inverse problem, which is then solved by an iterative method. The speed of convergence and the precision of the inversion process are controlled by a set of free parameters. Changing the relation between the accuracy of the approximation of the modeling functions and the size of the space discretization provides the freedom to compute a fast, quick-look solution or a very high quality model. The specificity of our approach is that these piecewise approximations of the modeling functions are equivalent to a reduction of the original dataset to a set of local moments in the unity-volumes.

The paper also presents two applications, one for the gravity field modeling and one for the magnetic field modeling. It is shown that our proposed method makes it possible: (i) to reduce the amount of data for computational and memory reason; (ii) to get precomputed datasets; (iii) to make available standardized data for databases. All these improvements are obtained without losing the quality of the final obtained model.

## 2. Space Discretization

A starting point in our considerations is tiling the dataset into smaller subsets, such that inside each one, the data

can be replaced by a suitably agglomerated derived quantities. Our space discretization is strongly related to the data distribution. We use an adaptive algorithm, which decomposes the space-domain, i.e. the shell at the satellite altitude, where the satellite revolves, into a sum of geometrical 3D-bodies—curvilinear triangular prisms (see Fig. 1). The shell is divided into $J$ unity-volumes $Q_j$, $\{j = 0, \dots, J-1\}$, where $J$ is the number of facets of an icosahedron centered at the Earth's center (see Chambodut *et al.*, 2005).

We note here, that a discretization based on subdivisions of the icosahedron does not provide optimal configuration, with respect to several relevant metrics (Katanforoush and Shahshahani, 2003). In this aspect other alternative distributions, such as the one based on polar coordinates, can be considered. On the other hand, for the purpose of our method, we do not need nearly perfect “even distribution” of the shell segments. We just need a practical way to split the data onto smaller subsets. Each unity-volume, assimilable to triangular prisms with two curved opposed faces, contains a certain amount of measurement points. If the size of a given prism is small enough, the data points inside are located in the vicinity of one another and the measured field appears smooth. Every shell segment $Q_j$ is defined such that it obtains the property of compact geometry: its thickness is chosen to be compatible with the size of the icosahedron edges at the satellite altitude. If the shell defined by the data is large, another layer of unity volumes is added in to the radial direction. The aim is to obtain unity-volumes with reasonable shapes, neither too flat and thin, nor too long and thin. For a given dataset, a balance between the number of prisms $J$ (directly related to their size), the number of data per prism $D_j$ and the order of multipole expansion used inside it should be determined. For example, in case of a few data in a large unity-volume, this configuration may allow to recover large wavelengths potential field, and thus leads to the use of low degree polynomials. Contrary, in case of a large number of data in a small unity-volume, this configuration may allow to recover small wavelengths potential field, and thus leads to the use of higher degree polynomials.
mials. Finally, there is a trade-off to find between a good approximation of the functions, that requests small enough cells, and a large enough number of observations within one cell. Indeed, the larger the number of observations in one prism is, the smaller the bias created by the approximation is.

Inside each prism, we define a special point called center of mass. We consider the following two definitions of such points:

- Case 1: Data independent—the definition is given by the geometrical center of the prism. Depending on the data distribution this case can lead to residuals with spatial structure linked to the altitude of the measurements.

- Case 2: Data dependent—the definition of the center of mass is given as the mean value of all the points inside. The advantage of this definition is related to a better spatial distribution of the residuals, the disadvantage comes from the fact that the local coefficient matrix (see next section) can not be precomputed independently from the dataset used in the model.

### 3. General Development into Local Multi-polar Expansions

We consider the following kind of linear model, where an arbitrary function, $s(x)$ is expanded into a family of modeling functions $g_n$, with coefficients $\gamma_n, n = 0, \dots, N-1$,

$$s(x) = \sum_{n=0}^{N-1} \gamma_n g_n(x). \qquad (1)$$

In practice the $g_n$ are typically spherical harmonics or wavelets.

We suppose that the measurements $f_k, k = 0, \dots, K-1$, at the observational points $x_k$, are the “true” values $m_k$ of the measured physical quantity to which some errors $e_k$ (noise due to instrumental errors and/or errors of positioning) are added:

$$f_k = m_k + e_k, \qquad (2)$$

where $e_k$ are Gaussian random variables with covariance matrix $\Sigma = (\Sigma_{k,k'})$ of size $K \times K$.

We also suppose a known a priori information about the model, expressed through a given quadratic form $\Gamma$. The entries of the $N \times N$ square matrix $\Gamma$ are then defined as $\Gamma_{n,n'} = \Gamma(g_n, g_{n'})$. This regularization matrix is used when the condition number of the normal, unregularized system is big. It can be suppressed for a well-conditioned system.

The coefficient vector $\gamma = (\gamma_0, \gamma_1, \dots, \gamma_{N-1})^T$ of the function $s(x)$ from (1), which minimizes the quadratic form $\Gamma$ under the constraint that the error has the supposed known covariance, can be computed as the solution of the following normal equation:

$$(\mathbf{F}^* \mathbf{\Sigma}^{-1} \mathbf{F} + \lambda \mathbf{\Gamma}) \mathbf{\gamma} = \mathbf{F}^* \mathbf{\Sigma}^{-1} \mathbf{f}, \qquad (3)$$

where $\mathbf{F} = (F_{k,n}) = (g_n(x_k))$ is the $K \times N$ system matrix of the modeling functions evaluated at the measurement points, $\mathbf{f} = (f_0, f_1, \dots, f_{K-1})^T$ is the observation vector and $\lambda \in \mathbb{R}$ is a Lagrange multiplier. In other words we

compute the coefficient vector $\gamma$ by solving a classical regularized weighted least-squares problem, with regularization parameter $\lambda$.

Often, since the modeling functions are smooth and the measurement points are dense, the functions $g_n$ are oversampled by the high number of measurement points. We take advantage of the intrinsic correlations between the columns of the system matrix of Eq. (3) and replace it by a numerically equivalent system. The central idea is to replace the information contained in the $f_k$ observation by a suitably agglomerated derived quantities, which from numerical point of view will have the same properties as the original data. In this way, we achieve a speed up in the computation of the normal equation.

In detail, we proceed as follows. Suppose that each function $g_n$ can be approximated, around a given point at satellite altitude, by a linear combination of simpler functions $h_\ell, \ell = 0, \dots, L-1$, where $L \ll K$ is a free parameter, which we discuss later. Such a representation is always possible. One can think of a local expansion into spherical harmonics, i.e. a linear combination of harmonic polynomials up to a certain degree, restricted to some suitably defined volumes at the satellite altitude. Inside each unity-volume, we approximate the values of all modeling functions $g_n$, around the center of mass point $x_0$, by the formula

$$g_n(x) \approx \sum_{\ell=0}^{L-1} C_{n,\ell}(x_0) h_\ell(x - x_0). \qquad (4)$$

Each local coefficient $C_{n,\ell}(x_0)$ depends on the type of the modeling functions used in the model, the order of the local approximation and the tiling of the measurement space. Explicit formulas for the local functions $h_\ell$ and the corresponding local coefficients $C_{n,\ell}$, when real and complex spherical harmonics are used as modeling functions, are given in Appendix A. Similar type of local approximations were also used, in a slightly different context, by Strykowski (2006) to speed-up the computation of the gravity field generated by an arbitrary distribution of mass, in forward gravity field modeling.

If the difference between the left and the right-hand side of Eq. (4) is a fraction of the measurement noise, a model based on the functions from the right-hand side cannot be distinguished from a model based on the functions from the left-hand side. The representation (4) allows us to factorize the system matrix $\mathbf{F}$ as follows:

$$\mathbf{F} \approx \mathbf{H} \mathbf{C}, \qquad (5)$$

where $\mathbf{H} = (H_{k,\ell}) = (h_\ell(x_k - x_0))$ is a block diagonal $K \times JL$ system matrix of local functions evaluated at the measurement points, $\mathbf{C}$ is a $JL \times N$ matrix of local coefficients and $J$ is the total number of volumes used in the space discretization (see Section 2). Based on the above factorization, the normal equations (3) can be rewritten as:

$$(\mathbf{C}^* \tilde{\mathbf{\Sigma}}^{-1} \mathbf{C} + \lambda \mathbf{\Gamma}) \tilde{\mathbf{\gamma}} = \mathbf{C}^* \tilde{\mathbf{\Sigma}}^{-1} \tilde{\mathbf{f}}, \qquad (6)$$

where $\tilde{\mathbf{\Sigma}}^{-1} = \mathbf{H}^* \mathbf{\Sigma}^{-1} \mathbf{H}$ is a $JL \times JL$ matrix and $\tilde{\mathbf{f}} = \tilde{\mathbf{\Sigma}} \mathbf{H}^* \mathbf{\Sigma}^{-1} \mathbf{f}$ is a vector of local moments of size $JL$.
The exact solution $\gamma$ of Eq. (3) can be computed from the solution of (6) by using an iterative method of stationary type (see Section 4.1). The advantage of such an approach is that we can interpret the matrix in the left hand side of Eq. (6) as a preconditioner to the original problem. In this sense it is not unusual to chose an approximation to the matrix of the underlying problem and use it to construct an iterative method. As long as the spectral radii of the approximation matrix is smaller than 1, the iterative method will convergence to the exact solution of the original problem. Clearly, the quality of the approximation, controlled by the free parameters: size of the prisms and degree of the local approximations, from one side and the speed of convergence of the method from other side are closely related. The better the quality of the approximation is the higher the speed of convergence should be. In the same time, since the original normal equation is usually not a well posted problem, using an approximation, which is not that precise will actually have the effect of an additional regularization. In Section 4, we have demonstrated numerically that 2-nd and 3-rd degree polynomial approximation of the basis function, together with the space discretization introduced in Section 2, provide satisfactory results for all considered examples. These choices are just one possibility, which proofs the applicability of our method. Clearly, better combinations might be possible. The optimal choice for all free parameters, which we have introduced in our method is currently under investigation.

### 3.1 Operation count

What do we gain by using the above factorization? To answer this question, let us first suppose that all local functions are restrictions to some spherical prisms $Q_j$, $j = 0, \dots, J-1$ of spherical harmonics of low degree and let $Q_{j(n)}$ denotes the support of $H_n$. The maximal number of $h_n$ that share the same support, we denote by $E$

$$E_j = \#\{n : \text{support } h_n = Q_j\}, \quad E = \max\{E_j\}. \quad (7)$$

Let us also denote by $D_j$ the maximal number of data points inside a given prism $Q_j$

$$D_j = \#\{k : x_k \in Q_j\}, \quad D = \max_j\{D_j\}. \quad (8)$$

The following relations are obviously satisfied

$$D \le K \le DJ, \quad \text{and} \quad E \le L \le EJ. \quad (9)$$

Note that, since a matrix multiplication of an $A \times B$ matrix with a $B \times C$ matrix takes $O(ABC)$ operations in general, for a full matrix $\Sigma$, the total number of operations needed to compute the normal matrix by means of factorization (5) may actually increase. However, if the errors are not correlated between different $Q_i$ (or only between neighboring $Q_j$), we may gain in computation time. In this case, the direct computation of the normal matrix requires

$$\mathrm{opcount}_1 = KN(N+1)/2 \quad (10)$$

operations.

On the other hand, we can take advantage from the special block diagonal structure of the matrix $\mathbf{H}$ (see Appendix A) and first compute the matrix $\mathbf{X} = \mathbf{H}^*\mathbf{\Sigma}^{-1}\mathbf{H}$,

followed by the computation of the matrices $\mathbf{Y} = \mathbf{XC}$ and $\mathbf{C}^*\mathbf{Y} = \mathbf{C}^*\mathbf{H}^*\mathbf{\Sigma}^{-1}\mathbf{H}\mathbf{C}$. This process takes at most

$$\mathrm{opcount}_2 \leq DJL(L+1)/2 + JL^2N + JLN(N+1)/2. \quad (11)$$

Clearly, the opcount$_2$ is dominated by the largest of its first and last terms. We can conclude that, in order to achieve computational savings, the choice of the parameters $L$ and $J$ should be done such that $R \ll K$, where: $R = JL$, is the total number of local coefficients used in the approximation (5), or in other words the rows number in $\mathbf{C}$. Since neither $L$ nor $J$ depends on $K$, the last requirement can be fulfilled by simply increasing the number of measurements used in the model.

### 3.2 Data storage

Additionally to the computational speed-up, achieved by using local multi-polar approximations of the modeling functions, the proposed approach leads to memory savings. This is due to the fact that instead of solving the original inverse problem with system matrix $\mathbf{F}$, we solve an equivalent problem with system matrix $\mathbf{C}$, which for $R \ll K$ has much smaller dimensionality. In practice, this allows to either increase the number of observations used in the model, or to increase the number of the modeling functions. We also note, that since the entries of the matrix $\mathbf{C}$ depend only on the number of the local functions used in the approximation and the total number of the modeling functions, for a space discretization (which uses a measurement independent definition of the center of mass, see Section 2), the matrix $\mathbf{C}$ can be precomputed once and for all. This fact can lead to additional computational savings, when the size of the unity-volumes used in the space discretization is small and the number of data inside each is large. Further more, one can also compute and store the SVD (or the QR) decomposition of the matrix $\mathbf{C}$ and then use it to find the least-squares solution of the original inverse problem. This is in fact the approach used in our numerical experiments.

## 4. Applications for Geopotential Field Modeling

In this section we test the proposed local multi-polar expansion approach for modeling the Earth's gravity (Section 4.1) and magnetic (Section 4.2) potential fields. As modeling functions, we use spherical harmonics and their approximations by local harmonic polynomials up to degree 2, $L = 9$ (see Appendix A). The space discretization is based on generation of six subdivisions of the icosahedron, with 20480 total number of facets at the Earth's surface.

An important question regarding the use of the multipolar expansion technique is: Depending on the size of the considered dataset, how should we choose the size of the unity-volumes (the number of data per unity-volume) and the maximal order of the local harmonic polynomials, in order to obtain optimal performance? It is clear that the higher the resolution of the space discretization is, the better the approximation will be. On the other hand, the maximal reduction of the size of the underlying problem, is obtained for the lowest possible resolution. Therefore, we need to find a balance between the resolution of the space discretization and the degree of the local expansion used in the model.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Ratio in % between the sizes of the original and the approximated problems versus the approximation relative error for $L = 4, 9, 16, 25$.</figcaption></figure>

Figure 2, shows the maximum relative error between all spherical harmonics up to degree 30 and their corresponding local multi-polar approximations, versus the ratio (in %) between the sizes of the underlying and the approximated problems (size reduction), for the dataset used in Section 4.2. For small and moderate size datasets, low orders local approximations, of first ($L = 4$) and second ($L = 9$) orders, are preferable. Higher order approximations would be beneficial, when the dataset is large enough, so the number of local coefficients per unity-volume can be compensated by the amount of data approximated in each cell. In this case we can benefit from the better quality of the higher order approximation.

The inversion codes, for the direct and the approximated problems, have been implemented in two programming languages: Fortran 90, for the scalar gravity model, and C++, for the vector magnetic model. The tests have been performed on a Intel Xeon™ 2 × 3.8 GHz CPU workstation, with 8 GB main memory and on a Intel Xeon™ 3.4 GHz CPU workstation, with 4 GB main memory, respectively.

## 4.1 Gravity field model

We first show an example of application of the method for gravity field modeling from GRACE. We consider the case of a monthly gravity field recovery up to degree 50 from GRACE simulated data sampled at 5 seconds intervals. It is indeed difficult to recover higher degrees from only one month of measurements, and consequently, the monthly geoid routinely provided by different teams are truncated at degree 50 or 60 (Biancale *et al.*, 2007; Bettadpur, 2008). In order to control the various sources of error and clearly show the performances of the multipolar approximations, we set up a simulation with very high accuracy requirements, closer to GRACE target accuracy than GRACE real accuracy. Moreover, our aim here is not to produce a gravity field from real datasets, that can be made available for users, but to show how the presented method works.

Our simulation data are inter-satellite potential differences as provided by the energy balance approach from the GRACE inter-satellite range rate (KBR) observables (Jekeli, 1999; Rowlands *et al.*, 2002; Han *et al.*, 2003). A realistic noise was added, taking into account the noise

in the KBR measurements, and errors arising from ocean tides mismodelling and from the continental hydrology signal, considered as a noise in this simulation. For that, we first generated one month of GRACE satellites orbits using the GINS software developed by GRGS (Toulouse, France). The GINS simulated measurements consist of GPS positions and velocitites of the GRACE satellites, accelerometer records accounting for the effects of non-conservative forces, and KBR satellite-to-satellite tracking measurements. The orbits are based on the EIGEN-GL04S static reference field up to degree 50 (Biancale *et al.*, 2007). The errors were introduced in the GINS simulated measurements as follows. First, an error on the intersatellites KBR range-rate measurements is introduced, using an empirical parameterization of the K-Band ranging system data at the period of revolution (i.e., 15 bias per day plus sinusoidal terms). The impact of this error on intersatellite potential differences (in rms) is about 0.001 m²/s². Second, a dynamical error was introduced, reflecting a mismodelling of the oceanic tides components of the time-varying gravity field and a continental hydrology loading. The ocean tide mismodelling error is taken as one tenth of the potential differences between the GOT2007 and FES2004 ocean tide models. It leads to an increase of the KBR errors by a factor of 2. The continental hydrology loading corresponds to the August 2003 total water contents from the WGHM model by Döll *et al.* (2003). We then simulate along the orbits the perturbations on the inter-satellite potential differences associated with these sources of noise by applying the time-integrated energy balance approach (see Jekeli, 1999; Rowlands *et al.*, 2002; Han *et al.*, 2003). The contribution of the static disturbing gravity potential from the EIGEN-GL04S model was directly sampled along the orbit. Figure 3 (left) shows the obtained 518400 synthetic data, and Fig. 3 (right) illustrates the synthetic noise. The noise rms is 0.0022 m²/s², corresponding to about 0.2 mm on the geoid.

From the least-squares adjustment of this synthetic dataset, we computed a spherical harmonics model of the gravity potential up to degrees 50, without applying a stabilization. Note that we do not use any a-priori gravity model (that is to say, our first guess of the solution is zero). We first perform an exact least-squares adjustment of the system. The map of residuals with respect to the synthetic data is represented on Fig. 4 (left), and corresponds to the KBR data noise and the tide model error. The rms of residuals amounts 0.002 m²/s², a little less since part of our synthetic noise has aliased into the solution. Figure 4 (right) shows the difference between the recovered spherical harmonics coefficients and the EIGEN-GL04S coefficients. The rms of residuals amounts 0.0023 m²/s², at the level of the synthetic noise. These residuals are clearly structured, reflecting the structure of the noise.

We then solved the normal system for different levels of approximation. Note that, in the case of inter-satellite potential differences, the entries of the system matrix are given by the difference between the bases functions at the two satellites positions for a given time $t$, namely ($F_{k,n}$) = ($g_n(x_k^{\text{Grace A}}(t))-g_n(x_k^{\text{Grace B}}(t))$). This leads to a band block diagonal $\tilde{\Sigma}^{-1}$ matrix.

To solve the approximated system, we apply an iteration.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Map of the simulated GRACE inter-satellite potential differences up to degree 50. Total disturbing gravity potential differences (**left**); Simulated noise (KBR measurements noise, error on the tide model and continental hydrology variations) (**right**).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Map of the residuals between the reconstructed spherical harmonics model up to degree 50 and the original dataset (**left**); Differences between the reconstructed spherical harmonics coefficients and the EIGEN-GL04S spherical harmonics coefficients (**right**). No size reduction is applied.</figcaption></figure>

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

tive approach for two reasons. The first one is that itera-
tive solvers are often considered for gravity field recovery
from satellite measurements, starting from a first guess of
the solution. In such context, exact inversions may be effi-
ciently replaced with approximate ones. The second reason
is that, by a single inversion, we may not be able to reach
the required level of accuracy, for a reasonable compres-
sion rate. However, using the set of spherical harmonics
coefficients obtained as the solution of the system, we may
build a first guess of the potential, and remove it from the
synthetic dataset. The residuals are then treated as input
data and adjusted during a second iteration, involving the
same approximated normal matrix and a new reduced data
vector **f**. Since at each step the approximated normal ma-
trix is unchanged, the extra computational cost involved,
for re-computing the vector **f** and the new solution of the
normal equation (6), is kept low. This process can be it-
erated until the desired level of precision is reached. Such
iterative approach allows to keep the approximation error
under control, and to ensure that it is smaller than the noise
level. Moreover, it provides a great flexibility, since both
quick-look solutions may be computed with high compres-

sion rates and few iterations, and precise solutions may be
obtained with lower compression rates and/or a larger num-
ber of iterations.

The approximation parameters are summarized in Ta-
ble 1. We consider local multi-polar expansions of spher-
ical harmonics of local orders 1—84% size reduction, and
2—64% size reduction—on a space discretization involving
one layer of 20480 prismatic cells centered at the average
data altitude.

Figure 5 (**left**) shows the spatial map of residuals in the 64% reduction case, after 4 iterations. The rms of the residuals is equal to 0.002 m²/s², and the residuals look very similar to the ones obtained without any approximation. The histograms of residuals for different levels of approximation are represented on Fig. 7, and one cannot distinguish between the result of an exact solution and the result of this approximated solution. Figure 5 (**right**) shows the difference between the recovered spherical harmonics coefficients and the ones obtained from the exact inversion. The rms of residuals amounts 0.0003 m²/s², and the discrepancies are mainly observed on the highest degrees spherical harmonics coefficients degrees (above 30), which make
<table><caption>Table 1. Local approximations applied in gravity modelling tests. A subdivision of space with one layer with total of 20480 prisms, all of which are filled, is used.</caption><thead><tr><th>Expansion order</th><th>1</th><th>2</th></tr></thead><tbody><tr><td>Number of reduced data/prism</td><td>4</td><td>9</td></tr><tr><td>Average number of original data/prism</td><td>23</td><td>23</td></tr><tr><td>Total number of reduced data</td><td>81920</td><td>184320</td></tr><tr><td>Total number of original data</td><td>518400</td><td>518400</td></tr><tr><td>Reduction in the size of the underlying problem</td><td>84%</td><td>64%</td></tr><tr><td>Nb of reduced data/SH coefficient (deg. 50)</td><td>31</td><td>71</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Map of the residuals between the reconstructed spherical harmonics model up to degree 50 and the original dataset (**left**); Differences between the reconstructed spherical harmonics coefficients and the spherical harmonics coefficients obtained from the exact inversion (**right**). A 64% size reduction is applied, and 4 iterations have been performed.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Map of the residuals between the reconstructed spherical harmonics model up to degree 50 and the original dataset (**left**); Differences between the reconstructed spherical harmonics coefficients and the spherical harmonics coefficients obtained from the exact inversion (**right**). A 84% size reduction is applied, and 9 iterations have been performed.</figcaption></figure>

sense since the lower degrees spherical harmonics are better approximated for a given size of the cells and local multipolar development order. The spherical harmonics components for the degrees lower than 30, that contain a lot

of signal, are almost perfectly recovered. Finally, this figure also shows that the multipolar approximation error after four iterations is lower than the noise level. If we continued iterating, the solution of the approximated system would
<figure><img src="image_9.png" /><figcaption>Fig. 7. Histograms of the spatial residuals between the simulated data and the reconstructed data in the 3 cases: exact inversion, 64% size reduction and 4 iterations, 84% size reduction and 9 iterations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Evolution of the rms of the spatial residuals between the simulated data and the reconstructed data as a function of the iteration number in the cases: 64% size reduction, 84% size reduction, no size reduction.</figcaption></figure>

converge towards the exact solution.

If we degrade the quality of approximation of the nor-
mal system, a larger number of iterations are requested to
meet the same accuracy requirements. Thus, in the 84%
reduction case, 9 iterations are necessary to reach the syn-
thetic noise level. Figure 6 (left) shows the spatial map of
residuals; the rms of the residuals is equal to 0.0027 m²/s².
There is a very slight difference between this solution and
the previous ones, as shown on the histogram plot of Fig. 7.
Accordingly, the difference between the recovered spherical
harmonics coefficients and the ones obtained from the exact
inversion is slightly larger, but still at the level of synthetic
noise, as shown on Fig. 6 (right). Again the largest discrep-
ancies are observed for the highest degrees. Note that the
recovery of the spherical harmonics coefficients for degrees
lower than 30 is again near perfect, even for a larger com-
pression rate. After adding 4 more iterations, the difference

between the recovered coefficients and the ones obtained
from the exact inversion falls below the noise level. Iter-
ating allows to introduce higher size reduction rates, while
keeping the quality of the approximation.

Finally, Fig. 8 shows how the rms of the spatial residuals
changes during the iterative resolution process. Our accu-
racy demands in this simulation were very high, but inter-
satellite potential differences that can be built from the real
GRACE data are known to have a lower accuracy, at the
level of 0.01 m²/s² (Shum *et al.*, 2004). Figure 8 shows that
we reach this accuracy at the second iteration step for the
64% reduction case, and at the third iteration step only for
the 84% reduction case. In addition, we did not use any
apriori guess of the spherical harmonics coefficients from a
background model in our simulations, contrary to what is
often done, especially when the system equations have to
be linearized (Bettadpur, 2007). Using a non-zero *a-priori*
<table><caption>Table 2. Local approximations applied in magnetic modelling tests. A subdivision of space with one layer with total of 20480 prisms, of which 7652 filled, is used.</caption><thead><tr><th>Expansion order</th><th>0</th><th>1</th></tr></thead><tbody><tr><td>Number of reduced data/prism</td><td>1</td><td>8</td></tr><tr><td>Number of original data/prism</td><td>9</td><td>9</td></tr><tr><td>Total number of reduced data</td><td>7652</td><td>22956</td></tr><tr><td>Total number of original data</td><td>70511</td><td>70511</td></tr><tr><td>Reduction in the size of the underlying problem</td><td>89%</td><td>68%</td></tr><tr><td>Nb of reduced data/SH coefficient (deg. 30)</td><td>8</td><td>24</td></tr><tr><td>Nb of reduced data/SH coefficient (deg. 60)</td><td>2</td><td>6</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Histograms of the residuals between the reconstructed spherical harmonics model and the original dataset after the second iteration step. Degree 30 model (**left**); degree 60 model (**right**). From top to bottom: X, Y and Z components of the geomagnetic field.</figcaption></figure>

on the solution, the perturbation on the spherical harmonics coefficients to be estimated would be much smaller and it is very likely that it could be performed at the level of GRACE present-day accuracy within one iteration only.

## 4.2 Magnetic field model

In the following, we present tests on the ability of the considered local multi-polar approximation technique to model the Earth's magnetic field. We built two different
<figure><img src="image_12.png" /><figcaption>Fig. 10. Maps of the residuals between the reconstructed spherical harmonics model up to degree 30 and the original dataset after the second iteration step. 89% size reduction (**left**); 68% size reduction (**right**). From **top** to **bottom**: X, Y and Z components of the geomagnetic field.</figcaption></figure>

sets of vector synthetic data at real positions of the CHAMP satellite for January 2007. Values of three components of the magnetic vector field were computed by using the spherical harmonics POMME 3.1 model (Maus *et al.*, 2006) up to degrees 30 and 60, respectively. In addition a Gaussian synthetic noise, with a standard deviation of 2 nT, for the degree 30, and 3 nT, for the degree 60 model was added to each of the field components. Due to the vector nature of the magnetic potential, the size of the synthetic dataset is taken smaller than the one used in the previous example. In order to keep the size reduction in approximately the same levels, as before, we considered constant ($L = 1$) and linear ($L = 4$) local harmonic polynomials approximations of the basis functions. Table 2 details the main parameters of the computations for the magnetic field model.

The results obtained from the least-squares fitting of the two synthetic datasets are based on the same orders of app

proximations for both magnetic datasets. The space discretization evolves on a layer of 20480 prismatic cells, of which only 7652 are filled. The center of mass of each prism is taken by averaging the data inside. The solution of the normal equation is found by using the iterative technique described in Section 4.1.

Figures 9 (**left**) and 10 show the histograms and the maps of the residuals between the reconstructed spherical harmonics model up to degree 30 and the original dataset, for X (northern), Y (eastern), Z (down vertical) components of the magnetic field. The results from the approximation models are obtained after two iteration steps. For comparison, the residuals computed by direct inversion, without using approximations of the modeling functions, are also given in Fig. 11. The fit obtained with the 68% size reduction of the underlying equation is at the level of the added synthetic noise, for all three components of the vector field.
<figure><img src="image_13.png" /><figcaption>Fig. 11. Maps of the residuals between the reconstructed spherical harmonics model up to degree 30 and the original dataset. No size reduction is applied. From *top* to *bottom*: X, Y and Z components of the geomagnetic field.</figcaption></figure>

When the size of the normal equation is reduced by 89%, the residuals are slightly larger from the introduced noise level.

Results for the three different components of the magnetic field, obtained from a degree 60 spherical harmonics model are shown in Figs. 9 (right) and 12. For the 68% reduction case, we reach the level of the applied synthetic noise, while for the 89% reduction case, we are above it. An other iteration step is required to improve the residuals in this case. The spatial distribution of the residuals is again related to the data distribution around the centers of mass of each prism. The errors in the approximated values of the spherical harmonics of high degrees (above 30) are more strongly related to the distance between the data points and the center of mass.

Finally, we comment the fact that the histograms for var-

ious compression rates of the magnetic field look more similar to each other, than those for the gravity field. This can be explained with the smoother nature of the magnetic field in comparison with that for the gravity field.

**5. Conclusion and Prospects**

In this paper, we have developed and tested a new method, which allows to compute gravity and magnetic field models from very large datasets in a fast and economic way. Our method uses a 3D tiling of the space and it is based on local harmonic polynomial approximations of the modeling functions inside each space segment. Introducing such local approximations is shown to be equivalent on replacing the original dataset with a set of local moments of various orders. The proposed approach allows to represent the original system matrix as a product of two smaller matrices: a local coefficient matrix and a matrix of local functions. Thus, the size of the resulting inverse problem can be significantly reduced. The local coefficient matrix C can be computed and stored for later uses before solving the least-squares problem. Due to the block-diagonal structure of the system matrix of local functions **H**, the computation of the modified covariance matrix $\tilde{\Sigma}^{-1}$ and the vector of local moments $\tilde{\tau}$ requires little computational effort and can be performed relatively fast. The considered approach generalizes the standard averaging technique and introduces the use of higher order local moments in potential field modeling. Depending on the chosen precision of the approximation, we can control the speed and the quality of the inversion. We may thus compute fast quick-look solutions or highly precise models with a great deal of flexibility. The tests performed in this paper show that we can practically control the amount of errors introduced by the approximations in the model. Distortions are small, and become negligible when the process is iterated. Our method is very general: it can be applied to any kind of modeling functions like wavelets or other radial basis functions. The local moments of the original data are computed around the centers of the data clusters inside the unity-volumes, but one may imagine in the future to design tilings of the space that could be more closely related to the data distribution instead of a fixed, uniform one. One may also consider local refinements of the tilings, in areas where a higher precision is requested, for instance because of the presence of important small scale features. This may be of primary interest for local modeling using wavelets, to increase the resolution, where needed. Finally, this paper opens broad ways in potential fields modeling from large and various datasets, at both global and local scales. The method can already be applied to the processing of the currently operating satellite missions such as the GRACE mission. In the future, we would like to apply this method to process larger datasets from several satellites, and/or include surface measurements. The method should be of primary interest for the forthcoming GOCE and Swarm data assimilation, using spherical harmonics or wavelets as modeling functions. In the case of the GOCE mission, the measurements will be less smooth, but the noise level will be larger, requiring an appropriate tuning of the approximation parameters and number of iterations. A major challenge is the treatment
<figure><img src="image_14.png" /><figcaption>Fig. 12. Maps of the residuals between the reconstructed spherical harmonics model up to degree 60 and the original dataset after the second iteration step. 89% size reduction (**left**); 68% size reduction (**right**). From **top** to **bottom**: X, Y and Z components of the geomagnetic field.</figcaption></figure>

of the coloured noise in the gradiometer measurements, es-
pecially in the space-wise approach. A possible way to
implement our method would be to pre-filter the measure-
ments before performing the space-wise inversion, similarly
to the approach proposed by Migliaccio *et al.* (2004), where
a Wiener filtering of the observation is performed along-
track before computing the gravity model. However, there
is a risk that the pre-filtering step affects not only the noise,
but also the signal, so it is necessary to assess the quality
of the filtering step based on the results of the gravity field
inversion, by comparing the *a-priori* and *a-posteriori* spec-
tra of both signal and noise, and iterate the process. A full
study of these aspects should be the topic of further work.
Finally, we indeed intend to apply this method for wavelet
modeling of potential fields, at both global and local scales.
Coupled with iterative resolution techniques for large sys-
tems, such as domain decomposition methods, the method

will lead to fast and accurate solvers. We may also think of
introducing a time-dependency in the local approximations,
in view of 4D-potential fields modeling. To conclude with,
the obtained models, constrained with the largest possible
amount of data at various altitude within a fully consistent
framework, will provide us with a high quality image of the
Earth's potential fields. When jointly analyzed with other
geophysical observables, those models will contribute to a
better understanding of the structure and dynamics of our
planet.

**Acknowledgments.** We thank Richard Biancale, from CNES/GRGS, Toulouse, France, for providing the GRACE simulated orbits data computed from the GINS software. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the framework of the projects BMBF/DFG "GEOTECHNOLOGIEN" and the DFG project KO 2870/3-1. Two of us (M. Holschneider and E. Schöll) also thank the
SFB 555 for financial support. I. Panet was supported by a CNES fellowship. This is IPGP contribution 2567.

## Appendix A. Local Multi-polar Expansions of Spherical Harmonics

Here we give explicit formulas for the complex and the real local multi-polar expansions of a spherical harmonic function in terms of harmonic polynomial. The real case formulas have been used in our numerical experiments. The following notations are used (see Fig. A.1): $\vec{x}_{0,j}$ is the centre of mass of the $j$-th unity-volume; $\vec{x}_{i,j}$ is the $i$-th measurement point in the $j$-th unity-volume; $\vec{x}_{i,j} - \vec{x}_{0,j}$ gives the local coordinates of the point $\vec{x}_{i,j}$ with respect to the coordinate system centred at the point $\vec{x}_{0,j}$. Obviously, for all unity-volumes the following relation holds $\|\vec{x}_{i,j} - \vec{x}_{0,j}\| > \|\vec{x}_{i,j} - \vec{x}_{0,j}\|$.

### A.1 Complex case

The value of each complex spherical harmonic at the point $\vec{x}_{i,j} \in Q_j$ is given by the following local multi-polar expansion around the point $\vec{x}_{0,j}$ (Epton and Dembart, 1995)

$$O_{\ell}^{m}(\vec{x}_{i,j}) = \sum_{\ell'=0}^{\infty} \sum_{m'=-\ell'}^{\ell'} I_{\ell'}^{m'}(\vec{x}_{i,j} - \vec{x}_{0,j}) O_{(\ell+\ell')}(\vec{x}_{0,j}), \quad (A.1)$$

where

$$I_{\ell}^{m}(\vec{x}_{i,j}) = i^{-|m|} A_{\ell}^{m} r^{\ell} Y_{\ell}^{m}(\theta, \phi), \quad (A.2)$$

$$O_{\ell}^{m}(\vec{x}_{i,j}) = \frac{(-1)^{\ell} i^{|m|} Y_{\ell}^{m}(\theta, \phi)}{A_{\ell}^{m} r^{\ell+1}}, \quad (A.3)$$

$$A_{\ell}^{m} = \frac{(-1)^{\ell}}{\sqrt{(\ell - m)!(\ell + m)!}}, \quad (A.4)$$

with $(r, \theta, \phi)$ and $(r_0, \theta_0, \phi_0)$ being the spherical coordinates of the points $\vec{x}_{i,j}$ and $\vec{x}_{0,j}$, respectively.

Equation (A.1) can also be written in the following more explicit form

$$\frac{Y_{\ell}^{m}(\theta, \phi)}{r^{\ell+1}} = \frac{A_{\ell}^{m}}{(-1)^{\ell} i^{|m|}} \sum_{\ell'=0}^{\infty} \sum_{m'=-\ell'}^{\ell'} \frac{A_{\ell'}^{m'}}{i^{|m'|}} h^{\ell'} Y_{\ell'}^{m'}(\alpha, \beta) \cdot \frac{(-1)^{(\ell+\ell')i|m-m'|} Y_{(\ell+\ell')}(\theta_0, \phi_0)}{A_{(\ell+\ell')}^{(m-m')}} r_0^{(\ell+\ell'+1)}, \quad (A.5)$$

where $(h, \alpha, \beta)$ are the spherical coordinates of the point $\vec{x}_{i,j} - \vec{x}_{0,j}$. The vector case formulas, used in Section 4.2, are obtained by computing the gradient of Eq. (A.5).

### A.2 Real case

Here we give an equivalent formulation of the local multi-polar expansion (A.5), which uses only real arithmetics. We denote by $\Psi_{\ell}^{m}$ the real spherical harmonics of degree $\ell$ and order $m$ defined as

$$\Psi_{\ell}^{m} = \begin{cases} \frac{Y_{\ell}^{m} + \bar{Y}_{\ell}^{m}}{\sqrt{2}}; & m > 0, \\ Y_{\ell}^{0}; & m = 0, \\ i \frac{Y_{\ell}^{m} - \bar{Y}_{\ell}^{m}}{\sqrt{2}}; & m < 0. \end{cases} \quad (A.6)$$

Before we present how a real spherical harmonic can be approximated, up to an arbitrary given order $L$ ($L = (L +$

<figure><img src="image_15.png" /><figcaption>Fig. A.1. Position vectors in global $(X, Y, Z)$ and local $(X', Y', Z')$ coordinate systems.</figcaption></figure>

1)²), in terms of local harmonic polynomials, we introduce the following notations. For each spherical harmonic $\Psi_{\ell}^{m}$ and for each $\ell' = 0, 1, \dots, L$ and $m' = -\ell', \dots, \ell'$, we define the coefficients

$$A_{\ell, \ell'}^{m, m'} (r_0, \theta_0) = 2^{\nu} (-1)^{\ell+\ell'} \frac{R_E^{\ell+1}}{r_0^{\ell+\ell+1}} \frac{(\ell+\ell' - |m-m'|)!}{(\ell' + |m'|)!} \cdot P_{(\ell+\ell')}^{m-m'} (\cos \theta_0), \quad (A.7)$$

where $(r_0, \theta_0, \phi_0)$ are the spherical coordinates of the centre of mass $\vec{x}_{0,j}$, $R_E$ is the radius of the Earth and $\epsilon$ is defined as

$$\epsilon = \begin{cases} m; & |m'| > |m| \quad mm' > 0, \\ m'; & |m'| \le |m| \quad mm' > 0, \\ 0; & 0 \le |m'| \le |m| \quad mm' < 0, \quad m = 0, \\ m' - m; & |m'| > |m| \quad mm' < 0. \end{cases} \quad (A.8)$$

The exponent $\nu$ is equal to 0 for $m = 0$ and to $\frac{1}{2}$ for all other $m$. We also introduce the related coefficients

$$\begin{aligned} & AC_{\ell, \ell'}^{m, m'} (r_0, \theta_0, \phi_0) = A_{\ell, \ell'}^{m, m'} (r_0, \theta_0) \cos((m - m')\phi_0), \\ & AS_{\ell, \ell'}^{m, m'} (r_0, \theta_0, \phi_0) = A_{\ell, \ell'}^{m, m'} (r_0, \theta_0) \sin((m - m')\phi_0). \end{aligned} \quad (A.9)$$

We are now ready to express the value of a real spherical harmonic at any given point $\vec{x}_{i,j}$, with spherical coordinates $(r, \theta, \phi)$, as a linear combination of real local harmonic polynomials. The following approximation is obtained from (A.5) by taking into account the definition (A.6)

$$\frac{\Psi_{\ell}^{m}(\theta, \phi)}{\left(\frac{r}{R_{\mathrm{E}}}\right)^{\ell+1}} = A_{\ell}^{m} \sum_{\ell'=0}^{\infty} \sum_{m'=-\ell'}^{\ell'} C_{\ell, \ell'}^{m, m'}(r_0, \theta_0, \phi_0) H_{\ell'}^{m'}(h, \alpha, \beta), \quad (A.10)$$

where $A_{\ell}^{m}$ is defined as

$$A_{\ell}^{m} = \frac{(-1)^{m}}{\sqrt{(\ell - m)!(\ell + m)!}}, \quad (A.11)$$
and coefficients $C_{\ell, \ell'}^{m, m'}$ are given by

$$C_{\ell, \ell'}^{m, m'}(r_0, \theta_0, \phi_0) = \begin{cases} AC_{\ell, \ell'}^{m, m'} + AC_{\ell, \ell'}^{m, -m'}; & m' > 0, \\ AC_{\ell, \ell'}^{m, 0}; & m' = 0, m \ge 0, \\ AS_{\ell, \ell'}^{m, m'} - AS_{\ell, \ell'}^{m, -m'}; & m' < 0, \\ -AS_{\ell, \ell'}^{m, m'} - AS_{\ell, \ell'}^{m, -m'}; & m' > 0, \\ -AS_{\ell, \ell'}^{m, 0}; & m' = 0, m < 0. \\ AC_{\ell, \ell'}^{m, m'} - AC_{\ell, \ell'}^{m, -m'}; & m' < 0, \end{cases} \quad (A.12)$$

The local harmonic functions $H_{\ell'}^{m'}$ are

$$H_{\ell'}^{m'}(h, \alpha, \beta) = h^{\ell'} P_{\ell'}^{m'|1}(\cos \alpha) \begin{cases} \cos(m'\beta) & m' \ge 0, \\ \sin(-m'\beta) & m' < 0, \end{cases} \quad (A.13)$$

with $(h, \alpha, \beta)$ being the local spherical coordinates of the point $\vec{x}_{i,j}$ with respect to the origin $\vec{x}_{0,j}$.

From (A.10) it is clear that the total number of coefficients used in a real multi-polar approximation is equal to the product of the number of all local functions times the number of all unity-volumes i.e. equal to $R = LJ$.

The elements of the coefficient matrix C and the values of the local functions, which enter in the matrix H from (5) can be computed directly from (A.12) and (A.13) respectively.

## References

Bettadpur, S., Level 2 gravity field product user handbook GRACE 327-734, http://podaac.jpl.nasa.gov/grace, 2007.

Bettadpur, S., Release notes for CSR RL04 GRACE L2 products, http://podaac.jpl.nasa.gov/grace, 2008.

Biancale et al., Five years of gravity variations from GRACE and LAGEOS data at 10-day interval over the period from July 29th 2002 to September 10th 2007, http://bgi.cnes.fr:8110/geoid-variations/README.html, 2007.

Chambodut, A., J. Schwarte, B. Langlais, H. Lühr, and M. Mandea, The selection of data in field modelling, OIST-4 Proceedings, Danish Meteorological Institute Scientific Report 03-09, 31–34, 2003.

Chambodut, A., I. Panet, M. Mandea, M. Diament, M. Holschneider, and O. Jamet, Wavelet frames: an alternative to spherical harmonic representation of potential fields, Geophys. J. Int., 163(3), 875–899, 2005.

Colombo, O. L., Numerical methods for harmonic analysis on the sphere, Reports of the Department of Geodetic Science, 310, Ohio State University, Columbus, 1981.

Colombo, O. L., The global mapping of gravity with two satellites, The Netherlands Geodetic Commission. Publications on Geodesy, 7, Delft, 1984.

Ditmar, P., R. Klees, and F. Kostenko, Fast and accurate computation of spherical harmonics coefficients from satellite gravity gradiometry data, J. Geod., 76, 690–705, 2003.

Döll, P. F., F. Kaspar, and B. Kaspar, A global hydrological model for deriving water availability indicators: model tuning and validation, J. Hydrol., 270, 105–134, 2003.

Epton, M. and B. Dembart, Multipole translation theory for the three-dimensional Laplace and Helmholtz equations, SIAM J. Sci. Comput., 16, 865–897, 1995.

Han, S.-C., Efficient determination of global gravity field from satellite-to-satellite tracking mission, Celestial Mech. Dyn. Astron., 88, 69–102, 2004.

Han, S.-C., C. Jekeli, and C. K. Shum, Static and temporal gravity field recovery using GRACE potential difference observables, Adv. Geosci., 1, 19–26, 2003.

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

Jekeli, C., The determination of gravitational potential differences from satellite-to-satellite tracking, Celestial Mech. Dyn. Astron., 75, 85–101, 1999.

Katanforoush, A. and M. Shahshahani, Distributing points on the sphere, I Exp. Math., 12, 199–209, 2003.

Keller, W., A wavelet approach for the construction of multi-grid solvers for large linear systems, in Proc. AG 2001 Scientific Assembly, 2–7 September 2001, Budapest, 2001.

Klees, R., R. Koop, R. Visser, and J. van der Ijssel, Efficient gravity field recovery from GOCE gravity gradient observations, J. Geod., 74, 561–571, 2000.

Kusche, J., Implementation of multigrid solvers for satellite gravity anomaly recovery, J. Geod., 74, 773–782, 2000.

Langel, R. A. and W. J. Hinze, The Magnetic Field of the Earth's Lithosphere: The Satellite Perspective, Cambridge Univ. Press, New York, 1998.

Lemoine, F. G., S. C. Kenyon, J. K. Factor, R. G. Trimmer, N. K. Pavlis, D. S. Chinn, C. Cox, S. M. Klosko, S. B. Luthcke, M. H. Torrence, Y. M. Wang, R. G. Williamson, E. C. Pavlis, R. H. Rapp, and T. R. Olson, The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96, NASA Technical Paper, 1998-206861, Greenbelt, Maryland, 1998.

Lesur, V., I. Wardinski, M. Rother, and M. Mandea, GRIMM—The GFZ Reference Internal Magnetic Model based on vector satellite and observatory data, Geophys. J. Int., 2007 (in press).

Maus, S., M. Rother, C. Stolle, W. Mai, S. Choi, H. Lühr, D. Cooke, and C. Roth, Third generation of the Potsdam Magnetic Model of the Earth (POMME), Geochem. Geophys. Geosyst., 7, Q07008, 2006.

Maus, S., T. Sazonova, K. Hemant, J. D. Fairhead, and D. Ravat, National Geophysical Data Center candidate for the World Digital Magnetic Anomaly Map, Geochem. Geophys. Geosyst., 8, Q06017, doi:10.1029/2007GC001643, 2007.

Migliaccio, F., M. Reguzzoni, and F. Sanso, Space-wise approach to satellite gravity field determination in the presence of coloured noise, J. Geod., 78, 304–313, 2004.

Migliaccio, F., M. Reguzzoni, F. Sanso, and N. Tselfes, On the use of grid-ded data to estimate potential fields coefficients, Proceedings of the 3rd International GOCE USsers Workshop, Frascati, Italy, 6–8 November 2006, 2006.

Moritz, H., Advanced Physical Geodesy, Herbert Wichmann, Karlsruhe, Germany, 1980.

Nakajima, K., Parallel multilevel iterative linear solvers with unstructured adaptive grids for simulations in Earth science, ACES 2001 Special Issue of Concurrency and Computation: Practice and Experience, 2001.

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

Reubelt, T., M. Gotzelmann, and W. Grafarend, A new CHAMP gravitational field model based on the GIS acceleration approach and two years of kinematic CHAMP data, Paper presented at the Joint CHAMP-GRACE Science Team meeting, GFZ Potsdam, Germany, 5–8 July 2004, 2004.

Rowlands, D. D., R. D. Ray, D. S. Chinn, and F. G. Lemoine, Short-arc analysis of intersatellite tracking data in a gravity mapping mission, J. Geod., 76, 307–316, 2002.

Sanso, F. and C. C. Tscherning, Fast spherical collocation: theory and examples, J. Geod., 77, 101–112, 2003.

Schuh, W.-D., Scientific data processing algorithms, in From Eotvos to milligal, final report, ESA report study, ESA/ESTEC contract 13392/98/LN/GD, edited by Sunkel, H., 2000.

Shum, C. K., S. C. Han, A. Braun, F. Schwartz, and M. Schmidt, Quantification of GRACE continental hydrology observations: aliasing and high-frequency signal recovery, GRACE-Hydrology workshop, UC-Irvine, March 22nd, 2004.

Strykowski, G., Outline of a new space-domain method for forward modelling, in Gravity Field of the Earth, Proceedings of the 1st International Symposium of the International Gravity Field Service (IGFS), 28th August–1st September 2006, Istanbul, Turkey, 2006 (accepted).

<footnote>B. Minchev (e-mail: minchev@math.uni-potsdam.de), A. Chambodut, M. Holschneider, I. Panet, E. Schöll, M. Mandea, and G. Ramillien</footnote>