Accurate Height Determination in Uneven Terrains with Integration of Global Navigation Satellite System Technology and Geometric Levelling: A Case Study in Lebanon

Abstract

The technology for determining a point’s coordinates on the earth’s surface using the global navigation satellite system (GNSS) is becoming the norm along with ground-based methods. In this case, determining coordinates does not cause any particular difficulties. However, to identify normal heights using this technology with a given accuracy, special research is required. The fact is that satellite determinations of geodetic heights (h) over an ellipsoid surface differ from ground-based measurements of normal height (H^N) over a quasi-geoid surface by a certain value called quasi-geoid height or height anomaly (ζ). In relation to determining heights of a certain territory, the concept of geoid height (N) is usually operated when dealing with a geoid model. In this work, geodetic and normal heights are determined for five control points in three different regions in Lebanon, where measurements are carried out using GNSS technology and geometric levelling. The obtained quasi-geoid heights are compared with geoid heights derived from the global Earth model EGM2008. The results obtained showed that, in the absence of gravimetric data, the combination of global Earth model data, geometric levelling for selected areas, and satellite determinations allows for the creation of a highly accurate altitude network for mountainous areas.

1. Introduction

Geometric leveling still remains the most reliable method for high-precision determination of heights. This method involves direct measurements of elevations, from which normal heights are determined. A normal height system relates heights to a hypothetical surface called the quasi-geoid that closely approximates mean sea level, while a geodetic height system is based on the reference ellipsoid which approximates the Earth’s shape (Figure 1).

The active introduction of satellite determinations into geodetic practice was not only an alternative to geometric leveling, but also in many ways became its replacement. The main disadvantage of satellite leveling is the indirect determination of values from geodetic heights. In this case, a geoid/quasi-geoid model should be used, which includes a set of geoid/quasi-geoid heights above the ellipsoid on a regular grid that makes it possible to determine the difference between geodetic and normal heights (or orthometric heights in case of using a geoid model) anywhere on the planet.

The different magnitudes of gravity on Earth are mainly due to the uneven distribution of mass (rocks), resulting in different gravitational forces that ultimately create the current shape of the geoid. The task arises of identifying the degree of this heterogeneity (height anomalies, differences in the shape of the quasi-geoid from the ellipsoid), which in particular can be achieved for a certain area of the earth’s surface by constructing local models of the quasi-geoid [1]. Local models of the quasi-geoid are essential for understanding the degree of heterogeneity in gravity, particularly in uneven terrains, where variations in mass distribution significantly impact gravitational forces and the shape of the geoid.

Uneven terrains refer to areas characterized by variations in elevation, rugged terrain features, and irregular topography. These areas often exhibit significant changes in slope, elevation, and relief, presenting challenges for accurate height determination and surveying. Uneven terrains can include mountainous regions, hillsides, valleys, and areas with complex geological formations. In such environments, traditional surveying methods may encounter difficulties due to the uneven distribution of gravitational forces and irregular surface features. Consequently, precise positioning and height accuracy become crucial considerations in engineering and geodetic applications conducted in these areas.

Currently, through GNSS technology, the geodetic coordinates (φ, λ, h) of a point of interest on Earth are effectively determined. However, the exact values of the normal heights (H^N) of these points should still be determined on the basis of geometric leveling [2,3]. The object of this study is the territory of Lebanon, where there is neither an official national geoid model, nor gravimetric data, and accordingly, there is no state altitude network. In this regard, the relevance of developing a method for determining normal heights is extremely high. Local quasi-geoid models can be determined by gravimetric, astro-geodetic, or geometric methods. Publications [4,5,6] propose various methods for constructing local quasi-geoid models and assessing their accuracy. In the present study, normal heights are determined by a geometric leveling method using the three-route levelling technique.

By default, GNSS is integrated with global models such as EGM2008. The specific choice of EGM2008 was mainly due to its wide application and high resolution [7,8,9]. Although the latest models of the global gravitational field, such as EIGEN 6C4, EIGEN-6C2, and EIGEN-6C3stat, are considered more accurate, the practical results of their use (K.I. Markovich, D.Sh. Fazilova, M. Szelachowska) show identical results when comparing models, or a difference within a few millimeters when determining height anomalies. In this regard, the use of EGM2008, due to its extensive global validation in various settings, seems appropriate [10,11,12].

