Search
Close
Search
Advanced Search
Search Menu
AI Discovery Assistant

SUMMARY

This study examines local geopotential field modelling over a mountainous region in Tanzania using vector airborne gravity data. We use the adjusted spherical cap and rectangular harmonic analyses. Both methods are based on expansion of gravitational potential into a series of orthogonal harmonic basis functions of local support in such a way that the expansion coefficients are determined by gravity observations. All three components of gravity vector are simultaneously inverted to derive the geopotential coefficients. In order to evaluate the accuracy of the local models, independent checkpoints are selected within the study region and around its boundary and the computed gravity vectors are compared with the independent gravity observations. The results show an excellent agreement with root mean square error (RMSE) of < 1.6 mGal over the study area. On the contrary, the RMSEs of global geopotential models against the checkpoints data are 7 mGal for the models up to the maximum degree of 2190 (a resolution of ∼9.1 km) and 5 mGal to 5399 (∼3.7 km). Our local models are significantly more accurate than the state-of-the-art global models and fully exploit the airborne vector data with the measurement error of ∼1 mGal. We also present the regional models constrained only by radial (vertical) or by lateral (horizontal) gravity observations. Those models are considerably less accurate than the one from 3-D gravity data inversion. Lastly, the regional models are validated against topography data. It is found that the gravity–topography correlation is 0.8–0.9 at 100 km, 0.5 at 20 km and higher than the correlations of the global models at all frequencies. The gravity–topography admittances estimated from our regional models indicate ∼130 mGal km−1 and imply the effective density of 2500 kg m−3 for topographical mass.

Geopotential theory, Lunar and planetary geodesy and gravity, Space geodetic surveys, Africa, Inverse theory, Numerical approximations and analysis

1 INTRODUCTION

The knowledge of geopotential field variation over local areas has important implication on geophysical exploration and tectonics. For the global field modelling, spherical harmonic functions (solutions of the Laplace equation) have been naturally used as basis functions. However, such ‘global’ support functions are not appropriate for local or regional modelling demanding a high spatial resolution. To synthesize high-frequency features of the geopotential field, alternative basis functions with local support such as spherical cap harmonic (SCH) functions or rectangular harmonic (RH) functions can be considered (Haines 1985).

Haines (1985) presented the theory of SCH functions for local geomagnetic field modelling in which, the boundary value problem of the Laplace equation is solved over the area of spherical cap. Accordingly, an arbitrary harmonic function with a boundary value on a spherical cap can be represented by the series expansion in terms of spherical harmonic with real degree and integer order, which is called the SCH. Haines (1985) demonstrated that these basis functions construct the orthogonal basis over the spherical cap area. Based on the author's knowledge, the pioneering application of this theory was Haines & Newitt (1986) for local magnetic modelling over Canada. Also De Franceschi et al. (1994) demonstrated the SCH method in ionosphere modelling.

De Santis (1992) developed the Haines theory and introduced the adjusted spherical cap harmonics (ASCH). In his method, a scale function is induced on the colatitude of spherical cap in such a way that the region of spherical cap is mapped into the hemispherical region. This mapping of boundary transformed the Legendre functions with real degrees into the Legendre functions with integer degrees and thus it reduces the computational costs related to SCH. He succeeded to recover smaller wavelengths of magnetic field from scattered magnetic data in comparison with the traditional SCH approach (De Santis 1991).

