Blažej Bucha

Senior Lecturer at the Department of Theoretical Geodesy, Slovak University of Technology in Bratislava.

This page offers output products from some of research papers that I co-authored. All codes are written in Matlab and were tested in Matlab R2015b and R2018b.

In case you have any comments, suggestions, questions or criticism, please feel contact me at blazej.bucha@stuba.sk.

GrafLab (GRAvity Field LABoratory) — Ultra-high-degree surface and solid spherical harmonic synthesis

GrafLab is a MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders, which allows:

• evaluation of 38 functionals of the geopotential up to ultra-high degrees and orders (e.g., the geoid, height anomaly, gravity anomalies/disturbances, deflections of the vertical, gravitational tensor),
• evaluation of commission error of 26 functionals using full variance-covariance matrix of spherical harmonic coefficients,
• depiction of the computed data on a map.

The Matlab source code accompanied by a manual and test data can be downloaded here.

GrafLab can be fully controlled through GUI (see the figure below) or as a function (without GUI). The latter option is beneficial especially when GrafLab needs to be called multiple times with, for instance, different geopotential model, maximum degree or evaluation points.

Further details can be found in the paper:

Bucha, B., Janák, J., 2013. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Computers & Geosciences 56, 186-196, http://doi.org/10.1016/j.cageo.2013.03.012.

Some sample outputs from GrafLab

• Geoid undulations evaluated and decipted via GrafLab from EGM96 and DTM2006.0 up to degree 360:

This output can be reproduced by using test data and running the following command:

GrafLab('OK',3986004.415E+8,6378136.3,0,360,1,'EGM96.mat',0,0,-90,0.25,90,0,0.25,360,0,[],[],[],[],'EGM96_Geoid_to360',0,10,1,'DTM2006.mat',0,1,1,2,6,1,60,600,1)

• Besides the evaluation of gravity field quantities such as the geoid, gravity anomalies, deflections of the vertical and many others, GrafLab can also be used to synthesize, for instance, the Earth's topography (look for the text "Some useful tips and tricks" in the list_of_changes_in_GrafLab.txt file):

This output can be reproduced by using test data and running the following command:

GrafLab('OK',1,1,0,360,[0 0 0 0 0],'DTM2006.mat',1,0,-90,0.25,90,0,0.25,360,0,[],[],[],[],'Topography',0,11,1,[],0,1,1,2,6,1,60,600,1)

• In fact, this command can be used to perform any surface spherical harmonic synthesis of the form $$f(\varphi,\lambda) = \sum_{n=0}^{n_{\max}} \sum_{m=0}^{n} \left( \bar{C}_{nm}\, \cos(m\, \lambda) + \bar{S}_{nm}\, \sin(m\, \lambda) \right) \bar{P}_{nm}(\sin\varphi)\,,$$ $$f(\varphi,\lambda) = \sum_{n=0}^{n_{\max}} \sum_{m=0}^{n} \Bigl( \bar{C}_{nm}\, \cos(m\, \lambda)$$ $$+ \bar{S}_{nm}\, \sin(m\, \lambda) \Bigr) \, \bar{P}_{nm}(\sin\varphi)$$ where $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are $4\pi$-fully-normalized (real) surface spherical harmonic coefficients of the function $f$, $n$ and $m$ are spherical harmonic degree and order, respectively, $\bar{P}_{nm}(\sin\varphi)$ are the $4\pi$-fully-normalized (real) associated Legendre functions of the first-kind, and, finally, $\varphi$ and $\lambda$ are the spherical latitude and longitude, respectively. This means that GrafLab can synthesize a wide range of (real) functions given on a sphere up to ultra-high harmonic degrees (tens of thousands and beyond), provided that the respective coefficients $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are available.

• After some simple tricks (look for the text "Some useful tips and tricks" in the list_of_changes_in_GrafLab.txt file), GrafLab can be used to synthesize gravity field of other celestial objects, such as the Moon, Mars and others. Below is shown an example: the gravitational vector (in the local north-oriented reference frame) induced by topographic masses of the Moon. The evaluation points reside on a Brillouin sphere (out of all gravitating masses).

