*The Gravity Effect of Topography: A Comparison among Three Different Methods DOI: http://dx.doi.org/10.5772/intechopen.97718*

and polar axis coinciding with the line connecting *O* to the observation point *P* (ECPTesseroid), for which it is possible to use the exact solution in Eq. (9) (**Figure 2**).

Mapping can be done by means of two procedures, depending on the spherical distance of the tesseroid from the observation point (**Figure 3a**).

The first procedure (ST procedure) involves a local second-order decomposition of each tesseroid into a number *NS* of small equal-area sectors of a spherical band, which develop along the local ECP meridians (dashed black lines) and parallels (dashed blue lines in **Figure 3b**). Sensitivity tests [11] show that the *optimum* value of *NS* for which the ST procedure converges varies only with the latitude and decreases from the equator to the North Pole. In the present study we use *NS* ¼ 4, in agreement with [11].

#### **Figure 2.**

one can compute the gravitational effect of a right parallelepiped assuming the computational point *P* at the origin of a Cartesian reference system as the closed

� � � �

This formula is implemented in e.g. the TC program of the GRAVSOFT package

The gravitational effect of spherical tesseroid has been studied in several papers. Relevant studies on this topic have been carried out by [3, 9–12]. In this study the

The UNIPOL approach grounds on the result that a closed formula of the volume Newtonian integral (Eq. (8)) exists when the observation point *P* is located along the polar axis (see e.g. [13]) and the tesseroids coincide with sectors of a spherical

In this case the gravitational contribution of spherical tesseroid of height *h* can

*<sup>φ</sup>*<sup>2</sup> � <sup>2</sup> � � <sup>þ</sup> *r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>2</sup>

*<sup>φ</sup>*<sup>1</sup> � <sup>2</sup> � � <sup>þ</sup> *r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>1</sup>

*zr*

� �,ð Þ *<sup>z</sup>*1, *<sup>z</sup>*<sup>2</sup> are the coordinates of the

*<sup>r</sup>*<sup>2</sup> <sup>Δ</sup>*λ*ð Þ *<sup>I</sup>*<sup>1</sup> <sup>þ</sup> *<sup>I</sup>*<sup>2</sup> <sup>þ</sup> *<sup>I</sup>*<sup>3</sup> <sup>þ</sup> *<sup>I</sup>*<sup>4</sup> <sup>þ</sup> *<sup>I</sup>*<sup>5</sup> (9)

q

q

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *<sup>R</sup>*<sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>2</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *<sup>R</sup>*<sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>1</sup>

1

CA

þ ð Þ� *R* þ *h r* cos *φ*<sup>2</sup>

þ *R* � *r* cos *φ*<sup>2</sup>

� � � *x*1 *x*2 *y*1 *y*2 � � � � � � *z*1 *z*2

(8)

*<sup>∂</sup><sup>z</sup>* <sup>¼</sup> *<sup>G</sup><sup>ρ</sup>* <sup>k</sup>*<sup>x</sup>* ln ð Þþ *<sup>y</sup>* <sup>þ</sup> *<sup>r</sup> <sup>y</sup>* ln ð Þ� *<sup>x</sup>* <sup>þ</sup> *<sup>r</sup> <sup>z</sup>* arctan *xy*

[8] that will be used in the computations described in the next section.

approach presented in [11], hereinafter referred to as UNIPOL, is used.

*<sup>x</sup>*<sup>2</sup> <sup>þ</sup> *<sup>y</sup>*<sup>2</sup> <sup>þ</sup> *<sup>z</sup>*<sup>2</sup> <sup>p</sup> and ð Þ *<sup>x</sup>*1, *<sup>x</sup>*<sup>2</sup> , *<sup>y</sup>*1, *<sup>y</sup>*<sup>2</sup>

*Geodetic Sciences - Theory, Applications and Recent Developments*

*∂V <sup>∂</sup><sup>r</sup>* <sup>¼</sup> *<sup>G</sup><sup>ρ</sup>*

<sup>2</sup> 3 cos <sup>2</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *<sup>R</sup>* <sup>þ</sup> *<sup>h</sup>* <sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>2</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *<sup>R</sup>* <sup>þ</sup> *<sup>h</sup>* <sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>1</sup>

