See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/316889552

A software package for computing a regional gravimetric geoid model by the

KTH method

Article in Earth Science Informatics · January 2015

DOI: 10.1007/s12145-014-0149-3

CITATIONS READS

22 681

2 authors:

Ramazan Alpay Abbak Aydın Üstün

Technical University of Konya Kocaeli University
36 PUBLICATIONS 274 CITATIONS 64 PUBLICATIONS 458 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Assesment of bathymetric maps View project

Improvement of geoid and orthometric heights using topographical density and vertical gravity gradient measurements View project

All content following this page was uploaded by Ramazan Alpay Abbak on 16 July 2022.

The user has requested enhancement of the downloaded file.

Earth Sci Inform (2015) 8:255–265
DOI 10.1007/s12145-014-0149-3

SOFTWARE ARTICLE

A software package for computing a regional gravimetric

geoid model by the KTH method
Ramazan Alpay Abbak · Aydin Ustun

Received: 18 September 2013 / Accepted: 30 January 2014 / Published online: 8 March 2014
© Springer-Verlag Berlin Heidelberg 2014

Abstract Nowadays, the geodetic community has aimed to Introduction

determine 1-cm accuracy gravimetric geoid model, which
satisfies the demands of most engineering applications. The geoid, which is a equipotential surface of the Earth’s
However, the gravimetric geoid determination is a diffi- gravity field, coincides with the mean sea level more or less
cult mission which needs an exclusive attention to obtain by extending inside the continents. Due to irregular distribu-
reliable results for this purpose. Today, Least-Squares Mod- tion of topographical masses and their densities, the geoidal
ification of Stokes (LSMS) formula which is so-called the heights vary from −105 m to +85 m with respect to the
KTH method (Swedish Royal Institute of Technology) has best fitting reference ellipsoid. Hence, the knowledge about
been performed in the regional geoid studies. Based upon the position of the geoid and/or its surface normals known
the earlier investigations, the KTH method provides more as vertical directions highly contributes to the studies of
reasonable results than the Remove Compute Restore tech- the Earth’s interior, oceanographical processes, geophysical
nique, especially in roughly terrain with sparse terrestrial interpretations etc. Furthermore, the precise geoid informa-
gravity data. Nevertheless, a compact and practical software tion plays an important role in the geodetic infrastructure
package is now not available for users and researchers in because the topographic elevations and depths of the sea are
geosciences. Thus, in this paper, a scientific software called referred from it. In other words, a variety of implementa-
“LSMSSOFT” is developed and presented by adding a new tions in geodesy, geophysics, oceanography and engineering
algorithm which speeds up the evaluation of Stokes’ inte- requires such kind of heights which are related to the geoid.
gral. Afterwards, the LSMSSOFT is applied to a case study For this purpose, the spirit levelling technique presents an
for the construction of a geoid model over the Auvergne accurate height determination. Nevertheless, this technique
test area in France. Consequently, the algorithm treated in is very time consuming, laborious and costly task for the
the software and its results imply that the LSMSSOFT is an geodesists to reach the reasonable results. During the two
alternative software package for modelling the gravimetric last decades, the Global Positioning System (GPS) provides
geoid by the KTH method. us an alternative technique to obtain the accurate heights
associated with the Earth’s gravity field. However, this alter-
Keywords Evaluation of Stokes’ integral · Gravimetric native method requests a refined gravimetric geoid model
geoid · Modification of Stokes’ formula · Software in order to supply a comparable accuracy with the spirit
package · Terrestrial gravity data levelling, i.e. a few centimetres.
The gravimetric geoid model can be precisely com-
Communicated by: H. A. Babaie puted by Stokes’ integral with the help of terrestrial grav-
ity anomalies distributed globally. For practical imple-
R. A. Abbak (!) · A. Ustun
mentations, however, the integration domain is spatially
Department of Geomatics Engineering, Selcuk University,
Konya, Turkey restricted owing to the limited availability of terrestrial
e-mail: aabbak@selcuk.edu.tr gravity measurements over the entire Earth. Molodensky
A. Ustun et al. (1962) proposed that the truncation bias (due to igno-
e-mail: austun@selcuk.edu.tr rance of gravity data at remote zone) can be reduced if the
256 Earth Sci Inform (2015) 8:255–265

