Combination of different observation types through a multi-resolution representation of the regional gravity field using the pyramid algorithm and parameter estimation

The optimal combination of different types of gravity observations is the key to obtaining a high-resolution and high-precision regional gravity model. Current studies based on spherical radial basis functions (SRBFs) majorly consider a single-level approach for data combination. Despite the promising results reported in numerous publications, it has been suspected that the single-level model might be biased towards high-resolution measurements. Instead, a multi-resolution representation (MRR) can be applied to further take into consideration the varying spectral sensitivities of different observation techniques. In this study, we develop a new MRR scheme based on the pyramid algorithm and sequential parameter estimation. We propose strategies to solve the challenges in the practical application of the pyramid algorithm, and this study represents its first successful realization in regional gravity field modeling. The modeling results based on both simulated and real gravity data show that either the single-level approach or the MRR without pyramid algorithm is able to capture gravity information from lower resolution measurements as sufficient as our newly developed MRR algorithm. In the simulated case, the RMS error w.r.t. the validation data obtained by the MRR based on the pyramid algorithm decreases by 50% and 35%, in comparison to that of the single-level model and the MRR without pyramid algorithm, respectively. In the real case, the improvement achieved by the MRR based on the pyramid algorithm is 35% and 23% in the onshore area, and it reaches 63% and 57% in the offshore area, compared to the single-level approach and the MRR without pyramid algorithm, respectively.


Introduction
Regional gravity field modeling is an important topic in geodesy; it plays an essential role for applications in geophysics and physical height system realization. A highresolution and high-precision geoid model derived from regional gravity field modeling is the key to the realization of the International Height Reference System (IHRS, Sánchez et al. 2021). In the last decades, spherical radial basis functions (SRBFs, see e.g., Freeden et al. 1998;Freeden and Michel 2004;Schmidt et al. 2007) have been widely used for regional gravity field modeling. Thanks to their localizing B Qing Liu qingqing.liu@tum.de 1 Deutsches Geodätisches Forschungsinstitut der Technischen Universität München (DGFI-TUM), Arcisstrasse 21, 80333 Munich, Germany features in both the spectral and the spatial domain, SRBFs are an appropriate approach to consider the heterogeneity of different gravity data types with varying spectral and spatial resolutions.
In regional gravity field modeling, local high-resolution gravity measurements (e.g., terrestrial, airborne, and shipborne gravity data) are usually combined with mediumresolution data (e.g., inferred from satellite altimetry missions) and low-resolution global data provided by dedicated satellite gravity missions, such as the Gravity Recovery and Climate Experiment (GRACE, Tapley et al. 2004) and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE, Rummel et al. 2002). Each observation technique has different spectral sensitivities and varying data distribution. Thus, a data combination method that extracts the maximum gravity information from different measurements is required to ensure the best possible high-precision regional gravity model. Typically, the data combination based on SRBFs is implemented using a single-level approach (see, e.g., Bentel et al. 2013a, b;Lieb et al. 2016;Bucha et al. 2016;Wu et al. 2017a, b;Liu et al. 2020a, b). In this approach, different observation techniques are combined at the maximum degree of expansion, and their contributions to the final gravity model are determined by the relative weighting between each data set. This approach is straightforward, and promising modeling results have been reported in the aforementioned publications. However, the contribution of measurements with medium to low resolution (e.g., altimetry data, satellite gravity data) could be understated, since they are not sensitive at high spectral degrees. For instance, Klees et al. (2018) demonstrate that the single-scale model lacks the flexibility to deal with data sets of significantly different bandwidths, while Wu et al. (2018) state that the single-level approach may fail to extract the full information contained in the gravity data.
To take the spectral sensitivity of different observation techniques into consideration and to combine heterogeneous data by spectral weights, a spectral combination can be applied (see, e.g., Sjöberg 1981;Kern et al. 2003;Denker 2013). The spectral combination based on SRBFs can be realized through a multi-resolution representation (MRR). The MRR of the gravity field was initially proposed by Freeden et al. (1998), Freeden (1999), and Haagmans et al. (2002), and its realization has been investigated in the last two decades. The MRR is often associated with wavelet functions. Kusche et al. (1998) conclude that the spherical wavelet functions show outstanding properties in the medium frequency band, and they are thus recommended for modeling the medium wavelength part of the gravity field. Beylkin and Cramer (2002) present several multi-resolution gravity models and show their performance results. Fengler et al. (2007) set up an MRR approach based on GRACE data using the cubic polynomial (CuP) wavelets and show that they are an appropriate tool to investigate regional and temporal variations of the Earth's gravity field. Panet et al. (2011) introduce the iterative domain decomposition by extending the Poisson wavelet modeling approach (Holschneider et al. 2003) to combine the surface (land, marine and altimetry) gravity data with satellite gravity data. Bolkas et al. (2016) state that wavelet decomposition is a useful tool for fusing terrestrial gravity data with satellite and airborne data and show that the multiscale fused model is able to fill data gaps. Schmidt et al. (2005Schmidt et al. ( , 2006Schmidt et al. ( , 2007 develop an MRR approach where the representation is decomposed into an expansion in terms of spherical harmonics for the longwavelength part, and a number of frequency-dependent detail signals in terms of wavelets for the medium and high frequency parts. This MRR was applied for regional gravity field modeling by Lieb (2017) and Wu et al. (2018), where the coefficients for calculating the detail signals are estimated simultaneously at each resolution level, using all or selected data groups. However, if all observation types are used at each level, the corresponding detail signals are strongly correlated (Lieb 2017). If only specific data sets are used, large data gaps may occur, and the prior information is insufficient for filling these data gaps at higher resolution levels. This leads to large erroneous effects in the output signals. Consequently, Lieb (2017) recommends for future work the implementation of the pyramid algorithm (Freeden et al. 1998) to consider all available information by connecting the different levels. Klees et al. (2018) show that a two-scale (or multi-scale) model needs to be applied in combination with a sequential estimation of the scale-dependent coefficients. The sequential estimation can be realized by applying the pyramid algorithm, which determines the coefficients of the lower resolution levels sequentially from the coefficients of the higher level, by successive low-pass filtering.
Although the pyramid algorithm was proposed nearly two decades ago, its practical application based on parameter estimation has never been realized in regional gravity field modeling, until this study. This is due to the challenges of applying the pyramid algorithm in the regional case, as stated by Lieb (2017): (1) it is difficult to set up a proper low-pass filter matrix in the regional case, and (2) the margin size has to be adapted appropriately at each level to minimize edge effects. We come up with methods in this study to solve these difficulties: (1) a low-pass filter is newly introduced in case of using the Reuter grid; (2) the choice of the margin size and the setting up of the estimation model at each resolution level are characterized; (3) strategies for reducing edge effects in the calculated detail signals are proposed.
A new MRR scheme based on the pyramid algorithm is developed, as visualized in Fig. 1, the text in blue highlights the novelty of this work. One main innovation is the realization of the sequential parameter estimation in the MRR scheme. The coefficients are firstly estimated at the highest chosen resolution level I using only the high-resolution data set(s) 1, and they are used to compute the detail signal G I by wavelet functions. Then, these estimated coefficients are transformed to the next lower level I − 1 by applying a lowpass filtering, i.e., the pyramid algorithm using the proposed low-pass filter matrix L I −1 . At the lower level I −1, the coefficient vector is updated by the lower-resolution data set(s) 2 introduced at this level. This updated coefficient vector is then used in combination with wavelet functions to calculate the detail signal G I −1 . Continuing this process until the lowest level i of the MRR, all data sets are introduced into the scheme at different resolution levels and the coefficient vectors as well as the detail signals are obtained. At each lower resolution level, the coefficients from the pyramid algorithm are combined directly with the new data set(s), through the parameter estimation procedure. In the two-scale model presented by Klees et al. (2018), the coefficients from the high-resolution level are used to generate a new pseudo Fig. 1 The developed MRR scheme based on the pyramid algorithm. The text in blue highlights the novelty of this work data set (with the same resolution as the low-resolution data set), which is then combined with the low-resolution data set to estimate the coefficients of the lower level. However, the authors point out that the high-resolution data set is required to have a larger area coverage than the low-resolution data set (due to edge effects), which frequently cannot be fulfilled when dealing with real data. In this study, we solve this limitation by the direct combination of the coefficient vector and the new data set(s), in analogy to a Kalman filter.
In this research, we additionally demonstrate that: (1) different observation types can be introduced into the evaluation model at the spectral level of their highest sensitivities, which makes it possible to benefit from the individual strengths of each data set optimally; (2) since all computations within the pyramid algorithm are based on linear equation systems (Schmidt et al. 2015), all covariance information can easily be calculated following the law of error propagation from the higher levels and serve as input for the lower levels; (3) as the number of required SRBFs decreases from the highest level to the lower levels, the design matrices of the lowerresolution data sets are now calculated with a smaller number of grid points, which reduces the computational effort significantly. We test the MRR based on the pyramid algorithm by using simulated data and then apply it to real data in different study areas. The modeling results are directly compared to those obtained from the single-level approach and the MRR without using the pyramid algorithm (where the coefficients are estimated independently at each level using all types of observations) in order to highlight the benefits of the MRR based on the pyramid algorithm.
This work is organized as follows: In Sect. 2, we introduce the fundamental concepts of the MRR based on SRBFs, the wavelet functions, and the pyramid algorithm. Section 3 explains the newly developed MRR procedure based on the sequential parameter estimation, and shows how the coefficient vector from the pyramid algorithm is updated by including new observations at the lower levels. Sections 4 and 5 present the performance of our approach based on simulated and real data, respectively. The model configuration is explained, and the developed strategy for reducing edge effects is demonstrated. The results are evaluated and discussed, and comparisons are made with both the single-level approach and the MRR without pyramid algorithm. Finally, Sect. 6 provides conclusions and outlook.