De Santis & Torta (1997) further applied SCH analysis for local gravity field modelling. This model also provided some functions of local gravity field such as geoid undulation and deflection of the vertical. De Santis et al. (1999) provided a strict relation between the associated Legendre functions used in global modelling and the spherical cap basis functions of the local support. Accordingly, they elegantly deduced the relationship between global and local coefficients. Thebault et al. (2004) discussed difficulties of upward continuation in SCH and derived a rigorous relationship between global and local harmonic coefficients. They also mentioned that the SCH have some troubles in a modelling of gravity field over small areas (for those regions where the curvature of the Earth's surface can be neglected) and they do not have a desirable accuracy in a radial direction. Younis (2013), with the aid of global gravity models, produced the grid data of gravity field over the part of Germany and then transferred this grid data into spherical cap coordinate system. Consequently, he used these observations to model the local gravity field by ASCH method. Hwang et al. (2012) estimated gravity anomalies and geoid undulations in South Korea using SCH and determined the thickness of sediment in Ulleung Basin. Younis (2015) utilized terrestrial vector gravity observations in a conjunction with orthometric height to model geoid undulation in Germany using SCH method. Hwang & Chen (1997) introduced two sets of fully normalized SCH to model local sea-level variations by satellite altimetry data on China Sea area. Torta (2020) provided a detailed review of the SCH literature, underlining various strengths and weaknesses of the technique. This study includes all the works concerned with the SCH technique and its alternatives, in different applications such as the analysis of external geomagnetic field variations, estimation of ionosphere parameters, application in geodesy and regional gravity field modelling and other general topics of research. He also mentioned that the SCH method has been quite rarely studied in geodesy compared with other disciplines. So the importance of further work in this field is well acknowledged.

Another method that can be used for local geopotential field modelling is RH analysis which was first introduced by Alldredge (1981) and it has been used in geomagnetic and gravity field modelling (Alldredge 1981, 1982; Malin et al. 1996). In this method, the Laplace equation is solved in the Cartesian coordinate system which results in basis functions of trigonometric nature. Though RH is theoretically less rigorous than SCH and the edge effects exist due to the periodic assumption of potential functions, the computation of the RH (trigonometric) basis functions is more stable and efficient at high degree expansion than that of the associate Legendre function (Jiang et al. 2014).

In this study, two local geopotential models over Tanzania region are constructed by analysing vector (3-D) airborne gravity observations using adjusted spherical cap and RH expansions. All three components of gravity vector data are simultaneously inverted to determine the spectral coefficients of the regional geopotential field models. We develop an approach to reduce the edge effect by introducing auxiliary gravity data covering the extended region beyond the study area to better constrain the basis functions at the boundary. Finally, we compare the accuracy and efficiency of two geopotential models and validate them against independent data.

2 SPHERICAL CAP HARMONIC EXPANSION OF GRAVITATIONAL POTENTIAL

The external gravitational potential of the Earth satisfies the following Laplace equation in the Earth centred-earth-fixed (ECEF) coordinate system:
(1)
where ∆ is the Laplace operator and V stands for the gravitational potential function. The solution of this equation in a spherical coordinate system with radial coordinate r, co-latitude θ and longitude λ is represented based on a series expansion in terms of spherical eigen functions. In the case where the boundary conditions are prescribed over a whole sphere, those eigen functions will be degenerated into the so-called solid spherical harmonics (with integer degree and order). However, if the boundary conditions are restricted to the part of a sphere in the form of a spherical cap, those eigen functions are collapsed on the basis functions known as SCH. To this end, one may write V as a representation of the external gravitational potential in terms of SCH (Haines 1985):
(2)
where R is the radius of the reference sphere, amk and bmk are SCH coefficients, nk(m) is real degree of harmonic function and |$P_{{n_k}(m)}^m(.)$| is the Legendre function of real degree and integer order.
Fig. 1 describes the coordinate frame (Xc–Yc–Zc) used in the above SCH function in the relationship to the usual global ECEF coordinate frame (Xs–Ys–Zs). The Zc axis passes through the centre of the spherical cap region. The cap is bounded by the polar coordinate of 0≤θ≤θ0 and azimuth of 0≤λ≤2π. In the above equation, |$P_{{n_k}(m)}^m(\cos\theta )$| are the Legendre functions of real but not necessarily integer degrees. The degree of Legendre functions which is denoted here by nk(m) can be determined from the solution of the following equations defined at the cap boundary (Haines 1985):
(3)

Since the different real values of n, as roots of the above equations, depend on m, they are denoted by nk(m) in which k being an integer subscript chosen in a way to order various roots n at a given m. Here, the index k has been chosen to start at m, analogous to the case of integral n. (Haines 1985; Torta 2020).

Spherical cap coordinate system defined by the position of spherical cap poles and angle θ0.
Figure 1.

Spherical cap coordinate system defined by the position of spherical cap poles and angle θ0.

Open in new tabDownload slide
In general, determination of nk(m) from the boundary condition of eq. (3) is relatively time consuming (Younis 2015). To overcome this drawback, De Santis (1992) introduced a modification to the Haines method by the following coordinate transformation.
(4)
where S is nothing but a scale factor of π/(2θ0). In this case, spherical cap coordinates (θ,λ,r) are transformed into the hemispherical coordinate system (θ′,λ',r′) and thus the Legendre functions with non-integer degree in eq. (2) are converted into the associated Legendre functions with integer degrees and orders. Consequently, one can show that eq. (2) can be modified as (De Santis 1992).
(5)
Where Pkm(cosθ) is the Legendre function of integer degree and order, Ckm and Skm are ASCH coefficients, also nk is real number that can be defined by (De Santis & Torta 1997)
(6)
Similar to the global spherical harmonic expansion, the wavelength of the ASCH technique can be approximately computed with the index k as follows (De Santis & Torta 1997).
(7)
The series like eq. (5) is usually truncated at Kmax which defines the spatial resolution (half of the minimum wavelength) of the geopotential field model as
(8)
After Haines (1991), the power spectral density suitable for ASCH can be computed by
(9)

It should be emphasized that the above terminology is only valid for the spherical cap area that covers the polar region, namely, the z-axis of the Earth-fixed coordinate system should pass through the North Pole of the spherical cap (Haines 1985). For arbitrary spherical cap which is not necessarily located in a polar configuration, the SCH should be rotated in a proper fashion and the coordinate transformation should be applied. In this case, the new coordinate system is defined which is obtained by the position of new pole (centre of the region) and the prime meridian which simultaneously passes through the North Pole, new Pole and South Pole (Liu et al. 2011; Younis 2015).

3 RECTANGULAR HARMONIC ANALYSIS OF GRAVITATIONAL POTENTIAL

The idea of RHA was first proposed by Alldredge (1981). RH arise from the solution of the Laplace equation in a local Cartesian coordinate system. This Cartesian coordinate is here defined in the same way as the local astronomic coordinate with the origin at the centre of the study area, the x- and y-axes, respectively, in the northward and eastward directions and the z-axis is vertically downward. Thus, one may write (Alldredge 1981; Jiang et al. 2014)
(10)

In the above equation, j = q-i + 1, Dij, Eij, Fij and Gij are the RH coefficients and Lx and Ly are the lengths of the study area in north and east directions, respectively. For local gravity field modelling, the upper limit of the above series will be limited to Qmax-1. The value of Qmax depends on the size of the study area and the amount of gravity data that are used for local gravity modelling.

The spatial resolution (half of the wavelength) of the RH can be computed approximately as (Jiang et al. 2014).
(11)
The power spectra of this model σj for each j can be expressed as
(12)
where |${\bar{D}_{ij}},\,\,{\bar{E}_{ij}},\,\,{\bar{F}_{ij}}$| and |${\bar{G}_{ij}}$| are normalized RH coefficients which are simply derived by dividing the expansion coefficients by mean square integrals of trigonometric basis functions over the region of study.

4 GRAVITATIONAL POTENTIAL FIELD MODELLING OVER THE TANZANIA REGION

4.1 Airborne vector gravity data

In this section, the gravity field modelling over Tanzania is considered as a case study to show the performance of each approach for local gravity field modelling and their accuracies are compared. We use the airborne vector gravity data on mountainous area located in the north of Tanzania, within 6.6° S < ϕ <4.9° S in latitude and 35.9 ° E < λ < 38.6° E in longitude (see Fig. 2). One of the applied local coordinate system in this study is a spherical cap coordinate system with the cap radius of 1.5° and its centre located at 5.75° S and 37.25° E (see Fig. 3). The selected spherical cap encircles the whole gravity survey area. On the other hand, in the local Cartesian coordinate system used for RH modelling, the area covers the rectangle region with the length of LY = 167.47 km, and width of LX = 279.49 km (see Fig. 4).

Sander Geophysics Tanzania Project Map. The marked areas indicate observation areas surveyed by SGL where the red line shows the study area in this paper.
Figure 2.

Sander Geophysics Tanzania Project Map. The marked areas indicate observation areas surveyed by SGL where the red line shows the study area in this paper.

Open in new tabDownload slide
Geometry of spherical cap and the study area.
Figure 3.

Geometry of spherical cap and the study area.

Open in new tabDownload slide
Topography of the study area (metre). The horizontal axes are in degrees (latitude and longitude).
Figure 4.

Topography of the study area (metre). The horizontal axes are in degrees (latitude and longitude).

Open in new tabDownload slide

The difference between maximum and minimum values of (measured) gravity in this region is approximately 217 mGal likely from large variability of rough topography over the region and thus the region may be an ideal place to test for different local gravity modellings. These data were gathered by AIRGrav (Airborne Initially Referenced Gravimeter) system, which is an economical alternative to ground and ship borne surveys, designed primarily for petroleum exploration (Sander & Ferguson 2010). The surveys were conducted with 68 flight lines with the medium height of about 1170 m from the reference ellipsoid WGS84. In addition, several extra lines were surveyed as control lines (see Fig. 5).

Main flight line (black line). Control flight line (red line).
Figure 5.

Main flight line (black line). Control flight line (red line).

Open in new tabDownload slide

Finally, a total of 37 951 vector gravity points were observed across all flight lines. To validate the accuracy of the geopotential models, 500 out of 37 951 gravity observation points uniformly distributed over the whole study area are saved and considered as independent validation data sets. These validation data sets are omitted in the process of gravity field modelling and are used as checkpoints.

The airborne gravity observations are in a local navigation frame. To use them for spherical cap and RH modelling, one should consider a coordinate transformation suitable for each of gravity modelling (Nakagawa & Yukutake 1985; Jekeli 2012; Hinze et al. 2005; Fiori 2011). In this study all three components of gravity vector are simultaneously used to derive the geopotential model coefficients using ASCH and RH analysis.

4.2 Linear models of geopotential modelling

The observation equations used for determining the geopotential models from the airborne vector gravity measurements are derived by taking analytic differentiation of the model given in eqs. (2) and (10). In ASCH approach, they are given by
(13)
(14)
(15)

Where gr, gθ and gλ are, respectively, radial, meridional and longitudinal components of gravity vector in the spherical cap coordinate frame.

In the same way, for RH approach, they can be written as
(16)
(17)
(18)
where gx, gy and gz are, respectively, north, east, and vertical components of gravity vector in the local Cartesian coordinate frame. Here we introduce variables like ν = 2π/Lx, w = 2π/Ly, u = [(iν)2+(iw)2]1/2.
The above observation equations define linear relationship between the measurements (vector gravity data) and the unknown parameters (harmonic coefficients) and they can be expressed by the following parametric equations:
(19)
where l is the vector of observations, x the vector of the unknown parameters, A is the design matrix and v is the vector of random errors of observations. This equation can be solved by usual least-squares inversion.
Since the local gravity field modelling is often an ill-posed problem characterized by the instability of normal equation matrix, some regularization methods should be used to tackle the ill-conditioning nature of the normal equations (Hansen 1994; Hansen 2001). The generalized Tikhonov regularization is one of the common methods used to solve for this kind of ill-conditioned problems. In this case, the solution is sought in such a way that the following objective function is minimized.
(20)
which results in
(21)
in which |$\alpha $| is called a regularization parameter and |$||.||$| denotes the Euclidean norm. Thus, both |${{\bf \hat{x}}}$| and |${{\bf \hat{v}}}$| are functions of regularization parameter. The key idea of all regularization methods is based on the compromise between the norm of the solution |$||{{\bf \hat{x}}}||$|⁠, and the norm of the residual |$||{{\bf A\hat{x}}} - {{\bf l}}||$| using the optimum regularization parameter. If too much regularization is imposed on the solution, then it will not fit the given data properly and the residual will be too large. On the other hand, if too little regularization is imposed, then the solution will fit to the data but will be dominated by the contributions of data errors. As such the regularization parameter should be chosen with a proper scheme (Hansen 2001).

One common criterion for choosing the regularization parameter is the L-curve method. In this method, the logarithm of the norm of the solution |$\log ||{{\bf \hat{x}}}||$| is plotted against the logarithm of the norm of residuals |$\log ||{{\bf A\hat{x}}} - {{\bf l}}||$|⁠. This graph would have a L-curved shape. (Hansen 2001). Fig. 6 shows a typical of L-curve graph and the behaviour of solutions and residuals as a function of regularization parameter. The optimum choice of parameter |$\alpha $| lies on the corner point of the L-curve. The rationale behind this choice is that the corner separates the flat and vertical parts of the curve where the solution is dominated by regularization errors and perturbation errors, respectively (Hansen 2001).

The L-curve function used for finding the regularization parameter.
Figure 6.

The L-curve function used for finding the regularization parameter.

Open in new tabDownload slide

The other challenge associated with regional modelling of potential field, particularly in RH modelling, is the edge effect. This is essentially caused by no information of the potential fields beyond the survey area while the potential field variations continue across the boundary of the region. The mathematical modelling basis functions have freedom around the margin of the study areas and are completely not constrained beyond the border (Jiang et al. 2014). This effect is small and can be neglected in the ASCH because the solution of the Laplace equation is sought over the spherical cap area and it is constrained by a boundary condition at the edge (see eqs 8 and 9 in Haines 1985). The boundary condition constrains the behaviour of the basis functions and the geopotential solutions at the edge of the cap (see eq. 3). In RH approach, however, we use the solution of Laplace equation in the rectangular Cartesian coordinate system and there is no boundary condition assumed.

Various mathematical methods have been suggested to mitigate such effect like zero padding, smooth extension, symmetric extension or other algorithms (Baniamerian et al. 2017). In this study, we use a priori geopotential models such as global models to produce long wavelength gravity information outside the study region. For the RH modelling, we first extend the flight paths beyond the border of the surveyed region and use heights of the flight lines changing based on the topography (see Fig. 7). The a priori global model is evaluated along these extra paths and the computed gravity data are used in addition to the measured data to determine high resolution geopotential field models. For the ASCH method, the synthetic data are used covering the outside of the rectangular survey region but the inside of the cap (as shown in Fig. 3). We expect the accuracy of geopotential field models to be degraded near the edge of the area. When the models are validated, we quantify the difference between the model and validation data separately inside the region and the rest of area as shown in two different colours in Fig. 8.

Height model of airborne observations (green colour), and their extrapolations (red colour).
Figure 7.

Height model of airborne observations (green colour), and their extrapolations (red colour).

Open in new tabDownload slide
Distribution of checkpoints. Inner checkpoint (red) and Edge checkpoint (blue).
Figure 8.

Distribution of checkpoints. Inner checkpoint (red) and Edge checkpoint (blue).

Open in new tabDownload slide

4.3 Data reduction

Before using the data in local modelling, the low-frequency parts of gravity signal are removed from the original gravity observations, and the resulting residual gravity data are used for local modelling (Alberts 2009). By removing the long wavelength of gravity signals computed with a global model such as EGM2008, we make the airborne gravity data localized and more suitable for modelling with the local support basis functions. Accordingly, based on a remove-compute-restore procedure, the vector gravity data have been constrained in processing using EGM2008 degree and order up to 300 (equivalently, wavelengths greater than 134 km). So long wavelength gravity signals are deducted from three components of vector airborne gravity data and the high-frequency part of gravity signal is analysed in order to construct the local feature of gravity in the study area.

4.4 Spatial resolution

The spatial resolution of the model is determined by the maximum degree of the series expansion, Kmax and Qmax in eqs (2) and (10). The number of potential field coefficients is, respectively, (Kmax + 1)2 and 2Qmax(Qmax−1)+3 for ASCH and RH. The maximum degree of harmonic expansion has a great influence on the accuracy of local gravity modelling. In fact, correct choice of this parameter guarantees appropriate fitting of the models to the data as well as computational efficiency. In this study, the optimal values of Kmax and Qmax for ASCH and RH will be found based on the comparison against the validation data sets. We will make different local models with different choices of Kmax and Qmax and find the models providing smaller difference against the validation data. To this end, higher weight for the internal checkpoints is considered. Namely, if we have two candidates for each of Kmax and Qmax, we choose the one yielding smaller RMSE for those checkpoints inside the region of study (Nakagawa & Yukutake 1985; Malin et al. 1996).

5 NUMERICAL RESULTS AND DISCUSSIONS

In this section, the results of local gravity field modelling over Tanzania region are presented. First of all, the values of Kmax and Qmax should be determined. For ASCH and RH, the computed values are respectively 280 and 190, which are determined by iterative scheme discussed above. Table 1 shows the optimal values of these parameters for each method along with the final RMSEs that can be reached for local gravity field modelling using these two methods. Beyond these values of Kmax and Qmax, the RMSEs remain unchanged. In fact, up to this point, increasing the number of parameters of local geopotentials would have no significant effect on the accuracy of two models.

Table 1.
Open in new tab

RMSEs of local gravity models (mGal) result in an optimum value of Kmax and Qmax.

MethodQmax/KmaxDirectionInner checkpointsEdge checkpoints
RMSE for each componentsRMSE-totalMEAN for each componentsMEAN-totalRMSE for each componentsRMSE-totalMEAN for each componentsMEAN-total
RH190gX0.7001.4690.3001.0060.5001.4630.3000.890
gY2.1211.3921.9661.306
gZ1.2171.0051.5190.761
ASHA280gθ1.8881.5461.2881.0771.9911.5811.3891.166
gλ1.6121.2081.6011.281
gr1.0030.6000.9840.712
MethodQmax/KmaxDirectionInner checkpointsEdge checkpoints
RMSE for each componentsRMSE-totalMEAN for each componentsMEAN-totalRMSE for each componentsRMSE-totalMEAN for each componentsMEAN-total
RH190gX0.7001.4690.3001.0060.5001.4630.3000.890
gY2.1211.3921.9661.306
gZ1.2171.0051.5190.761
ASHA280gθ1.8881.5461.2881.0771.9911.5811.3891.166
gλ1.6121.2081.6011.281
gr1.0030.6000.9840.712
Table 1.
Open in new tab

RMSEs of local gravity models (mGal) result in an optimum value of Kmax and Qmax.

MethodQmax/KmaxDirectionInner checkpointsEdge checkpoints
RMSE for each componentsRMSE-totalMEAN for each componentsMEAN-totalRMSE for each componentsRMSE-totalMEAN for each componentsMEAN-total
RH190gX0.7001.4690.3001.0060.5001.4630.3000.890
gY2.1211.3921.9661.306
gZ1.2171.0051.5190.761
ASHA280gθ1.8881.5461.2881.0771.9911.5811.3891.166
gλ1.6121.2081.6011.281
gr1.0030.6000.9840.712
MethodQmax/KmaxDirectionInner checkpointsEdge checkpoints
RMSE for each componentsRMSE-totalMEAN for each componentsMEAN-totalRMSE for each componentsRMSE-totalMEAN for each componentsMEAN-total
RH190gX0.7001.4690.3001.0060.5001.4630.3000.890
gY2.1211.3921.9661.306
gZ1.2171.0051.5190.761
ASHA280gθ1.8881.5461.2881.0771.9911.5811.3891.166
gλ1.6121.2081.6011.281
gr1.0030.6000.9840.712

As it is clear in Table 1, MEAN and root-mean-square error (RMSE) are used to evaluate the spatial resolution of local geopotential models. At all checkpoints (both inner and edge checkpoints), we compute gravity vectors from two models and compare the model values with actual gravity observations to determine the error. Using this error, we compute the RMSE and mean value for each gravity components. We have three RMSE and three mean values corresponding to each gravity vector component, we also define RMSE-total and MEAN-total as the RMSE and mean values over three components. These two metrics covey the overall accuracy of two models. Accordingly, the obtained mean and RMSE values can be used to analyse the spatial resolution of the regional geopotential models. As such, it can be deduced that both local geopotential models provide the accuracy of about 2 mGal.

It should be emphasized that without padding the extended area beyond the surveyed region using the synthetic gravity data, the modelling results get considerably worsened, particularly in the z-direction of the RH model. The edge effects in all directions are significantly reduced using the synthetic data computed from a priori geopotential models. According to eqs (7) and (8), the spatial resolution of the ASCH technique in this area would be
(22)
In the same way, the spatial resolution of the RH technique by eq. (11) in this area would be
(23)

Results of local gravity field models are shown in Figs 9 and 10. In this study, 43 730 gravity vector data are used, of which 5780 are simulated data generated around the study area by the global EGM2008 model (as described in Section 4.1). Therefore, a total of 131 190 number of observation equations are produced to determine the coefficient of geopotential model using vector modelling approach. The number of unknown parameters varies based on the maximum degree Kmax and Qmax. To this end for Kmax = 280, the number of unknown parameters of ASCH would be 79 242 and for Qmax = 190 this would be 71 821 for RH. By estimating coefficients of the local geopotential models and reconstructing the gravity potential from these harmonic coefficients, the gravity vector in the study area may be predicted at arbitrary locations with desirable accuracy (see Figs 9 and 10). To show the performance of ASCH, the modelling results are only presented for the inner circle completely located inside the study area with the actual gravity data.

ASCH modelling with Kmax = 280 in spherical cap coordinate system. (a) SCH model and (b) discrepancies between predicted and observed gravity. First row, ASCH modelling of gθ. Second row, ASCH modelling of gλ. Third row, ASCH modelling of gr.
Figure 9.

ASCH modelling with Kmax = 280 in spherical cap coordinate system. (a) SCH model and (b) discrepancies between predicted and observed gravity. First row, ASCH modelling of gθ. Second row, ASCH modelling of gλ. Third row, ASCH modelling of gr.

Open in new tabDownload slide
RH modelling with Kmax = 190 in rectangular coordinate system. (a) RH model and (b) discrepancies between predicted and observed gravity. First row, RH modelling of gX. Second row, RH modelling of gY. Third row, RH modelling of gZ.
Figure 10.

RH modelling with Kmax = 190 in rectangular coordinate system. (a) RH model and (b) discrepancies between predicted and observed gravity. First row, RH modelling of gX. Second row, RH modelling of gY. Third row, RH modelling of gZ.

Open in new tabDownload slide

To show the effects of using vector gravity observations on the accuracy and performance of local modelling, further investigations are conducted here. For ASCH, we separately model the local gravity using only the radial (gr) or lateral (gθ, gλ) components of gravity vector and derive the accuracy of modelling at inner and edge checkpoints. On the other hand, for RH, we separately model the local gravity using only the vertical (gz) or horizontal (gx, gy) components of gravity vector and derive the accuracy of modelling (see Fig. 11). We explore three models (for both RH and ASCH cases) estimated from three different cases of the data set to compare how accurately one can determine a geopotential field using only either vertical or horizontal gravity data, in comparison with the solution from the 3-D vector gravity data.

Comparison of vertical, horizontal and 3-D vector local gravity modelling. ASCH (left-hand panel) and RH (right-hand panel).
Figure 11.

Comparison of vertical, horizontal and 3-D vector local gravity modelling. ASCH (left-hand panel) and RH (right-hand panel).

Open in new tabDownload slide

The results of this comparison are graphically presented in Fig. 11 (left-hand panel for ASCH and right-hand panel for RH). In Fig. 11, each vertical rectangle in left- and right-hand panels illustrates three values of RMSE related to three components of gravity vector. (Three different colours are attributed to RMSEs in three components of gravity vector; green for RMSE in |$g_r^{ASCH}$| or |$g_z^{RH}$|⁠, orange for RMSE in |$g_\lambda ^{ASCH}$| or |$g_y^{RH}$|⁠, blue for RMSE in |$g_\theta ^{ASCH}$| or |$g_x^{RH}$|⁠). Moreover, the rectangles with no hatching indicate that the observation of radial component for ASCH or vertical component for RH are used for computation of the respective geopotential model coefficients from the radial or vertical gravity observations (Model I). The rectangles with single-line hatching refer to the case when the observations of lateral components for ASCH or horizontal components for RH are used in the modelling from lateral or horizontal gravity observations (Model II). Lastly, the rectangles with double-line hatching show the results for 3-D vector gravity field modelling in two approaches (Model III). For example, in the left-hand panel, the length of green rectangle with single hatching shows the value of RMSE for radial component, when the observations of only lateral components of gravity are used for gravity field modelling in ASCH approach (RMSE of Model II in radial direction). Moreover, each bar which comprises of the three coloured rectangles shows the sum of three RMSEs at specific value of Kmax. For example, the first bar in the left-hand panel, illustrates that for Kmax = 50, the sum of RMSEs is 15 mGal when only the radial gravity observation is used (Model I) for local modelling. It should be pointed out that since the numbers of observations are different for Models I, II and III, the values of Kmax and Qmax would be inevitably different for each model.

As can be seen in left-hand panel of Fig. 11, when the gravity observations of lateral component are used to construct the local model using ASCH basis function, the RMSE for radial components becomes large and it becomes worse when the maximum degree of expansion increases. However, the RMSE for horizontal components is still good and it becomes better when Kmax increases. This happens for both inner and edge checkpoints.

The same situation is occurred when the radial component of gravity vector is used. But the RMSE values in three components as well as the divergence rate are much lower than compared to the previous case. This reveals that radial component of gravity data are more critical in regional potential modelling. The reason for this is that in spherical coordinate, the radial component of gravity is much larger than the lateral components while the data noise is similar in all components (it should be noted that in Figs 9 and 10 the high-frequency residual gravity is depicted and in this residual scale the radial gravity is not much larger than the lateral components). When all three components are simultaneously used in the inversion, the RMSE in three directions may reach to 2 mGal and RMSEs are more or less consistent for all three components, unlike the case of the inversion from radial or lateral gravity data.

Similarly, right-hand panel of Fig. 11, illustrates that the RH gravity modelling using vertical gravity observation causes to improve the accuracy of model in vertical directions as it is expected. Besides, horizontal modelling using RH technique gives rise to improvement in the accuracy of the model in horizontal directions, but the consequence of horizontal modelling for vertical direction is not desirable. However, vector modelling leads to further significant improvement in three directions.

Thus, it can be concluded that the gravity field modelling using vertical or radial component leads to better spatial resolution than the gravity field modelling using lateral or horizontal components.

The results of local gravity field modelling are also compared with the global geopotential models to illustrate the improvement achieved by our new regional models. Vector gravity accelerations computed from global geopotential models are compared with the measured airborne gravity data at the locations of checkpoints (see Table 2). The global gravity models of EGM2008, GECO and EIGEN06C4 with the maximum degree of 2190 (spatial resolution of ∼9 km) present 6–7 mGal RMS differences against the validation data, while XGM2019e with the maximum degree of 5399 (spatial resolution of ∼4 km) shows 5 mGal. It demonstrates the significant improvement made by expanding the model to higher degrees and orders and reducing the omission error. In comparison, note that our regional models (spatial resolution of ∼1.5 km) present 1–2 mGal RMS difference (see Table 1).

Table 2.
Open in new tab

Error of Global gravity models (mGal) in the inner and edge regions.

ModelMax degreeDirectionInner checkpointEdge checkpoint
RMSE for each componentsRMSE-totalRMSE for each componentsRMSE-total
EGM 20082190gθ4.1207.3964.3186.742
gλ5.9544.441
gr10.5679.901
EIGEN06C42190gθ4.1297.1504.4336.483
gλ5.9374.543
gr10.0529.263
GECO2190gθ4.2627.3884.7946.834
gλ5.9054.542
gr10.5219.824
XGM2019e5399gθ4.2875.1024.2955.467
gλ3.2914.073
gr6.9907.175
ModelMax degreeDirectionInner checkpointEdge checkpoint
RMSE for each componentsRMSE-totalRMSE for each componentsRMSE-total
EGM 20082190gθ4.1207.3964.3186.742
gλ5.9544.441
gr10.5679.901
EIGEN06C42190gθ4.1297.1504.4336.483
gλ5.9374.543
gr10.0529.263
GECO2190gθ4.2627.3884.7946.834
gλ5.9054.542
gr10.5219.824
XGM2019e5399gθ4.2875.1024.2955.467
gλ3.2914.073
gr6.9907.175
Table 2.
Open in new tab

Error of Global gravity models (mGal) in the inner and edge regions.

ModelMax degreeDirectionInner checkpointEdge checkpoint
RMSE for each componentsRMSE-totalRMSE for each componentsRMSE-total
EGM 20082190gθ4.1207.3964.3186.742
gλ5.9544.441
gr10.5679.901
EIGEN06C42190gθ4.1297.1504.4336.483
gλ5.9374.543
gr10.0529.263
GECO2190gθ4.2627.3884.7946.834
gλ5.9054.542
gr10.5219.824
XGM2019e5399gθ4.2875.1024.2955.467
gλ3.2914.073
gr6.9907.175
ModelMax degreeDirectionInner checkpointEdge checkpoint
RMSE for each componentsRMSE-totalRMSE for each componentsRMSE-total
EGM 20082190gθ4.1207.3964.3186.742
gλ5.9544.441
gr10.5679.901
EIGEN06C42190gθ4.1297.1504.4336.483
gλ5.9374.543
gr10.0529.263
GECO2190gθ4.2627.3884.7946.834
gλ5.9054.542
gr10.5219.824
XGM2019e5399gθ4.2875.1024.2955.467
gλ3.2914.073
gr6.9907.175

In order to compare the ASCH and RH modelling in terms of convergence rate, computational costs and accuracy, we examine the RMSE-total as a function of K and Q specifying the number of geopotential coefficients or model parameters (Fig. 12). Increasing the number of coefficients results in better fit to the data as reflected in smaller RMSE. The ASCH modelling shows a faster convergence, while the RH modelling needs a less number of model parameters up on the convergence. Lastly, the histograms of residuals from two models (⁠|${{\bf \hat{v}}}$| in eq. 21) are shown in Fig. 13. In both cases, the histograms generally follow the normal distribution.

RMSE-total (mGal) versus the total number of coeffcients used in each of modelling.
Figure 12.

RMSE-total (mGal) versus the total number of coeffcients used in each of modelling.

Open in new tabDownload slide
Histogram of discrepancy between the models and observations (the absicissa is the model residual and the ordinate is the frequency of residual).
Figure 13.

Histogram of discrepancy between the models and observations (the absicissa is the model residual and the ordinate is the frequency of residual).

Open in new tabDownload slide

6 VALIDATION OF LOCAL MODELS

We validate our regional models against independent data such as topography assuming higher degree and order (shorter wavelength) gravity information is primarily caused by topographic mass change (see also McGovern et al. 2002). To this end, the correlation spectrum |${C_{\Psi \Gamma }}$| between gravity model (⁠|$\Gamma $|⁠) and topography model (⁠|$\Psi $|⁠) at degree l is computed as
(24)
where |${\sigma _{\Psi \Psi }}(l)$| and |${\sigma _{\Gamma \Gamma }}(l)$| are the power spectrum of topography and gravity models, respectively, and |${\sigma _{\Psi \Gamma }}(l)$| is the cross-power spectrum between gravity and topography. They are computed as follow:
(25)
(26)
where |${\Gamma _{nm}}$| and |${\Psi _{nm}}$| are the coefficients of gravity and topography model, respectively. Using them, we also compute the gravity–topography admittance as follow:
(27)

The gravity–topography correlation and admittance spectra are computed using our regional models and also for global geopotential models of EGM2008 versus topography model DTM2006

The residual topography data of the study area are computed using the DTM2006 model from degree 300 to the maximum resolution of the model. Then the grid data are expanded in terms of two harmonic basis functions (ASCH and RH) and their corresponding expansion coefficients are computed by the least-squares analysis. Finally, the correlation and admittance spectra between local gravity model and local topography model are calculated following eqs (24) and (27), respectively. For comparative analysis, the correlation and admittance between EGM2008 and DTM2006 as global gravity and topography models are also evaluated (see Fig. 14).

Left: correlation between gravity model and topography model in function of wavelength. ASCH (green), RH (red) and GGMs (blue). Right: admittance between gravity model and topography model in function of wavelength. ASCH (green), RH (red) and GGMs (blue).
Figure 14.

Left: correlation between gravity model and topography model in function of wavelength. ASCH (green), RH (red) and GGMs (blue). Right: admittance between gravity model and topography model in function of wavelength. ASCH (green), RH (red) and GGMs (blue).

Open in new tabDownload slide

Fig. 14 shows that improved correlations are found from the local models; they are 0.8–0.9 at 100 km and 0.5 at 20 km, higher than the global correlation at all frequencies. The admittance spectra of the regional model (ASCH) present ∼130 mGal km−1 for wavelengths greater than 30 km. This implies the effective density of topographical mass is 2500 kg m−3 which is much more realistic estimate of the upper crustal density than the one from the global models.

The local gravity solutions are also compared with the global geopotential field models in terms of power spectrum in Fig. 15. There is an unusual power decrease starting at 100 km from EGM2008, which is responsible for the low admittance estimates. The local models found significantly larger gravity signals in the region than the global model estimates over the wavelengths of 100–10 km.

Comparison of degree RMS of local models with those of global models. ASCH (green), RH (red) and GGMs (blue).
Figure 15.

Comparison of degree RMS of local models with those of global models. ASCH (green), RH (red) and GGMs (blue).

Open in new tabDownload slide

7 CONCLUSIONS

In this paper, the problem of local gravity field modelling over Tanzania mountainous region is investigated using ASCH and RH analysis of vector airborne gravity data. The vector gravity observations are expressed as a truncated series of adjusted spherical or RH functions of the regional geopotential field and the resulting normal equations are resolved using regularized least-squares methods compounded of the Tikhonov regularization algorithm. The approach of vector modelling enables robust estimation of expansion coefficients with uniform accuracy over all directions. To handle the edge effect inevitable in any local modelling, the EGM2008 geopotential model is used to generate auxiliary (long wavelength) gravity data beyond the boundary of the region. This operation significantly improves the accuracy of local gravity models at the interior and the border of the region. Furthermore, the maximum degree of each model is determined using the accuracy criteria controlled by RMSE at all independent checkpoints. The values of Qmax for RH and Kmax for ASCH equal to 190 and 280, respectively, approximately corresponding to the average spatial resolution of 1.17 and 1.57 km. The following general conclusion can be made:

  1. The adopted local harmonic models (RH and ASCH) have the ability to represent efficiently local geopotential field to a high spatial resolution (like a couple of km) while satisfying the Laplace's equation. The upward and downward continuation and computation of different gravimetric quantities are straightforward.

  2. The independent data evaluation of the ASCH and RH models found the modelling error of 1–2 mGal in both the central parts and the boundary of the region. The local models present much higher accuracy and resolution than any available global models.

  3. Compared to the solutions from either radial (vertical) or lateral (horizontal) components of the vector airborne gravity observations, the 3-D vector gravity modelling always leads to a significantly better spatial resolution and a higher convergence rate.

  4. It is shown that local geopotential models present higher topography correlation in the study region than the global correction of EGM2008 over 100–5 km. In addition, through the admittance computation, the local models determine an effective density of about 2500 kg m−3 in the study area. The regional geopotential models present higher power spectrum than the global estimate from the EGM2008 model in wavelengths between 100 and 10 km.

  5. Our new regional geopotential models will be one of the most comprehensive model in Tanzania which may be useful to geophysical studies and exploration demanding a higher spatial resolution of the gravity field.

ACKNOWLEDGEMENTS

The authors would like to thank the Ministry of Energy & Minerals of the United Republic of Tanzania and Sander Geophysics Ltd. (SGL) for providing the airborne vector gravity observations that have been used in this paper. The authors also would like to appreciate from IT Service Desk of the University of Newcastle for providing the High Performance Computer (HPC) to handle the numerical analysis of this paper.

DATA AVAILABILITY STATEMENT

The data sets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Author contribution statement

The contributions of the authors are as follows Mohsen Feisi (M.F) and Mehdi Raoofian-Naeeni (M.RN) conceived of the presented idea of this manuscript and discussed it with Shin-Chan Han (SH.H) for review of the work and further developments. M.F and M.RN performed the computations such as processing and analysing the data, solving the normal equations and preparing the plots. Then SH.H verified and checked the implemented methods and made various suggestions to improve the paper. M.F and M.RN wrote the draft of the manuscript. Then SH.H checked the English and the structure of the paper and made the revision. Finally, all authors discussed the results and contributed to the final version of the manuscript.

REFERENCES

Alberts
B.A.
,
2009
.
Regional Gravity Field Modeling using Airborne Gravity Data
,
Delft, 2009. 200 pagina's. ISBN: 978 90 6132 312 9

Google Scholar

OpenURL Placeholder Text

 

Alldredge
L.
,
1981
.
Rectangular harmonic analysis applied to the geomagnetic field
,
J. geophys. Res.: Solid Earth
,
86
,
3021
–
3026
.
10.1029/JB086iB04p03021

Google Scholar

Crossref
Search ADS

 

Alldredge
L.R.
,
1982
.
Geomagnetic local and regional harmonic analyses
,
J. geophys. Res.: Solid Earth
,
87
,
1921
–
1926
.
10.1029/JB087iB03p01921

Google Scholar

Crossref
Search ADS

 

Baniamerian
J.
,
Oskooi
B.
,
Fedi
M.
,
2017
.
Source imaging of potential fields through a matrix space-domain algorithm
,
J. appl. Geophys.
,
136
,
51
–
60
.
10.1016/j.jappgeo.2016.10.035

Google Scholar

Crossref
Search ADS

 

De Franceschi
G.
,
De Santis
A.
,
Pau
S.
,
1994
.
Ionospheric mapping by regional spherical harmonic analysis: new developments
,
Adv. Space Res.
,
14
,
61
–
64
.
10.1016/0273-1177(94)90240-2

Google Scholar

Crossref
Search ADS

 

De Santis
A.
,
1991
.
Translated origin spherical cap harmonic analysis
,
Geophys. J. Int.
,
106
,
253
–
263
.
10.1111/j.1365-246X.1991.tb04615.x

Google Scholar

Crossref
Search ADS

 

De Santis
A.
,
1992
.
Conventional spherical harmonic analysis for regional modelling of the geomagnetic field
,
Geophys. Res. Lett.
,
19
,
1065
–
1067
.
10.1029/92GL01068

Google Scholar

Crossref
Search ADS

 

De Santis
A.
,
Torta
J.
,
1997
.
Spherical cap harmonic analysis: a comment on its proper use for local gravity field representation
,
J. Geod.
,
71
,
526
–
532
.
10.1007/s001900050120

Google Scholar

Crossref
Search ADS

 

De Santis
A.
,
Torta
J.
,
Lowes
F.
,
1999
.
Spherical cap harmonics revisited and their relationship to ordinary spherical harmonics
,
Phys. Chem. Earth Part A.
,
24
,
935
–
941
.
10.1016/S1464-1895(99)00138-6

Google Scholar

Crossref
Search ADS

 

Fiori
R.
,
2011
.
Application of Spherical Cap Harmonic Analysis to Plasma Convection Mapping at High Latitudes
,
University of Saskatchewan
.

Google Scholar

OpenURL Placeholder Text

 

Haines
G.
,
1985
.
Spherical cap harmonic analysis
,
J. geophys. Res.: Solid Earth
,
90
,
2583
–
2591
.
10.1029/JB090iB03p02583

Google Scholar

Crossref
Search ADS

 

Haines
G.
,
1991
.
Power spectra of sub-periodic functions
,
Phys. Earth planet. Inter.
,
65
,
231
–
247
.
10.1016/0031-9201(91)90130-A

Google Scholar

Crossref
Search ADS

 

Haines
G.
,
Newitt
L.
,
1986
.
Canadian geomagnetic reference field 1985
,
J. Geomag. Geoelectr.
,
38
,
895
–
921
.
10.5636/jgg.38.895

Google Scholar

Crossref
Search ADS

 

Hansen
P.C.
,
1994
.
Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems
,
Numer. Algor.
,
6
,
1
–
35
.
10.1007/BF02149761

Google Scholar

Crossref
Search ADS

 

Hansen
P.C.
,
2001
.
The L-curve and its use in the numerical treatment of inverse problems
, in
Computational Inverse Problems in Electrocardiology
(pp.
119
–
142
.),
WIT Press
.

Google Scholar

OpenURL Placeholder Text

 

Hinze
W.J.
et al. ,
2005
.
New standards for reducing gravity data: the North American gravity database
,
Geophysics
,
70
,
J25
–
J32
.
10.1190/1.1988183

Google Scholar

Crossref
Search ADS

 

Hwang
C.
,
Chen
S.-K.
,
1997
.
Fully normalized spherical cap harmonics: application to the analysis of sea-level data from TOPEX/POSEIDON and ERS-1
,
Geophys. J. Int.
,
129
,
450
–
460
.
10.1111/j.1365-246X.1997.tb01595.x

Google Scholar

Crossref
Search ADS

 

Hwang
J.S.
,
Han
H.-C.
,
Han
S.-C.
,
Kim
K.-O.
,
Kim
J.-H.
,
Kang
M.-H.
,
Kim
C.H.
,
2012
.
Gravity and geoid model in South Korea and its vicinity by spherical cap harmonic analysis
,
J. Geodyn.
,
53
,
27
–
33
.
10.1016/j.jog.2011.08.001

Google Scholar

Crossref
Search ADS

 

Jekeli
C.
,
2012
.
Inertial Navigation Systems with Geodetic Applications
,
edn, Vol., pp. Pages
,
Walter de Gruyter
.
http://dx.doi.org/10.1515/9783110800234
.

Google Scholar

OpenURL Placeholder Text

 

Jiang
T.
,
Li
J.
,
Dang
Y.
,
Zhang
C.
,
Wang
Z.
,
Ke
B.
,
2014
.
Regional gravity field modeling based on rectangular harmonic analysis
,
Sci. China Earth Sci.
,
57
,
1637
–
1644
.
10.1007/s11430-013-4784-1

Google Scholar

Crossref
Search ADS

 

Liu
J.
,
Chen
R.
,
Wang
Z.
,
Zhang
H.
,
2011
.
Spherical cap harmonic model for mapping and predicting regional TEC
,
GPS Solut.
,
15
,
109
–
119
.
10.1007/s10291-010-0174-8

Google Scholar

Crossref
Search ADS

 

Malin
S.
,
Düzgit
Z.
,
Baydemir
N.
,
1996
.
Rectangular harmonic analysis revisited
,
J. geophys. Res.: Solid Earth
,
101
,
28205
–
28209
.
10.1029/96JB01885

Google Scholar

Crossref
Search ADS

 

McGovern
P.J.
et al. ,
2002
.
Localized gravity/topography admittance and correlation spectra on Mars: Implications for regional and global evolution
,
Journal of Geophysical Research: Planets
,
107
(
E12
),
19
–
1-19-25
.,
doi.org/10.1029/2002JE001854
.
10.1029/2002JE001854

Google Scholar

Crossref
Search ADS

 

Nakagawa
I.
,
Yukutake
T.
,
1985
.
Rectangular harmonic analyses of geomagnetic anomalies derived from MAGSAT data over the area of the Japanese Islands
,
J. Geomag. Geoelectr.
,
37
,
957
–
977
.
10.5636/jgg.37.957

Google Scholar

Crossref
Search ADS

 

Sander
L.
,
Ferguson
S.
,
2010
.
Advances in SGL AIRGrav acquisition and processing
, in
Proceedings of the Australian Society of Exploration Geophysicists Conference, Sydney, Australia
.

Google Scholar

OpenURL Placeholder Text

 

Thébault
E.
,
Schott
J.
,
Mandea
M.
,
Hoffbeck
J.
,
2004
.
A new proposal for spherical cap harmonic modelling
,
Geophys. J. Int.
,
159
,
83
–
103
.
10.1111/j.1365-246X.2004.02361.x

Google Scholar

Crossref
Search ADS

 

Torta
J.M.
,
2020
.
Modelling by spherical cap harmonic analysis: a literature review
,
Surv. Geophys.
,
41
,
201
–
247
.
10.1007/s10712-019-09576-2

Google Scholar

Crossref
Search ADS

 

Younis
G.
,
2013
.
Regional Gravity Field Modeling with Adjusted Spherical Cap Harmonics in an Integrated Approach
,
TUprints-TU Darmstadt publication service
.

Google Scholar

OpenURL Placeholder Text

 

Younis
G.
,
2015
.
Local earth gravity/potential modeling using ASCH
,
Arabian J. Geosci.
,
8
,
8681
–
8685
.
10.1007/s12517-014-1767-2

Google Scholar

Crossref
Search ADS

 
© The Author(s) 2021. Published by Oxford University Press on behalf of The Royal Astronomical Society.
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Subject
Gravity, Geodesy and Tides
Issue Section:
Gravity, geodesy and tides
Download all slides