- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications
- Imprint
- Data protection

Journal cover
Journal topic
**Geoscientific Model Development**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications
- Imprint
- Data protection

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications
- Imprint
- Data protection

**Development and technical paper**
03 Dec 2018

**Development and technical paper** | 03 Dec 2018

A multilayer approach and its application to model a local gravimetric quasi-geoid model over the North Sea: QGNSea V1.0

^{1}School of Earth Sciences and Engineering, Hohai University, Nanjing, China^{2}State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan, China^{3}MOE Key Laboratory of Fundamental Physical Quantities Measurement, School of Physics, Huazhong University of Science and Technology, Wuhan, China^{4}School of Geodesy and Geomatics, Wuhan University, Wuhan, China^{5}School of Civil and Transportation Engineering, Guangdong University of Technology, Guangdong, China

^{1}School of Earth Sciences and Engineering, Hohai University, Nanjing, China^{2}State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan, China^{3}MOE Key Laboratory of Fundamental Physical Quantities Measurement, School of Physics, Huazhong University of Science and Technology, Wuhan, China^{4}School of Geodesy and Geomatics, Wuhan University, Wuhan, China^{5}School of Civil and Transportation Engineering, Guangdong University of Technology, Guangdong, China

**Correspondence**: Yihao Wu (yihaowu@hhu.edu.cn) and Chuang Xu (chuangxu@hust.edu.cn)

**Correspondence**: Yihao Wu (yihaowu@hhu.edu.cn) and Chuang Xu (chuangxu@hust.edu.cn)

Abstract

Back to toptop
A multilayer approach is set up for local gravity field recovery within the framework of multi-resolution representation, where the gravity field is parameterized as the superposition of multiple layers of Poisson wavelets located at different depths beneath the Earth's surface. The layers are designed to recover gravity signals at different scales, where the shallow and deep layers mainly capture the short- and long-wavelength signals, respectively. The depths of these layers are linked to the locations of different anomaly sources beneath the Earth's surface, which are estimated by wavelet decomposition and power spectrum analysis. For testing the performance of this approach, a gravimetric quasi-geoid model over the North Sea, QGNSea V1.0, is modeled and validated against independent control data. The results show that the multilayer approach fits the gravity data better than the traditional single-layer approach, particularly in regions with topographical variation. An Akaike information criterion (AIC) test shows that the multilayer model obtains a smaller AIC value and achieves a better balance between the goodness of fit of data and the simplicity of the model. Further, an evaluation using independent GPS/leveling data tests the ability of regional models computed from different approaches towards realistic extrapolation, which shows that the accuracies of the QGNSea V1.0 derived from the multilayer approach are better by 0.4, 0.9, and 1.1 cm in the Netherlands, Belgium, and parts of Germany, respectively, than that using the single-layer approach. Further validation with existing models shows that QGNSea V1.0 is superior with respect to performance and may be beneficial for studying ocean circulation between the North Sea and its neighboring waters.

How to cite

Back to top
top
How to cite.

Wu, Y., Luo, Z., Zhong, B., and Xu, C.: A multilayer approach and its application to model a local gravimetric quasi-geoid model over the North Sea: QGNSea V1.0, Geosci. Model Dev., 11, 4797–4815, https://doi.org/10.5194/gmd-11-4797-2018, 2018.

1 Introduction

Back to toptop
Knowledge of the Earth's gravity field at the regional scale is crucial for a variety of applications in geodesy. It not only facilitates the use of the Global Satellite Navigation System to determine orthometric/normal heights in geodesy and surveying engineering but also plays a fundamental role in oceanography and geophysics.

Regional gravity field determination is typically conducted within the framework of the remove–compute–restore (RCR) methodology (Sjöberg, 2005), where long-wavelength signals are often recovered by satellite-only global geopotential models (GGMs) derived from dedicated satellite gravity missions such as the Gravity Field and Climate Experiment (GRACE; Tapley et al., 2004) and Gravity Field and Steady-State Ocean Circulation Explorer (GOCE; Rummel et al., 2002). Middle- and short-wavelength signals are extracted from locally distributed gravity-related measurements (Guo et al., 2010; Wang et al., 2012, 2018). Spherical radial basis functions (SRBFs) have become of great interest for gravity field modeling at the regional scale in recent years (Eicker et al., 2013; Naeimi et al., 2015). Typically, the most commonly SRBFs are implemented using the single-layer approach; i.e., the parameterization of the gravity field is based on only a single-layer of the SRBF grid (Wittwer, 2009; Bentel et al., 2013; Slobbe et al., 2014; Wu and Luo, 2016).

It has been suspected for long that the single-layer approach may fail to extract the full information contained in local gravity data; thus, the multi-resolution representation (MRR) method with SRBFs has been investigated in recent years (Freeden et al., 1998; Fengler et al., 2004, 2007). Freeden and Schreiner (2006) proposed a multiscale approach based on locally supported wavelets for determining regional geoid undulations from deflections of the vertical. Further, Freeden et al. (2009) demonstrated that a multiscale approach using spherical wavelets provided a powerful technique for the investigation of local fine-structured features such as those caused by plumes, which allowed the scale- and space-dependent characterization of this geophysical phenomenon. Schmidt et al. (2005, 2006, 2007) developed a multi-representational method for static and spatiotemporal gravitational field modeling using SRBFs, where the input gravity signals were decomposed into a number of frequency-dependent detail signals; they concluded that this approach could improve the spanning fixed time intervals with respect to the usual time-variable gravity fields. Chambodut et al. (2005) set up a multiscale method for magnetic and gravity field recovery using Poisson wavelets and created a set of hierarchical meshes associated with the wavelets at different scales, where a level of subdivision corresponded to a given wavelet scale. Panet et al. (2011) extended the approach developed by Chambodut et al. (2005) by applying a domain decomposition approach to defining the hierarchical subdomains of wavelets at different scales; this enabled the splitting of a large problem into smaller ones. The results of these studies show that the multiscale approach using SRBFs has a good potential for gravity field recovery. However, to the best of our knowledge, no direct comparisons have been made between the single-layer approach and the multiscale one regarding their performance in local gravity field recovery. Further, the existing multiscale methods mainly construct the multiscale framework in a mathematical sense, and no explicit geophysical meanings are investigated. In this study, inspired by the power spectral analysis of local gravity signals, we develop new parameterizations of the SRBF network within the MRR approach. In this approach, multiple layers are linked to the anomaly sources at different depths beneath the Earth's surface, and the aim is to recover the signals with different spectral contents. Moreover, the performance of the multilayer approach and traditional single-layer approach is directly compared, and the advantages and disadvantages of the two methods are analyzed.

The structure of the paper is as follows. The study area and data collection methods are described in Sect. 2. Then, the MRR method with SRBFs is introduced, where Poisson wavelets that have band-limited properties are chosen as the basis functions. Wavelet decomposition and power spectrum analysis are applied in constructing the network of Poisson wavelets. In addition, the function model based on this multilayer approach is derived and the method for estimating the unknown coefficients of Poisson wavelets is introduced. The construction of the multilayer model is described in Sect. 3. The performance of the two approaches (single-layer and multilayer) is also compared in this section. Finally, a gravimetric quasi-geoid over the North Sea, called QGNSea V1.0, is modeled using the multilayer approach and compared with other models for cross validation. We present the summary and the main conclusions of this study in Sect. 4.

2 Data and methods

Back to toptop
A region in Europe, from 49 to 61^{∘} N and −6 to 10^{∘} E,
covering the mainland of the Netherlands, Belgium, parts of the North Sea,
the UK, Germany, and France, is chosen as a case study. Data regarding point-wise
terrestrial and shipborne gravity anomalies are used in this study, which
were provided by different institutions. The details of the data
pre-processing procedures can be found in Wu et al. (2017c), where crossover
adjustment and low-pass filters were applied to remove systematic errors and
reduce high-frequency noise, respectively. Since the terrestrial and
shipborne gravity data were derived from different institutions over various
time spans, the horizontal and vertical reference systems need to be
unified. The
European Terrestrial Reference System 1989 (ETRS89) and European Vertical
Reference Frame 2007 (EVRF2007) are chosen as the horizontal and vertical
systems, respectively (Slobbe, 2013). Datum transformations are performed on
all of the data following the methods proposed by Wu et al. (2017c).
Moreover, the long-wavelength signal content in the data is reduced by
removing the contribution of the GOCO05s global geopotential model completely
to degree and order (d/o) of 280 (Mayer-Gürr et al., 2015). At the very
short wavelengths, a residual terrain model (RTM) is applied. The details of
the RTM reduction process and the residual gravity data can be found in Wu et
al. (2017c).

According to Schmidt et al. (2006, 2007), the MRR of the Earth's disturbing
potential *T*(z) at position *z* is expressed as

$$\begin{array}{}\text{(1)}& T\left(z\right)=\stackrel{\mathrm{\u203e}}{T}\left(z\right)+\sum _{i=\mathrm{1}}^{I}{t}_{i}\left(z\right)+\mathit{\delta}\left(z\right),\end{array}$$