<sup>2</sup> 3 cos <sup>2</sup>

<sup>2</sup> 3 cos <sup>2</sup>

q

*φ*2*r* 3

<sup>3</sup> ð Þ *<sup>R</sup>* <sup>þ</sup> *<sup>h</sup>* <sup>2</sup> <sup>þ</sup> *<sup>r</sup>*

h i

*<sup>φ</sup>*<sup>2</sup> � <sup>2</sup> � � <sup>þ</sup> *rR* cos *<sup>φ</sup>*<sup>2</sup> � � �

*<sup>φ</sup>*<sup>1</sup> � <sup>2</sup> � � <sup>þ</sup> *rR* cos *<sup>φ</sup>*<sup>1</sup> � � �

> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð Þ *<sup>R</sup>* <sup>þ</sup> *<sup>h</sup>* <sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>2</sup>

> > ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi *<sup>R</sup>*<sup>2</sup> <sup>þ</sup> *<sup>r</sup>*<sup>2</sup> � <sup>2</sup>*r R*ð Þ <sup>þ</sup> *<sup>h</sup>* cos *<sup>φ</sup>*<sup>2</sup>

In Eq. (9), *R* is the radius of the Earth and *φ* and *λ* colatitude and longitude, respectively. When the observation point *P* is not located along the polar axis, the UNIPOL approach maps each tesseroid defined in the Earth-Centered Rotational reference frame (ECRTesseroid) into a sector of a spherical zonal band of a spherical cap in the Earth-Centered *P*-Rotational reference frame, having the same origin *O*

h i

<sup>2</sup> 3 cos <sup>2</sup>

formula [2]:

*∂V*

parallelepiped edges.

where *<sup>r</sup>* <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

zonal band of a spherical cap.

<sup>3</sup> ð Þ *<sup>R</sup>* <sup>þ</sup> *<sup>h</sup>* <sup>2</sup> <sup>þ</sup> *<sup>r</sup>*

*<sup>R</sup>*<sup>2</sup> <sup>þ</sup> *<sup>r</sup>*

*<sup>R</sup>*<sup>2</sup> <sup>þ</sup> *<sup>r</sup>*

q

*<sup>I</sup>*<sup>5</sup> ¼ � cos *<sup>φ</sup>*<sup>2</sup> sin <sup>2</sup>

0

B@

be expressed as:

with

*<sup>I</sup>*<sup>1</sup> <sup>¼</sup> <sup>1</sup>

�

*<sup>I</sup>*<sup>2</sup> <sup>¼</sup> <sup>1</sup> 3

*<sup>I</sup>*<sup>3</sup> ¼ � <sup>1</sup>

q

�

*<sup>I</sup>*<sup>4</sup> <sup>¼</sup> <sup>1</sup> 3

� ln

**46**

q

�

*Geometry used to calculate the contribution of a spherical tesseroid (green rectangle) to the gravitational acceleration at a point P outside a spherically symmetric earth. a) Representation in the earth-Centered rotational reference frame (ECR). b) Representation in the earth-Centered P-rotational reference frame (ECP), defined with respect to the observation point P. R stands for the mean radius of the spherical earth; φ and λ stand for colatitude and longitude, respectively.*

#### **Figure 3.**

*Scheme of the UNIPOL approach. a) Set of tesseroids in the earth-centred rotational (ECR) reference frame. The yellow circlet indicates the observation point. The blue dashed circle around the observation point marks the area where tesseroids are at a distance* ≤*0.1° from the observation point and the ST procedure is required. Cyan and pink colors are used to indicate the tesseroids that require the application of the ST (cyan) and the RT (pink) procedure. b) Tesseroids mapped in the earth-Centered P-rotational (ECP) reference frame ECP and decomposition of some of them into sectors of a spherical band (light blue-colored tesseroids) in the local ECP reference frame. c) Tesseroids mapped in the earth-Centered P-rotational (ECP) reference frame ECP (pink areas) and re-orientation of some of them in the local ECP reference frame (red sectors). Modified from Marotta and Barzaghi [11].*

