Computing Local Geoid Model Using DTM and GPS Geodetic Points. Case Study: Mejez El Bab-Tunisia

Different methods have been deployed to compute the geoid, the altimetry reference for surveying applications. One of their main goals is to allow the use of GPS (Global Positioning System) or GNSS heights, which are related to an ellipsoid and therefore must be corrected. Some of these methods are accurate but quite heavy as developed by [1], but one of them is easy to use while giving very good results in a local system: some mm for a 10 × 10 km area developed by [2] [3]. In our study, we have used software called “Géoide Program”, previously used at the CERN in Switzerland and set up by [4], which they complete this software allowing a parameterization of general data to provide results in a general system. Then, tests have shown the way to optimize computations without any loss of accuracy. For our computations we use gridded of geodetic heights, from Lambert or WGS 84 datum’s, DTM (Digital Terrain Model) and leveled GPS points. To obtain these results, components of the vertical deflection are computed for every point on the grid, deduced from the attraction exerted by the mass Model. Then, geodetic heights are computed by an incremental way from an arbitrary reference. Once the calculation is performed, the geodetic height of any point located in the modelled area can be interpolated. The variations of parameters (mainly size and increments of the DTM and of the modeled area, and ground density) have shown that they do not play a significant role although DTM must be large enough to take into account an important area around a selected zone. However, the choice of the levelled GPS points is primordial. We have performed tests with real data concerning Mejez El Bab zone, in north of Tunisia. Nevertheless, for a few hundreds of square kilometers area, and just by using a DTM and a few levelled GPS points, this method provides results that look How to cite this paper: Rebaï, N., Zenned, O., Trabelsi, H. and Achour, H. (2018) Computing Local Geoid Model Using DTM and GPS Geodetic Points. Case Study: Mejez El Bab-Tunisia. International Journal of Geosciences, 9, 161-178. https://doi.org/10.4236/ijg.2018.93011 Received: December 29, 2017 Accepted: March 10, 2018 Published: March 13, 2018 Copyright © 2018 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access


Introduction
Knowledge about the shape of the Earth and the variations of its field of gravity is now crucial for any precise work in the spatial domain of artificial satellites trajectory, machines location and astronomy of position. It also allows the elaboration of geophysical hypotheses concerning the inside of the globe [5] [6] [7].
Since the geoid of a good approximation was a reference of the heights until the end of the 80's, the geoid has remained an object of scientific study. The fast development of the GPS system and the necessity of converting the ellipsoidal elevations into heights have made the geoid an important stake [3] [8].
Knowledge about the geoid is essential for the levelling because it is considered as the reference surface [9] [10] [11]. The shape of the geoid depends on the distribution of the masses inside the globe. Indeed, the fact that the Earth is not a homogeneous entity makes the geoid present undulations with regard to the ellipsoid. These undulations reflect the different aspects of heterogeneous densities [12].
The determination of a Tunisian geoid of centimetric accuracy turns out to be more and more essential. However, the gravimetric data, necessary for its determination, are very scattered. Their homogeneity, density and quality are unknown [13] [14].
Thus, the main objective of this work is the determination of local geoid model of Mejez El Bab using the DTM (Digital Terrain Model) of the sector and levelled GPS points. It is, in fact, about a quasi-geoid model. The purpose of this work is, also, to master the methodology of the determination of a geoïd in order to spread this calculation all over Tunisia.