where $\stackrel{\mathrm{\u203e}}{T}\left(z\right)$ represents a reference model,
e.g., a GGM computed from spherical harmonics; *δ*(z) represents unmodeled signals; *I* is the number of
levels (resolutions); *t*_{i}(z) is the detailed signal of
level *i*, and the higher the level value *i* is, the finer are the
structures extractable from the input data; *t*_{i}(z) is
computed as a linear combination of SRBFs (Schmidt et al., 2007).

$$\begin{array}{}\text{(2)}& {t}_{i}\left(z\right)=\sum _{n=\mathrm{1}}^{{N}_{i}}{\mathit{\beta}}_{i,n}{\mathrm{\Psi}}_{i}\left(z,{y}_{i,n}\right),\end{array}$$

where Ψ(*z*,*y*) is the SRBF, *N*_{i} and *β*_{i,k}
are the number and unknown coefficient of the SRBF at level *i*,
respectively, and *y*_{i,n} is the position of the SRBF at this level.

The reference GGM and RTM corrections are removed from the original data to
decrease the signal correlation length and smooth the data (Omang and
Forsberg, 2000). Then, only the residual disturbing potential *T*_{res}(z) is parameterized by the SRBF using the MRR approach.
Ignoring the unmodeled signals, the residual disturbing potential is
expressed as a series of detailed signals at different levels, combining
Eqs. (1) and (2):

$$\begin{array}{}\text{(3)}& {T}_{\mathrm{res}}\left(z\right)=\sum _{i=\mathrm{1}}^{I}\sum _{n=\mathrm{1}}^{{N}_{i}}{\mathit{\beta}}_{i,n}{\mathrm{\Psi}}_{i}\left(z,{y}_{i,n}\right),\end{array}$$

where Ψ_{i} is computed as the difference between the spherical scaling
functions with low-pass filter characteristics corresponding to consecutive
levels *i*+1 and *i*; Ψ_{i} can also be expressed as the SRBF with
band-limited properties in the frequency domain (Schmidt et al., 2007).
Ψ represents Poisson wavelets with band-limited properties (Chambodut et
al., 2005) in this study, the complete definition of which can be found in
Holschneider and Iglewska-Nowak (2007).

Poisson wavelets can also be identified as the multipoles inside the Earth,
and the scales of Poisson wavelets can be linked to their depths beneath the
Earth's surface. These depths are the key parameters in determining wavelet
properties in space and frequency domains (Chambodut et al., 2005). The
detailed signal at level *i* in Eq. (2) can be estimated using a linear
combination of Poisson wavelets located at a specific depth. Poisson wavelets
at various depths demonstrate different properties in the frequency domain.
At shallow depths, the scales decrease, and wavelet spectrums shift toward
the high degrees of the spherical harmonics (SH) and become more sensitive to
local signal features with high-frequency properties, and vice versa
(Chambodut et al., 2005). Moreover, Poisson wavelets at different depths can
be linked to the detailed signals at various levels and are sensitive to
various spectral contents of input signals. They can be used for
multi-resolution representation. These properties are crucial for local
gravity field modeling. The residual disturbing potential is typically a
band-limited signal within the RCR framework, and Poisson wavelets with
band-pass filter characteristics are preferable for band-limited signal
recovery (Bentel et al., 2013).

Rather than as an MRR, we interpret Eq. (3) as the multilayer approach that takes into consideration that Poisson wavelets at different depths have different characteristics, and different layers correspond to Poisson wavelets' grids at various depths. We place the Poisson wavelets in the Fibonacci grids under the Earth's surface and keep these grids parallel with the Earth's surface (Tenzer et al., 2012). Instead of associating the Poisson wavelets at different depths with the hierarchical meshes with various levels (Chambodut et al., 2005), we apply a wavelet analysis approach to estimate the depths of multiple layers. This approach is inspired by the power spectrum analysis of the local gravity signals, which shows that the gravity signals are superpositions of contributions from the anomaly sources at different depths, and the signals originating from different anomaly sources have heterogeneous spectral contents (Spector and Grant, 1970; Syberg, 1972; Xu et al., 2018). Since the Poisson wavelets at different depths are sensitive to signals with heterogeneous frequency characteristics, we place Poisson wavelet grids at locations where anomaly sources are situated. In this manner, the contributions of the anomaly sources at various depths can be estimated.

In order to separate the contributions of different anomaly sources, the
wavelet multiscale analysis is applied to decompose the gravity data Δ*g*(*φ*,*λ*) into wavelet approximation *A*_{W}(*φ*,*λ*) and a number of wavelet details ${D}_{\mathrm{w}}(\mathit{\phi},\mathit{\lambda})(w=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots},W)$ at different scales (Jiang et al., 2012; Audet, 2014;
Xu et al., 2017).

$$\begin{array}{}\text{(4)}& \mathrm{\Delta}g(\mathit{\phi},\mathit{\lambda})={A}_{\mathrm{W}}(\mathit{\phi},\mathit{\lambda})+\sum _{w=\mathrm{1}}^{W}{D}_{\mathrm{w}}(\mathit{\phi},\mathit{\lambda}),\end{array}$$

where (*φ*,*λ*) is the geodetic latitude and longitude, *W* is
the maximum order for decomposition, *A*_{W}(*φ*,*λ*) is
the regional anomaly caused by deep and large-scale geological bodies, and
*D*_{w}(*φ*,*λ*) is the local anomaly originating from
shallow and small-scale heterogeneous substances. Wavelet analysis generates
low-order wavelet details that are constant despite the decomposition order;
only high-order wavelet details and the corresponding wavelet approximation
change with decomposition order. Based on this, we can choose the proper
decomposition order to obtain desirable solutions.

According to the solution to the two-dimensional Laplace equation, each
*D*_{w}(*φ*,*λ*) in Eq. (4) can be expressed
as (Spector and Grant, 1970; Syberg, 1972; Cianciara and Marcak, 1976)

$$\begin{array}{}\text{(5)}& {D}_{\mathrm{w}}\left(\mathit{\phi},\mathit{\lambda}\right)=\sum _{\mathit{\phi}}\sum _{\mathit{\lambda}}{G}_{K}{e}^{i\mathrm{2}\mathit{\pi}({K}_{\mathit{\phi}}\mathit{\phi}+{K}_{\mathit{\lambda}}\mathit{\lambda})}{e}^{\mathrm{2}\mathit{\pi}KH},\end{array}$$

where *G*_{K} denotes the amplitude, $K=\sqrt{{K}_{\mathit{\phi}}^{\mathrm{2}}+{K}_{\mathit{\lambda}}^{\mathrm{2}}}$ is the wave number, and *H* is the elevation of *D*_{w}(*φ*,*λ*). Thus, *G*_{K} can be determined as

$$\begin{array}{}\text{(6)}& {G}_{K}=\sum _{\mathit{\phi}}\sum _{\mathit{\lambda}}{D}_{\mathrm{w}}\left(\mathit{\phi},\mathit{\lambda}\right){e}^{-i\mathrm{2}\mathit{\pi}({K}_{\mathit{\phi}}\mathit{\phi}+{K}_{\mathit{\lambda}}\mathit{\lambda})}{e}^{\pm \mathrm{2}\mathit{\pi}KH}.\end{array}$$

When *H*=0, Eq. (6) can be written as

$$\begin{array}{}\text{(7)}& ({G}_{K}{)}_{\mathrm{0}}=\sum _{\mathit{\phi}}\sum _{\mathit{\lambda}}{D}_{\mathrm{w}}\left(\mathit{\phi},\mathit{\lambda}\right){e}^{-i\mathrm{2}\mathit{\pi}({K}_{\mathit{\phi}}\mathit{\phi}+{K}_{\mathit{\lambda}}\mathit{\lambda})}.\end{array}$$

Inserting Eq. (7) into Eq. (6), *G*_{K} is rewritten as

$$\begin{array}{}\text{(8)}& {G}_{K}=({G}_{K}{)}_{\mathrm{0}}{e}^{\pm \mathrm{2}\mathit{\pi}KH}.\end{array}$$

Hence,

$$\begin{array}{}\text{(9)}& {P}_{K}=({P}_{K}{)}_{\mathrm{0}}{e}^{\pm \mathrm{4}\mathit{\pi}KH},\end{array}$$

where *P*_{K}=(*G*_{K})^{2} is the power. Then,

$$\begin{array}{}\text{(10)}& \mathrm{ln}{P}_{K}=\mathrm{ln}({P}_{K}{)}_{\mathrm{0}}\pm \mathrm{4}\mathit{\pi}KH,\end{array}$$

where ln*P*_{K} is the natural logarithm of *P*_{K}. Based on the linear
correlation between *K* and ln*P*_{K} in Eq. (10), the corresponding
average source depth *h*_{w} of *D*_{w}(*φ*,*λ*) can be estimated as (Spector and Grant, 1970; Xu et al.,
2018)

$$\begin{array}{}\text{(11)}& {h}_{\mathrm{w}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}{\displaystyle \frac{\mathrm{\Delta}\mathrm{ln}{P}_{K}^{\mathrm{w}}}{\mathrm{\Delta}{K}_{\mathrm{w}}}}\phantom{\rule{1em}{0ex}}w=\mathrm{1},\mathrm{2},\mathrm{\dots},W,\end{array}$$