The second procedure (RT procedure) is based on the rotation of an ECRTesseroid (pink areas in **Figure 3c**) around its center and on its resizing to form a sector of a spherical band (red lines in **Figure 3c**) that develops along the local ECP-meridian and has the ECP-longitudinal dimension such that the ECPTesseroid maintains the same area of the original ECRTesseroid.

The *l*1*pq* and *l*2*pq* terms are the three-dimensional distances between *P* and the end points of the segment *pq*, while *s*1*pq* and *s*2*pq* are the one-dimensional distances between the origin *P*" of a 1D local coordinate system defined on the segment *pq* and its two end points. Last term of Eq. (11), *sin gAp* , is the singularity term that appears for specific locations of *P*' with respect to the closed polygonal line *Gp*. It expresses the analytical solutions of the corresponding limiting values of the line integrals that are obtained from the partial application of the divergence theorem for a small circle containing the singularity point when its radius tends toward zero. This singularity terms yield finally the values �2*πhp* when *P*' lies inside *Sp*, �*πhp* when *P*' is located on *Gp* but not at any of its vertices, �*ϑhp* when *P*' is located at one of the vertices of *Gp* and *ϑ* is the angle defined by two subsequent segments that meet at the corresponding vertex, and 0 when *P*' is located outside *Sp*. In Eq. (11), the two coordinate transformations occurring in every face of the polyhedron are given by two nested summations, firstly over faces and secondly over segments, of the same transcendental expressions depending on the vertices of the polyhedron. This means properly managing the relative position of points *P*, *P*' and *P*" with respect to every surface *Sp* that can be done algorithmically by standard tools of

*The Gravity Effect of Topography: A Comparison among Three Different Methods*

*DOI: http://dx.doi.org/10.5772/intechopen.97718*

vector calculus. **Figure 4** presents a sketch of their relative positions.

*Q <sup>j</sup> φ <sup>j</sup>*

, *λ <sup>j</sup>*,*rj*

dron will be: *A φ<sup>j</sup>*

, *λ <sup>j</sup>*,*ri* , *<sup>E</sup> <sup>φ</sup><sup>j</sup>*

*D φ<sup>j</sup>* þ *Δφ<sup>j</sup>*

**Figure 4.**

*point* P*.*

**49**

, the grid spacings *Δφ<sup>j</sup>*

, *λ <sup>j</sup>*,*ri*

, *<sup>B</sup> <sup>φ</sup> <sup>j</sup>*, *<sup>λ</sup> <sup>j</sup>* <sup>þ</sup> *Δλ <sup>j</sup>*,*ri*

, *λ <sup>j</sup>*,*rj* , *<sup>F</sup> <sup>φ</sup><sup>j</sup>*

Once clarified that the linear integral approach requires the control of the relative position of the vertices of a polyhedral mass with respect to the computation point, in the following a brief description of the implementation of such approach in the framework of residual terrain correction algorithms assessment is provided. Let us consider a set of *n* points in space *Pi* (*i* = 1 to *n*) on which the gravitational effect of the terrain must be computed and a set of *m* gridded points *Q <sup>j</sup>* (*j* ¼ 1 to *m*) with grid spacings *Δφ<sup>j</sup>* in latitude and *Δλ <sup>j</sup>* in longitude describing the DTM considered to compute this effect. The coordinates of *Pi* and *Q <sup>j</sup>* are expressed in terms of spherical coordinates ð Þ *φ*, *λ*,*r* referring to the same reference system, as well as the angular grid spacings *Δφ<sup>j</sup>* and *Δλ <sup>j</sup>*. On the basis of the spherical coordinates of

a six-facet polyhedron can be defined. In particular, the eight vertices of this polyhe-

*Sketch of the relative position of points* P*,* P' *and* P" *and of the local reference system defined in the computation*

, *<sup>C</sup> <sup>φ</sup><sup>j</sup>* <sup>þ</sup> *Δφ <sup>j</sup>*

, *λ <sup>j</sup>* þ *Δλ <sup>j</sup>*,*rj*