Geographical and Geological Context of the Study Area
The choice of the study area was made according to several criteria in particular: its chaotic relief, the knowledge of its geology and a precise DTM obtained from the interpolation of the contour line of topographic map (scale 1:50,000).
The region of Mejez El Bab is situated in northern Tunisia (Figure 1). The studied sector is part of the "zone of domes" [15] [16] [17] [18]. This sector is of interest to the structures that become integrated into the front-country of the alpine chain, situated on the borders of the valley of Mejerda [19] [20] [21]. In this  zone the tectonic structure is characterized by Triassicoutcrops aligned according to a NE-SW direction [22] [23] [24]. On the paleogeographic plan, the "zone of domes", which is of about 50 km, corresponds to one bet of the qualified deep marine domain known as the "Tunisian groove". This domain presents a complex paleogeographic and structural evolution, bound to the multiple periods of opening and closing, accompanied with Triassic halokinetic movements [25] [26].
The studied structures in Mejez El Bab, Soughia sector correspond to a mega-anticline Plio-Quaternary resuming previous anticline and synclinal structures.
Indeed  [11]. In consequence, the corrections and the orthometric heights are not physically determined. Thus, they can be calculated from the moderate raw, measuring the differences of raw level, only by means of theoretically approached formula giving the gravity g in each point of the globe.
The obtained values are only approximate orthometric corrections. Resulting in the fundamental mark, the origin of the heights of the network is the mark of "Port de France", which was built 60 years ago, and having a new height of 7 m.
The 10 cm correction with regard to the former network was necessary for following the information provided by the tide gauge of la Goulette [33]. However, it is not an obvious solution because this mark is situated in a geologically unstable ground that is far away from the site of the tide gauge of la Goulette.
The Lambert conical conform projection used in the geodetic system. In order to avoid serious distortions, Tunisia has been divided into two parts: Lambert North of Tunisia (11 g, 40 g) and Lambert South of Tunisia (11 g, 37 g). The extension of the geodesic network southward requires a 3rd zone Lambert that can be characterized by (11 g, 34 g). In fact, it is important to draw the attention on the fact that we have two types of coordinates system in Lambert projection used by the Topographic Service of Tunisia, called "STT" coordinates are as follows: X has a positive sense northward and Y has a positive sense westward.
With X, Y counted from the origin of the projection axes in which X 0 = 0 and Y 0 = 0.
The formulas allowing the shift from "IGN" coordinate system to "STT" (shift of the origin in meters, 300,000, 500,000) are:

Methodology
For the calculation of the Mejez El Bab Local geoid model, we principally used two methods. The first principle method consists to the integration of the vertical deflection (ξ and η˚). The second method is based on the calculation of the potential engendered by the DTM of the area.
Our aim in this study is to elaborate "Géoide PROGRAM". It allows us to calculate geodetic heights in points of a grid defined by the user, using DTM de-  The attraction of a prism at a point of the grid is calculated by way of two formulas: The first (F1) consists of a direct application of Newton's law: where r is the distance between two bodies (m 1 , m 2 ), V the volume of the prism supposed constant and ρ is the density.
The second which is more rigorous, is obtained by integrating the first one on the volume of the prism: the integration gives (F3): ( ) The formula F1 is used until a critical distance beyond which differences of results between both formulae are unimportant, and where the formula F3 is thus applied. After several attempts this critical distance was fixed to 50,000 m because when it was considered superior, the calculation time increased strikingly without significant variation of the result.
If the zone has a more heterogeneous ground, the choice of the critical distance could have more consequences. Since the attraction is being inversely proportional to the distance, it is the closest prisms which have a dominating place in the calculation of the attraction and thus the vertical deflection. Yet, differences between the attraction's calculation, according to the rigorous formula or not, are more important as the relief is pronounced.

Calculation of the Components of the Vertical Deflection ξ and η˚
In the studied area, the attraction t of every prism is calculated on the point.
Dg is the sum of the attractions t: Dg t = ∑ The vertical deflection according to the X axis is given by: cos sin sin cos Yet, since the acceleration of the gravity g is hardly measurable, Comet and Maillard [4] replaced it by γ (The normal acceleration).

