Robust Positioning, Preliminary Orbit Determination, and Trajectory Prediction of Space Debris using In-Space Iterative-Bearing-Only Observations

Abstract

This paper is concerned with the preliminary localisation, orbit determination and model-based path forecasting of space debris based on a robust procedure. In this work, an in-orbit observer utilises only relative bearing observations iteratively. To this end, the problem is first formulated in order to calculate the distance vector between the space debris and any orbiting observer. Afterwards, the obtained position vector is corrected through an Extended Kalman Filter (EKF) for shrinking the sensor and process errors and increasing robustness of the computations in the presence of uncertainties. After preliminary positioning, the related classical orbital elements are acquired via the predicted position and velocity vectors using a hybrid technique. Extensive simulations demonstrate the efficacy and robustness of the aforementioned method, and in particular it is verified that the proposed scheme is capable of producing a suitable solution for preliminary localisation and orbit determination of space debris based on the presented space-based observation, which is practical in phasing and chasing manoeuvres of any grabber space robot.

1. INTRODUCTION

Space debris are man-made objects, which have existed in space since the first launch of artificial satellites. This includes satellites whose operational lifespan is over, components of rockets and other objects. Space debris is a major in-space operational hazard for both current and future satellites and manned missions. Therefore, it is rapidly becoming necessary to consider cleaning up potential orbits from any debris. The first stage of this is to precisely detect the orbits of debris (Flury, 1995). Orbit determination (Armellin et al., 2016; Curtis, 2013) is germane to one of the branches of astronomy, which observes, calculates and predicts space objects' orbits. One common application of orbit determination is supporting Global Positioning System (GPS) satellites (Montenbruck et al., 2005) Andres Johan Lexell and Friedrich Gauss initiated modern orbit determination (Sten and Aspaas, 2013). Furthermore, Gauss introduced a technique through which six orbital elements are obtained through only three observations at three different times (Chung, 2006). For orbit determination, using Gauss's Method, right ascension and declination should be specified during three different times. Moreover, local sidereal time is also specified at the observation times. This technique is specific to short arcs; in other words, the time difference between two observations should be shortened (Fadrique et al., 2012). One of the other methods for Angles-only Orbit Determination is the Laplace Method, which is a standard technique for initial orbit calculations. This method was used for the first time for the orbit determination of a meteor named P/1846 D1. The Laplace method is similar to the Gauss technique (Richard, 2003).

Another method for orbit determination is Gooding's (1993) method that considered various issues, the first of which is about the best position and time for observation of space objects. The second point pertaining to his work is utilisation of auxiliary sensory tools with the purpose of improving the observation conditions. Moreover, the effects of sensor type and their accuracy for initial distance determination were considered in the method. Further methods incorporate a myriad of corrections on the Laplace and Gauss methods. Briggs and Slowley (1959) presented a technique through which an iterative method was introduced utilising computer technologies. Escobal (1965) has presented an approach for angles-only orbit determination, namely the Double Iteration Method. This technique is similar to the classical orbit determination methods; however, there are some differences between them. For instance, Escobal's approach is far more accurate than previous classical methods, but compared to the previous ones, the method suffers from long computational time. Vallado (2010) has recently considered the three basic techniques including Gauss, Laplace and Escobal, and proposed these mentioned methods as the standard orbit determination techniques. All the Gauss, Laplace, Gooding and Escobal methods focus on angles. The Gooding method is an extended iterative technique that is specific to earth-orbiting satellites, where their angular data are based on several rotations around the earth. This method is not compatible with our objectives in this work. On the other hand, Escobal's technique, that is Double Iteration, is suitable for long arcs.

The Laplace approach will be optimal only if the distance between the grabber satellite and space debris is extremely long. Therefore, the best method for short arcs, which is the main objective of this paper, is the Gauss approach (Marsden, 1985). In addition, Xu et al. (2016) achieved some rendezvous by the proposed adaptive control based on the measurements of relative position and velocity between the chaser and target spacecraft. Gruntman (2014) researched some methods that use optical methods for detecting sub millimetre and millimetre size debris in low earth orbits. In addition, some further research on extensions of these methods has also been conducted in Armellin et al. (2016) and Kaufman et al. (2016) for further analysis of the uncertainties. Figure 1 exhibits the proposed plan of robust in-orbit localisation and orbit prediction of the space debris.

Figure 1. Diagram of robust in-orbit positioning, orbit determination, and path prediction.

2. ITERATIVE BEARING-ONLY POSITIONING OF SPACE DEBRIS

All the orbit determination methods are based on a common basis (Townsend, 2008). If $\vec{\bi r}$ is position vector of the space debris, $\vec{\bi{\rho}}$ is the distance vector between the space debris and grabber satellite and $\vec{\bi R}$ is position vector of grabber satellite, it leads to:

(1) $$\vec{r}=\vec{\rho }+\vec{R}=\rho \hat{L}+\vec{R}$$

