Preliminary orbits with line-of-sight correction for LEO satellites observed with radar

TL;DR

This paper introduces a method for preliminary orbit determination of LEO satellites using radar data that accounts for Earth's oblateness, improving initial orbit estimates by incorporating the J2 perturbation into the dynamical model.

Contribution

It extends previous Keplerian models by including Earth's oblateness (J2 effect) for more accurate orbit determination from radar observations.

Findings

Enhanced initial orbit accuracy for low Earth orbit satellites.

Improved convergence rates in orbit determination algorithms.

Effective correction of line-of-sight measurements considering Earth's oblateness.

Abstract

We propose a method to account for the Earth oblateness effect in preliminary orbit determination of satellites in low orbits with radar observations. This method is an improvement of the one described in (Gronchi et al 2015), which uses a pure Keplerian dynamical model. Since the effect of the Earth oblateness is strong at low altitudes, its inclusion in the model can sensibly improve the initial orbit, giving a better starting guess for differential corrections and increasing the chances to obtain their convergence. The input set consists of two tracks of radar observations, each one composed of at least 4 observations taken during the same pass of the satellite. A single observation gives the topocentric position of the satellite, where the range is very accurate, while the line of sight direction is poorly determined. From these data we can compute by a polynomial fit the values of the range and range rate at the mean epochs of the two tracks. In order to obtain a preliminary orbit we wish to compute the angular velocities, that is the rate of change of the line of sight. In the same spirit of (Gronchi et al 2015), we also wish to correct the values of the angular measurements, so that they fit the selected dynamical model if the same holds for the radial distance and velocity. The selected model is a perturbed Keplerian dynamics, where the only perturbation included is the secular effect of the $J_2$ term of the geopotential. The proposed algorithm models this problem with 8 equations in 8 unknowns.

Equations

$$(t_j,{\rho }_j,{\alpha }_j,{\delta }_j),j=1,\dots ,4$$

$$(\bar{t},\bar{\alpha },\bar{\delta },\rho ,\dot{\rho },\ddot{\rho }),$$

$${\mathcal{E}}_{att}=(\alpha ,\delta ,\dot{\alpha },\dot{\delta },\rho ,\dot{\rho }).$$

$$\alpha =\bar{\alpha }+\Delta \alpha ,\delta =\bar{\delta }+\Delta \delta .$$

$$\xi =\rho \dot{\alpha }\cos \delta ,\zeta =\rho \dot{\delta },$$

$${\mathcal{E}}_{att}^{\ast }=(\Delta \alpha ,\Delta \delta ,\xi ,\zeta ,\rho ,\dot{\rho }),$$

$$({\bar{t}}_1,{\bar{\alpha }}_1,{\bar{\delta }}_1,{\rho }_1,{\dot{\rho }}_1,{\ddot{\rho }}_1),({\bar{t}}_2,{\bar{\alpha }}_2,{\bar{\delta }}_2,{\rho }_2,{\dot{\rho }}_2,{\ddot{\rho }}_2).$$

$$\mathbf{r}=\mathbf{q}+\rho {\mathbf{e}}^{\rho },$$

$${\mathbf{e}}^{\rho }=(\cos \delta \cos \alpha ,\cos \delta \sin \alpha ,\sin \delta ).$$

$${\mathbf{e}}^{\alpha }=\left(-\cos \delta \sin \alpha ,\cos \delta \cos \alpha ,0\right),{\mathbf{e}}^{\delta }=\left(-\sin \delta \cos \alpha ,-\sin \delta \sin \alpha ,\cos \delta \right).$$

$$\dot{\mathbf{r}}=\xi {\mathbf{e}}^{\alpha }+\zeta {\mathbf{e}}^{\delta }+(\dot{\rho }{\mathbf{e}}^{\rho }+\dot{\mathbf{q}}).$$

$${\mathcal{E}}_{kep}$$

$${\mathcal{E}}_{car}$$

$${\mathcal{E}}_{att}$$

