Analysis of the movement model of a spacecraft in earth orbit

Full Text

Introduction

Satellite navigation (SN) has become an integral part of modern life. From mobile applications for city orientation to precise positioning in the maritime, aviation and space industries, SN provides reliable and accurate determination of location, speed and time [1]. These advantages are especially useful for autonomous navigation, automatic piloting, mapping and other applications that are just beginning to develop actively [2].

However, this technology is not perfect and has its own problems that can affect its efficiency and reliability:

1. signal interference: satellite signals transmitted by navigation systems can be subject to various interference such as multipath, noise and signalblocking. This can cause signal distortion and reduce the accuracy of positioning;
2. geographical limitations: in some places on Earth, such as deep valleys, mountainous areas or densely built-up cities, signals from satellites can be weakened or completely blocked. This leads to problems in obtaining a reliable navigation signal. A similar problem is observed in outerspace for spacecraft (SC) in geostationary orbit (GEO), when the available number of visible navigation spacecraft (VNS) is insufficient to accurately determine the state of the spacecraft;
3. dependence on the state of the navigation spacecraft: the SN is completely dependent on the functioning of global navigation satellite systems (GNSS) [3]. Any errors, technical failures, failures in satellite networks or deliberate acts of interference can lead to the unavailability or inadequate operation of the SN.

It should be noted that there are methods and technologies that help to cope with SN problems. For example, there are interference suppression methods that improve the quality of received signals, protective measures to reduce the vulnerability of the system, as well as alternative navigation systems, such as inertial ones. When used in combination with GNSS, the inertial system provides the consumer with coordinate and velocity data in the absence of GNSS signal reception.

Instantaneous (one-time) SWR solution for GNSS is a method that allows determining the position of an object (the state vector of the object – coordinates, speed and time) in real time. This method meets the requirements of most ground users, and together with the inertial system allows overcoming most of the problems of the SN.

Instantaneous SWR solution for objects in near-Earth orbits (NEO) is similar to the “ground” one, while the feature of orbital motion, which distinguishes it from the motion of ground objects, is high predictability – high accuracy of the motion forecast, determined through the parameters of gravity of the Earth, the Moon, the Sun, atmospheric resistance, light pressure, etc.

For navigation on the NEO, the forecast makes it possible to:

1. calculate in real time the motion at intervals of loss of tracking of the signals of the NSC (similar to the inertial navigation system), and also to forecast (predict) the motion in the future;
2. calculate from instantaneous coordinates and velocities (KS) the average orbit (SRO) on the interval, which can be more accurate than each of the initial instantaneous KS-solutions;
3. use recursive algorithms [4–7], allowing in real time to take into account various sources of errors, such as measurement errors, noise, interspatial interference, systematic errors in the coordinates of the navigation spacecraft.

Modern research in the field of modeling the orbital dynamics of spacecraft in near-Earth orbits is faced with the need to balance between the accuracy of forecasting and the computational complexity due to the specifics of applied problems. The papers [8–9] present simplified models aimed at short- term forecasting of trajectories of objects in low circular orbit (LCO) and medium-earth orbit (MEO), which is relevant for the tasks of pointing antennas with a wide coverage angle. In particular, [8] demonstrates the application of the SGP (Simplified General Perturbations) model for calculating the azimuth and elevation angle of the spacecraft basedon TLE (Two-LineElements) data, with an estimate of errors relative to a more accurate model. The study [9] proposes an algorithm for determining the position and orientation of the Earth relative to the spacecraft, which is critical for simulating docking processes. In contrast to the indicated works, the study [10] focuses on the trajectories of interplanetary missions (using the example of an expedition to the asteroid Apophis), where the set of disturbing factors includes the gravitational effects of remote celestial bodies, then on sphericity of the asteroid, and the pressure of sunlight. In the papers [11; 12] trajectory optimization methods are considered, including the pseudo-spectral Gaussian method and relativistic propagation models, but their computational complexity limits their applicability in real time for resource-constrainedsystems.

Complementing existing studies, a new approach to modeling orbital dynamics using machine learning is proposed in [13] to improve the accuracy of spacecraft trajectory prediction. The authors demonstrate that the integration of neural networks with classical orbital mechanics methods can significantly reduce the computational load while maintaining high prediction accuracy, but all processing occurs outside of real time.

It should be noted that most existing studies ignore GEO, limiting themselves to the navigation spacecraft and MEO. The objective of this paper is to comprehensively describe the spacecraft motion model implemented in the navigation modules produced by JSC KB NAVIS, with an emphasis on analyzing the impact of motion parameter errors on the prediction accuracy for various orbit classes. The results of numerical modeling of trajectory parameter deviations are presented, as well as flight test data confirming the need and possibility of taking into account dynamic disturbances in real time.

1. Physical model of forecasting the spacecraft state vector

WGS 84 and PZ-90 were chosen as coordinate systems for geodetic support of orbital flights and solving navigation problems. Both are characterized by the fact that the origin of the system is located at the center of mass of the Earth. The coordinate axis Z, in accordance with the recommendations of the International Earth Rotation Service (IERS), is directed toward the middle North Pole. The coordinate axis X lies in the plane of the earth’s equator of the same epoch, forming an intersection with the plane of the prime meridian established by the same IERS, and determines the position of the zero point of the adopted counting system. The Y axis complements the coordinate system to the right.

The difference between the two coordinate systems consists in the values​​ of the set of constants. WGS 84 is mainly used, so an example is given for this coordinate system:

• angular velocity of Earth rotation ${\mathrm{\omega }}_{\text{з}}$ = 7.292115×10^–5 rad/s;
• Earth radius $R_{\text{з}}$ (semi-major axis $a_{\text{e}}$) = 6378137 m;
• gravitational constant µ = 398600.4418×10⁹ m³/s²;
• eccentricity of the Earth’s orbit 6.69437999013×10^–3.

It is assumed that the spacecraft is not subject to control effects, i.e. it does not use control engines and has a constant mass. Thus, the position and velocity of the spacecraft at any given moment will be affected by the following components:

1. the force of attraction of a spherical and homogeneous model of the Earth;
2. the centrifugal force and the Coriolis force caused by the rotation of the Earth and the motion of the spacecraft;
3. the anomalous force of attraction of the Earth caused by the difference between the Earth and a spherical and homogeneous mass [14];
4. atmospheric resistance;
5. the attraction of the Sun;
6. the attraction of the Moon;
7. the pressure of sunlight.

2. Mathematical formulation of the problem of forecasting the state vector of the spacecraft

2.1. General system of equations

