Acessibilidade / Reportar erro

A comparative study of Runge-Kutta method orders for computing GLONASS broadcast orbits

Abstract:

The determination of the position and velocity of the GLONASS user’s begins with the GLONASS satellite’s orbit computation. In this paper, using the positions and velocities of GLONASS satellites given in the broadcast ephemerides file every 30 minutes as initial conditions, we study the effect of the order and the step size in the integration of the differential equation of satellite motion by the Runge-Kutta method to get the positions and velocities of GLONASS satellites at any time. This method consists of several orders, including order four recommended in the GLONASS Interface Control Document; order five; and order four and five. These three orders are tested in this study using two distinct step sizes (1 sec and 0.5 sec). In terms of the differences obtained between the forward (+15 min) and the backward (-15 min) integration processes and the runtime, order four is the most suitable for the determination of GLONASS orbits compared to the other orders employed in this study. The data used is the positions, velocities, and luni-solar acceleration of all GLONNAS satellites on April 20, 2021, given in the broadcast ephemerides file.

Keywords:
Broadcast Ephemerides; GLONASS Orbits; Runge-Kutta Order; Step Size; Runtime

1. Introduction

The acronym GLONASS stands for Globalnaya Navigazionnaya Sputnikovaya Sistema, meaning Global Navigation Satellite System. It is the Russian equivalent to the American Global Positioning System (GPS), created by the former Soviet Union and presently operated by the Russian Space Forces (Dach et al. 2015Dach, R., Lutz, S., Walser, P. and Fridez, P. 2015. Bernese GNSS Software Version 5.2. User manual. Astronomical Institute, University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297
https://doi.org/https://doi.org/10.7892/...
). The first GLONASS satellite was launched in 1982 (Hofmann-Wellenhof et al. 2008Hofmann-Wellenhof, B., Lichtenegger, H. and Wasle, E. 2008. GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more. Springer Wien New York. , Jeffrey et al. 2023Jeffrey, C. et al. 2023. An Introduction to GNSS. A primer in using Global Navigation Satellite Systems for positioning and autonomy. Third Edition. Published by Novatel Inc. Hexagon Calgary Campus. https://www.calameo.com/read/0019157965c9cf7056211
https://www.calameo.com/read/0019157965c...
). In 1993, the system was formally declared operational, and in 1995, it was brought to a fully operational constellation (24 GLONASS satellites of the first generation) (GLONASS 2023GLONASS. 2023. Global Navigation Satellite System. Available at: Available at: https://glonass-iac.ru/en/ . [Accessed: September 17, 2023]
https://glonass-iac.ru/en/...
).

Like GPS, GLONASS is a radio-based satellite navigation system that provides users worldwide with free, real-time positioning (P), navigation (N), and timing (T) services in a global geodetic reference system known under the name PZ-90 Parametri Zemli, or Parameters of the Earth (Dach et al. 2015Dach, R., Lutz, S., Walser, P. and Fridez, P. 2015. Bernese GNSS Software Version 5.2. User manual. Astronomical Institute, University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297
https://doi.org/https://doi.org/10.7892/...
; ESA 2011; Revnivykh 2009Revnivykh, S. 2009. Global Navigation Satellite System. 4-th meeting of International Committee on GNSS. 13-18 September 2009, St-Petersburg, the Russian Federation.). The PZ-90 is the system similar to the WGS84 (World Geodetic System) used as the nominal reference for the GLONASS navigation system (Kang et al. 2002Kang, J. et al. 2002. Application of GPS/GLONASS Combination to the Revision of Digital Map. Proceedings of FIG XXII International Congress. Washington. D.C. USA, April 19-26, 2002. ). The main parameters of the PZ-90 system are given in ICD-GLONASS (2016)ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. and Rossbach (2000Rossbach, U. 2000. Positioning and Navigation Using the Russian Satellite System GLONASS. Thesis to obtain the academic degree of doctor engineer. University of the Federal Armed Forces in Munich. https://d-nb.info/969738064/34
https://d-nb.info/969738064/34 ...
).

To enable full global coverage, the GLONASS space segment was planned to contain at least 24 satellites equally distributed in three orbital planes, with eight satellites per plane. The GLONASS constellation satellites are placed into nominally circular orbits with target inclinations of 64.8 degrees and an orbital radius of about 25510 km (Leick 2004Leick, A. 2004. GPS satellite surveying. Third edition, New Jersey: John Wiley & Sons, Inc. ; Revnivykh 2009Revnivykh, S. 2009. Global Navigation Satellite System. 4-th meeting of International Committee on GNSS. 13-18 September 2009, St-Petersburg, the Russian Federation.; Thirsk 2015). Each GLONASS satellite completes the orbit in approximately 11 hours 15 minutes 44 seconds, repeating the geometry every eight sidereal days (Dawidowicz et al. 2015Dawidowicz, K., Kazmierczak, R. and Swiatek, K. 2015. Short static GPS/GLONASS observation processing in the context of antenna phase center variation problem. Boletim de Ciências Geodésicas, 21(1), pp. 213-230. https://doi.org/10.1590/S1982-217020150001000014
https://doi.org/https://doi.org/10.1590/...
).