Calculation of the Geodetic Height
The difference of height between two successive points of the grid is obtained through multiplying the distance between both points by the value of the vertical deflection in the direction obtained from both points. Indeed, in a point the angle between the normal and the vertical line is equal to the one between the tangent in the ellipsoid and the tangent in the geoid [4].
In the consideration of the geoid as a line between two consecutive points of the grid the geodetic heights can be so calculated step by step. This method is applied for general leveling and not only in local leveling coordinate system, the case of "Géoïde" Soft developed by CERN (Figure 4). The modelling was made possible from a DTM with rectangular grid, in "Lambert IGN" or in WGS 84.
Once the components of the vertical deflection (ξ, η) are calculated, a wedging is made to add the undulations which are not caused by the relief. It consists in a translation and a rotation, and is made by leveled GPS points, for which the geodetic heights are known by subtraction [42].
The principle of the wedging is to overlay the geodetic height points calculated by the model and leveled GPS points in the same geodetic system.
For this, a fundamental point will be used; it is going to be made a vertical translation and a rotation considering it as chosen in the middle of the Calculated geoid model in order to avoid losing precision too important for boundaries zone.
Indeed, an imprecision in the determination of the rotation parameters has a proportional effect on the height of the point to the distance between it and to the fundamental one. If this point is taken in the middle, the maximal distance is reduced and therefore allows minimizing errors [4].
For wedging the geoid Model, we calculate the value of the vertical translation dN 0 and two rotation angles according to the direction of ξ and η˚: dξ 0 et dη 0 .
Thus, we calculate at first step as described previously, the geodetic heights of all the points of the grid and the wedging points.
During the calculating, we have chosen N 0 . The latter is the approached value: It is N of wedging point that is the closest to the fundamental point. Other parameters dξ 0 and dη 0 are fixed at first to 0.
Then, the least mean square method is applied with three unknown parameters (dN 0 , dξ 0 , dη 0 ) and a number of relations of observations equal to those of the wedging points. The relations of observations are explained in [4].

The Potential Calculation Method
We have based our calculation of the gravitational potential on a formula used by [42]. This method consists in the calculation of the potential of gravity. Indeed, a model made up of n prisms can result in a potential T. This allows the conclusion that a constant volume generates a constant potential.

Case Study: Mejez El Bab Zone, North of Tunisia
The present work leads to the elaboration of a "GEOIDE Soft" presented in (Figure 4). This program enables the design of a first prototype of local geoid model in Tunisia, and precisely in Mejez El Bab.
This program allows the modelling of a local surface of height 0 at a grid by use of a DTM, and the wedging of this latter by means of levelled GPS points. It was so expressed in the Tunisian parameters. Besides, we added to this program the second method based on the calculation of the potential of gravity [43] [44] [45].
Every method has its own specificities and scopes. For the case of Tunisia which has a system of orthometric height, we were able to make the wedging only for the method of integration of the components of the vertical deflection.
Indeed, the method of the Potential is based on normal heights, but it has allowed us to validate the model of geoid obtained and to make a comparison between two different methods ( Figure 5). In both cases, we obtain a geoid of a centimetric precision. This one reproduces globally the relief of Mejez El Bab.
In all what follows, the used data are plotted using the curve line functions of "Surfer Software" to modelling the relief.
What is represented in the ( Figure 6) is the entire zone covered by the DTM in "Lambert North IGN". What is worth noting is that Southwest half is made up of a rather strong relief that has an influence on the geoid.

