A methodology for least-squares local quasi-geoid modelling using a noisy satellite-only gravity field model

(1)

A methodology for least-squares local quasi-geoid modelling using a noisy satellite-only

gravity field model

Klees, R.; Slobbe, D. C.; Farahani, H. H. DOI

10.1007/s00190-017-1076-0 Publication date

2018

Document Version Final published version Published in

Journal of Geodesy

Citation (APA)

Klees, R., Slobbe, D. C., & Farahani, H. H. (2018). A methodology for least-squares local quasi-geoid modelling using a noisy satellite-only gravity field model. Journal of Geodesy, 92(4), 431–442.

https://doi.org/10.1007/s00190-017-1076-0 Important note

To cite this publication, please use the final published version (if applicable). Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

(2)

DOI 10.1007/s00190-017-1076-0

O R I G I NA L A RT I C L E

A methodology for least-squares local quasi-geoid modelling using

a noisy satellite-only gravity field model

R. Klees1 · D. C. Slobbe1 · H. H. Farahani1

Received: 20 July 2016 / Accepted: 1 October 2017

Abstract The paper is about a methodology to combine a

noisy satellite-only global gravity field model (GGM) with other noisy datasets to estimate a local quasi-geoid model using weighted least-squares techniques. In this way, we attempt to improve the quality of the estimated quasi-geoid model and to complement it with a full noise covariance matrix for quality control and further data processing. The methodology goes beyond the classical remove–compute– restore approach, which does not account for the noise in the satellite-only GGM. We suggest and analyse three dif-ferent approaches of data combination. Two of them are based on a local single-scale spherical radial basis func-tion (SRBF) model of the disturbing potential, and one is based on a two-scale SRBF model. Using numerical exper-iments, we show that a single-scale SRBF model does not fully exploit the information in the satellite-only GGM. We explain this by a lack of flexibility of a single-scale SRBF model to deal with datasets of significantly different band-widths. The two-scale SRBF model performs well in this respect, provided that the model coefficients representing the two scales are estimated separately. The corresponding methodology is developed in this paper. Using the statistics of the least-squares residuals and the statistics of the errors in the estimated two-scale quasi-geoid model, we demonstrate that the developed methodology provides a two-scale quasi-geoid model, which exploits the information in all datasets.

Keywords Local quasi-geoid modelling· Least-squares

approximation· Spherical radial basis functions · Poisson

B

R. Klees r.klees@tudelft.nl

1 _{Delft University of Technology, Stevinweg 1, 2628 CN Delft,} The Netherlands

wavelets· Noisy global gravity field model · Multi-scale analysis

1 Introduction

In this paper, we investigate the combination of a noisy satellite-only global gravity model (GGM) with noisy high-resolution datasets (e.g. terrestrial gravity anomalies) to estimate a local quasi-geoid model using weighted least-squares techniques. By considering the satellite-only GGM as one of the noisy datasets, we expect to improve the quality of the estimated local quasi-geoid model. By exploiting exist-ing information about the noise variances and covariances in combination with weighted least-squares techniques, we aim at making a step forward towards a comprehensive descrip-tion of the quality of the estimated quasi-geoid model in terms of a full noise covariance matrix for quality control and further data processing.