Each GLONASS satellite transmits a navigation message that contains the positions and velocities of all GLONASS satellites every 30 minutes in the PZ-90 system (Dach et al. 2015Dach, R., Lutz, S., Walser, P. and Fridez, P. 2015. Bernese GNSS Software Version 5.2. User manual. Astronomical Institute, University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297
https://doi.org/https://doi.org/10.7892/...
; Revnivykh 2009Revnivykh, S. 2009. Global Navigation Satellite System. 4-th meeting of International Committee on GNSS. 13-18 September 2009, St-Petersburg, the Russian Federation.). This information is an important parameter in the determination of the user’s position using the trilateration principle (Hofmann-Wellenhof et al. 2008Hofmann-Wellenhof, B., Lichtenegger, H. and Wasle, E. 2008. GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more. Springer Wien New York. ).

Since the launch of the first GLONASS satellite, navigation signals have changed significantly. Traditionally, GLONASS satellites transmit navigational signals on two frequency sub-bands, L1 ~ 1602 MHz and L2 ~ 1246 MHz, relying on the Frequency Division Multiple Access (FDMA) technique (Dawidowicz et al. 2015Dawidowicz, K., Kazmierczak, R. and Swiatek, K. 2015. Short static GPS/GLONASS observation processing in the context of antenna phase center variation problem. Boletim de Ciências Geodésicas, 21(1), pp. 213-230. https://doi.org/10.1590/S1982-217020150001000014
https://doi.org/https://doi.org/10.1590/...
). This is a different technique from the Code Division Multiple Access (CDMA) used by GPS and the other GNSS constellations (Jeffrey et al. 2023Jeffrey, C. et al. 2023. An Introduction to GNSS. A primer in using Global Navigation Satellite Systems for positioning and autonomy. Third Edition. Published by Novatel Inc. Hexagon Calgary Campus. https://www.calameo.com/read/0019157965c9cf7056211
https://www.calameo.com/read/0019157965c...
). The new K-GLONASS satellites were developed to broadcast new signals in the L3~1202.025 MHz band, relying on the CDMA and FDMA techniques (GLONASS 2023GLONASS. 2023. Global Navigation Satellite System. Available at: Available at: https://glonass-iac.ru/en/ . [Accessed: September 17, 2023]
https://glonass-iac.ru/en/...
).

The motion of GLONASS satellites orbiting around the Earth is described by differential equations of second order (Medjahed et al. 2021Medjahed, S. A. et al. 2021. Implementation of the variation of the luni-solar acceleration into GLONASS orbit calculus. Geodetski vestnik, 65(03), pp. 459-471. https://doi: 10.15292/geodetski-vestnik.2021.03.459-471
https://doi.org/https://doi: 10.15292/ge...
); this equation can be solved numerically by the Runge-Kutta integration method recommended in the GLONASS Interface Control Document (ICD-GLONASS 2016ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. , Góral and Skorupa 2015Góral, W. and Skorupa, B. 2015. Calculation of position and velocity of GLONASS satellite based on analytical theory of motion. Artificial Satellites, 50(3), pp. 105-114. https://doi: 10.1515/arsa-2015-0008
https://doi.org/https://doi: 10.1515/ars...
). The ICD-GLONASS is a technical document established by the Coordination Scientific Information Center that gives information about the GLONASS system and specifies parameters of the interface between the GLONASS space segment and the user equipment (Hofmann-Wellenhof et al. 2008Hofmann-Wellenhof, B., Lichtenegger, H. and Wasle, E. 2008. GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more. Springer Wien New York. ; Leick 2004Leick, A. 2004. GPS satellite surveying. Third edition, New Jersey: John Wiley & Sons, Inc. ).

Determining the position of a GLONASS user’s begins with determining the position of at least four satellites in orbit. In effect, the broadcast ephemerides are transmitted at 30-minute intervals. The integration of the differential equation of GLONASS satellite motion using the Runge-Kutta method permits determining at any given time the GLONASS satellite’s position and velocity.

The main objective of this paper is the study of the order and integration step size effect on the position and velocity of GLONASS satellites computed by the Runge-Kutta method and, consequently, on the GLONASS user’s position.

Three different orders of the Runge-Kutta (RK) method were used in this study: order four (RK04); order five (RK05); and order four and five (RK45), and two different step sizes: 1 sec and 0.5 sec, to compute the position and velocity of GLONASS satellites. The runtime or execution time of the scripts (programs) developed in this study is an important factor to take into consideration in the evaluation of the results.

The rest of this paper is organized as follows: In Section 2, the differential equation of GLONASS satellite motion is given, and we present in this section the principle of the Runge-Kutta integration method and the different orders of the Runge-Kutta method used in this study. Section 03 is divided into three sub-sections reserved for the description of the methodology followed in this study and the application, as follows:

Firstly, in Section 3.1, we present the broadcast ephemeris files (data used), the GLONASS satellites in the PZ90-ECEF (Earth-Centered Earth-Fixed) and the transformation from ECEF to the ECI (Earth-Centered Inertial) reference frame.

Secondly, in Section 3.2, we present the forward and backward integration processes used in this study.

Thirdly, in Section 3.3, we present the statistical analysis and discussions of the application results based on the differences between the forward and backward integration processes. The runtime of all programs is given in the form of graphs and discussions. Finally, a summary of conclusions is given.

2. Runge-Kutta Integration of GLONASS Orbits

