*Article* **A New Approach for Improving GNSS Geodetic Position by Reducing Residual Tropospheric Error (RTE) Based on Surface Meteorological Data**

**Mario Bakota 1, Serdjo Kos 2, Zoran Mrak <sup>2</sup> and David Brˇci´c 2,\***


**\*** Correspondence: david.brcic@pfri.uniri.hr

**Abstract:** Positioning error components related to tropospheric and ionospheric delays are caused by the atmosphere in positioning determined by global navigation satellite systems (GNSS). Depending on the user's requirements, the position error caused by tropospheric influences, which is commonly referred to as zenith tropospheric delay (ZTD), must be estimated during position determination or determined later by external tropospheric corrections. In this study, a new approach was adopted based on the reduction of residual tropospheric error (RTE), i.e., the unmodeled part of the tropospheric error that remains included in the total geodetic position error, along with other unmodeled systematic and random errors. The study was performed based on Global Navigation Satellite System (GLONASS) positioning solutions and accompanying meteorological parameters in a defined and harmonized temporal-spatial frame of three locations in the Republic of Croatia. A multidisciplinary approach-based analysis from a navigational science aspect was applied. The residual amount of satellite positioning signal tropospheric delay was quantitatively reduced by employing statistical analysis methods. The result of statistical regression is a model which correlates surface meteorological parameters with RTE. Considering the input data, the model has a regional character, and it is based on the Saastamoinen model of zenith tropospheric delay. The verification results show that the model reduces the RTE and thus increases the geodetic accuracy of the observed GNSS stations (with horizontal components of position accuracy of up to 3.8% and vertical components of position of up to 4.37%, respectively). To obtain these results, the Root Mean Square Error (RMSE) was used as the fundamental parameter for position accuracy evaluation. Although developed based on GLONASS data, the proposed model also shows a considerable degree of success in the verification of geodetic positions based on Global Positioning System (GPS). The purpose of the research, and one of its scientific contributions, is that the proposed method can be used to quantitatively monitor the dynamics of changes in deviations of X, Y, and Z coordinate values along coordinate axes. The results show that there is a distinct interdependence of the dynamics of Y and Z coordinate changes (with almost mirror symmetry), which has not been investigated and published so far. The resultant position of the coordinates is created by deviations of the coordinates along the *Y* and *Z* axes—in the vertical plane of space, the deviations of the coordinate X (horizontal plane) are mostly uniform and independent of deviations along the *Y* and *Z* axes. The proposed model shows the realized state of the statistical position equilibrium of the selected GNSS stations which were observed using RTE values. Although of regional character, the model is suitable for application in larger areas with similar climatological profiles and for users who do not require a maximum level of geodetic accuracy achieved by using Satellite-Based Augmentation Systems (SBAS) or other more advanced, time-consuming, and equipment-consuming positioning techniques.

**Keywords:** GNSS; tropospheric error; surface meteorological data; statistical position equilibrium; Saastamoinen model of zenith tropospheric delay

**Citation:** Bakota, M.; Kos, S.; Mrak, Z.; Brˇci´c, D. A New Approach for Improving GNSS Geodetic Position by Reducing Residual Tropospheric Error (RTE) Based on Surface Meteorological Data. *Remote Sens.* **2023**, *15*, 162. https://doi.org/ 10.3390/rs15010162

Academic Editor: Yunbin Yuan

Received: 1 December 2022 Revised: 16 December 2022 Accepted: 23 December 2022 Published: 27 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction and Background**

A total satellite positioning error budget can be decomposed in satellite, receiver, and propagation medium components. The latter refers to atmospheric layers that affect the path of satellite navigation signals, namely the ionosphere and the troposphere. In addition, other present errors with a measurable and significant impact on GNSS positioning include signal multipath errors, receiver noise error, satellite and receiver clock errors, and satellite orbital error [1–3].