where $\hat{L}$ is direction unit vector between the space debris and grabber satellite. A schematic view of the problem's geometry for orbit determination of space debris is displayed in Figure 2 For calculating the position vector of the space debris, it is necessary to present some assumptions in this paper. Since the main goal of this work is utilisation of a low-cost orbit determination based on the observations from an orbiting robot; this work is limited to utilising bearing-only methods. In other words, since finding the slant range between the space debris and grabber satellite needs expensive equipment, the Gauss technique, which is one of the bearing-only methods, is employed in this research (Curtis, 2013). Therefore, the azimuth and elevation angles between the observer and debris are specified. Note that the following angles are related to three different times over a short arc. For orbit determination using the Gauss method, the obtained azimuth and elevation angles are first converted into right ascension and declination, which are shown in Figure 3.

Figure 2. Schematic view of the geometry of orbit determination for a sample space debris.

Figure 3. Schematic view of the grabber satellite's position.

For this purpose, the unit vector of the sensor coordinate system, which is mounted on the viewer, is first obtained for transformation from sensor coordinate system to the right ascension and declination coordinate system as follows:

(2) $$\left[{\matrix{L_{hx} \cr L_{hy} \cr L_{hz}}} \right]=\left[{\matrix {-\cos {(AZ)}\cos {(El)} \cr \sin {(AZ)}\cos {(El)} \cr \sin {(El)}}} \right]$$

where L _h , AZ and EL are line of sight unit vectors between the space debris and grabber satellite, azimuth, and elevation angles, respectively. After calculating the line of sight unit vector, it is transformed to the right ascension and declination frame as noted below:

(3) $$\left[{\matrix {L_{hx} \cr L_{hy} \cr L_{hz} \cr}} \right] =\left[{\matrix {\sin {(\phi )}\cos {({\Omega })} & -\sin {({\Omega })} & \cos {(\phi )}\cos {({\Omega })} \cr \sin {(\phi )}\sin {({\Omega })} & \cos {({\Omega })} & \cos {(\phi )}\sin {({\Omega })} \cr -\cos {(\phi )} & 0 & \sin {(\phi )}}} \right]\!. \left[{\matrix {L_{hi} \cr L_{hj} \cr L_{hk} \cr}} \right]$$

(4) $$\Omega = cos^{-1}(\displaystyle{{\vec{R}_{ij}\kern-1pt{.}\,\vec{X}}\over{\mid\vec{R}_{ij} \parallel\vec{X}\mid}})$$

(5) $$\phi = cos^{-1}\left(\displaystyle{{\vec{R}_{ij}\kern-1pt{.}\,\vec{R}}\over{\mid\vec{R}_{ij} \parallel\vec{R}\mid}}\right)$$

where Ω, ϕ are right ascension and declination angles, respectively. Also $\vec{\hbox{X}}$ and $\vec{\hbox{R}}$ are x-direction in the Earth-Centred-Inertial (ECI) frame and position vector of observer, respectively. There are numerical non-zero values including a, b and c which are capable of creating the following relationship:

(6) $${a\vec{r}_{1}}+{b\vec{r}_{2}}+{c\vec{r}_{3}}=0$$

where $\vec{r}_{1}$ , $\vec{r}_{2}$ and $\vec{r}_{3}$ are position vector of space debris during the first, second and third observations. If the aforementioned relationship is solved in terms of ${\vec{r}}_{2}$ , the following relationship will be created:

(7) $$\vec{r}_{2}=c_{1}\vec{r}_{1}+c_{3}\vec{r}_{3}$$

For calculating the r ₁ and r ₃, the Lagrange coefficients, f and g are required through which the position vector for a specified time could be expressed in terms of position and velocity vectors during other times as noted below:

(8-1) $$\vec{r}_{1}=f_{1}\vec{r}_{2}+g_{1}\vec{v}_{2}$$

(8-2) $$\vec{r}_{3}=f_{3}\vec{r}_{2}+g_{3}\vec{v}_{2}$$

(8-3) $$f_{1}=1-\displaystyle{{\lpar t_{1}-t_{2} \rpar ^{2}}\over{2}}h_{2}$$

(8-4) $$f_{3}=1-\displaystyle{{\lpar t_{3}-t_{2} \rpar ^{2}}\over{2}}h_{2}$$

(8-5) $$g_{1}=\lpar t_{1}-t_{2} \rpar -\displaystyle{{\lpar t_{1}-t_{2} \rpar ^{3}}\over{6}}h_{2}$$

(8-6) $$g_{3}=\lpar t_{3}-t_{2} \rpar -\displaystyle{{\lpar t_{3}-t_{2} \rpar ^{3}}\over{6}}h_{2}$$

where $\vec{v}_{2}$ is the velocity vector of space debris during the second observation. Therefore, $\vec{r}_{1}$ and $\vec{r}_{3}$ are calculated in terms of $\vec{r}_{2}$ . In addition, h ₂ is angular momentum of space debris during the second observation and μ is the earth gravitational parameter. By combination the equations (Curtis, 2013), we have:

(9-1) $${c_{2}\rho }_{2}\hat{L}_{2}+c_{1}\rho_{1}\hat{L}_{1}+c_{3}\rho_{3}\hat{L}_{3}=-c_{1}\vec{r}_{site1}-c_{3}\vec{r}_{site3}-c_{2}\vec{r}_{site2}\equiv M$$

