ON THE INFLUENCE OF GEOPHYSICAL FACTORS ON SATELLITE MOTION

I. M. Yaglomsky

The first artificial Earth satellites will move in comparatively low orbits, since reaching great altitudes is associated with considerable expenditure of energy. These satellites must for the most part be in a region where the effect of air resistance is still noticeable; as a result of it they will gradually descend and ultimately fall to Earth or burn up in the atmosphere. In addition to the force of air resistance, the satellite will be acted upon by a disturbing force caused by the difference between the gravitational field of the spheroidal Earth and a central field. Finally, anomalies of gravity, due to the difference between the true gravitational field and the field of the terrestrial spheroid, also have a certain significance.

The influence of all the forces listed above, which we shall call geophysical factors, must be taken into account in calculating the orbit for the purpose of predicting the satellite’s motion. Clearly, the accuracy with which these forces are known will determine the accuracy of prediction of the orbit for given initial conditions of motion.

On the other hand, one may pose the problem of refining the coefficients entering into the mathematical expressions for the forces of interest to us on the basis of measurements of the satellite’s coordinates. Such a determination of the field of acting forces from satellite observations, if it proves possible, will subsequently permit more accurate prediction of satellite orbits.

Below we give some estimates of the accuracy currently possible in taking into account the influence of geophysical factors on the motion of a satellite, and also present preliminary considerations on refining the field of the indicated forces from the results of measurements of its orbit.

In calculating the motion of a satellite we shall use the system of differential equations, well known in astronomy, for the elliptic osculating elements[^6]. These equations are also given in papers[^8][^9] with certain changes which are not essential for the problems considered here.

The elliptic osculating elements $\Omega, i, \omega, p, e, \tau$ completely determine the position of the satellite in space and, for each instant of time, for known perturbations $S, T, W$ can be computed by direct integration of the indicated system. In astronomy, however, the method of successive approximations is usually used. This method gives good results for comparatively small perturbations, precisely such as are the perturbations caused by the action of geophysical factors.

I. M. YATSUNSKII

§ 1. NONSPHERICITY OF THE EARTH

As is known, the Earth in its shape is close to a spheroid (or to an ellipsoid of rotation). The gravitational potential for the spheroidal Earth can, with accuracy sufficient for our purposes, be expressed as follows:

\[ V_{\mathrm e}=\frac{\mu}{r}-\varepsilon\,\frac{\sin^2\psi-\frac{1}{3}}{r^3}. \tag{1} \]

The constant \(\varepsilon\), which depends on the parameters of the general terrestrial ellipsoid, has the form

\[ \varepsilon=\mu a^2\left(\alpha-\frac{m}{2}\right). \tag{2} \]

Here \(a\) is the semimajor axis of the general terrestrial ellipsoid, \(\alpha=\dfrac{a-b}{a}\) is the flattening of the Earth; \(m=\dfrac{\tilde{\omega}^{\,2}a}{g_a}\), \(b\) is the semiminor axis of the terrestrial ellipsoid, \(\tilde{\omega}\) is the angular velocity of the Earth’s rotation about its axis, \(g_a\) is the acceleration of gravity at the equator of the general terrestrial ellipsoid, and \(\mu\) is the Earth’s gravitational constant.

Differentiating \(V_{\mathrm e}\) in the directions \(S, T, W\), one readily obtains the following expressions for the disturbing accelerations:

\[ S=\frac{3\varepsilon}{r^4}\left(\sin^2 i\sin^2 u-\frac{1}{3}\right),\quad T=-\frac{\varepsilon}{r^4}\sin^2 i\sin 2u,\quad W=-\frac{\varepsilon}{r^4}\sin 2i\sin u. \tag{3} \]

Substituting the values of \(S, T, W\) into the right-hand sides of the system of equations of the osculating elements and integrating each of the equations separately, with the quantities \(\Omega, i,\ldots\) in the right-hand sides taken as constant and equal to their initial values, we obtain, in the first approximation, the changes in the elements of the elliptical orbit over the time of integration.

Without dwelling on the derivation of these formulas, since similar derivations have been made earlier, we shall give the results of integration with respect to \(\vartheta\) from \(0\) to \(2\pi n\) (where \(n\) is the number of revolutions in the angle \(\vartheta\)), i.e., we single out the so-called secular terms:

\[ \Omega-\Omega_0=\Delta\Omega_{\mathrm v} =-\frac{\varepsilon\cos i_0}{\mu p_0^2}\,2\pi n, \tag{4} \]

\[ i-i_0=\Delta i_{\mathrm v}=0, \tag{5} \]

\[ p-p_0=\Delta p_{\mathrm v}=0, \tag{6} \]

\[ e-e_0=\Delta e_{\mathrm v}=0, \tag{7} \]

\[ \omega-\omega_0=\Delta\omega_{\mathrm v} =\frac{\varepsilon}{\mu p_0^2}\left(4-5\sin^2 i_0\right)\pi n, \tag{8} \]

\[ \tau-\tau_0=\Delta\tau_{\mathrm v} =-\frac{\varepsilon}{\mu\sqrt{p_0\mu}}\,(1+3e_0)\left(3\sin^2 i_0\sin^2\omega_0-1\right)2\pi n^*). \tag{9} \]