Spherical radial basis functions
In general, SRBFs are isotropic localizing functions centered at grid points P k on a sphere R with radius R (Freeden et al. 1998;Schmidt et al. 2007), and can be defined by the Legendre series wherein x = r · r is the position vector of an arbitrary point P(ϕ, λ, r ), with latitude ϕ, longitude λ, and radial distance r = |x| = R + h ; h is the height of P over the sphere R ; r = [cos ϕ cos λ, cos ϕ sin λ, sin ϕ] T is the corresponding unit vector. x k = R · r k is the position vector of a grid point P k . P n is the Legendre polynomial of degree n, which is a function of the spherical distance between the point P and the grid point P k . n max is the maximum degree of the expansion. B n is the Legendre coefficients, which specify the shape of the SRBFs. In case of B n = 1 for all degree values n = 0, 1, . . . , n max , the SRBF is defined as the reproducing kernel K rep (x, x k ), which has the ability of the unique reproduction of a harmonic function. To be more specific, the convolution of a harmonic function F(x) (band-limited with n = 0, 1, . . . , n ≤ n max ) with the reproducing kernel K rep is equal to the original function (Schmidt et al. 2007), i.e., Any gravity observation y(x) can be represented as a series expansion in terms of the SRBFs, i.e., where K and d k are the number of grid points and the corresponding series coefficients, respectively. s(x) is the truncation error, and e(x) is the observation error. It is worth mentioning that the general expression (1) of the basis functions B(x, x k ) needs to be adapted for describing different gravitational functionals (e.g., the gravity disturbance δg or the gravity anomaly g). A list of basis functions adapted to different functionals can be found in Koop (1993) and Liu et al. (2020a).

Multi-resolution representation
The fundamental idea of the MRR is to split a given input signal into a smoothed version and a number of detail signals by successive low-pass filtering (Schmidt et al. 2007). Following Schmidt et al. (2006), the MRR of the gravity functional F(x) can be expressed as whereF(x) is a reference model, i.e., usually the longwavelength component from a global gravity model (GGM), and it is used as the background model within the removecompute-restore (RCR, e.g., Forsberg 1993) procedure.
is the detail signal of resolution level i, and s(x) is the truncation error. To set up the resolution level i, the frequency domain needs to be discretized. Each resolution level defines a frequency band, and thus, includes a certain number of degree values. In this study, the upper boundary (maximum degree value) n i of each level is defined as It is worth mentioning that the resolution levels can also be adapted differently, i.e., the base 2 in Eq. (5), which specifies the range of the frequency bands, can be chosen as a different number (see Schmidt et al. 2007 for more details). Figure 2 (cf. Lieb et al. 2016) visualizes the discretization of the frequency domain, and the relation between the resolution level i and its corresponding spectral degree n i as well as spatial resolution ρ i . The corresponding maximum spatial resolution of level i is ρ i = π R/n i , with R the Earth's radius. Gravity data obtained from different observation techniques can then be classified according to their spectral resolution. Terrestrial, airborne, and shipborne observations are high-resolution regional data, thus, they cannot be used to estimate the low-frequency part. Therefore, they should be combined with the medium-resolution satellite altimetry measurements, and the low-resolution satellite gravimetry data, such as GRACE and GOCE.

