Characterization and stabilization of the downward continuation problem for airborne gravity data

Abstract

In this study, we compare six commonly used methods for the downward continuation of airborne gravity data. We consider exact and noisy simulated data on grids and along flight trajectories and real data from the GRAV-D airborne campaign. We use simulated and real surface gravity data for validation. The methods comprise spherical harmonic analysis, least-squares collocation, residual least-squares collocation, least-squares radial basis functions, the inverse Poisson method and Moritz’s analytical downward continuation method. We show that all the methods perform similar in terms of surface gravity values. For real data, the downward continued airborne gravity values are used to compute a geoid model using a Stokes-integral-based approach. The quality of the computed geoid model is validated using high-quality GSVS17 GPS-levelling data. We show that the geoid model quality is similar for all the methods. However, the least-squares collocation approach appears to be more flexible and easier to use than the other methods provided that the optimal covariance function is found. We recommend it for the downward continuation of GRAV-D data, and other methods for second check.

Appendices

Appendix A

1.1 The derivation of the equations of downward continuation of airborne gravity data using various approaches

1.1.1 Spherical harmonic analysis (SHA)

In spherical approximation, band-limited (full/residual) gravity disturbances from airborne gravimetry can be modelled as

$$ \begin{aligned} & \delta g\left( {r~,\phi ~,\lambda ~} \right) \approx - \frac{{\partial T}}{{\partial r}} \\ & \quad = \frac{{{\text{GM}}}}{{R^{2} }}\mathop \sum \limits_{{n = n1}}^{{n2}} \left( {\frac{R}{{r~}}} \right)^{{n + 2}} \left( {n + 1} \right)\\ &\quad \times \mathop \sum \limits_{{m = 0}}^{n} \left( {\bar{c}_{{{\text{nm}}}} ~\cos \left( {m\lambda } \right) + \bar{s}_{{{\text{nm}}}} ~\sin \left( {m\lambda } \right)} \right)\bar{p}_{{{\text{nm}}}} ~\left( {\sin \phi } \right) \\ \end{aligned} $$

(13)

where $T$ is the gravity disturbing potential; ${\text{GM}}$ is the product of the Universal gravitational constant and the assigned reference mass; $R$ is the mean Earth’s radius; $\overline{p}_{{{\text{nm}}}}$ are the fully normalized associated Legendre functions of degree $n$ and order $m$; $\left\{ {\overline{c}_{{{\text{nm}}}} ,\overline{s}_{{{\text{nm}}}} } \right\}$ are the Stokes’s coefficients to be solved from the airborne gravity disturbances ${\updelta }g$; and the limits of $n1$ to $n2$ define the bandwidth of the airborne gravity disturbances.

However, solving for all coefficients $\left\{ {\overline{c}_{{{\text{nm}}}} ,\overline{s}_{{{\text{nm}}}} } \right\}$ from airborne data given inside a local area is not possible; in fact, there are infinitely many solutions to this problem. To foster one particular solution, one may assume zero gravity data outside the local area.

This treatment will definitely reduce the power in the spherical harmonic coefficients because of applying the global function for a representation of space-limited data. This set of coefficients are no longer representing the true gravity field of the Earth, especially outside the study region (Jekeli 2016). An intuitive explanation is that all the data of the world implies that $\left\{ {\overline{c}_{{{\text{nm}}}} ,\overline{s}_{{{\text{nm}}}} } \right\}$ are zeros, except only a small amount of regional data says they are not zeros. Note that the spherical harmonic functions, $\overline{p}_{{{\text{nm}}}} \left( {\sin \phi } \right)\cos \left( {m\lambda } \right)$ and $\overline{p}_{{{\text{nm}}}} \left( {\sin \phi } \right)\sin \left( {m\lambda } \right)$, are global functions; using regional data, it is impossible to estimate all spherical harmonic coefficients.

In addition to the theoretical shortcomings of using SHA to downward continue airborne data, numerically, Eq. (13) cannot be established directly on the original airborne survey point $i\left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)$ in order to utilize the EGM208 block diagonal engine to quickly solve the 4 –million spherical harmonic coefficients. In practice, the point airborne data has to be reduced from the flight altitude onto a grid on the surface of a larger (in terms of semi-major axis) reference ellipsoid that is located in the mean flight height in order to solve the spheroidal harmonic coefficients quickly by using the well-known block diagonal method (Pavlis 1998), which are then transformed into spherical harmonics by using Jekeli’s method (1988). Finally, the spherical harmonic coefficients are scaled back to the desired system. The details of spheroidal harmonic analysis and their transformations to the spherical harmonic are well-known and will not be repeated here.

In the GRAV-D airborne gravity data, the flight heights can vary by several thousand meters in a typical survey campaign. The real flights are not on a mean flight height as assumed by some algorithms. A novel iterative scheme (Smith et al. 2013) described by the following steps was developed to downward continue the real data from varying flight heights on to an inflated ellipsoid whose mean radius is $R^{{{\text{big}}}}$, and its surface is located at the numerical mean flight height. Then the spherical harmonic coefficients are scaled back by using the ratio between $R^{{{\text{big}}}}$ and $R_{{}}^{{}}$.

