2.1. The Utilized Solution
The GGI method was developed as a quasigeoid modeling method based on the geophysical gravity data inversion technique. Generally, it consists of building a local disturbing potential model represented by three components [
17,
18]:
The potential produced by the local topographic masses included in Volume .
The potential produced by the mass anomalies located between the geoid and the Moho surface, which are included in Volume . Volume horizontally covers the same area as Volume .
The component (external disturbing potential), which represents the disturbing potential not covered by the and components.
Volumes
and
are horizontally limited, slightly beyond the study area, and do not cover all the masses producing the disturbing potential (
Figure 1). Hence, the
component covering the influence of the masses located further away is introduced; a long wavelength nature is assumed for this component. Based on this, this component can be represented in the form of low-degree polynomials:
where
are polynomial coefficients,
is the height of the point
P, and
are the coordinates of the point in the local Cartesian coordinate system. The origin of the system is located in the middle of the study area, at the geoid level. The
X and
Y axes lie on the horizontal plane and are directed toward the north and east, respectively. The
Z axis is directed towards the geodetic zenith.
The components
and
are provided by Newton’s integrals:
where
G is Newton’s gravitational constant,
and
are density distribution functions,
and
are elements of volumes, and
l is the distance between the computational point and the attracting masses.
Finally, the disturbing potential on the terrain surface can be defined as follows:
The unknown parameters of the model (Equation (6)) are the 3D density-distribution functions
and
which are defined in Volumes
and
, respectively, and the polynomial coefficients
. The model parameters are determined using the least squares method based on the gravity anomalies or gravity disturbances provided in a dense network of gravity points and the disturbing potential values provided in a sparse network of points with known GNSS/levelling height anomalies, which are converted into disturbing potential values using the Bruns formula [
1].
The density-distribution functions
and
are determined using the procedure for the linear inversion of gravity data [
19] in its discrete form. Therefore, Volumes
and
were divided into finite volume blocks of constant density values. Until now, a very simple division into constant density blocks was used in these calculations; this method assumed only one layer of lateral density variation for Volume
and one layer of lateral density variation for Volume
. In applications, Volume
is defined by digital elevation model (DEM) blocks, which are grouped into zones of constant density (see [
18],
Figure 2). The constant density blocks defining Volume
in the horizontal plane correspond to the constant density zones of Volume
and extend from the geoid to the Moho surface. Both Volumes
and
are represented by rectangular prisms. Equations (4) and (5) can now be written in the following form [
20]:
where
is the determined constant density of Zone
k,
n is the number of zones,
is the number of rectangular prisms of the DEM in Zone
k,
is the determined density of the rectangular prism
j defining the volume
, and
s is the number of prisms.
The coefficients
and
are Newton’s integrals for the rectangular prisms
i of the DEM and
j of the volume
, respectively:
where
and
are the coordinates defining the rectangular prism
i and
j respectively;
;
; and
.
The integrals in Equations (9) and (10), and their derivatives, are then replaced with exact solutions for prism ([
21], Equations (4) and (8)).
In the calculations, the sphericity of the Earth is also considered. The height coordinates of the DEM blocks and the prisms of Volume κ are changed by the correction
, where
R is the mean Earth radius and
is the horizontal distance between the point
and the particular prism centre [
16].
An important step of the calculations is the adoption of the initial density models
for Volume
and
for Volume
. The model
is usually assumed to be constant for all topography. The model
is defined as a kind of topographic-isostatic model, whereby the masses of Volume
balance the masses of Volume
. According to this, in the discrete form of the
model, the constant density value of Block
j of Volume
located exactly below the constant density zone
i of Volume
is defined by the following equation:
where
is the height of Block
j of Volume
κ and
is the mean height of Zone
i of Volume
.
The discrete initial density model defined above can be written in the following form:
The final density model can be written as follows:
Now, we can define the vector of unknown parameters of the model (Equation (6)) in the following form:
where
is the vector of residual densities and
.
Considering the defined vector of unknowns (Equation (14)), the observation equation for the disturbing potential
can be written as follows [
16]:
where
is the vector of known coefficients resulting from the equation of the disturbing potential model (Equation (6)), with its components defined by Equations (3)–(5).
For the gravity disturbance (
), the appropriate equation will have the following form:
where
is the known vector resulting from taking the
derivative of Equation (6).
The values
and
in Equations (15) and (16) are adjustment errors, and the approximate observation quantities are defined as:
and
are determined based on the vector
, where
is a five-dimensional zero vector.
Equations (15) and (16) for the series of observations form a system of equations, which we will write as follows:
where
is the design matrix of known coefficients,
is the vector of observations, and
is the vector of adjustment errors.
To eliminate the ambiguity of the gravity data inversion, the approach suggested by [
22] was used. An additional condition is imposed on the determined parameters:
where
,
is the zero weighting matrix assigned to the vector
, and
is the density model weighting matrix. A detailed description of the definition and determination of the model-weighting matrix can be found in [
18].
With the defined weight-observation matrix
P and considering Equation (20), the least squares objective function can be written as follows:
The condition defined in Equation (21) leads to the solution of the system of equations defined in Equation (19):
The determined parameters of the model (the densities
and
of Volumes
and
respectively, as well as the coefficients
) allow for the calculation of various functionals of the disturbing potential at any point in the considered area based on Equation (6). In order to determine the quasigeoid model, on the basis of the DEM, the heights of the regular grid nodes are determined; the disturbing potential values are calculated from the GGI model. These values are converted using the Bruns formula to the height anomalies (
Figure 1):
where
is the disturbing potential calculated at Point
P on the terrain surface and
is the normal gravity at the telluroid.
Let us also add that calculations can be performed with the use of global geopotential models (GGMs). In this case, the remove–compute–restore procedure is used (e.g., [
16]), and the quasigeoid models determined are slightly more accurate than models built without the use of GGMs ([
23,
24]).
The quasigeoid model determined according to the procedure presented above (with or without the use of GGM) is fitted to GNSS/leveling height anomalies. In [
16], a modification of the GGI solution was proposed that allows for the determination of gravimetric quasigeoid models, which are fitted to the particular GGM being used. In this case, the GGI model is developed without GNSS/levelling data (only GGM and gravity data are used). Both quasigeoid modeling solutions using the GGI method (fitted to GNSS/levelling data and the gravimetric solution) mentioned were used to build the current quasigeoid model for Poland ([
25]).
Thus far, the disturbing potential values have been determined by the GGI approach at points on the terrain surface, allowing for the determination of the height anomalies (Equation (23)). However, the disturbing potential can be determined from the model (Equation (6)) at any point, including inside Volumes
and
or at points on the geoid (
Figure 1). By determining the disturbing potential at Point
on the geoid (
), the geoid undulation can also be determined based on the Bruns formula:
where
is the value of the normal gravity at Point
on the ellipsoid.
Hence, we can determine the GQS value directly as a difference:
where
In Equations (23)–(25), the values and are determined using the same GGI model provided by Equation (6), and the parameters are estimated according to the standard GGI modeling procedure described above.
Note also that Equations (23)–(25) represent the basic form of the GGI model i.e., without a GGM. In this form of the model, the two main components of the GQS values defined in Equation (25) correspond to the main components of the complete GQS values defined by the classical approaches (e.g., [
14] Equation (3)). The
component depends on the density distribution of the masses lying below the geoid, which corresponds to the Bouguer gravity anomaly component in the classical approach, whereas the
component is defined almost the same as the topographic potential correction in the classical approach. Because of these similarities, we will consider and analyse this basic form of the GGI model in the following.
The points with known GNSS/leveling height anomalies necessary for the GQS values modelling can be replaced by a regular grid of height anomalies obtained from an existing quasigeoid model. This approach is recommended because GNSS/levelling points are usually arranged in an irregular grid, often with low resolution. The use of the existing quasigeoid model allows for the optimal selection of these points. This does not present any problem when we note that modeling of the GQS values with the GGI method is an extension of the functionality of this method and can be considered as a complementary step in quasigeoid modeling. In this case, much of the preparatory and computational work completed at the quasigeoid modeling can be used directly in the GQS and geoid determination. Of course, if a quasigeoid model has already been developed (using any approach) in the area being considered, the solution can only be used to implement the GQS value determination.
If the constant density of the topographic masses is assumed, or if a known model of these densities is used and will be not corrected in the GGI modelling process, the following equality will hold true:
. The determined density model will, therefore, only include the densities of Volume
and will take the following form:
Hence, the vector of the determined model parameters (Equation (6)) will be written in the following form:
where
.
In this case, constant density zones of volume are defined only to estimate the initial density model (Equation (11)).