Tropospheric error is caused by the propagation of a radio navigation signal through the lowest layer of the atmosphere. In general, the troposphere is divided into two layers: the wet (non-hydrostatic) layer (with a height of up to about 10 km above the Earth's surface) and the dry (hydrostatic) layer (which is 10–40 km above the Earth's surface) [4,5]. These layers cause delays in satellite navigation signals. The troposphere is a non-dispersive medium; therefore, the magnitude of this error component does not depend on signal frequency and cannot be determined as is the case with ionospheric delay. Tropospheric error (which is the reduction of radio navigation signal propagation speed and its deviation from the geometric path and is also commonly called tropospheric delay) can be modeled based on the fundamental meteorological parameters of the troposphere: temperature, humidity, and atmospheric pressure [6–10]. Various models have been developed to predict and reduce tropospheric delay with different scopes of application, including the following: the Two-Quartic Hopfield model (n/a) ((n/a)—not specified or valid for any elevation angle.) [11], the Saastamoinen model (10◦ and above) [12], the Modified Hopfield model (n/a) [13], the Marini model (10◦ and above) [14], the Davis et al. model (C*f*A) (5◦ and above) [15], the Ifadis model (2◦ and above) [16], and the Askne and Nordius model (n/a) [17], etc. The above models allow the value of the tropospheric delay to be estimated with varying degrees of accuracy depending on the input components of the model (including dry, wet, or total delay and with or without mapping function); therefore, they have different elevation angles of mapping function (given in the brackets). Tropospheric delay caused by the wet component is a consequence of the presence of water vapor in all its forms in the upper layer of the troposphere (up to 10 km). Tropospheric delay is partly caused by non-hydrostatic causes; as such, the zenith wet delay (ZWD) is much smaller in absolute terms (several millimeters) than the zenith hydrostatic delay (ZHD—which reaches several tens of centimeters) [4]. The ZWD cannot be accurately modeled using surface meteorological data due to the extreme space and time volatile features of water vapor. With existing tropospheric models, the error caused by ZWD can be reduced by up to 10–20% of its actual value [18]. The causes of the wet component are not in hydrostatic equilibrium; therefore, models based on the partial pressure of water vapor or relative surface humidity do not provide sufficient accuracy. They require empirical constants that vary spatially and temporally [19], although approaches and models for ZWD estimation are being developed based on surface meteorological parameters [20].

The causes of the hydrostatic component of the zenith delay (or atmospheric dry gasses and the non-dipole component of water vapor refraction [18]) are in hydrostatic equilibrium [19] and are therefore determined relatively simply and precisely using the Saastamoinen model. Both components form the ZTD. Approaches of existing tropospheric delay models for ZHD and ZWD differ, and certain models [21–26] of ZTD usually include a non-hydrostatic zenith component which depends primarily on the temporal-spatial distribution of water vapor and the height of the distribution, with the influence of water vapor being the most important. Considering that the insufficient modeling capabilities of the non-hydrostatic causes of tropospheric delay in the zenith direction limit the accuracy of the mapping function for other signal elevation angles, the separation of the mapping function from the zenith delay allowed the development of a number of new models of mapping functions [27–32]. Such models can combine hydrostatic and non-hydrostatic causes and can be combined with tropospheric delay in the zenith direction. This results in hybrid models separated by hydrostatic and non-hydrostatic causes. In general, tropospheric errors in GNSS measurement range from 2–2.4 m in the zenith direction and 25 m at horizontal elevation angles [33,34].

The new approach used in this study investigates the overall effect of tropospheric delay on the accuracy of the GNSS geodetic position of the selected area. Several studies [10,35] show the relationship between tropospheric delays based on radiosonde signal measurements and deviations from position accuracy where the main effect was found to be the atmospheric refraction expressed by the number of N units, which is a value that varies greatly in time and space [36–38]. In this way, it is possible to combine the influence of the tropospheric delays in radio signals (usually with peak values up to several tens of millimeters) with the amounts of slant tropospheric errors in the GNSS system which can reach several tens of meters.