where $\mathrm{\Delta}\mathrm{ln}{P}_{K}^{\mathrm{w}}$ and Δ*K*_{w} are the
change rates of $\mathrm{ln}{P}_{K}^{\mathrm{w}}$ and *K*_{w}, respectively. In
this manner, the corresponding average source depths ${h}_{\mathrm{w}}\left(w=\mathrm{1},\mathrm{2},\mathrm{\dots},W\right)$ of all decomposed wavelet details ${D}_{\mathrm{w}}\left(\mathit{\phi},\mathit{\lambda}\right)\left(w=\mathrm{1},\mathrm{2},\mathrm{\dots},W\right)$ and
wavelet approximation *A*_{W}(*φ*,*λ*) can
be estimated.

In this study, terrestrial and shipborne gravity anomalies are merged for
modeling. Gravity anomalies, Δ*g*, and quasi-geoid heights, *ζ*,
are related to the disturbing potential based on the multilayer approach as
follows:

$$\begin{array}{}\text{(12)}& \begin{array}{l}\mathrm{\Delta}g\left(z\right)\approx -{\displaystyle \frac{\mathrm{2}}{\left|z\right|}}{T}_{\mathrm{res}}\left(z\right)-{\displaystyle \frac{\partial {T}_{\mathrm{res}}\left(z\right)}{\partial \left|z\right|}}\\ =\sum _{i=\mathrm{1}}^{I}\sum _{n=\mathrm{1}}^{{N}_{i}}{\mathit{\beta}}_{i,n}\left(-{\displaystyle \frac{\partial}{\partial \left|z\right|}}{\mathrm{\Psi}}_{i}\left(z,{y}_{i,n}\right)-{\displaystyle \frac{\mathrm{2}}{\left|z\right|}}{\mathrm{\Psi}}_{i}\left(z,{y}_{i,n}\right)\right)\\ \mathit{\zeta}\left(z\right)={\displaystyle \frac{{T}_{\mathrm{res}}\left(z\right)}{\mathit{\gamma}\left(z\right)}}=\sum _{i=\mathrm{1}}^{I}\sum _{n=\mathrm{1}}^{{N}_{i}}{\mathit{\beta}}_{i,n}{\displaystyle \frac{{\mathrm{\Psi}}_{i}\left(z,{y}_{i,n}\right)}{\mathit{\gamma}\left(z\right)}},\end{array}\end{array}$$

where *γ* is the normal gravity value.

We assume the observational errors to be white noise with zero mean, and the gravity field model using the multilayer approach is expressed as the standard Gauss–Markov model:

$$\begin{array}{ll}{\displaystyle}{\mathit{l}}_{j}-{\mathit{e}}_{j}& {\displaystyle}=\phantom{\rule{0.25em}{0ex}}{\mathbf{A}}_{j}\mathit{x},\phantom{\rule{0.25em}{0ex}}E\mathit{\left\{}{\mathit{e}}_{j}\mathit{\right\}}=\mathrm{0},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}D\mathit{\left\{}{\mathit{e}}_{j}\mathit{\right\}}={\mathbf{C}}_{j}={\mathit{\sigma}}_{j}^{\mathrm{2}}{\mathbf{Q}}_{j}\\ \text{(13)}& {\displaystyle}& {\displaystyle}={\mathit{\sigma}}_{j}^{\mathrm{2}}{\mathbf{P}}_{j}^{-\mathrm{1}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}j=\mathrm{1},\mathrm{2},\mathrm{\dots},J,\end{array}$$

where *l*_{j} is the *m*_{j}×1 corresponding observation vector of group *j*; *e*_{j} is
the *m*_{j}×1 vector of observational errors; **A**_{j} is the
*m*_{j}×*N* design matrix of group *j*; ** x** is the

The data in different groups are assumed to be independent, and the weight
matrix **P**_{j} is assumed as the scaled diagonal matrix with white
noise properties since it is usually difficult to acquire a realistic full
error variance–covariance matrix in real-life measurements. Point-wise data
can be directly combined for modeling through the functions described above.
However, the heterogeneous characteristics of the data, in terms of spatial
coverages and noise properties, may result in an ill-conditioned normal
matrix (Panet et al., 2011). We apply the first-order Tikhonov regularization
for tackling the problem of the ill-conditioned matrix (Kusche and Klees,
2002; Wu et al., 2017a). For a given *α* (regularization parameter) and
*κ* (regularization matrix), the least squares solution of Eq. (13) is
(Klees et al., 2008)

$$\begin{array}{}\text{(14)}& \widehat{\mathit{x}}={\left(\sum _{j=\mathrm{1}}^{J}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma}}_{j}^{\mathrm{2}}}}{\mathbf{A}}_{j}^{T}{\mathbf{P}}_{j}{\mathbf{A}}_{j}\right)+\mathit{\alpha}\mathit{\kappa}\right)}^{-\mathrm{1}}\left(\sum _{j=\mathrm{1}}^{J}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma}}_{j}^{\mathrm{2}}}}{\mathbf{A}}_{j}^{T}{\mathbf{P}}_{j}{\mathit{l}}_{j}\right)\right).\end{array}$$

Furthermore, we use Monte Carlo variance component estimation (MCVCE) to estimate the appropriate variance factors for different observation groups and the regularization parameter (Koch and Kusche, 2002; Kusche, 2003; Wu et al., 2017c).

3 Numerical results and discussion

Back to toptop
The network design of the multilayer model contains several key parameters such as the number of layers, the depth of each layer, and the number of Poisson wavelets in each layer. Since the different layers are linked to the anomaly sources located at different depths, wavelet decomposition is used to separate and extract the contributions of the different anomaly sources. Moreover, the signal analysis is applied to determine the number of multiple layers based on background knowledge of local gravity field signals, while the power spectrum analysis is used to estimate the average depth of each layer. Then, we use a trial-and-error approach to determine the optimal number of Poisson wavelets in each layer. A flowchart representing the design of the multilayer model is shown in Fig. 1, and the details will be discussed in Sect. 3.1 and 3.2.

To determine the depths of the different layers, the residual gravity data
are decomposed into signals at the different scales based on wavelet
analysis. Spline interpolation is used to compute the gridded data, and Coif3
basis functions are chosen for wavelet decomposition (Xu et al., 2017). The
preliminary maximum order for wavelet decomposition is arbitrarily chosen to
some extent. However, since low-order wavelet details are constant despite
change in the decomposition order, we can preliminarily choose a predefined
order and implement wavelet decomposition; we then analyze the derived
details in order to choose the optimal order. For signals useful for
constructing the multilayer model that remain unseparated, we change the
decomposition order until all the useful signals have been extracted.
Otherwise, we truncate to a specific order and compute the wavelet details and
approximation to conduct the multilayer network design. By trial and
error, the preliminary order for decomposition is chosen as nine. Figure 2
shows the derived wavelet details (the corresponding statistics are provided
in Table 1). With the increase in the decomposition order, more
long-wavelength features occur. Specifically, the low-order wavelet details
indicate high-frequency signals stemming from shallow and small-scale
substances, while high-order wavelet details with long-wavelength patterns
reflect anomalies caused by deep and large-scale geological bodies. The first-
and second-order details (i.e., *D*_{1} and *D*_{2}) seem to be dominated by
high-frequency signals that correlate strongly with the local topography (the
local digital terrain model (DTM) can be seen in Fig. 1 in Wu et al., 2017c).
We attribute this to the uncorrected topographical signals in the RTM
corrections; these remain due to the inaccuracy of the density parameters in
RTM and the limitations of the DTM in terms of spatial resolution and precision.
As a result, high-frequency signals originating from local topographical
variations cannot be thoroughly recovered from RTM reduction, and
consequently, the uncorrected signals leak into the first- and second-order
details. To avoid these high-frequency errors propagating into the final
solution, we neglect these two wavelet details in designing the network of
multilayer model. With nine layers, we observe that *D*_{9} reveals large-scale
signals with wavelengths of hundreds of kilometers. Given that the mean
distance between measured gravity data in this target area is approximately
6–7 km and the spatial resolution of the applied GGM (i.e., GOCO05S) is
roughly 72 km, the spectral contents of the residual signals to be recovered
is roughly between several kilometers and tens of kilometers within the RCR
framework, i.e., approximately between degrees 250 and 3000 in terms of
spherical harmonics representation. The spectral contents of the ninth-order
wavelet details exceed the frequency bands of the signals to be modeled;
thus, the maximum order for wavelet decomposition is truncated to eight. In
this manner, the third- to eighth-order (*D*_{3}–*D*_{8}) wavelet details
and the final approximation (*A*_{8}) (see the information in Fig. 3 and
Table 2) are applied for constructing the multilayer model; the model then
consists of seven layers at various depths.