(9-2) $$\left[\matrix{L_{1x} & L_{2x} & L_{3x} \cr L_{1y} & L_{2y} & L_{3y} \cr L_{1z} & L_{2z} & L_{3z}} \right]\!.\left[\matrix{ c_{1}\rho_{1}\cr {c_{2}\rho }_{2}\cr c_{3}\rho_{3}} \right]=M=[L] \left[\matrix{ c_{1}\rho_{1}\cr {c_{2}\rho }_{2}\cr c_{3}\rho_{3}} \right]$$

where c ₁ and c ₃ are functions of $\vert \vec{r}_{2} \vert $ , t ₁, t ₂ and t ₃. Moreover, $\hat{L}$ ₁, $\hat{L}$ ₂ and $\hat{L}$ ₃ are direction unit vectors between the debris and observer during the first, second and third observations. By inversing, the unit vectors of the line of sight matrix are:

(10-1) $$\left[\matrix{c_{1}\rho_{1} \cr {c_{2}\rho }_{2} \cr c_{3}\rho_{3}}\right]=[L]^{-1}M$$

(10-2) $$[L]^{-1}\equiv [A] =\left[\matrix{ a_{11} & a_{12} & a_{13} \cr a_{21} & a_{22} & a_{23} \cr a_{31} & a_{32} & a_{33}} \right]$$

Afterwards, in terms of new variables the following relations are revealed:

(11-1) $$c_{1}\approx {\tau_{1} \over \tau_{2}}\left[ 1+{\lpar \tau_{2}^{2}-\tau_{1}^{2} \rpar \over {6}}h_{2} \right]=A_{1}+B_{1}h_{2}$$

(11-2) $$c_{3}\approx {\tau_{3}\over \tau_{2}}\left[ 1+{\lpar \tau_{2}^{2}-\tau_{3}^{2} \rpar \over {6}}h_{2} \right]=A_{3}+B_{3}h_{2}$$

(11-3) $$\vec{A}^{T}=\left[\matrix{A_{1} & -1 & A_{3} } \right]$$

(11-4) $$\vec{B}^{T}=\left[ \matrix{B_{1} & 0 & B_{3}} \right]$$

For obtaining ρ₁, ρ₂ and ρ₃ the following relations are exploited:

(12-1) $$\rho_{2}=A_{2}^{\ast }+B_{2}^{\ast }h_{2}$$

(12-2) $$\rho_{1}={A_{1}^{\ast }+B_{1}^{\ast }h_{2}\over c_{1}}$$

(12-3) $$\rho_{3}={A_{3}^{\ast }+B_{3}^{\ast }h_{2}\over c_{3}}$$

(12-4) $$A_{2}^{\ast }=a_{21}\vec{A}{.}\,\vec{X}+a_{22}\vec{A}{.}\,\vec{Y}+a_{23}\vec{A}{.}\,\vec{Z}$$

(12-5) $$B_{2}^{\ast }=a_{21}\vec{B}{.}\,\vec{X}+a_{22}\vec{B}{.}\,\vec{Y}+a_{23}\vec{B}{.}\,\vec{Z}$$

(12-6) $$A_{1}^{\ast }=-(a_{11}\vec{A}{.}\,\vec{X}+a_{12}\vec{A}{.}\,\vec{Y}+a_{13}\vec{A}{.}\,\vec{Z})$$

(12-7) $$B_{1}^{\ast }=-(a_{11}\vec{B}{.}\,\vec{X}+a_{12}\vec{B}{.}\,\vec{Y}+a_{13}\vec{B}{.}\,\vec{Z})$$

(12-8) $$A_{3}^{\ast }=-(a_{31}\vec{A}{.}\,\vec{X}+a_{32}\vec{A}{.}\,\vec{Y}+a_{33}\vec{A}{.}\,\vec{Z})$$

(12-9) $$B_{3}^{\ast }=-(a_{31}\vec{B}{.}\,\vec{X}+a_{32}\vec{B}{.}\,\vec{Y}+a_{33}\vec{B}{.}\,\vec{Z})$$

Moreover, position vectors of the observer during three observations are expressed as noted below:

(13) $$[X \vdots Y \vdots Z]^{T}=\left[ \matrix{ \vec{r}_{site1x} & \vec{r}_{site2x} & \vec{r}_{site3x}\cr \vec{r}_{site1y} & \vec{r}_{site2y} & \vec{r}_{site3y}\cr \vec{r}_{site1z} & \vec{r}_{site2z} & \vec{r}_{site3z}} \right]$$

where $\vec{X}$ , $\vec{Y}$ and $\vec{Z}$ are three directions in the ECI frame. Finally, for computation of the position vector of space debris during the second observation, the following equation is solved:

(14) $$r_{2}^{8}-r_{2}^{6}[r_{site2}^{2}+A_{2}^{\ast 2}+2\lpar \hat{L}_{2}.\,\vec{r}_{site2}\rpar A_{2}^{\ast } \rsqb -r_{2}^{3}[2A_{2}^{\ast }B_{2}^{\ast } \mu +2\lpar \hat{L}_{2}.\, \vec{r}_{site2} \rpar \!B_{2}^{\ast } \mu \rsqb -B_{2}^{\ast 2} \mu^{2}=0$$

By solving and obtaining eight possible answers, two cases occur. In the first case, if the space debris during the observation time lies in higher orbit relative to the satellite, that is $\vert \vec{r}_{2} \vert >\vert \vec{r}_{site2} \vert $ , the largest root is the desired solution. In the second case, if the satellite lies in higher orbit during the observation time relative to the space debris, that is $\vert \vec{r}_{2} \vert <\vert \vec{r}_{site2} \vert $ , the second-largest root is the solution. Therefore, the position vector of the debris is eventually obtained as an input for the EKF for increasing the robustness and accuracy of the Gauss Method, which is described in detail in the next Section.

3. EKF AND ROBUST POSITIONING

In this work, the dynamics of the system are nonlinear. Hence, both the EKF and Unscented Kalman Filter (UKF) may be used. UKF is similar to the EKF but the linearization process is not applied in UKF. Because EKF is one of the most widely used filters for nonlinear systems, specifically because of its speed, this filter is employed in the current research (Giannitrapani et al., 2011). In fact, this filter converts the nonlinear system into a state-dependent time varying linear system; however, this transformation leads to some modelling errors. The significance of the Kalman filter is that all calculations for the system's state estimations (regardless of accuracies) are done using proper weighting. Relative motion between the earth and satellite is a two-body problem without any noise according to the law of universal gravitation.

3.1. Filtering and Estimation

Now, EKF is considered in detail. This filter is undoubtedly one of the most widely used filters for Gaussian nonlinear systems over the past few decades. By employing this filter, the system is first linearized through sampling times and then, systems' states are estimated for the next steps. Generally, the solving process for EKF is divided into two parts including predicting and updating times. During the estimation step, the measurement error covariance is first considered as follows (Wu et al., 2013; Kaufman et al., 2016; Armellin et al., 2016):

(15) $$p_{k}^{-}=F_{k-1}p_{k-1}^{+}F_{k-1}^{T}+Q_{k-1}$$

where plus and minus indices are germane to the next and previous times, respectively. State is estimated as noted below:

(16) $$\hat{x}_{k}^{-} = f_{k-1}(\hat{x}_{k-1}^{+}u_{k-1},0)$$

It should be noted that a fourth-order Runge-Kutta method has been utilised in this paper for calculating the $\hat{x}_{k}^{-}$ . Afterwards, the updating step is started by computation of the Kalman gain matrix as:

(17) $$k_{k}=p_{k}^{-}H_{k}^{T}{{(}H_{k}\ p_{k}^{-}\ H_{k}^{T})}^{-1}$$

where k _k is the gain matrix. State variables and covariance matrices during the update process are calculated by means of observation data as follows:

(18-1) $$\hat{x}_{k}^{+}=\hat{x}_{k}^{-}+k_{k}[y_{k}-h_{k}\lpar \hat{x}_{k}^{-},0 \rpar \rsqb $$

(18-2) $$p_{k}^{+}=(I-k_{k}H_{k}^{T})p_{k}^{-}$$

3.2. Range Prediction Model

A diagramatic snapshot of orbits of observer space robot and space debris is shown in Figure 4. Here the state vector of the system is defined as (Wu et al., 2013):

Figure 4. Supposed orbits of the grabber satellite and space debris.

(19) $$X = {({r_x},\,{r_y},\,{r_z},\,{\upsilon_x},\,{\upsilon_y},\,{\upsilon_z})}^T$$

Equations of motion in a Keplerian Orbit, comprising perturbations based on the law of universal gravitation are expressed as:

(20) $$\ddot r =-\displaystyle{\mu\over {r^2}}{\hat{\bi e}}_r + {\bi F}_p$$

where μ is Earth gravitational parameter and ${\bi F}_{\bi p}$ model uncertainties due to the perturbations. Moreover, the motion state differential equation is defined as:

(21) $$\dot{X_{k}}= \lpar \dot{r_{x}},\dot{r_{y}},\dot{r_{z}},\ddot{r_{x}},\ddot{r_{y}},\ddot{r_{z}} \rpar _{k}^{T}=f(X_{k}, k)$$

The second-order Taylor expansion for discrete form of the aforementioned relationship is:

(22) $$X_{k+1}=X_{k}+f\!\lpar X_{k} \rpar \!h+F\!\lpar X_{k} \rpar \!f\!\lpar X_{k} \rpar \!{h^{2}\over{2}}+O^{2}{}_{k}\lpar h \rpar =\phi \lpar X_{k} \rpar +O^{2}{}_{k} \lpar h \rpar $$

where $O^{2}{}_{k}\lpar h \rpar $ is model uncertainty (Big-O-2) which is modelled as a Gaussian white noise with zero average and $F\!\lpar X_{k} \rpar $ , Jacobian of F, is defined as below:

(23) $$F\!\lpar X_{k} \rpar ={\partial f\!(X_{k}) \over \partial X_{k}}$$

State vectors of the satellite and space debris during k-th step are respectively calculated as follows:

(24-1) $$\dot{x}_{s,k}=\lpar \dot{r}_{xs},\dot{r}_{ys},\dot{r}_{zs},\ddot{r}_{xs},\ddot{r}_{ys},\ddot{r}_{zs} \rpar _{k}^{T}$$

(24-2) $$\dot{x}_{b,k}=\lpar \dot{r}_{xb},\dot{r}_{yb},\dot{r}_{zb},\ddot{r}_{xb},\ddot{r}_{yb},\ddot{r}_{zb} \rpar _{k}^{T}$$

Consequently, the Jacobian matrix of the state equation is derived as:

(25) $${F(x}_{d})=\left( \matrix{ 0 & 0 & 0 & 1 & 0 & 0 \cr 0 & 0 & 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & 0 & 0 & 1 \cr -{\mu \left( r^{{3 \over 2}}-3x_{b}^{2}r^{{1 \over 2}} \right)\over r^{3}} & -3\mu x_{b}y_{b}r^{-{5\over 2}} & -3\mu x_{b}z_{b}\lpar k \rpar r^{-{5 \over 2}} & 0 & 0 & 0 \cr -3\mu x_{b}y_{b}r^{-{5 \over 2}} & -{\mu \left( r^{3\over 2}-3y_{b}^{2}r^{1\over 2} \right) \over r^{3}} & -3\mu z_{b}y_{b}r^{-{5 \over 2}} & 0 & 0 & 0 \cr -3\mu x_{b}z_{b}r^{-{5 \over 2}} & -3\mu y_{b}z_{b}r^{-{5 \over 2}} & -{\mu \left( r^{{3 \over 2}}-3z_{b}^{2}r^{1\over 2} \right)\over r^{3}} & 0 & 0 & 0} \right); at\ k^{th} step$$

where the exhibited elements are the resultant of the Jacobian operator on f. In addition, r denotes the squared range from the debris and is defined as noted below:

(26) $$r^{2}=r_{xb}^{2}+r_{xb}^{2}+r_{xb}^{2}$$

3.3. Bearing-Only Observation Model

In this work, measurement equation is considered as below:

(27) $$Y=\left[ \matrix{ AZ \cr EL} \right]=h\lpar X_{s},X_{b} \rpar +V$$

where V resembles measurement noises. Moreover, the measured azimuth and elevation angles through mounted sensors on the observer robot are expressed as follows:

(28) $$AZ=\tan^{-1}{\left({r_{yb}-r_{ys}\over r_{xb}-r_{xs}}\right)}$$

(29) $$EL=\left({r_{zb}-r_{zs}\over\sqrt {\lpar r_{xb}-r_{xs} \rpar^{2}+\lpar r_{yb}-r_{ys} \rpar^{2}}}\right)$$

In the following, the Jacobian matrix of the measurement equation is derived as noted below:

(30) $$\eqalign{H &= \displaystyle{{\partial h} \over {\partial X_b}};at\;k^{th}\;step \cr & = \left[ {\matrix{{ - \displaystyle{{r_{y,bs}} \over {r_{x,bs}^2 + r_{y,bs}^2 }}} & {\displaystyle{{r_{x,bs}} \over {r_{x,bs}^2 + r_{y,bs}^2 }}} & 0 & 0 & 0 & 0 \cr { - \displaystyle{{r_{z,bs}r_{x,bs}} \over {{\left( {r_{x,bs}^2 + r_{y,bs}^2 } \right)}^{{\textstyle{1 \over 2}}}d^2}}} & { - \displaystyle{{r_{z,bs}r_{y,bs}} \over {{\left( {r_{x,bs}^2 + r_y^2 } \right)}^{{\textstyle{1 \over 2}}}d^2}}} & {\displaystyle{{{\left( {r_{x,bs}^2 + r_{y,bs}^2 } \right)}^{{\textstyle{1 \over 2}}}} \over {d^2}}} & 0 & 0 & 0}} \right]} $$

in which r _x , r _y , r _z are elements of the range vector along x, y, and z directions while d is magnitude of the range vector at k ^th step:

(31-1) $$r_{x,bs}=r_{xb}-r_{xs}$$

(31-2) $$r_{y,bs}=r_{yb}-r_{ys}$$

(31-3) $$r_{z,bs}=r_{zb}-r_{zs}$$

(31-4) $$d^{2}= r^{2}_{x,bs} + r^{2}_{y,bs} + r^{2}_{z,bs}$$

4. SOLUTION OF LAMBERT PROBLEM AND ORBIT DEFINITION

The first stage in cleaning space debris from an orbit is determination of the debris' orbit through a space/ground observer. Since the position of the space debris is specified during a bearing-only process, the velocity vector should be calculated through some techniques for accomplishment of the requirements of orbit determination procedure. Since the Lambert approach (Curtis, 2013) is far more accurate than the Gibbs Method it is utilised in our preliminary orbit determination process. It is indicated that the Lambert Problem has a shorter computation time relative to the Gibbs one. In this method, orbital elements are calculated based on only two position observations and the time difference between them. It is supposed that the calculated position vectors during two different times are ${\bi r}_{\bf{1}}$ and ${\bi r}_{\bf{2}}$ and the time difference between them is indicated. Afterwards, it is determined whether the orbit is retrograde or prograde. The angle between the two position vectors in the Lambert Problem is calculated as follows:

(32) $$\alpha ={\cos^{-1}}\displaystyle{{r_1 \cdot r_2} \over{r_1r_2}}$$

Some auxiliary variables are employed for the sake of simplification, one of which is defined as:

(33) $$A = \sin \Delta \alpha \sqrt {\displaystyle{{r_1r_2} \over {1 - \cos \Delta \alpha}}} $$

Correspondingly, another auxiliary variable z is defined and obtained iteratively from numerical methods as follows:

(34) $$z_{i+1}=z_{i}- {F(z_{i})\over \acute{F}(z_{i})}$$

The mentioned equation is a recursive one and the calculation process is continued until the value of z approaches its amount in the prior step. The auxiliary function of $y\lpar z \rpar $ is obtained from the following functions:

(35-1) $$y(z)= r_{1}+r_{2}+A{zS\lpar z \rpar -1 \over \sqrt {C(z)}}$$

(35-2) $$S(z)=\left\{\eqalign{{\sqrt Z -\sin \sqrt Z \over \sqrt Z^{3}} \quad \quad\quad (z>0)\cr {\sinh \sqrt {-Z} -\sqrt {-Z} \over \sqrt {-Z}^{3}} \quad (z< 0)\cr {1\over6}\quad\quad\quad\quad\quad\quad\quad\quad(z=0)}\right.$$

(35-3) $$C(z)= \left\{\eqalign{{1-\cos \sqrt Z \over \sqrt Z^{3}} \quad \quad\quad (z>0) \cr {\cosh \sqrt {-Z} -1 \over \sqrt {-Z}^{3}} \quad \quad (z< 0) \cr {1\over 2} \quad \quad\quad\quad\quad\quad\quad(z=0)}\right.$$

where C(z) and $S\lpar z \rpar $ are called Stumpff functions, Afterwards, Lagrange coefficients and their derivatives are obtained as noted below:

(36-1) $$ f=1- {y\over r_{1}}$$

(36-2) $$g = A\sqrt {\displaystyle{y \over \mu}} $$

(36-3) $$\dot{g}= 1- {y\over r_{2}}$$

Eventually, velocity vectors at the two observation times are calculated as follows:

(37-1) $$v_{1}= {1\over g}(r_{2}-fr_{1})$$

(37-2) $$v_{2}=\ {1\over g}({\dot{g}r}_{2}-r_{1})$$

As noted above, the six classical orbital elements are obtained from a union of the precomputed position and velocity vectors. It is mentioned that the Lambert Approach is far easier than the Gibbs Method. To elucidate further, the space object in the Gibbs Method is observed three times, whereas two observations as well as the time difference between them are sufficient in the Lambert Approach. Note that as per our experiments, the required solving time in the Gibbs Method is nearly five times greater than in the Lambert Problem. After computing the position and velocity vectors, angular momentum, inclination angle, the eccentricity vector is obtained as noted below:

(38-1) $$h=r \times v$$

(38-2) $$i= \cos^{-1}{h_{z}\over h}$$

(38-3) $$e= {1\over \mu }\left[ \left(v^{2}-{\mu \over r} \right)r-rv_{r}v \right]$$

(38-4) $$v_{r}= {v\,{.}\,r \over r}$$

Afterwards, right ascension, argument of perigee, and the true anomaly angle are also attained as follows:

(39-1) $$\Omega = \cos^{-1}{N_{x} \over N}$$

(39-2) $$w = \cos^{-1}{N\,{.}\,e \over Ne}$$

(39-3) $$\theta = \cos^{-1}{r\,{.}\,e \over re}$$

(39-4) $$N=\hat{K} \times h$$

where $\hat{K}$ is a unit vector along z. Moreover, perigee radius, apogee radius and semi-major axis are calculated from relations exhibited in Curtis (2013).

5. MODEL-BASED ORBIT PREDICTION AND MECHANISATION

Prediction of space debris' orbit, during the observations, is an essential task; due to the perturbations such as the earth anomalies like oblateness and bulge, atmospheric and geomagnetic forces and moments, solar and galactic winds and so on. Any geocentric orbit suffers from dynamical changes of orbital elements. This lack of symmetry means that the force of gravity on a revolving object is not focused along the centre of the earth. While the gravitational field of a plainly spherical sphere depends on the distance from its centre, oblateness causes a variation with latitude, which is called zonal variation, as well. The second zonal harmonic $\lpar {J}_{{2}}\rpar \comma \; $ is a dimensionless factor, which indicates the effects of oblateness on orbits. In the following, the gravitational specific force (per unit mass) in an oblate globe is given by:

(40) $$\ddot r = -\displaystyle{\mu\over{r^2}}{\hat{\bi e}_{\bi r}} + {\bi F}_{\bi p}$$

in which ${\bi F}_{\bi p}$ stands for the perturbing acceleration that is due to the earth spherical defects. This perturbing acceleration is comprised of the following noted components:

(41) $${\bi F}_{\bi p}=p_{r}\hat{\bi e}_{\bi r} + p_{\bot}\hat{\bi e}_{{\bot }} + p_{\bi h}\hat{\bi h}$$

where $\hat{\bi e}_{r}\comma \; \hat{\bi e}_{\bot }\comma \; \hat{h}$ are respectively, radial, transverse and perpendicular unit vectors attached to the satellite. Furthermore, p _r , p _⊥and p _h all depend on j ₂ that are computed from the following relations:

(42-1) $$p_{r} = -{\mu \over r^{2}}{3\over 2}J_{2}\left( {R\over r} \right)^{2}[{1-3}{sin}^{2}i\ {\sin}^{2}(\omega +\theta)]$$

(42-2) $$p_{\bot } = -{\mu \over r^{2}}{3\over 2}J_{2}\left( {R\over r} \right)^{2}{sin}^{2}i\ {\sin}[2(\omega +\theta)]$$

(42-3) $$p_{\bf h} = -{\mu \over r^{{2}}}{{3}\over{2}}J_{{2}}\left( {R\over r} \right)^{{2}}{[1-3}{sin}{2}isin{(}\omega + \theta{)]}$$

Prussing and Conway (1993) showed how p _r , p _⊥ and p _h affect time rates of change for Right Ascension and Argument of Perigee, respectively.

(43) $$\dot{\Omega }={h\over \mu }{{\sin}(\omega +\theta )\over \sin {i\ (1+ecos\theta )}}p_{\bf{h}}$$

(44) $$\dot{\omega }=-{rcos\theta \over eh} p_{r}+{(2+ecos\theta )sin\theta \over eh}p_{\bot }-{rsin(\omega +\theta )\over h\tan i}p_{\bf{h}}$$

Evidently, the time variation of the Right Ascension depends merely on the component of the instantaneous perturbing force, normal to the orbital plane, while the rate of change of the Argument of Perigee is subject to all the perturbation components.

5.1. Mechanisation in Inertial and Earth-Fixed Frames

To compute orbital position and velocity ( $X\comma \; Y\comma \; Z\comma \; \dot{X}\comma \; \dot{Y}\comma \; \ \hbox{and}\ \dot{Z}$ ) at time t, given, $a\comma \; e\comma \; i\comma \; \omega\comma \; \Omega\comma \; \ \hbox{and}\ T$ , at the first step the mean anomaly, MA, the eccentric anomaly, EA, and the true anomaly, ν, are computed respectively using the following relations (Chobotov, 2002):

(45-1) $$M =\sqrt{{\mu\over a^{3}}}(t-T)$$

(45-2) $$MA =EA -e\sin EA$$

(45-3) $$\tan{\nu \over 2} = \left[{(1+e)\over (1-e)}\right]^{1/2} \tan{EA \over 2}$$

Afterwards, the radius, r, the specific angular momentum, h, the position and the velocity vectors' components, are computed as below:

(46-1) $$r= a(1-e\cos EA) ={p\over 1+e\cos \nu} = {a(1-e^{2})\over 1+e \cos \nu}, \hbox{h} = (\mu a(1-e^{2}))^{1/2}$$

$$X = r(\cos\Omega \cos(\omega + \nu) - \sin \Omega \sin(\omega + \nu) \cos i)$$

(46-2) $$Y = r(\sin\Omega \sin(\omega + \nu) - \cos \Omega \sin(\omega + \nu) \cos i)$$

$$Z = r(\sin i \sin(\omega + \nu))$$

$$\dot{\hbox{X}} = {\hbox{Xhe}\over \hbox{p}} \sin \nu- {\hbox{h}\over\hbox{r}}(\cos\Omega \sin(\omega + \nu) + \sin\Omega \cos(\omega + \nu)\cos \hbox{i})$$

(46-3) $$\dot{\hbox{Y}} = {\hbox{Yhe}\over\hbox{p}} \sin \nu -{\hbox{h}\over \hbox{r}}(\sin\Omega \sin(\omega + \nu) - \cos\Omega \cos(\omega + \nu)\cos \hbox{i})$$

$$\dot{\hbox{Z}} = {\hbox{Zhe}\over \hbox{p}} \sin \nu +{\hbox{h}\over \hbox{r}}\sin \hbox{i}\ \cos(\omega + \nu)$$

To resolve coordinates in ECEF (Earth-Centred-Earth-Fixed) frame, the following relation is exploited:

(47) $$\left[{\matrix{X \cr Y \cr Z \cr}}\right]_{ECEF} = \left[{\matrix{{\cos \theta_G} & {\sin \theta_G} & 0 \cr {\sin \theta_G} & {\cos \theta_G} & 0 \cr 0 & 0 & 1 \cr } } \right]\left[{\matrix{{\dot X + \omega_EY} \cr {\dot Y + \omega_EX} \cr {\dot X} \cr}}\right]_{ECI}$$

where θ _G is the Greenwich sidereal time, and ω _E is the Earth spin rate (7.2921158553x10⁻⁵ rad/sec). Furthermore, in order to reckon the space debris' ground track, which involves the geographic latitude and longitude, the following relations are utilised:

(48-1) $$l = \hbox{asin}\ (Z/r)$$

(48-2) $$\lambda = \hbox{if}\ m>0 ; \hbox{acos}\left( {l\over cos\lpar l \rpar } \right)\!, else; 360-\hbox{acos}\left( {l\over cos\lpar l \rpar } \right)$$

SIMULATION RESULTS

In the current work, we assume that the observer observes the debris' initial orbit from an unchanging orbit. This type of orbit is normally achieved via the Automatic Orbit Control System (AOCS) of the observer. The debris is also supposed to be the Switzerland CubeSat which encounters perturbations. Moreover, initial orbital elements for observer and space debris and orbital characteristics of the space debris obtained from the Lambert Method are presented in Table 1. By comparison with the reference data, it is clear that the employed hybrid procedure in this paper has sufficient accuracy during the preliminary orbit determination.

Table 1. Orbital elements of observer, space debris, and determined orbit.

The simulated orbits of the space debris and the observer robot are shown in Figure 5. In the next step, the variations of the space debris' position vector obtained using EKF are considered during the specified period. For observing the position vector variations, three different figures are considered based on the variations along x, y, and z. A view of simulated relative trajectory of space debris and the observer are shown in Figure 6, which demonstrates the complexity of the scenario. Moreover, the history of the reference orbit together with the estimated orbit by EKF are shown in Figure 6(b), where one might obtain a visual insight on the intimacy between them of the results.

Figure 5. A view of simulated orbits of space debris and the observer.

Figure 6. (a) A 3D/2D view of simulated relative trajectory of space debris and the observer. (b) History of the reference orbit and the estimated orbit by EKF.

The variations of the space debris' position vector along the x-axis are illustrated in Figure 7(a). As exposed, these variations due to different harmonies in the orbits of space debris and the observer are sinusoidal. Moreover, the distance between two peaks indicates the complete orbit cycle. The variations of the elements of space debris' position vector, along y-and-z-axes, are displayed in Figures 7(b), 7(c), respectively. Also, the variations of the space debris' velocity components, along x, y and z-axes, are exhibited in Figures 8(a), 8(b), and 8(c), in turn. The error components of the space debris' position vector are portrayed in Figures 9(a), 9(b), and 9(c), respectively. The error components of the space debris' velocity vector, along x, y and z-axes, are demonstrated in Figure 10(a), 10(b), and 10(c). As per the outcomes, error rates in the space debris' position and velocity vectors after the specified time increase a shallow slope because of the application of EKF. In other words, this issue demonstrates the efficacy of EKF, which is applied in this research. It is stated that these slight errors are created due to the uncertainties as well as produced errors because of the high simulation time. Note that after a specific time, errors will increase, and this issue indicates that it is necessary to utilise other auxiliary methods or integration of them for a precise orbit determination.

Figure 7. (a) Variations of the space debris' position vector along x-axis, (b) y-axis, and (c) z-axis.

Figure 8. (a) Variations of the space debris' velocity vector along x-axis, (b) y-axis, and (c) z-axis.

Figure 9. (a) Error variations of space debris' position vector along x-axis, (b) y-axis, and (c) z-axis.

Figure 10. (a) Error variations of space debris' velocity vector along x-axis, (b) y-axis, and (c) z-axis.

The determined orbital elements, including angular momentum, eccentricity, right ascension and inclination, true anomaly and argument of perigee are portrayed in Figures 11(a) to 11(b), respectively. The results of model-based prediction of space debris' orbit are demonstrated in Figures 12(a), 12(b), and 12(c). At the end, for better understanding of the orbital characteristics and the effects of perturbations, the ground track of the space debris and the observer are depicted in Figure 13. It is indicated that, due to considering a 98^° for inclination angle and high sensitivity of the perturbations on it, the simulation scenario experiences the uppermost possible anomalies in the orbit of the space debris.

Figure 11. Estimated (a) Angular momentum, (b) Eccentricity (c) Right ascension, (d) Inclination, (e) True anomaly, and (f) Argument of perigee of the space debris' orbit.

Figure 12. Predicted time history of the space debris orbit, (a) 3D position (b) Right Ascension, (c) Angular Momentum, (d) Argument of Perigee, and (e) Eccentricity.

Figure 13. Ground track of space debris.

CONCLUSION

The main objective of this paper has been preliminary positioning, orbit determination, and model-based trajectory prediction of space debris, based on a robust procedure using iterative bearing-only in-orbit observations. In the first place, the problem was formulated with the purpose of calculating the distance vector between the space debris and an observer satellite. In addition, the assumption employed in this work was that the orbit of the observer was specified right from the beginning. Afterwards, the obtained position vector of the debris, by Gauss method, was corrected through EKF to diminish the errors. The Lambert Method was applied in this work for calculating the velocity vectors, and eventually, orbital elements of the space debris were obtained via acquired position and velocity vectors of the space debris. For future practical applications, a model-based method was exploited for prediction of the position and orbital characteristics of the debris. The efficacy and robustness of the presented method were demonstrated by running desired simulations over a practical mission. Particularly, it was verified that the suggested scheme was capable of producing a suitable solution for preliminary localisation and orbit determination of space debris based on the presented space-based observations, which was practical in preliminary rendezvous manoeuvres such as phasing and chasing exercises of any debris grabber space-machine.