• (1) First, let’s assume the residual airborne gravity disturbances are given on the surface of an inflated ellipsoid and solve the spherical harmonic coefficients:

  $$ \begin{aligned}& \left\{ {\overline{c}_{{{\text{nm}}}}^{T0} ,\overline{s}_{{{\text{nm}}}}^{T0} } \right\}\\ &\quad= {\mathcal{F}}_{0}^{ - 1} \left\{ {{\updelta }\underline {g}_{0} \left( {r_{i} - h,\phi_{i} ,\lambda_{i} } \right) = {\updelta }\underline {g} \left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)} \right\}, \end{aligned}$$

  (14)

  where ${\updelta }\underline {g}_{0} \left( {r_{i} - h_{i} ,\phi_{i} ,\lambda_{i} } \right)$ are the pseudo gravity disturbances on the reference ellipsoid; ${\mathcal{F}}_{0}^{ - 1} \left\{ \cdot \right\}$ represent spherical harmonic analysis operator for the block diagonal method (for simplicity of description the steps to transform the ellipsoid harmonics into spherical harmonics are neglected, but it is applied in the computation);

• (2) Then, the following two terms are generated for the next iteration.

The first term is called the misclosure:

$$ {\updelta }\underline {g}^{M1} = {\mathcal{F}}_{0}^{{\left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{T0} ,\overline{s}_{nm}^{T0} } \right\} - {\updelta }\underline {g} , $$

(15)

where ${\mathcal{F}}_{0}^{{}} \left\{ \cdot \right\}$ is the spherical harmonic synthesis operator (SHS) at the original airborne survey points at the real altitudes based on $\left\{ {\overline{c}_{{{\text{nm}}}}^{T0} ,\overline{s}_{{{\text{nm}}}}^{T0} } \right\}$. Basically, this term shows the errors in the spherical harmonic coefficients due to the assumption made in step (1).

The second term is the DC term:

$$ {\updelta }\underline {g}^{{{\text{DC1}}}} = {\mathcal{F}}_{0}^{{\left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{T0} ,\overline{s}_{{{\text{nm}}}}^{T0} } \right\} - {\mathcal{F}}_{0}^{{\left( {r_{i} - h_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{\delta 0} ,\overline{s}_{{{\text{nm}}}}^{\delta 0} } \right\} $$

(16)

• (3) Based on the DC term, a new observation as shown in Eq. (17) is assumed to be given on the surface of the bigger ellipsoidal surface for the next round of SHA as shown in Eq. (18).

  $$ {\updelta }\underline {g}_{1} \left( {r_{i} - h_{i} ,\phi_{i} ,\lambda_{i} } \right) = {\updelta }\underline {g} \left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right) - {\updelta }\underline {g}^{{{\text{DC1}}}} , $$

  (17)
  $$ \left\{ {\overline{c}_{{{\text{nm}}}}^{T1} ,\overline{s}_{{{\text{nm}}}}^{T1} } \right\} = {\mathcal{F}}_{1}^{ - 1} \left\{ {{\updelta }\underline {g}_{1} \left( {r_{i} - h,\phi_{i} ,\lambda_{i} } \right)} \right\}, $$

  (18)

Another misclosure term is obtained by:

$$ {\updelta }\underline {g}^{M2} = {\mathcal{F}}_{1}^{{\left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{T1} ,\overline{s}_{{{\text{nm}}}}^{T1} } \right\} - {\updelta }\underline {g} , $$

(19)

• (4) In this step, the DC term of the misclosure term is computed according to:

  $$ \left\{ {\overline{c}_{{{\text{nm}}}}^{T2} ,\overline{s}_{{{\text{nm}}}}^{T2} } \right\} = {\mathcal{F}}_{2}^{ - 1} \left\{ {{\updelta }\underline {g}^{M2} } \right\}, $$

  (20)
  $$\begin{aligned} {\updelta }\underline{\Delta g}^{{{\text{DC2}}}}& = {\mathcal{F}}_{2}^{{\left( {r_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{T2} ,\overline{s}_{{{\text{nm}}}}^{T2} } \right\} \\& \quad - {\mathcal{F}}_{2}^{{\left( {r_{i} - h_{i} ,\phi_{i} ,\lambda_{i} } \right)}} \left\{ {\overline{c}_{{{\text{nm}}}}^{T2} ,\overline{s}_{{{\text{nm}}}}^{T2} } \right\}, \end{aligned}$$

  (21)

• (5) Finally, the original airborne gravity is reduced from the flight trajectories onto the surface of the inflated ellipsoid by:

  $$ {\updelta }\underline {g}_{3} \left( {r_{i} - h_{i} ,\phi_{i} ,\lambda_{i} } \right) = {\updelta }\underline {g}_{1} - \left( {{\updelta }\underline {g}^{M2} - {\updelta }\underline {g}^{{{\text{DC2}}}} } \right), $$

  (22)

Finally, the spherical harmonic coefficients for the original airborne data are given as

$$ \left\{ {\overline{c}_{{{\text{nm}}}}^{T} ,\overline{s}_{{{\text{nm}}}}^{T} } \right\} = {\mathcal{F}}_{3}^{ - 1} \left\{ {{\updelta }\underline {g}_{3} } \right\}, $$

(23)

Note that a gridding step is still needed to interpolate the scattered airborne gravity data onto the nodes of a grid (one data grid and one DEM grid) to start the iterations described above.

1.1.2 Least squares collocation (LSC) method

Forsberg’s (1987) covariance model is based on a power spectral model for a planar residual gravity field, matching the Kaula rule up in the main band frequency band up to a “depth” parameter D (equivalent to the use of a Bjerhammar sphere, but twice the magnitude). This gives rise to a gravity covariance function ${C}_{\delta g,\delta g}$ of shape − log(z + r), where z = z₁ + z₂ + D is a height parameter, and r the slant range between the two covariance points. To avoid singularities in the geoid terms for a fully self-consistent model, several log-terms are needed to be combined through a deeper “compensation depth” T, giving a relatively complicated expression of form

$$ \begin{aligned} & C_{{\delta g,\delta g}} \left\{ {P\left( {\phi _{i} ,\lambda _{i} ,h_{i} } \right),Q\left( {\phi _{j} ,\lambda _{j} ,h_{j} } \right)} \right\}\\ &\quad = \frac{{\sigma ^{2} }}{{3\log \left( {D + T} \right) - 3\log \left( {D + 2T} \right) + \log \left( {D + 3T} \right) - \log \left( D \right)}}\\ &\qquad \times \left\{ \begin{gathered} 3\log \left( {D + T + h_{i} + h_{j} + \sqrt {s + \left( {D + T + h_{i} + h_{j} } \right)^{2} } } \right) \hfill \\ - 3\log \left( {D + 2T + h_{i} + h_{j} + \sqrt {s + \left( {D + 2T + h_{i} + h_{j} } \right)^{2} } } \right) \hfill \\ + \log \left( {D + 3T + h_{i} + h_{j} + \sqrt {s + \left( {D + 3T + h_{i} + h_{j} } \right)^{2} } } \right) \hfill \\ - \log \left( {D + h_{i} + h_{j} + \sqrt {s + \left( {D + h_{i} + h_{j} } \right)^{2} } } \right) \hfill \\ \end{gathered} \right\} \end{aligned} $$

(24)

where $\sigma_{{}}^{2}$ is the data variance; $D$ is the depth to the high frequency attenuation layer; $T$ is the thickness to bottom layer, low frequency attenuation; and $s = \sqrt {x_{{}}^{2} + y_{{}}^{2} }$ the planar distance (approximated here with $x = 111.11{\text{km}} \times \left( {\lambda_{i} - \lambda_{j} } \right)\cos \overline{\phi }$ and $y = 111.11{\text{km}} \times \left( {\phi_{i} - \phi_{j} } \right)\cos \overline{\phi }$; $h_{i}^{{}}$ and $h_{j}^{{}}$ are the flight height in km).

The covariance between gravity disturbance, $\delta g,$ and height anomaly, ${\upzeta }$, based on the covariance model of gravity anomalies is (Forsberg, 1987):

$$ \begin{gathered} C_{\delta g,\zeta } \left\{ {P\left( {\phi_{i} ,\lambda_{i} ,h_{i} } \right),Q\left( {\phi_{j} ,\lambda_{j} ,h_{j} } \right)} \right\} = - \frac{1}{\gamma }\frac{{\sigma^{2} }}{{3\log \left( {D + T} \right) - 3\log \left( {D + 2T} \right) + \log \left( {D + 3T} \right) - \log \left( D \right)}} \hfill \\ \quad \times\left\{ \begin{gathered} \left( {\sqrt {s + \left( {D + h_{i} + h_{j} } \right)^{2} } - \left( {D + h_{i} + h_{j}^{{}} } \right)\log \left( {D + h_{i} + h_{j} + \sqrt {s + \left( {D + h_{i} + h_{j} } \right)^{2} } } \right)} \right) \hfill \\ - 3\left( {\sqrt {s + \left( {D + T + h_{i} + h_{j} } \right)^{2} } - \left( {D + T + h_{i} + h_{j} } \right)\log \left( {D + T + h_{i} + h_{j} + \sqrt {s + \left( {D + T + h_{i} + h_{j} } \right)^{2} } } \right)} \right) \hfill \\ + 3\left( {\sqrt {s + \left( {D + 2T + h_{i} + h_{j} } \right)^{2} } - \left( {D + 2T + h_{i} + h_{j} } \right)\log \left( {D + 2T + h_{i} + h_{j} + \sqrt {s + \left( {D + 2T + h_{i} + h_{j} } \right)^{2} } } \right)} \right) \hfill \\ - \left( {\sqrt {s + \left( {D + 3T + h_{i} + h_{j} } \right)^{2} } - \left( {D + 3T + h_{i} + h_{j} } \right)\log \left( {D + 3T + h_{i} + h_{j} + \sqrt {s + \left( {D + 3T + h_{i} + h_{j} } \right)^{2} } } \right)} \right) \hfill \\ \end{gathered} \right\} \hfill \\ \end{gathered} $$

(25)

Given the covariance function, the prediction of gravity anomalies and height anomalies including error variances using LSC can be implemented by the well-known collocation equations:

$$ \widehat{\delta g}^{{{\text{pred}}}} = C_{{\delta g^{{{\text{pred}}}} ,\delta g^{{{\text{obs}}}} }} \left( {C_{{\delta g^{{{\text{obs}}}} ,\delta g^{{{\text{obs}}}} }} + D_{n,n} } \right)^{ - 1} \left[ {\begin{array}{*{20}c} {\delta g_{1}^{{{\text{obs}}}} } \\ {\begin{array}{*{20}c} {\delta g_{2}^{{{\text{obs}}}} } \\ \vdots \\ \vdots \\ \end{array} } \\ {\delta g_{n}^{{{\text{obs}}}} } \\ \end{array} } \right], $$

(26)

$$ \begin{aligned} D\{ \widehat{\delta g}^{{{\text{pred}}}} \} &= C_{{\delta g,\underline{\delta g} }} - C_{{\delta g^{{{\text{pred}}}} ,\delta g^{{{\text{obs}}}} }} \left( {C_{{\delta g^{{{\text{obs}}}} ,\delta g^{{{\text{obs}}}} }} + D_{n,n} } \right)^{ - 1}\\ & \quad C_{{\delta g^{{{\text{pred}}}} ,\underline{\delta g}^{{{\text{obs}}}} }}^{T} ,\end{aligned} $$

(27)

$$ {\hat{\zeta }}^{{{\text{pred}}}} = C_{{{\upzeta }^{{{\text{pred}}}} ,\underline{\delta g}^{{{\text{obs}}}} }} \left( {C_{{\delta g^{{{\text{obs}}}} ,\delta g^{{{\text{obs}}}} }} + D_{n,n} } \right)^{ - 1} \left[ {\begin{array}{*{20}c} {\delta g_{1}^{{{\text{obs}}}} } \\ {\begin{array}{*{20}c} {\delta g_{2}^{{{\text{obs}}}} } \\ \vdots \\ \vdots \\ \end{array} } \\ {\delta g_{n}^{{{\text{obs}}}} } \\ \end{array} } \right], $$

(28)

$$ \begin{aligned} D\{ {\hat{\zeta }}^{{{\text{pred}}}} \} &= C_{{{\upzeta },{\upzeta }}} - C_{{{\upzeta }^{{{\text{pred}}}} ,\underline{\delta g}^{{{\text{obs}}}} }} \left( {C_{{\underline{\delta g}^{{{\text{obs}}}} ,\underline{\delta g}^{{{\text{obs}}}} }} + D_{n,n} } \right)^{ - 1} \\& \quad C_{{{\upzeta }_{{}}^{{{\text{pred}}}} ,\underline{\delta g}^{{{\text{obs}}}} }}^{T} ,\end{aligned} $$

(29)

where $C_{{{\upzeta },{\upzeta }}} = \frac{1}{{\gamma^{2} }}\left[ {\frac{3}{4}zr + \left( {\frac{1}{4}r^{2} - \frac{3}{4}z^{2} } \right)\log \left( {z + r} \right)} \right].$

In the planar approximation, the covariance model is described only by three parameters, i.e., gravity variance $\sigma_{{}}^{2}$, high-frequency attenuation depth parameter $D$, and low-frequency attenuation depth parameter $T$. Due to its simplicity, it is widely used in processing airborne gravity data (Forsberg et al. 1996, 1999, 2000, 2003). However, the simplified covariance function cannot always yield a good fit of the data in the remove-compute-restore scenario, especially for the geoid part, as the covariance functions have difficulties to approximate the negative covariances coming from removing a reference field (which again is a consequence of the planar approximation).

1.1.3 A.3 Residual LSC

Residual least-squares collocation is a modified form of LSC (Willberg et al. 2019) in the context of remove-compute-restore (RCR). In RLSC, the covariance modeling is taken from the global gravity field modeling (GGM) and forward topographic modeling instead of empirical parameterization. The general form of RLSC reads as (Willberg et al. 2019):

$$ {\mathbf{s}} = \left( {{\mathbf{C}}_{{\hat{s}\hat{l}}}^{m} + {\mathbf{C}}_{{\hat{s}\hat{l}}}^{t} } \right)\left( {{\mathbf{C}}_{ll} + {\mathbf{C}}_{{\hat{l}\hat{l}}}^{m} + {\mathbf{C}}_{{\hat{l}\hat{l}}}^{t} } \right)^{ - 1} \left( {{\mathbf{l}} - {\hat{\mathbf{l}}}^{m} - {\hat{\mathbf{l}}}^{t} } \right) + {\hat{\mathbf{s}}}^{m} + {\hat{\mathbf{s}}}^{t} , $$

(30)

where ${\mathbf{s}}$ is the signal, ${\mathbf{l}}$ is the observable, and the super-scripts $m$ and $t$ are representing information from GGM and topography, respectively. The circumflexes indicates the corresponding quantities are a-priori estimates from GGM and topographic modeling. ${\hat{\mathbf{l}}}^{m} , {\hat{\mathbf{l}}}^{t}$ are linked to the remove step, ${\hat{\mathbf{s}}}^{m} ,{\hat{\mathbf{s}}}^{t}$ are linked to the restore step.

From Eq. (30), we can see that RLSC uses separated error covariance matrices ${\mathbf{C}}$ to describe signal correlations from every input source. In particular, ${\mathbf{C}}_{{\widehat{ \cdot }\widehat{ \cdot }}}^{m}$ includes the full covariance information from a high-resolution GGM, which is not included in LSC. Accordingly, it considers accuracies and correlations of the high-resolution GGM and models the gravity field of the Earth as anisotropic and location-dependent (up to the maximum degree of the GGM). This more realistic representation of the gravity field should generally improve the quality of a regional gravity field model.

A more detailed differentiation between LSC and RLSC is explained in Willberg et al. (2019), which also highlights corresponding similarities to specific realizations of LSC that already include accuracy information from a GGM in different forms (e.g., Haagmans & van Gelderen, 1991; Pail et al. 2010).

1.1.4 Radial basis function (RBF)

The (residual) gravity disturbing potential is written as

$$ \begin{aligned}\delta T\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)& = \frac{{{\text{GM}}}}{{R_{{}}^{{}} }}\mathop \sum \limits_{q = 1}^{{Q_{L} }} \alpha_{q} \mathop \sum \limits_{n = n1}^{n2} \left( {\frac{R}{{r_{p} }}} \right)^{n + 1}\\ & \quad \times b_{n} \sqrt {2n + 1} \overline{p}_{n} \left( {\frac{{\vec{r}_{p} \cdot \vec{r}_{q} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q} } \right|}}} \right)\end{aligned} $$