Wavelet functions
SRBFs (1) can act as high-pass, low-pass or band-pass filters, depending on the chosen Legendre coefficients B n , and a harmonic function F(x) can be filtered by it through a spherical convolution (Schmidt et al. 2007;Liu et al. 2020b). In case of using a band-limited SRBF, e.g., a spherical scaling function B(x, x k ) =: i (x, x k,i ), which means the Legendre coefficient B n =: φ n,i > 0 for degree n = 0, 1, . . . , n max =: 2 i −1 and φ n,i = 0 for degree n > 2 i − 1, the SRBF acts as a low-pass filter. Based on the scaling function, a spherical wavelet function i (x, x k,i ), which can be interpreted as a band-pass filter, is defined as where x k,i is the position vector of the grid point P k,i at resolution level i, and the Legendre coefficients ψ n,i are computed as The detail signal G i (x) of level i in Eq. (4) can then be defined in terms of the spherical wavelet function as where K i and d k,i are the number of grid points and the corresponding series coefficients, respectively, at the level i. As explained in Liu et al. (2020b), the coefficients d k,i do not depend on the choice of the SRBFs, as soon as the SRBFs are band-limited up to the same degree, i.e., their Legendre coefficients are equal to 0 for all degree values n > 2 i − 1. Thus, the same set of unknown coefficients d k,i at level i can be used with both the spherical scaling function and the spherical wavelet function of the same level. Also, it is possible to use different types of SRBFs in the analysis step (in which the unknown coefficients are estimated) and in the synthesis step (in which the estimated coefficients are used to calculate the output gravity models), respectively. In this study, we use the following spherical scaling functions (see Schmidt et al. 2007;Lieb et al. 2016): 1. The Shannon function, which could be interpreted as the reproducing kernel (see Sect. 2.1), its Legendre coefficients are given by 2. The Blackman function, its Legendre coefficients are given by where The Shannon scaling function is used in the analysis step to estimate the unknown coefficients, and the Blackman wavelet function (with Legendre coefficients ψ Bla n,i = φ Bla n,i − φ Bla n,i−1 ) is applied in the synthesis step to calculate the detail signals at each level. The reason for using the Shannon function in the analysis step is to avoid the loss of spectral information. An SRBF with smoothing features, such as the Blackman function, is used in the synthesis step to reduce erroneous effects (Lieb et al. 2016;Liu et al. 2020b). Figure 3 visualizes the characteristics of the Blackman wavelet function in the spatial (upper plot) and the spectral domain (bottom plot) for different resolution levels. As we can see, the Blackman wavelet function has smoothing features. At the same time, it has strict band-pass features, i.e., only the Legendre coefficients within the spectral band [2 i−2 , 2 i ) are not zero. The spectral weight, which is defined by the Legendre coefficients, increases from degree 2 i−2 to 2 i−1 and decreases from degree 2 i−1 to 2 i . In the spectral domain, the corresponding frequency band becomes wider at higher resolution levels. In the spatial domain, with increasing resolution level i, the peak becomes sharper. Freeden (1999) shows that the coefficients of neighboring resolution levels depend on each other linearly, and thus, can be computed successively. Therefore, the pyramid algorithm can be set up to determine the coefficients of the lower resolution levels from the coefficients of the higher level by a low-pass filtering. This procedure is based on a down sampling strategy, as the number of coefficients decreases at each level. With the pyramid algorithm, the coefficients at level i (i = i , i + 1, . . . , I − 2, I − 1) can be expressed as

The pyramid algorithm
is the covariance matrix of d i|i+1 obtained from the covariance matrix d i+1 following the law of error propagation. According to Schmidt et al. (2007), the low-pass filter matrix L i can be decomposed as where the integration weights associated with the grid points P k,i of level i, and it depends on the implemented type of grid. Although the Reuter grid (Reuter 1982) is one of the most commonly used grids in regional gravity field modeling, there are no previous studies about how to define W i for this type of grid. In this study, we develop a strategy to set up the corresponding integration weights, which read where Z i is the total number of the Reuter grid points on the sphere at level i. Reuter grids are regarded as equidistributed point systems on the sphere (Fengler et al. 2004;Eicker 2008), and the corresponding integration weights are constant for each grid point.