*) In deriving the formula for \(\tau\), an expansion in powers of \(e\) was used, retaining terms of order \(e^2\).

Of the six osculating elements, only three \((\Omega, \omega, \tau)\) undergo secular perturbations caused by the influence of the ellipticity of the Earth. The parameter \((p)\) and the eccentricity \((e)\) of the orbit, as well as the inclination \((i)\), have no secular perturbations, at least in the first approximation.

Owing to the action of the secular terms \(\Omega\) and \(\omega\), the orbit will, with time, significantly change its orientation in space. The maximum rate of rotation of the node \((\Omega)\) and of the perigee \((\omega)\) for low orbits is \(\Delta \Omega_{\mathrm{v}} = 26'\), \(\Delta \omega_{\mathrm{v}} = 52'\) per revolution. In one day (14–16 revolutions) the displacement of the node reaches \(7^\circ\) and that of the perigee \(14^\circ\).

Let us determine with what accuracy it is possible to compute the secular deviations of the orbital parameters. As for the periodic terms, only the prediction within the limits of one revolution will depend on the possible accuracy with which they are taken into account. The error is not connected with the number of the revolution. However, in some cases, as will be seen from what follows, this error is not negligibly small.

Differentiating \(\Delta \Omega_{\mathrm{v}}\) with respect to \(a\), \(\alpha\), and \(m\), we shall have

\[ \delta(\Delta \Omega_{\mathrm{v}})_a = 2\Delta \Omega_{\mathrm{v}}\frac{\delta a}{a}, \tag{10} \]

\[ \delta(\Delta \Omega_{\mathrm{v}})_\alpha = -\frac{\Delta \Omega_{\mathrm{v}}}{a-\dfrac{m}{2}}\,\delta\alpha, \tag{11} \]

\[ \delta(\Delta \Omega_{\mathrm{v}})_m = -\frac{\Delta \Omega_{\mathrm{v}}}{2\left(a-\dfrac{m}{2}\right)}\,\delta m. \tag{12} \]

Similarly, for \(\Delta \omega_{\mathrm{v}}\) we obtain

\[ \delta(\Delta \omega_{\mathrm{v}})_a = 2\Delta \omega_{\mathrm{v}}\frac{\delta a}{a}, \tag{13} \]

\[ \delta(\Delta \omega_{\mathrm{v}})_\alpha = -\frac{\Delta \omega_{\mathrm{v}}}{a-\dfrac{m}{2}}\,\delta\alpha, \tag{14} \]

\[ \delta(\Delta \omega_{\mathrm{v}})_m = -\frac{\Delta \omega_{\mathrm{v}}}{2\left(a-\dfrac{m}{2}\right)}\,\delta m. \tag{15} \]

According to present-day data¹ the semi-major axis \(a\) may be regarded as known with an error of the order of \(100\) m, the flattening \(\alpha\) with an error of no more than two units in the denominator, and the quantity \(|\delta m| < 10^{-7}\). Substituting the values of the indicated errors, we shall have
\(|\delta \Omega_a| < 1''\), \(|\delta \Omega_\alpha| < 6'\), \(|\delta \Omega_m| < 1''\),
\(|\delta \omega_a| < 1''\), \(|\delta \omega_\alpha| < 12'\), \(|\delta \omega_m| < 2''\).
Thus the largest value is assumed by the possible error in the value of the Earth’s flattening \(\alpha\). For the orbits considered (with \(i\) close to zero), in one day an error of \(6'\) may accumulate in the position of the node and about \(12'\) in the position of the perigee.

The value of the gravitational constant \(\mu\) is not reflected directly in the motion of the node and the perigee, but the time \(\tau\) of the satellite’s passage through perigee depends on \(\mu\).

Varying formula (9) with respect to \(a\), \(\alpha\), \(m\), and \(\mu\), we obtain*)

\[ \delta(\Delta \tau_{\mathrm{v}})_a = 2\Delta \tau_{\mathrm{v}}\frac{\delta a}{a}, \tag{16} \]

\[ \delta(\Delta \tau_{\mathrm{v}})_\alpha = -\frac{\Delta \tau_{\mathrm{v}}}{a-\dfrac{m}{2}}\,\delta\alpha, \tag{17} \]

\[ \delta(\Delta \tau_{\mathrm{v}})_m = -\frac{\Delta \tau_{\mathrm{v}}}{2\left(a-\dfrac{m}{2}\right)}\,\delta m, \tag{18} \]

\[ \delta(\Delta \tau_{\mathrm{v}})_\mu = -\frac{1}{2}\Delta \tau_{\mathrm{v}}\frac{\delta\mu}{\mu}. \tag{19} \]

*) Here the relation between the quantities \(a\), \(\alpha\), \(m\), and \(\mu\) is not taken into account.

Taking, as before, \(\delta a=100\) m, \(\delta a=\dfrac{a}{450}\), \(\delta m<10^{-7}\), and for \(\delta\mu\) a value equal to two units in the fifth decimal place,\(^1\) we obtain the following errors (for \(n=15\)):