(31)

where $b_{n}$ is a set of parameters that control the type of RBF (Schmidt et al. 2007), which determine the spectrum of the RBF, $\vec{r}_{q}$ denotes the location of the center of the RBF,$\overline{p}_{n}$ is the $4\pi$-normalized Legendre polynomial of degree n, and $\alpha_{q}^{{}}$ are the parameters to be estimated from the data, and $Q_{L}$ is the total number of coefficients. Then, in spherical approximation, the residual gravity disturbance is given by:

$$\begin{aligned} &\delta \delta g\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right) = - \frac{\partial \delta T}{{\partial r_{p} }} = \frac{{{\text{GM}}}}{{R^{2} }}\mathop \sum \limits_{q = 1}^{{Q_{L} }} \alpha_{q} \\& \quad \times\mathop \sum \limits_{n = n1}^{n2} \left( {\frac{R}{{r_{p} }}} \right)^{n + 2} b_{n} \left( {n + 1} \right)\sqrt[{}]{2n + 1}\overline{p}_{n} \left( {\frac{{\vec{r}_{p} \cdot \vec{r}_{q} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q} } \right|}}} \right),\end{aligned} $$

(32)

Eq. (32) provides the functional model for the estimation of the RBF parameters using ordinary least-squares techniques. Note that the computation of residual gravity disturbances at the Earth’s surface requires an extra synthesis step once the RBF parameters have been computed. Conceptually, using least-squares to estimate the RBF parameters allows the use of the original data along the flight directory without any spatial interpolation prior to the estimation of the RBF coefficients. This convenience avoids all of the complicated data pre-processing needed in steps (1) to (5) in A.1. Furthermore, the spatial concentration of a RBF allows for a parameterization of a regional gravity field requiring much less coefficients compared to the use of spherical harmonics (e.g., Klees et al. 2008; Lieb et al. 2016; Li 2018a).