The problem is timely. The quality and spatial resolu-tion of the most recent satellite-only GGMs, which are mainly based on data of the Gravity Recovery and Climate Experiment (GRACE) and Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite missions, have improved dramatically compared to pre-mission models. Moreover, the GGM’s spherical harmonic coefficients are now complemented with a full noise covariance matrix. For models such as GOCO05s (Mayer-Gürr et al. 2015), the quality of the noise covariance matrix benefits from (i) a post-fit residual analysis, which provides more real-istic models of the data noise (e.g. Farahani et al. 2013), and (ii) numerically efficient algorithms to propagate the full data noise covariance matrices into the noise covariance matrix of the estimated spherical harmonic coefficients (e.g.

(3)

regarding high-resolution datasets frequently used in local quasi-geoid modelling. An example isFarahani et al.(2017) who derived coloured noise models for radar altimeter-based along-track quasi-geoid height differences, or Slobbe (2013) who successfully accounted for long-wavelength errors in terrestrial gravity anomalies (e.g.Heck 1990) by augmenting the functional model with additional parame-ters.

Until now, a GGM with full noise covariance matrix has not been used as one of the noisy datasets in the compu-tation of a local quasi-geoid model. The standard approach is remove–compute–restore approach, i.e. the GGM is used to facilitate a local approach to quasi-geoid modelling by removing most of the energy in the data at the long wave-lengths. From a theoretical point of view, the combination of a local set of terrestrial gravity anomalies with a noisy GGM has been considered in a number of publications since the 1980s. They all use a modification of Stokes’ formula and go back to the early work of Sjöberg and Wenzel (Sjöberg 1980,

1981; Wenzel 1981). At that time, only error degree vari-ances of the GGM were used. The modification of Stokes’ kernel was formulated as a global optimization problem, in which either the variance or the mean square error of the quasi-geoid heights was minimized. Since then, the stochas-tic spectral combination methods were applied routinely when computing local quasi-geoid models. An example is the European quasi-geoid EGG08 (Denker et al. 2008;Denker 2013), which uses the spectral combination method of Wen-zel (1981). Sjöberg (2005) was the first who derived the formalism of a local least-squares modification of Stokes’ formula, which uses, among others, the full noise covariance matrix of the GGM to compute the weights per spherical har-monic degree. In (Sjöberg 2011), the restriction to weights per spherical harmonic degree was given up. To our knowl-edge, this method has not been investigated yet in detail. Some numerical aspects were studied inÅgren(2004) and

Ellmann(2004), in particular numerical instabilities when estimating the spectral weights, which naturally arise in local applications.

Here we follow another approach to local quasi-geoid modelling, which uses least-squares techniques to estimate the parameters of a local spherical radial basis function (SRBF) model of the disturbing potential from the available datasets. Least-squares techniques have many advantages. In particular, they allow us to include noise covariance matrices, improve them using variance component estimation, pro-vide variance–covariance information about the estimated parameters and linear functionals of them, and allow the use of statistical hypothesis testing to test the validity of the mathematical model. Least-squares local quasi-geoid mod-elling using SRBFs has been intensively studied by various authors, see Klees et al. (2008) for a literature overview until 2007. Since then, a number of studies investigated

various aspects of the use of SRBFs in local quasi-geoid modelling. This comprises the choice of the type of SRBFs (e.g.Tenzer and Klees 2008;Bentel et al. 2013a,b), SRBF network design and numerical optimization (e.g. Wittwer 2009), regularization issues (e.g. Naeimi 2013), and the optimization of the location of the SRBFs (e.g. Lin et al. 2014). Some aspects related to the combination of data with different bandwidths have been discussed in Panet et al. (2011); Naeimi(2013); Bentel and Schmidt(2016);

Lieb et al. (2016); Lieb (2017). However, they do not cover numerical studies about the combination of a GGM with full noise covariance matrix with high-resolution noisy datasets.

An alternative to the use of a GGM as one of the noisy datasets in local quasi-geoid modelling is to complement the high-resolution local datasets with the original satellite data at altitude, e.g. satellite gravity gradients from the GOCE mission, low–low satellite-to-satellite tracking data from the GRACE mission, and high–low satellite-to-satellite track-ing data from GRACE, GOCE, and other low Earth orbiters. Examples in the context of local quasi-geoid modelling with SRBFs are Lieb et al.(2015) andLieb (2017). The major drawback of this approach is the complexity of the func-tional model for the low–low satellite-to-satellite tracking data and the huge amount of data. This may be the reason why numerical studies published so far (e.g.Lieb 2017) do not use high–low satellite-to-satellite tracking data, limit to a subset of the available GOCE data (e.g. the second radial derivative of the gravitational potential), or use GRACE-based along-track gravity potential differences as pseudo-observations instead of the original K-band ranging data, often in combi-nation with simplified noise models. Overall, this approach does not offer a significant advantage compared to the use of a GGM that is based on the same data and, therefore, is not pursued in this study.

The study addresses to main research questions: (i) What is a suitable functional model for the satellite-only GGM? (ii) How to combine the satellite-only GGM with high-resolution datasets to obtain a quasi-geoid model that optimally exploits the information content in all datasets?

The remainder of this paper is organized as follows. First, we introduce local single-scale and two-scale SRBF models of the disturbing potential and suggest three func-tional models to be used in a least-squares estimation of the quasi-geoid model from the satellite-only GGM and the high-resolution datasets. Following this, we describe the set-up of the numerical experiments, which were designed to investigate the performance of the SRBF models and the functional models. Thereafter, we present and discuss the results of the numerical experiments. We conclude by empha-sizing the main findings and identifying topics for future research.

(4)

2 Parameterization and functional models

2.1 Local parameterization of the disturbing potential

A prerequisite for local gravity field modelling is that the involved datasets do not contain much long-wavelength energy, where “long” relates to the size of the data area. One possibility to achieve this is by reducing all datasets for the contribution of a global model of the disturbing poten-tial. Then, the disturbing potential to be parameterized in local modelling is a residual quantity with little, though nonzero energy at the long wavelengths. Moreover, there is an upper limit of the highest attainable spatial resolution, which depends on the data distribution and signal-to-noise ratio. This allows us to consider a local model of the dis-turbing potential which is band-limited to some maximum degree, say Lmax. From now on, we call this residual and

band-limited disturbing potential simply “disturbing poten-tial”, denoted T .

In this study, we use SRBFs to model T over an area of interest. Basically, two models will be used. The first one is a single-scale model, i.e.

T(x) = I i=1 ciΦ(x, zi), (1) where Φ(x, zi) = R |x| Lmax l=0 φl|z i| |x| l Ql( ˆx · ˆzi), x∈ extσR, zi ∈ intσR, (2)

is a SRBF located at zi,{ci} are the SRBF coefficients, which

are to be estimated from the data using least-squares tech-niques, Qlis the reproducing kernel of the space of harmonic

functions of degree l,φlis the Legendre coefficient of degree l,ˆx and ˆziare points on the unit sphere, andσRis the surface

of a sphere of radius R. The model of Eq. (1) is referred to as a single-scale model.

Alternatively, we may use a multi-scale model involving several sets of SRBFs representing different scales, i.e.

T(x) = J j=1 Ij i=1 cj iΨj(x, zj i), (3)

where j indicates the scale,{cj·} are the SRBF coefficients

at scale j , andΨj(·, zj i) is a SRBF of scale j centred at the

point zj i. In the context of a multi-scale analysis, the SRBF Ψjmay be defined as (e.g.Lieb 2017)

Ψj(x, zj i) = ⎧ ⎪ ⎨ ⎪ ⎩ R |x| ll1=0φl(1) _|z 1i| |x| l Ql( ˆx · ˆz1i) for j=1 R |x| lj l=0(φl( j)− φl( j−1)) _|z j i| |x| l Ql( ˆx · ˆzj i) for j =2 . . . J. (4) Frequently, the relation lj = 2j−1 is used to relate the scale

index j to the maximum spherical harmonic degree lj, which

is resolved at scale j , though other choices are possible.

2.2 Functional models

In the framework of this study, we assume that there are basically two datasets, i.e. a low-resolution dataset and a high-resolution dataset. The low-resolution dataset d1is

syn-thesized from the spherical harmonic coefficients of the GGM as d1(x1k) = L1 n=0 2n+1 m=1 ˆcnm− c(nmref) (F1Hnm)(x1k), k= 1 . . . K1, (5)

where {ˆcnm} are the spherical harmonic coefficients of the

GGM,{c(nmref)} are the spherical harmonic coefficients of the

reference GGM, and Hnm is a solid spherical harmonic of

degree n. The low-resolution dataset is band-limited to a degree L1 ≤ LGGM, where LGGM is the maximum degree

of the GGM. We assume that the high-resolution dataset {d2(x2k) : k = 1 . . . K2} allows the resolution of

wave-lengths up to a maximum degree L2 ≤ Lmax, where L2

depends on the point density and the signal-to-noise ratio. Defining a kernel δL(x, y) = L n=0 1 4π R2 _R |x| n+1 _R |y| n+1 Qn( ˆx · ˆy), x, y ∈ extσR, (6)

a spherical convolution of T withδLas (δL∗ T )(x) =

σR

δL(x, y)T (y) dσR(y), (7)

and linear functionals F1and F2of the disturbing potential

T , we may relate the datasets d1 and d2 to the disturbing

potential T as E{d1}(x1k) = F1(δL1 ∗ T ) (x1k), k = 1 . . . K1, (8) E{d2}(x2k) = F2(δL2 ∗ T ) (x2k), k = 1 . . . K2, (9)

where E{·} denotes mathematical expectation. We will inves-tigate three functional models to estimate a local quasi-geoid

(5)

model by least squares from the low-resolution dataset d1

and the high-resolution dataset d2.

Functional model no. 1 uses the single-scale model of the

disturbing potential, Eq. (1), and reads:

E{d1}(x1k)= I i=1 ci F1(δL1∗ Φ) (x1k, zi), k = 1 . . . K1, (10) E{d2}(x2k)= I i=1 ci F2(δL2 ∗ Φ) (x2k, zi), k = 1 . . . K2. (11) The coefficients{ci} are estimated simultaneously from the

two noisy datasets using weighted least-squares techniques. The weight matrix of each dataset is the inverse of the noise cofactor matrix.

Functional model no. 2 uses the single-scale model of the

disturbing potential, Eq. (1), and reads

(P ∗ E{d1})(x1k)= I i=1 ci F1(P ∗ Φ) (x1k, zi), k =1 . . . K1, (12) E{d2}(x2k) = I i=1 ci F2(δL2∗ Φ) (x2k, zi), k = 1 . . . K2. (13) The kernel P of Eq. (12) is defined as

P(x, y) = ∞ n=0 1 4π R2 _R |x| n+1 _R |y| n+1 hnQn( ˆx · ˆy), x, y ∈ extσR. (14)

The Legendre coefficients{hn: n = 1, 2, . . .} are equal to 1

for degrees n≤ p1, taper off between degrees p1< n < p2,

and are zero for all degrees n≥ p2. An example is a cosine

taper, hn= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1, n < p1 0.5 + 0.5 cos π n−p1 p2−p1 , p1≤ n ≤ p2≤ L2 0, n > p2 . (15) This taper will be used in the numerical experiments of Sect.3. The coefficients{ci} are estimated simultaneously

from the two noisy datasets using weighted least-squares techniques. The weight matrix of each dataset is propor-tional to the inverse of the noise covariance matrix. The noise covariance matrix of P∗ d1is computed from the full noise

covariance matrix of d1using the law of covariance

propa-gation.

The difference between the functional models no. 1 and no. 2 is in the functional model of the low-resolution dataset. Functional model no. 2 uses a tapered SRBF, whereas functional model no. 1 uses a truncated SRBF. Moreover, functional model no. 2 applies the same taper to the dataset, whereas functional model no. 1 uses the original dataset.

Functional model no. 3 uses a two-scale model of the

dis-turbing potential, i.e. Eq. (3) with J = 2:

T(x) = I1 i=1 c1iΨ1(x, z1i) + I2 i=1 c2iΨ2(x, z2i). (16)

The first term on the right-hand side is a low-resolution model of T comprising degrees from 0 to L1, i.e. its

res-olution is identical to the resres-olution of dataset d1. The

second term on the right-hand side complements the low-resolution model to the maximum low-resolution L2 of dataset

d2. In the context of a multi-resolution analysis, it represents

a detail space comprising wavelengths from degrees L1+ 1

to L2.

The basis functionsΨ1andΨ2of Eq. (16) are defined as

Ψ1(x, z) = (P ∗ Φ)(x, z), (17) Ψ2(x, z) = (δL2− P) ∗ Φ (x, z), (18) withΦ(x, z) of Eq. (2). Inserting the last two equations into Eq. (16), the two-scale model of the disturbing potential T is written as T(x) = I1 i=1 c1i(P ∗ Φ)(x, z1i) + I2 i=1 c2i (δL2 − P) ∗ Φ (x, z2i), (19)

withΦ of Eq. (2). The coefficients{c1i} and {c2i} are

esti-mated in two steps. First, we use the functional model

E{d2}(x2k)= I2 i=1 c2i F2(δL2 ∗ Φ) (x2k, z2i), k =1 . . . K2, (20) and estimate the coefficients {c2i} using weighted least

squares. Then, we define a new dataset

d3(x1k) :=

I1

i₌₁

(6)

where{ˆc2i} denotes the least-squares estimate of {c2i}. The

resolution of the dataset d3is identical to the resolution of

the dataset P∗ d1. In that sense, d3and P∗ d1are spectrally

consistent. Then, we use the functional model (P ∗ E{d1})(x1k) E{d3}(x1k) = I1 i=1 c1i(F1Ψ1)(x1k, z1i) = I1 i=1 c1i (F1(P ∗ Φ)) (x1k, z1i), k = 1 . . . K1, (22)

and compute an estimate{ˆc1i} of the coefficients {c1i}, using

weighted least-squares techniques. The noise covariance matrix of dataset d3is computed from the noise covariance

matrix of the estimated coefficients{ˆc2i} using the law of

covariance propagation. It is a full matrix. The least-squares estimate of the disturbing potential is then given by Eq. (19), with{c1i} and {c2i} replaced by the estimates {ˆc1i} and {ˆc2i},

respectively.

Remarks

1. The motivation of using the functional model of Eq. (10) is the following. Dataset d1and its full noise covariance

matrix are band-limited to a degree L1 ≤ LGGM.

There-fore, the right-hand side of the functional model must also be band-limited to the same degree. To achieve this, we consider the signal 1·Φ, with Φ of Eq. (2). We expand this signal on the sphereσRin spherical harmonics and

truncate the expansion at degree L1. The result is

iden-tical toδL1 ∗ Φ. If the right-hand side of the functional

model would not be band-limited to degree L1, the

least-squares estimate of the coefficients{c1i} would be biased

towards zero for the wavelengths above degree L1. This

is due to the fact that a band-limited noise covariance matrix is equivalent to zero noise and noise correlations for degrees above L1.

2. The motivation to use the functional model of Eq. (12) is the result of numerical experiments which are described in Sect. 3 and discussed in Sect. 4. There, we will show that|d1(·)− _iI₌₁1 ci

F1(δL1∗ Φ)

(·, zi)| is much

larger than the noise in the dataset d1, i.e. the functional

model of Eq. (10) is not accurate enough. Compared to this, the error of the functional model of Eq. (12) can be made much smaller than the data noise standard deviation if P of Eq. (14) is chosen as in Eq. (15).

3. The functional model of Eqs. (20), (22) is different from the model suggested in Lieb(2017), which in our nota-tion is E{d1}(·) = I1 i=1 c1i F1(δL1 ∗ Φ) (·, z1i), (23) E{d2}(·) = I2 i=1 c1i F2(δL1∗ Φ) (·, z1i) + I2 i=1 c2i F2 (δL2− δL1) ∗ Φ (·, z2i). (24) Moreover, Lieb (2017) suggests to estimate the coef-ficients {c1i} and {c2i} simultaneously using weighted

least-squares techniques. Some preliminary experiments with this model andΦ set equal to the Abel–Poisson ker-nel (Freeden et al.1998) point to a sub-optimal quality of the estimated quasi-geoid model at the resolution of the dataset d1, which is likely caused by the simultaneous

estimation of the two sets of coefficients{c1i} and {c2i}.

However, additional numerical experiments are neces-sary to support these preliminary results. They are out of the scope of this study.

3 Numerical experiments

The parameterizations and functional models of Sect.2will be analysed using numerical experiments. Though from a practical point of view, working with real data may be desired, we decide to use a state-of-the-art combined GGM and a satellite-only GGM to generate the (noise-free) high-resolution and low-high-resolution datasets, respectively. The main motivation for us to prefer GGMs to real datasets is that some problems and limitations of the functional models of Sect.2.2would be masked by deficiencies in real datasets, e.g. unmodelled signal and noise and data gaps. This would make a proper interpretation of the results impossible. Gen-erating the exact data from GGMs and adding noise which is consistent with the corresponding noise covariance matrix provides a complete error control and facilitates a proper interpretation of the results.

The datasets are generated from the GGMs using a spherical harmonic synthesis and thereafter reduced for the contribution of a long-wavelength GGM, which serves as the reference model. Here we use EIGEN-6C4 (Förste et al. 2014) to generate the (noise-free) high-resolution dataset and the regularized version of GOCO05s (Mayer-Gürr et al. 2015) to generate the (noise-free) low-resolution dataset. The latter is also used as the reference model, though up to a smaller maximum degree. The low-resolution dataset con-sists of a set of height anomalies. This is a logic choice as the target quantity is a quasi-geoid model. The noise covari-ance matrix of the low-resolution dataset is obtained from the full noise covariance matrix of the spherical harmonic coeffi-cients of the unregularized version of GOCO05s by applying the law of covariance propagation. A logic choice for the

(7)

Table 1 Experimental set-ups used in Sect.4

No. 1 High-resolution dataset Gravity disturbances from EIGEN-6C4 complete to degree 500 - GOCO05s complete to degree 150; Reuter gridReuter(1982), control parameter 875 (width 1221,≈ 23 km along the meridians); zero-mean, white Gaussian, 2 mGal standard deviation

Low-resolution dataset Height anomalies from GOCO05s complete to degree 200 - GOCO05s complete to degree 150; Reuter grid, control parameter 350 (width 3051,≈ 57 km along the meridians); full noise covariance matrix GOCO05s complete to degree 200, propagated using the law of covariance propagation

Data area 44◦–68◦N and 11◦W–15◦E

Parameterization Poisson wavelets of order 3, depth = 30 km, Fibonacci grid with a mean node distance of 25 km, max absolute parameterization error = 4.7 mm, SD of parameterization error = 0.9 mm (over area of interest) No. 2 High-resolution dataset Gravity disturbances from EIGEN-6C4 complete to degree 500 - GOCO05s complete to degree 150;

Reuter grid, control parameter 875 (width 1221,≈ 23 km along the meridians); zero-mean, white Gaussian, 2 mGal standard deviation

Low-resolution dataset Height anomalies from GOCO05s complete to degree 230 - GOCO05s complete to degree 150; Reuter grid, control parameter 402.5 (width 2650,≈ 50 km along the meridians); full noise covariance matrix GOCO05s complete to degree 230, filtered with cosine taper, propagated using the law of covariance propagation

Data area High-resolution dataset: 39◦–73◦N and 16◦W–20◦E; low-resolution dataset: 44◦–68◦N and 11◦W–15◦E Parameterization Φ of Eq. (20): Poisson wavelets of order 3, depth = 30 km, Fibonacci grid with a mean node distance of 25 km; max absolute parameterization error = 4.5 mm, SD of parameterization error = 0.9 mm (over area of interest)

Ψ1of Eq. (22): Poisson wavelets of order 3, depth = 60 km, Fibonacci grid with a mean node distance of 60 km; max absolute parameterization error = 0.01 mm (over area of interest)

The area of interest is 49◦–63◦N and 6◦W–10◦E in all set-ups The parameterization area is identical to the data area

high-resolution dataset would be gravity anomalies. Here, we use gravity disturbances for simplicity reasons. Noise in gravity disturbances is zero-mean white Gaussian, i.e. the noise covariance matrix of the high-resolution dataset is a scaled unit matrix. More details about the datasets are pro-vided in Table1.

The datasets are generated at the Earth’s surface and cover an area, which is referred to as “the data area”. The Earth’s surface is represented by the digital elevation model Euro-DEM v1.0 (Hovenbitzer 2008) with 2grid width. In areas

where this model is not available, we use SRTM version 2.1 (Farr et al. 2007) with 3grid width. For the remaining areas, we use ASTER GDEM v2 (Tachikawa et al. 2011) with 1 grid width.

The SRBF of Eq. (2) is a Poisson wavelet of order 3 (Holschneider and Iglewska-Nowak 2007). Different from the Shannon kernel, which is frequently used in local quasi-geoid modelling, the Legendre spectrum of a Poisson wavelet relative to a sphere of radius R has a peak at degree _R3R_−|z|, where z< R is the location of the Poisson wavelet (cf. Fig1).

(8)

100 101 102 103 104 degree 0 0.2 0.4 0.6 0.8 1

Fig. 1 Normalized Legendre spectrum of Poisson wavelets of order 3

at a depth of (from right to left) 20, 40, 80, 160, and 320 km, respectively. Note the logarithmic scale of the horizontal axis

The Legendre spectrum of a Poisson wavelet may give the impression that a single-scale Poisson wavelet model is not able to accurately represent a quasi-geoid with a resolution one typically encounters in practice. Therefore,Chambodut et al. (2005) suggest to use Poisson wavelets of different scales to guarantee that the space of spherical harmonic complete to a degree L2is sufficiently well covered.

How-ever,Slobbe(2013) successfully used a single-scale Poisson wavelet model to compute a quasi-geoid model for the Netherlands mainland, continental shelf, and Wadden Islands with an accuracy of about 1.5 cm standard deviation using real data. The only prerequisite is that the energy in the data at the lowest and highest frequencies is reduced by using a refer-ence GGM and a digital terrain model, respectively. The Pois-son wavelets are located at a constant depth below the Earth’s surface and cover the data area. Their horizontal positions correspond to the points of a Fibonacci grid (Gonzalez 2010). Whenever a new set of Poisson wavelets is chosen in the numerical experiments, we have to determine the optimal depth and the optimal mean distance between the Poisson wavelets. This is done using noise-free datasets generated from the corresponding GGMs on grids dense enough to pre-serve the information content in the GGM. We define a set of candidate depths and candidate mean distances, and estimate the model coefficients by least squares using the correspond-ing dataset (i.e. gravity disturbances when lookcorrespond-ing for a high-resolution model and height anomalies when looking for a low-resolution model). The depth and mean distance that pro-vide the model with the smallest RMS difference to a height anomaly control dataset are selected. The fit of this model to the control dataset is referred to as the “parameterization error”. Note that the parameterization error is always defined in terms of height anomalies, no matter whether the dataset comprises gravity disturbances or height anomalies. Other control datasets are generated to assess the quality of the estimated quasi-geoid models in Sect.4. They always com-prise height anomalies on grids different from the data grids and are computed using a spherical harmonic synthesis of the GGMs from which the noise-free datasets were generated.

Table 2 Statistics of the least-squares residuals using functional model

no. 1 of Sect.2.2

Dataset Units # Of points Min Max Mean SD

d1 cm 1297 − 34.20 30.75 − 0.22 9.86

d2 mGal 8238 − 9.90 9.42 0.01 2.02

The statistics are computed over the area bounded by 44.5◦–67.5◦N and 11.5◦W–14.5◦E

# Of points refers to the whole dataset

When computing quasi-geoid models using weighted least-squares techniques, we calculate the normal equa-tions explicitly and apply Tikhonov regularization (Tikhonov 1963) with a unit regularization matrix. The regularization parameter is fixed using the method in (Wittwer 2009). The normal equations are solved using a parallelized QR-decomposition with column pivoting. This solver is preferred to a Cholesky decomposition due to its much better stability for ill-conditioned linear systems at the benefit of a smaller bias in the least-squares estimate due to a smaller regulariza-tion parameter.

Table1 summarizes the set-up of the numerical experi-ments, which will be used in Sect.4.

4 Results

4.1 Functional model no. 1

We use experimental set-up no. 1 of Table1. Table2shows some statistics of the least-squares residuals for the estimated quasi-geoid model. The standard deviation (SD) of the resid-uals is 9.86 cm for the low-resolution dataset and 2.02 mGal for the high-resolution dataset. The latter corresponds to the standard deviation of the noise in the high-resolution dataset. The former, however, is much larger than the noise. From this we conclude that the estimated quasi-geoid model fits the high-resolution dataset within noise, but gives a poor fit to the low-resolution dataset.

Table3shows the statistics of the errors in the estimated quasi-geoid model. They are computed over the area of inter-est.

The errors range from−10.68 to 12.76 cm; the error SD is 3.46 cm. After applying a low-pass filter at the cut-off degree 200, the estimated quasi-geoid model error SD is 7.62 and 7.68 cm depending on what low-resolution signal is taken as the reference.

To get more insight into the reason why the estimated quasi-geoid model does not fit the low-resolution dataset within noise, we repeat the experiment with noise-free data. The error SD of the estimated quasi-geoid model reduces from 3.46 to 0.09 cm. This is identical to the SD of the param-eterization error. Hence, when using noise-free data, the

(9)

esti-Table 3 Error statistics (in units of cm) of the quasi-geoid model which

has been estimated using the functional model no. 1 of Sect.2.2

Control dataset # Of points Min Max Mean SD

F1(δ200∗ T1) 1152 − 21.52 21.65 0.07 7.68 F1(δ500∗ T2) 6776 − 10.68 12.76 0.05 3.46 F1(δ200∗ T2) 1152 − 22.31 21.66 0.15 7.62 The statistics are computed over the area of interest, which is bounded by 49◦–63◦N and 6◦W–10◦E, on an equal-angular grid of width 1048 (1st row) and 2329(2nd and 3rd rows)

T1 = GOCO05s minus GOCO05s complete to degree 150; T2 = EIGEN-6C4 minus GOCO05s complete to degree 150. F1 = height anomaly functional

mated quasi-geoid model perfectly fits the high-resolution dataset. This does not apply, however, to the low-pass-filtered quasi-geoid model; the error SD is 7.12 and 7.46 cm, respec-tively, i.e. comparable to the results using noisy data. Hence, the poor fit of the estimated quasi-geoid model to the low-resolution dataset cannot be explained by the noise in this dataset. Additional numerical experiments (not shown here) reveal that the fit to the low-resolution dataset can only be improved by further increasing the size of the data area. The 5◦extension beyond the area of interest in all directions as used here is already a challenge in real quasi-geoid modelling as access to data of neighbouring countries is not guaranteed. Moreover, we found that the fit to the low-resolution dataset improves slowly when enlarging the data area. From this we conclude that the poor fit of the estimated quasi-geoid model to the low-resolution dataset is caused by the hard truncation of the Poisson wavelets. This introduces strong spatial-domain oscillations, which are cut off at the border of the data area when computing the elements of the design matrix. This introduces errors in the functional model, which exceed by far the noise in the low-resolution dataset, as shown in Table3.

We use experimental set-up no. 2 of Table1. The kernel P, which according to Eq. (12) is used in the functional model of the low-resolution dataset, is chosen according to Eqs. (14) and (15). The cosine taper parameters are set equal to p1=

150 and p2= 230. Hence, the filtered low-resolution dataset

P∗ d1of Eq. (12) is band-limited to degree L1= 230.

The choice of p1 and p2 is a trade-off between loss of

information in the low-resolution dataset by filtering (i.e. nonzero d1 − P ∗ d1), and a reduction in the area under

the side lobes of the cosine taper, which cause oscillations of the filtered Poisson wavelets extending beyond the data area. Moreover, the difference p2− p1determines how fast

the oscillations roll off. The difference p2− p1 = 80 has

been fixed after some numerical experiments. Note that the maximum degree of the reference GGM (which is Lref= 150

in our experiments) and the maximum degree of the noisy

Table 4 Statistics of the least-squares residuals using functional model

no. 2 of Sect.2.2

P∗ d1 cm 1751 − 15.45 13.29 0.00 3.42

d2 mGal 8238 − 7.90 8.52 0.00 1.94

GGM (which is LGGM = 280 for GOCO05s) impose lower

and upper bounds, respectively, on the choice of p1and p2,

i.e. p1≥ 150 and p2≤ 280.

The maximum possible value of p2is equal to the

max-imum degree of the GGM, LGGM. The GGM used in this

study is GOCO05s. However, for GOCO05s the cumulative height anomaly commission error increases exponentially with increasing degree. It is 1.5 cm at degree 200, but already 3.6 cm at degree 230, and 6.8 cm at degree 250. The height anomaly signal and noise degree variances intersect at degree 257. Hence, when assuming that the noise standard devia-tion in the high-resoludevia-tion dataset does not exceed 1–2 mGal (which applies to good terrestrial gravity anomaly datasets), it does not make sense to use a low-resolution dataset com-plete to the maximum possible degree of LGGM = 280.

Another reason in favour of a choice L1 < LGGMis the fact

that the condition number of the noise covariance matrix of the low-resolution dataset, which is propagated from the full noise covariance matrix of the spherical harmonic coeffi-cients of the GGM, increases with increasing L1. This makes

the computation of the least-squares estimator numerically challenging.

Table4shows the statistics of the least-squares residuals for the estimated quasi-geoid model. The SD is 1.94 mGal for the high-resolution dataset. This is close to the SD of the superimposed zero-mean white Gaussian noise of 2.0 mGal. From this we conclude that the model fit is within noise. The situation is different for the low-resolution dataset. The SD of the residuals is 3.42 cm. This is a factor of 2 larger than the average SD of the data noise. Obviously, for some reason, the estimated quasi-geoid model does not fit the low-resolution dataset as one may expect.

Table5shows the statistics of the errors in the estimated quasi-geoid model. They are computed over the area of inter-est.

The error SD is 3.01 cm. The error SD of the cosine-tapered quasi-geoid model is 1.89 and 1.76 cm, respectively, depending on what control dataset is used. We may compare this with a quasi-geoid model which is estimated using only the high-resolution dataset. The corresponding error SDs are 3.01, 1.58, and 1.71 cm, respectively. Hence, adding the low-resolution dataset does not improve the accuracy of the estimated quasi-geoid model at low frequencies (i.e. below

(10)

Table 5 Error statistics (in units of cm) of the quasi-geoid model which

has been estimated using the functional model no. 2 of Sect.2.2with cosine taper parameters p1= 150 and p2= 230

Control dataset # Of points Min Max Mean SD F1(P ∗ T1) 1400 − 5.60 5.40 0.00 1.89 F1(δ500∗ T2) 6776 − 11.28 11.33 0.01 3.01 F1(P ∗ T2) 1400 − 4.90 5.89 0.01 1.76 The statistics are computed over the area of interest, which is bounded by 49◦–63◦N and 6◦W–10◦E, on an equal-angular grid of width 1048 (1st row) and 2329(2nd and 3rd rows)

degree 230). This is unexpected, because the low-resolution dataset has a high quality and should improve the estimated quasi-geoid model at the low frequencies. This, and the large least-squares residuals of the dataset P∗ d1, which by far

exceed the noise, implies that for some reason, the informa-tion content in the low-resoluinforma-tion dataset is not fully exploited in the combined least-squares adjustment.

To help understand this result, we run two additional experiments. First of all, we compute the error of the functional model of the (noise-free) low-resolution dataset, Eq. (12), i.e.(P ∗ d1)(·) −

I

i=1ci (F1(P ∗ Φ)) (·, zi). We

found that it does not exceed 0.01 cm over the area of inter-est. In a second experiment, we compute a low-resolution quasi-geoid model using the noisy low-resolution dataset and the functional model of Eq. (12). The SD of the resid-uals is 0.13 cm. This is significantly smaller than for the solution which uses both datasets (SD = 3.42 cm, cf. 1st row in Table4). Moreover, the error SD of the low-resolution quasi-geoid model is 1.58 cm when evaluated over the area of interest. This is also smaller than the error we obtain when using both datasets (SD = 3.89 cm, cf. 1st row in Table5).

Our interpretation of the results of these experiments is that a single-scale model is not able to fit two datasets of sig-nificantly different bandwidths. Consequently, the weighted least-squares principle forces the solution to match the high-resolution dataset (because it comprises many more observations than the low-resolution dataset) at the price of a larger mismatch to the low-resolution dataset.

To support this interpretation, we choose experimental set-up no. 3 of Table1, which is similar to experimental set-up no. 2, but involves a high-resolution dataset and a low-resolution dataset with a much larger bandwidth difference of 335% compared to 117% of experimental set-up no. 2. The SD of the least-squares residuals of the low-resolution dataset increases from 3.42 cm (cf. Table4, 2nd row) to 5.13 cm, whereas the SD of the least-squares residuals of the high-resolution dataset does not change. Hence, when increasing the bandwidth difference between the high-resolution and the resolution datasets, the fit of the model to the

low-Table 6 Statistics of the least-squares residuals of the datasets used

to estimate the two-scale SRBF model using the functional models of Eqs. (20), (22)

d2 mGal 15967 − 6.01 7.30 − 0.01 1.74

P∗ d1 cm 1751 − 4.85 4.52 − 0.02 1.47

d3 cm 1751 − 8.04 6.90 − 0.33 2.69

resolution dataset becomes worse. This provides evidence that our interpretation is correct.

We use experimental set-up no. 4 of Table1. Note that this set-up uses a high-resolution dataset which extends over a larger area than the experimental set-ups no. 1–3. The rea-son is the following. When using the functional model of Eq. (22) to estimate a low-resolution quasi-geoid model, the dataset d3must be available over the data area, which in all

experiments extends by 5◦in all directions beyond the area of interest. If we would use the same data area for dataset d2

when estimating the coefficients{c21} using the functional

model of Eq. (20), the dataset d3 of Eq. (21) would suffer

from edge effects. Therefore, the high-resolution dataset d2

must be available over an area, which is larger than the data area of the dataset d3. The additional extension must be

cho-sen to reduce the edge effects in dataset d3below the noise

level. In the experimental set-up no. 4, we use an additional extension by 5◦in all directions. Test computations reveal that this choice causes edge effects in dataset d3, which are

negligible compared to the noise. We expect that the addi-tional extension can be chosen much smaller. To find the minimum extension may be the subject of another study.

Table6shows the statistics of the least-squares residuals for the model of Eqs. (20) and (22), respectively. The fit of the high-resolution dataset d2to the model of Eq. (20) has a

SD of 1.74 mGal. This is close to the superimposed noise of SD = 2 mGal; a similar fit has also been observed when using the functional model no. 2 (cf. Table4). However, compared to the functional model no. 2, the fit of the dataset P ∗ d1

to the model has improved dramatically: from 3.42 cm (1st row of Table4) to 1.47 cm (2nd row of Table6). A value of 1.47 cm is consistent with the noise in the dataset P∗d1. The

fit of the dataset d3to the model is 2.69 cm, i.e. dataset P∗d1

has a larger contribution to the model compared to dataset d3.

(Both datasets are evenly large.) This is also consistent with the expectations based on an analysis of the noise covariance matrices of the two datasets (not shown here).

The two-scale model appears to have a much higher qual-ity than the single-scale model of Sect.4.2. This follows from

(11)

Table 7 Error statistics (in units of centimetres) of the quasi-geoid

model which has been estimated using the functional model no. 3 of Sect.4.3

Control dataset # Of points Min Max Mean SD F1(P ∗ T1) 1400 − 2.45 2.73 0.00 0.86 F1(δ500∗ T2) 6776 − 6.12 6.65 0.02 1.97 F1(P ∗ T2) 6776 − 2.27 2.42 0.01 0.77 The low-pass filter is a cosine taper with parameters p1 = 150 and p2= 230. The statistics are computed over the area of interest, which is bounded by 49◦–63◦N and 6◦W–10◦E, on an equal-angular grid of width 1048(1st row) and 2329(2nd and 3rd rows)

the statistics of the differences at the control datasets, which are shown in Table7. For instance, the fit of the two-scale model to the control dataset F1(δ500 ∗ T2) improves from

SD = 3.01 cm (2nd row of Table5) to SD = 1.97 cm (2nd row of Table7). The fit to the low-resolution control data improves dramatically, too: from 1.89 cm (F1(P ∗ T1), 1st

row of Table5) and 1.76 cm (F1(P ∗ T2), 3rd row of Table5)

to 0.86 and 0.77 cm, respectively. From this we conclude that the two-scale model in combination with the functional model of Eqs. (20), (22) performs better at all wavelengths than any of the two single-scale models. The improvement is a factor of 2.2 for the wavelengths common to the high-and the low-resolution dataset, high-and a factor of 1.5 for the wavelength not resolved by the low-resolution dataset. The former is due to the fact that the suggested approach which uses the two-scale model fully exploits the higher accuracy of the low-resolution dataset, which is not the case if any of the single-scale models is used.

5 Summary and conclusions

In this study, we investigated different approaches to esti-mate a local SRBF model of the disturbing potential using weighted least squares from a high-resolution dataset and a low-resolution dataset. In practice, the low-resolution dataset represents a satellite-only spherical harmonic model of the global gravity field equipped with a full noise covari-ance matrix. Considering the latter as one of the noisy datasets in local quasi-geoid modelling is considered as a significant improvement to the traditional remove–compute– restore approach. It improves the quality of the estimated quasi-geoid model and paves the way to a complete quality description of the estimated quasi-geoid model in terms of a full noise covariance matrix.

Two approaches investigated in this study use a single-scale SRBF model, but differ in the functional model for the low-resolution dataset. The third one uses a two-scale

SRBF wavelet model and estimates the coefficients per scale independently of each other.

We showed that the functional model of the low-resolution dataset has to be chosen with care. A hard truncation of the SRBFs at the maximum degree of the low-resolution dataset is the right choice in global quasi-geoid modelling, but provides a wrong functional model in local quasi-geoid modelling. This is in line with the results in (Slobbe et al. 2012). Applying a taper to both the low-resolution dataset and the SRBF model solves this problem.

We also showed that a single-scale SRBF model cannot deal with datasets of different bandwidths. The estimated quasi-geoid model is biased towards the high-resolution dataset at the cost of a poor fit to the low-resolution dataset. The latter appeared to be much worse than the noise in this dataset suggested, which indicates that the information con-tent of the low-resolution dataset is not fully exploited.

We suggested the use of a two-scale SRBF model in com-bination with a sequential estimation of the scale-dependent coefficients. The latter differs from what has been suggested in the literature in the context of a multi-scale analysis. In this way, we ensure that the two datasets are weighted in line with their accuracy, the information content in the low-resolution and high-resolution datasets is fully exploited, and the mis-fit of the estimated quasi-geoid model is consistent with the noise in the datasets.

A challenge of the suggested approach in applications involving real datasets is the additional extension of the data area for the high-resolution dataset. In this study, a safe choice has been made to make edge effects insignificant. In applica-tions involving real datasets, access to high-quality terrestrial gravity anomaly datasets of neighbouring countries is not guaranteed. How much the data area needs to be extended and whether data with reduced accuracy can be used in the additional area without introducing distortions in the esti-mated quasi-geoid model has to be investigated.

It would be interesting to compare the two-scale approach suggested in this study with a multi-scale approach, which estimates the coefficients at the two scales simultaneously as suggested in Lieb (2017). Some preliminary experi-ments (not shown here) indicate that such a multi-scale approach provides a sub-optimal low-resolution quasi-geoid model compared to a sequential estimation as suggested here. Whether this may be corrected for by a further optimization of the multi-scale approach, for example, by introducing constraints between the model coefficients asso-ciated with different scales may be the subject of a future study.

Acknowledgements This study was performed in the framework of the

Netherlands Vertical Reference Frame (NEVREF) project, funded by the Netherlands Technology Foundation STW. This support is gratefully acknowledged.

(12)

Open Access This article is distributed under the terms of the Creative

Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

References

Ågren J (2004) Local geoid determination methods for the era of satellite gravimetry. Doctoral Dissertation, Royal Institute of Technology, Stockholm, p 246

Bentel K, Schmidt M, Gerlach C (2013a) Different radial basis func-tions and their applicability for local gravity field representation on the sphere. Int J Geomath 4:67–96.https://doi.org/10.1007/ s13137-012-0046-1

Bentel K, Schmidt M, Denby CR (2013) Artifacts in local gravity representations with spherical radial basis functions. J Geod Sci 3:173–187.https://doi.org/10.2478/Jogs-2013-0029

Bentel K, Schmidt, M (2016) Combining different types of gravity observations in regional gravity modeling in spherical radial basis functions. In: Sneeuw N, Novak P, Crespi M, Sanso F (eds) VIII Hotine–Marussi symposium on mathematical geodesy, IAG sym-posia 142, pp 115–120,https://doi.org/10.1007/1345_2015_2 Chambodut A, Panet I, Mandea M, Diament M, Holschneider M, Jamet

O (2005) Wavelet frames: an alternative to spherical harmonic representation of potential fields. Geophys J Int 163:875–899 Denker H, Barriot JP, Barzaghi R, Fairhead D, Forsberg R, Ihde J,

Keny-eres A, Marti U, Sarrailh M, Tziavos I (2008) The development of the European gravimetric geoid model EGG07. In: Sans F, Sideris MG (eds) Observing our changing earth, international associa-tion of geodesy symposia, vol 133. Springer, Berlin, pp 177–185. https://doi.org/10.1007/978-3-540-85426-5_21

Denker H (2013) Regional gravity field modeling: theory and practical results. In: Xu G (ed) Sciences of geodesy—II. Springer, Berlin, pp 185–291.https://doi.org/10.1007/978-3-642-28000-9 Ellmann A (2004) The geoid for the Baltic countries determined by the

least squares modification of Stokes’ formula. Doctoral Disserta-tion, Royal Institute of Technology, Stockholm, p 82

Farahani HH, Ditmar P, Klees R, Liu X, Zhao Q, Guo J (2013) The static gravity field model DGM-1S from GRACE and GOCE data: computation, validation and an analysis of GOCE mis-sion’s added value. J Geod 87:843–867.https://doi.org/10.1007/ s00190-013-0650-3

Farahani HH, Slobbe DC , Klees R , Kurt Seitz (2017) Impact of accounting for coloured noise in radar altimetry data on a regional quasi-geoid model. J Geod 91: 97–122

Farr TG, Rosen PA, Caro E, Crippen R, Duren R, Hensely S, Kobrick M, Paller M, Rodrigues E, Roth L, Seal D, Shaffer S, Shimada J, Umland J, Werner M, Oskin M, Burbank D, Alsdorf D (2007) The shuttle radar topography mission. Rev Geophys 45(2):33.https:// doi.org/10.1029/2005RG000183

Förste C, Bruinsma SL, Abrikosov O, Lemoine JM, Marty JC, Flechtner F, Balmino G, Barthelmes F, Biancale R (2014) EIGEN-6C4 the latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse. GFZ Data Serv.https://doi.org/10.5880/icgem.2015.1

Freeden W, Gervens T, Schreiner M (1998) Constructive approxima-tions on the sphere with applicaapproxima-tions to geomathematics. Oxford Science Publications, Clarendon Press, Oxford

Gonzalez A (2010) Measurement of areas on a sphere using Fibonacci and latitude-longitude lattices. Math Geosci 42:49–64.https://doi. org/10.1007/s11004-009-9257-x

Heck B (1990) An evaluation of some systematic error sources affecting terrestrial gravity anomalies. Bull Geod 64:88–108.https://doi. org/10.1007/BF02530617

Holschneider M, Iglewska-Nowak I (2007) Poisson wavelets on the sphere. J Fourier Anal Appl 13:405–419

Hovenbitzer M (2008) The European DEM (EuroDEM). Int Arch Pho-togramm Remote Sens Spat Inf Sci 37:1853–1856

Klees R, Tenzer R, Prutkin I, Wittwer T (2008) A data-driven approach to local gravity field modelling using spherical radial basis functions. J Geod 82:457–471. https://doi.org/10.1007/ s00190-007-0196-3

Lieb V, Bouman J, Dettmering D, Fuchs M, Schmidt M (2015) Com-bination of GOCE gravity gradients in regional gravity field modelling using radial basis functions. In: Sneeuw N, Novk P, Crespi M, Sans F (Eds) VIII Hotine–Marussi symposium on math-ematical geodesy, IAG symposia, vol 142, pp 101–108,https://doi. org/10.1007/1345_2015_71

Lieb V, Schmidt M, Dettmering D, Börger K (2016) Combination of various observation techniques for local modeling of the grav-ity field. J Geophys Res Solid Earth. https://doi.org/10.1002/ 2015JB012586

Lieb V (2017) Enhanced regional gravity field modeling from the com-bination of real data via MRR. Ph.D. thesis, DGK series C, no 795, München

Lin M, Denker H, Müller J (2014) Local gravity field modelling using free-positioned point masses. Stud Geophys Geod 58:207–226. https://doi.org/10.1007/s11200-013-1145-7

Mayer-Gürr T, Kvas A, Klinger B, Maier A, (2015) The new combined satellite only model GOCO05s. EGU General Assembly, (2015) Vienna. Austria.https://doi.org/10.13140/RG.2.1.4688.6807 Naeimi M (2013) Inversion of satellite gravity data using spherical radial

base functions. Doctoral Dissertation, Leibniz University Han-nover, Deutsche Geodätische Kommission, Reihe C, Heft Nr. 711, p 130

Panet I, Kuroishi Y, Holschneider M (2011) Wavelet modelling of the gravity field by domain decomposition methods: an example over Japan. Geophys J Int 184:203–219.https://doi.org/10.1111/ j.1365-246X.2010.04840.x

Reuter R (1982) Über Integralformeln der Einheitssphäre und harmonis-che Splinefunktionen. Veröff Geod Inst RWTH Aaharmonis-chen 33:1982 Sjöberg LE (1980) Least squares combination of satellite harmonics

and integral formulas in physical geodesy. Gerlands Beiträge zur Geophysik 89:371–377

Sjöberg LE (1981) Least squares combination of satellite and terrestrial data in physical geodesy. Ann Geophys 37:25–30

Sjöberg LE (2005) A local least-squares modification of Stokes’ for-mula. Stud Geophys Geod 49:23–30

Sjöberg LE (2011) Local least squares spectral filtering and combina-tion by harmonic funccombina-tions on the sphere. J Geod Sci 1:355–360. https://doi.org/10.2478/v10156-011-0015-x

Slobbe DC (2013) Roadmap to a mutually consistent set of offshore vertical reference frames. Ph.D. thesis, Delft Uni-versity of Technology, p 233. https://doi.org/10.4233/uuid: 68e4e599-51ab-40df-8fb4-918f5f54a453

Slobbe DC, Simons FJ, Klees R (2012) The spherical Slepian basis as a means to obtain spectral consistency between mean sea level and the geoid. J Geod 86:609–628.https://doi.org/10.1007/ s00190-012-0543-x

Tachikawa T, Hato M, Kaku M, Iwasaki A (2011) Characteristics of ASTER GDEM version 2. In: IEEE international geoscience and remote sensing symposium IGARSS (2011) Vancouver. BC, Canada,https://doi.org/10.1109/IGARSS.2011.6050017 Tenzer R, Klees R (2008) The choice of the spherical radial basis

functions in local gravity field modelling. Stud Geophys Geod 52:287–304

(13)

Tikhonov AN (1963) Solution of incorrectly formulated problems and the regularization method. Dokl Akad Nauk SSSR 151: 501–504 = Soviet Math Dokl 4: 1035–1038

Wenzel HG (1981) Zur Geoidbestimmung durch Kombination von Schwereanomalien und einem Kugelfunktionsmodell mit Hilfe von Integralformeln. Zeitschrift für Vermessungswesen 106:102– 111

Wittwer T (2009) local gravity field modelling with radial basis func-tions. Doctoral Dissertation, Delft University of Technology, Delft, The Netherlands, p 191