Calculating Global Navigation Satellite System Satellite Velocities and Accelerations by Utilizing the Orbit Fitting and Orbit Integration Methods

: The high-precision satellite velocities and accelerations calculated by the Global Navigation Satellite System (GNSS) are essential for tasks such as airborne gravity data processing. Users generally interpolate satellite positions in the precise ephemeris to calculate satellite velocity and acceleration. However, due to the edge effect, the accuracy of the interpolation is relatively low near day boundaries. In this study, a method for calculating GNSS satellite velocity and acceleration based on orbit fitting and orbit integration was proposed, and the high-precision transformation relationship between satellite velocity and acceleration in the Earth-Centered Inertial (ECI) coordinate system and the Earth-Centered, Earth-Fixed (ECEF) coordinate system was derived. The experimental results show that the satellite velocity accuracy is 1.5 × 10 − 6 m/s and the acceleration accuracy is 1.0 × 10 − 8 m/s 2 according to the proposed method. Thus, the proposed method improves the accuracy of calculating satellite velocity and acceleration near day boundaries, and helps GNSS users to obtain satellite velocity and acceleration information with consistent precision throughout the day.


Introduction
The satellite velocities and accelerations calculated by the Global Navigation Satellite System (GNSS) provide important information about satellite orbits [1][2][3] and are also an important part of GNSS data services.High-precision GNSS satellite velocities and accelerations are crucial for tasks such as airborne gravity measurement data processing [4][5][6], gravity field reproduces with Low-Earth-Orbit (LEO) satellites [7], and motion carrier state monitoring [8][9][10].Currently, the methods commonly utilized [11] to calculate GNSS satellite velocities and accelerations include the analytical method based on broadcast ephemeris and the numerical interpolation method based on precise ephemeris.
The broadcast ephemeris is an extrapolated forecast ephemeris, whose orbit accuracy is much lower than that of the precise ephemeris, and its parameters ignore the second-order derivatives of orbit inclination and ascending node longitude.Therefore, the accuracy of satellite velocity and acceleration calculated with the broadcast ephemeris is relatively low.For instance, the accuracy of satellite velocity calculated based on the broadcast ephemeris in the literature [12] is about 1.0 × 10 −3 m/s and the acceleration accuracy is about 1.0 × 10 −4 m/s 2 , which cannot meet the application requirements of higher accuracy.
The precise ephemeris is generally obtained with post-processing and contains daily tabular satellite positions and satellite clock offsets.The final or the rapid precise ephemerides provided by the International GNSS Service (IGS) undoubtedly contain the most accurate Remote Sens. 2024, 16, 2366 2 of 14 satellite velocity and acceleration information.The National Geodetic Survey (NGS) recommends utilizing sliding-window Lagrange interpolation on the precise ephemeris to compute the satellite positions, velocities, and accelerations at the required epochs [13].However, the use of Lagrange interpolation on the precise ephemeris is affected by the Runge phenomenon [14].Namely, the interpolation error at the edges of the interpolation arc is relatively large [15].Therefore, because the IGS provides precise ephemerides on a daily basis, the interpolation accuracy is relatively low near day boundaries.In the above situation, users can generally choose two methods for the interpolation near the day boundaries.One is to accept the interpolation results at the edge of the interpolated arc and extrapolate them if necessary (for the precise ephemeris that does not provide the 24:00 record).The other is to combine the precise ephemerides of adjacent days, so that the calculated time is roughly located in the middle of the interpolation arc.However, the Runge phenomenon and the extrapolation error of orbit interpolation are both large [16].Furthermore, as they are affected by the GNSS orbit determination error, the precise ephemerides of the adjacent days are not continuous.The discontinuities can reduce the accuracy of orbit interpolation [17].In summary, the accuracy of the satellite velocity and acceleration calculated with the precise ephemeris utilizing the current method is not uniform within a single day, and the calculation accuracy near the day boundary is much lower than that of the other periods, which affects the high-precision application of the satellite velocity and acceleration throughout the day.Therefore, it is necessary to study the GNSS satellite velocity and acceleration calculation method, which can realize uniform accuracy within a single day.
The orbit integration method [18,19] adopts rigorous force models and a numerical integration algorithm, which can realize single-day precise ephemeris fitting and shortterm orbit extrapolation with a high accuracy (millimeter level), and this method has been widely utilized in Precise Point Positioning (PPP) [20,21] and precise satellite clock offset estimation [22,23].However, except for Precise Orbit Determination (POD), the current applications of orbit integration are mainly utilized for obtaining continuous precise satellite positions in a single day, while the applications and studies of calculating satellite velocity and acceleration utilizing orbit integration are insufficient.Therefore, considering the better orbit fitting and extrapolation performance of orbit integration, this study investigated the use of orbit integration to accomplish the calculation of GNSS satellite velocity and acceleration with uniform precision within a single day.
This article is organized as follows: the Lagrange method to interpolate satellite velocities and accelerations is presented in Section 2, and the method for calculating GNSS satellite velocities and accelerations utilizing orbit fitting and orbit integration is introduced in Section 3. Section 4 describes the experiment, including the utilized data and processing strategies.Section 5 presents the experimental results and analysis.Finally, a summary of the work and the conclusions are provided in Section 6.