The performance of the RBF least-squares approach critically depends on the definition of the RBF network, i.e., the location, including the depth (for some types of RBFs), and distribution of the RBFs. There are many ways to find a suitable RBF network (e.g., Schmidt et al. 2007, Klees et al. 2008, Eicker et al. 2013, Naeimi 2013, Lin et al. 2019, Foroghi et al. 2018). In this study, we use the nodes of a Reuter grid to position the RBFs horizontally (Eicker et al. 2013, Naeimi 2013, Lieb et al. 2016, Klees et al. 2018, Klees et al. 2019).

After setting up the system, the observation equations can be written as a Gauss-Markov Model (Schaffrin 2006):

$$ {\mathbf{E}}\left( {\mathbf{y}} \right) = A{\varvec{\xi}}, D\left( y \right) = \sigma_{0}^{2} P^{ - 1} , $$

(33)

where ${\mathbf{y}}$ is a $m \times 1$ column vector that contains the residual gravity disturbances, $A$ is a $m \times k$ design matrix with $k = Q_{L}^{{}}$ number of basis functions used, and ${\varvec{\xi}}$ is a $k \times 1$ vector comprising the RBF parameters. The regularized ordinary least-squares solution is given by (Eq. 78 Schmidt et al. 2007):

$$ \hat{\user2{\xi }} = \left( {A^{T} PA + {{\varvec{\upkappa}}}^{{{\mathbf{om}}}} {\mathbf{I}}} \right)^{ - 1} (A^{T} {\mathbf{Py}}), $$