terrestrial gravity data is combined with the long- GPS/levelling data in an absolute and relative sense. Finally,
wavelength information of the gravity field derived from a short summary concludes this paper.
Global Geopotential Model (GGM). This method is called
modification of Stokes’ formula. Advancing in technol-
ogy and geodetic theory, nowadays, this combination has The KTH Method with additive corrections
been very effective method for reaching 1-cm accuracy
gravimetric geoid model. The KTH method based on the Least-Squares Modifica-
The modifications of Stokes’ formula can be broadly tion of Stokes’ (LSMS) Formula has been reviewed in
split into two categories; deterministic and stochastic several geodetic literature (e.g. Ellmann 2004; Kiamehr
approaches. While the deterministic approaches aim to min- 2006; Ulotu 2009; Abbak 2012b). In this section, we briefly
imize the truncation error caused by the lack of terrestrial explain the computation of geoidal height by the KTH
gravity data in remote zone by assuming that GGM coeffi- method and its additive corrections.
cients are errorless, stochastic methods are used to reduce According to the KTH approach, the geoidal height
global mean square error of the truncation error as well as N̂ can be summarized by the following formula (Sjöberg
random errors of GGM and terrestrial gravity data (Sjöberg 2003a):
2003a). Some other researchers such as Wong and Gore
(1969), Meissl (1971), Vincent and Marsch (1974), and N̂ = Ñ + δNcomb
Top
+ δNDWC + δNcomb
Atm
+ δNell (1)
Vanı́ček and Kleusberg (1987) applied different determin-
istic methods. Then, Sjöberg (1981, 1984, 1991, 2003a) where Ñ is the fundamental term of (approximate) geoid
Top
suggested a stochastic method (the KTH method) for the undulation, δNcomb is the combined topographic correction
modification. His proposals consider not only the trunca- including the total effect of the direct and indirect topo-
tion error but also errors of the GGM and terrestrial data, graphical effects of gravity anomaly to the sea level from
and GGM and local gravity dataset are combined in the the earth surface point, δNDWC is the downward continuation
spectral domain (i.e. frequency domain) to determine a pre- (DWC) of gravity anomaly, δNcomb Atm
is the combined atmo-
cise regional geoid model. Sjöberg and Hunegnaw (2000) spheric correction containing (direct and indirect) atmo-
compared deterministic and stochastic methods in the geoid spherical effects on the geoid and δNell is the ellipsoidal
modelling, and they demonstrated that the reasonable way correction for the spherical approximation of the geoid in
to obtain a more accurate model is to combine GGMs and Stokes’ formula to the reference ellipsoidal surface.
local data via the stochastic methods. A numerical study The KTH method combines terrestrial gravity data and
comparing some deterministic and stochastic modification long-wavelength components of the gravity field in a
methods can also be found in Ellmann (2005b). Also, some stochastically modified Stokes’ kernel. Accordingly, the
comparisons show that the KTH method is the most pre- stochastic KTH method uses the least squares principle
cise method in the gravimetric geoid modelling (Ågren et al. to minimize the expected global mean square error of the
2009; Abbak et al. 2012a; Yildiz et al. 2012). modified Stokes’ formula (Sjöberg 1984, 1991, 2003a).
Although the KTH is known as a successful method in By taking advantage of the orthogonality of spherical har-
the gravimetric geoid determination, still a compact and monics over the sphere, the original stokes integral can be
practical software package is not available for users and rewritten by two terms,
researchers studying in geodesy and related disciplines, e.g. !! M
geophysics. Therefore, in this paper, a scientific software R R "
Ñ = S L (ψ)&gdσ0 + bn &gn (2)
package which is called “LSMSSOFT” is developed in 4π γ 2γ
σ0 n=2
the C/C++ programming language by adding a new algo-
rithm and a range of options. By this way, the LSMSSOFT where bn is the set of spectral parameters that enables to use
is faster, simpler and more user-friendly than the KTH- the observed (original) gravity anomalies instead of residual
GEOLAB developed in the FORTRAN programming lan- anomalies in Eq. 2. R is the mean Earth’s radius, γ is the
guage by Kiamehr and Sjöberg (2010). normal gravity at computation point, S L (ψ) is the modified
The current paper starts with a brief review of the com- Stokes’ function, L is maximum degree of the modification,
putational scheme of the KTH method with additive cor- &gn is the gravity anomaly computed from the GGM, M is
rections. Then, we introduce the LSMSSOFT that uses an the upper limit of the spherical harmonic expension.
algorithm of the KTH method and its corrections. After- The combined topographic correction term in Eq. 1 can
wards, we explain the input data and procedure of grid- be computed by Sjöberg (2007),
ding terrestrial gravity anomalies over the target area. We # $
apply the LSMSSOFT for the construction of a gravimet- 2π GρH 2 2H
δNcomb = δNdir + δNind = −
Top Top
1+ (3)
ric geoid model of the Auvergne region by comparing with γ 3R
Earth Sci Inform (2015) 8:255–265 257