This study presents the determination of quasi-geoid height using geometric levelling and satellite measurements in three zones having different topographic characteristics in Lebanon. GNSS measurements in the static mode provide accuracy in the calculation of geodetic heights (within a few millimeters), which allows for an assessment of geoid waviness, and therefore for more accurate calculations of normal heights [13,14,15]. It should be noted that in practical activities (for example, in construction), GNSS measurement technologies are widely used today, which significantly implies its high accuracy in determining coordinates [16,17,18].

2. Materials and Methods

The geometric method was utilized in three distinct zones in Lebanon to determine height anomalies. Zone A (Bekaa governorate), Zone B (Mount Lebanon governorate), and Zone C (North governorate). These sites were selected based on their contrasting topographies. In each site, a total of five points were identified. Four of these points were arranged to form a quadrant, while the fifth point was allocated as a central point. These points were strategically positioned to encompass the surrounding area and capture any potential variations in the deviation of the plumb line.

GNSS measurements were first conducted in this work, where control points from different orders were taken as reference points. The triangulation networks in Lebanon include various orders, starting with the first-order polynomial network, where the side of a triangle ranges from 20 to 30 km. The second-order polynomial network relies on the first-order network, with triangle sides of about 20 km. The third-order polynomial network was created based on the first two networks, with triangle sides approximately 5 km in length. Finally, the fourth-order polynomial network, which is created based on the three previous networks, is convenient for identifying survey points throughout the country (Figure 2).

Static mode GNSS positioning was chosen over the real-time kinematic (RTK) mode due to its capacity to achieve superior accuracy over extended observation periods. Unlike RTK, which relies on continuous correction data from nearby base stations, static mode enables longer observation times, allowing receivers to accumulate signals from multiple satellites and attain highly precise positioning solutions, particularly for applications requiring centimeter-level accuracy across larger areas.

Topcon Hiper-V GNSS receivers with their controllers were used for measuring geodetic heights. Tripods and tribraches have been used to set up the receivers over the points. The instrument height was measured twice at the beginning and the ending of each session using a regular measuring tape. During each session, four satellite receivers were positioned and operated for a duration of one hour. This duration allows for sufficient satellite signal acquisition and data collection for precise positioning. Figure 3 shows the distribution of points in each zone. All GNSS measurements started from a chosen reference point from second- and third-order networks in each study area. In (A), the reference point 8F3Z from the third-order network was allocated. A_1–4 are the ancillary points that formed a quadrant with a central point CP. The created baselines between these points were measured starting from 8F3Z. The same procedure was followed in zone (B), where three available control points in the region were allocated, and measurements were carried out based on them. The control points 7MR5, 7JLB, and 7PRP are from the third-order network. However, in zone (C), the chosen control point was QDF7 from the second-order network.

The raw data obtained in the three zones have been processed and adjusted to obtain the geodetic coordinates and variances of the stations. The processing interval was set as automatic to determine the best interval setting to use based on the length of the baseline and the duration period. This option balances the highest processing efficiency with the highest results quality.

After processing the baselines in the three zones, the control points were used to perform the fully constrained adjustment. Fixing the control points shifts the observations to the correct location within the chosen datum. Then, the adjusted coordinates for all other points in the network can be determined with respect to the project datum. The adjusted geodetic coordinates are shown in

Table 1. Geodetic coordinates determined by GNSS static measurements.

Point	φ	λ	h (m)
Zone (A)
A1	N 33°31′04.23088″	E 35°40′19.02209″	879.954
A2	N 33°31′01.24131″	E 35°40′19.23139″	881.054
A3	N 33°31′00.71578″	E 35°40′05.80162″	885.348
A4	N 33°31′05.25177″	E 35°40′08.02764″	885.476
CP	N 33°31′01.69385″	E 35°40′13.54181″	882.446
Zone (B)
A1	N 33°43′39.48678″	E 35°28′28.34842″	276.318
A2	N 33°43′10.03531″	E 35°28′17.11964″	242.384
A3	N 33°43′04.36187″	E 35°27′42.56322″	152.402
A4	N 33°43′37.40257″	E 35°27′53.33982″	193.289
CP	N 33°43′22.26857″	E 35°28′00.94609″	192.591
Zone (C)
A1	N 34°36′52.51739″	E 36°07′53.08944″	212.957
A2	N 34°36′39.53158″	E 36°08′13.63695″	193.228
A3	N 34°36′12.38280″	E 36°08′00.49453″	236.286
A4	N 34°36′34.61689″	E 36°07′08.98025″	163.769
CP	N 34°36′34.27561″	E 36°07′55.86680″	164.500

