*3.2. Verification for Zadar Location*

The correction values generated by the proposed model are shown in Table 5.


**Table 5.** Values of RTE corrections based on the proposed model for Zadar GNSS location (in m).

It can also be seen that the values of generated RTE corrections for the GNSS station in Zadar were negative for the y and z axes in the observed period and ranged from −0.85 m to 0.9 m for the *x* axis. The absolute corrections' values ranged from 0.055 m (*y*-axis in 2014) to 1.74 m (*x*-axis in 2015). The statistical parameters RMSE, STD, range, and median for the Zadar GNNSS position are shown in Table 6.

**Table 6.** Comparison of the achieved values of RMSE, STD, minimum, maximum, and median values for the Zadar GNSS position obtained using the Saastamoinen model and the proposed model (in m).


Comparison of the parameters showed that the proposed model did not have a significant effect on the reduction of the final values of the accuracy deviations, and the range of the minimum and maximum values remained approximately the same. However, there was a significant change in the grouping values of the measurement frequencies for 2014 in all coordinate axes. For 2015, a positive shift of the median was obtained for the x and z axes, and an additional deviation from the zero value was registered for the *y* axis (0.26 m versus 0.10 m). The 2014 results showed an improvement in the accuracy of the RMSE value of the position at the same time for all axes. Accuracy improvements were obtained for the *x*-axis (3.33 cm/0.78%), the *y*-axis (0.9 cm/0.37%), and the *z*-axis (22.71 cm/4.33%). The parameter STD improved for the *x*-axis, and the value STD decreased for the *y*-axis (−0.03 cm/−0.01%) and the *z*-axis (−0.006 cm/−0.001%). For 2015, the results of the predictive model showed the RMSE values of the *x*-axis (11.21 cm/2.6%) and the *z*-axis (1.45 cm/0.35%), while the *y*-axis had a decrease in the RMSE value (−1.65 cm/−0.65%). The same tendency was seen in the parameter STD: the *x*-axis and *z*-axis obtained improvements, and a deterioration of the *y*-axis STD value (−0.03cm/−0.01%) was recorded.

#### *3.3. Verification for Dubrovnik Location*

The correction values generated by the proposed model for the Dubrovnik GNSS location are shown in Table 7.

**Table 7.** Values of the RTE corrections based on the proposed model for the Dubrovnik GNSS location (in m).


The values shown indicated that the model generated negative correction amounts of RTE for the *y*- and *z*-axes over the observed two-year period and values within the limits of −0.96 m to 1.69 m for the *x*-axis. The values of the absolute correction amounts ranged from 0.047 m (*y*-axis in 2014) to 2.33 m (*x*-axis in 2014). The values of the RMSE, standard deviation (STD) and associated statistical parameters for Dubrovnik GNSS position are shown in Table 8.

**Table 8.** Comparison of the realized values of RMSE, STD, range, and medians for the Dubrovnik GNNSS position obtained using the Saastamoinen model and the proposed model (in m).


Analysis of the presented accuracy parameters showed that the proposed model contributed to a deterioration of the RMSE value in 2014 for the *x*-axis (−1.35 cm/−0.3%), while it improved the accuracy of the other two axes, especially the *z*-axis (23.87 cm/4.37%), and also the *y*-axis (2.5 cm/0.92%). As for the two previously observed positions, there was no significant change in the value for the deviations range. An improvement in accuracy was observed for the median parameter in 2014 in the *y* and *z* axes, and an increase in the deviation was recorded for the *x* axis (from -0.17 m to 0.54 m). The parameter STD improved in the *x*-axis (1.78 cm/0.39%) and worsened in the *y*-axis (−0.004 cm/−0.001%) and *z*-axis (−0.12 cm/−0.02%).