The order of Poisson wavelets is fixed at three to achieve a good compromise
between localization in the space and frequency domains (Panet et al., 2011).
In addition, the depth and number of Poisson wavelets are crucial factors
affecting the quality of the regional solution (Klees et al., 2008). Poisson
wavelets belonging to different layers are placed on the Fibonacci grids at
various depths beneath the Earth's surface, and the power spectrum analysis
is applied to estimate the depths. In Fig. 4, the green curves show the
radially averaged logarithm power spectrums of signals at different scales,
and the red straight lines represent the slopes of the spectrums, indicating
the depths of the corresponding layers. The red lines represent the rates of
change for logarithmic power relative to the wave number, estimated by an
autoregressive method. The initial and terminal points of the red lines are
the inflection points of the curves, recognized according to the trend of the
curves (Xu et al., 2018). Table 4 provides the estimated depths of the
different layers, limited between 4 and 60 km. The shallowest layer is
located 4.5 km underneath the Earth's surface, while the deepest layer is
estimated to be approximately 59.2 km below the Earth's surface. The
thickness of the sediments in the study area is approximately
2–4 km, and the thickness of the upper–middle crust is roughly
15–20 km (Artemieva and Thybo, 2013). Thus, the first four
layers (layer 1, layer 2, layer 3, and layer 4) are located between the
sediments and upper–middle crust. The corresponding wavelet details (*D*_{3}, *D*_{4}, *D*_{5}, and *D*_{6}) comprise small-scale patterns due to the
highly heterogeneous structure of the crust. *D*_{3} and *D*_{4} correspond
to the tectonic structure in the upper crust. The distributions of *D*_{3}
and *D*_{4} (at the average depths of 4.5 and 9.2 km, respectively) on
land are more dispersed than those in the ocean, and the tectonic structure
underneath the land is found to be more complex than that beneath the ocean
in the upper crust. Moreover, the gravity anomalies in the northern part of
North Sea are more dispersed than those in the central and southern parts of
North Sea, which is consistent with the fact that the Viking Graben and two
basins (i.e., Forth Approaches Basin and Norwegian–Danish Basin) are located
in the northern and southern parts of North Sea, respectively (e.g., see
Fichler and Hospers, 1990, and Blundell et al., 1991). The mean source
depths of *D*_{5} and *D*_{6} are 13.7 and 19.6 km, respectively; they
correspond to the depth of the middle crust. The gravity anomalies in these
two layers present apparent positive–negative alternating patterns, which may
be interpreted as the crustal shearing and extrusion (Blundell et al., 1991;
Ziegler and Dèzes, 2006). The last three layers (layer 5, layer 6, and
layer 7) are assumed to be located between the Moho surface and upper mantle,
considering that the Moho depth in the region is approximately
25–30 km (Grad and Tiira, 2009), and the corresponding details
(*D*_{7}, *D*_{8}, and *A*_{8}) become smoother and more long-wavelength
signals occur. *D*_{7}, with a mean source depth of 27.0 km, primarily
reflects the Moho undulation. The distribution of positive–negative
alternating gravity anomalies in *D*_{7} is nearly south–north oriented,
which is in agreement with the features of the Moho relief in the area
(Fichler and Hospers, 1990; Ziegler and Dèzes, 2006). The average source
depths of *D*_{8} and *A*_{8} are 32.3 and 59.0 km, respectively,
corresponding to the depth of the upper mantle; this indicates that the
density distribution of the upper mantle is relatively smooth. Overall, these
decomposed gravity anomalies can reveal the tectonic structure of the study
area at different depths.

A trial-and-error approach is used to estimate the number of Poisson wavelets at each layer (Wittwer, 2009). For a specific layer with a fixed depth, we predefine different numbers of Poisson wavelets to form a certain number of Fibonacci grids. Then, the signals reconstructed from these grids are compared with the true values, i.e., those derived from wavelet decomposition. The parameter that obtains the smallest difference between the modeled and true signals is considered as the optimal parameter. By trial and error, the spatial resolutions of the Fibonacci grids (mean distance between Poisson wavelets) are changed from 20 to 14 km with a step of 1 km. Table 4 shows the accuracies of the solutions derived from the different Fibonacci grids of the multiple layers, and we take the situation of the first layer, for instance. With increase in the number of Poisson wavelets, the standard deviation of the differences between the reconstructed and true signals decreases gradually to 0.12 mGal when the spatial resolution of the grid increases to 16 km. Beyond this point, no significant changes occur on incorporating more Poisson wavelets. Moreover, introducing more Poisson wavelets increases the overlap between them, which may lead to highly ill-conditioned normal matrices, and the associated heavy regularization may decrease the solution quality (Wu et al., 2017b). The optimal mean distance between Poisson wavelets of the first layer is estimated as 16 km. Similarly, the spatial resolutions for the remaining layers can be determined in this way (see Table 4).

For regional gravity field recovery, point-wise terrestrial and shipborne gravity anomalies are combined. Since there is no accurate information on the accuracies of terrestrial and shipborne data, we assume an accuracy of 2 mGal for both of these types of data, and the posterior variance factors of different data are estimated from MCVCE. The weights of different data indicate their relative contributions and play a key role in data combination. The estimated variance factors for terrestrial and shipborne gravity data are approximately 1.45 and 1.30 mGal, respectively, when we model the local gravity field based on the multilayer approach. For terrestrial data, the estimated accuracy agrees with that derived by Klees et al. (2008), i.e., 1.48 mGal for parts of the Netherlands. However, it is difficult to judge whether this estimate is realistic in other regions because of a lack of accurate information on data precision. For shipborne data, the computed value of 1.30 mGal is smaller than the results of crossover adjustments, where the standard deviation for the residuals at the crossovers was estimated to be approximately 2.0 mGal (Slobbe, 2013). However, this value may be too optimistic, considering that much of the shipborne data were collected decades ago without GPS navigation. The first-order Tikhonov regularization is used to tackle the ill-conditioned problem, and the convergent regularization parameter is estimated to be approximately $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{5}}$ using the MCVCE method. Details on regularization parameter estimation and comparisons with different methods can be found in Wu et al. (2017b).

The performance of the single-layer approach is also investigated for
comparison. The parameterization of the local gravity field based on the
single-layer approach has been described in detail by, e.g., Klees et
al. (2008) and Slobbe (2013). By trial and error, the single layer of the
Poisson wavelets' grid is found to be located 40 km beneath the Earth's
surface, and the mean distance between the Poisson wavelets is defined as 8.7 km (Wu et al., 2016). Figure 5 shows the normalized spectrums for different
approaches. Considering the frequency range of the signals to be recovered in
the target area is approximately between degrees 250 and 3000 in the spherical
harmonics representation, we note the single-layer approach is only sensitive
to a part of the signal spectrum. It is sensitive approximately between
degrees 300 and 1200 if we assume the criterion for determining whether it is
sensitive or not within a specific frequency band to be half of the maximum
value of the normalized spectrum. However, for the high-frequency band
between degrees 1200 and 3000, this approach is less sensitive. On the
contrary, the multilayer approach effectively covers the spectrum of the
local gravity signals, and is sensitive to both the low- and high-frequency
bands. The residuals of data after least squares adjustment using different
methods are displayed in Fig. 6 (the boundary limits for this area are
contracted by 0.5^{∘} in all the directions to reduce edge
effects). The residuals derived from the multilayer approach are significantly
lower throughout the region compared with those obtained from the
single-layer approach, especially in the western parts of the UK, south of
Norway, and southwest parts of Germany. In these regions, high-frequency
signals that are correlated with local topography dominate the regional
gravity field. We also note improvements in the oceanic parts, especially in
the waters around the English Channel, Irish Sea, northwest of the North Sea,
and the Atlantic Ocean close to the northwest UK. Table 5 displays the standard
deviation (SD) value for the residuals of terrestrial (shipborne) gravity
anomalies, which decreases by 0.37 mGal (0.34 mGal) when the multilayer
approach is used. These results are reasonable, since the multilayer model
contains several layers shallower than 40 km, and the spectrums of these
layers shift to the high-frequency bands. Thus, the spectrum of the
multilayer approach is more sensitive to high-frequency signals, and
consequently, the local high-frequency signals can be better fitted by the
multilayer approach. It is also worth mentioning that the analysis of data
residuals cannot be treated as the only criterion to justify the performance
of different approaches. There are two major reasons for this. First, these
gravity data have been used for modeling purposes, and the SD values of data
residuals should be regarded as the internal agreement. Additionally, due to
the limitation of the accuracies of gravity data, we cannot arrive at firm
conclusions based only on the analysis of data residuals. It is also possible
that lower data residuals can be derived if we place the Poisson wavelets'
grid at a shallower depth when the single-layer approach is used. However, we
think a shallower single grid may reduce the data residuals but may not
derive a better solution when validated against the independent control data,
as described in detail by Wu et al. (2016). In the following part, we
introduce another high-quality independent data set for external validation,
i.e., GPS/leveling data, which gives us more confidence with respect to the
performance of different methods.