.

Normal heights are usually determined by the geometric levelling method, which is verified by taking measurements in two directions: forward and backward. However, the uneven distribution of the Earth’s mass leads to an uneven gravitational field, which affects the waviness of the geoid [19,20,21]. This may appear in different values during the geometric levelling measurements along different routes. If the values turn out to be the same, then the gravitational field in the area of study can be considered negligible. To determine whether there are noticeable changes due to gravitational forces, measurements are taken between the baselines along three different routes, and then the elevation differences are compared. The process of measuring heights through the three different routes (one straight route and two curvilinear routes) occurs between two main points, passing through an intermediate point in each followed curvilinear route (Figure 4). For example, the points A1 and CP are measured three times, where back sight and fore sight readings are taken along the first route (straight: baseline-B), then along the second route (curvilinear: baselines 2 and 5) passing through the intermediate point IP2, and then along the third route (curvilinear: baselines 21 and 24) passing through IP1. The same process is followed for other points. This process helps in estimating the gravitational influence in the area by comparing the values of the measured heights in each route.

Systematic errors such as collimation error, non-standard temperature error, earth curvature and atmospheric refraction, and disclosure error are eliminated [22,23,24]. The corrections of these errors are applied to the raw observations in each zone to indicate the adjusted elevation difference between the central point and the surrounding control points.

• Collimation error:

$$e=\frac{\left(\mathrm{B}\mathrm{S}2-\mathrm{F}\mathrm{S}2\right)-(\mathrm{B}\mathrm{S}1-\mathrm{F}\mathrm{S}1)}{2},$$

(1)

where e is the total vertical error in a horizontal sight, BS is the backsight reading, and FS is the foresight reading.

• Non-standard temperature error:

$$C_t=\left(t_m-t_s\right)\times D\times CE,$$

(2)

where $C_t$ is the rod temperature correction, $t_m$ is the mean observed temperature of the rod, $t_s$ is the standardization temperature of the rod, $D$ is the observed difference of elevation between the benchmarks, $CE$ is the mean coefficient of thermal expansion per unit length per degree temperature of the rod.

• Earth curvature and atmospheric refraction:

$$C_{C\&R}=0.0673\times D^2,$$

(3)

where D is the distance from the instrument to the staff station in kilometers.

• Disclosure error:

$$e_D=\sum BS-\sum FS,$$

(4)

Assuming that the height of the starting point in each zone is 100 m, the elevation differences are calculated between all points, and adjusted based on the errors’ corrections (see

Table 2. Correction of systematic errors and applying LSM in Zone (B).

Section	Route	Correction of Systematic Errors (mm)			Disclosure (m)	ΔHN (m)	LSM (m)
		Collimation	Earth Curvature and Atmospheric Refraction	Non-Standard Temperature
CP-A1	1	−0.268	0.312	−0.010	−0.001	−48.557	−48.549
	2	0.175	1.137	−0.009	−0.001
	3	−0.238	1.170	−0.007	0.000
CP-A2	1	−0.280	0.358	−0.004	−0.001	−28.732	−28.737
	2	0.243	1.190	−0.005	−0.001
	3	0.331	1.575	−0.007	0.000
CP-A3	1	0.050	0.315	−0.025	−0.002	−71.820	−71.816
	2	0.231	1.349	−0.019	0.005
	3	0.140	1.049	0.009	0.005
CP-A4	1	−0.686	0.517	0.000	0.001	0.576	0.573
	2	0.173	1.148	0.001	−0.004
	3	0.289	1.385	0.001	0.001

). After determining the adjusted height difference, the least squares method (LSM) was applied in order to minimize the errors as much as possible. The LSM is expressed in the Equations (5)–(11).

The weight matrix for an equivalent observation is as follows:

$$W_e=[A.Q.A^T{]}^{-1},$$

(5)

where A is the Jacobian matrix, and Q is the cofactor diagonal matrix representing the relative variances of the observations.

The unknown parameters are calculated by the following:

$$\Delta ={{[B}^T{.W}_e.B]}^{-1}{[B}^T{.W}_e.f],$$

(6)

where B is the design matrix, and f is the disclosure vector.

Since this is a nonlinear equation system, the corrections in matrix $(\Delta )$ are applied to the initial approximations, and the method is repeated until the system converges.