The gravitational vector was computed up to degree 2160 at 5 arc-min grid from the model "STU_Moon_topography_to720_gravity_to2160.gfc" (see the section on divergence-free gravity models of the Moon's gravity field below) using the following command:

GrafLab('OK',0.4902805823000e+13,0.1737999e+07,0,2160,[0 0 0 0 0],'STU_Moon_topography_to720_gravity_to2160.gfc',1,0,-90,0.08,90,0,0.08,360,11000,[],[],[],[],'Moon_gravitational_vector',0,16,3,[],0,1,1,2,6,1,60,600,1)

Computation time

Below are shown computation times required by GrafLab to synthesise the Earth's topography as represented by Earth2014.TBI2014.degree10800 up to harmonic degrees nmax = 90, 180, 360, 720, 1080, 2160, 4320, 6480, 8640 and 10800. The position of evaluation points is defined by the Gauss—Legendre quadrature (see Section Ultra-high-degree surface spherical harmonic analysis). The dimensions of the grids are (nmax+1)*(2*nmax+2) in terms of the latitude and the longitude, respectively, and Legendre's functions were evaluated after Fukushima (2012). The tests were conducted on a PC with Intel® CoreTM i7-6800K CPU, 128 GB of RAM (no more than a few GBs were needed in this example though) and a 250 GB SSD drive. The same grids and PC were used to demonstrate the computational times required by the software package to perform spherical harmonic analysis in Section Ultra-high-degree surface spherical harmonic analysis, so that the two figures can be mutually compared.

Go to top

isGrafLab (Irregular Surface GRAvity Field LABoratory) — Efficient ultra-high-degree solid spherical harmonic synthesis at regular grids residing on irregular surfaces

isGrafLab is a modified version of GrafLab that allows accurate and fast computation of functionals of the geopotential at dense grids on irregular surfaces such as the Earth's surface. It employs the highly efficient lumped coefficients approach for surface spherical harmonic synthesis of the functionals and their radial derivatives at regular surfaces (sphere or ellipsoid of revolution), and the Taylor series expansions to continue the functionals to the irregular surface. All the other options available in GrafLab, such as the possibility to use three different approaches to compute the fully normalized associated Legendre functions, the graphical user interface or data plotting, are also available in the new software.

The Matlab source code accompanied by a manual and test data can be downloaded here.

Further details can be found in the paper:

Bucha, B., Janák, J., 2014. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: Efficient computation at irregular surfaces. Computers & Geosciences 66, 219-227, http://doi.org/10.1016/j.cageo.2014.02.005.

Go to top

GOCE-only global gravity field models based on the kinematic orbit

The GOCE-only global spherical harmonic gravity field models models developed by Bucha et al. (2015) can be downloaded here (~1.6 MB).

The models were estimated by inverting the kinematic orbit of the GOCE satellite (time period from 01 November 2009 to 11 January 2010 which is equal to the R1 period used by the official ESA's models). Some further specifications:

• the models were originally derived by means of spherical radial basis functions (Shannon SRBF) and were transformed into spherical harmonics after the evaluation; the transformation is exact without any approximations,
• the acceleration approach was used to invert the kinematic orbit into gravity field models,
• the problematic near-zonal coefficients (polar gap) are stabilized by prior information obtained from our own low-degree GOCE-only models, so that the final models are still GOCE-only,
• the models are estimated up to relatively high harmonic degrees (75, 100, 130) to mitigate the aliasing, but it is not recommended to use them beyond degree, say, 70 due to the adverse signal-to-noise ratio,
• the models are provided in the gfc format as defined by ICGEM (International Centre for Global Earth Models).

The stabilization of the near-zonal coefficients is shown in terms of difference degree amplitudes in the figure below (Bucha et al., 2015).

When using the models, please use the following reference:

Bucha, B., Bezděk, A., Sebera, J., Janák, J., 2015. Global and regional gravity field determination from GOCE kinematic orbit by means of spherical radial basis functions. Surveys in Geophysics 36, 773-801, http://doi.org/10.1007/s10712-015-9344-0

Go to top

Divergence-free spherical harmonic models of the gravity field implied by the Moon's topographic masses

These models approximate the gravitational field induced by topographic masses of the Moon and mitigate the divergence effect of spherical harmonics in the vicinity of the Moon's topography. The models rely on the Runge—Krarup theorem and enable generally a more accurate evaluation of gravity field quantities close to the lunar topography as compared to the models from spectral gravity forward modelling. This is because the latter ones may suffer from the divergence effect when evaluating the spherical harmonic series on or below the limit sphere encompassing all the gravitating masses (that is, also on the topography). In the figure below, the two approaches are validated with respect to a third (reference) one (Bucha et al., 2019). The upper map shows the reference signal, the middle one depicts the differences related to the Runge—Krarup solution and the bottom one validates the spectral gravity forward modelling technique; note the different colour bars.

We provide in total of 12 models that approximate the gravity field implied by the lunar topographic masses. The topographic masses are expanded up to degrees 90, 180, 360 and 720, and the maximum degree of the gravity field models varies from 360 up to 2160. All models (~290 MB) are available for download here in the gfc format as defined by ICGEM. One of the models, "STU_Moon_topography_to720_gravity_to2160", can also be accessed from ICGEM, where it can be find under a shortened name "STU_MoonTopo720".

The corresponding models that may suffer from the divergence effect were derived by Hirt and Kuhn (2017) and can be downloaded here, including grids of reference gravity disturbances from spatial-domain Newtonian integration.

Please use the following reference when using the models in your works:

Bucha, B., Hirt, C., Kuhn, M., 2019. Divergence-free spherical harmonic gravity field modelling based on the Runge—Krarup theorem: a case study for the Moon. Journal of Geodesy 93, 489-513, https://doi.org/10.1007/s00190-018-1177-4

Go to top

Ultra-high-degree surface spherical harmonic analysis

These routines are designed for accurate and fast FFT-based ultra-high-degree (tens of thousands) surface spherical harmonic analysis. As their key features, the exact Gauss—Legendre quadrature is employed and Legendre's functions are evaluated after Fukushima (2012). Two versions are available:

• non-parallel, which is designed to harmonically analyze an input grid, while storing the grid, output spherical harmonic coefficients, Legendre's functions, etc., in RAM during the entire computation (requires tens of GBs of RAM when the maximum degree is equal to a few tens of thousands),
• parallel, which performs ultra-high-degree surface spherical harmonic analysis of integer powers of input function represented by spherical harmonic coefficients; this code is specifically tailored to HPC clusters, where only a few GBs of RAM per worker may be available.

As an example, shown below are accuracy levels achieved in a closed-loop test environment, where signals up to degrees 90, 180, 360, 720, 1080, 2160, 4320, 6480, 8640 and 10800 were harmonically analysed. These results were obtained by the "SHA_GL.m" routine (non-parallel computation) and the Earth's topography was used as the input signal (Earth2014.TBI2014.degree10800). The first figure shows the reference signal up to degree 10,800 and the second one reveals the accuracy achieved (recovered minus reference coefficients), both in terms of (difference) degree amplitudes (square root of degree variances). The third figure shows computation times required by the "SHA_GL.m" routine to perform the analysis. The tests were conducted on a PC with Intel® CoreTM i7-6800K CPU, 128 GB of RAM (no more than a few GBs were needed in this example though) and a 250 GB SSD drive.

When using the routines, please use the following reference:

Bucha, B., Hirt, C., Kuhn, M. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Submitted to Journal of Geodesy on 14 March 2019.

Go to top

Molodensky's truncation coefficients for cap-modified spectral gravity forward modelling (near- and far-zone gravity effects implied by topographic masses)

Cap-modified spectral gravity forward modelling is a spectral technique to deliver near- and/or far-zone gravity effects implied by topographic masses, the shape of which is given as a surface spherical harmonic series. From the numerical point of view, one of its most difficult parts is the accurate computation of the associated Molodensky's truncation coefficients. This is because this evaluation may necessiate to extend the number of significant digits from 16 (double precision) to, say, 64 or even 256. Here, this was achieved via the ADVANPIX toolbox, which can be used in Matlab.

Some pre-computed ready-to-use values of the near-zone and far-zone truncation coefficients that were computed within the Bucha et al. study (see below) are available here (~1.6 GB). The coefficients are evaluated for

• the spherical distance of psi_0 = 100000 m/6378137 m (100 km integration radius from the evaluation point),
• the reference sphere having the radius R = 6378137 m,
• the radius of the evaluation point r = 6378137 m + 7000 m,
• harmonic degrees n = 0,...,21600,
• topography powers p = 1,...,30,
• radial derivatives k = 0,...,40,
• the first- and second-order horizontal derivatives.

The coefficients were computed using 256 significant digits, ensuring 24-digit accuracy or better. After the evaluation, the coefficients were converted to double precision with 16 significant digits. Importantly, in some cases, the loss of significance errors may be encountered during the spherical harmonic sythesis when using the coefficients (see the reference below).

Please use the following reference when using either the codes or the coefficients:

Bucha, B., Hirt, C., Kuhn, M. Cap integration in spectral gravity forward modelling up to the full gravity tensor. Submitted to Journal of Geodesy on 14 March 2019.

Go to top

List of publications

• Hirt, C., Bucha, B., Yang, M., Kuhn, M., 2019. A numerical study of residual terrain modelling (RTM) techniques and the harmonic correction using ultra-high-degree spectral gravity modelling. Journal of Geodesy, https://doi.org/10.1007/s00190-019-01261-x
• Hirt, C., Yang, M., Kuhn, M., Bucha, B., Kurzmann, A., Pail, R., 2019. SRTM2gravity: An Ultrahigh Resolution Global Model of Gravimetric Terrain Corrections. Geophysical Research Letters 46, 4618-4627, https://doi.org/10.1029/2019GL082521
• Bucha, B., Hirt, C., Kuhn, M., 2019. Divergence-free spherical harmonic gravity field modelling based on the Runge—Krarup theorem: a case study for the Moon. Journal of Geodesy 93, 489-513, https://doi.org/10.1007/s00190-018-1177-4
• Bucha, B., Hirt, C., Kuhn, M., 2019. Cap integration in spectral gravity forward modelling: near- and far-zone gravity effects via Molodensky's truncation coefficients. Journal of Geodesy 93, 65-83, https://doi.org/10.1007/s00190-018-1139-x
• Rexer, M., Hirt, C., Bucha, B., Holmes, S., 2018. Solution to the spectral filter problem of residual terrain modelling (RTM). Journal of Geodesy 92, 675-690, https://doi.org/10.1007/s00190-017-1086-y
• Bucha, B., Janák, J., Papčo, J., Bezděk, A., 2016. High-resolution regional gravity field modelling in a mountainous area from terrestrial gravity data. Geophysical Journal International 207, 949-966, http://doi.org/10.1093/gji/ggw311
• Bucha, B., Bezděk, A., Sebera, J., Janák, J., 2015. Global and regional gravity field determination from GOCE kinematic orbit by means of spherical radial basis functions. Surveys in Geophysics 36, 773-801, http://doi.org/10.1007/s10712-015-9344-0
• Bucha, B., Janák, J., 2014. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders: Efficient computation at irregular surfaces. Computers & Geosciences 66, 219-227, http://doi.org/10.1016/j.cageo.2014.02.005
• Bucha, B., Janák, J., 2013. A MATLAB-based graphical user interface program for computing functionals of the geopotential up to ultra-high degrees and orders. Computers & Geosciences 56, 186-196, http://doi.org/10.1016/j.cageo.2013.03.012

Preprints can be found at my ResearchGate profile.

Go to top

Go to top