The orbital motion of a satellite is a result of the Earth’s gravitational attraction; mathematically, the equation of the satellite’s motion is a differential equation of second order given in (Leick 2004Leick, A. 2004. GPS satellite surveying. Third edition, New Jersey: John Wiley & Sons, Inc. ; Lin et al. 2009Lin, Y., Guo, H. and Yu, M. 2009. A comparison for GLONASS satellite coordinate calculation. International Conference on Information Engineering and Computer Science. IEEE. https://DOI: 10.1109/ICIECS.2009.5365110
https://doi.org/https://DOI: 10.1109/ICI...
; Son et al. 2019Son, E. et al. 2019. Comparison of Numerical Orbit Integration between Runge-Kutta and Adams-Bashforth-Moulton using Global Navigation Satellite System Broadcast Ephemeris. Journal of Positioning, Navigation, and Timing, 8(4), pp. 201-208. https://doi.org/10.11003/JPNT.2019.8.4.201
https://doi.org/https://doi.org/10.11003...
) by Equation (1):

γ x = - η 1 X - η 2 1 - η 3 X + w e 2 X + 2 w e V y + γ x - L S γ y = - η 1 Y - η 2 1 - η 3 Y + w e 2 Y - 2 w e V x + γ y - L S γ z = - η 1 Z - η 2 3 - η 3 Z + γ z - L S (1)

In Equation (1), we have the following notation:

  • η1=GM/r3, G is the gravitational constant, M is the mass of Earth and ris the orbital radius =X2+Y2+Z2.

  • η2=(3/2)J2GMa2/r5, J2 is the second order of the harmonic coefficient and a is the equatorial radius.

  • η3=5z2/r2, we is the rotation rate of Earth, x y z ) are the accelerations of satellites, (X,Y,Z) are the coordinates of satellites, (Vx,Vy) are the velocities of satellites and x-LS y-LS z-LS ) are the luni-solar accelerations (sun and moon perturbations).

The constants GM, J 2 , a and w e are given in the ICD-GLONASS (2016)ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. , Son et al. (2019Son, E. et al. 2019. Comparison of Numerical Orbit Integration between Runge-Kutta and Adams-Bashforth-Moulton using Global Navigation Satellite System Broadcast Ephemeris. Journal of Positioning, Navigation, and Timing, 8(4), pp. 201-208. https://doi.org/10.11003/JPNT.2019.8.4.201
https://doi.org/https://doi.org/10.11003...
) and Rossbach (2000Rossbach, U. 2000. Positioning and Navigation Using the Russian Satellite System GLONASS. Thesis to obtain the academic degree of doctor engineer. University of the Federal Armed Forces in Munich. https://d-nb.info/969738064/34
https://d-nb.info/969738064/34 ...
) documents.

To determine the position (Y) and velocity (V) vectors of a GLONASS satellite at time (Tc) as given in Figure 4, we solved the differential Equation (1) of a second-order transformed to a first-order given in (ICD-GLONASS 2016ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. ; Rossbach 2000Rossbach, U. 2000. Positioning and Navigation Using the Russian Satellite System GLONASS. Thesis to obtain the academic degree of doctor engineer. University of the Federal Armed Forces in Munich. https://d-nb.info/969738064/34
https://d-nb.info/969738064/34 ...
; Subirana et al. 2013Subirana, J. S. et al. 2013. GNSS Data Processing: Vol. I Fundamentals and Algorithms. ESA Communications. ISBN 978-92-9221-886-7. https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
https://gssc.esa.int/navipedia/GNSS_Book...
) by:

d X t = V x d Y t = V y d Z t = V z d V x t = - η 1 X - η 2 1 - η 3 X + w e 2 X + 2 w e V y + γ x - L S d V y t = - η 1 Y - η 2 1 - η 3 Y + w e 2 Y + 2 w e V x + γ y - L S d V z t = - η 1 Z - η 2 3 - η 3 Z + γ z - L S (2)

There are many numerical methods used for the integration of the Equation (2); the Runge-Kutta method has been designed to solve first-order differential equations (Es-hagh 2005Es-hagh, M. 2005. Step variable numerical orbit integration of a low earth orbiting satellite. Journal of the Earth & Space Physics, (1), pp. 1-12. ). This method (the R-K method) was developed around 1900 by the German mathematicians Carl Runge (1856-1927) and Martin Wilhlem Kutta (1867-1944) in the context of the resolution of differential equations in the field of atomic spectra (Arora et al. 2020Arora, G., Joshi, V. and Garki, I.S. 2020. Developments in Runge-Kutta method to solve ordinary differential equations. In : Mangey Ram, ed.2020,ed. Recent advances in mathematics for engineering. CRC Press, Taylor & Francis Group pp. 193-202. https://doi.org/10.1201/9780429200304-9
https://doi.org/https://doi.org/10.1201/...
; Vallado 2013Vallado, D. A. 2013. Fundamentals of Astrodynamics and Applications. Fourth Edition. Space Technology Library. Microcosm Press. ). With the numerical development, other special and modified Runge Kutte methods were developed, such as R-K order three and four proposed by Abebe and Gustaf; R-K Fehlberg; R-K order four and five developed by Dormand and Prince and the R-K Merson (Arora et al. 2020Arora, G., Joshi, V. and Garki, I.S. 2020. Developments in Runge-Kutta method to solve ordinary differential equations. In : Mangey Ram, ed.2020,ed. Recent advances in mathematics for engineering. CRC Press, Taylor & Francis Group pp. 193-202. https://doi.org/10.1201/9780429200304-9
https://doi.org/https://doi.org/10.1201/...
; Bradley 2015Bradley, B. K. 2015. Numerical Algorithms for Precise and Efficient Orbit Propagation and Positioning. Ph.D Thesis, Aerospace Engineering Sciences, University of Colorado at Boulder. ; Vallado 2013Vallado, D. A. 2013. Fundamentals of Astrodynamics and Applications. Fourth Edition. Space Technology Library. Microcosm Press. ).