It is also of interest to implement an Akaike information criterion (AIC)
test for different models. Although, the multilayer model fits the gravity
observations better, it also increases the level of estimated parameters. AIC
rewards the goodness of fit of data but also includes a penalty as the
number of estimated parameters increases. In other words, it deals with the
trade-off between the goodness of fit of the data and the simplicity of the
model. The AIC value is an estimator of the relative quality of statistical
models for a given set of data, providing a means for model selection. The
model that yields the minimum AIC value may be more preferable (Akaike, 1974;
Burnham and Anderson, 2002). The definition for the AIC value can be seen in
Eq. (A1) in the Appendix. Since we model the gravity field in the
framework of least squares system, we can simply take AIC $=\mathrm{2}k+m\mathrm{ln}(\mathrm{RSS}/m)$ for
model comparison, where *k* is the number of estimated parameters in the model,
*m* is the number of observations, and RSS is the residual sum of squares.
For details, see the Appendix. In this study, 894 649 point-wise
gravity observations are used for modeling, and 47 504 and 19 477 parameters
are estimated in the multilayer and single-layer models, respectively. The
RSS values for the multilayer and single-layer model are computed as
8.8527×10^{5} and 1.3296×10^{6} mGal^{2}, respectively, based on the data residuals after the least
squares adjustment. Then, the AIC values for the multilayer and single-layer
models are estimated as 85 581 and 393 400, respectively. Based on these
statistics, we note that the multilayer model yields a smaller AIC value,
which may be more preferable because it achieves a better balance between the
goodness of fit of the data and the simplicity of the model.

To test the ability of realistic extrapolation of different regional models recovered from various methods, we introduce GPS/leveling data in the Netherlands (534 points), Belgium (2707 points), and parts of Germany (213 points) as the independent validation data. This is a comparison of the predicted values derived from the regional model (e.g., model computed from the multilayer or single-layer approach) and those derived from independent survey/measurements. These data were provided in terms of geometric quasi-geoid heights derived from the high-quality GPS measurements and leveling surveys. The overall estimated accuracy of these observed quasi-geoid heights was approximately at the 1 cm level. It is worth mentioning that these GPS/leveling data have not been combined for modeling, and their three-dimensional coordinates do not coincide with the positions of gravity data. For validating purposes, it is necessary to reconstruct the regional model based on the estimated Poisson wavelets' coefficients and coordinates of GPS/leveling points (see Eq. 12), and compute the gravimetric quasi-geoid heights at these predicted points. We compute the standard deviation of the point-wise difference between GPS/leveling data and the gravimetric quasi-geoid height derived from the regional approach. This serves as an external validation.

The validation results demonstrate that the discrepancies between the GPS/leveling points and quasi-geoid heights derived from the multilayer approach decrease substantially compared with those computed from the single-layer approach (Fig. 7). The most prominent improvements occur in the northwest of Belgium, west of Germany, and eastern parts of the Netherlands, which are in good agreement with the results for data residuals' analysis demonstrated in Fig. 6. As shown in Table 6, the accuracies of gravimetric quasi-geoids derived from the multilayer approach improve by 0.4, 0.9, and 1.1 cm in the Netherlands, Belgium, and parts of Germany, respectively. Moreover, the mean values indicate that the solution computed from the multilayer approach further reduces the biases between the gravimetric solution and local GPS/leveling data, with magnitudes of 0.8, 0.7, and 1.1 cm in these three regions, respectively, compared to the those modeled from the single-layer approach. From these results, we can see that the multilayer approach not only leads to a reduction for the data residuals but also generates a better solution assessed by the independent control data. To construct the multilayer model, we consider that the gravity signals are the sum of the contributions generated from the anomaly sources, and different layers are designed for recovering these contributions with heterogeneous spectral contents. As a result, the spectrum of the multilayer approach is sensitive to the frequency bands of local gravity signals, both in the low- and high-frequency bands, and the local signals may be better recovered. We also notice that there are still biases between the regional gravimetric solutions and local GPS/leveling data (see the mean values in Table 6), which are mainly due to the commission errors in the GGM and uncorrected systematic errors in the local gravity data and leveling systems (Fotopoulos, 2005). Generally, corrector surface (Fotopoulos, 2005; Nahavandchi and Soltanpour, 2006) or more complicated algorithms, like least squares collocation (Tscherning, 1978), boundary-value methodology (Klees and Prutkin, 2008; Prutkin and Klees, 2008), and a direct approach (Wu et al., 2017a), can be applied to reduce the systematic errors and properly combine GPS/leveling data and gravimetric solutions. However, since the objective of this study is to develop a multilayer approach for gravimetric quasi-geoid modeling that may be served as a basis for further geophysical applications, the derived quasi-geoid is not purely gravimetric with implementing the data merging approach. Furthermore, we only have well-distributed GPS/leveling data in a limited region, i.e., in the Netherlands, Belgium, and parts of Germany, while in other regions, no high-quality control data are available. If we use the locally distributed GPS/leveling data to remove these systematic errors and compute the combined quasi-geoid, the final solution may be distorted in other regions, especially around the ocean, because no control data exist in these regions. Thus, we do not implement the methods mentioned above for computing the combined quasi-geoid. We use the gravimetric model derived from the multilayer approach for the following study, which is hereafter denoted as QGNSea V1.0 (quasi-geoid over the North Sea version 1.0).

QGNSea V1.0 is compared with a regional model called EGG08 (Denker, 2013),
and four other recently published high-order GGMs, i.e., EGM2008 with a full
d/o of 2190 and 2159 (Pavlis et al., 2012), EIGEN-6C4 (d/o
2190) (Förste et al., 2014), GECO (d/o 2190) (Gilardoni et al., 2016),
and SGG-UGM-1 (d/o 2159) (Liang et al., 2018). The reason for choosing these
four GGMs for comparisons is that these models have relatively higher spatial
resolutions and better accuracies compared to most other available GGMs (see
the information in http://icgem.gfz-potsdam.de/home
last access: 19 November 2018). EGG08 is a regional
gravimetric quasi-geoid model in Europe, which was recovered by Stokes'
integral based on locally distributed gravity data. This model is provided in
terms of gridded data instead of spherical harmonics, and its spatial
resolution is 1 arcmin in latitude and 1.5 arcmin in longitude,
respectively (Denker, 2013). The other four models are global geopotential
models provided in terms of spherical harmonics, and EGM2008 was computed by
merging GRACE measurements, terrestrial, altimetry-derived, and airborne
gravity data. Since no GOCE data have been incorporated for developing
EGM2008, and the recently published GGMs have been developed by combining
GOCE data, which are supposed to improve the gravity field in the frequency
bands approximately from degrees 30 to 220 in the spherical harmonics
representation (Gruber et al., 2014). EIGEN-6C4 was computed by combining
GRACE, GOCE, and terrestrial gravity data and other data sets; GECO was
computed by incorporating the GOCE-only TIM R5 (d/o 250) solution into
EGM2008; and SGG-UGM-1 was computed by the combination of EGM2008 gravity
anomalies and GOCE gravity gradients and satellite-to-satellite tracking
data. The differences between QGNSea V1.0 and other models are shown in
Fig. 8 (the boundary limits for the area are reduced by 0.5^{∘} in
all directions to reduce edge effects), the magnitude of which reaches the
decimeter level. For EGG08, we note the most prominent differences appear in
the eastern parts of the Irish Sea and center of Germany. Different data
pre-processing procedures and methods for parameterization partly account for
these differences. For example, QGNSea V1.0 is recovered from the multilayer
approach using Poisson wavelets, and proper weights for different observation
groups are estimated through MCVCE, while the spectral combination technique
and spectral weights were implemented in EGG08 for merging heterogeneous data
(Denker, 2013). Larger differences are observed between QGNSea V1.0 and these
four GGMs, and remarkable differences are seen in southern Norway, northern
parts of the North Sea, eastern parts of the Irish Sea, and northwest parts
of Germany. These differences are interpreted as resulting from the different
modeling techniques, and the additional signals introduced by QGNSea V1.0,
stemming from the incorporation of more high-quality gravity data. The
evaluation results with GPS/leveling data displayed in Fig. 9 and Table 7
show that the gravimetric quasi-geoid inverted from the multilayer approach
has the best quality, especially in the north of the Netherlands and western
and eastern parts of Belgium. Note that we removed the mean values between
the gravimetric model (both for the regional models and GGMs) and local
GPS/leveling data, since these GGMs deviate from the local GPS/leveling data
by tens of centimeters or even more in this area, due to the commission
errors and uncorrected systematic errors in gravity data and inconsistencies
among different height systems. Thus, if the mean biases are not removed,
these differences can become dominated by the systematic errors, which is
undesirable for model comparison. The SD value of the misfit between the
GPS/leveling data and QGNSea V1.0 is 1.5 cm, while this value increases to
2.2 cm for EGG08. In contrast, the accuracies of the four GGMs,
approximately at 2.6 cm levels, are slightly worse than that of EGG08.
Compared to the GGMs, the added values introduced by local high-quality data
led to the primary improvements in QGNSea V1.0. We find that the four GGMs
have comparable accuracies. However, those developed by combining GOCE data
and EGM2008 (i.e., GECO and SGG-UGM-1) do not demonstrate better performance
than EGM2008 alone, with SGG-UGM-1 even showing a slightly worse performance
than EGM2008. This is particularly prominent in the eastern parts of Belgium.
However, the possible reasons require further investigation. A new European
gravimetric quasi-geoid model, EGG2015, is also observed to have been
computed, where the GOCE-derived GGMs were used as reference models (Denker,
2015). However, this model is not publicly available, and its performance
cannot be assessed in this local region. Systematic errors can be seen in the
results presented in Fig. 9. These errors remain because they cannot be
thoroughly removed by simply removing the mean differences. However, as
mentioned previously, the target of this study is to develop a multilayer
approach for gravimetric quasi-geoid modeling. Implementing the data merging
approach for combining local GPS/leveling and gravimetric model may lead to a
distorted solution. Thus, a detailed discussion regarding the removal of
these systematic errors is out of the scope of this study.