The equivalent residuals $V_e$ are as follows:

$$V_e=B\Delta -f,$$

(7)

The observational residuals are as follows:

$$V=Q{A}^T{W}_e{V}_e,$$

(8)

The reference variance for the adjustment can be computed using the equivalent residuals and weight matrix employing the equation:

$${\sigma }_0^2=\frac{V_e^T.W_e.V_e}{r},$$

(9)

where r is the number of redundancies in the system and defined as $r=n-m$.

The covariance matrix:

$${\mathsf{\Sigma }}_{\Delta \Delta }={\sigma }_0^2.{{[B}^T.W_e.B]}^{-1}={\sigma }_0^2.Q_{\Delta \Delta }$$

(10)

The standard deviation of the individual adjusted parameters is obtained from the following:

$${\sigma }_i={\sigma }_0.\sqrt{q_{\mathrm{i},\mathrm{i}}},$$

(11)

where $q_{\mathrm{i},\mathrm{i}}$ is ith the diagonal element of the covariance matrix $Q_{\Delta \Delta }$.

The results of the adjusted elevation differences (ΔH^N) and the applied LSM for zone B are shown in Table 2. The same procedure was followed for both zones A and C.

In order to evaluate the achieved accuracy using the proposed method (three-route levelling), the same adjustment process took place for the levelling measurements taken at the first route only (direct route between two points), excluding the two curvilinear routes, which is exactly the known classic levelling method. After correcting all the systematic errors, the LSM was applied. The obtained standard deviations in both methods are compared (classic levelling vs. three-route levelling), as shown in Figure 5. The three-route levelling technique shows consistent standard deviation values in the three zones, unlike the classic method. A consistent standard deviation implies the reliability of the model’s performance, and thus reliable measurements are less likely to contain systematic errors or biases.

Adjusted normal heights derived by the three-route levelling method and geodetic heights derived from GNSS measurements were used in order to calculate the quasi-geoid height (ζ) for each zone using Equation (12) [25,26,27].

$$\zeta =h-H^N$$

(12)

The values of quasi-geoid heights are compared with the geoid heights derived from EGM2008 in order to estimate its accuracy in the Lebanese territory. By comparing these values, it can be noticed that the EGM2008 can be used in flat areas (zone A) as it shows high similarity in the values. However, in uneven terrains (zones B and C), the EGM2008 shows high discrepancies, which ensures the actual gravitational influence on the height accuracy in those areas (Figure 6). In fact, these discrepancies are due to the absence of gravimetric data in the country, and thus the proposed geometric method is recommended for use in engineering applications that require high accuracy in the vertical positioning [28,29,30]. The model’s accuracy can be achieved in other areas having similar topographic characteristics in the country.

3. Results

According on the varied topography in Lebanon, it is classified into zones based on elevation ranges starting from sea level and ending with the highest elevation of 3088 m above sea level. This classification enables us to find a convenient coefficient for each zone that can be used independently for converting geodetic heights into normal heights based on the developed method for determining normal heights which is the three-route levelling technique.

In general, coefficients serve to scale, adjust, or modify the magnitude, direction, or behavior of the associated variable or term within the context of the equation or model in which they appear. A coefficient in the context of quasi-geoid height analysis represents a numerical value that adjusts the relationship between the observed quasi-geoid heights and the underlying geophysical or topographic factors influencing them. It acts as a multiplier or a scaling factor that modifies the magnitude or shape of the anomaly distribution curve, thereby enhancing the accuracy of the model. Through careful calibration and optimization, coefficients enable researchers to develop robust models that capture the underlying dynamics of the Earth’s surface with greater accuracy [31,32,33].

The results of the determined quasi-geoid heights in the three selected zones in Lebanon are generalized to cover the surroundings of each zone taking into consideration the proximity of elevation ranges in these surroundings. A coefficient was calculated by averaging the absolute difference between quasi-geoid heights and geoid heights. Considering this coefficient in each zone, an enhancement of the derived normal heights from geodetic heights will occur.

Geoid heights derived from EGM2008 at each point in the three zones, and the calculated quasi-geoid heights from normal and geodetic heights, are shown in

Table 3. Differences between quasi-geoid heights and geoid heights in each zone.

Point	ζ (m)			N (EGM2008, m)
	Zone A	Zone B	Zone C	Zone A	Zone B	Zone C