$${\mathcal{E}}_{att}^{\ast }$$

$${\mathcal{E}}_{kep}\stackrel{{\varphi }_1}{\to }{\mathcal{E}}_{car},{\mathcal{E}}_{car}\stackrel{{\varphi }_2}{\to }{\mathcal{E}}_{att},{\mathcal{E}}_{att}\stackrel{{\varphi }_3}{\to }{\mathcal{E}}_{att}^{\ast },$$

$$\Phi ={\varphi }_3\circ {\varphi }_2\circ {\varphi }_1$$

$$\ddot{\mathbf{r}}=\nabla U(\mathbf{r}),$$

$$U({\bf r})=\frac{\mu}{r}\Bigl{[}1-J_{2}\Bigl{(}\frac{R_{\oplus}}{r}\Bigr{)}^{2}P_{2}(\sin\delta)\Bigr{]}.$$

$$P_2(\sin \delta )=\frac{3}{2}{\sin }^2\delta -\frac{1}{2}=\frac{3}{2}\frac{z^2}{r^2}-\frac{1}{2}.$$

$$\begin{cases}\quad\dot{a}=0\\ \quad\dot{e}=0\\ \quad\dot{I}=0\\ \quad\dot{\Omega}=\displaystyle-\frac{3}{2}J_{2}\frac{R_{\oplus}^{2}}{p^{2}}n_{0}\cos{I}\\ \quad\stackrel{{\scriptstyle}}{{\dot{\omega}}}=\displaystyle\frac{3}{4}J_{2}\frac{R_{\oplus}^{2}}{p^{2}}n_{0}(4-5\sin^{2}{I})\\ \quad\stackrel{{\scriptstyle}}{{\dot{\ell}}}=n=\displaystyle n_{0}\Bigl{[}1+\frac{3}{2}J_{2}\frac{R^{2}_{\oplus}}{p^{2}}\Bigl{(}1-\frac{3}{2}\sin^{2}I\Bigr{)}\sqrt{1-e^{2}}\Bigr{]}\end{cases}$$

$${\dot{\mathcal{E}}}_{kep}={\mathbf{X}}_{kep}({\mathcal{E}}_{kep}).$$

$${\dot{\mathcal{E}}}_{car}=\mathbf{Y}({\mathcal{E}}_{car}),$$

$$\displaystyle{\bf Y}=\Bigl{(}\frac{\partial\phi_{1}}{\partial\mathscr{E}_{kep}}{\bf X}_{kep}\Bigr{)}\circ\phi_{1}^{-1}$$

$$\ddot{{\bf r}}=\Bigl{(}\frac{\partial\dot{\bf r}}{\partial\mathscr{E}_{kep}}{\bf X}_{kep}\Bigr{)}\circ\phi_{1}^{-1}=:\tilde{\bf y}({\bf r},\dot{\bf r}).$$

$$\mathcal{K}=0,$$

$$\mathcal{K}=\Bigl{(}\ddot{\bf r}-\tilde{\bf y}\Bigr{)}\cdot{\bf e}^{\rho}=\ddot{\rho}-\rho\eta^{2}+\ddot{{\bf q}}\cdot{\bf e}^{\rho}-\tilde{\bf y}\cdot{\bf e}^{\rho},$$

$$|\dot{\mathbf{r}}{|}^2={\xi }^2+{\zeta }^2+2\dot{\mathbf{q}}\cdot {\mathbf{e}}^{\alpha }\xi +2\dot{\mathbf{q}}\cdot {\mathbf{e}}^{\delta }\zeta +|\dot{\rho }{\mathbf{e}}^{\rho }+\dot{\mathbf{q}}{|}^2,$$

$$\dot{\mathbf{r}}\cdot \mathbf{r}=\mathbf{q}\cdot {\mathbf{e}}^{\alpha }\xi +\mathbf{q}\cdot {\mathbf{e}}^{\delta }\zeta +(\dot{\rho }{\mathbf{e}}^{\rho }+\dot{\mathbf{q}})\cdot \mathbf{r}.$$