Lagrange Interpolation Method
To calculate the GNSS satellite positions at the desired epochs utilizing the precise ephemeris, an interpolation method is generally adopted on a set of data points.The interpolation methods include Lagrange interpolation, trigonometric polynomial interpolation, and Chebyshev polynomial interpolation.Among these algorithms, the Lagrange polynomial accurately passes through the data points and is recommended by the NGS [13].By differentiating the Lagrange interpolation formula, the first-order and second-order derivatives can be obtained: Remote Sens. 2024, 16, 2366 ..
where (t i , x(t i )) represents the data points utilized for interpolation; x(t i ) denotes the known satellite positions (X, Y, or Z component) in the precise ephemeris; .
x(t) represent the satellite velocity and acceleration of a certain component at epoch t, respectively.By utilizing Equations ( 1) and ( 2), the satellite velocity and acceleration can be interpolated.

Orbit Fitting and Orbit Integration Methods
Currently, solving the satellite motion equations with orbit integration is the main method utilized to obtain high-precision satellite orbits [18].In order to realize more accurate orbit integration, the known information is needed to determine the orbit state parameters at the initial moment, i.e., orbit fitting.In the following, the principles of orbit integration and orbit fitting are described.

Orbit Integration
Orbit integration mainly refers to the process of numerical integration of the orbit utilizing a set of orbit state parameters and force models.The motion equation of a GNSS satellite in the Earth-Centered, Inertial (ECI) coordinate system [18] can be expressed as where x = (r, v, p) T represents the state vector composed of satellite position, satellite velocity and force model parameters; x 0 represents the state vector at the initial moment of integration.The satellite motion equation is a first-order differential equation, which can be solved by utilizing the numerical integration method after the initial moment state vector and force model parameters are determined.The integration of the satellite motion equation can only obtain the satellite state vector at other moments, i.e., x * the reference orbit, and the following formula should also be integrated if POD is to be carried out: .
When the orbit integration is implemented, the numerical integration of extended state vectors is generally carried out to realize the integration of Equations ( 3) and (4) at the same time.According to the integration step size and the number of known data points, the numerical integration can be generally categorized into a single-step method with a variable step size and a multi-step method with a fixed step size [18].The single-step method does not need to record information about data points beyond the integration moment and is mostly utilized as an initial step in the multi-step method.The multi-step method has a high computational efficiency and is mostly utilized for orbit calculation.In this study, we mainly utilize the seventh (eighth) order variable step-size Runge-Kutta-Fehlberg (RKF) single-step method and the tenth order fixed-step-size Adams linear multi-step method for the numerical integration of orbits.