(34)

where the optimal regularization parameter ${{\varvec{\upkappa}}}_{{{\mathbf{om}}}}$ is computed using Algorithm B in (Xu 1992) and

$$ {\varvec{D}}\{ \widehat{{{\varvec{\xi}}\} }} = \sigma_{y}^{2} \left( {A^{T} A + {{\varvec{\upkappa}}}_{{{\mathbf{om}}}} {\mathbf{I}}} \right)^{ - 1} \left( {A^{T} A} \right)\left( {A^{T} A + {{\varvec{\upkappa}}}_{{{\mathbf{om}}}} {\mathbf{I}}} \right)^{ - 1} , $$

(35)

The downward continued gravity disturbance is simply given by:

$$ {\hat{\mathbf{y}}} = A^{dc} \hat{\user2{\xi }};\user2{ }D\left\{ {{\hat{\mathbf{y}}}} \right\} = A^{dc} D\left\{ {\hat{\user2{\xi }}} \right\}A^{dcT} , $$

(36)

1.1.5 Inverse Poisson method

The formula for the upward continuation of a harmonic function $r_{p} {\updelta }g\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)$ at flight altitude P ($r_{p} ,\phi_{p} ,\lambda_{p}$) from a spherical surface ($\sigma_{R}^{{}} \left( {Q_{{}}^{^{\prime}} } \right)$) is given by the Poisson integral as (e.g., Alberts and Klees 2004):

$$\begin{aligned} &r_{p} {\updelta }g\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right) = \frac{1}{{4\pi R^{2} }}\iint\limits_{{\sigma_{R}^{{}} }} {R{\updelta }g\left( {r_{q}^{^{\prime}} ,\phi_{q}^{^{\prime}} ,\lambda_{q}^{^{\prime}} } \right)}\\ & \quad {\wp \left( {R,\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}},r_{p} } \right)d\sigma_{R} \left( {Q^{^{\prime}} } \right),} \end{aligned}$$

(37)

where $\wp \left( {R,\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}},r_{p} } \right)$ is the Poisson kernel given by

$$ \begin{aligned}&\wp \left( {R,\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}},r_{p} } \right) = \mathop \sum \limits_{n = 0}^{\infty } \left( {\frac{R}{{r_{p} }}} \right)^{n + 1} \sqrt {2n + 1} \overline{p}_{n} \left( {\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}}} \right) \\&\quad = \mathop \sum \limits_{n = 0}^{\infty } 2n + 1\left( {\frac{R}{{r_{p} }}} \right)^{n + 1} P_{n} \left( {\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}}} \right) = R\frac{{r_{p}^{2} - R^{2} }}{{l^{3} }}, \end{aligned}$$

(38)

The distances between $P\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)$ and $Q^{\prime}\left( {r_{q}^{^{\prime}} ,\phi_{q}^{^{\prime}} ,\lambda_{q}^{^{\prime}} } \right)$ is

$$ l = \sqrt {r_{p}^{2} + R^{2} - 2rR\frac{{\vec{r}_{p} \cdot \vec{r}_{q}^{^{\prime}} }}{{\left| {\vec{r}_{p} } \right|\left| {\vec{r}_{q}^{^{\prime}} } \right|}}} = \sqrt {r_{p}^{2} + R^{2} - 2\vec{r}_{p}^{{}} \cdot \vec{r}_{q}^{^{\prime}} } $$

(39)

Thus, the relationship between the observed gravity disturbances at altitude and the counterparts on the geoid in a spherical approximation is given by:

$${{\updelta }}g\left( P \right) = \frac{{R^{2} }}{{4\pi r_{p} }}\iint\limits_{{\sigma _{R} }} {{{\updelta }}g\left( {Q^{'} } \right)\frac{{r_{p}^{2} - R^{2} }}{{l^{3} }}{\rm{d\sigma }}_{R} \left( {Q^{'} } \right)}$$

(40)

Solving $\updelta g\left({Q}^{^{\prime}}\right)$ from $\updelta g\left(P\right)$ is called the DC of the gravity disturbance or the inverse Poisson problem. In general, it belongs to the Fredholm integral equation of the first kind in mathematics, which is ill-posed due to the amplification of noise when dealing with a high-resolution field. It is essential to stabilize the Poisson DC using regularization or smoothing (see e. g. Martinec 1996).