A model of empirical character is proposed which combines the value of the nonmodeled geodetic position RTE GNSS with the meteorological surface parameters of the observed positions. Rather than determining the appropriateness of a particular model of tropospheric delay (either in the zenith or slant direction), the goal of the presented research was to determine the existence of structural dynamics of deviations of x, y, and z coordinates based on the observed GNSS stations in the function of reducing the nonmodeled tropospheric error based on relevant and real surface meteorological data.

The basis for determining tropospheric errors within GNSS position accuracy errors was the Saastamoinen model. The proposed model, together with the previously modeled part of the tropospheric error, reduces the tropospheric error of the GNSS position by correlating the unmodeled portion of the tropospheric error with real meteorological parameters, thereby increasing the overall geodetic accuracy of the determined position.

The following section presents the methodology used in the development of the proposed model, including data collection, statistical regression, and the proposed model validation. The obtained verification results for each location and time period are presented in the third section. The performance of the proposed model and the periodicity characteristics of RTE are presented and discussed in the fourth section, including its potential suitability with GPS. Concluding remarks and possible directions for further research are given in the final section.

#### **2. Methodology**

The development of the presented model is based on combined GLONASS positioning solutions and meteorological data from GNSS reference stations in the mid-latitudes of the Republic of Croatia. Positioning solutions were calculated based on position data in the Receiver Independent Exchange (RINEX) 2.0 compressed format [39]. Initial theories regarding the development of the model included the assumption that the propagation medium error, the user segment, and the microenvironment errors have a constant time– space character. Initial limitations of the proposed model include:


#### *2.1. Time Frame of the Study*

The time frame used for the creation and validation of the proposed model covers the year 2019. The model verification was performed using data from 2014 and 2015. Data from the same locations were analyzed since the model verification required the mutual compatibility of meteorological and GLONASS data at all stages of the study.

#### *2.2. Meteorological Data Collection*

The set of meteorological input data was determined based on the relationship between the propagation of radio signals through the neutral atmosphere and tropospheric dynamic processes. Hydrostatic and non-hydrostatic causes of tropospheric delay were identified. The independent input variables for the model were: pressure *P* (hPa), temperature *T* ( ◦C), precipitation water *PWV* (mm), precipitation *Pr* (mm), and relative humidity *Rh* (percent). Meteorological data can be interpolated from existing numerical weather models (NWM), which is acceptable for analyzing existing tropospheric models, but not for developing a new model according to the selected regional model development. Therefore, in this study, meteorological data collected using automatic meteorological instruments at selected GNSS locations were used by the State Hydrometeorological Institute of the Republic of Croatia (DHMZ) [40] as a source of meteorological data. The time resolution of the meteorological data was 10 min.

Given the current limitations of available data sources, the study was based on the wellknown approach of determining the tropospheric error by determining pseudoranges [10]. The used data delimit the area from 42.60◦ to 46.38◦ North latitude and from 15.22◦ to 18.11◦ East longitude with an altitude range from 64.3 m to 457.9 m above sea level (Figure 1).

**Figure 1.** The wide geographic area of selected GNSS measuring stations.

The climatological profiles of the observed locations differ. According to the Köppen climate classification, which is determined based on the average annual course of air temperature and precipitation, Cakovec is classified as ˇ *Cb* (which is a moderately warm, rainy climate with an average monthly temperature of the coldest month between −3 ◦C and 18 ◦C) [41]. At the same time, the warmest month of the year has an average temperature of less than 22 ◦C [40]. There are no particularly dry months in the year, and it is the area with the least precipitation in the cold season. Cakovec has a humid climate ˇ according to Thornthwaite's climate classification which is based on the ratio of the amount of water needed for potential evapotranspiration and expressed by the humidity coefficient IP/E [41]. Zadar and Dubrovnik belong to areas with a temperate climate with long and hot summers (*Ca*-mark according to the Köppen climate classification), while Zadar belongs to an area with a subhumid (semi-humid) climate in terms of the humidity coefficient IP/E [42]. Dubrovnik has a humid climate due to heavier precipitation. Water vapor and precipitation amounts are the main meteorological input parameters for determining the non-hydrostatic zenith component of the GNSS position. The above measurement locations do not have significantly pronounced differences in their humidity coefficients; as such, they represent a suitable choice for creating the proposed model, which is also applicable to larger geographical areas with similar climate profiles.