Orbit Fitting
Orbit fitting mainly refers to the process of improving the state vector at an initial moment by utilizing the known time series of satellite positions.The improved state vector can accurately describe the satellite motion within the fitted orbit arc and the short-term extrapolated orbit arc.Utilizing a numerical integration of the extended states, the reference orbit x * (t) and the state transfer matrix Ψ(t, t 0 ) at time t can be obtained.Let r(t) represent the known satellite positions in the ECI coordinate system at time t; then, the orbit fitting observation equation can be described as where Ψ r (t, t 0 ) is the first three rows of Ψ(t, t 0 ), i.e., the partial derivatives of the satellite position vector with respect to the state vector at the initial moment, and ε(t) denotes the observation noise.The observed value in Equation ( 6) is the ECI satellite position (X, Y, Z components) minus the reference orbit.Among them, the reference orbit is generated by orbit integration, and the ECI satellite position is obtained through coordinate transformation from the Earth-Centered, Earth-Fixed (ECEF) satellite position in the precise ephemeris.
Considering the above observations to be independently and equally weighted, the normal equation of orbit fitting can be expressed as where . After superim- posing the normal equations of multiple epochs, the unknown parameters ∆x 0 can be solved with the least squares method, and the state vector at the initial moment can be improved.
Figure 1 shows detailed steps for satellite orbit integration and orbit fitting.It can be seen that orbit integration and orbit fitting complement each other.Orbit integration is a prerequisite of orbit fitting, and by utilizing orbit fitting, the improved initial moment state vector required by orbit integration can be obtained.

Transformation of ECI and ECEF for Velocity and Acceleration
After orbit fitting, the reference orbit generated by the orbit integration contains the exact ECI satellite positions and velocities, and certain force model parameters are optimized together.However, what the user needs is the satellite velocity and acceleration in the Earth-Centered, Earth-Fixed (ECEF) coordinate system.Therefore, the satellite velocity and acceleration need to be transformed from ECI to ECEF after the orbit integration is completed.

Transformation of ECI and ECEF for Velocity and Acceleration
After orbit fitting, the reference orbit generated by the orbit integration contains the exact ECI satellite positions and velocities, and certain force model parameters are optimized together.However, what the user needs is the satellite velocity and acceleration in the Earth-Centered, Earth-Fixed (ECEF) coordinate system.Therefore, the satellite velocity and acceleration need to be transformed from ECI to ECEF after the orbit integration is completed.
Since most of the commonly utilized POD software, including GAMIT 10.71 [24], and PANDA [25], adopts an approximate method to transform the ECI and ECEF for satellite velocity and acceleration, only the rate of the Earth Rotation Angle (ERA) is considered, and its transformation accuracy is limited.Moreover, technical documents such as the Standards of Fundamental Astronomy (SOFA) [26] and the International Earth Rotation and Reference Systems Service (IERS) convention [27] do not provide transformation methods.Therefore, several formulas are derived below to achieve high-precision ECI-to-ECEF transformation of satellite velocity and acceleration.
According to the IERS convention [27], at the observation time t , the transformation relationship between the ECI position vector ECI r and the ECEF position vector ECEF r can be described as where is the Earth rotation matrix, and is the polar motion matrix.
By taking the first-order and second-order derivatives of the above equation, Since most of the commonly utilized POD software, including GAMIT 10.71 [24], and PANDA [25], adopts an approximate method to transform the ECI and ECEF for satellite velocity and acceleration, only the rate of the Earth Rotation Angle (ERA) is considered, and its transformation accuracy is limited.Moreover, technical documents such as the Standards of Fundamental Astronomy (SOFA) [26] and the International Earth Rotation and Reference Systems Service (IERS) convention [27] do not provide transformation methods.Therefore, several formulas are derived below to achieve high-precision ECI-to-ECEF transformation of satellite velocity and acceleration.
According to the IERS convention [27], at the observation time t, the transformation relationship between the ECI position vector r ECI and the ECEF position vector r ECEF can be described as where Q(t) is the precession-nutation matrix, R(t) is the Earth rotation matrix, and W(t) is the polar motion matrix.By taking the first-order and second-order derivatives of the above equation, ..