, *<sup>G</sup> <sup>φ</sup><sup>j</sup>* <sup>þ</sup> *Δφ<sup>j</sup>*

, *Δλ <sup>j</sup>* and the radius of the computation point *ri*,

, *λ <sup>j</sup>* þ *Δλ <sup>j</sup>*,*ri* ,

> , *λ <sup>j</sup>* þ *Δλ <sup>j</sup>*,*rj* ,

In agreement with [11], in the present study we use the ST procedure for angular distance between the observational point and the center of the tesseroid less than 0.1°. For angular distance greater than 0.1° we use the RT procedure.

Finally, the Flat Tesseroid (FT) approximation is considered. In this case, the bottom and the top of the tesseroid are flattened thus obtaining a particular polyhedron. A polyhedron is by definition a three-dimensional body with flat polygonal faces, straight edges and sharp corners or vertices and, through proper determination of its vertices, it can realize volumes approximating both right and spherical prisms.

By assuming the FT approximation, one can apply the linear integral algorithm for the computation of the gravitational effect implied by the topographic masses as proved in [4]. The implementation of such approach is based on the double application of the divergence theorem of Gauss to volume integrals of the general expressions for the gravity field resulting from a polyhedral constant density mass and evaluated at an arbitrary point in space. Given the potential of the body *B* in Eq. (7), the three components of the attraction of this body in a Cartesian orthogonal frame (*x, y, z*) can be expressed as

$$\mathcal{V}\_{\mathcal{X}\_i}(P) = \mathcal{G}\rho \iint\limits\_{\mathcal{B}} \frac{\partial}{\partial \mathbf{x}\_i} \left(\frac{\mathbf{1}}{l}\right) d\boldsymbol{v}(\mathbf{Q}) \, i = \mathbf{1}, 2, 3 \tag{10}$$

where *Vxi* ð Þ *P* is the partial derivative of the potential *V P*ð Þ along the *xi* direction.

The double application of the divergence theorem transforms these expressions first into an equal set of surface integrals (of the same number of faces of the polyhedral source) and subsequently each of them into a set of line integrals defined for each individual segment belonging on that face of the polyhedron. The solution of each line integral produces the final analytical formulas that, in terms of firstderivatives of the potential, are given by [4]:

$$V\_{\chi\_i}(P) = \mathbf{G}\rho \sum\_{p=1}^{n} \cos\left(\mathbf{N}\_p, \mathbf{e}\_i\right) \left[\sum\_{q=1}^{m} \sigma\_{pq} h\_{pq} LN\_{pq} + h\_p \sum\_{q=1}^{m} \sigma\_{pq} AN\_{pq} + \sin\left(\mathbf{g}\_{A\_p}\right)\right] \tag{11}$$

In Eq. (11), *p* defines one of the polygonal surfaces *Sp*, *q* defines one of the segments delimitating the *p* polygonal surface, *N<sup>p</sup>* is the outer unit normal of the polyhedral plane *p*, *e<sup>i</sup>* is the unit vector situated in the computation point *P*, *σpq* is equal to �1 when the normal of segment *q* lying on the plane of polygon *Sp* points to the half-plane that contains the projection *P*' of *P* on *Sp*, *σpq* is equal to +1 otherwise, *hpq* is the distance between *P*' and the segment *q*, *hp* is the distance between *P* and the plane *p*. *LNpq* and *ANpq* are transcendental expressions:

$$LN\_{pq} = \ln \frac{s\_{2pq} + l\_{2pq}}{s\_{1pq} + l\_{1pq}} \tag{12}$$

$$AN\_{pq} = \arctan \frac{h\_p s\_{2pq}}{h\_{pq} l\_{2pq}} - \arctan \frac{h\_p s\_{1pq}}{h\_{pq} l\_{1pq}} \tag{13}$$

## *The Gravity Effect of Topography: A Comparison among Three Different Methods DOI: http://dx.doi.org/10.5772/intechopen.97718*