#### *2.3. Geographic GNSS Data Collection*

The study was carried out based on data from locations equipped with GNSS stations that provided predicted meteorological parameters with adequate time resolution. Therefore, measuring stations in the Regional Reference Frame Sub-Commission for Europe (EUREF) in Cakovec, Zadar, and Dubrovnik were selected. The basic stations' data are ˇ listed in Table 1 [43].



<sup>1</sup> Earth centered, earth fixed coordinate system (The European Terrestrial Reference System, 1989).

#### *2.4. Model Development*

The development of the model involves the harmonization of two parallel input components: the value of the geodetic deviation of the user's position caused by the tropospheric error and the temporal-spatial harmonization of the corresponding meteorological input data.

#### 2.4.1. Determination of the Size of the Tropospheric Error

Basically, the determination of the tropospheric error by measuring the pseudorange can be conducted in two ways:

(i) by determining and removing all systematic and random errors; and

(ii) by determining the deviation from the known position using a preselected tropospheric model and estimating the deviation residuals as unknown parameters [25].

In the first method, the pseudorange is calculated according to the general formula [10]:

$$R\_A^l = \rho\_A^l + \delta\rho\_{mul} + \delta\rho\_{rel} + c\delta\_A - c\delta^l + I\_A^l + T\_A^l + e\_A^l \tag{1}$$

where *R<sup>i</sup> <sup>A</sup>* is the pseudorange from position *<sup>A</sup>* to satellite *<sup>i</sup>*; *<sup>ρ</sup><sup>i</sup> <sup>A</sup>* is the geometric distance; *c* is l speed of light; *I<sup>i</sup> <sup>A</sup>*, *<sup>T</sup><sup>i</sup> <sup>A</sup>* is ionospheric and tropospheric delay; *<sup>c</sup>δA*, *<sup>c</sup>δ<sup>i</sup>* is the satellite orbit and clock error; *δρmul* is the multipath trajectories error; *δρrel* is relativistic error; and *e<sup>i</sup> <sup>A</sup>* is the random error.

The multipath error *δρmul* can be determined programmatically, usually based on the receiver and antenna equipment manufacturer. Sources of ionospheric delay, relativistic errors, satellite orbit errors, and clock errors are included in the navigation message or, in the case of post-processing, from the appropriate ground truth data. In any case, it is necessary to perform software processing to ensure appropriate data sources are used for additional corrections. For additional verification, the values of the isolated and determined tropospheric delays of a known position were compared with the values determined using the radiosonde signals or another system [10].

The model development was based on the determination of the user's position in accordance with the selected modeling approach and its initial limitations. It was obtained by determining the pseudoranges and isolating the accuracy deviation from the user's position caused by the tropospheric error within the total positioning error. A position determined in this way can be defined in a simpler form as the difference between the signal reception time which is determined based on its clock and time (determined by the satellite clock), which can be represented with the following expression [4,44]:

$$P\_{r,i}^S = \mathcal{c}\left(\overline{\mathfrak{f}}\_r - \overline{\mathfrak{f}}^s\right) \tag{2}$$

where *P<sup>S</sup> <sup>r</sup>*,*<sup>i</sup>* is the pseudorange of the i-th satellite; *tr* is the time of the signal reception determined by the receiver clock (s); and *t <sup>s</sup>* is the time of signal transmission determined by the satellite clock (s).

The misalignment of the satellite and receiver clocks (with input data contained in the navigation message), the ionospheric, tropospheric, and measurement error, and the expression (2) takes the form (3) by introducing the parameters of the geometric distance between the antennas of the satellite and the receiver:

$$\begin{split}P\_{r,i}^{\mathbb{S}} &= c\left(\left(t\_{r} + dt\_{r}(t\_{r})\right) - \left(t^{\mathbb{S}} + dT^{\mathbb{S}}\left(t^{\mathbb{S}}\right)\right)\right) + \varepsilon\_{P} \\ &= c\left(t\_{r} - t^{\mathbb{S}}\right) + c\left(dt\_{r}(t\_{r}) - dT^{\mathbb{S}}\left(t^{\mathbb{S}}\right)\right) + \varepsilon\_{P} \\ &= \left(\rho\_{r}^{\mathbb{S}} + I\_{r,i}^{\mathbb{S}} + T\_{r}^{\mathbb{S}} + c\left(dt\_{r}(t\_{r}) - dT^{\mathbb{S}}\left(t^{\mathbb{S}}\right)\right) + \varepsilon\_{P}\right) \\ &= \rho\_{r}^{\mathbb{S}} + c\left(dt\_{r}(t\_{r}) - dT^{\mathbb{S}}\left(t^{\mathbb{S}}\right)\right) + I\_{r,i}^{\mathbb{S}} + T\_{r}^{\mathbb{S}} + \varepsilon\_{P} \end{split} \tag{3}$$

where *I<sup>S</sup> <sup>r</sup>*,*<sup>i</sup>* is the geometric distance between the satellite and receiver antenna; *dtr*, *dT<sup>S</sup>* is the receiver and satellite clock offset; *I<sup>S</sup> <sup>r</sup>*,*<sup>i</sup>* is the ionospheric error; *<sup>T</sup><sup>S</sup> <sup>r</sup>* is the tropospheric error; and *ε<sup>P</sup>* is the measurement error.

The conversion of the geodetic position into a data set in ECEF ETRS89 format with values expressed in meters is expressed according to the following expressions:

$$el^2 = f(2-f)\tag{4}$$

$$\nu = \frac{a}{\sqrt{1 - e l^2 \sin \phi\_r^2}}\tag{5}$$

$$r\_r = \begin{pmatrix} (\nu + h)\cos\phi\_l \cos\lambda\_r \\ (\nu + h)\cos\phi\_l \sin\lambda\_r \\ \nu \left(1 - \varepsilon l^2\right) \sin\phi\_f \end{pmatrix} \tag{6}$$

where *λ<sup>r</sup>* is the geodetic longitude; *φ<sup>r</sup>* is the geodetic latitude; *h* is the ellipsoidal height; *a* is the semi-major axis of the reference Earth ellipsoid (6,378,137 m); *el*<sup>2</sup> is the first numerical eccentricity of the ellipsoid; and *f* is the flattening coefficient of the reference Earth ellipsoid.

The Saastamoinen model of tropospheric zenith correction was used, and the present parameters were transmitted in the navigation message as a common source of input program corrections for all measuring stations and periods and for standardizing other program settings. The input program parameters for pressure, absolute temperature, and partial pressure are determined by the expressions of the standard atmosphere model [44]:

$$P = 1013.25 \cdot \left(1 - 2.2557 \cdot 10^{-5} h\right)^{5.2568} \tag{7}$$

$$T = 15 - 6.5 \cdot 10^{-3} h + 273.15\tag{8}$$

$$\sigma = 6.108 \cdot \exp\left(\frac{17.15 \ T - 4684}{T - 38.45}\right) \cdot \frac{h\_{rel}}{100} \tag{9}$$

where *P* is the total air pressure (hPa); *T* is the temperature in degrees Kelvin; *e* is the partial atmospheric pressure (hPa); *h* is the geodetic altitude above sea level; and *hrel* is the relative humidity. The applied Saastamoinen model uses a constant relative humidity value of 70%.

The tropospheric correction in this configuration was calculated using the following algorithms [45,46]. For the mapping function in the selected program setting, the Niell model is calculated according to the expression [47–49]:

$$m(\varepsilon) = m\_{\text{\textquotedblleft}}(\varepsilon) \{ 1 + \cot \varepsilon \left( G\_N \cos z + G\_E \sin z \right) \}\tag{10}$$