For 2015, the proposed model improved the accuracy (parameter RMSE) of the *x*-axis (2.2 cm/0.57%), *y*-axis (0.3 cm/0.13%), and *z*-axis (11.15 cm /2.52%). The value of STD also improved for the *x*-axis, while it decreased for the *y*-axis (−0.03 cm/−0.01%) and the *z*-axis (−0.14 cm/−0.03%). At the same time, there were shifts in the median value toward the central value for all axes, with shifts on the *y*-axis (from −1.679 to 0.03 m) and the *z*-axis (from −3.84 to 0.53 m).

#### **4. Discussion**

Summary results of the verification procedure for the observed locations are presented in Table 9. The table contains data obtained by comparing the RMSE and STD parameters of the position accuracy obtained using the Saastamoinen model and the proposed model for all observed locations within the observed time frame.


**Table 9.** Total results of RMSE and STD values achieved using the proposed model.

Positive values represented an increase in RMSE and STD values, i.e., a decrease in position accuracy, while negative values represented a decrease in the aforementioned statistical parameters, i.e., an increase in geodetic position accuracy. The graphical results of the performance of the proposed model are shown in Figures 3 and 4 where the curve represents the relative value (in percentages) and the columns represent the absolute amount (in cm).

**Figure 3.** Movement of RMSE deviation values using the proposed model from the exact geodetic position of the observed locations expressed in absolute (right *y*-axis) and relative (left *y*-axis) values.

**Figure 4.** Movement of the parameter STD using the proposed model of observed locations expressed in absolute (right *y*-axis) and relative (left *y*-axis) values.

The presented results of RMSE as the main performance parameter showed that the proposed model improved the overall geodetic accuracy of the observed positions. The *x*-axis position deviations decreased for five measurements, while the deviations in the *y*-axis and especially in the *z*-axis decreased simultaneously (Dubrovnik, 2014).

The model for the *y*-axis achieved a decrease in the deviation in four measurements, while where there was an increase in the deviation in the remaining measurements (Cakovec, 2014 and Zadar, 2015), and there improvement in the accuracy of the position ˇ along the second horizontal *x*-axis was achieved in addition to the improvement achieved along the vertical *z*-axis.

The results of the obtained STD values as measures of the variability of the obtained results were half-hearted (partial). Improvement was obtained in nine measurements and worse results were obtained in the other nine measurements, as shown in Figure 4. It is important to point out that the absolute amounts of the obtained worse results were within the limits of 0.0014 cm to 0.1448 cm, although they were nominally worse than the STD value obtained using the Saastamoinen model. Considering the mentioned amounts, the obtained worse results can be ignored in further interpretation of the results and evaluation of the study's success.

#### *4.1. Model Suitability for Application within the GPS*

The proposed model shows a certain degree of success in the application within the GPS, although it was developed on the basis of GLONASS position records. The results showing the movement of the RMS parameter are shown in Figure 5.

The proposed model generally increases geodetic accuracy throughout the verification period. The RMS value of the *x*-axis deviation for all cities in the entire verification period reduced with absolute reduction amounts ranging from 3.29 cm to 32.26 cm (1.32–11.08%). At the same time, the geodetic accuracy decreased along the *y*-axis for the absolute values from 1.64 cm to 4.77 cm (1.09% to 5.54%). In the *z* axis, the model showed a variable result: it achieved an improvement in four of the six measurements (in absolute amounts from 0.65 cm to 2.49 cm, i.e., 0.34% to 1.23%). In the remaining two measurements, the proposed model degraded geodetic accuracy (varying from 3.55 cm to 16.32 cm, i.e., 1.78% to 7.95%). Differences in the model's ability to reduce the residual tropospheric delay

of the satellite signal when compared to GPS are due to a number of reasons, including the different spatial constellation of systems, differences in modulations of the satellite navigation signals, etc. For adequate results within the GPS, a regression analysis should be performed as a function of the interdependence of meteorological input predictors and RTE output variables based on GPS positional data, which further observes and analyzes the dynamics and variations of RTE movements.

**Figure 5.** Movement of the RMSE deviation value using the proposed model from the exact geodetic position of the observed locations obtained using a GPS system and expressed in absolute (right *y*-axis) and relative (left *y*-axis) values.