R(t)
. (10) in which v and a denote the velocity and acceleration vectors, respectively;

R(
Similarly, the second-order derivatives of (UT1 − UTC) can also be calculated with the interpolation method, and then, E

W(t) can be obtained by interpolating the corresponding precession-nutation matrices and polar motion matrices.
Utilizing Equations ( 9) and ( 10) enables the high-precision transformation of ECEF to ECI for satellite velocity and acceleration measurement.The transformations of ECI to ECEF and ECEF to ECI are the inverse transformations of each other.After obtaining the transformation matrix of ECEF to ECI, the transformation matrix of ECI to ECEF can be obtained by performing inverse operations on the transformation matrix.
Equations ( 11) and (12) give the approximate transformation methods commonly utilized in POD software: In contrast to the transformation method proposed in this study, the approximate method ignores the first-order derivatives of the precession-nutation matrix and the polar motion matrix for velocity transformation, and further ignores the first-order derivatives of the Earth rotation matrix as well as all second-order derivative terms for acceleration transformation.

Data and Processing Strategies
In order to assess the effectiveness of the satellite velocity and acceleration calculation method based on orbit integration, experimental validation is carried out in terms of orbit fitting and orbit extrapolation, respectively.The experimental data and processing strategies of the two types of experiments are presented below.

Orbit Fitting
The BeiDou orbit data in the Multi-GNSS Experiment (MGEX) [28] precise ephemeris of Wuhan University from 001 days to 120 days in 2022 were selected for the experiment, with a product sampling rate of 15 min and a time span of 00:00~23:45 (GPS Time, GPST).The reason for choosing BeiDou is that the system adopts three kinds of orbit satellites: Geostationary Earth Orbit (GEO), Inclined Geosynchronous Orbit (IGSO), and Medium Earth Orbit (MEO).Other GNSS systems only adopt one kind of orbit satellite: MEO.Therefore, the data used in this experiment are universal.
In data processing, the single-day precise ephemeris is firstly utilized as a virtual observation to improve the initial state parameters, solar radiation pressure model parameters, and velocity pulse parameters for each satellite.Then, the improved satellite parameters are utilized to perform the numerical integration, and the satellite velocity and acceleration in the ECEF coordinate system are output.The force models utilized in orbit fitting and orbit integration are shown in Table 1.Finally, the satellite velocity and acceleration calculated by the 10th-order sliding-window Lagrange interpolation are utilized as the reference to evaluate the calculation accuracy of the orbit integration.To ensure the interpolation accuracy and the density of check points, the check points are at 1 min intervals and are located in the middle of each sliding window.

Orbit Extrapolation
The time span of the precise ephemeris provided by IGS and its analysis centers is mostly 00:00~23:45 or 00:00~23:55.In order to calculate the satellite velocity and acceleration during the time period of 23:45~24:00 or 23:55~24:00, the orbit integration should be extrapolated for 15 min or 5 min.The orbit extrapolation experiment adopts the same experimental data as those used in Section 3.1.[31] Every two hours Relativity [18] Considering Antenna thrust [32] Considering Earth orientation IERS C04 (http://celestrak.com/SpaceData/EOP-format.asp (accessed on 10 April 2024))

Results and Analysis
This section presents and analyzes the experimental results of orbit fitting and orbit extrapolation.By comparing the performances of different strategies, the effectiveness of calculating satellite velocity and acceleration by orbit fitting and orbit integration is verified.

Orbit Fitting Results
Figures 2 and 3 show the single-day root mean square (RMS) of satellite velocities and accelerations calculated by utilizing orbit fitting and orbit integration, respectively.The results were further classified according to GEO, IGSO, and MEO.As can be seen from Figure 2, the satellite velocity discrepancy calculated by utilizing orbit integration and sliding window Lagrange interpolation is small, and each component of the three types of satellites is less than 1.5 × 10 −6 m/s.Furthermore, the calculation accuracy fluctuates between different days.This is mainly caused by the satellite orbit accuracy of that day, the corresponding orbit fitting accuracy and the transformation accuracy between ECI and ECEF.
Remote Sens. 2024, 16, x FOR PEER REVIEW 8 of 14    To verify the correctness and accuracy of the transformation relationship between ECI and ECEF for satellite velocity and acceleration in this study, the multi-day average RMS for satellite velocity and acceleration calculated utilizing the approximate transformation model and the transformation model are as presented in Tables 2 and 3, respectively.For comparison, Table 4 shows the calculation accuracy when the first-order and second-order derivatives of ( )  From Figure 3, it can be seen that the acceleration accuracy of each component for GEO satellites is less than 1.0 × 10 −8 m/s 2 , and the corresponding accuracies for IGSO and MEO satellites are less than 0.5 × 10 −8 m/s 2 .Therefore, the satellite velocity and acceleration calculated by the orbit fitting and orbit integration methods in this study have a high accuracy and can meet the requirements of high-precision applications.
To verify the correctness and accuracy of the transformation relationship between ECI and ECEF for satellite velocity and acceleration in this study, the multi-day average RMS for satellite velocity and acceleration calculated utilizing the approximate transformation model and the transformation model are as presented in Tables 2 and 3, respectively.For comparison, Table 4 shows the calculation accuracy when the first-order and second-order derivatives of (UT1 − UTC), Q(t), and W(t) are ignored in the new transformation model.In this case, the mathematical interpolation operation for the above three items is not required when calculating the transformation matrix.
As can be seen from Tables 2 and 3, the satellite velocity accuracy calculated utilizing the ECI and ECEF transformation model in this study is greatly improved compared with that of the approximate transformation model.Except for the v x and v y components of the GEO satellites, which are improved by about 70%, and the other satellites, the enhancement of each component is by almost close to or more than one order of magnitude.The accuracy of satellite acceleration calculated utilizing different transformation models varies greatly, and for the a z components, the model in this study improves the accuracy by one order of magnitude compared with the approximate model.For the a x and a y components, the computational error of the approximate model is large and reaches the 10 −2 m/s 2 level, while the computational error was reduced to the 10 −9 m/s 2 level when the transformation model in this study was adopted.This result is mainly caused by the neglect of the first-and second-order derivatives of the Earth rotation matrix, namely .R(t) and ..

R(t).
Furthermore, due to the Earth's rotation around the Z-axis, neglecting the above two items mainly affects the x and y components.As a result, the acceleration of the x and y components in Table 2 is much larger than that of the z component.As can be seen from Tables 3 and 4, when the first and second order derivatives of (UT1 − UTC), Q(t), and W(t) are neglected, the calculated satellite velocity and acceleration accuracy is lower than that of the rigorous transformation model in this study.Among them, the accuracies of the GEO satellite v x and v y components are reduced by about four times, and the accuracies of the v z components and the velocities of the IGSO and MEO satellites are reduced by at least one order of magnitude.For the acceleration accuracies, except for the GEO a x and a y components, which are basically unchanged, the a z component and IGSO and MEO accuracies are reduced by a factor of about three.However, since the mathematical interpolation operation is avoided by neglecting the three first-and second-order derivatives, the above simplified model can also be utilized when meeting the accuracy requirements of practical applications.
In addition, in order to improve the orbit fitting accuracy, the segmented velocity pulse was added in this study.To analyze whether the inclusion of the velocity pulse helps to improve the calculation accuracy of satellite velocity and acceleration, the multi-day average RMS values of the satellite velocity and acceleration for the non-eclipse seasons [33] and the eclipse seasons were calculated utilizing the transformation model of this study, and the results are presented in Tables 5 and 6, respectively.
As can be seen from Tables 5 and 6, the multi-day average RMS in the non-eclipse season is lower than that in the eclipse season, and the satellite velocity accuracy measured by utilizing the segmented velocity pulse is greatly improved compared to that without the velocity pulse, especially for the GEO satellites and the IGSO satellites in the eclipse season.Therefore, the inclusion of the velocity pulse parameter should be considered in practical applications.However, the addition of the velocity pulse does not improve the acceleration accuracy significantly.This may be because small changes in satellite position and velocity have little effect on acceleration calculations.

Orbit Extrapolation Results
In order to avoid the edge effect of Lagrange interpolation, the precise ephemeris from 00:00 to 22:45 for each day is utilized for orbit fitting and then extrapolated for 15 min with orbit integration to calculate the satellite velocity and acceleration between 22:45 and 23:00.The initial moment of orbit integration is chosen to be 00:00.Since the extrapolation time period is far away from the edge of the interpolation arc at 21:15~23:45, the satellite velocity and acceleration can still be calculated reliably utilizing Lagrange interpolation, and the interpolation result can be utilized as a reference to evaluate the extrapolation accuracy of orbit integration.Meanwhile, the extrapolation results from the orbit integration are compared with the extrapolation results from the Lagrange fitting, which is different from using the Lagrange interpolation as a reference.The specific orbit extrapolation experiment scheme is shown in Figure 4.  Figure 5 shows the 3D RMS of the satellite velocity and acceleration extrapolated utilizing orbit fitting and Lagrange fitting, respectively.It can be clearly seen that the satellite velocity and acceleration calculated by utilizing Lagrange fitting are affected by the edge effect, and the interpolation error increases earlier than in the interpolation of satellite positions [16,17].Furthermore, the satellite velocity and acceleration error extrapolated utilizing Lagrange fitting is large compared to orbit fitting.As shown in Table 7, the 15 min extrapolation error for velocity can reach 1~2 mm/s, and the acceleration error can reach 10 −6 m/s 2 .By contrast, the extrapolation accuracy of orbit integration is high, and there is no significant increase in errors within 15 min for both satellite velocity and acceleration.Therefore, by utilizing orbit fitting, orbit integration, and orbit extrapolation, high-precision satellite velocities and accelerations can be obtained, albeit near the day boundary.Figure 5 shows the 3D RMS of the satellite velocity and acceleration extrapolated utilizing orbit fitting and Lagrange fitting, respectively.It can be clearly seen that the satellite velocity and acceleration calculated by utilizing Lagrange fitting are affected by the edge effect, and the interpolation error increases earlier than in the interpolation of satellite positions [16,17].Furthermore, the satellite velocity and acceleration error extrapolated utilizing Lagrange fitting is large compared to orbit fitting.As shown in Table 7, the 15 min extrapolation error for velocity can reach 1~2 mm/s, and the acceleration error can reach 10 −6 m/s 2 .By contrast, the extrapolation accuracy of orbit integration is high, and there is no significant increase in errors within 15 min for both satellite velocity and acceleration.Therefore, by utilizing orbit fitting, orbit integration, and orbit extrapolation, high-precision satellite velocities and accelerations can be obtained, albeit near the day boundary.
15 min extrapolation error for velocity can reach 1~2 mm/s, and the acceleration error can reach 10 −6 m/s 2 .By contrast, the extrapolation accuracy of orbit integration is high, and there is no significant increase in errors within 15 min for both satellite velocity and acceleration.Therefore, by utilizing orbit fitting, orbit integration, and orbit extrapolation, high-precision satellite velocities and accelerations can be obtained, albeit near the day boundary.Considering the high extrapolation accuracy of orbit fitting and orbit integration methods, we take the result as a reference to further demonstrate the effect of day boundary discontinuity on the satellite velocity and acceleration interpolated utilizing the Lagrange interpolation method.The experimental scheme is shown in Figure 6.As shown in the figure, the precise ephemeris of the adjacent day is combined for the Lagrange interpolation.Considering the high extrapolation accuracy of orbit fitting and orbit integration methods, we take the result as a reference to further demonstrate the effect of day boundary discontinuity on the satellite velocity and acceleration interpolated utilizing the Lagrange interpolation method.The experimental scheme is shown in Figure 6.As shown in the figure, the precise ephemeris of the adjacent day is combined for the Lagrange interpolation.

Conclusions
Undoubtedly, the GNSS precise ephemeris contains accurate satellite velocity and acceleration information.Users generally calculate satellite velocity and acceleration by interpolating satellite positions in the precise ephemeris.However, the interpolation method is easily affected by the edge effect, leading to a decrease in computational accuracy near the day boundary.To address the above issues, a calculation method for GNSS satellite velocity and acceleration based on orbit fitting and orbit integration was proposed.The satellite velocity and acceleration calculated utilizing orbit integration are in the ECI coordinate system.To obtain the satellite velocity and acceleration in the ECEF

Conclusions
Undoubtedly, the GNSS precise ephemeris contains accurate satellite velocity and acceleration information.Users generally calculate satellite velocity and acceleration by interpolating satellite positions in the precise ephemeris.However, the interpolation method is easily affected by the edge effect, leading to a decrease in computational accuracy near the day boundary.To address the above issues, a calculation method for GNSS satellite velocity and acceleration based on orbit fitting and orbit integration was proposed.The satellite velocity and acceleration calculated utilizing orbit integration are in the ECI coordinate system.To obtain the satellite velocity and acceleration in the ECEF coordinate system, a high-precision transformation method of ECI and ECEF for satellite velocity and acceleration was developed.
The accuracy of satellite velocity and acceleration calculated utilizing the method proposed in this study is relatively high.The satellite velocity RMS values for GEO, IGSO, and MEO are all less than 1.5 × 10 −6 m/s.The corresponding satellite acceleration RMS for GEO is less than 1.0 × 10 −8 m/s 2 , and the values for IGSO and MEO satellites are less than 0.5 × 10 −8 m/s 2 .In addition, the accuracy of extrapolating satellite velocity and acceleration by one-day orbit fitting arc is high, and there is basically no significant error growth within 15 min.Thus, the orbit fitting and orbit integration methods implement the high-precision calculation of GNSS satellite velocity and acceleration in a single day with uniform accuracy.Moreover, the above accuracy values can meet the requirements of high-precision satellite velocity and acceleration applications.
The method proposed in this study was validated utilizing GNSS ephemeris data.Its applicability to LEO satellites has not been verified yet.Compared with GNSS satellites, LEO satellites have higher operating speeds and more complex force models.Therefore, utilizing the proposed method to calculate the velocity and acceleration of LEO satellites will be a future research topics.

Figure 1 .
Figure 1.Flow chart of orbit fitting and orbit integration.

Figure 1 .
Figure 1.Flow chart of orbit fitting and orbit integration.
−ERA), where R 3 denotes the rotation matrix around the Z-axis, and ERA = 2π • (0.7790572732640 + 1.00273781191135448 • T u ), where T u is the simplified Julian day at UT1.Because T u = UTC + (UT1 − UTC), therefore .T u = 1 + UT1 • − UTC ; the first-orderderivatives of (UT1 − UTC) can be calculated by interpolating the corresponding Earth orientation parameters, and then, E .RA can be calculated;..

Figure 2 .
Figure 2. Satellite velocity accuracy calculated utilizing orbit fitting and orbit integration.The RMS values were calculated based on the differences between the result calculated by orbit fitting and orbit integration and that calculated by sliding window Lagrange interpolation.

Figure 2 .
Figure 2. Satellite velocity accuracy calculated utilizing orbit fitting and orbit integration.The RMS values were calculated based on the differences between the result calculated by orbit fitting and orbit integration and that calculated by sliding window Lagrange interpolation.

Figure 2 .
Figure2.Satellite velocity accuracy calculated utilizing orbit fitting and orbit integration.The RMS values were calculated based on the differences between the result calculated by orbit fitting and orbit integration and that calculated by sliding window Lagrange interpolation.

Figure 3 .
Figure 3. Satellite acceleration accuracy calculated utilizing orbit fitting and orbit integration.

Figure 3 .
Figure 3. Satellite acceleration accuracy calculated utilizing orbit fitting and orbit integration.

Figure 5 .
Figure 5. Three-dimensional RMS of satellite velocity and acceleration extrapolated utilizing orbit integration and Lagrange interpolation.The left side of the red dividing line represents fitting, and the right side represents extrapolation.

Figure 5 .
Figure 5. Three-dimensional RMS of satellite velocity and acceleration extrapolated utilizing orbit integration and Lagrange interpolation.The left side of the red dividing line represents fitting, and the right side represents extrapolation.

Figure 6 .
Figure 6.Experimental scheme for demonstrating the effect of day-boundary discontinuity.

Figure 7
Figure 7 presents the experimental results.It is clearly shown that the day-boundary discontinuity affects the accuracy of interpolating velocity and acceleration for GEO, IGSO, and MEO satellites.These results further demonstrate that utilizing the Lagrange interpolation method by means of combining the adjacent ephemeris cannot obtain highprecision satellite velocity and acceleration near the day boundary.

Figure 6 .
Figure 6.Experimental scheme for demonstrating the effect of day-boundary discontinuity.

Figure 7
Figure 7 presents the experimental results.It is clearly shown that the day-boundary discontinuity affects the accuracy of interpolating velocity and acceleration for GEO, IGSO, and MEO satellites.These results further demonstrate that utilizing the Lagrange interpolation method by means of combining the adjacent ephemeris cannot obtain high-precision satellite velocity and acceleration near the day boundary.

Figure 7
Figure7presents the experimental results.It is clearly shown that the day-boundary discontinuity affects the accuracy of interpolating velocity and acceleration for GEO, IGSO, and MEO satellites.These results further demonstrate that utilizing the Lagrange interpolation method by means of combining the adjacent ephemeris cannot obtain highprecision satellite velocity and acceleration near the day boundary.

Figure 7 .
Figure 7. Three-dimensional RMS of satellite velocity and acceleration interpolated utilizing the Lagrange method by means of combining adjacent ephemeris.The time ranges from 23:30 to 24:00 on 1 January 2022.

Figure 7 .
Figure 7. Three-dimensional RMS of satellite velocity and acceleration interpolated utilizing the Lagrange method by means of combining adjacent ephemeris.The time ranges from 23:30 to 24:00 on 1 January 2022.

Table 1 .
Orbit fitting and orbit integration force models.

Table 2 .
The multi-day average RMS values utilizing the approximate transformation model.The RMS values were calculated based on the differences between the result calculated utilizing orbit fitting and orbit integration and that calculated utilizing sliding window Lagrange interpolation.Note: the unit of acceleration is m/s 2 .

Table 3 .
The multi-day average RMS values utilizing the transformation model proposed in this study.

Table 4 .
The multi-day average RMS values when the first-and second-order derivatives of (UT1-UTC), Q(t), and W(t) are ignored in the transformation model proposed in this study.

Table 5 .
The multi-day average RMS values during non-eclipse seasons.Values before and after the slash represent the results of inclusion and exclusion of velocity pulses, respectively.

Table 6 .
The multi-day average RMS values during eclipse seasons.Values before and after the slash represent the results of inclusion and exclusion of velocity pulses, respectively.During the experimental data collection period in this study, no GEO satellites were utilized in the eclipse season.

Table 7 .
Satellite velocity and acceleration accuracy extrapolated utilizing Lagrange interpolation for 15 min.

Table 7 .
Satellite velocity and acceleration accuracy extrapolated utilizing Lagrange interpolation for 15 min.