CC BY
8
Russian Journal of Earth Sciences
Downward Continuation of Airborne Gravimetry Data by Means of Spherical Radial Basis Functions
Leyla S. Sugaipova1,2 and Yury M. Neyman*'1,2
1 Moscow State University of Geodesy and Cartography, Moscow, Russia 2Center of Geodesy, Cartography and SDI, Moscow, Russia
Received 21 June 2021; accepted 21 July 2022; published 10 August 2022.
The problem of downward continuation of airborne gravimetry data is discussed. Use of spherical radial basis functions (SRBF) to solve this ill-posed problem is proposed. Gravity disturbances observed at flight height are continued downward to disturbing potential. The SRBF method is numerically tested using synthesised data for flight heights 2000 m, 4600 m and 6000 m and grid steps 1 arcmin and 2.5 arcmin in area bounded by colatitudes 40°, 43° and longitudes 153°, 157° (spherical coordinates). The experiments prove that the SRBF method can provide stable and accurate results. Moreover, as a result of this procedure one have an approximator in the form of a linear combination of SRBF which allows to determine the values of different transforms of potential by applying the corresponding operators to this expression.
Keywords: airborne gravimetry, disturbing potential, downward continuation, gravity disturbance, spherical radial basis functions.
Citation: Sugaipova, L. S., and Yury M. Neyman, (2022), Downward Continuation of Airborne Gravimetry Data by Means of Spherical Radial Basis Functions, Russian Journal of Earth Sciences, Vol. 22, ES4004, doi: 10.2205/2022ES000804.
1 Introduction
The main problem of physical geodesy is to determine disturbing gravity potential of Earth. It requires gravimetric data at different scales covering homogenously the whole planet. Traditional terrestrial gravimetry data covers gravity field band from medium to higher frequencies. But there are many places that are difficult to access or even unaccessible for it. Dedicated satellite missions provide measurements in lower part of needed bandwidth over Earth except the poles. Airborne gravimetry proved to be fast and economically efficient method to collect homogeneous gravimetric measurements in medium to high bandwidth that can be used also in remote areas.
One of the key steps in airborne gravimetry data processing is its downward continuation. Many publications are dedicated to this problem. Usually gravimetric data continued down to ground, geoid or some mean surface using inverse Poisson integral [Alberts and Klees, 2004; Kingdon and Vanicek, 2011; Novak and Heck, 2002; Novak et al., 2003; Martinec, 1996; Liu et al., 2017; Forsberg et al., 2000; Forsberg, 2003], to name a few. Since it is an ill-posed problem some kind of reducing its effect on results is needed. There were proposed different solution methods such as iteration pro-
'Corresponding author: yuney@miigaik.ru
cedures, least squares collocation, regularization and so on. In [Müller et al., 2008] space localizing spline functions is used to downward continue airborne gravimetry and gradiometry data. This task can be solved also by means of radial basis functions [Li et al., 2022].
Since, as a result, we are interested in derivation of the disturbing potential at ground, ellipsoid or some mean topographical surface, we propose here to use spherical radial basis functions for downward continuation of gravity disturbances Sg from flight level to disturbing potential T at ground level. Potential is smoother then gravity disturbances, so this procedure may provide better results in terms of stability and accuracy [Novak and Heck, 2002].
Usually the remove-calculate-restore technique is used for processing airborne gravimetry data. Gravity measurements collected at flight height are band-limited to around L » 4000 [Novak and Heck, 2002] in terms of harmonic analysis. By means of one of the modern satellite Earth gravity models (EGM) we can remove low-frequency part (l < 360) of the data before processing.
Now, let's formulate the problem we are going to face here. We will consider it in spherical approximation. The band-limited reduced for topographic and atmospheric effects gravity disturbances Sgb are given at flight height. The problem is to deter-
mine the band-limited disturbing gravity potential Tb at some set of points. More rigorously [Sansd, 1995; Novak and Heck, 2002]:
AT(r,0,X) = 0for r>R, ôg(r,0,X) = -dT{r:9,X) for r = R + H,
dr
(1)
T(r,0X) = O(^) for r
^ CO,
where R is the radius of the geocentric sphere (spherical approximation of geoid); H is a known constant flight height over sphere R; (r,0,X) are spherical coordinates. The last equation in (1) reminds us that as 5gb are reduced by satellite model to degree l we cannot recover the lower part of Tb.
The harmonic band-limited disturbing gravity potential Tb can be expanded into series of spherical harmonics on the surface of the sphere of radius R as follows
u GM L
T (R,0,X) = GRMY_Tn(0,X),
n=l
f (P) « £Xj O (P,Qj ),
j=1
The equations relating the measurements and the coefficients x are
Ax = u + v; A = (a,j) = (O(P„ Qj));
i = 1,2,...,n; j = 1,2,...,k; k < n,
where u is the observation vector, v is the vector of corrections to observations, Pi are the measurement nodes.
SRBF are functions of the spherical distance ty between points P and Q. Spectral properties of SRBF are defined by the non-negative numerical sequence (n) called a symbol, and the parameter j called a scale:
R2 \n+1 rPrQI
c j r2 \n+1
Oj(P,Q) = Yj2n + 1) - Qj(n)Pn(cos #
n=0 VprQ'
(4)
(2)
where Tn (6,X) are the Laplace spherical harmonics of degree n and GM denotes the product of the gravitational constant and the Earth's mass.
2 Theoretical Basics
We propose here to use spherical radial basis functions (SRBF) to solve the downward continuation problem for airborne gravimetric data.
In [Neyman and Sugaipova, 2016; Sugaipova, 2018; Neyman et al., 2021] it is shown, that approximation by means of SRBF of a required function
f (P):
here rp and rQ are radial distances of points P and Q respectively.
The symbol and the scale of SRBF allow to take into account the structure of the field under consideration in the studied local area. More information on SRBF, regular grids and methods for optimal poles selection can be found, for example, in [Michel, 2013; Reuter, 1982; Neyman et al., 2021; Sugaipova, 2018].
In practice, we will limit the summation in (4) to j = 180/A, where A is the average value of the step of the grid with the measurements.
If the function under study is the gravity potential or any of its transformants, it is advisable to associate the symbol (n) with the spectrum of this function, in particular, with degree variances ffn of its spherical harmonic expansion
llOjl 1^2 = ^^ (2n +1)x
n=0
$j(n) = c2 ^ Qj(n) =
(3)
can be used as convenient substitute for integration of this function over local region.
In (3) O (P,Qj) is SRBF, Qj are so called SRBF poles, Xj stands for weight coefficients to be determined.
In the simplest case, the number of poles coincides with the number n of data nodes. But if n is too large it can lead to situation when we have to solve large system of linear equations with an ill-conditioned matrix of coefficients. Therefore, it is recommended to construct some regular grid of poles in the area of interest. In this case, the poles may not coincide with the data nodes, and their number k may be noticeably less than the points with measurements. Furthermore, if we select the most optimal poles, it can even more reduce number of used ones, and, consequently, number of equations.
V2n + 1
Thus, the SRBF (4) for our problem takes the form
\n
\rPrQ /
J
Oj (P,Q)=Y^ V2
n + 1an
n=0
i n+1
■ Pn(cos # (5)
Degree variances can be calculated directly using one of the EGM, or by means of the degree variance models. One of the most recent degree variance models can be found in [Sansd and Sideris, 2013; Rexer and Hirt, 2015]:
2 = / GM \2 A • Bn ;
°n =\(n - 1)(n - 2)(n + 4)(n + 17); A = 5.0 x 10-8; B = 0.999845, n ^ 3.
Gravity signal strengths at the surface of the topography is well approximated by degree variance model [Rexer and Hirt, 2015]:
2
1.79 x 10-7 • 0.999995
n
ff2 = l M _
n \ R ) (n - 1)(n- 2)(n + 4)(n + 17)'
Its use can be recommended when real data are continued down to the topography.
n
3 Numerical tests and results
Simulated data were used in order to test accuracy and stability of approximation by means of SRBF described in the previous section. Test area is limited by colatitudes 0m¡n = 40°, 0max = 43° and longitudes Am¡n = 153.5°, Amax = 156.5° (spherical coordinates). The inner zone bounded by colatitudes 0 = 40.5°, 0 = 42.5° and longitudes A = 153°, A = 157° was also considered in order to check if there are any edge effects. Reference system used is WGS84.
Earth gravity model XGM2019e [Zingerle et al., 2020] complete to degree 5540 and order 5480 were used to synthesize model data. Two grids of band-limited gravity disturbances
— M %L ■ í R \n+2
bgb (R + H,0,A) = -^Yjn +1) ( RRH) Tn(0, A)
n=l
at height H = 4600 m were calculated: 1) with grid step A = 2.5' and band limits l = 241, L = 4320; 2) with grid step A = 1' and band limits l = 241, L = 5400.
Band-limited disturbing potential Tb (2), generated from XGM2019e at height H = 0 m and the same grid spacings A and band limits l, L as for gravity disturbance were used as reference data for comparing results of downward continuation.
The numerical tests consisted of two steps. At the first step the initial data - gravity disturbances bg at flight height H = 4600 m - were approximated by linear combination of SRBF:
1 _
j(P, Q) = V2n + 1an • Pn(cos^).
n=0
]
Lgr = (180 • 60)/A,
Output of the first stage is the set of coordinates of SRBF poles RP and their weight coefficients x. This set serves as input of the second stage when approximating the potential T on the sphere R (H = 0 m).
In this case, equation (5) takes a form
o T
J , / — + — \n + 1
(P,Q) = ^V2n +1 0n(-R-) • Pn(cos^), (7)
n=0
where the sphere of radius R + H plays now role of reference sphere, and degree variances correspond to potential T:
rn = (C2 + Sn) • I
GM \2 4 4
- m /c .
\ — /
This expression is obtained from (5) with R = rP = tq = R + H.
Degree variances of bg were derived from XGM2019e with appropriate dimensioning:
an = (C„2 + S2)(n + 1)2 ■ (GM)2 X 1010mGal2,
where C2, S2 are the EGM spherical harmonic coefficients. It is taken into account that 1 mGal = 10-5 m/c2.
As already noted in previous section, the poles may not coincide with the observation nodes. It can be recommended to use some regular grid for their location. The Reuter's grid [Reuter, 1982] was used here. Grid density Lgr, namely, the number of parallels included in it, was selected depending on the spacing between the observations
The SRBF 0;T (P, Q) (7) is plotted in Figure 1 as a function of spherical harmonics degree n for different values of flight height H: 2000, 4600, 6000, 8000 and 10000 m. The value of spherical distance is fixed to ty = 1°. It can be clearly seen from the graph that (P, Q) has the divergent nature for L ^ to. With increasing elevation H the oscillations become more severe. It can mean that with an increase in H the procedure of downward continuation of gravitational information using the SRBF (7) can lose stability and lead to a deterioration of the results. But for the considered values H = 2000, 4600 and 6000 m, as it will be seen later from the computational experiments, the accuracy of field recovery remains quite high.
The graph in the Figure 2 shows the SRBF (P, Q) as a function of the spherical distance ty for l = 241 and L = 4320.
Tables 1, 2 and 3 provide calculation statistics for the residual Tres potential field, which is the difference between the reference potential values T and the approximation result Tcal
cal
of gravity distur-
cal
(6)
where the numerator and the denominator are dimensioned in arcmin. Selection of the most optimal poles was not performed.
Tres = T - T
and for the residual field bances
bgres = bg-
with bg denoting reference values, and bgcal standing for approximated ones.
The minimum, maximum, mean and standard deviation values of Tres and bgres are given in the Table 1 for flight level H = 4600 m. In the column with number of points, in parentheses the number of poles used is shown.
Recall that at the first stage the gravity disturbance bg is approximated and its output sets of poles RP and coefficients x are used as input for the second stage to approximate T. Thus, between x, on the one hand, and T, on the other hand, there is only an indirect relationship, since the measure
1500 2000 2500 3000 3500 4000
Degree of the series expansion, n
Figure 1: Behavior of SRBF as a function of harmonic series expansion degree n for different flight heights H.
Figure 2: SRBF as a function of spherical distance ty for l = 241 and L = 4320.
for calculation of x is the accuracy of recovering of Sg, and not T.
As can be seen from Table 1, approximation by SRBF demonstrates high accuracy of recovery even for the whole test area but for inner zone the results are an order of magnitude better.
To obtain a more complete picture regarding the SRBF method, additional calculations were made.
1. In order to determine whether approximators constructed for both T and Sg can be generalized, 100 points were randomly selected in the test area; 49 of them fell into the inner zone (Figure 3). In these points values of T and Sg were calculated at the heights H = 0 m and H = 4600 m, respectively. The statistics of results are displayed in the Table 2. Two
computation options were considered. In one case, sets of poles RP and coefficients x obtained for grid spacing A = 2.5' were used, and in the other - for A = 1'. It is natural to expect that the use of a set corresponding to a smaller step should give better accuracy, since more detailed observation information is used. It is the case for Sg, but calculations for the potential T show no improvement. Perhaps the reasons for this are the smoother nature of the potential and the fact mentioned above that there is only indirect connection between RP, x and T.
2. In order to determine the robustness of the estimates, the approximation was made at different heights, namely H = 2000 m and H =
Table 1: Statistics of approximation results for bg and T at a height level H = 4600 m
Number Frequency Step min max mean std
Region of points band A
(number l -L of poles)
bgres (mGal)
7081 241-4320 2.5' -0.4083 0.4951 -0.0002 0.0251 (4647)
40°-43°, 43621 241-5400 1' -0.1660 0.1622 -3.0e-05 0.0050 153°-157° (28770)
43621 241-5400 1' -1.4061 0.9187 0.0003 0.0624 (7241)
3577 241-4320 2.5' — — — 0.0090 (4647)
40.5°-42.5°, 21901 241-5400 1' — — — 0.0024 153.5°-156.5° (28770)
21901 241-5400 1' — — — 0.0265 (7241)
Tres (m2/ c2)
7081 241-4320 2.5' -1.4423 1.7844 0.0078 0.2026 (4647)
40°-43°, 43621 241-5400 1' -3.3988 4.9013 0.0045 0.8899 153°-157° (28770)
43621 241-5400 1' -1.8811 2.0591 0.0079 0.2002 (7241)
3577 241-4320 2.5' -0.0595 0.1068 0.0004 0.0182 (4647)
40.5°-42.5°, 21901 241-5400 1' -1.8121 2.3117 -0.0478 0.7807 153.5°-156.5° (28770)
21901 241-5400 1' -0.0885 0.1884 0.0009 0.0246 (7241)
6000 m. A grid with step A = 2.5' was considered. The Table 2 provides the statistics of results. In general, it can be seen that the signal weakens with the height, which is natural to expect, but the accuracy of the approximation does not deteriorate.
3. In order to see the influence of number of used SRBF poles two scenarios of tests for A = 1' were performed. In one case the grid density Lgr was used as in (6) and it led to 28,770 poles in the area of interest (see Table 1). In another case we took 2 times less value of Lgr providing 7241 poles. Comparing results one can see that more poles lead to increase in accuracy of recovering of bg by an order of magnitude. But for disturbing potential T accuracy of results deteriorated. One reason for this is the influence of instability of the down-
ward continuation problem. Further, recall that approximation of gravity disturbances is performed in the plane of observations. In contrast, approximation of disturbing potential is carried out in the orthogonal plane and, consequently, has a nature of extrapolation. This fact can be assumed among other reasons of degradation of approximation accuracy for T.
4 Conclusions
We discussed here the determination of the disturbing gravity potential on the surface of the sphere of radius R (spherical approximation of the geoid) from band-limited airborne gravity data at flight level. Simulated gravity data, generated
Figure 3: Location of test points (red). The points used in the calculation of x are shown in blue dots.
The inner zone is limited by black rectangle.
Table 2: Statistics of approximation results for Sg and T on the set of test points
Region Number of points Frequency band l - L RPand X for A min max mean std
8gres (mGal)
40°-43°, 100 241-4320 2.5' -0.3138 0.1921 -0.0088 0.0662
153°-157° 241-5400 1' -0.2013 0.1299 0.0002 0.0441
40.5°-42.5°, 49 241-4320 2.5' -0.1341 0.0910 -0.0048 0.0415
153.5°-156.5° 241-5400 1' -0.0467 0.0360 -0.0039 0.0197
Tres (m2/ c2)
40°-43°, 100 241-4320 2.5' -0.8263 0.5797 0.0030 0.1492
153°-157° 241-5400 1' -0.8409 0.5860 0.0051 0.1567
40.5°-42.5°, 49 241-4320 2.5' -0.0232 0.0045 0.0010 0.0177
153.5°-156.5° 241-5400 1' -0.0491 0.0440 0.0010 0.0178
from the high-frequency EGM, were used for numerical testing of the proposed method. These tests have shown that approximation by SRBF provides highly accurate and stable results of downward continuation of gravity signal from flight levels at least as high as H = 6000 m.
The SRBF method provides an approximating construction as a linear combination of SRBF (3) not only for the initial field of gravity disturbance bg, but also for the disturbing potential T. On the one hand, it allows to calculate the values of these functions at any points in the approximation area,
Table 3: Statistics of approximation results for bg and T at different flight levels
Region Height, m min max mean std bgres (mGal)
40°-43°, 2000 -0.4883 0.5922 -0.0003 0.0310
153°-157° 4600 -0.4083 0.4951 -0.0002 0.0251
6000 -0.3710 0.4499 -0.0002 0.0226
40.5°-42.5°, 2000 — — — 0.0111
153.5°-156.5° 4600 — — — 0.0090
6000 — — — 0.0081
Tres (m2/ c2)
40°-43°, 2000 -2.1812 2.5425 0.0097 0.2473
153°-157° 4600 -1.4423 1.7844 0.0078 0.2026
6000 -1.6990 1.9237 0.0064 0.1805
40.5°-42.5°, 2000 -0.0664 0.1113 0.0005 0.0199
153.5°-156.5° 4600 -0.0595 0.1068 0.0004 0.0182
6000 -0.0688 0.1179 0.0004 0.0186
on the other hand, it becomes possible to determine the values of different transforms of potential by applying the corresponding operators to the expression (3).
Further studies would be advisable on the next questions:
• to research influence of noise on the accuracy and stability and its propagation to the results of approximation and downward continuation by using simulated data with added random noise or real data;
• to investigate influence of a number and location of the SRBF poles on the results of approximation and downward continuation. Methods for optimal poles selection developed, for example, in [Neyman et al., 2021; Sugaipova, 2018] can be recommended for this task.
Acknowledgments. This work was partly conducted in the framework of budgetary funding, adopted by the Ministry of Science and Higher Education of the Russian Federation (project No. 0708-2020-001).
References
Alberts, B., and R. Klees, A comparison of methods for the inversion of airborne gravity data, Journal of Geodesy, 78(1-2), 55-65, doi:10.1007/s00190-003-0366-x, 2004.
Forsberg, R., Downward continuation of airborne gravity data, in Gravity and Geoid 2002, 3rd Meeting of the International Gravity and Geoid Commission (IGGC), edited by I. N. Tziavos, pp. 51-56, ZITI Editions, 2003.
Forsberg, R., A. Olesen, L. Bastos, A. Gidskehaug, U. Meyer, and L. Timmen, Airborne geoid determination, Earth, Planets and Space, 52(10), 863-866, doi:10.1186/bf03352296, 2000.
Kingdon, R., and P. Vanicek, Poisson downward continuation solution by the jacobi method, Journal of Geodetic Science, 1(1), 74-81, doi:10.2478/v10156-010-0009-0, 2011.
Li, X., J. Huang, R. Klees, et al., Characterization and stabilization of the downward continuation problem for airborne gravity data, J. Geod. 96, 96(4), doi:10.1007/s00190-022-01607-y, 2022.
Liu, X., Z. Sun, K. Xu, and M. Ouyang, Downward continuation of airborne gravimetry data based on poisson integral iteration method, Geodesy and Geodynam-ics, 8(4), 273-277, doi:10.1016/j.geog.2017.03.009, 2017.
Martinec, Z., Stability investigations of a discrete downward continuation problem for geoid determination in the canadian rocky mountains, Journal of Geodesy, 70(11), 805-828, doi:10.1007/bf00867158, 1996.
Michel, V., Lectures on Constructive Approximation. Fourier, Spline, and Wavelet Methods on the Real Line, the Sphere, and the Ball, 324 pp., Springer Science, New York, doi:10.1007/978-0-8176-8403-7, 2013.
Müller, F., T. Mayer-Gürr, and A. Makhloof, Downward Continuation of Airborne Gravimetry and Gradiom-etry Data Using Space Localizing Spline Functions, in Observing our Changing Earth, vol. 133, edited by M. G. Sideris, pp. 143-153, Springer, International Association of Geodesy Symposia, doi:10.1007/978-3-540-85426-5_17, 2008.
Neyman, Y. M., and L. S. Sugaipova, Basics of multiscale approximation of geopotential, MIIGAiK, Moscow, p. 218, (In Russian), 2016.
Neyman, Y. M., L. S. Sugaipova, and S. Patalov, From Hotine-Koch integral to approximation by means of spherical radial basis functions, Geodesy and Aeropho-tosurveying, 65(1), 5-12, doi:10.30533/0536-101x-2021-65-1-5-12, (In Russian), 2021.
Novak, P., and B. Heck, Downward continuation and geoid determination based on band-limited airborne gravity data, Journal of Geodesy, 76(5), 269-278, doi:10.1007/s00190-002-0252-y, 2002.
Novak, P., M. Kern, K.-P. Schwarz, et al., On geoid determination from airborne gravity, Journal of Geodesy, 76(9-10), 510-522, doi:10.1007/s00190-002-0284-3, 2003.
Reuter, R., Über Integralformeln der Einheitssphäre und harmonische Splinefunktionen, RWTH, Aachen, 1982.
Rexer, M., and C. Hirt, Spectral analysis of the Earth's topographic potential via 2D-DFT: a new data-based degree variance model to degree 90,000, Journal of Geodesy, 89(9), 887-909, doi:10.1007/s00190-015-0822-4, 2015.
Sanso, F., The long road from measurements to boundary value problems in physical geodesy, Manuscripta geodaetica, 20(5), 326-344, 1995.
Sanso, F., and M. G. Sideris, The local modelling of the gravity field by collocation, in Geoid Determination: Theory and Methods, edited by F. Sanso and M. G. Sideris, pp. 203-258, Springer, Berlin, Heidelberg, doi:10.1007/978-3-540-74700-0_5, 2013.
Sugaipova, L. S., Development and study of methods for multiscale modelling of geopotential, Ph.D. thesis, Moscow, (In Russian), 2018.
Zingerle, P., R. Pail, T. Gruber, and X. Oikonomidou, The combined global gravity field model XGM2019e, Journal of Geodesy, 94(7), 66, doi:10.1007/s00190-020-01398-0, 2020.