Forecasting the parameters of the spacecraft motion on the orbit from the moment of time t0 to the moment t is performed by numerically integrating the following differential equations of motion in a rectangular geocentric coordinate system (RCS) to determine the position coordinates (x, y, z) and velocity components (V_x, V_y, V_z)

$\frac{dx}{dt}=V_x,$ (1)

$\frac{dy}{dt}=V_y,$ (2)

$\frac{dz}{dt}=V_z,$ (3)

$\frac{d{V}_x}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot x+{\omega }_3^2\cdot x+2\cdot {\omega }_3\cdot V_y+{\ddot{x}}_a+{\ddot{x}}_s+{\ddot{x}}_c+{\ddot{x}}_m+{\ddot{x}}_l,$ (4)

$\frac{d{V}_y}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot y+{\omega }_3^2\cdot y+2\cdot {\omega }_3\cdot V_x+{\ddot{y}}_a+{\ddot{y}}_s+{\ddot{y}}_c+{\ddot{y}}_m+{\ddot{y}}_l,$ (5)

$\frac{d{V}_z}{dt}=-\frac{\mathrm{\mu }}{r^3}\cdot z+{\ddot{z}}_a+{\ddot{z}}_s+{\ddot{z}}_c+{\ddot{z}}_m+{\ddot{z}}_l,$ (6)

where $r=\sqrt{x^2+y^2+z^2}$ is the radius vector of the spacecraft position; ${\ddot{x}}_a,{\ddot{y}}_a,{\ddot{z}}_a$ are the components of the acceleration vector caused by the difference between the Earth model and a spherical one; ${\ddot{x}}_s,{\ddot{y}}_s,{\ddot{z}}_s$ are the components of the acceleration vector caused by aerodynamic braking of the spacecraft in the Earth’s atmosphere; ${\ddot{x}}_c,{\ddot{y}}_c,{\ddot{z}}_c$ are the components of the acceleration vector caused by the attraction of the Sun; ${\ddot{x}}_m,{\ddot{y}}_m,{\ddot{z}}_m$ are the components of the acceleration vector caused by the attraction of the Moon; ${\ddot{x}}_l,{\ddot{y}}_l,{\ddot{z}}_l$ are the components of the acceleration vector due to the pressure of sunlight.

The force of attraction of a spherical and homogeneous model of the Earth corresponds to the first term in equations (4)–(6). The centrifugal force corresponds to the second, and the Coriolis force to the third term of formulas (4)–(5).

2.2. Determination of spacecraft accelerations caused by geopotential

The components of the spacecraft acceleration vector caused by the difference between the Earth model and a spherical one are first determined in the topocentric coordinate system (TCS) through the geopotential coefficients [14], after which they are transformed into RCS components ${\ddot{x}}_a,{\ddot{y}}_a,{\ddot{z}}_a$:

$\left(\begin{array}{c}{\dot{x}}_a\\ {\dot{y}}_a\\ {\dot{z}}_a\end{array}\right)=\left(\begin{array}{ccc}\frac{-x}{r}& -\frac{x\cdot z}{r\cdot r_1}& -\frac{y}{r_1}\\ -\frac{y}{r}& \frac{y\cdot z}{r\cdot r_1}& \frac{x}{r_1}\\ -\frac{z}{r}& \frac{r_1}{r_1}& 0\end{array}\right)\times \left(\begin{array}{c}{\dot{U}}_r\\ {\dot{U}}_{\varphi }\\ {\dot{U}}_{\lambda }\end{array}\right)$, (7)

where ${\dot{U}}_r,{\dot{U}}_{\varphi },{\dot{U}}_{\lambda }.$ are the components of the acceleration vector in the TCS: by radius vector r, by latitude ϕ and longitude λ; $r_1=\sqrt{x^2+y^2}$.

The components of the acceleration vector ${\ddot{U}}_r,{\ddot{U}}_{\varphi },{\ddot{U}}_{\lambda }$ are determined as the sums of the trigonometric series:

$\begin{array}{l}{\ddot{U}}_r=\mathrm{\mu }/r^2\sum _{n=2}^N{(a_e/r)}^n(n+1)\sum _{m=0}^n{A}_{nm}{\overline{P}}_{nm},\\ {\ddot{U}}_{\varphi }=\mathrm{\mu }/r^2\sum _{n=2}^N{(a_e/r)}^n\sum _{m=0}^n{A}_{nm}\cdot {{\overline{P}}^{\text{'}}}_{nm},\\ {\ddot{U}}_{\lambda }=\frac{\mathrm{\mu }}{r^2\cdot \cos \varphi }\sum _{n=2}^N{(a_e/r)}^n\sum _{m=0}^{n}m{B}_{nm}{\overline{P}}_{nm},\end{array}$

where

$\begin{array}{l}\mathrm{sin\varphi }=z/r,\\ \mathrm{cos\varphi }=r_1/r,\\ A_{nm}={\overline{C}}_{nm}\cdot \cos (m\cdot \lambda )+{\overline{S}}_{nm}\cdot \sin (m\cdot \lambda ),\\ B_{nm}={\overline{S}}_{nm}\cdot \cos (m\cdot \lambda )-{\overline{C}}_{nm}\cdot \sin (m\cdot \lambda ).\end{array}$

${\overline{C}}_{nm},{\overline{S}}_{nm}$ are the normalized coefficients of the expansion of the Earth’s potential.

The trigonometric functions of angles that are multiples of l are calculated using the formulas

$\begin{array}{l}\sin (m\cdot \lambda )=\sin \left(\lambda \right)\cdot \cos \left(\right[m-1]\cdot \lambda )+\cos \left(\lambda \right)\cdot \sin \left(\right[m-1]\cdot \lambda ),\\ \cos (m\cdot \lambda )=\cos \left(\lambda \right)\cdot \cos \left(\right[m-1]\cdot \lambda )-\sin \left(\lambda \right)\cdot \sin \left(\right[m-1]\cdot \lambda ),\\ \cos \lambda =x/r_1,\sin \lambda =y/r_1.\end{array}$

The normalized Legendre polynomials ${\overline{P}}_{nm}$ are functions of $\sin \left(\mathrm{\varphi }\right)$, $\cos \left(\mathrm{\varphi }\right)$ and are determined using the following recurrent dependencies:

$\begin{array}{l}{\overline{P}}_{nm}=\sqrt{\frac{2\cdot n+1}{2\cdot n}}\cdot \cos \left(\varphi \right)\cdot {\overline{P}}_{n-\mathrm{1,}m-1},\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}n}=\text{m}>\text{1,}\\ {}_{\text{nm}}=\sqrt{\frac{2\cdot n+1}{n^2-m^2}}[\sqrt{2\cdot n-1}\cdot \sin (\varphi )\cdot {\overline{P}}_{n-\mathrm{1,}m}-\sqrt{\frac{(n-m-1)\cdot (n+m-1)}{2\cdot n-3}}\cdot {\overline{P}}_{n-\mathrm{2,}m}],\text{\hspace{1em}n}>\text{m,}\\ {}_{\text{nm}}=0,\text{\hspace{1em}m}>\text{n\hspace{0.33em}.}\end{array}$

Initial values: ${\overline{\mathrm{P}}}_{00}=\mathrm{1,}{\overline{\mathrm{P}}}_{10}=\sqrt{3}\mathrm{sin\varphi },{\mathrm{P}}_{11}=\sqrt{3}\mathrm{cos\varphi }$.

The derivatives of the normalized Legendre polynomials ${\overline{P}}_{nm}$ with respect to are determined by the following recurrent dependencies:

$\begin{array}{l}{{\overline{P}}^{\text{'}}}_{nm}=\sqrt{\frac{2\cdot n+1}{2\cdot n}}\cdot \left(\cos \right(\varphi )\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{1,}m-1}-\sin (\varphi )\cdot {\overline{P}}_{n-\mathrm{1,}m-1}),\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}n}=\text{m}>\text{1,}\\ {}_{\text{nm}}=\sqrt{\frac{2\cdot n+1}{n^2-m^2}}\cdot \left[\sqrt{2\cdot n-1}\right(\sin \left(\varphi \right)\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{1,}m}+\cos \left(\varphi \right)\cdot {\overline{P}}_{n-\mathrm{1,}m})-\sqrt{\frac{(n-m-1)\cdot (n+m-1)}{2\cdot n-3}}\cdot {{\overline{P}}^{\text{'}}}_{n-\mathrm{2,}m}],\text{\hspace{1em}n}>\text{m,}\\ {}_{\text{nm}}=0,\text{\hspace{1em}m}>\text{n\hspace{0.33em}.}\end{array}$