The principle of the R-K method consists in the computation of the approximate value of (Y i+1 ) from the previous value of (Y i ) (Medjahed et al. 2021Medjahed, S. A. et al. 2021. Implementation of the variation of the luni-solar acceleration into GLONASS orbit calculus. Geodetski vestnik, 65(03), pp. 459-471. https://doi: 10.15292/geodetski-vestnik.2021.03.459-471
https://doi.org/https://doi: 10.15292/ge...
). This integration method is efficient, accurate, relatively simple, and easy to implement in programming software such as MATLAB; the integration step size can be easily changed (Son et al. 2019Son, E. et al. 2019. Comparison of Numerical Orbit Integration between Runge-Kutta and Adams-Bashforth-Moulton using Global Navigation Satellite System Broadcast Ephemeris. Journal of Positioning, Navigation, and Timing, 8(4), pp. 201-208. https://doi.org/10.11003/JPNT.2019.8.4.201
https://doi.org/https://doi.org/10.11003...
). The disadvantages of this method include the numerical errors, the instability in solving some problems, the computation time, and the convergence difficulties.

In this study, Equation (2) is solved using the following orders of the R-K method (Arora et al. 2020Arora, G., Joshi, V. and Garki, I.S. 2020. Developments in Runge-Kutta method to solve ordinary differential equations. In : Mangey Ram, ed.2020,ed. Recent advances in mathematics for engineering. CRC Press, Taylor & Francis Group pp. 193-202. https://doi.org/10.1201/9780429200304-9
https://doi.org/https://doi.org/10.1201/...
; Dormand and Prince 1980Dormand, J. R. and Prince, P. J. 1980. A family of embedded Runge-Kutta formulae. Journal of computational and applied mathematics, 6(1), pp. 19-26. https://doi.org/10.1016/0771-050X(80)90013-3
https://doi.org/https://doi.org/10.1016/...
; Es-hagh 2005Es-hagh, M. 2005. Step variable numerical orbit integration of a low earth orbiting satellite. Journal of the Earth & Space Physics, (1), pp. 1-12. ; Lin et al. 2009Lin, Y., Guo, H. and Yu, M. 2009. A comparison for GLONASS satellite coordinate calculation. International Conference on Information Engineering and Computer Science. IEEE. https://DOI: 10.1109/ICIECS.2009.5365110
https://doi.org/https://DOI: 10.1109/ICI...
):

1. The R-K Method’s Fourth Order (RK04)

Y i + 1 = Y i + K 1 + 2 K 2 + 2 K 3 + K 4 6

K 1 = f ( t n , Y n ) K 2 = f ( t n + h / 2 , Y n + h K 1 / 2 ) K 3 = f t n + h 2 , Y n + h K 2 2 K 4 = f ( t n + Y n + h K 3 ) (3)

2. The R-K Method’s Fifth Order (RK05)

Y i + 1 = Y i + ( 7 K 1 + 32 K 3 + 12 K 4 + 32 K 5 + 7 K 6 ) / 90 )

K 1 = f ( t n , Y n ) K 2 = f ( t n + h / 2 , Y n + h K 1 / 2 ) K 3 = f t n + h / 4 , Y n ( ( 3 h K 1 + h K 2 / 16 ) K 4 = f t n h 2 , Y n + h K 3 2 K 5 = f t n + 3 h 4 , Y n + ( - 3 h K 2 + 6 h K 3 + 9 h K 4 / 16 ) K 6 = f t n + h , Y n + ( K 1 + 4 h K 2 + 6 h K 3 - 12 h K 4 + 8 h K 5 / 7 ) (4)

3. The R-K Method’s Four and Five Order (RK45):

Y i + 1 = Y i + ( 35 K 1 / 384 + 500 K 3 / 113 + 125 K 4 / 192 - 2187 K 5 / 6784 + 11 K 6 ) / 84 )

K 1 = f ( X n , Y n ) K 2 = f ( X n + h / 5 , Y n + h K 1 / 5 ) K 3 = f X n + 3 h 10 , Y n + 3 h K 1 / 40 + 9 h K 2 / 40 K 4 = f X n + 4 h / 5 , Y n + 44 h K 1 / 45 - 56 h K 2 / 15 + 32 h K 3 / 9 K 5 = f X n + 8 h 9 , Y n + 19732 h K 1 + 6561 - 25360 h K 2 2187 + 64448 h K 3 / 6561 - 212 h K 4 / 729 K 6 = f X n + h , Y n + 9017 h K 1 3168 - 355 h K 2 33 + 46732 h K 3 5247 + 49 h K 4 176 - 5103 h K 5 / 18656 (5)