Parameter estimation
In this paper, a new MRR scheme is developed based on the pyramid algorithm using sequential parameter estimation. The initial step of the MRR procedure is to estimate the unknown coefficient vector at the highest resolution level I using parameter estimation (Koch 1999;Schmidt et al. 2007). Assuming that P groups of observations are used at level I , the estimated coefficient vector d I and the corresponding covariance matrix D( d I ) = d I can be obtained as (Liu et al. 2020a, b) y p,I is the observation vector of the p th gravity data set used at level I , P p,I is its positive definite weight matrix, and A p,I is the design matrix, which contains the corresponding (adapted) scaling functions. μ is the expectation vector of the coefficient vector d I , and P μ is its given positive definite weight matrix. σ 2 p,I and σ 2 μ are the corresponding variance factors for the observations y p,I and the prior information μ, respectively.
After estimating the scaling coefficient vector and its covariance matrix, they are used to calculate the scaling coefficient vector d I −1|I and the corresponding covariance matrix d I −1|I of the next lower level, by applying the pyramid algorithm according to Eqs. (12) and (13). Then, instead of transforming the coefficient vector d I −1|I directly to the next lower level I − 2, which is the usual procedure in previous studies about the pyramid algorithm, it is updated by the gravity observations introduced at this level I − 1. The updated coefficient vector d I −1 is then transformed to the level I − 2 following Eqs. (12) and (13). Assuming that Q groups of observations y q,I −1 are introduced at level I − 1, the combination of y q,I −1 and d I −1|I is realized through the parameter estimation: The updated coefficient vector of level I − 1 is estimated as with the covariance matrix The variance factors σ 2 1,I −1 , σ 2 2,I −1 , . . . , σ 2 Q,I −1 of the data sets y 1,I −1 , y 2,I −1 , . . . , y Q,I −1 are estimated by the variance component estimation (VCE, Koch and Kusche 2002).
The combination of can also be solved in analogy to the Kalman filter (Kalman 1960), where d I −1|I and d I −1|I can be regarded as the predicted state vector and the related predicted covariance matrix, respectively. Then the corrected state vector d I −1 as well as its covariance matrix d I −1 are computed by incorporating the involved measurements y I −1 at level I − 1 where K I −1 is the gain matrix ] T , and y I −1 being the covariance matrix of the observation vector y I −1 . After taking into account the matrix identities, the solution delivered by Eqs. (21) and (22) ends up identical to Eqs. (19) and (20). We refer to Koch (1999) and Erdogan et al. (2020) for the details of the matrix identities.
In the synthesis step, the estimated coefficient vector d I −1 and its covariance matrix d I −1 are used to calculate the estimated detail signal G I −1 as well as its covariance matrix where T is the vector of the estimated detail signal values at the computation points x 1 , x 2 , . . . , x C of the output gravity model. B I −1 is the design matrix, which contains the Blackman wavelet functions I −1 (x c , x k,I −1 ), as defined in Eq. (6), between the computation points of the gravity model and the grid points of level I − 1.
The developed MRR procedure is summarized in Fig. 4: after determining the estimated coefficient vector d I at the highest level I , it can be used to compute the detail signal G I and to start the pyramid algorithm, i.e., to calculate the coefficient vector d I −1|I of the next lower level. Then, this coefficient vector d I −1|I is updated by introducing the observations at level I −1. Continuing this process until the lowest level of the MRR, the scaling coefficients and the detail signals of each level can be obtained, and thus, the final gravity functional is obtained according to Eq. (4). In this process, the input gravity data obtained from different observation techniques can be introduced into the estimation model at the resolution level of their highest sensitivities. Hence, they are able to contribute to the final model with maximum gravity information. Typically, the terrestrial data are used at the highest level. Then the shipborne or airborne data can be introduced at a level lower, followed by the altimetry data and the satellite gravity data, if applicable (see Fig. 2). Another advantage of including data by levels is that the high frequencies of the airborne data, which have large noise (Childers et al. 1999), could be excluded to guarantee a stable solution (Jiang and Wang 2016). Moreover, since the data are now introduced at the lower levels which require less grid points, the size of the design matrices, and consequently, the computation time is significantly reduced.

Data
In the following, the MRR based on the pyramid algorithm is first realized and evaluated using simulated data to benefit from the availability of an accurate validation data set serving as the "truth". The study area is between 10 • and 20 • longitude and between 39 • and 45 • latitude (Fig. 5), covering parts of South Europe, the Adriatic Sea, and the Tyrrhenian Sea. Five types of gravity observations are used, namely terrestrial, airborne, and altimetry data, as well as satellite data from GOCE and GRACE. These data are simulated from the global gravity model GECO (Gilardoni et al. 2016), with the position of the observations provided by the IAG-ICCT (International Association of Geodesy -Inter Commission Committee on Theory) Joint Study Group (JSG) 0.3 ("Comparison of current methodologies in regional gravity field modeling"), running from 2011 to 2015. All observations Fig. 4 The multi-resolution representation (MRR) based on the pyramid algorithm are simulated in the sense of disturbing gravity field quantities, i.e., functionals of the disturbing potential. A detailed data description can be found in Liu et al. (2020a). The terrestrial observations (yellow dots in Fig. 5) are generated on a regular grid at the surface of the topography (DTM2006.0, Pavlis et al. 2006) with a grid spacing of 5 (which corresponds to a spatial resolution of 10 km). They are simulated up to spherical harmonic (SH) degree and order (d/o) 2190, in terms of the first order radial derivatives of the disturbing potential. The airborne data (red dots in Fig. 5) are located over the Adriatic Sea along synthetic flight tracks with an altitude of 2.5 km, generated in terms of the first order radial derivatives of the disturbing potential up to SH d/o 1600. The altimetry data (green dots in Fig. 5) are simulated up to SH d/o 1000 in terms of geoid height N , based on the real ground track of the altimetry mission Envisat (Extended Mission). The GOCE (grey dots in Fig. 5) and GRACE (blue dots in Fig. 5) data are simulated up to SH d/o 250, based on real satellite orbits with a time span of 61 days and one month, respectively. They are used in terms of the second-order radial derivatives of the disturbing potential for GOCE and the disturbing potential differences between the two satellites for GRACE. Observation noise is also generated and added to the gravity data, while the noise level is chosen according to the assumptions of the JSG 0.3. White noise with standard The study area (observation area ∂ O ) and simulated gravity data, including terrestrial (yellow dots), airborne (red dots), altimetry (green dots), GOCE (grey dots), and GRACE (blue dots) data. The black rectangle represents the investigation area ∂ I , where the final gravity model is calculated deviations of 0.01 mGal, 1 mGal, and 0.03 m is added to the terrestrial, airborne, and altimetry data, respectively. Colored noise (Austen and Grafarend 2004;Naeimi 2013) is added to the satellite data of GOCE and GRACE, with standard deviations of 10 mE and 8 · 10 −4 m 2 /s 2 , respectively. Validation data are disturbing potential values simulated from GECO with a spatial resolution of 5 × 5 and a maximum degree of 2190.