For the DC of airborne data, the inversion is done for a regional area using a regional dataset. It is well-known that the missing data outside the data area (the far-zone) has limited impact on the solution over the area of interest if the data are reduced for the contribution of a global gravity model. Hence, it is not necessary to modify the Poisson kernel in the case of the remove-compute-restore (RCR) scheme. In the case of RCR, a one-degree cap has been shown to be sufficient (Huang 2002).

After selecting an inner zone cap size, it remains to discretize the integral and solve the linear system of equations. The inner zone ($\psi <{\psi }_{0}$) contribution of the (residual) gravity disturbance at point $p$ can be expressed by numerical discretization of the integral in Eq. (40) as.

$$ {\updelta }g_{i}^{p} = \mathop \sum \limits_{j = 1}^{M} B_{ij} {\updelta }g_{j}^{q^{\prime}} $$

(41)

where

$$ \begin{aligned} & B_{ij} \\ & = \left\{ {\begin{array}{ll} {\frac{R}{{4\pi r_{i} }}\sigma_{j} \wp \left( {R,{{\cos}}(\psi_{ij} ),r_{p} } \right)} & {{\text{if}} \psi_{ij} \le \psi_{0} ;{\text{and}} i \ne j} \\ {\frac{R}{r}d^{l} \left( {r_{p} ,\psi_{0} ,R} \right) - \mathop \sum \limits_{j = 1,j \ne i}^{K} B_{ij} } & {{\text{if}} \psi_{ij} \le \psi_{0} ;{\text{and }}i = j} \\ 0 & {{\text{if}}\; \psi_{ij} > \psi_{0} } \\ \end{array} } \right.\end{aligned} $$

(42)

$$ d^{l} \left( {r_{p} ,\psi_{0} ,R} \right) = \frac{1}{2}\left[ {\frac{r + R}{r}\left( {1 - \frac{r - R}{l}} \right)} \right] $$

(43)

and the integral area on a unit sphere is ${\sigma }_{j}={\Delta }_{\phi }{\Delta }_{\lambda }\mathrm{cos}(\phi )$.

Eq. (40) can be expressed in a matrix representation as $={\varvec{g}}$, which is then solved by direct computations, or the following iterations (Huang 2002):

$$ {\varvec{f}} = A{\varvec{f}} + {\varvec{g}} $$

(44)

where $A=I-B$, I being an identity matrix. Eq. (44) shows the sums of the observed downward continued gravity disturbances and correction terms. The iteration in Eq. (44) can start from

$$ {\varvec{f}}_{{}}^{0} = {\varvec{g}} $$

(45)

Followed by, ${{\varvec{f}}}^{1}={\varvec{A}}{{\varvec{f}}}^{0}+{\varvec{g}}$, ${{\varvec{f}}}^{2}={\varvec{A}}{{\varvec{f}}}^{1}+{\varvec{g}}$,…,${{\varvec{f}}}^{{\varvec{i}}}={\varvec{A}}{{\varvec{f}}}^{{\varvec{i}}-1}+{\varvec{g}}$, and ${{\varvec{f}}}^{{\varvec{i}}+1}={\varvec{A}}{{\varvec{f}}}^{{\varvec{i}}}+{\varvec{g}}$.

The dimension of the linear system in Eq. (44) depends on the discretization step and data coverage, and can become very large. In order to avoid solving a huge system of linear equations, a block-wise approach is adopted to solve the DC for a block of 1° × 1° at a time without overlapping. No tile pattern is found in the DC results. Each side of the block is extended at least 1° in spherical distance to avoid edge effects. The iteration stops when the RMS differences between two consecutive solutions are smaller than a specified threshold value. The RMS differences correspond to the solution residuals of the previous iteration, as shown by Kingdon and Vaníček (2011). In this study, the threshold RMS difference and residual are set to 0.001 and 1.0 mGal, respectively, for the convergence of the inverse Poisson method.

1.1.6 Moritz’s analytical downward continuation

Moritz’s ADC method (Moritz 1980, Sect. 45) was originally derived to continue the gravity anomaly from the telluroid to the point level for the determination of the height anomaly by Molodensky’s theory. The DC correction is expressed as the sum of a series of terms. The G1 term represents the first order approximation, i.e., a linear gradient DC correction. When the DC of the data is well-posed, the higher-order terms attenuate rapidly to lead to a convergent correction. Otherwise the series becomes numerically divergent when the DC problem is numerically ill-posed, and must be truncated empirically to produce an approximate result.

In this study, we apply it to the DC of airborne gravity disturbance. The computational equations (Huang 2002, Sect. 5.2) are shown below.

$$ \delta g\left( {r_{g} ,\phi_{p} ,\lambda_{p} } \right) = \mathop \sum \limits_{n = 0}^{\infty } G_{n} $$

(46)

$$ G_{0} \left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right) = \delta g\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right) $$

(47)

$$ G_{n} \left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right) = - \mathop \sum \limits_{r = 1}^{n} H^{r} L_{r} \left( {G_{n - r} \left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)} \right) $$