As expected, the results do not achieve an improvement as with the GLONASS system. Given the time and space limitations of the research, GLONASS was chosen as an affirmed part of the GNSS system since the main goal of the research was the development and adoption of a methodology based on statistical analysis. However, the verification of the proposed model using GPS data showed the potential of the adopted approach and the possibility of further development based on the input of multi-GNSS positioning data for a wider GNSS application.

#### *4.2. Periodic Effect of the Proposed Model on the Positioning Accuracy*

The proposed model showed the possibility of a quantitative influence on the reduction of the tropospheric error which increases the accuracy of GLONASS positions. The tropospheric error is a stochastic phenomenon due to natural causes; however, the conducted research proves the possibility of developing a statistically significant correlation between tropospheric dynamics (expressed as a set of meteorological parameters) and the improvement of GNSS position accuracy. The regularity of the movement of the absolute differences of the values of RTE according to the proposed and Saastamoinen models can be observed; thus, indirectly, the movement of the geodetic accuracy achieved via the GLONASS system can also be observed. The annual motions shown in Figures 6–8 (see Appendix A for additional Figures A1–A7) show the annual oscillation for each position axis. The value for the *x*-axis is shown in the upper part of the graph, while the middle part of the graph shows the values for the *z*-axis, and the lower part of the graph shows the realized values for the *z*-axis. The property of periodicity of the motion manifests itself at smaller temporal resolutions—for example, at quarterly or monthly time frames—regardless of the selected GNSS position (see Appendix A: GNSS Position Zadar 2015).

**Figure 6.** Movement of the absolute difference of the parameter RTE at the GNSS position Cakovec ˇ in 2014.

**Figure 7.** Movement of the absolute difference of the parameter RTE at the GNSS position Zadar in 2014.

**Figure 8.** Movement of the absolute difference of the parameter RTE at the GNSS position Dubrovnik in 2014.

The annual trends in the difference of shown absolute RTE deviations clearly indicate the seasonal nature of the tropospheric error. On the *x*-axis, the difference in absolute RTE deviations typically reaches a maximum in the winter months, although the upper values of the deviations may also occur in the later periods of spring or autumn, i.e., in summer in the case of climatological deviations for the mentioned seasons.

The same pattern of correlation of the differences in absolute deviations from RTE along the y and z axes is evident in the other two observed GNSS positions. There are differences in the absolute values (amounts) of the RTE deviation difference, but there are no differences in the general distribution of the measured RTE deviation difference.

#### *4.3. Observations on the Specifics of the Conducted Research*

The dynamic specificities of the presented relationship between the positioning accuracy and the absolute values of the RTE difference are shown in the vertical level in the direction of the *z*-axis (regarding the relationship with the zenith angle of the incoming radio navigation signal) and in the horizontal plane in the direction of the *y*-axis. The vertical plane is defined by the *y*- and *z*-axes; therefore, improvement in the positional accuracy is expressed in such a way that the increase in deviation along the *y*-axis is followed by a decrease in deviation along the *z*-axis and vice versa. At the horizontal level, the deviation along the *x*-axis is somewhat proportionally related to the value achieved along the *y*-axis. Simultaneous deviations along the *y*- and *z*-axis show the mutual relationship achieved at the vertical level which they jointly define in the observed space as the starting point.

It is possible to define the RTE vector using the origin in the center O of the coordinate space (*Oxyz*) based on analytical interpretation of the accuracy dynamics of GNSS positions (observed based on the movement of the difference of the absolute values of RTE realized by the Saastamoinen model and the proposed model where the reduction of the parameter RTE improves the geodetic accuracy of the GNSS position). The motion of the RTE vector within the stationary *Oxyz* system describes the translational motion of the RTE parameter/point within the available (six) degrees of freedom. The proposed model shows improvement of the RTE parameter within the observed *Oxyz* system of taking the position (or the reduction of the RTE parameter) independent of time where the specified property of the *Oxyz* system is also independent of time [51]. Therefore, it is justified to conclude that the observed GNSS system has reached a state of statistical position equilibrium. The original assumption was that a decrease in the RTE parameter leads to a decrease in the tropospheric error which increases the geodetic accuracy of the position. The observed GNSS GLONASS system also reaches this state due to the realized movement of the RTE parameter.