$$\mathbf{c}=\mathbf{A}\xi +\mathbf{B}\zeta +\mathbf{C},$$

$$\mathcal{E}=\frac{1}{2}|\dot{\mathbf{r}}{|}^2-\frac{\mu }{|\mathbf{r}|},$$

$$\displaystyle\hskip 28.45274pt\mu{\bf L}(\rho,\dot{\rho})=\dot{\bf r}\times{\bf c}-\mu\frac{{\bf r}}{|{\bf r}|}=\Bigl{(}|\dot{\bf r}|^{2}-\frac{\mu}{|{\bf r}|}\Bigr{)}{\bf r}-(\dot{\bf r}\cdot{\bf r})\dot{\bf r},$$

$$\mathbf{A}=\mathbf{r}\times {\mathbf{e}}^{\alpha },\mathbf{B}=\mathbf{r}\times {\mathbf{e}}^{\delta },\mathbf{C}=\mathbf{r}\times \dot{\mathbf{q}}+\dot{\rho }\mathbf{q}\times {\mathbf{e}}^{\rho }.$$

Full text

Keywords: orbit determination; Earth satellites; radar observations; Earth oblateness.

1 Introduction

We investigate a method of preliminary orbit determination that takes into account the Earth oblateness effect. The method is conceived for the computation of orbits of Earth satellites at low altitudes (LEO), with radar observations. This is an improvement of the algorithm described in [9], where a two-body approximation is employed. The Earth oblateness has been already included in the dynamical model for the computation of preliminary orbits, see e.g. [5], where radar observations were assumed to measure the range, the range rate, the right ascension and the declination. In this paper we include the secular effect of the $J_{2}$ term of the geopotential in the dynamical model, and we aim to compute preliminary orbits where also corrections to the values of the angular measurements are applied. To write the equations for this problem we shall use the Keplerian integrals evolution, the equations of motion projected onto the line of sight, and Lambert’s theorem. We assume that each radar measurement at epoch $t$ is composed by a precise value of the range $\rho$ (with standard deviation $\approx 10$ m) and poorly determined values of right ascension $\alpha$ and declination $\delta$ (with standard deviation $\approx 0.2$ deg), as we did in [9]. The available data are radar tracks of the form

[TABLE]

where $\Delta t=t_{j+1}-t_{j}\approx 10$ s. Given a radar track we can derive the vector

[TABLE]

where $\bar{t}$, $\bar{\alpha}$ and $\bar{\delta}$ are the mean values and $\rho$, $\dot{\rho}$, $\ddot{\rho}$ are obtained through a cubic fit.

To describe the osculating two-body orbit of the satellite we use spherical coordinates (known as attributable coordinates)

[TABLE]

Given the data in (2), we miss the values of $\dot{\alpha},\dot{\delta}$ to have an orbit, which are the unknowns of our orbit determination problem. Moreover, to improve the quality of the preliminary orbit we wish to correct the values of $\bar{\alpha},\bar{\delta}$. For this purpose we introduce the quantities $\Delta\alpha,\Delta\delta$, that are unknown small deviations from the mean values $\bar{\alpha}$, $\bar{\delta}$:

[TABLE]

We call infinitesimal angles the deviations $\Delta\alpha,\Delta\delta$ . Moreover, in place of the unknowns $\dot{\alpha},\dot{\delta}$ we use

[TABLE]

that are the components of the topocentric velocity of the satellite orthogonal to the line of sight. The orbit at time $\bar{t}-\rho/c$, where $c$ is the velocity of light, is completely determined by the coordinates

[TABLE]

that we call modified attributable coordinates. Our method consists in correlating two radar tracks to determine the values of the unknowns and compute preliminary orbits, see [10]. Therefore, the input data set is

[TABLE]