where *ε* is the signal elevation angle; *z* is the signal zenith angle; *GN* is the tropospheric gradient in the north direction; *GE* is the tropospheric gradient in the east direction; and *mw* is the non-hydrostatic mapping function of the non-hydrostatic Niell (New) Mapping Functions (NMF).

The total tropospheric delay was calculated according to the following expression:

$$d\_{tro} = m\_d(\varepsilon)d\_d^z + m(\varepsilon)(d\_{tot}^z - d\_d^z) \tag{11}$$

where *dtro* is the total tropospheric delay; *d<sup>z</sup> <sup>d</sup>* is the hydrostatic component of zenith tropospheric delay (in meters and determined by the Saastamoinen model); *d<sup>z</sup> tot* is the tropospheric total zenithal retardation; and *md* is the (NMF) hydrostatic mapping function. The parameters of tropospheric gradients and total tropospheric delay in zenith direction were estimated using the extended Kalman filter (EKF) [45]. The parameters for the correction of the ionospheric delay (A—the numerical coefficient of the maximum total free electron content in the ionospheric layer F2; F—the index of solar activity; and Ap—the daily index of geomagnetic activity) included in the GLONASS navigation message (broadcast ionosphere model) [47] have the following form:

$$P\_{ion} = \left(\mathfrak{a}\_{0\prime}\mathfrak{a}\_{1\prime}\mathfrak{a}\_{2\prime}\mathfrak{a}\_{3\prime}\mathfrak{a}\_{0\prime}\mathfrak{b}\_{1\prime}\mathfrak{b}\_{2\prime}\mathfrak{b}\_{3}\right)^{T} \tag{12}$$

The determination of ionospheric delay *I<sup>S</sup> <sup>r</sup>* ˆ (m) was performed using Klobuchar's model due to single-frequency computing. Despite the efficiency of Klobuchar's model, more effective reducing of the ionospheric delay would be achieved using an iono-free combination or dual-frequency receiver. However, considering the format of available input data, a single-frequency receiver with the Klobuchar model was used. Input data regarding clock parameters and ephemerides were sent in the navigation message (broadcast ephemerides and clock parameters) and used in calculations in the following form [44,46,50]:

$$P\_{\rm cph}(t\_b) = \begin{pmatrix} \mathbf{x}\_{\prime} y\_{\prime} z\_{\prime} \boldsymbol{\nu}\_{\mathbf{x}\prime} \ \boldsymbol{\nu}\_{y\prime} \boldsymbol{\nu}\_{z\prime} a\_{y\prime} a\_{y\prime} a\_{z\prime} \boldsymbol{\pi}\_{n\prime} \boldsymbol{\gamma}\_{n\prime} \end{pmatrix} \tag{13}$$

In addition, the initial software setup included a single positioning mode and a 3◦ elevation mask value (due to the scope of applicability of the NMF mapping function) to isolate the deviation in the geodetic accuracy of the user's position caused by the tropospheric error component.

### 2.4.2. Statistical Analysis and Model Development

The remaining value of the tropospheric error was determined based on the difference in the geodetic accuracy of the position determined both with and without tropospheric correction (using the Saastamoinen model). It is important to emphasize that, in addition to the tropospheric component, such determined deviations still contain several unmodeled systematic and random errors, including the residual ionospheric error, the satellite position and clock error, errors related to multipath, and solid tides. However, the only input difference when comparing the accuracy of the geodetic position was the application of the Saastamoinen model of tropospheric correction (as the other applied algorithms were identical); therefore, the effect of the resulting final difference between the two models can be deterministically attributed to the tropospheric component within the total geodetic position error.

The difference in geodetic accuracy of the user's position can be analyzed as a function of the influence of the applied Saastamoinen model of tropospheric delay correction since all other parameters were set identically in both cases. This provided the theoretical basis for quantifying the effect of the Saastamoinen model on improving the geodetic accuracy of the observed GNSS positions.

The resulting RTE value of the known position caused by the non-modeled part of the tropospheric delay was subjected to a statistical regression procedure to determine the correlation with real meteorological surface parameters. Statistical correlation with its positive parameters are the basis for developing a mathematical model to reduce RTE as a function of surface meteorological parameters (as independent input variables).