Model configuration
In this study case, the highest resolution level of MRR is chosen as I = 11, considering the spatial resolution of the data (see Fig. 2). The terrestrial data are used at the highest level to calculate the unknown coefficients of this level following Eq. (16), with the expectation vector μ set to the zero vector and the weight matrix P p,I and P μ set to be the identity matrix I (Lieb et al. 2016;Liu et al. 2020a, b). As listed in Table 1, the airborne data are included at level 10, the altimetry data are then added at level 9, followed by the GOCE observations at level 8 and GRACE observations at level 7. The long-wavelength component (up to level 6) is modeled by the GGM within the RCR procedure, i.e., the background model GECO up to spherical harmonic degree n 6 = 2 6 − 1 = 63 is removed from each observation, and then restored to the estimated model in the synthesis step. The single-level model and the MRR without the pyramid algorithm (see Sect. 1) are also calculated to serve as comparison scenarios. For the single-level model, all observations are combined at level 11, with the long-wavelength component (up to degree 63) modeled by GECO. For the MRR without pyramid algorithm, the unknown coefficients are estimated at each resolution level by combining all types of observations, and these coefficients are used to calculate the detail signals of each level. Wu et al. (2018) point out that the existing publications lack comparisons between the MRR and the single-level approach. Thus, the direct comparison of the single-level approach, the MRR without pyramid algorithm, and the MRR based on the pyramid algorithm, presented in this study fills this gap in current literature. The Reuter grid is used, which generates a homogeneous coverage of equidistributed grid points on the sphere. The total number Z of Reuter grid points on the sphere is decided by a control parameter γ , and γ + 1 denotes the number of points along the meridian (Eicker 2008). In this study, we choose the parameter γ to be equal to the maximum spectral degree n i of the expansion at each resolution level i (Wittwer 2009;Bentel et al. 2013b). In regional gravity field modeling, the computation area ∂ C , where the SRBFs are located, needs to be chosen larger than the observation area ∂ O , where the observations are given (Fig. 5), and ∂ O needs to be larger than the investigation area ∂ I , where the final gravity models are computed, i.e., ∂ C ⊃ ∂ O ⊃ ∂ I . This hierarchy is necessary to mitigate edge effects. The margin size η C,O between the computation area ∂ C and the observation area ∂ O is determined following (Lieb et al. 2016), with where ϕ max is the maximum latitude value. The margin size is influenced by the shape of the SRBFs; they become wider at the lower resolution levels (i.e., when n max is smaller, see Fig. 3), and thus, a larger margin size has to be chosen to reduce edge effects (Liu et al. 2020b). In this case, the η C,O is chosen as 0.3 • , 0.6 • , 1.2 • , 2.4 • , 4.8 • for the levels 11, 10, 9, 8, 7, correspondingly (as shown in Fig. 6). Consequently, the number of unknown coefficients is decided by the generated Reuter grid points that are located inside the computation area ∂ C of each resolution level (which amounts to K 11 = 7759, K 10 = 2238, K 9 = 856, K 8 = 326, and K 7 = 157). In  Figure 5 (black rectangle) presents the corresponding investigation area ∂ I . In case of the MRR based on the pyramid algorithm, only specific observation groups are used at the higher resolution levels, which means the involved observations do not have full coverage over the observation area ∂ O . This will cause strong edge effects at the border of the high-resolution observations in the calculated detail signals, and further contaminate the final gravity model. For example, at level 11, the terrestrial data are available only in the onshore area. Strong edge effects in G 11 thus show up at the border of the terrestrial data coverage, i.e., near the coastal lines (as will be shown in Sect. 4.3). As pointed out by Lieb (2017), it is a challenging task to handle the edge effects properly in case of MRR, and it is one of the major difficulties in the practical realization of the pyramid algorithm. To address this issue, we develop a strategy in this study to reduce the edge effects in the calculated detail signals of each level. Besides the observation area ∂ O (Fig. 5) and the investigation area ∂ I (black rectangle in Fig. 5) for the whole study area, we also define ∂ O i and ∂ I i for each resolution level i when calculating the detail signals. To be more specific, ∂ O i depends on the data coverage of the observation groups involved at this level, and ∂ I i is adapted to ∂ I i = ∂ O i ∩ ∂ I . As an example, ∂ O 11 is defined as the onshore areas in Fig. 5 since only terrestrial observations are involved at level 11. Consequently, ∂ I 11 (see Fig. 7) is adapted to the onshore areas within the investigation area ∂ I . The detail signals G i of level i are then calculated within ∂ I i . Fig. 6 The estimated scaling coefficients at level 11 (first row), and levels 10, 9, 8, 7 (second to fifth row). From the second to the fifth row, the left column represents the coefficients d i|i+1 estimated directly from the pyramid algorithm, and the right column represents the updated coefficients d i after including the new data at this level (following the procedure explained in Sect. 3). The black rectangle inside each plot shows the observation area ∂ O  Fig. 7 The detail signal G i of the MRR based on the pyramid algorithm before (left column) and after (mid column) adapting the investigation area ∂ I i , as well as the modeled gravity signal (in terms of disturbing potential) F I =F + I i=7 G i (right column) from level 11 (first row) to level 7 (fifth row), withF (right column, last row) modeled from GECO