Results and Discussion
The obtained model represents a geoid that reproduces globally the relief, as shown in the (Figure 7). Indeed, we discern in its South West part a bump of direction NE-SW where the geodetic height reaches 0.03 meters, bounded by two hollows.
This bump would correspond to Jebel Bou Mous, Jebel Mourra, Jebel Jebs and Jebel Kechtilou (Figure 1), whereas the two hollows on both sides would represent the plain of Mejerda to the North (between 0.05 m and 0.03 m) and that of Goubellat in the South (between −0.03 m and 0.01 m). Also, in the southeast corner, we notice a light rise of the geoid reaching 0.02 m corresponding to Jebel Khamfir and Jebel Ed Deridjah. Whereas in the extreme North-West, the geoid rises again but achieves only −0.01 m in 0 m. On the map this rise corresponds to Jebel Jedidi and Jebel Chaabane.
The model obtained by the potential method ( Figure 8) is a geoid which slightly reproduces the relief but by being too much inflated. Hence, we realized   that is necessary to subtract the average value from the potential generated by the whole DTM in the middle point. We obtain then the geoid of the Figure 5 with a minimum of 0.07 m and a maximum of 0.06, corresponding at 13 cm of undulations.
The tests cover the influence of the various parameters (size and stitch of the DTM and the zone). In every test, a parameter is adapted to a standard, and the differences on the obtained heights are then studied.
The objective of these experiences is to compare the various models with what is in priori the best. The comparison is a quadratic gap average, and must be For all the parameters, we note that their choice is bound, for the majority of them, to be bound to the relief of the considered region. A precise analysis of the reliefs surrounding the zone must thus be made. The size of the DTM is fundamental as far as a bigger DTM is needed, but it is not enough. It is useful that the limits of the DTM are fixed by taking into account the relief with a maximum integration of the stressed parts [4]. The stitch of the DTM does not seem, as far as it is concerned, to be the most sensitive parameter. However, the realized tests converge on the choice of 200 m ( Table 1; Table 2). For the stitch of the zone to be modelled, better results are obtained with a reduced stitch ξ and η. Seen the limit of the calculation, the choice of 500 m seems suitable.
Regarding the points of wedging, their choice is determining for the final precision (Table 3). Seen their limited number, we are not able to test their influence on the modelling of the geoid. However, the more this number decreases, the more the position in the zone and the quality of points become more important during the wedging of the modeled surface. A compromise between the minimization of the number of points and the precision of the result must thus be reached.

Conclusions
The

North of Tunisia.
A prospective improvement of the model would be possible by certain techniques. First, varying the density on the DTM while considering a chaotic terrain has important implications. Second, having precise DTM is helpful. Third, the size of the DTM plays a crucial role. Therefore, it could be interesting to increase its size without too much increasing the mass of the calculations by using a much bigger stitch in periphery. The "Géoïde Soft", being for the moment limited in capacity, can't handle a DTM of more than 178*178 knots. These limits could be exceeded by finding the way to develop more the programs without amplifying the calculation. Also, it would be necessary to arrange more leveled GPS points with the maximum of accuracy, which requires a precise leveling.
Through this example, by using the data in the best way, it appears that rms on the obtained geodetic heights is better than 3 cm, which is close to GPS altimetry precision in z. Therefore, the total rms might be perhaps better, which could be confirmed with more accurate GPS measurements.
Nevertheless, for a few hundreds of square kilometers area, and just by using a DTM and a few levelled GPS points, this method provides results that look extremely promising, at least for surveying activities, as it shows a good possibility to use GPS for coarse precision levelling, and as DTM is now widely available in many countries.
The Terrestrial geoids are so difficult to be interpreted as they represent the effect resulting from the heterogeneity existing in the various internal layers of the Earth having complex structure and dynamics, which interact with each other through mechanical or thermal couplings, etc. However, by setting more data and ways, an alliance geoid-gravimetric would be interesting [41]. It could help us to emit various hypotheses and to interpret the heterogeneity of the interior of the earth.
Indeed, we consider in this study that the attraction value or the integration of vertical deflection components is directly proportional to the density. Therefore, an increase in density accentuates the relief. In our case, we have considered that the density is constant and equal to the continental crust (2.67). In order to have a better definition of the geological structures, it is necessary to refine the gravimetric models according to the Regional scale. These models will contribute, in a significant way, to the realization of a Tunisian geoid of high accuracy. Indeed, it would be necessary to integrate various methods to reach the required precision (gravimetric collocation) [44] [45].
Today, the essential results are known by the global models of the geopotential. Such models integrate all the information collected from artificial satellites, measures of gravity and, for several years, even altimetry spatial measures.