We want to determine the values of the vectors $(\Delta\alpha_{1},\Delta\delta_{1},\xi_{1},\zeta_{1})$ and $(\Delta\alpha_{2},\Delta\delta_{2},\xi_{2},\zeta_{2})$ at epochs $\bar{t}_{1}-\rho_{1}/c$ and $\bar{t}_{2}-\rho_{2}/c$.

2 Notation

Let us denote by ${\bf e}^{\rho}$ the unit vector corresponding to the line of sight, and by ${\bf q}$ the geocentric position of the observer. Then the geocentric position of the observed body is

[TABLE]

where $\rho$ is the range. Using as angular coordinates the topocentric right ascension $\alpha$ and declination $\delta$ in an equatorial reference frame (e.g. J2000), we have

[TABLE]

We introduce the unit vectors

[TABLE]

The set $\{{\bf e}^{\rho},{\bf e}^{\alpha},{\bf e}^{\delta}\}$ is an orthonormal system. Denoting by $\dot{\bf r}$ the geocentric velocity of the satellite, we have

[TABLE]

We will use the following different sets of coordinates for the orbits:

[TABLE]

that are respectively Keplerian, Cartesian, attributable and modified attributable coordinates.

We also consider the coordinate changes

[TABLE]

and the composite transformation

[TABLE]

from $\mathscr{E}_{kep}$ to $\mathscr{E}^{*}_{att}$.

3 The equations of motion

Let us consider Newton’s equation

[TABLE]

for the motion of a point mass in the Earth gravity field where the potential $U$ is truncated at the $J_{2}$-term, that is

[TABLE]

Here, $r=|{\bf r}|$ is the geocentric distance, $R_{\oplus}$ is the equatorial radius of the Earth and $P_{2}$ is the Legendre polynomial of second degree

[TABLE]

The problem defined by equation (7) is non-integrable, see [3]. If we average out the short period term in (7) we obtain an integrable system (see [11]) given by

[TABLE]

with $p=a(1-e^{2})$ the parameter of the two-body trajectory and $n_{0}=\sqrt{\mu/a^{3}}$ the mean motion. Note that, in the dynamics defined by (9), the elements $a,e,I$ remain constant while the ascending node $\Omega$, the argument of perigee $\omega$, and the mean anomaly $\ell$ change uniformly with time. Equations (9) can be written shortly as

[TABLE]

In the following we shall assume that the observed body is moving according to the integrable dynamics defined by equations (9), (10), and we shall call oblateness effect (or $J_{2}$ effect) the deviation from the pure Keplerian motion defined by these equations.

To solve our problem, we express the equations of motion in terms of the coordinates $\mathscr{E}^{*}_{att}$. First we express the equation (10) in Cartesian coordinates $\mathscr{E}_{car}=({\bf r},\dot{\bf r})$. We obtain

[TABLE]

where

[TABLE]

is the transformed vector field. From the expression above we obtain the acceleration $\ddot{\bf r}$ as a function of the Cartesian coordinates ${\cal E}_{car}$ along the solutions of (11):

[TABLE]

As done in [9] for the pure Keplerian dynamics, we project the perturbed equation of motion (12) along the line of sight ${\bf e}^{\rho}$ and obtain the equation

[TABLE]

with

[TABLE]

where $\eta=\sqrt{\dot{\alpha}^{2}\cos^{2}\delta+\dot{\delta}^{2}}$ is the proper motion.

Equation (13) can be expressed as a function of the unknown modified attributable coordinates $(\Delta\alpha,\Delta\delta,\xi,\zeta)$ using the expressions of ${\bf r},\dot{\bf r}$ given in (5), (6).

4 The $J_{2}$ effect on the two-body integrals

We recall the expressions of the conserved quantities in the Keplerian dynamics, i.e. the angular momentum ${\bf c}$, the energy $\mathcal{E}$ and the Laplace-Lenz vector ${\bf L}$, as a function of ${\bf r},\dot{\bf r}$. These quantities can be read as functions of the $\mathscr{E}_{att}$ attributable coordinates using (5), (6), and