The R-K method order four and five (RK45) is available in MATLAB, where it is known as ODE45, Ordinary Differential Equation (Bradley 2015Bradley, B. K. 2015. Numerical Algorithms for Precise and Efficient Orbit Propagation and Positioning. Ph.D Thesis, Aerospace Engineering Sciences, University of Colorado at Boulder. ; Lin et al. 2009Lin, Y., Guo, H. and Yu, M. 2009. A comparison for GLONASS satellite coordinate calculation. International Conference on Information Engineering and Computer Science. IEEE. https://DOI: 10.1109/ICIECS.2009.5365110
https://doi.org/https://DOI: 10.1109/ICI...
).

3. GLONASS Orbits Computation: Methodology and Application

3.1 Data Used in GLONASS Satellite Computation

The broadcast ephemeris are part of the navigation message (Hofmann-Wellenhof et al. 2008Hofmann-Wellenhof, B., Lichtenegger, H. and Wasle, E. 2008. GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more. Springer Wien New York. ). The orbital information of GLONASS satellites is provided in the form of position, velocity, and luni-solar acceleration every 30 minutes in RINEX (Receiver Independent Exchange) format (Góral and Skorupa 2015Góral, W. and Skorupa, B. 2015. Calculation of position and velocity of GLONASS satellite based on analytical theory of motion. Artificial Satellites, 50(3), pp. 105-114. https://doi: 10.1515/arsa-2015-0008
https://doi.org/https://doi: 10.1515/ars...
; Kang et al. 2002Kang, J. et al. 2002. Application of GPS/GLONASS Combination to the Revision of Digital Map. Proceedings of FIG XXII International Congress. Washington. D.C. USA, April 19-26, 2002. ; Krzyżek and Skorupa 2015Krzyżek, R. and Skorupa, B. 2015. The influence of application a simplified transformation model between reference frames ECEF and ECI on to prediction accuracy of position and velocity of GLONASS satellites. Reports on Geodesy and Geoinformatics, 99(1), pp. 19-https://doi.org/10.2478/rgg-2015-0009
https://doi.org/https://doi.org/10.2478/...
). The description of the broadcast ephemeris files is given in Gurtner (2007Gurtner, W. 2007. RINEX : The Receiver Independent Exchange Format. Astronomical Institute, University of Bern. https://files.igs.org/pub/data/format/rinex300.pdf.
https://files.igs.org/pub/data/format/ri...
) and Subirana et al. (2013Subirana, J. S. et al. 2013. GNSS Data Processing: Vol. I Fundamentals and Algorithms. ESA Communications. ISBN 978-92-9221-886-7. https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
https://gssc.esa.int/navipedia/GNSS_Book...
). In Table 1, we present an example of the broadcast ephemerides of GLONASS satellite number 9.

Table 1:
Broadcast ephemeris of GLONASS satellite number 9. (April 20, 2021, at 13.45 p.m.)

The data used in this study is the broadcast ephemeris file downloaded from the International GNSS Service (IGS 2021) and available in the Crustal Dynamics Data Information System (CDDIS 2021) database.

The broadcast orbits of all GLONASS satellites on April 20, 2021, as provided in the broadcast ephemerides file, are shown in Figure 1.

Figure 1:
GLONASS orbits in PZ90-ECEF reference frame.

The left figure shows all GLONASS satellites in the PZ90-ECEF reference frame (April 20, 2021), and the right figure shows the orbits of GLONASS satellites numbers 1 (Plane I), 9 (Plane II) and 17 (Plane III).

Figure 2 shows the broadcast orbits of GLONASS satellites on April 20, 2021, transformed from the ECEF to the ECI frame by Equation (6) given in (ICD-GLONASS 2016ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. ; Krzyżek and Skorupa 2015Krzyżek, R. and Skorupa, B. 2015. The influence of application a simplified transformation model between reference frames ECEF and ECI on to prediction accuracy of position and velocity of GLONASS satellites. Reports on Geodesy and Geoinformatics, 99(1), pp. 19-https://doi.org/10.2478/rgg-2015-0009
https://doi.org/https://doi.org/10.2478/...
; Subirana et al. 2013Subirana, J. S. et al. 2013. GNSS Data Processing: Vol. I Fundamentals and Algorithms. ESA Communications. ISBN 978-92-9221-886-7. https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
https://gssc.esa.int/navipedia/GNSS_Book...
) by:

X a = X e cos θ G e - Y e sin θ G e Y a = X e sin θ G e - Y e cos θ G e Z a = Z e V X a = V X e cos θ G e - V Y e sin θ G e - w e Y a V Y a = V X e sin θ G e - V Y e cos θ G e - w e X a V z a = V Z e γ X a = γ X e cos θ G e - γ y e sin θ G e γ Y a = γ X e sin θ G e - γ y e cos θ G e γ Z a = γ Z e (6)