Modeling results and discussions
The estimated coefficients at each level are plotted in Fig. 6, and the black rectangle inside each plot represents the observation area ∂ O (Fig. 5). From the first row (level 11) to the fifth row (level 7), the margin size η C,O increases and the number of unknown coefficients decreases. The top row shows the scaling coefficients d 11 at level 11, which are estimated from the terrestrial observations only. Comparing the plot of d 11 to Fig. 5, it is clear that the values of the coefficients are almost zero in the area where no terrestrial data exist, and larger absolute values are observed where these data are available. This shows that additional gravity information with respect to the background model is only captured at the locations with terrestrial data coverage, which is reasonable (Liu et al. 2020b). The left plot at the second row shows the scaling coefficients d 10|11 calculated directly from the pyramid algorithm, and the right one displays the updated coefficients d 10 after including the airborne data at level 10. Comparing the two plots at level 10, we can see that the airborne observations insert additional gravity information in the region where they are located. They fill parts of the data gaps from the terrestrial data, and the gravity signals captured from the previous level (level 11) are preserved at the same time. At level 9 (third row), the same behavior as at level 10 can be observed, i.e., the altimetry data fill data gaps at the left bottom corner of the observation area, and meanwhile, the gravity signals from level 10 are kept. At level 8 (fourth row) and level 7 (fifth row), the satellite data which have even distributions are used. It can be observed that the right plots show darker colors (larger values) than the left ones inside the observation area ∂ O , which indicate the contribution of the GOCE and GRACE data at level 8 and level 7, respectively. Figure 7 (left column) visualizes the detail signals G i of the MRR based on the pyramid algorithm at each level, which show the spectral information contained in the corresponding frequency ranges (see Fig. 2). At level 11, the detail signal G 11 captures gravity information only in onshore area where the terrestrial data are located. However, at the border of the terrestrial data, strong edge effects show up. The same is valid in the detail signal G 10 of level 10, large edge effects occur at the border of the data coverage, i.e., in the coastal area of the Tyrrhenian Sea. We thus apply the strategy explained in Sect. 4.2 to reduce these edge effects by defining ∂ I i , i.e., by setting ∂ I 11 = ∂ O 11 ∩∂ I and ∂ I 10 = ∂ O 10 ∩∂ I . After level 9, the involved observations have full coverage over the observation area ∂ O (i.e., ∂ I i = ∂ I ), and no edge effects are visible within ∂ I . Correspondingly, the new detail signals of level 11 and 10 after adapting ∂ I i are presented in Fig. 7 (mid column), and the edge effects are significantly reduced. Figure 7 (right column) shows the gravity signal F I =F + I i=7 G i of each level. The gravity signal (in terms of disturbing potential) of level 6 (F, right column, last row) is the long wavelength component from the global gravity model GECO, which only contains very smooth information. When the resolution level increases from level 7 (fifth row) to level 11 (first row), more and more fine structures show up in both the detail signals and the gravity signals.
The final modeling result ( Fig. 7 right column, first row) is evaluated by the validation data, and their difference is visualized in Fig. 8d. For comparison, the single-level model, the MRR without pyramid algorithm, and the MRR based on pyramid algorithm before adapting the ∂ I i are also computed, and their differences to the validation data are shown in Fig. 8a, b, and c, respectively. The corresponding statistics are listed in Table 2. As shown in Fig. 8a, the single-level model delivers small differences compared to the validation data in onshore regions, where the terrestrial data are available. However, in offshore regions with no terrestrial data coverage, the differences are quite large. This result demonstrates that the single-level approach majorly recovers gravity information from the terrestrial observations, and the contribution of other measurements which are sensitive to lower spectral bands is not captured sufficiently. It indicates that the single-level model is not able to benefit from all the observation types, as mentioned in the Introduction, and the MRR is necessary, especially in cases where the terrestrial data do not have large coverage over the study area. After applying the MRR, the differences w.r.t. the validation data in offshore regions decrease (see Fig. 8b), which indicates that the gravity information in lower-resolution observations are better extracted. The MRR without pyramid algorithm delivers an RMS error of 4.21 m 2 /s 2 , which is 23% smaller than the one given by the single-level approach. However, Fig. 8b still shows the same pattern as the single-level model, i.e., larger differences show up in the offshore regions compared to the onshore regions. Although the MRR (without pyramid algorithm) already gives better results than the single-level approach, it is not optimal. Indeed, Lieb (2017) points out that the detail signals of different levels become correlated when all the observation groups are used at each level and recommends the implementation of the pyramid algorithm as further development of the MRR approach.
In case of the MRR based on the pyramid algorithm, as shown in Fig. 8d, the differences between the calculated gravity model and the validation data do not show dependency on the distribution of certain types of observations, i.e., the offshore area does not show larger differences than the area with terrestrial data. It suggests that each observation type makes a contribution to the final result, and the MRR benefits from all the measurements. This statement is supported by the fact that the RMS value delivered by the MRR based on the pyramid algorithm decreases by 50% w.r.t. the one delivered by the single-level approach, and 35% w.r.t. that of the MRR without pyramid algorithm. The comparison between Fig. 8d and Fig. 8c shows the benefit of applying the proposed strategy for reducing edge effects, i.e., by adapting the investigation area ∂ I i at each level. The edge effects at the border (outside the coverage) of the terrestrial data are significantly reduced in Fig. 8d. The improvement achieved by applying this strategy is 21% in terms of RMS, w.r.t. the validation data. However, at the border (inside the coverage) of the terrestrial observations, the edge effects remain the same after adapting ∂ I i , as also shown in the calculated detail signals at level 11 and level 10 (Fig. 7, mid column). Thus, a main challenge to be faced in future studies regarding the MRR based on the pyramid algorithm is to further reduce these edge effects.