\[ \begin{aligned} \left|\delta(\Delta\tau_{\mathrm{B}})_a\right|&<0.004\ \text{sec.},& \left|\delta(\Delta\tau_{\mathrm{B}})_\alpha\right|&<1.7\ \text{sec.},\\ \left|\delta(\Delta\tau_{\mathrm{B}})_m\right|&<0.004\ \text{sec.},& \left|\delta(\Delta\tau_{\mathrm{B}})_\mu\right|&<0.002\ \text{sec.} \end{aligned} \]

In this case as well, the error in the flattening \(\delta\alpha\) has the largest value.

Over the course of a day the time of passage through perigee can acquire a secular change of up to \(1.7\) sec. as a result of the inaccurate knowledge, chiefly, of the quantity \(\alpha\). A very favorable factor from the point of view of the possibility of refining the constants associated with the Earth from observations of the satellite’s motion is the predominant influence of one parameter, namely the magnitude of the flattening \(\alpha\)*). The latter makes it possible to hope to refine, above all, this least accurately known parameter. However, such a problem proves not to be as simple as it may seem at first glance. The complexity of this problem will become clear after an analysis of the action of the anomalies of the force of gravity.

Let us estimate the magnitude of the error in the periodic terms that appear when the equations for the osculating elements are integrated from \(0\) to \(\vartheta\), if for \(S, T, W\) the expressions (3) are adopted. For these periodic terms one may write the estimates

\[ \left|\Delta\Omega_{\mathrm{per}}\right|<\frac{\varepsilon\cos i}{2\mu p^2}, \tag{20} \]

\[ \left|\Delta i_{\mathrm{per}}\right|<\frac{\varepsilon\sin 2i}{4\mu p^2}, \tag{21} \]

\[ \left|\Delta p_{\mathrm{per}}\right|<\frac{\varepsilon\sin^3 i}{\mu p}, \tag{22} \]

\[ \left|\Delta e_{\mathrm{per}}\right|<\frac{2\varepsilon}{\mu p^2}, \tag{23} \]

\[ \left|\Delta\omega_{\mathrm{per}}\right|<\frac{2\varepsilon}{\mu p^2 e}, \tag{24} \]

\[ \left|\Delta\tau_{\mathrm{per}}\right|<\frac{\varepsilon}{e\mu\sqrt{p\mu}}. \tag{25} \]

Inequalities (20)–(25) show that errors in the Earth’s constants will have the most substantial influence on the displacement of the perigee \((\omega)\) and the time of passage through perigee \((\tau)\) for orbits with small eccentricities.

Let us calculate the maximum values of the errors \(\delta(\Delta\Omega_{\mathrm{per}})\), \(\delta(\Delta i_{\mathrm{per}}),\ldots\), taking into account only the influence of the possible error in \(\alpha\), since it is the largest. In doing so we shall take

\[ p=7000\ \text{km}, \qquad e=0.05. \]

We obtain

\[ \delta(\Delta\Omega_{\mathrm{per}})_\alpha < \frac{a^2\cos i}{2p^2} =2'' \quad \text{for } i=0, \]

\[ \delta(\Delta i_{\mathrm{per}})_\alpha < \frac{a^2\sin 2i}{4p^2}\,\delta\alpha =1'' \quad \text{for } i=45^\circ, \]

\[ \delta(\Delta p_{\mathrm{per}})_\alpha < \frac{a^2\sin^2 i}{p}\,\delta\alpha =120\ \text{m} \quad \text{for } i=90^\circ, \]

\(^1\) Note that the error in \(\mu\) is reflected in the first term of expression (1) and can be revealed from observations of the satellite’s period of revolution.

\[ \delta(\Delta e_{\mathrm{per}})_a < \frac{2a^2\sin^2 i}{p^2}\,\delta a = 0.3\cdot 10^{-4} \quad \text{for } i=90^\circ, \]