CP	24.827	23.136	25.861	25.245	23.189	25.248
A1	24.836	23.529	25.769	25.254	23.260	25.271
A2	24.836	23.174	25.852	25.253	23.256	25.319
A3	24.816	23.092	25.831	25.234	23.142	25.341
A4	24.819	23.079	25.703	25.239	23.140	25.231

.

The average of the absolute difference was determined in each zone and considered as an additive coefficient for the surrounding area to be used by the GNSS user community. Taking into consideration the comparison between geoid heights and quasi-geoid heights, negative and positive signs were assigned to the coefficient values in order to bring the geoid heights as close as possible to the quasi-geoid heights (

Table 4. Determined coefficient for each zone.

	Zone (A)	Zone (B)	Zone (C)
Coefficient	−0.418	−0.103	0.521

).

To verify the improved accuracy of the geoid heights using the obtained coefficients, a total of twelve checkpoints with known normal and geodetic heights are selected within the surrounding areas of the three zones. These checkpoints are carefully chosen to represent a range of elevations and terrain characteristics in each site. The measured geodetic coordinates of the checkpoints (φ, λ, h) are used to extract the corresponding geoid height (N) value from the EGM2008 model for each point [34,35]. The coefficient is then added to the (N) value, by which the normal heights can be determined (Equation (12)).

Table 5. The measured geodetic and normal heights (H^N) of the checkpoints and their derived geoid heights (N) and orthometric heights (H) from EGM2008.

Checkpoint	φ	λ	h (m)	HN (m)	N (m)	H (m)
1	35.66919	33.51608	883.009	857.965	25.66	857.34
2	35.66930	33.51757	884.935	859.637	25.66	859.27
3	35.67127	33.51815	879.181	853.828	25.65	853.53
4	35.67133	33.51651	879.411	854.252	25.66	853.76
5	35.47342	33.72284	265.147	238.731	26.68	238.47
6	35.46775	33.71645	201.866	175.286	26.73	175.14
7	35.45836	33.72328	80.887	54.162	26.81	54.08
8	35.47109	33.72965	222.225	195.670	26.69	195.53
9	36.13131	34.6029	236.914	213.032	23.44	213.48
10	36.12046	34.607	141.587	117.808	23.42	118.17
11	36.12918	34.61443	208.496	184.433	23.42	185.08
12	36.14072	34.61285	206.589	182.840	23.44	183.15

shows the measured geodetic and normal heights, the derived geoid height (N) from the EGM2008 model, and its corresponding orthometric height. This comprehensive approach ensures the reliability and applicability of the adjusted geoid heights across diverse elevation ranges and terrain types. By incorporating a range of elevations and terrain characteristics in the selection of checkpoints, the effectiveness of the coefficients in enhancing height determination accuracy can be thoroughly evaluated.

After applying the corresponding coefficients to the geoid heights for each checkpoint, the resulting (N) values are subtracted from the geodetic heights in order to calculate new orthometric heights. Figure 7 shows the comparison between the checkpoints’ known normal heights and the orthometric heights before and after adding the coefficient. D1 corresponds to the differences between heights before adding the coefficients, while D2 corresponds to the differences after adding the coefficients. It is noted that the obtained orthometric heights by the help of the coefficients serve a higher accuracy through being closer to the zero level when compared to the actual normal heights which ensures the validity of the proposed method.

4. Conclusions

This study highlights the importance of considering gravitational influences when determining precise heights for engineering applications and many other disciplines. The presence of a local geoid/quasi-geoid model is crucial in accurately accounting for the irregularities in the Earth’s gravitational field, which can significantly impact height calculations. This study’s findings revealed that the uneven distribution of mass within the Earth is affecting the determination of normal heights in Lebanon. This necessitates the creation of a local quasi-geoid model in the country to establish a reliable reference surface for height measurements. By calculating height anomaly using the GNSS and three-route levelling measurements, and developing coefficients for each zone, a more refined quasi-geoid model can be created, enhancing the accuracy of height determinations, as well as the accuracy of geoid heights derived from EGM2008, which was found unsuitable to be used independently for all the regions of the country, particularly for areas with uneven terrain. For applications requiring high accuracy in height determination, it is absolutely recommended to use gravimetric data; otherwise, the proposed combined geometric approach of GNSS and the three-route levelling technique can be used to enhance the accuracy and reliability of height determinations in Lebanon and similar countries, ultimately improving construction, infrastructure development, mine surveys, and urban planning in the country in the absence of a local geoid model.