[TABLE]

We have

[TABLE]

where

[TABLE]

Including the $J_{2}$ effect in the dynamics the angular momentum and the Laplace-Lenz vectors are not conserved anymore. However, the following relations hold true:

[TABLE]

where

[TABLE]

Here we denote by $R_{\phi}^{{\bf v}}$ the rotation of an angle $\phi$ around the axis defined by the vector ${\bf v}$. The angular variations $\Delta\Omega$ and $\Delta\omega$ are obtained as

[TABLE]

using equation (9) and the epochs $\tilde{t}_{i}=\bar{t}_{i}-\rho_{i}/c\quad i=1,2$, corrected by aberration, with $c$ the velocity of light.

Remark 1

We can also write

[TABLE]

and

[TABLE]

with

[TABLE]

5 Lambert’s theorem with the $J_{2}$ effect

Let us denote by $\mathcal{L}$ the expression defining Lambert’s equation. In the dynamics defined by (9), the mean motion evolves linearly, thus Lambert’s equation can be written as

[TABLE]

with $n$ given by the last equation in (9). Moreover, $k\in\mathbb{N}$ is the number of revolutions in the time interval $[\tilde{t}_{1},\tilde{t}_{2}]$. The angles $\beta$, $\gamma$ are defined by

[TABLE]

with $r_{1},r_{2}$ the distances from the center of force. In (18) the distance

[TABLE]

is the length of the chord joining the two positions of the body at epochs $\tilde{t}_{1},\tilde{t}_{2}$ after rotating the osculating ellipse at epoch $\tilde{t}_{1}$ so that it overlaps with the osculating ellipse at epoch $\tilde{t}_{2}$. The rotation $\tilde{R}$ is given explicitly by

[TABLE]

with $\tilde{R}_{1}$ and $\tilde{R}_{2}$ the transformations from the selected equatorial reference to the orbital reference frame at the epochs $\tilde{t}_{1}$ and $\tilde{t}_{2}$ respectively (${\bf\hat{n}_{1}}$ and ${\bf\hat{n}_{2}}$ are the directions of the lines of nodes):

[TABLE]

For a fixed number of revolutions $k$ we have 4 different choices for the pairs $(\beta,\gamma)$, see [9] and [1] for the details.

6 Linkage

We wish to link two sets of radar data of the form (2), with mean epochs $\bar{t}_{i}$, $i=1,2$, and compute one (or more) preliminary orbits. In the following we use labels 1, 2 for the quantities introduced in the previous sections, according to the epoch.

Moreover, let us define ${\bf v}_{2}={\bf e}^{\rho}_{2}\times{\bf q}_{2}$. Taking into account the $J_{2}$ effect we consider the system

[TABLE]

of 8 equations in the 8 unknowns $({\bf X},\bm{\Delta})$, with

[TABLE]

Note that the unknowns are divided into 2 sets so that $\bm{\Delta}$ is the vector of infinitesimal angles.

We search for solutions of equation

[TABLE]

where

[TABLE]

and ${\bf X}(\bm{\Delta})$ is a solution of

[TABLE]

In our scheme we use a double-iteration method. To solve system (20) we use the following procedure. Let us consider the function

[TABLE]

Solutions of (21) correspond to absolute minimum points of (23). To find them, we search for stationary points of ${\cal F}$, i.e. solutions of

[TABLE]

Applying Newton-Raphson method would yield the iterations

[TABLE]

where the Hessian matrix

[TABLE]

has components

[TABLE]

We drop the terms with the second derivatives of ${\cal G}$, that have quite complicated expressions. This variation of the algorithm is usually called pseudo-Newton method, and can converge if we start close to a solution. In fact, the second derivatives of ${\cal G}$ are multiplied by the components of ${\cal G}$, that are small in this situation.