\[ \delta(\Delta\omega_{\mathrm{per}})_a < \frac{2a^2}{p^2 e}\,\delta a = 2', \]

\[ \delta(\Delta\tau_{\mathrm{per}})_a < \frac{a^2}{e\sqrt{p\mu}}\,\delta a = 0.3 \text{ sec.} \]

Obviously, the errors in the position of the node and in the value of the inclination, as well as in \(p\) and \(e\), are very insignificant. The remaining two errors are more noticeable for the orbits under consideration.

§ 2. GRAVITY ANOMALIES

The field of gravity anomalies for the entire Earth can be represented in mathematical form as an expansion in spherical functions.

Let us write the expression for the potential of gravity anomalies in the form

\[ V_a=\gamma_m \sum_{n=2}^{\infty}\sum_{m=0}^{n} \left(\frac{R}{r}\right)^{n+1} (\alpha_{nm}\cos m\lambda+\beta_{nm}\sin m\lambda)P_{nm}(\sin\varphi). \tag{26} \]

Here \(\gamma_m\) is the mean value of gravity for the geoid, \(R\) is the mean radius of the Earth, \(r\) is the radius vector of the satellite, \(\varphi,\lambda\) are the geographic latitude and longitude of the satellite, and \(\alpha_{nm}, \beta_{nm}\) are the coefficients of the expansion, whose numerical values are determined by appropriate processing of gravimetric data.

Formula (26) is easily obtained from the general expression for the disturbing potential \(V_a\) in the entire external field\({}^{2}\), if instead of \(\Delta g\) (anomalies) one substitutes into it the corresponding expansion in spherical functions\({}^{1}\).

Differentiating (26) in the directions \(S,T,W\), we obtain

\[ S=-\frac{\gamma_m}{R}\sum_{n=2}^{\infty}(n+1) \left(\frac{R}{r}\right)^{n+2} \sum_{m=0}^{n} (\alpha_{nm}\cos m\lambda+\beta_{nm}\sin m\lambda)P_{nm}(\sin\varphi), \tag{27} \]

\[ T=\frac{\gamma_m}{R}\sum_{n=2}^{\infty}\left(\frac{R}{r}\right)^{n+2} \sum_{m=0}^{n} (\alpha_{nm}\cos m\lambda+\beta_{nm}\sin m\lambda) P'_{nm}(\sin\varphi)\cos\delta+ \]

\[ +\frac{\gamma_m}{R\cos\varphi}\sum_{n=2}^{\infty}\left(\frac{R}{r}\right)^{n+2} \sum_{m=0}^{n} (-\alpha_{nm}\sin m\lambda+\beta_{nm}\cos m\lambda) mP_{nm}(\sin\varphi)\sin\delta, \tag{28} \]

\[ W=\frac{\gamma_m}{R}\sum_{n=2}^{\infty}\left(\frac{R}{r}\right)^{n+2} \sum_{m=0}^{n} (\alpha_{nm}\cos m\lambda+\beta_{nm}\sin m\lambda) P'_{nm}(\sin\varphi)\sin\delta- \]

\[ -\frac{\gamma_m}{R\cos\varphi}\sum_{n=2}^{\infty}\left(\frac{R}{r}\right)^{n+2} \sum_{m=0}^{n} (-\alpha_{nm}\sin m\lambda+\beta_{nm}\cos m\lambda) mP_{nm}(\sin\varphi)\cos\delta. \tag{29} \]

Here \(\delta\) is the current azimuth; a prime denotes differentiation with respect to \(\varphi\).

The number of terms of the expansion (harmonics) \(n\) must be sufficiently large if we wish to have any satisfactory accuracy. At present, however, the absence of detailed gravimetric material for the entire Earth’s surface\({}^{1}\) makes the use of higher harmonics meaningless. Hence the accuracy of such expansions

very small. In estimating the influence of anomalies of the force of gravity on the motion of the satellite, we use the materials¹, bearing in mind, however, that the error in allowing for the action of the anomalies may exceed 100%. Thus, at present it is possible only approximately to establish the qualitative character and the order of magnitude of the deviations in the satellite’s motion caused by the action of anomalies of the force of gravity.

Changes in the osculating elements were calculated for one of the possible orbits (mean altitude 500 km) over the time of five revolutions. In view of the great computational difficulties, the calculations were carried out only with allowance for the second harmonic, i.e., with terms constituting the greater part of the perturbing accelerations. Analysis of the results of these calculations makes it possible to draw the following conclusions. The changes in the osculating elements are, for the most part, oscillatory in character; moreover, with the passage of time (within five revolutions) the amplitude of the oscillations increases. In addition, for some elements a gradually accumulating component is observed, i.e., with the passage of time the level about which the oscillations take place is displaced. This phenomenon can be explained by the presence of long-period oscillations connected with the period of rotation of the Earth about its axis.

In absolute value, some of the deviations reach appreciable magnitudes. Thus, for example, the changes in the orbital parameter at the end of the period considered become equal to 1–1.2 km. The changes in eccentricity reach \(2 \cdot 10^{-4}\), the displacement of perigee—to several minutes of arc, \(\tau\)—to one tenth of a second. The changes in \(\Omega\) and \(i\) remain within 0.5.

These estimates are only qualitative in character, and the true field of anomalies may cause somewhat different deviations. Nevertheless, it appears possible to judge the order of magnitude of the deviations.

§ 3. AIR RESISTANCE

The force of air resistance, unlike the forces of gravitation, is not conservative. The work performed by this force reduces the total energy of the satellite, if it is regarded as a moving material point. As a result, the dimensions of the orbit gradually decrease; the satellite enters the dense layers of the atmosphere and, depending on the heat resistance of its shell, may burn up or fall to Earth.

Let us estimate the possible accuracy of allowing for the forces of air resistance, depending on the accuracy of our knowledge of the satellite’s drag coefficient and of the distribution of atmospheric density with altitude.

First of all, we have quite insufficient data on the determination of density as a function of altitude for altitudes above 150 km. It is true that in recent years some fragmentary data have appeared on the density up to 200 km and higher, obtained on the basis of rocket investigations. However, these results are isolated and cannot be extended to all regions where artificial satellites will fly. Moreover, these data differ noticeably from the data of density determinations by other methods. As will be seen below, for satellites close to the Earth a special role is played by the characteristics of atmospheric resistance at altitudes approximately from 200 to 500 km. Concerning possible values of the density at these altitudes, let us cite one statement by L. Spitzer³: “... neither the study of the ionosphere nor spectroscopic observations permit an unambiguous estimate of the temperature and density of particles at an altitude of 300 km. The data given above indicate only that the particle density must be, at least, so

same as the electron density, equal to \(10^6\ \mathrm{cm}^{-3}\), and cannot greatly exceed \(10^{10}\ \mathrm{cm}^{-3}\)... It seems probable that the particle density is considerably greater than \(10^8\ \mathrm{cm}^{-3}\),... we shall assume that the true value of the density differs little from \(10^{10}\ \mathrm{cm}^{-3}\). Although at present this value seems the most plausible, we cannot regard as completely excluded the possibility that the density at an altitude of 300 km may for a large part of the time be equal to \(10^6\ \mathrm{cm}^{-3}\)." According to other data\(^7\), obtained from an analysis of the peculiarities of the behavior of the \(F_2\) layer, the concentration of particles in the \(F_2\) layer is of the order of \(10^8\ \mathrm{cm}^{-3}\). Approximately the same value of the particle concentration is adopted in \(^4\) for an altitude of 300 km \((\sim 7 \cdot 10^8\ \mathrm{cm}^{-3})\).

We have deliberately cited data on the density at an altitude of 300 km from various sources in order to show that, in choosing any of the given density values, one can easily be in error by more than an order of magnitude.

The density of the atmosphere is known somewhat more accurately at altitudes below 150 km, where it has repeatedly been determined from meteor observations. Spitzer adopts for an altitude of 150 km a particle concentration of \(10^{12}\ \mathrm{cm}^{-3}\), Mitra—\(2.5 \cdot 10^{-11}\ \mathrm{cm}^{-3}\). Below, in analyzing the probable accuracy of taking air resistance into account, we shall use two atmospheres—the model given by Mitra\(^4\), and a model specially constructed from Spitzer’s data, taking in the latter case, at an altitude of 300 km, the particle concentration to be \(10^{10}\ \mathrm{cm}^{-3}\). The latter atmosphere is regarded by us as the densest, i.e. as having the strongest effect on the satellite’s orbit.

The second factor on which the trajectory of a satellite in the atmosphere depends is the drag coefficient \(c_x\). In the expression for the drag force, \(c_x\) enters as a multiplier, just as does the density. Therefore the relative error in \(c_x\) will affect the motion of the satellite to the same extent as the relative error in the density. At present the drag coefficient \(c_x\) for simple bodies (sphere, cylinder) is calculated theoretically with comparative accuracy\(^5\). Despite the fact that the error of such a calculation may be large in absolute value, it is in any case many times smaller than the uncertainty in determining the air density.

Before proceeding to the derivation of formulas for calculating the orbit with allowance for drag, we shall make two remarks. First, calculations show that oscillations in the magnitude of the radius vector of the orbit due to the noncentrality of the gravitational force field may reach 10–15 km. At the same time, oscillations of the height above the level of the terrestrial spheroid are about half as large, i.e. only 5–8 km. Thus the flight altitude and, consequently, the air density can be computed with sufficient accuracy if the Earth is taken to be a sphere and the orbit to be an exact ellipse. Therefore, when calculating the orbit with allowance for air resistance, we shall neglect the noncentrality of the gravitational force field and the ellipticity of the Earth, assuming that the changes in the orbital parameters in this case will be close to the true ones.

The second remark is that, in deriving the computational formulas, we shall regard the eccentricity of the orbit as a small quantity \((e<0.1)\). Indeed, for \(e=0.1\) and a perigee altitude \(h_\pi = 300\) km, we have an apogee altitude \(h_a = 1770\) km. Thus such a restriction on eccentricities is quite permissible. It is easy to see that, as a result of the assumptions made, the case of a plane orbit is considered. Consequently, perturbations of the node and inclination will not occur. Secular perturbations of the perigee \(\omega\) will likewise be absent. The matter reduces only to the calculation of the secular perturbations of the parameter \(p\) and the eccentricity \(e\) of the orbit, i.e. to calculating the change of these quantities over an integer number of revolutions.

Indeed, for the acceleration due to the force of air resistance one may write

\[ R_x=\frac{C_x S_m}{2m}\rho v^2=\frac{1}{2}b\rho V^2, \]

where \(\rho\) is the air density, \(v\) is the velocity of motion along the orbit, \(S_m\) is the area of the midsection, \(m\) is the mass of the satellite,

\[ b=\frac{C_x S_m}{m}. \]

For the perturbing accelerations we shall have\(^9\)

\[ S=-\frac{R_x}{v}\sqrt{\frac{\mu}{p}}\,e\sin\vartheta =-\frac{1}{2}b\rho\frac{\mu}{p}e\sin\vartheta\left(1+2e\cos\vartheta+e^2\right)^{\frac12}, \]

\[ T=-\frac{R_x}{v}\sqrt{\frac{\mu}{p}}(1+e\cos\vartheta) =-\frac{1}{2}b\rho\frac{\mu}{p}(1+e\cos\vartheta)\left(1+2e\cos\vartheta+e^2\right)^{\frac12}, \]

\[ W=0. \]

After substituting the expressions \(S, T\) into the equations for the osculating elements, for \(p\) and \(e\) we obtain

\[ \left. \begin{aligned} \frac{dp}{d\vartheta} &=-b\rho p^2\, \frac{\sqrt{1+2e\cos\vartheta+e^2}}{(1+e\cos\vartheta)^2},\\ \frac{de}{d\vartheta} &=-b\rho p\, \frac{(e+\cos\vartheta)\sqrt{1+2e\cos\vartheta+e^2}}{(1+e\cos\vartheta)^2}. \end{aligned} \right\} \tag{30} \]

The time of passage through perigee \(\tau\) will not interest us, since we confine ourselves to considering complete revolutions. Before proceeding to the integration of the expressions obtained, it is necessary to choose an analytic dependence for the air density. One of the simplest and most commonly used expressions for \(\rho\) is an exponential function of the form

\[ \rho=Ae^{-kh}, \tag{31} \]

where \(h\) is the current altitude, and \(A\) and \(k\) are constants.

A comparison of the values of \(\rho\) according to (31) with the available data shows that the discrepancies may reach large values if one attempts, with the same constants \(A\) and \(k\), to describe the density over the whole range of working altitudes (from 100 to 500 km). Satisfactory results are obtained if formula (31) is applied over intervals of length 90–100 km.

For the Mitra atmosphere, such an approximation gives a maximum error of up to 20% in the middle of some intervals, which, given the existing data on density, is a small error. Replacing the altitude \(h\) by the formula

\[ h=r-R=\frac{p}{1+e\cos\vartheta}-R, \]

for the density we obtain

\[ \rho=Ae^{-k\left(\frac{p}{1+e\cos\vartheta}-R\right)}. \tag{32} \]

Replacing further \(\rho\) and \(r\) by their expressions in terms of \(p, e, \vartheta\), for the derivatives \(\dfrac{dp}{d\vartheta}\), \(\dfrac{de}{d\vartheta}\) we finally obtain

\[ \frac{dp}{dv} = -Abp^2 e^{-k\left(\frac{p}{1+e\cos\vartheta}-R\right)} \cdot \frac{\sqrt{1+2e\cos\vartheta+e^2}}{(1+e\cos\vartheta)^2}, \tag{33} \]

\[ \frac{de}{dv} = -Abpe^{-k\left(\frac{p}{1+e\cos\vartheta}-R\right)} \cdot \frac{(e+\cos\vartheta)\sqrt{1+2e\cos\vartheta+e^2}}{(1+e\cos\vartheta)^2}. \tag{34} \]

Now let us make use of the smallness of the eccentricity \(e\). Applying expansion in a power series, instead of (33) and (34) we obtain

\[ \frac{dp}{d\vartheta} = -Abp^2 e^{-k(p-R)} e^{\xi\cos\vartheta} \left( 1-\xi e\cos^2\vartheta+\frac{\xi^2}{2}e^2\cos^4\vartheta \right) (1-e\cos\vartheta), \tag{35} \]

\[ \frac{de}{d\vartheta} = -Abpe^{-k(p-R)} e^{\xi\cos\vartheta} \left( 1-\xi e\cos^2\vartheta+\frac{\xi^2}{2}e^2\cos^4\vartheta \right) \times \]

\[ \times(\cos\vartheta+e\sin^2\vartheta), \tag{36} \]

where \(\xi=kpe\).

In expanding the exponential function into a series, terms of order \(e^2\) were retained; in other cases, terms of order \(e\). Integrating (35) and (36) with respect to \(\vartheta\) over the limits from \(0\) to \(2\pi n\), we shall have, in the first approximation,

\[ p-p_0 = -2\pi nAbp_0^2 e^{-k(p_0-R)} \left[ \left(1-\xi e_0+\frac{\xi^2}{2}e_0^2\right)I_0(\xi) + \frac{e_0^2}{2}I_2(\xi) \right], \tag{37} \]

\[ e-e_0 = -2\pi nAbp_0 e^{-k(p_0-R)} \left[ \left( 1+\frac{e_0}{\xi}-\xi e_0-e_0^2+\frac{\xi^2}{2}e_0^2 \right)I_1(\xi) + \right. \]

\[ \left. + \left( e_0+3\frac{e_0^2}{\xi}-\xi e_0^2 \right)I_2(\xi) + \frac{3}{2}e_0^2 I_3(\xi) \right]. \tag{38} \]

Here \(I_0(\xi),\ I_1(\xi),\ I_2(\xi),\ I_3(\xi)\) are Bessel functions of imaginary argument.

Fig. 1.

The calculation by formulas (37) and (38) is carried out as follows. Given some number of revolutions \(n\), from the initial data \(p_0\) and \(e_0\) we compute the corresponding changes \((p-p_0)\) and \((e-e_0)\) and the new values of the parameter \(p\) and the eccentricity \(e\). Taking these \(p\) and \(e\) as initial values, we again prescribe \(n\) and find the following values of \(p\) and \(e\), and so on.

The number \(n\) is chosen so that the decrease in the height of the apogee is \(30 \div 40\) km.

The error of the calculation can be no more than 20% in comparison with numerical integration in the case where the coefficients \(k\) and \(A\) are chosen constant over altitude intervals of up to 50 km, and the transition to new \(k\) and \(A\) is made according to the height of the perigee.

Fig. 2.

In order to obtain an estimate of the accuracy of accounting for air resistance, the changes in \(p\) and \(e\) for a certain orbit in the atmospheres of Mitra and Spitzer were calculated from the formulas obtained. The results are given in Figs. 1, 2. The changes in the heights of the apogee and perigee are shown in Figs. 3, 4.

Fig. 3.

Fig. 4.

From these graphs it is evident that, depending on which of the two atmospheres is chosen, the changes in the heights of the apogee and perigee may differ by a factor of \(5 \div 7\). Consequently, the lifetime of a satellite may vary by the same factor.

§ 4. POSSIBILITY OF DETERMINING AIR DENSITY,
THE CONSTANTS OF THE TERRESTRIAL ELLIPSOID, AND ANOMALIES
OF GRAVITY FROM SATELLITE OBSERVATIONS
===========================================

As the estimates given show, the influence of geophysical factors on the motion of a satellite can be taken into account with a certain error. This error is especially significant for low-flying satellites, where air resistance is noticeable, i.e., for orbits with a perigee height of \(150 \div 300\) km. Determining the elements of these orbits from the results of measurements over a long period of time, we shall find a systematic decrease in the parameter and eccentricity caused by air resistance. Then, using formulas (37) and (38), from the known number of revolutions and from the measured \(p\) and \(e\), one can find the constants \(A\) and \(k\). These constants will correspond to the air density in a certain range of heights near perigee.

Passing to the next number \(n\) and to the elements \(p, e\), we find \(A\) and \(k\) for the next height interval. Such a method of determining the density function cannot be very accurate, all the more so because the coefficient \(c_x\), even for the simplest satellite shapes, is known with a large error. But even these results will be very useful, since values of the air density in the height range under consideration are practically absent.

Let us consider the possibility of refining the Earth’s gravitational field from the results of orbit measurements. This problem is considerably more complicated than the determination of air density. The point is that the action of air resistance is easily separated from the action of other forces if the satellite’s flight altitude is sufficiently low.

As for the gravitational field, however, it seems almost impossible to separate the influence of the various constants of this field. It is relatively easy to exclude only the action of air resistance. For this it is necessary to raise the perigee of the orbit above \(300 \div 350\) km, so that the satellite motion (it is assumed that the mean density of the satellite is not too small) is almost entirely determined by the Earth’s gravitational field.

As was clarified above, among the Earth’s constants the magnitude of the flattening \(\alpha\) is of principal importance for the satellite. If anomalies of gravity could not cause secular displacements of the node and perigee, then from observations of these displacements the flattening could be obtained almost unambiguously. However, secular displacements of the node and perigee due to anomalies are entirely probable, and in order of magnitude they may be the same as those due to \(\delta\alpha\). To separate the influence of \(\delta\alpha\) from anomalies it is expedient to have several orbits with different inclinations. In this case the part of the secular displacement of the node and perigee that depends on \(\delta\alpha\) will vary in a quite definite way (see formulas (4) and (8)). The other part, connected with anomalies, must vary according to some other, unknown law. Thus one can attempt to isolate the influence of \(\delta\alpha\), if it is seen that the main part of the displacements is due precisely to \(\delta\alpha\), and not to anomalies, and thereby refine the value of \(\alpha\) for the terrestrial ellipsoid. In principle, it is also possible to pose the problem of determining the coefficients of a series of spherical functions for the anomaly field, i.e., \(a_{nm}, \beta_{nm}\) (at least for the first harmonics), from observations of satellite motion.

Indeed, integrating from \(0\) to \(t\) each of the equations for the osculating elements with constant \(\delta\), \(i\), \(\omega\), \(p\), \(e\), \(\tau\) in the right-hand sides, and \(S, T, W\) according to (27), (28), (29), and varying the resulting expressions, we obtain—

reduce to the following system of algebraic equations:

\[ l_i=a_{i1}a_{00}+a_{i2}a_{20}+a_{i3}a_{21}+a_{i4}\beta_{21}+\ldots+a_{ik-1}a_{nn}+a_{ik}\beta_{nn} \tag{39} \]

\[ (i=1,2,\ldots,6). \]

Here the unknowns are the expansion coefficients \(a_{00}\) (included among the unknowns in order to determine (see below) the constants \(\alpha\) and \(\gamma_e\)), \(a_{20}\), \(a_{21}\), \(\beta_{21},\ldots,\beta_{nn}\). In all there are \(k=n^2+2n-2\) unknowns, where \(n\) is the order of the highest harmonic1. The quantities \(l_1,l_2,\ldots,l_6\) are the differences between the computed values of the osculating elements (respectively \(\Omega, i, \omega, p, e, \tau\)) for the instant of time \(t\) and their values obtained as a result of measurements of the coordinates for the same instant of time. The coefficients \(a_{ij}\) are calculated from the data of the computed orbit. The simplest method of computing these coefficients is the method of quadratures. Thus, for example,

\[ a_{13}= \frac{\gamma_m}{R\sqrt{p^4}\sin i} \int_0^t r\sin u \left(\frac{R}{r}\right)^4 \left[ \cos\lambda\, P'_{21}(\sin\varphi)\sin\delta + \frac{\sin\lambda}{\cos\varphi}\, P_{21}(\sin\varphi)\cos\delta \right]dt. \]

To obtain the required number of equations (\(\geq n^2+2n-2\)), expressions (39) are compiled for a number of values of the instants of time \(t\). The redundant number of such equations (which in this case are regarded as conditional) makes it possible to apply the method of least squares and to determine the most probable values of the coefficients \(a_{nm}, \beta_{nm}\).

To equations of the form (39), with the corresponding weights, one may also add the following:

\[ \begin{aligned} a_1a_{00}+a_2a_{20}+a_3a_{21}+\ldots+a_k\beta_{nn}&=\Delta g_1,\\ b_1a_{00}+b_2a_{20}+b_3a_{21}+\ldots+b_k\beta_{nn}&=\Delta g_2,\\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \end{aligned} \tag{40} \]

On the right-hand sides of equalities (40) stand certain values, averaged over preselected zones, of the anomalies of the force of gravity on the surface of the geoid. The coefficients \(a_i, b_i,\ldots\) are certain spherical functions specified in the same zones1. Joint processing by the method of least squares of the results of measurement of the orbit (equations (39)) and of the available gravimetric material (equations (40)) can increase the reliability of the determination of the sought coefficients. To improve the accuracy it is also expedient to use another path1, namely, from the values \(a_{nm}, \beta_{nm}\) found from the solution of equations (39) and (40) (or, what is the same,

\[ A_{nm}=\frac{\gamma_m}{R}(n-1)a_{nm},\qquad B_{nm}=\frac{\gamma_m}{R}(n-1)\beta_{nm} \]

) to determine these coefficients a second time, using the formulas

\[ \begin{aligned} A_{n0}&=\frac{2n+1}{4\pi}\int \Delta g\,P_{n0}(\sin\varphi)\,d\sigma,\\[4pt] B_{n0}&=0,\\[4pt] A_{nm}&=\frac{2n+1}{2\pi}\frac{(n-m)!}{(n+m)!} \int \Delta g\,P_{nm}(\sin\varphi)\cos m\lambda\,d\sigma,\\[4pt] B_{nm}&=\frac{2n+1}{2\pi}\frac{(n-m)!}{(n+m)!} \int \Delta g\,P_{nm}(\sin\varphi)\sin m\lambda\,d\sigma. \end{aligned} \tag{41} \]

In expressions (41) the integration is carried out over the surface of the unit sphere \(\sigma\). For the quantities \(\Delta g\), the anomaly values are taken as those computed according to the previously obtained coefficients \(\alpha_{nm}\) and \(\beta_{nm}\).

Equating the values found, \(A_{00}=a_{00}\), \(A_{20}=a_{20}\), we can also determine the flattening \(\alpha\) and the acceleration of gravity at the equator \(\gamma_e\) for the new general terrestrial ellipsoid. Indeed, for the distribution of the accelerations of gravity on the ellipsoid we have the expression1

\[ \gamma = a_{00} + a_{20} P_{20} + a_{40} P_{40}, \]

where the expansion coefficients \(a_{00}\), \(a_{20}\), and \(a_{40}\) are known functions of \(\alpha\) and \(\gamma_e\), namely

\[ \alpha = \frac{\frac{5}{2}m-\beta}{1+\frac{17}{14}m}, \qquad \beta = \left(\frac{3}{2}a_{20}+\frac{5}{8}a_{40}\right)\frac{1}{\gamma_e}, \]

\[ \gamma_e = a_{00} - \frac{1}{2}a_{20} + \frac{3}{8}a_{40}\,{}^*). \]

Conclusion

The analysis carried out of the influence of geophysical factors on the motion of a satellite close to the Earth, and the estimate of the possible accuracy with which they can be taken into account, show that, at the present level of knowledge in this field, errors in orbit prediction may prove to be large. On the basis of prolonged measurements of the orbit, it is quite possible to refine substantially the data on atmospheric density at altitudes from 150 to 300 km. It also appears possible to refine the value of the flattening of the terrestrial spheroid; however, reliable data can apparently be obtained only if there are several satellites whose orbits have sharply different inclinations.

The most difficult problem is the determination of the field of anomalies of gravity, since here it is necessary to determine, on the basis of measurements, the periodic changes of the osculating orbital elements. Obviously, the larger the territory over which measurements of the orbit are made, and the longer the interval of time during which they continue, the greater the confidence in the results will be.

References Cited

1. I. Zhongolovich, Proceedings of the Institute of Theoretical Astronomy, issue 3, 1952.
2. N. I. Idelson, Potential Theory with Applications to the Theory of the Figure of the Earth and Geophysics, 2nd ed., 1936.
3. The Atmospheres of the Earth and the Planets. Collection of articles edited by G. P. Kuiper, IL, 1951.
4. S. K. Mitra, The Upper Atmosphere, IL, 1955.
5. Gas Dynamics. Collection of articles edited by S. G. Popov and S. V. Fal’kovich, IL, 1950.
6. G. N. Duboshin, Introduction to Celestial Mechanics, Moscow–Leningrad, ONTI, 1938.
7. Ya. L. Al’pert, ZhETF 18, 995 (1948).
8. D. E. Okhotsimsky, T. M. Eneev, G. P. Taratynova—in this issue of UFN.
9. G. P. Taratynova—in this issue of UFN.

1. The quantities \(a_{40}\) and \(m\) may be regarded as known exactly. ↩↩↩↩