The *l*1*pq* and *l*2*pq* terms are the three-dimensional distances between *P* and the end points of the segment *pq*, while *s*1*pq* and *s*2*pq* are the one-dimensional distances between the origin *P*" of a 1D local coordinate system defined on the segment *pq* and its two end points. Last term of Eq. (11), *sin gAp* , is the singularity term that appears for specific locations of *P*' with respect to the closed polygonal line *Gp*. It expresses the analytical solutions of the corresponding limiting values of the line integrals that are obtained from the partial application of the divergence theorem for a small circle containing the singularity point when its radius tends toward zero. This singularity terms yield finally the values �2*πhp* when *P*' lies inside *Sp*, �*πhp* when *P*' is located on *Gp* but not at any of its vertices, �*ϑhp* when *P*' is located at one of the vertices of *Gp* and *ϑ* is the angle defined by two subsequent segments that meet at the corresponding vertex, and 0 when *P*' is located outside *Sp*. In Eq. (11), the two coordinate transformations occurring in every face of the polyhedron are given by two nested summations, firstly over faces and secondly over segments, of the same transcendental expressions depending on the vertices of the polyhedron. This means properly managing the relative position of points *P*, *P*' and *P*" with respect to every surface *Sp* that can be done algorithmically by standard tools of vector calculus. **Figure 4** presents a sketch of their relative positions.

Once clarified that the linear integral approach requires the control of the relative position of the vertices of a polyhedral mass with respect to the computation point, in the following a brief description of the implementation of such approach in the framework of residual terrain correction algorithms assessment is provided. Let us consider a set of *n* points in space *Pi* (*i* = 1 to *n*) on which the gravitational effect of the terrain must be computed and a set of *m* gridded points *Q <sup>j</sup>* (*j* ¼ 1 to *m*) with grid spacings *Δφ<sup>j</sup>* in latitude and *Δλ <sup>j</sup>* in longitude describing the DTM considered to compute this effect. The coordinates of *Pi* and *Q <sup>j</sup>* are expressed in terms of spherical coordinates ð Þ *φ*, *λ*,*r* referring to the same reference system, as well as the angular grid spacings *Δφ<sup>j</sup>* and *Δλ <sup>j</sup>*. On the basis of the spherical coordinates of

*Q <sup>j</sup> φ <sup>j</sup>* , *λ <sup>j</sup>*,*rj* , the grid spacings *Δφ<sup>j</sup>* , *Δλ <sup>j</sup>* and the radius of the computation point *ri*, a six-facet polyhedron can be defined. In particular, the eight vertices of this polyhedron will be: *A φ<sup>j</sup>* , *λ <sup>j</sup>*,*ri* , *<sup>B</sup> <sup>φ</sup> <sup>j</sup>*, *<sup>λ</sup> <sup>j</sup>* <sup>þ</sup> *Δλ <sup>j</sup>*,*ri* , *<sup>C</sup> <sup>φ</sup><sup>j</sup>* <sup>þ</sup> *Δφ <sup>j</sup>*, *<sup>λ</sup> <sup>j</sup>* <sup>þ</sup> *Δλ <sup>j</sup>*,*ri* , *D φ<sup>j</sup>* þ *Δφj*, *λ <sup>j</sup>*,*ri* , *<sup>E</sup> <sup>φ</sup><sup>j</sup>* , *λ <sup>j</sup>*,*rj* , *<sup>F</sup> <sup>φ</sup><sup>j</sup>* , *λ <sup>j</sup>* þ *Δλ <sup>j</sup>*,*rj* , *<sup>G</sup> <sup>φ</sup><sup>j</sup>* <sup>þ</sup> *Δφ<sup>j</sup>* , *λ <sup>j</sup>* þ *Δλ <sup>j</sup>*,*rj* ,

#### **Figure 4.**

*Sketch of the relative position of points* P*,* P' *and* P" *and of the local reference system defined in the computation point* P*.*

The second procedure (RT procedure) is based on the rotation of an ECRTesseroid (pink areas in **Figure 3c**) around its center and on its resizing to form a sector of a spherical band (red lines in **Figure 3c**) that develops along the local ECP-meridian and has the ECP-longitudinal dimension such that the ECPTesseroid maintains the

In agreement with [11], in the present study we use the ST procedure for angular distance between the observational point and the center of the tesseroid less than