(48)

where the $L$ operator is given by:

$$\begin{aligned}& L\left( {G\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)} \right) = \frac{{R^{2} }}{2\pi }\iint\limits_{{\sigma_{R}^{{}} }} {\frac{G\left( Q \right) - G\left( P \right)}{{(2R\sin \left( {\frac{\psi }{2}} \right))^{3} }}d\sigma_{R} \left( Q \right)} \\ & \quad {- \frac{1}{R}G\left( {r_{p} ,\phi_{p} ,\lambda_{p} } \right)} \end{aligned}$$

(49)

Empirically the DC is truncated to order 5 so that a potential divergence manifest itself, and can be controlled numerically for an ill-posed case.

Appendix B

2.1 Method for computing a geoid model from band-limited downward continued airborne gravity data

Using the remove-compute-restore (RCR) Stokes–Helmert scheme, a geoid model is computed from the gravity disturbances given on the surface of the reference ellipsoid by (Huang and Véronneau 2013):

$$ N\left( {\Omega } \right) = N_{0} \left( {\Omega } \right) + N_{1} \left( {\Omega } \right) + N_{{\text{H}}} \left( {\Omega } \right) + \delta N_{{{\text{PITE}}}} \left( {\Omega } \right) $$

(50)

where ${N}_{0}\left(\Omega \right) \mathrm{and} {N}_{1}\left(\Omega \right)$ are degree-0 and -1 terms; the last term stands for the primary indirect topographical effect on the geoid, and the Helmert co-geoid components above degree-1 can be expressed as

$$ N_{{\text{H}}} \left( {\Omega } \right) = N_{{{\text{HRef}}}} \left( {\Omega } \right) + dN\left( {\Omega } \right) $$

(51)

$$ N_{{{\text{HRef}}}} \left( {\Omega } \right) = \zeta_{{{\text{Ref}}}} \left( {\Omega } \right) - \delta N_{{{\text{DTE}}}} \left( {\Omega } \right) $$

(52)

$$ dN\left( {\Omega } \right) = \frac{R}{{4\pi \gamma \left( {\Omega } \right)}}\int\limits_{{{\upsigma }_{{0}} }} {S_{{{\text{MDB}}}} \left( \psi \right)dg_{{\text{H}}} \left( {{\Omega }^{^{\prime}} } \right)d\sigma } $$

(53)

When the residual geoid height is determined to the same spectral band as the reference field, the direct topographical effect (DTE) on the airborne data are computed to the same band as the reference field. Then we have

$$\begin{aligned} dg_{{\text{H}}} \left( {\Omega } \right) &= \left[ {\Delta g_{{{\text{Airborne}}}} \left( {\Omega } \right) - \delta \Delta g_{{{\text{DTE}}}} \left( {\Omega } \right)} \right] \\&- \left[ {\Delta g_{{{\text{Ref}}}} \left( {\Omega } \right) - \Delta g_{{{\text{DTE}}}} \left( {\Omega } \right)} \right] \end{aligned}$$

(54)

or

$$ dg_{{\text{H}}} \left( {\Omega } \right) = \Delta g_{{{\text{Airborne}}}} \left( {\Omega } \right) - \Delta g_{{{\text{Ref}}}} \left( {\Omega } \right) $$

(55)

Therefore, Eq.(53) can be re-written as

$$ \begin{aligned} dN\left( {\Omega } \right) &= d\zeta \left( {\Omega } \right) = \frac{R}{{4\pi \gamma \left( {\Omega } \right)}}\int\limits_{{{\upsigma }_{0} }} {S_{{{\text{MDB}}}} \left( \psi \right)}\\ & \quad{ \times \left[ {\Delta g_{{{\text{Airborne}}}} \left( {\Omega } \right) - \Delta g_{{{\text{Ref}}}} \left( {\Omega } \right)} \right]d\sigma } \end{aligned}$$

(56)

Considering Eq. (51–56) and the degree-0 and degree-1 terms for the geoid height are equal to the same terms for the height anomaly on these reference ellipsoid, Eq. (50) can be re-written as

$$ N\left( {\Omega } \right) = \zeta_{0} \left( {\Omega } \right) + \zeta_{1} \left( {\Omega } \right) + \zeta_{{{\text{Ref}}}} \left( {\Omega } \right) + d\zeta \left( {\Omega } \right) + C_{{\text{T}}} \left( {\Omega } \right) $$

(57)

where

$$ C_{{\text{T}}} \left( {\Omega } \right) = - \delta N_{{{\text{DTE}}}} \left( {\Omega } \right) + \delta N_{{{\text{PITE}}}} \left( {\Omega } \right) $$

(58)

The ${C}_{\text{T}}\left(\Omega \right)$ term is given in (Huang and Véronneau 2015, equations A9–10), which converts the height anomaly evaluated on the reference ellipsoid to the geoid height. Note that Sjöberg (Eq. 23, 2007) derived an equivalent equation to (10) to the 3rd order approximation, and named it the topographic bias with a negative sign.