Data and model configuration
After testing the performance of the MRR based on the pyramid algorithm with simulated data successfully, we apply this method to real data sets. Figure 9 shows the study area, located between 5.9 • and 14.3 • longitude and between 53.2 • and 55.3 • latitude, covering Northern Germany, parts of the North Sea and the Baltic Sea, and a small part of the Netherlands and Denmark. The gravity observations are taken from Lieb et al. (2016), where a detailed data description can be found. The terrestrial data (yellow dots in Fig. 9) are provided by the Federal States Schleswig-Holstein, Mecklenburg West-Pomerania, and Lower Saxony, with 23,465 observations covering Northern Germany. This high-resolution data set is given in terms of absolute gravity g and used in terms of gravity anomalies g. Two airborne data sets (orange flight tracks in Fig. 9) are provided by the Federal Agency of Cartography and Geodesy (BKG), one over the North Sea, and another one over the Baltic Sea, collected in 2007/2008 and 2006, respectively. They have been preprocessed and are provided in terms of gravity disturbance δg; the data accuracy is estimated by BKG to be 2 to 3 mGal. Fig. 9 The observation area ∂ O and the distribution of real gravity data, including terrestrial (yellow dots), shipborne (red dots), airborne (orange dots), and altimetry (green dots) data. The black rectangle represents the investigation area ∂ I , where the final gravity model is calculated The flight campaign over the North Sea contains 5,651 observations, with an average flight altitude of 592 m; the flight campaign over the Baltic Sea contains 6,508 observations, with an average flight altitude of 832 m. The 1,183 shipborne measurements (red dots in Fig. 9) in the Baltic Sea are preprocessed and provided by the Federal State Mecklenburg West-Pomerania in terms of gravity anomalies g. Satellite altimetry data (green dots in Fig. 9) are provided by the DGFI-TUM altimetry group in terms of geoid height N , which is derived from the measured sea surface heights (SSH) and the instantaneous dynamic ocean topography (iDOT) (Bosch et al. 2013). The altimetry data include measurements from multiple altimetry missions, namely the ERS-1e/f (Geodetic Mission phase, 1994-1995, Envisat (Extended Mission phase, 2010-2012, cycles Nr. 096-113), Jason-1 (Geodetic Mission phase, 2012-2013, cycles Nr. 500-537), and Cryosat (2010-2013, cycles Nr. 011-035), with an average spatial resolution of 10 km in the North Sea. In the Baltic Sea, we use a sparse altimetry data distribution deliberately in order to test our approach also in areas with poor data coverage. For each altimetry measurement, corrections derived from a multi-mission cross-calibration (Bosch et al. 2014) have been applied. The satellite data from GOCE and GRACE are included in this study case as the satellite-only global gravity model (GGM), instead of direct observations. Lieb et al. (2016) point out that the study area is too small (especially in the north-south direction) in comparison with the spatial resolution of the satellite data, and thus, does not allow resolving reliable long-wavelength information. In this way, we present the performance of the MRR based on the pyramid algorithm in two different scenarios of using the satellite data, i.e., as direct observations and as a satellite-only GGM, in the simulated (Sect. 4) and real case, respectively. The highest resolution level of the MRR is chosen as I = 12 according to the spatial resolution of the data, and only the terrestrial data are used at the highest level to calculate the unknown coefficients of this level and to start the pyramid algorithm. The shipborne and the airborne data are included at level 11 and level 10, respectively, and the altimetry data are added at level 9. The long-wavelength component up to level 8 (degree 255) is modeled by using the RCR procedure with GOCO06s (Kvas et al. 2021), which enhances an optimal combination of the GOCE and GRACE data. Same as in the simulated case (see Sect. 4), we also calculate the single-level model and the MRR without the pyramid algorithm for comparison. The number of Reuter grid points at each level is determined in the same manner as for the simulated study case, and the margin size η C,O between the computation area ∂ C and the observation area ∂ O is determined following Eq. (26). In this study case, the η C,O is chosen as 0.15 • , 0.3 • , 0.6 • , 1.2 • for levels 12, 11, 10, 9, correspondingly. Consequently, the number of unknown coefficients amounts to K 12 = 6638, K 11 = 2111, K 10 = 684, and K 9 = 315 at each resolution level. The margin size η O,I between the observation area ∂ O and the investigation area ∂ I (black rectangle in Fig. 9) is chosen as the median of the applied η C,O , i.e., η O,I = 0.45 • . As discussed in Sect. 4.2, strong edge effects show up in the calculated detail signals of the high resolution levels, due to the fact that the involved observations at these levels do not have full coverage over ∂ O . Again, we apply the strategy proposed in Sect. 4.2, and define ∂ I i = ∂ O i ∩ ∂ I for calculating the detail signals at each level. In the following, we always refer the MRR based on the pyramid algorithm to the one after adapting ∂ I i .

Modeling results and discussions
The elements of the estimated coefficient vectors d i and their standard deviations at each level i are plotted in Fig. 10; the black rectangle inside each plot represents the observation area ∂ O . As we can see, the estimated coefficients at level 12 only contain additional gravity information from the terrestrial data. Consequently, their standard deviations are much larger in regions without terrestrial observations, which is reasonable. At level 11, level 10, and level 9, the shipborne, Fig. 10 The estimated scaling coefficients (left column) and their standard deviations (right column) from level 12 (first row) to level 9 (fourth row). The black rectangle inside each plot represents the observation area ∂ O airborne, and altimetry data insert additional information, respectively. The new observations introduced at the lower levels fill the data gaps from the terrestrial observations, and the features from the highest level are preserved at the same time. The behavior of the standard deviations coincides with the coefficients, i.e., from level 12 to level 10, the standard deviations at the locations with shipborne and airborne data decrease, and from level 10 to level 9, those within the altimetry data coverage decrease.
The computed regional quasi-geoid models from the single-level approach, the MRR without pyramid algorithm, and the MRR based on the pyramid algorithm (after adapting ∂ I i ) are firstly validated using GPS/leveling data in Northern Germany (Gruber et al. 2011). 53 data points are available within the investigation area ∂ I (black rectan-gle in Fig. 9), which are derived from GPS-based ellipsoidal heights and leveling-based normal heights. The differences between the computed height anomaly results and those from GPS/leveling are shown in Fig. 11, and the corresponding statistics are listed in Table 3. The mean value of the difference between the gravimetrically determined height anomalies and those from GPS/leveling amounts to around 33 cm. It is consistent with the value reported in Gruber et al. (2011)), and is caused by differences in the height system definition, i.e., the local normal heights in Germany refer to the vertical datum of the European Vertical Network (EUVN), which is defined as the equipotential surface of the Earth's gravity field passing through the "Normaal Amsterdams Peil" (NAP; fundamental tide gauge in Amsterdam, the Netherlands). Figure 11 shows clearly that the MRR based on Fig. 11 Differences between the calculated quasi-geoid model and the GPS/leveling data, delivered by a the single-level approach, b the MRR without pyramid algorithm, and c the MRR based on the pyramid algorithm. Note that the mean values of the differences are removed the pyramid algorithm delivers the smallest difference w.r.t. the GPS/leveling data, with an RMS error of 2.23 cm, which is 23% smaller than the one given by the MRR without pyramid algorithm, and 35% smaller than that of the single-level model. Such large improvements demonstrate the benefits of applying the MRR based on the pyramid algorithm. Figure 11a (single-level approach) and Fig. 11b (MRR without pyramid algorithm) show much larger differences w.r.t. the GPS/leveling data in comparison to Fig. 11c (MRR based on the pyramid algorithm), even in regions with very dense terrestrial data (between 11 • and 13.85 • longitude).
In the offshore area where no GPS/leveling data are available, the computed gravity models are validated with the 2 × 2 altimetric gravity anomaly grid DTU17 (Andersen and Knudsen 2019). The differences between the computed gravity anomaly results and the DTU17 are shown in Fig. 12, and the corresponding statistics are listed in Table   3. Again, the largest differences are delivered by the singlelevel approach (Fig. 12a), with an RMS error of 7.22 mGal. In the single-level model, much smaller differences show up in regions with shipborne data coverage (see Fig. 9), which suggests that it majorly recovers gravity information from the high-resolution shipborne data, and information from other measurement types are not captured sufficiently. This result agrees with the conclusion drawn from the simulated study case (see Sect. 4.3). In comparison to the single-level model, applying the MRR (without pyramid algorithm) improves the modeling results by 13%. Comparing Fig. 12b with Fig. 12a reveals significant improvements in regions where the altimetry data are located (between 6.35 • and 7 • longitude). This indicates that the MRR extracts gravity information from the lower resolution altimetry data better than the single-level approach. However, it can still be seen from Fig. 12b that larger differences occur in regions without shipborne data Fig. 12 Differences between the calculated gravity anomaly results and the DTU17 grid in the offshore area, delivered by a the single-level approach, b the MRR without pyramid algorithm, and c the MRR based on the pyramid algorithm coverage. The differences w.r.t. the DTU17 are significantly reduced when the MRR based on the pyramid algorithm is applied, giving an RMS of 2.67 mGal, which is 57% smaller than the one delivered by the MRR without pyramid algorithm, and 63% smaller than that of the single-level approach. The improvement achieved by applying the MRR based on the pyramid algorithm is larger in the offshore area than in the onshore area, where high-resolution terrestrial data are available. It demonstrates that the single-level approach cannot represent the lower-resolution data in an optimal way, and the MRR based on the pyramid algorithm is beneficial, especially in regions where high-resolution gravity data are not available. In the MRR based on the pyramid algorithm (Fig. 12c), larger differences again occur at the border of the higher-resolution gravity data (shipborne and airborne), which is caused by edge effects, as discussed in Sect. 4.3.
It is worth mentioning that validation in the offshore area was also made w.r.t. the NKG2015 gravimetric quasi-geoid model (Ågren et al. 2016) to rule out possible conflicts when validating our results with external models also based on satellite altimetry data as the DTU17. The differences between the calculated quasi-geoid results and the NKG2015 model show the same pattern as those w.r.t. the DTU17, i.e., Fig. 12. Thus, this comparison is not shown in detail here due to the length of the manuscript. In comparison to the NKG2015, the improvement achieved by applying the developed MRR scheme based on the pyramid algorithm is 39% w.r.t. the MRR without pyramid algorithm, and it reaches 55% w.r.t. the single-level approach, in terms of RMS value.

