Integration of GNSS Technology and Geometric Levelling for Accurate Height Determination in Uneven Terrains: A Case Study in Lebanon

1. Introduction

Geometric leveling still remains the most reliable method of high-precision determination of heights. This method involves direct measurements of elevations, from which normal heights are determined. Normal height system relates heights to a hypothetical surface called the quasi-geoid that closely approximates mean sea level, while 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].

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 (HN) of these points should still be determined on the basis of geometric leveling [2,3]. The object of the study is the territory of Lebanon, where there is neither 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 geometric leveling method using 3-route levelling technique.

Figure 1. Comparison of different height systems and their base surface.

Figure 1. Comparison of different height systems and their base surface.

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, EIGEN-6C3stat are considered more accurate, the practical results of their use (K.I. Markovich, D.Sh. Fazilova, etc.) show identical results when comparing models, or a difference within 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 characteristic 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, 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 topography. 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.

Four satellite receivers were used in static mode, and carried out for 1 hour/session. Figure 2 shows the distribution of points in each zone. All GNSS measurements started based on a chosen control point from classes II or III in each study area. In (A), the control point 8F3Z from class-III was allocated. A1-A4 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 based on them. The control points are 7MR5, 7JLB, and 7PRP are from class-III. However, in zone (C), the chosen control point was QDF7 from class-II.

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.

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, then the elevation differences are compared. The process of measuring heights through the three different routes (1 straight route and 2 curvilinear routes) occurs between 2 main points, passing through an intermediate point in each followed curvilinear route (Figure 3). For example, the points A1 and CP are measured three times, where back sight and fore sight readings are taken along the 1st route (straight: baseline-B), then along the 2nd route (curvilinear: baselines 2 and 5) passing through the intermediate point IP1, then along the 3rd route (curvilinear: baselines 21 and 24) passing through IP4. 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). 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:

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

(5)

The unknown parameters are calculated by:

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

(6)

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

The equivalent residuals $V_e$ are:

$$V_e=B∆-f,$$

(7)

The observational residuals are:

$$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 }}_{∆∆}={\sigma }_0^2.{{[B}^T.W_e.B]}^{-1},$$

(10)

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

$$\sigma ={\sigma }_0.\sqrt{a_{\mathrm{1,1}}},$$

(11)

where $a_{\mathrm{1,1}}$ is the diagonal element of the covariance matrix.

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.

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

Section	Route	Correction of Systematic errors (mm)			Disclosure (m)	ΔHN	LSM
		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

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

Section	Route	Correction of Systematic errors (mm)			Disclosure (m)	ΔHN	LSM
		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

In order to evaluate the achieved accuracy using the proposed method (3-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 deviation in both methods are compared (classic levelling vs. 3-route levelling) as shown in Figure 4. The 3-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 3-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 8 [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 5). 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 the 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 “the 3-route levelling technique”.

In general, coefficients serve to scale, adjust, or modify the magnitude, direction, or behaviour 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.

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 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.

The average difference was determined in each zone and considered as an additive coefficient for the surrounding area to be used by GNSS user community (Table 4).

To verify the improved accuracy of the geoid heights using the obtained coefficients, a total of 12 check points with known normal and geodetic heights are selected within the surrounding areas of the three zones. These check points are carefully chosen to represent a range of elevations and terrain characteristics in each site. The measured geodetic coordinates of the check points (ϕ, λ, h) are used to extract the corresponding geoid height (N) value from EGM2008 model for each point. The coefficient is then added to the (N) value, by which the normal heights can be determined (equation 12). Table 5 shows the measured geodetic and normal heights, the derived geoid height (N) from EGM2008 model and its corresponding orthometric height.

After applying the corresponding coefficients to the geoid heights for each check point, the resulting (N) values are subtracted from geodetic heights in order to calculate new orthometric heights. Figure 6 shows the comparison between the check points’ 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 being closer to the 0 level when compared to the actual normal heights which ensures the validity of the proposed method.

5. Conclusions

This study highlights the importance of considering the 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. The study 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 3-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 3-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.