The stages in the creation of the proposed model are shown in Figure 2.

**Figure 2.** The proposed model development.

Meteorological independent input variables are as follows: *T* is the temperature (◦C); *Pr* is the precipitation (mm); *Rh* is the relative humidity (percentage); *PWV* is the precipitation water (mm); *P* is the pressure (hPa); and Δ*XM*, Δ*YM*, and Δ*ZM* are deviations along the *x*, *y*, and *z* axes in the proposed model (m).

The proposed model contains a mathematical expression for each axis as the observed deviation from the exact geodetic position followed the ECEF coordinate system. The final form of the proposed model is shown as follows. The amounts of the coefficients of the input predictors are the result of the regression analysis that describes the form and intensity of the mutual connection between the meteorological input predictors and the output variable RTE. For the *X* axis, the model component Δ*XM* has the form:

$$
\Delta X\_M = \sqrt{-4.66541 - 0.01108Rh + 0.07062P + 0.11806Pr} - 0.05187T + 0.1246PNV \tag{14}
$$

For the *Y* axis, the model component Δ*YM* has the form:

$$
\Delta Y\_M = -\left(\sqrt{0.06656 - 0.00004Rh - 0.00003P + 0.00041Pr + 0.000047T - 0.00345PNV}\right) \tag{15}
$$

For the *Z* axis, the model component Δ*ZM* has the form:

$$
\Delta Z\_M = -\frac{1}{7.17152 - 0.01049R\_h - 0.00407P + 0.12603P\_r - 0.03837T + 0.01937PWh} \tag{16}
$$

The final form of the proposed model (*PSi*) is based on an extension of the existing Saastamoinen model and represents the sum of the corrections made by the Saastamoinen model and the proposed model for each axis:

$$P\_{Si} = M\_S + \begin{cases} \Delta X\_M\\ \Delta Y\_M\\ \Delta Z\_M \end{cases} \tag{17}$$

where *MS* is the correction value realized by the Saastamoinen model (in m).

#### 2.4.3. Validation of the Proposed Model

Validation was performed using the cross-validation method with a 50:50 train-test ratio. Common statistical indicators of significance of the model components for a single coordinate axis include p value, multiple correlation coefficients, multiple determination, and adjusted determination coefficient. The results are shown in Table 2.


**Table 2.** Statistical indicators of the proposed model validation.

The obtained results show statistically significant correlations (upper limit *p* = 0.05), i.e., there is a statistically significant correlation between the independent input variables (the meteorological parameters) and the dependent output variable (RTE) for each coordinate axis, expressed by the Equations (14)–(16). The obtained correlation coefficients show: (i) the presence of a statistically significant correlation between the input predictors and the dependent variable in all components of the proposed model; (ii) a statistically significant prediction of the criterion variable by the input predictors in all components of the proposed model; and (iii) a statistically significant relationship between the obtained correlation, the number of samples (the input records), and the number of input variables (the predictors) of all components of the proposed model. Although the existence of a statistically significant connection between the input predictors and the output RTE variable is shown, the realized adjusted determination coefficients show different connection intensities, which is evident considering the associated low values of the coefficients of the input predictors in the model component for the *y* axis.

#### **3. Results and Findings**

The verification of the proposed model was performed using independent input data and meteorological data from the same locations for the period of 2014–2015. The coordinates of the geodetic positions for the period of 2014–2015 were calculated in an identical way and with the same input algorithms as the geodetic positions for 2019. The proposed model was created on this basis. The fundamental parameter that highlights the success of the proposed model is the reduction of the RMSE deviation of the positional accuracy compared to the positional accuracy achieved exclusively by the Saastamoinen model. The verification process was performed separately for each location and period. The movement of the coordinates along the coordinate axes in the horizontal and vertical planes and the inversion of the deviation of individual coordinates were observed using the proposed validated and verified empirical model; for example, y and z are complexly and graphically shown on a series of diagrams provided in Appendix A of the paper.