Finally, the Flat Tesseroid (FT) approximation is considered. In this case, the bottom and the top of the tesseroid are flattened thus obtaining a particular polyhedron. A polyhedron is by definition a three-dimensional body with flat polygonal faces, straight edges and sharp corners or vertices and, through proper determination of its vertices, it can realize volumes approximating both right and spherical

By assuming the FT approximation, one can apply the linear integral algorithm for the computation of the gravitational effect implied by the topographic masses as proved in [4]. The implementation of such approach is based on the double application of the divergence theorem of Gauss to volume integrals of the general expressions for the gravity field resulting from a polyhedral constant density mass and evaluated at an arbitrary point in space. Given the potential of the body *B* in Eq. (7), the three components of the attraction of this body in a Cartesian

0.1°. For angular distance greater than 0.1° we use the RT procedure.

*Geodetic Sciences - Theory, Applications and Recent Developments*

same area of the original ECRTesseroid.

orthogonal frame (*x, y, z*) can be expressed as

*Vxi*

derivatives of the potential, are given by [4]:

cos *Np*, *e<sup>i</sup>*

� � <sup>X</sup>*<sup>m</sup>*

the plane *p*. *LNpq* and *ANpq* are transcendental expressions:

*ANpq* ¼ arctan

*q*¼1

X*n p*¼1

ð Þ¼ *P Gρ*∭

*B*

*∂ ∂xi*

first into an equal set of surface integrals (of the same number of faces of the polyhedral source) and subsequently each of them into a set of line integrals defined for each individual segment belonging on that face of the polyhedron. The solution of each line integral produces the final analytical formulas that, in terms of first-

1 *l* � �

The double application of the divergence theorem transforms these expressions

*σpqhpqLNpq* þ *hp*

In Eq. (11), *p* defines one of the polygonal surfaces *Sp*, *q* defines one of the segments delimitating the *p* polygonal surface, *N<sup>p</sup>* is the outer unit normal of the polyhedral plane *p*, *e<sup>i</sup>* is the unit vector situated in the computation point *P*, *σpq* is equal to �1 when the normal of segment *q* lying on the plane of polygon *Sp* points to the half-plane that contains the projection *P*' of *P* on *Sp*, *σpq* is equal to +1 otherwise, *hpq* is the distance between *P*' and the segment *q*, *hp* is the distance between *P* and

*LNpq* <sup>¼</sup> ln *<sup>s</sup>*2*pq* <sup>þ</sup> *<sup>l</sup>*2*pq*

*hps*2*pq hpql*2*pq*

*s*1*pq* þ *l*1*pq*

� arctan

*hps*1*pq hpql*1*pq*

X*m q*¼1

� � " #

ð Þ *P* is the partial derivative of the potential *V P*ð Þ along the *xi* direction.

*dv Q*ð Þ *i* ¼ 1, 2, 3 (10)

*σpqANpq* þ sin *gAp*

(11)

(12)

(13)

prisms.

where *Vxi*

*Vxi*

**48**

ð Þ¼ *P Gρ*

*H φ <sup>j</sup>* þ *Δφ <sup>j</sup>*, *λ <sup>j</sup>*,*rj* . Note that the two planar surfaces identified by the closed polygons *ABCD* and *EFGH* have their vertices whose radiuses depend in the first case on the radius of the computation point *Pi* while in the second case on the radius of the DTM point *Q <sup>j</sup>* . **Figure 5** illustrates in graphical form how the single polyhedron is built.

consistency between the different *Vr* computed by the same DTM point *Q <sup>j</sup>* with

*The Gravity Effect of Topography: A Comparison among Three Different Methods*

**3. The numerical tests on the three proposed approaches**

(called AREA\_1):

inner area:

**Figure 6.**

**51**

*The DTM (AREA\_1) and the points for TC computation.*

The three different mathematical models presented in Section 2 have been applied in two *TC* computation tests. Both tests have been carried out in the Alpine area. In the first test, the SRTM3 DTM (see [14]) have been selected in the area

46° ≤φ≤47° 11° ≤λ≤12°

46*:*25°≤φ ≤46*:*75° 11*:*25°≤ λ≤11*:*75°