where G is the Newtonian gravitational constant, ρ is the The LSMSSOFT for computing a gravimetric geoid
Earth’s crust density, and H is the elevation of the topogra- model
phy at the computation point, P .
The DWC correction term in Eq. 1 can be considered by This section comprehensively explains the general structure
Sjöberg (2003b), of the software encoded in the GNU C/C++ platform. It is
(1) L1,Far L2 named “The Least-Squares Modification of Stokes’ formula
δNDWC = δNdwc + δNdwc + δNdwc (4) SOFTware (LSMSSOFT)”. Subsequently, with the aim of
where speeding up the integration procedures, a new approach
designed by the authors is rigorously treated in this soft-
&gP Ñ 1 ∂&g ware. The next subsection introduces the new approach
(1)
δNdwc = HP + 3 HP − |P HP2 (5)
γ rp 2γ ∂r briefly.
and
%# $ & Functions in the software
M
L1,Far R " R n+2
δNdwc = bn − 1 &gn (6) The software contains one main function and 6 sub-
2γ rP
n=2
functions (sub-routines) which are CHLS, SVD, COEF-
also FICIENT, LEVELELL, LEGENDRE and HELP. Short
!! ' ( descriptions of all functions are given as follows:
L2 R L ∂&g
δNdwc = S (ψ) |Q (HP − HQ ) dσo (7)
4π γ ∂r • MAIN takes the names of data files and user’s options
σo
from the command line. In this context, there must be
where rP = R + HP is the spherical radius of the point P , totally three data files: GGM, anomaly and elevation.
HP is the orthometric height of point P , Q represents the The GGM file is publicly available from the Interna-
moving integration point. tional Centre for Global Earth Models (ICGEM) web-
The combined atmospheric correction in Eq. 1 can be site (http://icgem.gfz-potsdam.de/ICGEM). The soft-
approximated by Sjöberg (1999), ware has the ability to read all available GGM files in
!! the ICGEM website directly as an input file without
GRρ a
δNcomb = −
Atm
S L (ψ)HP dσo (8) any change. The anomaly file in which data is stored
γ
σo in the point format (ϕ, λ, &g), must cover the target
area and at least 2 degrees the nearest vicinity of the
where ρ a is the density of the atmosphere at sea level. Gen- target area. The elevation file must coincide with the
erally, it is selected that ρ a is equal to 1.23 kg/m3 for the anomaly file in respect to the resolution and coverage
computations (Sjöberg 1999; Ecker and Mittermayer 1969). as well. Some parameters (M, ψ, C0) and boundary of
The ellipsoidal correction in Eq. 1 can be simply calcu- the target area with the resolution should be entered to
lated by Ellmann and Sjöberg (2004), program. After all, the computations according to the
theory given in the previous section, are completed by
δNell ≈ [(0.0036 − 0.0109 sin2 ϕ)&g + 0.0050Ñ cos2 ϕ]QL
0 considering all parameters and options, then, the grid-
(9) ded geoid heights (and all segments if desired) in the
where QL target area are printed to the stdout sequentially. The
0 denotes the Molodensky truncation coefficient.
Regarding the computational efficiency, the formulae flowchart that shows the input files, their formats and
above can be rearranged as follows, outputs are illustrated in Fig. 1.
!! • CHLS generates the inversion of the design matrix
N̂ =
R
S L (ψ)[&g + &g ∗ − 4π Gρ a HP ]dσ0 A of unknown modification coefficients by using the
4π γ method of Cholesky Decomposition. Cholesky decom-
σ0
position for the least-square solution estimates the mod-
# $n+2
R "
M
R (10) ification parameter vector s from the linear system of
+ bn &gn
2γ rP equations,
n=2
+ δNcomb
Top (1)
+ δNdwc + δNell As = h (11)

where &g ∗ = ∂&g ∂r |Q (HP − HQ ). In the next section, where A and h contain the coefficients of the lin-
considering the compact (10) the computational scheme of ear sytem and constant terms, respectively (Ellmann
the KTH approach and its components as subroutines are 2005a). Because the algorithm of Cholesky Decompo-
introduced step by step. sition is more suitable for a symmetric and positive-
258 Earth Sci Inform (2015) 8:255–265

Fig. 1 The flowchart of the

LSMSSOFT

definite matrix, this algorithm is employed for the using the recursion formulae. The parameters of this
“biased” version of the KTH method. function are the spherical latitude and maximum degree
• SVD performs the singular value decomposition of the of the polynomial.
design matrix A. Since the design matrix becomes ill- • HELP supplies the general information about the usage
conditioned in the “unbiased” and “optimum” version of the software. In the case where “-H” or unexpected
of the KTH method, the solution for s becomes numer- options are used, the HELP function becomes active.
ically unstable. Thus, there is an alternative method It shows how the LSMSSOFT works together with
which is the Singular Value Decomposition (SVD) that the options, and what parameter is needed for each
avoids the unstability of the inversion of a symmet- option.
ric matrix. The algorithm in this software was adopted
from Press et al. (2007). The software should be compiled in the Linux based
• COEFFICIENT estimates the modification parameters server or personal computer, or in Windows operating sys-
sn and bn by the least-squares method. The function tem by installing free like-linux platforms, e.g. Cygwin.
considers three stochastic modifications, i.e. the biased Cygwin can be downloaded from http://www.cygwin.com.
(Sjöberg 1984), unbiased (Sjöberg 1991) and optimum
(Sjöberg 2003a) version of the KTH method. In this The new algorithm
function, the whole computations are realized according
to Ellmann (2005a), and results converges his studies. The main topic of this subsection is the efficient evaluation
• LEVELELL generates the normal gravity field by con- of the convolution integral on the sphere, i.e. the Stoke-
cerning the reference ellipsoid. As a default, the nor- sian integration, by means of the new approach. With the
mal geopotential coefficients are calculated based upon intention of expediting the computational procedures of the
GRS80 reference ellipsoid. However, any geopotential integration in the Stokes’ formula, a new algorithm designed
model with the different GM value can be used in the by the authors was adopted to the LSMSSOFT properly.
software as it rescales the coefficients only. A similar philosophy for the algorithm was succesfully
• LEGENDRE calculates the fully normalized associated adapted by Huang et al. (2000). Our algorithm solely con-
legendre polynomials of the first and second kind by siders a compartment for looking a running point Q, which
Earth Sci Inform (2015) 8:255–265 259

is varying according to the grid and integration cap sizes. The performance of the new strategy was tested against
This structure will make it faster than earlier approach. to the classical pointwise method for a gravity data set in
As it is in other gravimetric geoid methodologies, the central France (i.e current study area). For the area with a
Stokesian integration is the most time consuming part in maximum latitude difference of 2 degrees, the new method
the estimation of the geoid model in the KTH scheme, as can be up to 0.17 millimetres bias compared to the clas-
well. This is because the modification part of Stokes’ for- sical pointwise integration, which should be comfortably
mula is repeatedly calculated in every compartment using neglected in the practice. The mean and standard deviation
recursive algorithm (see Fig. 2). In order to avoid these of the differences are 0.00 mm and 0.04 mm, recpectively.
repeated computations, at the beginning of the longitude, Related with the time difference, the new method finishes
a sample computation is done by encoding the values in the integration in 0.53 seconds on a computer with 2.0 GHz
respect to the positions of the running point. Then, in every Intel processor, whereas the pointwise method is up to 29.00
compartment, the values of the modification parts are sim- seconds on the same processor. This means that the classi-
ply recalled from these encoded computations. Additionally, cal method is 55 times worse than the new method in point
the half of the compartment is enough to be evaluated since of view of CPU time.
it has longitudinal symmetry.
The new algorithm can be briefly defined as follows: The
curser starts the integration from the minimum latitude and A numerical example in the Central France
longitude introduced by the user. The research area (here-
after compartment) is rectangular shaped according to the In this section, we test the potential of the LSMSSOFT by
cap and grid sizes (Fig. 2). Considering the cap size ψ0 using an example data from Auvergne test area. At first, we
(specified by the user), the curser begins to look for the make an overview of the datasets, which will be utilized in
running point Q. If the radial distance between the compu- the LSMSSOFT for geoid modelling procedures. Then, we
tation point P and running point Q is less than or equal to explain our scheme of gridding the gravity anomalies in the
the cap size ψ0 , the modification part of Stokes’ function target area. These gravity anomalies are used in the con-
S M (ψ) is calculated, otherwise it is passed. This calculated struction of a gravimetric geoid model by applying the KTH
modification value is encoded with respect to its position approach via the LSMSSOFT. Finally, we evaluate this new
with x and y coordinate pairs in the compartment (i.e x model by means of the GPS/levelling data in an absolute
and y are the numbers of grid from the lower left corner of and relative sense.
the compartment). All integrations in the east-west direction
are made with the help of previously computed Stokesian Input data
modification parts. When we move to the next latitude, the
modification parts are recomputed at the minimum longi- The target area is the Auvergne area located in the central
tude. Similarly, all integrations can be repeated for every part of France. Its total area is approximately 40 000 km2
parallel where the geoidal heights are computed. The new including inner waters. The lowest and the highest points in
algorithm is portrayed in Fig. 2. the target area are 0 m at the southeastern part of the area
and 2800 m at the peak of Alps. The area has the average
height about 500 m. Its geographical limits are from 45◦ to
47◦ northern latitude and from 2◦ to 4◦ eastern longitude
(Fig. 3). This region was selected because it is one of the
most complicated areas in France from the point of view of
rough topography and geoidal height variation.
The gravity data used in the current research are obtained
from Duquenne (2007). The gravity points are related to the
International Gravity Standardization Net 1971 (IGSN71),
and also its geographical datum is the World Geodetic Sys-
tem 1984 (WGS84). The geographical distribution of all
gravity data points in the target area is shown in Fig. 3b.
The gravity data coverage is satisfactory. The number
of gravity points within the target area is about 190 000,
which produce a density of one point per 2 km2 , approx-
imately. As it is well known, the density of the gravity
Fig. 2 Evaluating the integration within compartment centered at the data directly affects the accuracy of any geoid model.
computation point P Also the procedure of the LSMSSOFT needs the gravity
260 Earth Sci Inform (2015) 8:255–265

Fig. 3 Topography and

distribution of the gravity
surveys (black points) and
GPS/levelling benchmarks (red
triangles)

data outside of the edges of study area. In this sense, the formula. Satellite-only models may avoid this correlation
gravity data in the nearest vicinity of the target area is when they are used as a reference GGM.
available. The EGM Development Team of the National
In 2000, the National Aeronautics and Space Adminis- Geospatial-Intelligence Agency (NGA) in USA, has been
tration (NASA) organized the Shuttle Radar Topography released the Earth Gravitational Model 2008 (EGM2008)
Mission (SRTM) project to generate the most complete on the internet freely (Pavlis et al. 2012). EGM2008
high-resolution digital topographic database of the Earth to presents remarkable improvement compared to previous
date. Thus, a Digital Elevation Model (DEM) with 3×3 arc- combined GGMs due to the release of new gravity data
second resolution was created and made publicly available from more sources, such as terrestrial gravity, GRACE
on the Internet. The DEM will fulfill the elevation needs in tracking and satellite altimetry data. Therefore, EGM2008
the current study for the topographical, atmospherical and is the best combined model all over the world, comparing
downward continuation corrections. And also the DEM will via external data (for details see all papers in Newton’s Bul-
be utilized to predict the mean gravity anomalies in our letin Nr. 4 entitled “External Quality Evaluation Reports of
gridding procedure. EGM2008”). Thus, in the current study EGM2008 will be
GGMs are the representation of Earth’s gravitational used for the filling the gravity gaps in both the water bodies
potential in terms of spherical harmonic coefficients and mountainous regions where terrestrial gravity data is
obtained from space-borne and terrestrial data sources. unsurveyed.
Every model is different from the other with respect to the In the target area, there are 56 GPS/levelling bench-
resolution, data source, analysis method, quality and density marks which validate the new geoid model in an absolute
of data. Thus the selection of any GGM in geoid deter- and relative sense. The accuracy of the ellipsoidal heights
mination directly affects the regional geoid solution (e.g. determined by GPS is estimated to be ±3 cm, whereas the
Abbak 2012b). accuracy of orthometric heights can be estimated to be ±2
ITG-GRACE2010S will be used for the regional geoid cm in the study area (Duquenne 2007).
modelling with the LSMSSOFT in this study. This model is
a satellite-only model derived from 84 months of GRACE Gridding the gravity anomalies
tracking data, which is complete to degree and order 180
(Mayer-Guerr et al. 2010). We choose this satellite-only The observed gravity data on physical surface of Earth can
model because this is an appropriate satellite-only GGM for be generally represented as Free-air or Bouguer gravity
long-wavelength components of gravity knowledge, com- anomalies in geodesy and geophysics. The free-air anoma-
paring with the spectral tools (Ustun and Abbak 2010). lies are used for geoid modelling, whereas the Bouguer
Additionally, the GGM and ground gravity data are corre- anomalies are better fitted for interpolation due to their
lated especially when the combined GGM is used in Stokes’ smoothness.
Earth Sci Inform (2015) 8:255–265 261

There are two types of Bouguer gravity anomalies: sim- The KTH geoid model
ple and refined/complete. Janak and Vanı́ček (2005) suggest
that the complete Bouguer anomaly is inevitable in high This subsection presents the computation processes of
mountains with the sparse gravity points while in the low- the new geoid model which is named as Auvergne Geoid
elevations the easier simple Bouguer approach could be Model via the KTH method 2013 (AG13). The LSMSSOFT
sufficient in the interpolation. Since our target area is par- was used for the whole numerical evaluations. Firstly, the
tially mountainous and is of overdetermined gravity data, in Stokesian integration with gridded free-air gravity anoma-
this research the simple Bouguer anomalies were regarded lies was implied in the target area to obtain approximate
enough to predict the mean gravity anomalies. Nonetheless, geoid undulations. Subsequently, other corrections were
the complete Bouguer anomalies should be preferably used separately added to these approximate undulations to get
in rough terrain for the gravity anomaly prediction. the accurate ones.
The gridding gravity anomalies were made as follows: Due to the shortage of terrestrial gravity data out of the
the gravity observations distributed randomly were directly target area, the integration cap size (ψ0 ) around the com-
reduced to the simple Bouguer gravity anomalies. Then putation point is limited to only few degrees while the
Bouguer gravity anomalies were interpolated to grid centers modification of Stokes’ formula is mainly trying to reduce
according to Bjerhammer’s deterministic rule (Bjerhammar the truncation error (Kiamehr 2006). Also the selection of
1973). During the nearest neighbouring interpolation pro- upper bound of Stokes’ function L plays an essential role
cess, the minimum and maximum number of the points for in geoid modelling procedure in terms of the computa-
every grid point were chosen 4 and 6, respectively, and also tional efficiency. Whereas higher degree GGMs (M) may
the maximum interpolation radius was selected as 12 arc- significantly balance the lack of gravity data, the errors
minutes. Finally, simple Bouguer anomalies were converted of the potential coefficients of GGM increase gradually.
to free-air gravity anomalies in grids (3×3 arc-minute reso- Therefore, a reasonable compensation should be determined
lution) by restoring the mean Bouguer plate effects (Fig. 4). among these parameters. For example, this procedure can be

Fig. 4 Free-air gravity anomaly

grid
262 Earth Sci Inform (2015) 8:255–265

done such that the value of a parameter is changing whilst from a simple linear regression to a 7-parameter similarity
the others are keeping constant. The value of the parameter transformation (Kotsakis and Sideris 1999). In this research,
which gives the best result compared to the GPS/levelling the 7-parameter model was preffered since it gives the min-
data is fixed up (e.g. Abbak 2012b). imum standard deviation for gravimetric geoid model of the
Checking the parameters with the GPS/levelling data, target area. While the one parameter model is utilized for the
the final geoid solution of AG13 was computed by the comparison, the standard deviation of AG13 and EGM08
biased version of the KTH method while considering the models are 99.4 mm and 47.2 mm, respectively.
parameters based on the free-air surface gravity anoma- Another way for assesment of the real potential of the
lies within 3×3 arc-minute resolution, ITG-GRACE2010S, geoid models is the test of their fit to the GPS/levelling data
the cap size ψ0 = 1.0◦ , M = L = 120, σ&g 2 = 25 in a relative sense. In this sense, the discrepancies between
2
mGal (Fig. 5). Additionally, the additive corrections over two geoidal height differences &Nij , which are derived
the territory of the Auvergne were calculated by using same from geometric and gravimetric methods are simply cal-
parameters (Fig. 6). culated for the whole baselines. These differences can be
represented in the relative form in parts per million (ppm):
Validation of the AG13 model ) Gra )
) N − N Gra − (N Geo − N Geo ) )
) i j i j )
&Nij = ) ) (12)
There are two ways to evaluate the accuracy of any gravi- ) Dij )
metric geoid model via the GPS/leveling data: the absolute
and relative techniques. The absolute technique is estima- where Dij is the distance between i th and j th points (unit:
tion of the RMS errors of the discrepancies between the km), N Gra and N Geo are the gravimetric and geometric undu-
gravimetric and geometric geoid models. A plenty of math- lations of these points (unit: mm).
ematical models can be used to absorb the long-wavelengths In Table 1, global and regional gravimetric geoid heights
and systematic errors in the geoid model which are ranging are evaluated by using the GPS/levelling data in an