For further comparison, we compute the local mean dynamic topography (MDT), which illustrates the departure of the mean sea surface (MSS) from the quasi-geoid/geoid (Becker et al., 2014; Bingham et al., 2014). We compute the MDTs in a geodetic manner, with raw MDTs computed as the differences between MSS and local geoid/quasi-geoid models. The derived MDTs are further smoothed with a Gaussian filter to suppress the small-scale signals from the MSS or local geoid/quasi-geoid that cannot be resolved (Andersen et al., 2013). The DTU13MSS from 1993 to 2012 is chosen as the MSS, and this model is provided as the gridded data with a spatial resolution of 1 arcmin×1 arcmin (Andersen et al., 2013). Considering that QGNSea V1.0 and EGG08 have better performance than other models when validated against local GPS/leveling data, we only compute local MDTs based on these two gravimetric quasi-geoid models. DTU13MSS and QGNSea V1.0/EGG08 are directly combined to obtain the raw MDT. Then, a Gaussian filter with a correlation length of 6 km is further applied to smooth the derived MDT, considering the signals at very short scales cannot be recovered from the local gravity data due to the limited spatial resolution of the gravimetric measurements.

The MDTs modeled based on QGNSea V1.0 and EGG08 are denoted as
MDTNS_QGNSea and MDTNS_EGG08, respectively
(Fig. 10). The results of these models agree with each other in most
regions over the North Sea. Prominent signals such as the Norwegian coastal
currents can be seen in the MDTs; e.g., see Idžanović et al. (2017).
The signals observed in MDTNS_QGNSea do not provide a full
picture of Norwegian coastal currents due to the limited data coverage in
Norway and its neighboring ocean areas. In most areas of the North Sea, the
MDTs show considerably smooth patterns, indicating a small change in the sea
surface topography; this result is consistent with Hipkin et al. (2004).
However, extreme values are observed surrounding most offshore areas; e.g.,
see the features over the offshore regions close
to the Wash (around
53^{∘} N, 0.5^{∘} W) and Thames
(around 51.5^{∘} N, 1^{∘} W) estuaries in England, and along the coastal areas of
France, the Netherlands, and Germany. MDT signals in these areas are traditionally
difficult to model and are frequently identified as errors (Hipkin et al.,
2004). The problems for computing geodetic MDTs in offshore regions are
twofold. First, the quasi-geoid/geoid is poorly modeled in coastal areas due
to unfavorable data coverage, and data inconsistencies are usually observed
when combining land and marine gravity surveys. Further, the quality of
altimetry data is dramatically reduced near offshore areas, and associated
errors in the derived MSS propagate into the final MDT (Andersen et al.,
2013). However, airborne gravimetric survey provides a seamless method of
gravity measurements over land and oceans, which may improve this situation
(Andersen and Knudsen, 2000).

4 Conclusions

Back to toptop
A multilayer approach for gravity field recovery at the regional scale, within the framework of multi-resolution representation, is developed, where the residual gravity field is parameterized as the superposition of the multiple layers of Poisson wavelets located at different depths beneath the Earth's surface. Since the gravity signals are the sum of the contributions of the anomaly sources at different depths, we place the multiple layers of the model at the locations of the different anomaly sources. Further, wavelet decomposition and power spectrum analysis are applied to estimate the depths of the different layers.

To test the performance of this multilayer approach, we model a local gravimetric quasi-geoid model, QGNSea V1.0, over the North Sea and validate this model against independent control data. Based on wavelet decomposition and power spectrum analysis, multiple layers located between 4.5 and 59.2 km underneath the Earth's surface are constructed to capture signals at different scales. The numerical results show that the multilayer approach is sensitive to the spectrum of signals, both in the low- and high-frequency bands, while the traditional single-layer approach is only sensitive to parts of the signals' spectrum. Comparisons with the single-layer approach show that the multilayer approach fits the gravity observations better, especially in the regions where the gravity signals show strong correlations with the variation of local topography. Moreover, an AIC test, which estimates the relative quality of the statistical models for a given set of data, is introduced for model selection in view of the statistical test. The associated results demonstrate that the multilayer model obtains a smaller AIC value and achieves a better balance between the goodness of fit of data and the simplicity of the model. Evaluation using independent GPS/leveling data tests the ability of regional models recovered from different methods towards realistic extrapolation, and shows that QGNSea V1.0 using the multilayer approach fits the local GPS/leveling data better than that using the single-layer approach, by the magnitudes of 0.4, 0.9, and 1.1 cm in the Netherlands, Belgium, and parts of Germany, respectively. Further comparisons with the existing models show that QGNSea V1.0 is superior in terms of performance and may be beneficial for investigating ocean circulation in the North Sea and surrounding oceanic areas.

Future work should focus on further improving the QGNSea V1.0. First, a data-adaptive algorithm may be developed for designing the optimal network in the multilayer approach, such as an algorithm for choosing the order for wavelet decomposition and determining the number of multiple layers, since human interventions are currently needed for estimating these key parameters. Moreover, satellite data (e.g., K-band range rate data and gravity gradients) from the GRACE and GOCE missions can be combined with ground-based gravimetry and altimetry data through the multilayer approach. Doing so can further improve the quality of local gravity field recovery, especially in the long-wavelength bands. However, deeper layers than those used in this study to combine surface data may be implemented to incorporate satellite observations, since these data mainly contribute to low-frequency bands of the gravity field. In addition, the stochastic model may need to be refined. For instance, the effects of the GGM errors on the solutions can be quantified if the full error variance–covariance matrix of the spherical coefficients is incorporated into the stochastic model. Thus, the different data may be more properly weighted and the solutions may be further improved.

Code and data availability

Back to toptop
Code and data availability.

Model code is available at https://github.com/yihaowu/QGNSea (last access: 19 November 2018). Gravity data were provided by the British Geological Service; the Geological Survey of Northern Ireland; the Nordic Geodetic Commission; Bundesamt für Kartographie und Geodäsie (Germany); Institut für Erdmessung (Germany); Bureau Gravimétrique International IAG service (France); Banque de données Gravimétriques de la France; and Bureau de Recherches Géologiques et Minières (France). GPS/leveling data were provided by geoinformation and ICT of Rijkswaterstaat (RWS-AGI) and the GPS Kernnet of the Kadaster, National Geographic Institute (NGI) and the Royal Observatory (ROB), and Bundesamt für Kartographie und Geodäsie. The global geopotential models can be publicly accessed from http://icgem.gfz-potsdam.de/home (last access: 19 November 2018).

Appendix A: Akaike information criterion

Back to toptop
Suppose that we have a statistical model of some data, and the AIC value of the model is (Burnham and Anderson, 2002)

$$\begin{array}{}\text{(A1)}& \mathrm{AIC}=\mathrm{2}k-\mathrm{2}\mathrm{ln}\left(\widehat{L}\right),\end{array}$$

where *k* is the number of estimated parameters in the model, and $\widehat{L}$ is
the maximum value of the likelihood function for the model (Akaike, 1974;
Burnham and Anderson, 2002).

For least squares fitting, the maximum likelihood estimate for the variance of a model's residual distributions is

$$\begin{array}{}\text{(A2)}& {\widehat{\mathit{\theta}}}^{\mathrm{2}}=\mathrm{RSS}/m,\end{array}$$

where RSS is the residual sum of squares, and *m* is the number of
observations.

Then, the maximum value of a log-likelihood function of the least squares model is (Burnham and Anderson, 2002)

$$\begin{array}{ll}{\displaystyle}& {\displaystyle}-{\displaystyle \frac{m}{\mathrm{2}}}\mathrm{ln}\left(\mathrm{2}\mathit{\pi}\right)-{\displaystyle \frac{m}{\mathrm{2}}}\mathrm{ln}\left({\widehat{\mathit{\theta}}}^{\mathrm{2}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\widehat{\mathit{\theta}}}^{\mathrm{2}}}}\\ \text{(A3)}& {\displaystyle}& {\displaystyle}\mathrm{RSS}=-{\displaystyle \frac{m}{\mathrm{2}}}\mathrm{ln}(\mathrm{RSS}/m)+C,\end{array}$$

where *C* is a constant independent of the model.

Combining Eqs. (A1) and (A3), for the least squares model, the AIC value is expressed as

$$\begin{array}{}\text{(A4)}& \mathrm{AIC}=\mathrm{2}k+m\mathrm{ln}(\mathrm{RSS}/m)+C.\end{array}$$

Since only differences in AIC are meaningful, the constant *C* can be ignored,
and we can conveniently take $\mathrm{AIC}=\mathrm{2}k+m\mathrm{ln}(\mathrm{RSS}/m)$ for model comparisons.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-11-4797-2018-supplement.