Initial values: ${{\overline{\mathrm{P}}}^{\text{'}}}_{00}=\mathrm{0,}{{\overline{\mathrm{P}}}^{\text{'}}}_{10}=\sqrt{3}\cdot \cos \left(\mathrm{\varphi }\right),{{\mathrm{P}}^{\text{'}}}_{11}=-\sqrt{3}\cdot \sin \left(\mathrm{\varphi }\right).$

2.3. Determining spacecraft accelerations due to atmospheric braking

The components of the acceleration vector ${\ddot{x}}_s,{\ddot{y}}_s,{\ddot{z}}_s$ caused by aerodynamic braking of the spacecraft in the Earth’s atmosphere [14] are calculated using the formulas

$\begin{array}{l}{}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{x}}\text{, }\\ {}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{y}}\text{,}\\ {}_{\text{s}}=-{\text{S}}_{\text{b}}\cdot \mathrm{\rho }\cdot \text{V}\cdot {\text{V}}_{\text{z}}\text{,}\end{array}$ (8)

where ρ is the density of the upper atmosphere from the national standard model (GOST R 25645.166–2004); $V=\sqrt{{V_x}^2+{V_y}^2+{V_z}^2}$ is the magnitude of the spacecraft velocity vector; ${\text{S}}_{\text{b}}=C_x{S}_m/2m$ is the ballistic coefficient; m is the spacecraft mass; $C_x$ is the aerodynamic drag coefficient; $S_m$ is the spacecraft midsection area.

The values ​​of the aerodynamic coefficient $C_x$ and the midsection area $S_m$ are determined by the size and shape of a specific spacecraft.

2.4. Determining spacecraft accelerations due to solar gravitation

The components of the acceleration vector ${\overrightarrow{\ddot{r}}}_c$ = (${\ddot{x}}_c$, ${\ddot{y}}_c$, ${\ddot{z}}_c$), caused by disturbances from the Sun [14], are determined using the following formulas:

${\overrightarrow{\ddot{\mathrm{r}}}}_{\mathrm{c}}=\frac{{\mathrm{\mu }}_{\mathrm{c}}\cdot ({\overrightarrow{\mathrm{r}}}_{\mathrm{c}}-\overrightarrow{\mathrm{r}})}{{\left|{\overrightarrow{\mathrm{r}}}_{\mathrm{c}}-\overrightarrow{\mathrm{r}}\right|}^3}-\frac{{\mathrm{\mu }}_{\mathrm{c}}\cdot {\overrightarrow{\mathrm{r}}}_{\mathrm{c}}}{{\left|{\overrightarrow{\mathrm{r}}}_{\mathrm{c}}\right|}^3}$,

where ${\mathrm{\mu }}_c$ is the product of the gravitational constant and the mass of the Sun; ${\overrightarrow{r}}_c$ is the vector of the Sun’s position relative to the Earth; $\overrightarrow{r}$ is the vector of the object’s position relative to the Earth.

The calculation of the ecliptic and geocentric coordinates of the Sun [15] uses the mean elements of the Sun’s orbit as constants, determined for the J2000 epoch. The estimated calculation error is no more than 10 arc minutes.

The rectangular equatorial coordinates of the Sun are calculated using the radius vector, longitude, and inclination of the ecliptic to the equator. When calculating the coordinates of the Sun on the current date, its elliptical motion is taken into account. The radius vector of the Sun and the longitude of the Sun on the current date are determined through series by the mean anomaly. The mean anomaly on the current date is determined from the mean anomaly for the J2000 epoch, taking into account the mean motion for the J2000 epoch and the time in Julian centuries from the J2000 epoch. The inclination of the ecliptic to the equator on the current date is determined from the inclination of the ecliptic at J2000 and the change ininclination over time in Julian centuries from the J2000 epoch.

The coordinates of the Sun in the rectangular geocentric coordinate system on the current date are determined from the equatorial coordinates by rotating with the Greenwich hour angle. The Greenwich hour angle on the current date is determined without taking into account the movement of the poles, precession and nutation from the Greenwich hour angle at the J2000 epoch, taking into account its change for the time after the J2000 epoch.

The equatorial coordinates of the Sun are determined by the following formulas:

$\begin{array}{c}x=r\cdot \text{cos}\left(L\right),\\ y=r\cdot \text{sin}\left(L\right)\cdot \text{cos}\left(EPS\right),\\ z=r\cdot \text{sin}\left(L\right)\cdot \text{sin}\left(EPS\right),\end{array}$

where r is the distance to the Sun:

$r=[\mathrm{149,619}-\mathrm{2,499}\cdot \text{cos}(M)-\mathrm{0,021}\cdot \text{cos}(2\cdot M\left)\right]\cdot {10}^9,$

M is the mean anomaly:

$M=357°\mathrm{,528}+35999°\mathrm{,05}\times T$,

L is the longitude of the Sun taking into account precession:

$L=282°\mathrm{,94}+M+[6892"\cdot \text{sin}(M)+72"\times \text{sin}(2\cdot M\left)\right]/3600+1°3972\cdot T$,

EPS is the inclination of the ecliptic to the equator:

$EPS=23°\mathrm{,43929111}–46"\mathrm{,8150}\cdot T–0"\mathrm{,00059}\cdot T^2+0"\mathrm{,001813}\cdot T^3$,

T is the time in Julian centuries from the J2000 epoch.

The coordinates of the Sun (the Sun vector ${\overline{R}}_{\text{S}}=\left[x_S{y}_S{z}_S\right]$) in the rectangular geocentric coordinate system are determined by the formulas

$\begin{array}{c}{x}_S=x\cdot \text{cos}\left(Tgr\right)+y\cdot \text{sin}\left(Tgr\right),\\ y_S=y\cdot \text{cos}\left(Tgr\right)–x\cdot \text{sin}\left(Tgr\right),\\ z_S=z,\end{array}$

where $Tgr$ is the Greenwich hour angle from the J2000 epoch:

$Tgr={280}^{\circ }\mathrm{,46061837504}+{360}^{\circ }\mathrm{,9856473662862}\cdot d$,

$i=Tgr/360$,

$RAD=180/\pi$,

$Tgr=(Tgr–i\cdot 360)/RAD$, if $Tgr<0$, then $Tgr=Tgr+2\pi$,

$d=\left(MJD–51544\right)+(FMJD\cdot 86400–dUTC)/86400–\mathrm{0,5};$

MJD is the number of days from the beginning of the J2000 epoch,

FMJD is the number of seconds from the beginning of the epoch,

dUTC is the time correction to the UTC scale.

The calculation of the equatorial coordinates of the Sun was verified using the Astronomical Year book of the USSR for 1983 and control points from the NASA website. Examples of the calculation are presented in Table 1.

 

Table 1

Evaluation of the error in calculating the coordinates of the Sun

Calculation of the equatorial coordinates of the Sun on 21.11.2010 12:00:00
Equatorial coordinates of the Sun	X, km	Y, km	Z, km
NASA website	–7.6318610 E+07	–1.1610898 E+08	–5.0335183 E+07
Calculation by formulas	–7.6038571 E+07	–1.1626342 E+08	–5.0402949 E+07
Calculation error	0.0280039 E+07	0.0015444 E+08	0.0067766 E+07
Max error	0.3269033 E+06 km=0.00218=7.5 arc minutes
Calculation of the equatorial coordinates of the Sun on 21.06.2011 12:00:00
Equatorial coordinates of the Sun	X, km	Y, km	Z, km
NASA website	9.9620932 E+05	1.3947997 E+08	6.0467823 E+07
Calculation by formulas	6.4872767 E+05	1.3948265 E+08	6.0468820 E+07
Calculation error	3.4748165 E+05	0.0000268 E+08	0.0000997 E+07
Max error	0.347493 E+06 km = 0.00232 = 8.0 arc minutes

 

2.5. Determining the spacecraft accelerations due to the attraction of the Moon

The components of the acceleration vector ${\overrightarrow{\ddot{r}}}_m$ = (${\ddot{x}}_m$, ${\ddot{y}}_m$, ${\ddot{z}}_m$), caused by disturbances from the Moon [14], similar to the Sun, are determined by the following formulas:

${\overrightarrow{\ddot{r}}}_m=\frac{{\mathrm{\mu }}_m\cdot ({\overrightarrow{r}}_m-\overrightarrow{r})}{{\left|{\overrightarrow{r}}_m-\overrightarrow{r}\right|}^3}-\frac{{\mathrm{\mu }}_m\cdot {\overrightarrow{r}}_m}{{\left|{\overrightarrow{r}}_m\right|}^3}$,

where ${\mathrm{\mu }}_m$ is the product of the gravitational constant and the mass of the Moon; ${\overrightarrow{r}}_m$ = ($x_m$, $y_m$, $z_m$) is the vector of the Moon’s position relative to the Earth; $\overrightarrow{r}$ is the vector of the spacecraft’s position relative to the Earth.

To calculate the coordinates of the Moon [15], it is necessary to know the current time in the format of year, month, day, hour, minute, second, correction to the UTC (US) time zone.

The calculation procedure and formulas are as follows:

conversion of the current date and time to time in Julian centuries T from the epoch J2000 similar to the Sun;

calculation of the average longitude of the Moon el0 (in degrees):

$el0=218.31617+481267.88088\cdot T–1.3972\cdot T^2$;

average anomaly of the Moon el (indegrees):

$el=134.96292+477198/86753\cdot T$;

average anomaly of the Sun elp (in degrees):

$elp=357.528+35999.04944\cdot T$;

mean angular distance f of the Moon from the ascending node (in degrees):

$f=93.27283+483202.01873\cdot T;$

difference between mean longitudes of the Sun and Moon (in degrees):

$d=297.85027+445267.11135\cdot T$;

true ecliptic longitude of the Moon (epoch J2000) differs from mean longitude by periodic terms Lon and dlon (in degrees):

$Lon=el0+dlon$

$\begin{array}{l}dlon=[22640\cdot \text{sin}\left(el\right)\text{ + }769\cdot \text{sin}\left(2\cdot el\right)–4586\cdot \text{sin}\left(el–2\cdot d\right)+\\ +2370\cdot \text{sin}\left(2\cdot d\right)–668\cdot \text{sin}\left(elp\right)–412\cdot \text{sin}\left(2\cdot f\right)–212\cdot \text{sin}\left(2\cdot el–2\cdot d\right)\\ –206\cdot \text{sin}\left(el+elp–2\cdot d\right)+192\cdot \text{sin}\left(el+2\cdot d\right)-165\cdot \text{sin}\left(elp–2\cdot d\right)+148\cdot \text{sin}\left(el-elp\right)\\ –125\cdot \text{sin}\left(d\right)–110\cdot \text{sin}\left(el+elp\right)-55\cdot \text{sin}\left(2\cdot f–2\cdot d\right)]/3600\end{array}$;

true ecliptic latitude of the Moon Lat (epoch J2000) (in degrees):

$\begin{array}{l}Lat=[18520\cdot \text{sin}\left(f+Lon–el0+q\right)–526\cdot \text{sin}\left(f–2\cdot d\right)+44\cdot \text{sin}\left(el+f–2\cdot d\right)–\\ –31\cdot \text{sin}\left(–el+f–2\cdot d\right)–25\cdot \text{sin}\left(–2\cdot el+f\right)–23\cdot \text{sin}\left(elp+f–2\cdot d\right)\\ +21\cdot \text{sin}\left(–el+f\right) +11\cdot \text{sin}\left(–elp+f–2\cdot d\right)]/3600,\end{array}$

where $\begin{array}{l}q=412/3600\cdot \text{sin}\left(2\cdot f\right)+541/3600\cdot \text{sin}\left(elp\right)\end{array}$;

distance Rse from the center of the Earth to the Moon (in meters):

$\begin{array}{l}Rse=[38500–20905\cdot \text{cos}\left(el\right)–3699\cdot \text{cos}\left(2\cdot d–el\right)–2956\cdot \text{cos}\left(2\cdot d\right)–\\ -570\cdot \text{cos}\left(2\cdot el\right)+46\times \text{cos}\left(2el–2d\right)–20\times 5\text{cos}\left(elp–2d\right)-\\ -171\times \text{cos}\left(el+2d\right)+2–152\times \text{cos}\left(el+elp–2d\right)]\times 1000;\end{array}$

true ecliptic longitude of the Moon taking into account precession (indegrees):

$Lon=Lon+1.3972\times T ;$

inclination of the ecliptic to the equator Obe (in degrees):

$Obe=23.43929111 ;$

ecliptic coordinates of the Moon xse, yse, zse:

$\begin{array}{l}xe=Rse\cdot \text{cos}\left(Lon\right)\cdot \text{cos}\left(Lat\right),\\ ye=Rse\cdot \text{sin}\left(Lon\right)\cdot \text{cos}\left(Lat\right),\\ ze=Rse\cdot \text{sin}\left(Lat\right);\end{array}$

equatorial coordinates of the Moon xse, yse, zse:

$\begin{array}{l}xse=xe\\ yse=ye\cdot \text{cos}\left(Obe\right)–ze\cdot \text{sin}\left(Obe\right), \\ zse=ze\cdot \text{cos }\left(Obe\right)+ye\cdot \text{sin}\left(Obe\right);\end{array}$

coordinates of the Moon in the rectangular geocentric coordinate system:

$\begin{array}{l}{x}_m=xse\cdot \text{cos(}Tgr)+yse\cdot \text{sin}(Tgr),\\ y_m=yse\cdot \text{cos}\left(Tgr\right)–xse\cdot \text{sin}\left(Tgr\right),\\ z_m=zse,\end{array}$

where the Greenwich hour angle Tgr from the J2000 epoch is calculated using formulas similar to the Sun.

2.6. Determining the spacecraft accelerations due to solar radiation pressure

To calculate the spacecraft accelerations due to the action flight pressure [14], it is necessary to know the effective reflection coefficient Crefl, which depends on the cross-sectional area, mass and other properties of the surface. The components of the acceleration vector ${\overrightarrow{\ddot{r}}}_l$ = (${\ddot{x}}_l$, ${\ddot{y}}_l$, ${\ddot{z}}_l$), caused by solar radiation pressure are determined by the following formula:

${\overrightarrow{\ddot{r}}}_l=Crefl\cdot {\left(\frac{1.4959787\cdot {10}^8}{\left|{\overrightarrow{r}}_c\right|}\right)}^2\cdot \frac{(\overrightarrow{r}-{\overrightarrow{r}}_c)}{\left|{\overrightarrow{r}}_c\right|}$, (9)

where $\left|{\overrightarrow{r}}_c\right|$ is the distance from the Sun to the spacecraft.

2.7. The fourth-order Runge-Kutta method

The numerical integration of the system of ordinary nonlinear differential equations of spacecraft motion (1)–(6) of the form ${y^{\text{'}}}_j=f_j(t,y_1,y_2,y_3,y_4,y_5,y_6)$ is carried out by the fourth-order Runge-Kutta method with a constant step [16] (j = 1,…, 6 is the number of equations describing the spacecraft motion).

For the k-th step and the j-th equation, it denote

$y_{jk}=y_j\left(t_k\right),f_{jk}=f_j(t_k,{y_1}_k,{y_2}_k,{y_3}_k,{y_4}_k,{y_5}_k,{y_6}_k)$.

At the k+1 integration step, the value of the sought functions is calculated using the formula:

$y_{jk+1}=y_{jk}+h/6\cdot (k_{1j}+2\cdot k_{2j}+2\cdot k_{3j}+k_{4j})$,

where h is the integration step,

$\begin{array}{l}{k}_{1j}=f_j(t_k,{y_1}_k,{y_2}_k,{y_3}_k,{y_4}_k,{y_5}_k,{y_6}_k),\\ k_{2j}=f_j(t_k+\frac{h}{2},{y_1}_k+\frac{h}{2}\cdot k_{11},{y_2}_k+\frac{h}{2}\cdot k_{12},{y_3}_k+\frac{h}{2}\cdot k_{13},{y_4}_k+\frac{h}{2}\cdot k_{14},{y_5}_k+\frac{h}{2}\cdot k_{15},{y_6}_k+\frac{h}{2}\cdot k_{16}),\\ k_{3j}=f_j(t_k+\frac{h}{2},{y_1}_k+\frac{h}{2}\cdot k_{21},{y_2}_k+\frac{h}{2}\cdot k_{22},{y_3}_k+\frac{h}{2}\cdot k_{23},{y_4}_k+\frac{h}{2}\cdot k_{24},{y_5}_k+\frac{h}{2}\cdot k_{25},{y_6}_k+\frac{h}{2}{k}_{26}),\\ k_{4j}=f_j(t_k+h,{y_1}_k+h\cdot k_{31},{y_2}_k+h\cdot k_{32},{y_3}_k+h\cdot k_{33},{y_4}_k+h\cdot k_{34},{y_5}_k+h\cdot k_{35},{y_6}_k+h\cdot k_{36}).\end{array}$

With numerical integration over a large number of steps, there is a significant loss of accuracy due to the accumulation of rounding errors. It is believed that the accumulation of rounding errors in coordinates is proportional to the number of integration steps raised to the power of 3/2:

${\epsilon }_n^i\approx 3n^{3/2}\cdot \epsilon \left(0\right)$,

where $\epsilon \left(0\right)$ is the calculation accuracy at each step.

3. Errors in predicting the spacecraft motion along the NEO in the presence of errors

The above prediction equations were programmed in the MATLAB environment and the calculations were performed:

• for NEO with a period of revolution around the Earth of 1.5 hours and an orbital altitude above sea level of 500 km;
• SRO with a period of 12 hours and an altitude of 20,000 km;
• GEO with a period of 24 hours and a naltitude of 35,777 km.

The results of the prediction from the harmonics of the Earth’s potential to the 64^th order, taking into account accelerations from the Sun and the Moon, with a zero ballistic coefficient and a zero effective reflection coefficient, were used as reference data.

The motion parameters (including the starting point of the prediction) were varied, and the prediction error was estimated. The set of data obtained allows, in the first approximation, based on the permissible prediction error, to estimate the accuracy with which the parameters used in the prediction should be known.

3.1. Forecast errors from the error in the initial position

In order to estimate the effect of the error in coordinates and velocities along different axes, a uniform distribution of the displacement of the initial state vector with the specified modulus was specified separately for coordinates (Table 2) and velocities (Table 3) and then the maximum result of the forecast error was selected.

 

Table 2

Forecast error from errors in the initial coordinates

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Error modulus, m	Forcast interval, h
	0.5	1.0	1.5			3.0		6.0
0.5	2.1	5.5	8.0			15.7		30.8
2	9.7	27.3	36.6			73.1		145.6
10	39.2	101.7	144.6			286.5		566.9
30	118.5	267.2	301.0			619.4		1256.9
MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Error modulus, m	Forcast interval, h
	0.5	1.0		1.5	3.0		6.0	12.0	24.0
0.5	0.5	0.5		0.7	1.8		5.6	8.7	17.7
2	2.0	2.2		2.5	6.0		18.2	38.6	76.3
10	10.4	12.1		15.7	38.3		120.9	202.6	407.0
30	32.2	39.4		51.8	122.9		399.0	730.2	1461.4
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
Error modulus, m	Forcast interval, h
	0.5	1.0		1.5	3.0		6.0	12.0	24.0
0.5	0.5	0.5		0.6	0.9		2.6	10.0	18.7
2	2.0	2.0		2.2	2.9		7.3	29.4	55.9
10	10.1	10.6		11.3	16.3		43.3	166.2	310.6
30	30.6	32.4		35.5	53.9		145.1	555.2	1036.8

 

Table 3

Forecast error from initial errors in velocities

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Error modulus, m/s:	Forcast interval, h
	0.5		1.0			1.5		3.0		6.0
0.0001	0.3		1.4			1.7		3.4		6.8
0.0005	1.9		6.9			8.6		17.0		33.6
0.0020	7.7		28.1			34.7		68.7		135.8
0.0050	19.8		67.2			86.1		170.2		336.3
MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Error modulus, m/s:	Forcast interval, h
	0.5	1.0		1.5	3.0		6.0		12.0		24.0
0.0001	0.1	0.3		0.6	1.4		6.7		12.4		24.5
0.0005	0.9	1.9		3.1	9.3		34.5		65.7		130.1
0.0020	3.6	7.7		12.7	34.3		131.4		261.4		517.3
0.0050	9.1	19.2		30.6	94.5		352.4		621.9		1232.9
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
Error modulus, m/s:	Forcast interval, h
	0.5	1.0		1.5	3.0		6.0		12.0		24.0
0.0001	0.1	0.3		0.5	1.2		3.6		14.7		25.8
0.0005	0.9	1.8		2.8	6.4		19.6		74.6		129.4
0.0020	3.6	7.3		11.1	24.9		76.5		268.1		507.2
0.0050	9.0	18.3		28.3	64.9		180.2		625.7		1277.6

 

In real conditions, position and velocity errors act in a complex manner. Based on the initial errors and the required forecast interval, the data below can be used to estimate the need to take into account other forecast parameters and their accuracy.

3.2. Forecast errors from different numbers of considered geopotential harmonics

The anomalous force of gravity of the Earth (7) depends on the number of considered geopotential harmonics. The error in calculating coordinates taking into account different numbers of harmonics of the Earth’s potential from the initial coordinates and velocities is presented in Table 4.

 

Table 4

Forecast errors from the number of considered geopotential harmonics

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Number of harmonics:	Forcast interval, h
	0.5		1.0		1.5		3.0			6.0
32	1.1		1.4		1.4		4.4			14.4
16	5.3		14.8		24.4		19.0			25.2
8	20.8		65.1		66.7		140.2			197.2
2	100.2		236.6		236.6		451.9			626.1
MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Number of harmonics:	Forcast interval, h
	0.5	1.0		1.5		3.0		6.0	12.0		24.0
32	0.0	0.0		0.0		0.0		0.0	0.0		0.0
16	0.0	0.0		0.0		0.0		0.0	0.0		0.0
8	0.0	0.0		0.0		0.0		0.1	0.3		0.8
2	0.1	0.7		1.5		5.6		12.4	27.6		32.0
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
Number of harmonics:	Forcast interval, h
	0.5	1.0		1.5		3.0		6.0	12.0		24.0
32	0.0	0.0		0.0		0.0		0.0	0.0		0.0
16	0.0	0.0		0.0		0.0		0.0	0.0		0.0
8	0.0	0.0		0.0		0.0		0.0	0.0		0.1
2	0.1	0.1		0.1		0.1		0.6	3.3		21.3

 

To achieve minimal forecast errors for NKO, it is necessary to take into account the largest number of harmonics of the Earth’s potential. With increasing orbital altitude, the required number of harmonics decreases, but 8 harmonics is the necessary minimum.

3.3. Atmospheric drag influence

The atmospheric state was taken into account according to GOST R 25645.166–2004. The variation of the ballistic coefficient Sb in (8) has a stronger effect in low orbit [17] due to the denser atmosphere. Therefore, several NKOs were additionally considered (Table 5). Sb is linearly related to acceleration, which means thati n (1)–(6) the coordinate error will have the same dependence, despite the orbit nonlinearity, which is shown for one NKO.

 

Table 5

Atmospheric drag influence

Sb	Forcast interval, h
	0.5	1.0	1.5	3.0	6.0
NKO. Orbital period: Т = 1.45 ч. Altitude 200 km. Forecast error, m:
6.25е-05	1.3	11.8	36.6	143.5	570.6
6.25е-04	13.6	118.9	366.5	1437,7	5722.7
6.25е-03	136.2	1191.2	3683.4	>10000	>10000
6.25е-02	1364.6	>10000	>10000	>10000	>10000
NKO. Orbital period: Т = 1.49 h. Altitude 288 km. Forecast error, m:
6.25е-03	7.8	67.8	190.0	766.8	3098.1
NKO. Orbital period: Т = 1.50 h. Altitude 500 km. Forecast error, m:
6.25е-03	0.1	0.8	3.1	12.2	48.9
NKO. Orbital period: Т = 1.68 h. Altitude 600 km. Forecast error, m:
6.25е-03	0.1	0.2	0.7	2.7	10.9
NKO.Orbital period: Т = 1.76 h. Altitude 1000 km. Forecast error, m:
6.25е-03	0.0	0.0	0.1	0.1	0.4

 

The weakening of the influence of the ballistic coefficient with increasing orbital altitude is shown in Table 5. It is also noticeable that an increase in the ballistic coefficient by an order of magnitude leads to an increase in the forecast error by an order of magnitude.

Due to the absence of an atmosphere above 3,000 km, there is no effect on the forecast error for MEO and GEO. However, for NKO taking into account the ballistic coefficient is necessary, and for elliptical orbits combining the properties of NKO and MEO, it is necessary within the specified altitude.

3.4. Forecast errors due to the Sun’s gravity

In assessing the influence of the Sun’s gravity on the forecast of the spacecraft position, two parameters are considered (Table 6):

• the error in determining the angular position of the Sun. In order to set the angular error in the coordinates of the Sun (the error in the influence of gravity), the coordinates were shifted by a uniformly distributed random value so that the shift angle was equal to the specified value, then the maximum result of the forecast error was selected among equal values ​​of the shift angle;
• the error in determining the distance to the Sun due to the periodic (yearly) change in the radius of the Earth’s orbit from 147.098 to 152.098 thousand km.

The forecast error is also given in the case when the influence of the Sun’s gravity is not taken into account.

It follows from the data in Table 6 that the variation in the object-Sun distance due to the position of the object in orbit can practically be ignored.

 

Table 6

Ошибки прогноза от притяжения Солнца

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Error parameter		Forcast interval, h
		0.5	1.0	1.5	3.0	6.0
Angle, degrees	0.1	0.0	0.0	0.0	0.0	0.1
	0.5	0.0	0.1	0.1	0.1	0.1
	2	0.1	0.2	0.2	0.2	0.3
	10	0.5	1.0	1.0	1.4	2.5
Distance, thousand km	20	0.0	0.0	0.0	0.0	0.0
	100	0.0	0.0	0.0	0.0	0.0
	500	0.0	0.0	0.0	0.1	0.1
	2.500	0.0	0.1	0.1	0.2	0.4
Excluding ${\overrightarrow{\ddot{r}}}_c$		0.6	2.2	2.3	4.1	7.6

 

End of Table 6

MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Error parameter		Forcast interval, h
		0.5	1.0	1.5	3.0	6.0	12.0	24.0
Angle, degrees	0.1	0.0	0.0	0.1	0.2	1.0	2.2	2.6
	0.5	0.1	0.2	0.5	1.1	4.6	8.5	14.9
	2	0.2	0.6	1.4	4.9	18.5	33.7	58.3
	10	1.0	3.9	8.5	25.2	74.2	173.4	222.3
Distance, thousand km	20	0.0	0.0	0.0	0.0	0.1	0.2	0.3
	100	0.0	0.0	0.0	0.2	0.6	0.8	1.6
	500	0.0	0.1	0.2	0.8	3.2	3.9	7.9
	2.500	0.1	0.4	0.9	4.1	16.4	19.8	40.6
Excluding ${\overrightarrow{\ddot{r}}}_c$		2.0	7.9	17.6	78.9	314.4	379.3	778.3
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
Error parameter		Forcast interval, h
		0.5	1.0	1.5	3.0	6.0	12.0	24.0
Angle, degrees	0,1	0.021	0.084	0.2	0.7	2.2	9.9	43.3
	0,5	0.114	0.455	1.0	4.0	13.4	52.4	210.3
	2	0.401	1.578	3.5	13.1	39.9	189.2	755.1
	10	1.756	6.871	15.4	62.6	216.0	940.6	3986.0
Distance, thousand km	20	0.0	0.0	0.0	0.1	0.2	0.9	0.976
	100	0.0	0.0	0.1	0.3	0.8	4.5	4.9
	500	0.0	0.2	0.4	1.3	4.5	22.8	24.6
	2.500	0.2	0.8	1.8	6.5	21.6	109.5	118.1
Excluding ${\overrightarrow{\ddot{r}}}_c$		4.5	17.4	38.2	135.2	451.0	2284.5	2464.0

 

3.5. Forecast errors due to the Moon’s gravity

In assessing the influence of the Moon’s gravity on the forecast of the spacecraft’s position, two parameters are considered (Table 7):

• the error in determining the angular position of the Moon, similar to the Sun;
• error in determining the distance to the Moon due to periodic changes in the radius of the lunar orbit from 356.410 to 740 km.

 

Table 7

Forecast errors due to the Moon’s attraction

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Error parameter		Forcast interval, h
		0.5	1.0	1.5	3.0	6.0
Angle, degrees	0.1	0.0	0.0	0.0	0.0	0.0
	0.5	0.1	0.1	0.1	0.1	0.2
	2	0.2	0.3	0.3	0.4	0.6
	10	1.2	1.6	1.9	3.6	6.5
Distance, thousand km	200	0.0	0.0	0.0	0.0	0.0
	1.000	0.0	0.1	0.1	0.1	0.2
	5.000	0.1	0.3	0.3	0.5	0.9
	25.000	0.3	1.4	1.5	2.6	4.9
Excluding ${\overrightarrow{\ddot{r}}}_m$		1.3	5.9	6.1	10.8	20.4

 

End of Table 7

MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Error parameter			Forcast interval, h
			0.5	1.0	1.5		3.0		6.0		12.0		24.0
Angle, degrees	0.1		0.0	0.1	0.2		0.7		3.0		5.4		6.0
	0.5		0.2	0.6	1.1		3.0		11.9		21.9		31.8
	2		0.6	2.1	4.4		14.0		69.4		121.1		133.1
	10		2.6	9.7	20.3		73.0		255.4		570.3		570.3
Distance, thousand km	200		0.0	0.0	0.1		0.4		1.2		1.4		3.0
	1.000		0.0	0.2	0.4		1.8		6.2		6.9		15.2
	5.000		0.2	0.8	1.9		8.6		30.2		34.0		74.5
	25.000		1.2	5.0	11.4		51.3		178.9		200.1		442.6
Excluding ${\overrightarrow{\ddot{r}}}_m$			5.1	20.4	46.6		210.8		758.7		885.9		1855.5
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
Error parameter			Forcast interval, h
			0.5	1.0	1.5	3.0		6.0		12.0		24.0
Angle, degrees	0.1		0.0	0.1	0.2	0.7		3.0		18.6		24.4
	0.5		0.1	0.3	0.8	2.7		12.6		92.4		140.7
	2		0.5	2.0	4.3	15.8		71.1		452.0		590.0
	10		2.8	11.1	24.5	89.3		413.5		2373.7		4052.2
Distance, thousand km	200		0.0	0.0	0.1	0.3		1.2		5.0		15.7
	1.000		0.0	0.2	0.4	1.5		6.1		25.1		78.0
	5.000		0.2	0.9	1.9	7.7		30.7		127.2		394.8
	25.000		0.3	1.2	2.6	10.3		41.2		170.8		530.3
Excluding ${\overrightarrow{\ddot{r}}}_m$			5.9	22.9	50.6	204.8		831.2		3355.3		>1.0000

 

The forecast error is also given in the case where the influence of the Moon’s gravity is not taken into account.

The influence of the Sun and Moon on the spacecraft depends on the type of orbit – low or high. For NKO, the accuracy of calculating the position of the Moon and the Sun may not be critical, minimal accuracy is required to achieve optimal results. For high orbits, the influence of the secelestial bodies becomes much more significant, requiring more accurate models and calculations. When moving in orbits that have the properties of CEO and NKO, for example, high elliptical, it is necessary to take into account the features of each orbit depending on the phase of the object’s motion.

3.6. Effect of solar radiation pressure

Varying the effective reflection coefficient Crefl relative to the zero value gives the following results (Table 8).

 

Table 8

Crefl forecast errors

NKO. Orbital period: Т = 1.5 h. Altitude 500 km. Forecast error, m:
Crefl	Forcast interval, h
	0.5	1.0		1.5		3.0		6.0
10–8	0.0	0.1		0.3		0.6		0.9
10–7	0.2	0.9		3.0		5.7		9.3
10–6	1.6	9.4		30.4		57.0		93.1
10–5	16.0	94.3		303.8		570.3		931.0
MEO. Orbital period: Т = 12 h. Altitude 20,000 km. Forecast error, m:
Crefl	Forcast interval, h
	0.5	1.0	1.5		3.0		6.0	12.0	24.0
10–8	0.0	0.1	0.1		0.55		2.6	8.4	16.9
CEO. Orbital period: Т = 24 h. Altitude 35,777 km. Forecast error, m:
10–8	0.0	0.1	0.1		0.6		2.0	13.2	61.2

 

When forecasting motion along NKO, it should not be forgetten to exclude accelerations due to light pressure in areas where the Sun is shadowed by the Earth.

As in the case of the ballistic coefficient, an order of magnitude increase in the effective reflection coefficient error increases the forecast error by an order of magnitude, which is also confirmed by the linear dependence in (9) and does not depend on the orbit nonlinearity.

4. Using the forecast model when processing real data

During the flight tests of the object on the NKO, a set of instantaneous KS-solutions was obtained using satellite navigation. Due primarily to the errors in the ephemeris and time support of GNSS, the sequence of instantaneous KS solutions should be a sawtooth broken line around the real trajectory. Therefore, the standard was the SRO calculated from instantaneous KS-solutions by the least squares method using the forecast (to reduce individual KS solutions to a common moment) and subsequent “reproduction” of the calculated SRO by all the original instantaneous KS-solutions.

The forecast parameters were selected based on the data given above in Section 3 for a forecast interval of 4.5 hours, an expected accuracy of the instantaneous position of 15 m (3σ) and an acceptable forecast error of no more than 2 m:

• harmonics of the Earth’s potential up to the 32^nd order;
• taking into account accelerations from the gravity of the Sun and the Moon without errors;
• ballistic coefficient 0.000625;
• effective reflection coefficient 10^-7.

The graph of deviations of instantaneous KS-solutions from the calculated SRO is shown in the figure.

 

График отклонений реальных мгновенных КС-решений объекта на НКО от рассчитанной СрО (3 витка – 4,5 ч)

Graph of deviation of change of instantaneous KS-solutions of the object on NKO from the calculated SRO (3 seconds – 4.5 hours)

 

The figure shows that the nature of instantaneous KS-solutions corresponds to the expected; there is no obvious deterioration inaccuracy overtime, which indicates the correct choice of motion forecast parameters. The error of real KS-solutions relative to the calculated SRO is 15 m and 4 cm/s (3 σ).

Conclusion

The results presented in the work demonstrate that the task of developing a spacecraft motion model for navigation modules operating in real time has been successfully solved. The proposed model provides an effective compromise between the accuracy of trajectory prediction and computational complexity, which is especially important for devices with limited resources.

The main results of the work are the following:

• development of a model based on kinematic equations with adjustments that take into account the main dynamic effects (gravitational anomalies, atmospheric drag, the influence of the Sun and Moon gravity, sunlight pressure);
• successful testing of the model in real conditions, confirming its applicability for navigation tasks, especially for small spacecrafts;
• results obtained by varying the model parameters, which demonstrate the effect of incomplete data on the parameters and motion conditions of the spacecraft, the navigation of which is carried out using satellite navigation equipment.

This article serves as a basis for further study of methods and algorithms aimed at determining the location based on current navigation data. It should be noted that the motion model is the basis on which the Kalman filter can be effectively used as a navigation algorithm [18], allowing to improve the accuracy and reliability of navigation systems.

In the future, it is possible to further improve the model and optimize the computational processes, for example, through the study and transition to regular quaternion equations in other variables [19].

About the authors

Ivan A. Sokolov

JSC “Design Bureau of Navigation Systems”; St. Petersburg University

Author for correspondence.
Email: ivansokolof1997@mail.ru
ORCID iD: 0009-0008-2709-7080

postgraduate student, St. Petersburg State University; leading engineer of the laboratory of basic software for navigation equipment for space users, JSС “Design Bureau of Navigation Systems”

Russian Federation, 3, Kulneva St., Moscow, 121170; 7-9, Universitetskaya Embankment, St. Petersburg, 199034

Gennady S. Tsekhanovich

JSC “Design Bureau of Navigation Systems”

Email: ggsstt49@mail.ru
ORCID iD: 0009-0007-3748-7704

Cand. Sc., head of the laboratory of basic software for navigation equipment for space users

Russian Federation, 3, Kulneva St., Moscow, 121170

References

Supplementary files