Conclusion and outlook
This study focuses on the spectral combination of different types of gravity observations through the MRR based on the pyramid algorithm, which is realized successfully for the first time in connection with sequential parameter estimation, in regional gravity field modeling. We address in this paper the challenges regarding the practical realization of the pyramid algorithm. Firstly, a successive low-pass filtering for transforming the estimated coefficients of the highest resolution level to lower resolution levels is proposed. Furthermore, we develop an innovative MRR scheme where the coefficients for calculating the detail signals at each lower resolution level are not only determined from the pyramid algorithm but also updated by a direct combination with observation groups included according to their spectral resolution. The main contribution of our approach is that the final gravity model is able to benefit from the individual strengths of each observation type. The settings of the MRR are characterized, including the type of the SRBFs, the location of the SRBFs, and the margin size applied at each level. When the MRR based on the pyramid algorithm is applied, only specific data sets are used at the higher resolution levels, resulting in strong edge effects at the border of the high-resolution observations. Therefore, we additionally develop a suitable strategy to adapt the investigation area ∂ I i at each level i according to the coverage of involved observations. This is also a remarkable innovation of our approach as it reduces edge effects in the calculated detail signals (and therefore, the final gravity model) significantly.
The performance of the MRR based on the pyramid algorithm is evaluated in comparison to the conventional single-level approach and the MRR without pyramid algorithm using both simulated and real gravity data. In the simulated study case, the RMS obtained by the MRR based on the pyramid algorithm is 50% smaller than the one delivered by the single-level approach, and 35% smaller than the one given by the MRR without pyramid algorithm, comparing to the validation data. Moreover, the single-level model shows very large differences to the validation data in offshore regions, which indicates that the contribution of other types of observations is not captured as sufficiently as the terrestrial observations. This shows that the single-level approach is biased towards the terrestrial measurements, and the gravity information from measurements with medium to low resolution is not extracted sufficiently. Applying the MRR (without pyramid algorithm) improves the modeling results in offshore regions. However, it still shows larger differences w.r.t. the validation data in the offshore area in comparison to the onshore area. Thus, it is important and beneficial to apply the MRR based on the pyramid algorithm, especially when the high-resolution terrestrial data do not have full coverage over the study area.
In the real data case, the MRR based on the pyramid algorithm is applied for regional gravity field modeling in Northern Germany. The terrestrial, shipborne, airborne, and altimetry observations are used at the resolution levels 12, 11, 10, and 9, respectively. At the lower levels, the gravity information obtained from the highest level is preserved, and at the same time, the new observations introduced at this specific level contribute with additional information and fill data gaps. Such features in the estimated coefficients are observed in both the simulated and real data cases. The computed gravity models are validated in terms of quasigeoid and gravity anomaly with GPS/leveling data in land and the DTU17 in offshore area, respectively. In comparison to the single-level approach, the improvement achieved by the MRR based on the pyramid algorithm is 35% and 63% in terms of RMS value, w.r.t. the GPS/leveling data and DTU17, respectively. Compared with the MRR without pyramid algorithm, the RMS error obtained by the MRR based on the pyramid algorithm decreases by 23% and 57% w.r.t. the GPS/leveling data and DTU17, respectively. Such significant improvements further demonstrate the benefits of the MRR based on the pyramid algorithm. Results in the real case again show that both the single-level approach and the MRR without pyramid algorithm cannot recover gravity information from the lower-resolution observations as sufficient as the MRR based on the pyramid algorithm.
In both the simulated and the real cases, larger differences w.r.t. the validation data in the MRR based on the pyramid algorithm occur at the border of the high-resolution data, due to edge effects. Thus, a major concern for future work is to develop strategies for further reducing the edge effects in the calculated detail signals. We also plan to include the full variance-covariance matrix of the GGM. In current studies, the weight matrix of the prior information (used at the highest level) is defined as the identity matrix, which is computationally easy. However, better and more realistic modeling results might be obtained by considering the full covariance matrix of the GGM instead of a simple identity matrix. In addition, it is also planned for future work to use real GOCE gravity gradients and K-band range-rate (KBRR) data from GRACE as direct observations in our developed MRR scheme.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.