Where X e , Y e , Z e , Vx e , Vy e , Vz e xe , γ ye , γ ze are the positions, velocities, and luni-solar accelerations of GLONASS satellites and θ Ge is the sidereal time at epoch t given in the Appendix APPENDIX Computation of sidereal time for the transformation from ECEF to ECI (Subirana et al. 2013; Vallado 2013). We have: Ge=G0+wete; (G0: the Greenwich meridian at midnight. G 0 = 24110.54841 + 8640184.812866 T u + 0.093104 T u 2 - 6.2 × 1 0 - 6 T u 3 Tu: Julian Centuries; Tu=JD- 2451545/36525; JD: Julian Date JD=1721013.5+367Y-int7/4Y +intM+9/12+int275M/9+D+60h+m+s/60/1440 Y = 2021 Y e a r ; M = 04 M o n t h ; D = 20 D a y ; m = 15 o r 45 M i n u t e s ; s = 0 S e c o n d The sidereal time computed for all GLONASS satellites (April 20, 2021): .

Figure 2:
GLONASS orbits in ECI reference frame.

The left figure shows the GLONASS constellation (24 satellites) in the ECI reference frame (April 20, 2021), and the right figure shows the orbits of GLONASS satellites numbers 1 (Plane I), 9 (Plane II) and 17 (Plane III).

A satellite’s ground track is its orbit-projected perpendicular to the Earth’s surface; Figure 3 depicts the ground track of GLONASS satellites on September 17, 2023 (GLONASS 2023GLONASS. 2023. Global Navigation Satellite System. Available at: Available at: https://glonass-iac.ru/en/ . [Accessed: September 17, 2023]
https://glonass-iac.ru/en/...
).

Figure 3:
GLOANSS ground tracks of satellites.

3.2 Integration Processes of GLONASS Orbits

In this study, the GLONASS broadcast orbits are computed by integration within an interval of ±15 minutes in the forward integration (+15 min) and backward integration (-15 min), as shown in Figure 4, using the data (GLONASS broadcast ephemerides) presented above.

Figure 4:
GLONASS orbits integration processes.

  1. In the forward integration process, the computation of Y = (X, Y, Z) and V = (Vx, Vy, Vz) at time Tc = H: 30 min or Tc = H: 00 min from Tb1 with Tb1  Tc.

  2. In the backward integration process, the computation of Y = (X, Y, Z) and V = (Vx, Vy, Vz) at time Tc = H: 30 min or Tc = H: 00 min from Tb2 with Tb2  Tc.

The forward and backward integration processes are given in Kang et al. (2002Kang, J. et al. 2002. Application of GPS/GLONASS Combination to the Revision of Digital Map. Proceedings of FIG XXII International Congress. Washington. D.C. USA, April 19-26, 2002. ), Medjahed et al., (2022), and Rossbach (2000Rossbach, U. 2000. Positioning and Navigation Using the Russian Satellite System GLONASS. Thesis to obtain the academic degree of doctor engineer. University of the Federal Armed Forces in Munich. https://d-nb.info/969738064/34
https://d-nb.info/969738064/34 ...
).

In the document ICD-GLONASS (2016)ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. , a simplified algorithm is given (Algorithm J-2, page 53), where the luni-solar acceleration is considered constant on an interval of 15 min (ICD-GLONASS 2016ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. , Leick 2004Leick, A. 2004. GPS satellite surveying. Third edition, New Jersey: John Wiley & Sons, Inc. ). It can also be considered variable during the integration interval if we apply the precise algorithm given in the document ICD-GLONASS (2016)ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow. (Algorithm J-1, Page 45).

The steps of GLONASS orbit computation followed in this study are represented in the following diagram:

Figure 5:
GLONASS orbits integration steps.

3.3 Results and Discussions of GLONASS Orbits Integration

In this section, we present the results of integration. Figure 6 shows the integrated orbit (in forward and backward by RK05 and h = 1 sec) and the orbit provided in the data file (broadcast ephemerides) every 30 minutes.

Figure 6:
Orbits of GLONASS satellites number 9 (April 20, 2021).

The statistics of differences between the forward and backward integration processes of GLONASS satellite number 9 on the X, Y, and Z-axes in position and velocity are given in Table 2.

Table 2:
Statistic of GLONASS satellite number 9 during 24 hours (20/04/2021).

On the X-axis, the minimum value was obtained with RK05 (0.006 m) and the maximum value was obtained with RK04 (5.128 m). On the Y-axis, the minimum value was obtained with RK45 (0.000 m) and the maximum value was obtained with RK05 (2.640 m). On the Z-axis, the same minimum (0.05 m) and maximum (1.505 m) were found with RK04, RK05 and RK45. On the velocity components, the same minimum and maximum values were found with RK04, RK05, and RK45 on the X, Y and Z-axes.

The results of the 3D-difference (Min and Max) in position of all GLONASS satellites between forward and backward integration are shown in Figure 7; where the 3D-Diff in position =ΔX2+ΔY2+ΔZ2 and Δ is the absolute difference between forward and backward integration on the X, Y and Z-axes.

Figure 7:
Difference in position (all GLONASS satellites).

The results of the 3D differences in position between forward and backward integration are visually identical during the integration interval (24 hours) for all GLONASS satellites. Table 3 shows the statistics of the differences (forward and backward) presented in Figure 7.

Table 3:
Statistics of all GLONASS satellites (20/04/2021).

The minimum difference was obtained by RK45 (0.041 m); this is due to the precision of the Equation (5) and the number of K-coefficients. The maximum differences obtained are metric which is acceptable for the computation of the broadcast orbits.

The differences obtained are due to several factors, such as the precision of the broadcast orbits and the precision of the R-K method. The luni-solar acceleration is taken constant during 15 minutes of integration. Considering this perturbation as a variable can improve the position and velocity of GLONASS satellites.

The results of RK05 and RK45 are almost identical, with a slight priority for RK45. This is because the equations (4) and (5) of these two orders are similar.

The 3D difference between RK04 recommended in ICD-GLONASS and the other orders does not exceed a few centimeters; this difference is negligible compared to the dimensions of GLONASS orbits.

According to the results of Table 3 (average results), the orders are classified as follows: firstly, order four (RK04); secondly, order four and five (R45); and thirdly, order five (RK05).

Figure 8 displays the runtime (computation or execution time) of the scripts developed in this study of 1 epoch (15 minutes) and 1000 iterations.

Figure 8:
The runtime of 1000 iterations in second.

The runtime is in function of the R-K order and the integration step size; order four (RK04) is faster in execution compared to the other orders because the simplicity of their equations and the number of K-coefficients (4 per iteration).

Due to the similarity of equations and the number of K-coefficients (6 per iteration), the RK05 and RK45 are almost identical, with a slight priority for RK05, the order four and five of R-K method (RK45) is slower because of the complexity of their equation.

The mean execution time is almost twice when the integration step is divided as follows:

  • Step size = 1 sec is 0.020 sec, 0.076 sec, and 0.081 sec for RK04, RK05, and RK45, respectively.

  • Step size = 0.5 sec is 0.040 sec, 0.160 sec, and 0.179 sec for RK04, RK05, and RK45, respectively.

According to the results of the runtime, the orders are classified as follows: firstly, order four (RK04); secondly, order five (RK05); and thirdly, order four and five (RK45).

Finally, in terms of 3D difference and runtime results, RK04 offers the best results in integration and runtime, followed by RK05 and then RK45.

4. Conclusions

In this paper, we studied the effect of the order and step size on the integration of GLONASS orbits using the Runge-Kutta method. The runtime is an important factor to take into consideration.

The order effect of the R-K method is less important and slightly changes the results and the effect of the step size is very important on the runtime. Reducing the integration step significantly increases the runtime.

To profit from the advantages of the fourth order, including ease of implementation in software due to the simplicity of the equations and the reduced execution time (runtime), it is recommended to use order four (RK04) of the Runge-Kutta integration method for the smaller integration steps (1 sec, 0.5 sec, 0.25 sec, etc.) . The order five (RK05) and the order four and five (RK45) of the Runge-Kutta integration method can be used for the big integration steps (10 sec, 15 sec, 20 sec, etc.).

The comparison between the broadcast GLONASS orbits computed by the R-K method and the precise GLONASS orbits can evaluate the results of our study. However, this comparison requires, if we use the IGS products, the transformation of the precise GLONASS orbits in the PZ-90 reference frame and taking into account the difference between GLONASS time the reference of the broadcast ephemerides and GPS time the reference of the precise ephemerides. The precise ephemerides are available at 15-minute intervals, whereas the broadcast ephemerides are available at 30-minute intervals. This requires the use of an interpolation method such as Lagrange interpolation to get the precise positions and velocities of the GLONASS satellites at any given time (time of broadcast ephemerides).

REFERENCES

  • Arora, G., Joshi, V. and Garki, I.S. 2020. Developments in Runge-Kutta method to solve ordinary differential equations. In : Mangey Ram, ed.2020,ed. Recent advances in mathematics for engineering CRC Press, Taylor & Francis Group pp. 193-202. https://doi.org/10.1201/9780429200304-9
    » https://doi.org/https://doi.org/10.1201/9780429200304-9
  • Bradley, B. K. 2015. Numerical Algorithms for Precise and Efficient Orbit Propagation and Positioning Ph.D Thesis, Aerospace Engineering Sciences, University of Colorado at Boulder.
  • Dach, R., Lutz, S., Walser, P. and Fridez, P. 2015. Bernese GNSS Software Version 5.2. User manual Astronomical Institute, University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297
    » https://doi.org/https://doi.org/10.7892/boris.72297
  • Dawidowicz, K., Kazmierczak, R. and Swiatek, K. 2015. Short static GPS/GLONASS observation processing in the context of antenna phase center variation problem. Boletim de Ciências Geodésicas, 21(1), pp. 213-230. https://doi.org/10.1590/S1982-217020150001000014
    » https://doi.org/https://doi.org/10.1590/S1982-217020150001000014
  • Dormand, J. R. and Prince, P. J. 1980. A family of embedded Runge-Kutta formulae. Journal of computational and applied mathematics, 6(1), pp. 19-26. https://doi.org/10.1016/0771-050X(80)90013-3
    » https://doi.org/https://doi.org/10.1016/0771-050X(80)90013-3
  • ESA. 2011. European Space Agency GLONASS general introduction. Available at: Available at: https://gssc.esa.int/navipedia/index.php [Accessed: September 04, 2023]
    » https://gssc.esa.int/navipedia/index.php
  • Es-hagh, M. 2005. Step variable numerical orbit integration of a low earth orbiting satellite. Journal of the Earth & Space Physics, (1), pp. 1-12.
  • GLONASS. 2023. Global Navigation Satellite System Available at: Available at: https://glonass-iac.ru/en/ [Accessed: September 17, 2023]
    » https://glonass-iac.ru/en/
  • Góral, W. and Skorupa, B. 2015. Calculation of position and velocity of GLONASS satellite based on analytical theory of motion. Artificial Satellites, 50(3), pp. 105-114. https://doi: 10.1515/arsa-2015-0008
    » https://doi.org/https://doi: 10.1515/arsa-2015-0008
  • Gurtner, W. 2007. RINEX : The Receiver Independent Exchange Format Astronomical Institute, University of Bern. https://files.igs.org/pub/data/format/rinex300.pdf
    » https://files.igs.org/pub/data/format/rinex300.pdf
  • Hofmann-Wellenhof, B., Lichtenegger, H. and Wasle, E. 2008. GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and more Springer Wien New York.
  • ICD-GLONASS. 2016. General Description of Code, Interface Control. Document .Global Navigation Satellite System. Edition 1.0. Moscow.
  • Jeffrey, C. et al. 2023. An Introduction to GNSS. A primer in using Global Navigation Satellite Systems for positioning and autonomy Third Edition. Published by Novatel Inc. Hexagon Calgary Campus. https://www.calameo.com/read/0019157965c9cf7056211
    » https://www.calameo.com/read/0019157965c9cf7056211
  • Kang, J. et al. 2002. Application of GPS/GLONASS Combination to the Revision of Digital Map. Proceedings of FIG XXII International Congress Washington. D.C. USA, April 19-26, 2002.
  • Krzyżek, R. and Skorupa, B. 2015. The influence of application a simplified transformation model between reference frames ECEF and ECI on to prediction accuracy of position and velocity of GLONASS satellites. Reports on Geodesy and Geoinformatics, 99(1), pp. 19-https://doi.org/10.2478/rgg-2015-0009
    » https://doi.org/https://doi.org/10.2478/rgg-2015-0009
  • Leick, A. 2004. GPS satellite surveying Third edition, New Jersey: John Wiley & Sons, Inc.
  • Lin, Y., Guo, H. and Yu, M. 2009. A comparison for GLONASS satellite coordinate calculation. International Conference on Information Engineering and Computer Science IEEE. https://DOI: 10.1109/ICIECS.2009.5365110
    » https://doi.org/https://DOI: 10.1109/ICIECS.2009.5365110
  • Medjahed, S. A. et al. 2021. Implementation of the variation of the luni-solar acceleration into GLONASS orbit calculus. Geodetski vestnik, 65(03), pp. 459-471. https://doi: 10.15292/geodetski-vestnik.2021.03.459-471
    » https://doi.org/https://doi: 10.15292/geodetski-vestnik.2021.03.459-471
  • Revnivykh, S. 2009. Global Navigation Satellite System. 4-th meeting of International Committee on GNSS 13-18 September 2009, St-Petersburg, the Russian Federation.
  • Rossbach, U. 2000. Positioning and Navigation Using the Russian Satellite System GLONASS Thesis to obtain the academic degree of doctor engineer. University of the Federal Armed Forces in Munich. https://d-nb.info/969738064/34
    » https://d-nb.info/969738064/34
  • Son, E. et al. 2019. Comparison of Numerical Orbit Integration between Runge-Kutta and Adams-Bashforth-Moulton using Global Navigation Satellite System Broadcast Ephemeris. Journal of Positioning, Navigation, and Timing, 8(4), pp. 201-208. https://doi.org/10.11003/JPNT.2019.8.4.201
    » https://doi.org/https://doi.org/10.11003/JPNT.2019.8.4.201
  • Subirana, J. S. et al. 2013. GNSS Data Processing: Vol. I Fundamentals and Algorithms ESA Communications. ISBN 978-92-9221-886-7. https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
    » https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
  • Vallado, D. A. 2013. Fundamentals of Astrodynamics and Applications Fourth Edition. Space Technology Library. Microcosm Press.

DATA

  • CDDIS. 2021. Crustal Dynamics Data Information System. GNSS Orbit Products. Available at: https://cddis.nasa.gov/Data_and_Derived_ Products/GNSS/orbit_products.html. [Accessed: April 30, 2021].
  • IGS. 2021. International GNSS Service. Available at: https://igs.org/. [Accessed: April 30, 2021]

APPENDIX

Computation of sidereal time for the transformation from ECEF to ECI (Subirana et al. 2013; Vallado 2013). We have: Ge=G0+wete; (G0: the Greenwich meridian at midnight.

G 0 = 24110.54841 + 8640184.812866 T u + 0.093104 T u 2 - 6.2 × 1 0 - 6 T u 3

Tu: Julian Centuries; Tu=JD- 2451545/36525; JD: Julian Date JD=1721013.5+367Y-int7/4Y +intM+9/12+int275M/9+D+60h+m+s/60/1440

Y = 2021 Y e a r ; M = 04 M o n t h ; D = 20 D a y ; m = 15 o r 45 M i n u t e s ; s = 0 S e c o n d

The sidereal time computed for all GLONASS satellites (April 20, 2021):

Publication Dates

  • Publication in this collection
    30 Aug 2024
  • Date of issue
    2024

History

  • Received
    23 Mar 2023
  • Accepted
    22 July 2024
Universidade Federal do Paraná Centro Politécnico, Jardim das Américas, 81531-990 Curitiba - Paraná - Brasil, Tel./Fax: (55 41) 3361-3637 - Curitiba - PR - Brazil
E-mail: bcg_editor@ufpr.br