Fig. 5 Geoid heights of the new

gravimetric model (AG13).
Contour interval is 0.25 m
Earth Sci Inform (2015) 8:255–265 263

Fig. 6 All corrections of the

AG13 model in the target area
[m]

absolute and relative sense. Considering both comparisons Conclusions

in the table, the AG13 model has the better statistics than
to EGM08 in the tests. These results also confirm the ear- The present paper summarized the theory of the KTH
lier investigation which was carried out via KTH-GEOLAB method with the additive corrections to a precise regional
software in the same region by Yildiz et al. (2012). Thus, it geoid determination. Then, a scientific software package
is concluded that the LSMSSOFT is fairly successful in the that is constructed according to the aforementioned the-
gravimetric geoid modelling by the KTH method. ory of the KTH method was presented by adding a new

Table 1 Evaluating
gravimetric geoid models by Geoid model Absolute statistics (mm) Relative statistics (ppm)
means of GPS/levelling data in
the absolute and relative senses Min Max Mean RMS Min Max Mean

AG13 −64.4 47.8 0.5 25.6 0.0 3.1 0.4

EGM08 −64.5 88.5 3.6 34.3 0.0 5.7 0.5
264 Earth Sci Inform (2015) 8:255–265

algorithm with a wide number of the options. Finally, the Ellmann A (2004) The geoid for the baltic countries determined by the
paper was followed by a sample numerical application in least squares modification of Stokes formula. PhD thesis. Royal
Institute of Technology (KTH). Division of Geodesy, Stockholm
partially mountainous part of France where the terrestrial
Ellmann A (2005a) Computation of three stochastic modifications
data is overdetermined. The numerical analyses show that Stokes formula for regional geoid determination. Comput Geosci
the software yields us more reasonable results in the target 31(/6):742–755
area, comparing with GPS/leveling data. Ellmann A (2005b) Two deterministic and three stochastic modifica-
tions of Stokes’ formula: a case study for the Baltic countries. J
This study contributes to the important computational
Geodesy 79:11–23
aspects of the KTH method with the additive corrections. In Ellmann A, Sjöberg LE (2004) Ellipsoidal correction for the modified
order to find the best solution in the computation of a gravi- Stokes formula. Boll Geod Sci Aff 63:3
metric geoid model by the LSMSSOFT, the interested users Huang J, Vanı́ček P, Novak P (2000) An alternative algorithm to
FFT for the numerical evaluation of Stokes’s integral. Studia
are strongly advised to experiment with their own data and
Geophysica et Geodaetica 44:374–380
some alternative modification parameters, and then compare Janak J, Vanı́ček P (2005) Mean free-air gravity anomalies in the
them, finally draw conclusions from the test results. mountains. Studia Geophysica et Geodaetica 49:31–42
On the other hand, the software is intentionally designed Kiamehr R (2006) Precise gravimetric geoid model for iran based on
GRACE and SRTM data and the least squares modification of the
to be composed of a range of various sub-functions. This
stokes formula with some geodynamic interpretations. PhD thesis,
structure lets users to replace their own sub-functions (e.g. royal institute of technology (KTH). Department of transport and
an alternative singular value decomposition) without any economics, Stockholm
negative influence on the remained part of the software. Kiamehr R, Sjöberg EL (2010) KTH-GEOLAB scientific software
Ultimately, the software is completely structural and user for precise geoid determination based on the least-squares modi-
fication of Stokes’ formula. Royal institute of technology (KTH),
friendly, and can be simply applied for the academic pur- Division of geodesy, Stockholm. Technical manual
poses. We hope that the availability of the LSMSSOFT Kotsakis C, Sideris MG (1999) On the adjustment of combined
increases the eagerness of the geoscience researchers to con- GPS/levelling/geoid networks 73(8):412–421
duct the KTH method for their own data in the different Mayer-Guerr T, Kurtenbach E, Eicker A (2010) ITG-Grace2010
gravity field model. http://www.igg.uni-bonn.de/apmg/index.php?
target areas. id=itg-grace2010
Meissl P (1971) A study of covariance functions from discrete mean
Acknowledgments This study was financially supported by The value. Technical report, Dept. Geod. Sci. Ohio State University,
Research Fund of Selcuk University of Turkey under grant 09-101- Columbus
009. The first author developed some parts of the software while he Molodensky MS, Eremeev VF, Yurkina MI (1962) Methods for study
was at the Royal Institute of Technology in Sweden under supervision of the external gravity field and figure of the Earth. Israel program
of Prof. Sjöberg and at the Tallinn Technology University in Estonia for scientific translation
under supervision of Prof. Ellmann. The authors cordially appreciate Pavlis NK, Holmes SA, Kenyon SC, Factor JK (2012) The devel-
the constructive remarks made by two anonymous reviewers in the first opment and evaluation of the earth gravitational model 2008
version of this paper. The authors would like to give a special thanks (EGM2008). J Geophys Res Solid Earth 117:38
to Prof. Ramin Kiamehr for his valuable comments on the original Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2007) Numer-
manuscript. ical recipes: the art of scientific computing. Cambridge University
Press, New York
Sjöberg LE (1981) Least squares combination of satellite and terres-
trial data in physical geodesy. Ann Geophys 37(1):25–30
Sjöberg LE (1984) Least squares modification of stokes and venning-
References meinesz formulas by accounting for errors of truncation, potential
coeffcients and gravity data. Technical report, Department of
Abbak RA, Erol B, Ustun A (2012a) Comparison of the KTH and Geodesy, Institute of Geophysics, University of Uppsala
remove-compute-restore techniques to geoid modelling in a moun- Sjöberg LE (1991) Refined least squares modification of Stokes
tainous area. Comput Geosci 48:31–40 formula. Manuscripta Geodaetica 16:367–375
Abbak RA, Sjoberg LE, Ellmann A, Ustun A (2012b) A precise gravi- Sjöberg LE (1999) The iag approach to the atmospheric geoid correc-
metric geoid model in mountainous areas with scarce gravity data tion in Stokes’ formula and a new strategy. J Geodesy 73:362–
– a case study in Central Turkey. Studia Geophysica et Geodaetica 366
56:909–927 Sjöberg LE (2003a) A general model of modifying Stokes formula and
Ågren J, Sjöberg EL, Kiamehr R (2009) The new gravimetric quasi- its least squares solution. J Geodesy 77:790–804
geoid model KTH08 over Sweden. J Appl Geodesy 3:143–153 Sjöberg LE (2003b) A solution to the downward continuation effect on
Bjerhammar A (1973) Theory of errors and generalized matrix the geoid determination by Stokes formula. J Geodesy 77:94–100
inverses. Elsevier, Amsterdam Sjöberg LE (2007) The topographic bias by analytical continuation in
Duquenne H (2007) A data set to test geoid computation methods. physical geodesy. J Geodesy 87:345–350
Harita Dergisi. Special Issue. pp 61–65 Sjöberg LE, Hunegnaw A (2000) Some modifications of Stokes’ for-
Ecker E, Mittermayer E (1969) Gravity corrections for the influence mula that account for truncation and potential coefficient errors. J
of the atmosphere. Boll Geofis Teor Appl 11:70–80 Geodesy 74:232–238
Earth Sci Inform (2015) 8:255–265 265

Ulotu PE (2009) Geoid model of tanzania from sparse and varying Vincent S, Marsch J (1974) Gravimetric global geoid. In: on the use of
gravity data density by the KTH method. PhD thesis, royal insti- artifical satellites for geodesy and geodynamics
tute of technology (KTH), Division of Transport and Economics, Wong L, Gore R (1969) Accuracy of geoid heights from the modified
Stockholm Stokes kernels. J Royal Astron Soc 18:81–91
Ustun A, Abbak RA (2010) On global and regional spectral evaluation Yildiz H, Forsberg R, Ågren J, Tscherning CC, Sjöberg LE (2012)
of global geopotential models. J Geophys Eng 7(4):369–379 Comparison of remove-compute-restore and least squares modifi-
Vanı́ček P, Kleusberg A (1987) The Canadian geoid – stokes approach. cation of Stokes’ formula techniques to quasi-geoid determination
Manuscripta Geodatica 12:86–98 over the auvergne test area. J Geodetic Sci 2(1):53–64
View publication stats