The computation points are thus in an inner area which have 0.25° from the extents of the outer one containing the DTM. So, the terrain correction in points

The statistics of the height data in this area (see **Figure 6**) are given in **Table 1**. Points for the *TC* computation have been selected on a 3<sup>0</sup> � 3<sup>0</sup> regular grid in the

The three different models adopted in computing the *TC* account for a different geometry of the single terrain element. The tesseroid formula (Eq. (9)) considers the radial convergence component of the vertical edges of each terrain element and, besides that, gives a spherical approximation of the top and the bottom of this element. The flat tesseroid has the same geometry in the radial component but top and bottom are planar surfaces. Thus, its geometry is, in principle, less accurate than the one of the tesseroid. Finally, the computation based on the right parallelepiped disregards also the radial convergence of the vertical edges. Hence, this model describes the topography geometry of each element in a way that does not adhere properly the two main features of the given DTM element. However, since in computing the *TC* effect only the differences between the height of the computation point and the heights of the DTM needs to be considered, the above-mentioned differences between the geometries of the elements used in computing the *TC* should have a limited impact on the computed values. To prove this, we have devised a test in the Alpine area where the largest discrepancies are expected in the *TC* effect evaluated with the three different geometries of each single DTM element.

respect to the different computation points *Pi*.

*DOI: http://dx.doi.org/10.5772/intechopen.97718*

This procedure defines the spherical coordinates of the vertices of the polyhedron that are at the height of the terrain and at the height of the computation point, i.e. the level at which the Bouguer plate is computed. The computation of the firstorder derivative along the direction of *ri* of the gravitational effect of such polyhedron corresponds to the terrain correction to be applied at *Pi* as contribute of *Q <sup>j</sup>* . Such value is obtained by running the code *polyhedron.f* made available by the author [4] implementing the linear integral approach. As input, the relative cartesian coordinates of the polyhedron vertices with respect to the computation point and their topological relationships are required. These are obtained applying a change of reference system to the polyhedron vertices. In particular, their coordinates were roto-translated into a local reference system having the origin at the computation point *Pi*, the *z* axis pointing up along the direction of *ri* and the *x* and *y* axes parallel to the local East and North directions, respectively. Regarding the topological relationships defining the outer normal direction of the six planes of the polyhedron, they are defined by a topology matrix containing the counterclockwise sequence of the vertices as seen from outside. As output, the absolute value of the computed *Vr* is taken. This procedure contemplates two nested loops over all the *m* DTM points *Q <sup>j</sup>* and the *n* computation points *Pi*. Within the loop over the computation points, different local reference systems are defined. This leads to slight changes in the directions of the *x* and *y* axes but not on the *z* axis, always normal and pointing outside the reference sphere defined on *Pi*, then maintaining

#### **Figure 5.**

*Sketch of the polyhedron vertices building procedure on the basis of the DTM point* Qj *and the computation point* Pi *.*

*The Gravity Effect of Topography: A Comparison among Three Different Methods DOI: http://dx.doi.org/10.5772/intechopen.97718*

consistency between the different *Vr* computed by the same DTM point *Q <sup>j</sup>* with respect to the different computation points *Pi*.

The three different models adopted in computing the *TC* account for a different geometry of the single terrain element. The tesseroid formula (Eq. (9)) considers the radial convergence component of the vertical edges of each terrain element and, besides that, gives a spherical approximation of the top and the bottom of this element. The flat tesseroid has the same geometry in the radial component but top and bottom are planar surfaces. Thus, its geometry is, in principle, less accurate than the one of the tesseroid. Finally, the computation based on the right parallelepiped disregards also the radial convergence of the vertical edges. Hence, this model describes the topography geometry of each element in a way that does not adhere properly the two main features of the given DTM element. However, since in computing the *TC* effect only the differences between the height of the computation point and the heights of the DTM needs to be considered, the above-mentioned differences between the geometries of the elements used in computing the *TC* should have a limited impact on the computed values. To prove this, we have devised a test in the Alpine area where the largest discrepancies are expected in the *TC* effect evaluated with the three different geometries of each single DTM element.