Author contributions

Back to toptop
Author contributions.

YW and CX developed the algorithm. YW and ZL initiated the study and designed the numerical experiments. YW wrote the manuscript with the contributions from CX and BZ. All authors were involved in discussions throughout the development.

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements

Back to toptop
Acknowledgements.

The authors would like to give sincerest thanks to the three anonymous reviewers
and Cornelis Slobbe for the beneficial suggestions and comments, which
were of great value for improving the manuscript. We also thank the executive
editor, Lutz Gross, for the kind assistance and constructive comments. We
are grateful for the kind support from the editorial office. We acknowledge funding
from the Netherlands Vertical Reference Frame project. We thank Roland
Klees and Cornelis Slobbe from Delft University of Technology for kindly
providing the original software. This study was supported by the Fundamental
Research Funds for the Central Universities (no. 2018B07314), the National
Natural Science Foundation of China (nos. 41830110, 41474061, 41504015,
41474109), the State Scholarship Fund from Chinese Scholarship Council
(no. 201306270014), the Open Research Fund Program of the State Key
Laboratory of Geodesy and Earth's Dynamics (nos. SKLGED2018-1-2-E and
SKLGED2018-1-3-E), and Key Laboratory of Geospace Environment and Geodesy,
Ministry of Education, Wuhan University (no. 17-01-09).

Edited by: Lutz Gross

Reviewed by: three anonymous referees

References

Back to toptop
Akaike, H.: A new look at the statistical model identification, IEEE T. Automat. Contr., 19, 716–723, https://doi.org/10.1109/TAC.1974.1100705, 1974.

Andersen, O. B. and Knudsen, P.: The role of satellite altimetry in gravity field modeling in coastal areas, Phys. Chem. Earth., 25, 17–24, https://doi.org/10.1016/S1464-1895(00)00004-1, 2000.

Andersen, O. B., Knudsen, P., and Stenseng, L.: The DTU13 global mean sea surface from 20 years of satellite altimetry, in: Ocean Surface Topography Science Team Meeting, Boulder, Colo., USA, 8–11 October 2013.

Artemieva, I. M. and Thybo, H.: EUNAseis: A seismic model for Moho and crustal structure in Europe, Greenland, and the North Atlantic region, Tectonophysics, 609, 97–153, https://doi.org/10.1016/j.tecto.2013.08.004, 2013.

Audet, P.: Toward mapping the effective elastic thickness of planetary lithospheres from a spherical wavelet analysis of gravity and topography, Phys. Earth. Planet. In., 226, 48–82, https://doi.org/10.1016/j.pepi.2013.09.011, 2014.

Becker, S., Brockmann, J. M., and Schuh, W. D.: Mean dynamic topography estimates purely based on GOCE gravity field models and altimetry, Geophys. Res. Lett., 41, 2063–2069, https://doi.org/10.1002/2014GL059510, 2014.

Bentel, K., Schmidt, M., and Rolstad, D. C.: Artifacts in regional gravity representations with spherical radial basis functions, Journal of Geodetic Science, 3, 173–187, https://doi.org/10.2478/jogs-2013-0029, 2013.

Bingham, R. J., Haines, K., and Lea, D. J.: How well can we measure the ocean's mean dynamic topography from space?, J. Geophys. Res.-Oceans, 119, 3336–3356, https://doi.org/10.1002/2013JC009354, 2014.

Blundell, D. J., Hobbs, R. W., Klemperer, S. L., Scott-Robinson, R., Long, R. E., West, T. E., and Duin, E.: Crustal structure of the central and southern North Sea from BIRPS deep seismic reflection profiling, J. Geol. Soc. London, 148, 445–457, 1991.

Burnham, K. P. and Anderson, D. R.: Model Selection and Multimodel Inference: A practical information-theoretic approach, 2nd Edn., Springer-Verlag, New York, 2002.

Chambodut, A., Panet, I., Mandea, M., Diament, M., Holschneider, M., and Jamet, O.: Wavelet frames: an alternative to spherical harmonic representation of potential fields, Geophys. J. Int., 163, 875–899, https://doi.org/10.1111/j.1365-246X.2005.02754.x, 2005.

Cianciara, B. and Marcak, H.: Interpretation of gravity anomalies by means of local power spectra, Geophys. Prospect., 24, 273–286, https://doi.org/10.1111/j.1365-2478.1976.tb00925.x, 1976.

Denker, H.: Regional gravity field modeling: theory and practical results, in: Sciences of Geodesy – II, edited by: Xu, G., Springer, Berlin, Heidelberg, 185–291, 2013.

Denker, H.: A new European gravimetric (quasi)geoid EGG2015, International Union of Geodesy and Geophysics General Assembly, Prague, Czech Republic, 22 June–2 July 2015.

Eicker, A., Schall, J., and Kusche, J.: Regional gravity modeling from spaceborne data: case studies with GOCE, Geophys. J. Int., 196, 1431–1440, https://doi.org/10.1093/gji/ggt485, 2013.

Fengler, M. J., Freeden, W., and Michel, V.: The Kaiserslautern multiscale geopotential model SWITCH-03 from orbit perturbations of the satellite CHAMP and its comparison to the models EGM96, UCPH2002_02_0.5, EIGEN-1s and EIGEN-2, Geophys. J. Int., 157, 499–514, https://doi.org/10.1111/j.1365-246X.2004.02209.x, 2004.

Fengler, M. J., Freeden, W., Kohlhaas, A., Michel, V., and Peters, T.: Wavelet Modeling of Regional and Temporal Variations of the Earth's Gravitational Potential Observed by GRACE, J. Geodesy, 81, 5–15, https://doi.org/10.1007/s00190-006-0040-1, 2007.

Fichler, C. and Hospers, J.: Deep crustal structure of the northern North Sea Viking Graben: results from deep reflection seismic and gravity data, Tectonophysics, 178, 241–254, https://doi.org/10.1016/0040-1951(90)90150-7, 1990.

Förste, C., Bruinsma, S. L., Abrikosov, O., Lemoine, J. M., Schaller, T., Götze, H. J., Ebbing, J., Marty, J. C., Flechtner, F., Balmino, G., and Biancale, R.: EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse, The 5th GOCE User Workshop, Paris, France, 2014.

Fotopoulos, G.: Calibration of geoid error models via a combined adjustment of ellipsoidal, orthometric and gravimetric geoid height data, J. Geodesy, 79, 111–123, https://doi.org/10.1007/s00190-005-0449-y, 2005.

Freeden, W. and Schreiner, M.: Local multiscale modelling of geoid undulations from deflections of the vertical, J. Geodesy, 79, 641–651, https://doi.org/10.1007/s00190-005-0017-5, 2006.

Freeden, W., Gervens, T., and Schreiner, M.: Constructive Approximation on the Sphere (With Applications to Geomathematics), Clarendon Press, Oxford, 1998.

Freeden, W., Fehlinger, T., Klug, M., Mathar, D., and Wolf, K.: Classical globally reflected gravity field determination in modern locally oriented multiscale framework, J. Geodesy, 83, 1171–1191, https://doi.org/10.1007/s00190-009-0335-0, 2009.

Gilardoni, M., Reguzzoni, M., and Sampietro, D.: GECO: a global gravity model by locally combining GOCE data and EGM2008, Stud. Geophys. Geod., 60, 228–247, https://doi.org/10.1007/s11200-015-1114-4, 2016.

Grad, M. and Tiira, T.: The Moho depth map of the European Plate, Geophys. J. Int., 176, 279–292, https://doi.org/10.1111/j.1365-246X.2008.03919.x, 2009.

Gruber, T., Rummel, R., Abrikosov, O., and van Hees, R.: GOCE Level 2 Product Data Handbook, GO-MA-HPF-GS-0110, Issue 4.2, available at: https://earth.esa.int/documents/10174/1650485/GOCE_Product_Data_Handbook_Level-2 (last access: 19 November 2018), 2014.

Guo, J., Gao, Y., Hwang, C., and Sun, J.: A multi-subwaveform parametric retracker of the radar satellite altimetric waveform and recovery of gravity anomalies over coastal oceans, Sci. China-Earth Sci., 53, 610–616, https://doi.org/10.1007/s11430-009-0171-3, 2010.

Hipkin, R. G., Haines, K., Beggan, C., Bingley, R., Hernandez, F., Holt, J., and Baker, T.: The geoid EDIN2000 and mean sea surface topography around the British Isles, Geophys. J. Int., 157, 565–577, https://doi.org/10.1111/j.1365-246X.2004.01989.x, 2004.

Holschneider, M. and Iglewska-Nowak, I.: Poisson wavelets on the sphere, J. Fourier. Anal. Appl., 13, 405–419, https://doi.org/10.1007/s00041-006-6909-9, 2007.

Idžanović, M., Ophaug, V., and Andersen, O. B.: The coastal mean dynamic topography in Norway observed by CryoSat-2 and GOCE, Geophys. Res. Lett., 44, 5609–5617, https://doi.org/10.1002/2017GL073777, 2017.

Jiang, W., Zhang, J., Tian, T., and Wang, X.: Crustal structure of Chuan-Dian region derived from gravity data and its tectonic implications, Phys. Earth Planet. In., 212, 76–87, https://doi.org/10.1016/j.pepi.2012.07.001, 2012.