At each iteration of the pseudo-Newton method we use again pseudo-Newton iterations to compute ${\bf X}(\bm{\Delta})$ from system (22). Precisely, for a given value of $\bm{\Delta}=\bm{\Delta}_{k}$ we search for ${\bf X}_{k}={\bf X}(\bm{\Delta}_{k})$ by looking for minimum points of the function

[TABLE]

where

[TABLE]

Taking advantage of the assumed smallness of the solutions $\bm{\Delta}$, we consider $\bm{\Delta}={\bf 0}$ as starting guess.

Remark 2

Equations (22) are not polynomial in ${\bf X}$, like the corresponding equations in [9].

7 Computing ${\bf X},\bm{\Delta}$

In equation (21) the components of the vector ${\bf G}({\bf X},\bm{\Delta})$ are similar to the ones of the corresponding vector in [9]. However, the following differences occur:

• i)

in $\mathcal{K}$ at epochs $\tilde{t}_{1},\tilde{t}_{2}$, the term $\tilde{\bf y}\cdot{\bf e}^{\rho}$ replaces the radial component of the Keplerian force $-\mu{\bf r}\cdot{\bf e}^{\rho}/|{\bf r}|^{3}$;

• ii)

in place of the Laplace-Lenz conservation law we have equation (16);

• iii)

in $\mathcal{L}$ we have a different expression of the mean motion $n$, as a function of the constant elements $a,e,I$, coming from the dynamical model (9). Moreover, we have a different expression for the length of the chord, see (19).

To search for the values of $\bm{\Delta}$ that solve equation (21) we have to compute the first derivatives of ${\cal G}(\bm{\Delta})$ with respect to $\bm{\Delta}$, appearing in (26). Since these computations are similar to the ones described in [9, Sect.7], we describe below only the differences coming from the adopted dynamical model.

7.1 The derivatives of $(R_{L}{\bf L}_{1}-{\bf L}_{2})\cdot{\bf v}_{2},\mathcal{L}$

Differentiating $(R_{L}{\bf L}_{1}-{\bf L}_{2})\cdot{\bf v}_{2}$ with respect to $\bm{\Delta}$ we have to compute the derivatives of the rotation $R_{L}$. For this purpose we use Euler-Rodrigues formula for the rotation matrices (see [6])

[TABLE]

where $\hat{{\bf v}}$ is the skew-symmetric matrix associated to ${\bf v}=(v_{1},v_{2},v_{3})$ by

[TABLE]

that is

[TABLE]

We have

[TABLE]

The derivatives of $\mathcal{L}$ differ from the ones written in ([9]) only for the different expressions of $n$ and $d_{L}$. We have

[TABLE]

[TABLE]

7.2 The derivatives of $\mathcal{K}_{1}$, $\mathcal{K}_{2}$

We also need to compute the following derivatives:

[TABLE]

The derivatives of

[TABLE]

are given by

[TABLE]

with

[TABLE]

In (29) we made a little abuse of notation: ${\bf e}^{\rho}$ stands both for a function of $(\alpha,\delta)$ and for a function of $(\Delta\alpha,\Delta\delta)$.

Then, we introduce ${\bf y^{*}}$ that is the vector $\tilde{\bf y}$ as a function of the $\mathscr{E}^{*}_{att}$ coordinates:

[TABLE]

The components of ${\bf y^{*}}$ are

[TABLE]

and their derivatives are given by

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

Denoting by $(\dot{\Omega},\dot{\omega},\dot{\ell})$ the last three components of the vector field ${\bf X}_{kep}$, see (9), we have

[TABLE]

where ${\cal O}_{3}$ is the null matrix of order 3 and the components of $\frac{\partial(\dot{\Omega},\dot{\omega},\dot{\ell})}{\partial(a,e,I)}$ are given by

[TABLE]

with

[TABLE]

and

[TABLE]

Acknowledgments

This work is partially supported by the Marie Curie Initial Training Network Stardust, FP7-PEOPLE-2012-ITN, Grant Agreement 317185.