#### **5. Conclusions**

The proposed model showed its success by reducing deviations from exact GNSS geodetic positions by acting on the non-modeled part of the tropospheric error. The basis of the proposed model was the existing Saastamoinen model, and correction values (obtained using the proposed model) for each axis should have been added to this model. The proposed model needs meteorological input parameters which can be interpolated from a given standard weather model depending on availability; alternatively, standard predicted input values can be used (as is the case of the absence of meteorological parameters in the Saastamoinen model).

The proposed model was verified at the same GNSS stations within a two-year period. The proposed model showed an improvement in position accuracy achieved by reducing the residual tropospheric error when compared to the Saastamoinen model. The model did not achieve a simultaneous improvement in all axes for all locations during the entire verification period, but it demonstrated superiority over the Saastamoinen model. Improvements in the horizontal axes of the position up to a maximum of 3.87% were achieved (14.26 cm), while the accuracy of the second horizontal axis reduced by 0.65% (1.65 cm) for two measurements. At the same time, the accuracy of the height component of the position improved in all measurements to a maximum of 4.37% (23.87 cm).

The proposed model must be used programmatically (software) as a complement to the Saastamoinen model, although there are certain limitations in terms of the geographical area of application, i.e., the possibility of application in areas with similar climate profiles. At the same time, the optimal application areas of the proposed model are found in stationary and dynamic systems to determine the position of the user in real time with lower accuracy. Therefore, it is clear that the model can be used within the existing application areas of the Saastamoinen model with all existing limitations and advantages.

An additional result of the applied methodology and the use of the RTE parameter is the statistical position equilibrium of the observed GNSS GLONASS positions defined by the oscillation around the central position values within the *Oxyz*-space. The model was developed based on GLONASS position data; however, it also shows a certain level of success in the verification of geodetic positions based on GPS position. The possible continuation of this research includes applying this suggested approach to developing a model based on multi-GNSS positioning data and verifying its effectiveness with other available GNSS systems. Future research of this type will be focused on determining the positional statistical balance of the X, Y, and Z coordinates along the x, y, and z coordinate axes as a function of various satellite navigation arguments.

**Author Contributions:** Conceptualization, S.K. and M.B.; methodology, S.K. and M.B.; software, M.B.; validation, S.K., Z.M. and D.B.; formal analysis, M.B. and Z.M.; investigation, M.B.; resources, D.B.; data curation, M.B.; writing—original draft preparation, S.K. and M.B.; writing—review and editing, D.B.; visualization, M.B.; supervision, S.K. and D.B.; project administration, D.B.; funding acquisition, S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The GNSS observation data are available at http://www.epncb.oma. be/index.php (accessed 29 August 2021). Meteorological data are available at https://meteo.hr/ proizvodi.php?section=katalog\_zahtjevi&param=zahtjev\_podaci\_usluge upon request (accessed 29 August 2021).

**Acknowledgments:** This work was fully supported by the University of Rijeka under the project uniri-tehnic 18-66: *Research of environmental impacts on satellite navigation systems*.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **Appendix A**

**Figure A1.** Movement of the absolute difference of the RTE GNSS Cakovec position parameter ˇ in 2015.

**Figure A2.** Movement of the absolute difference of the RTE GNSS Zadar position parameter in 2015.

**Figure A3.** Movement of the absolute difference of the RTE GNSS Zadar position parameter in 2015.

**Figure A4.** Movement of the absolute difference of the RTE GNSS Zadar position parameter in the period July–September 2015.

**Figure A5.** Movement of the absolute difference of the RTE GNSS position parameter in Zadar, July 2015.

**Figure A6.** Movement of the absolute difference of the RTE GNSS position in Zadar, August 2015.

**Figure A7.** Movement of the absolute difference of the RTE GNSS position parameter in Zadar, September 2015.