Klees, R. and Prutkin, I.: The combination of GNSS-levelling data and gravimetric (quasi-) geoid heights in the presence of noise, J. Geodesy, 84, 731–749, https://doi.org/10.1007/s00190-010-0406-2, 2008.

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

Koch, K. R. and Kusche, J.: Regularization of geopotential determination from satellite data by variance components, J. Geodesy, 76, 259–268, https://doi.org/10.1007/s00190-002-0245-x, 2002.

Kusche, J.: A Monte-Carlo technique for weight estimation in satellite geodesy, J. Geodesy, 76, 641–652, https://doi.org/10.1007/s00190-002-0302-5, 2003.

Kusche, J. and Klees, R.: Regularization of gravity field estimation from satellite gravity gradients, J. Geodesy, 76, 359–368, https://doi.org/10.1007/s00190-002-0257-6, 2002.

Liang, W., Xu, X., Li, J., and Zhu, G.: The determination of an ultra-high gravity field model SGG-UGM-1 by combining EGM2008 gravity anomaly and GOCE observation data, Acta Geodaeticaet Cartographica Sinica, 47, 425–434, https://doi.org/10.11947/j.AGCS.2018.20170269, 2018.

Mayer-Gürr, T., Pail, R., Gruber, T., Fecher, T., Rexer, M., Schuh, W.-D., Kusche, J., Brockmann, J.-M., Rieser, D., Zehentner, N., Kvas, A., Klinger, B., Baur, O., Höck, E., Krauss, S., and Jäggi, A.: The combined satellite gravity field model GOCO05s, EGU General Assembly, Vienna, Austria, 12–17 April 2015.

Naeimi, M., Flury, J., and Brieden, P.: On the regularization of regional gravity field solutions in spherical radial base functions, Geophys. J. Int., 202, 1041–1053, https://doi.org/10.1093/gji/ggv210, 2015.

Nahavandchi, N. and Soltanpour, A.: Improved determination of heights using a conversion surface by combining gravimetric quasi-geoid/geoid and GNSS-levelling height differences, Stud. Geophys. Geod., 50, 165–180, https://doi.org/10.1007/s11200-006-0010-3, 2006.

Omang, O. C. D. and Forsberg, R.: How to handle topography in practical geoid determination: three examples, J. Geodesy, 74, 458–466, https://doi.org/10.1007/s001900000107, 2000.

Panet, I., Kuroishi, Y., and Holschneider, M.: 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, 2011.

Pavlis, N. K., Holmes, S. A., Kenyon, S. C., and Factor, J. F.: The development and evaluation of Earth Gravitational Model (EGM2008), J. Geophys. Res., 117, B04406, https://doi.org/10.1029/2011JB008916, 2012.

Prutkin, I. and Klees, R.: On the non-uniqueness of local quasi-geoids computed from terrestrial gravity anomalies, J. Geodesy, 82, 147–156, https://doi.org/10.1007/s00190-007-0161-1, 2008.

Rummel, R., Balmino, G., Johannessen, J., Visser, P., and Woodworth, P.: Dedicated gravity field missions-Principle and aims, J. Geodyn., 33, 3–20, https://doi.org/10.1016/S0264-3707(01)00050-3, 2002.

Schmidt, M., Fabert, O., and Shum, C. K.: On the estimation of a multi-resolution representation of the gravity field based on spherical harmonics and wavelets, J. Geodyn., 39, 512–526, https://doi.org/10.1016/j.jog.2005.04.007, 2005.

Schmidt, M., Han, S. C., Kusche, J., Sanchez, L., and Shum, C. K.: Regional high- resolution spatiotemporal gravity modeling from GRACE data using spherical wavelets, Geophys. Res. Lett., 33, L08403, https://doi.org/10.1029/2005GL025509, 2006.

Schmidt, M., Fengler, M., Mayer-Gürr, T., Eicker, A., Kusche, J., Sánchez, L., and Han, S. C.: Regional gravity modeling in terms of spherical base functions, J. Geodesy, 81, 17–38, https://doi.org/10.1007/s00190-006-0101-5, 2007.

Sjöberg, L. E.: A discussion on the approximations made in the practical implementation of the remove-compute-restore technique in regional geoid modelling, J. Geodesy, 78, 645–653, https://doi.org/10.1007/s00190-004-0430-1, 2005.

Slobbe, D. C.: Roadmap to a mutually consistent set of offshore vertical reference frames, PhD thesis, Delft University of Technology, the Netherland, 2013.

Slobbe, D. C., Klees, R., and Gunter, B. C.: Realization of a consistent set of vertical reference surfaces in coastal areas, J. Geodesy, 88, 601–615, https://doi.org/10.1007/s00190-014-0709-9, 2014.

Spector, A. and Grant, F. S.: Statistical models for interpreting aeromagnetic data, Geophysics, 35, 293–302, https://doi.org/10.1190/1.1440092, 1970.

Syberg, F. J. R.: A Fourier method for the regional-residual problem of potential fields, Geophys. Prospect., 20, 47–75, https://doi.org/10.1111/j.1365-2478.1972.tb00619.x, 1972.

Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C.: The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res. Lett., 31, L09607, https://doi.org/10.1029/2004GL019920, 2004.

Tenzer, R., Klees, R., and Wittwer, T.: Local Gravity Field Modelling in Rugged Terrain Using Spherical Radial Basis Functions: Case Study for the Canadian Rocky Mountains, in: Geodesy for Planet Earth, edited by: Kenyon, S., Pacino, M., and Marti, U., International Association of Geodesy Symposia, Springer, Berlin, Heidelberg, 401–409, 2012.

Tscherning, C. C.: Introduction to functional analysis with a view to its application in approximation theory, in: Approximation methods in geodesy, edited by: Moritz, H. and Sünkel, H., Wichmann H, Karlsruhe, Germany, 157–192, 1978.

Wang, J., Guo, J., Liu, X., Shen, Y., and Kong, Q.: Local oceanic vertical deflection determination with gravity data along a profile, Mar. Geod., 41, 24–43, https://doi.org/10.1080/01490419.2017.1380091, 2018.

Wang, Y., Saleh, J., Li, X., and Roman, D. R.: The US Gravimetric Geoid of 2009 (USGG2009): model development and evaluation, J. Geodesy, 86, 165–180, https://doi.org/10.1007/s00190-011-0506-7, 2012.

Wittwer, T.: Regional gravity field modelling with radial basis functions, PhD thesis, Delft University of Technology, the Netherlands, 2009.

Wu, Y. and Luo, Z.: The approach of regional geoid refinement based on combining multi-satellite altimetry observations and heterogeneous gravity data sets, Chinese J. Geophys.-Ch., 59, 1596–1607, https://doi.org/10.6038/cjg20160505, 2016.

Wu, Y., Luo, Z., and Zhou, B.: Regional gravity modelling based on heterogeneous data sets by using Poisson wavelets radial basis functions, Chinese J. Geophys.-Ch., 59, 852–864, https://doi.org/10.6038/cjg20160308, 2016.

Wu, Y., Luo, Z., Chen, W., and Chen, Y.: High-resolution regional gravity field recovery from Poisson wavelets using heterogeneous observational techniques, Earth Planets Space, 69, 34, https://doi.org/10.1186/s40623-017-0618-2, 2017a.

Wu, Y., Zhong, B., and Luo, Z.: Investigation of the Tikhonov regularization method in regional gravity field modeling by Poisson wavelets radial basis functions, J. Earth. Sci., 9, 1–10, https://doi.org/10.1007/s12583-017-0771-32016, 2017b.

Wu, Y., Zhou, H., Zhong, B., and Luo, Z.: Regional gravity field recovery using the GOCE gravity gradient tensor and heterogeneous gravimetry and altimetry data, J. Geophys. Res.-Sol. Ea., 122, 6928–6952, https://doi.org/10.1002/2017JB014196, 2017c.

Xu, C., Liu, Z., Luo, Z., Wu, Y., and Wang, H.: Moho topography of the Tibetan Plateau using multi-scale gravity analysis and its tectonic implications, J. Asian. Earth Sci., 138, 378–386, https://doi.org/10.1016/j.jseaes.2017.02.028, 2017.

Xu, C., Luo, Z., Sun, R., Zhou, H., and Wu, Y.: Multilayer densities using a wavelet-based gravity method and their tectonic implications beneath the Tibetan Plateau, Geophys. J. Int., 213, 2085–2095, https://doi.org/10.1093/gji/ggy110, 2018.

Ziegler, P. A. and Dèzes, P.: Crustal evolution of western and central Europe, Geological Society, London, Memoirs, 32, 43–56, https://doi.org/10.1144/GSL.MEM.2006.032.01.03, 2006.

Short summary

A multilayer approach is parameterized for model development, and the multiple layers are located at different depths beneath the Earth’s surface. This method may be beneficial for gravity/manget field modeling, which may outperform the traditional single-layer approach.

A multilayer approach is parameterized for model development, and the multiple layers are...

Geoscientific Model Development

An interactive open-access journal of the European Geosciences Union