Tracking an Untracked Space Debris After an Inelastic Collision Using
  Physics Informed Neural Network
With the sustained rise in satellite deployment in Low Earth Orbits, the
collision risk from untracked space debris is also increasing. Often
small-sized space debris (below 10 cm) are hard to track using the existing
state-of-the-art methods. However, knowing such space debris' trajectory is
crucial to avoid future collisions. We present a Physics Informed Neural
Network (PINN) - based approach for estimation of the trajectory of space
debris after a collision event between active satellite and space debris. In
this work, we have simulated 8565 inelastic collision events between active
satellites and space debris. Using the velocities of the colliding objects
before the collision, we calculate the post-collision velocities and record the
observations. The state (position and velocity), coefficient of restitution,
and mass estimation of un-tracked space debris after an inelastic collision
event along with the tracked active satellite can be posed as an optimization
problem by observing the deviation of the active satellite from the trajectory.
We have applied the classical optimization method, the Lagrange multiplier
approach, for solving the above optimization problem and observed that its
state estimation is not satisfactory as the system is under-determined.
Subsequently, we have designed Deep Neural network-based methods and Physics
Informed Neural Network (PINN )based methods for solving the above optimization
problem. We have compared the performance of the models using root mean square
error (RMSE) and interquartile range of the predictions. It has been observed
that the PINN-based methods provide a better prediction for position, velocity,
mass and coefficient of restitution of the space debris compared to other
methods.
Harsha M., Gurpreet Singh, Vinod Kumar, Arun Balaji Buduru, Sanat K. Biswas
Tracking an Untracked Space Debris After an Inelastic Collision Using Physics Informed Neural Network

###### Abstract

With the sustained rise in satellite deployment in Low Earth Orbits, the collision risk from untracked space debris is also increasing. Often small-sized space debris (below 10 cm) are hard to track using the existing state-of-the-art methods. However, knowing such space debris' trajectory is crucial to avoid future collisions. We present a Physics Informed Neural Network (PINN) - based approach for estimation of the trajectory of space debris after a collision event between active satellite and space debris. In this work, we have simulated 8565 inelastic collision events between active satellites and space debris. To obtain the states of the active satellite, we use the TLE data of 1647 Starlink and 66 LEMUR satellites obtained from space-track.org. The velocity of space debris is initialized using our proposed velocity sampling method, and the coefficient of restitution is sampled from our proposed Gaussian mixture-based probability density function. Using the velocities of the colliding objects before the collision, we calculate the post-collision velocities and record the observations. The state (position and velocity), coefficient of restitution, and mass estimation of un-tracked space debris after an inelastic collision event along with the tracked active satellite can be posed as an optimization problem by observing the deviation of the active satellite from the trajectory. We have applied the classical optimization method, the Lagrange multiplier approach, for solving the above optimization problem and observed that its state estimation is not satisfactory as the system is under-determined. Subsequently, we have designed Deep Neural network-based methods and Physics Informed Neural Network (PINN )based methods for solving the above optimization problem. We have compared the performance of the models using root mean square error (RMSE) and interquartile range of the predictions. It has been observed that the PINN-based methods provide a better prediction for position, velocity, mass and coefficient of restitution of the space debris compared to other methods.

## Introduction

Since the launch of Sputnik in 1957, approximately 13,000 satellites have been deployed in the Low Earth orbit [1]. With collisions, explosions and fragmentation, the number of space debris sized less than 1 cm is now estimated to be more than 1,000,000 [1]. The US Space Surveillance Network has a catalogue of about 30,000 resident space objects only [1]. Since 97% of space debris is not tracked, they pose a greater threat to active satellites as well as future space missions.

With this large number of space debris, the risk of in-orbit collision has increased, and the number of collision events has seen a rise in the past decade. Collision in space can be destructive, for example, Iridium-Cosmos collision [2] of 2009, or non-destructive, for example, Canadarm2 of the International Space Station hit by a small piece of space debris in 2021 [3]. Another incident which was catalogued as a non-destructive collision was in 2013 when Blitz satellite was possibly hit by un-tracked space debris [4]. While there has been exhaustive research on modelling destructive collisions [5, 6], on risk assessment due to small space debris [7], and reconnecting fragments in space to their parent satellite body [8], the extraction of trajectory information of untracked space debris observing a non-destructive collision has not been studied extensively.

Recently, Harsha et al. formulated a space debris position, velocity and mass estimation problem considering the position and velocity deviation of an active satellite as observation with non-destructive and elastic collision assumption [9]. The performance of the classical estimation methods as well as machine learning (ML) based approaches, were examined for this problem. It was observed that the position, velocity and mass estimation performance of the Ensemble Neural Network-based technique are similar to the classical methods. It should be noted that the elastic collision assumption in the above preliminary study was idealistic. The above estimation problem becomes complex when a more rigorous inelastic collision model is considered. Hence it is of interest to examine the trajectory estimation performance of both classical and ML-based methods under the inelastic collision assumption.

Deep Neural Network (DNN)-based techniques have already been proposed for asteroid exploration, spacecraft rendezvous and terrain navigation [10]. The convolutional neural network (CNN) has been proposed for pose estimation of uncooperative
spacecraft[11] and for object detection and tracking in space[12]. NASA's Double asteroid redirection test (DART) mission successfully used ML techniques to navigate the spacecraft autonomously and collide with Didymos by real-time analysis and processing of data from the spacecraft's onboard cameras and range-finding instruments. Further, the ML algorithm optimized trajectories and evaluated various impact scenarios, considering factors such as the asteroid's shape, composition, and spin[13].

Supervised ML methods have found applications across many fields of engineering and science, for example object detection and recognition[14, 15], recommendation systems[16], text-to-speech conversion[17], machine translation[18], and information retrieval[19]. These methods achieve satisfactory performance when trained with large amounts of data. However, these standard neural network-based methods inevitably face the challenge of drawing conclusions and making decisions under partial information, where the data acquisition is costly and scarce while dealing with complex and non-linear physical systems[20]. In addition, the solution of the traditional neural network does not guarantee adherence to the underlying physical properties in problems involving physical systems. To address these issues, the Physics informed Neural network (PINN) was proposed[20]. PINN[21, 22] uses the laws governing the system dynamics in the loss function as prior information. Note that the training of a neural network essentially implies finding the set of weights and biases for all the nodes in the network, which minimises the loss function. Hence, the inclusion of the system dynamics based on physical laws in the loss function act as a constraint for the output while training the neural network using the training data set. The PINNs can be used as surrogates in learning models used for autonomous navigation, uncertainty quantification and any real-time application that need inference[23].

In this article, we have studied the problem of tracking unknown space debris, observing an inelastic and non-destructive collision event. Additionally, we assumed the space debris was in a stable orbit before the collision. Under the inelastic collision assumption, the coefficient of restitution[24] is an unknown parameter in addition to the variables that need to be estimated under the elastic collision assumption. We have examined the performance of the classical estimation technique and various DNN and PINN-based methods for estimating the position, velocity, and mass of the space debris and the coefficient of restitution for the inelastic collision for the above-mentioned problem using 8235 collision simulations for training and 330 collision simulations for testing. The inelastic collisions were simulated using the model described by Schwager et al.[25]. The performance of the DNN and PINN-based methods are compared with the classical method for the inelastic and non-destructive collision. Towards tracking unknown space debris, we have made the following key contributions:

1. Formulation of space debris tracking problem observing a non-destructive and inelastic collision: We have posed the unknown space debris tracking problem after an inelastic event as a position, velocity, mass, and the coefficient of restitution estimation problem considering the active satellite position and velocity deviation as observations.
2. Formulation of an appropriate physics loss function and application of PINN: We have formulated a physics loss function for the above problem and trained a PINN using the developed loss function to solve the above estimation problem.
3. A velocity sampling method for simulating random in-orbit collision events: For the training of DNN and PINN-based models as well for testing the performance of the ML-based approaches and the classical approach, we simulated a total of 8565 inelastic collision events. Note that, for a collision event, the active satellite position and the space debris position will be approximately the same at the time of the collision, while the velocities will be different. However, any random velocity vector at the given position does not necessarily result in a stable elliptical orbit. We have proposed a velocity sampling method to generate the velocity vector of space debris for a given active satellite position at the time of the collision, which guarantees an elliptical orbit with a periapsis above a given minimum allowable altitude from the mean sea level.

The rest of the article is organised as follows: We have formally defined the problem in the next section, and then discussed the classical estimation method and various DNN-based approaches and provided the physics loss function for the PINN for solving the problem. Then we have discussed the collision data generation method, including the formulation of the velocity sampling algorithm and simulation of collision in space. Next, we have compared the performance of the classical estimation technique as well as DNN and PINN-based methods for solving the estimation problem. We conclude the article by consolidating our observations and delineating the future research direction.

## Problem Definition

Consider a satellite and space debris undergoing an inelastic collision in space. Our objective is to find the trajectory of the space debris by observing the position and velocity deviation of the active satellite after the collision. Let \(\mathbf{r}_{sat}\) and \(\mathbf{v}_{sat}\) be the position and velocity of the satellite, and \(\mathbf{r}_{d}\) and \(\mathbf{v}_{d}\) be the position and velocity of the space debris. We assume satellite and space debris are spherical for ease of the analysis. Let \(t_{c}^{-}\) and \(t_{c}^{+}\) denote the time before collision and time after collision respectively. At the time of the collision, the position of the satellite and the position of debris can be written as

\[\mathbf{r}_{sat}(t_{c}^{+})=\mathbf{r}_{sat}(t_{c}^{-}) \tag{1}\]
\[\mathbf{r}_{d}(t_{c}^{+})=\mathbf{r}_{d}(t_{c}^{-})=\mathbf{r}_{sat}(t_{c}^{-})-(\rho_{sat}+ \rho_{d})\mathbf{\hat{v}}_{rel} \tag{2}\]

where \(\rho_{sat}\) and \(\rho_{d}\) are radii of the satellite and debris, respectively and \(\mathbf{\hat{v}}_{rel}\) is the unit vector along the relative velocity of debris with respect to the satellite and is given by

\[\mathbf{\hat{v}}_{rel}=\frac{\mathbf{v}_{sat}(t_{c}^{-})-\mathbf{v}_{d}(t_{c}^{-})}{||\mathbf{ v}_{sat}(t_{c}^{-})-\mathbf{v}_{d}(t_{c}^{-})||} \tag{3}\]

The velocity of the satellite and debris after an inelastic collision is[24]

\[\mathbf{v}_{sat}(t_{c}^{+})=\frac{m_{sat}\mathbf{v}_{sat}(t_{c}^{-})+m_{d}\mathbf{v}_{d} (t_{c}^{-})+\epsilon m_{d}(\mathbf{v}_{d}(t_{c}^{-})-\mathbf{v}_{sat}(t_{c}^{-}))}{m_{ sat}+m_{d}} \tag{4}\]

\[\mathbf{v}_{d}(t_{c}^{+})=\frac{m_{d}\mathbf{v}_{d}(t_{c}^{-})+m_{sat}\mathbf{v}_{sat}(t_ {c}^{-})+\epsilon m_{sat}(\mathbf{v}_{sat}(t_{c}^{-})-\mathbf{v}_{d}(t_{c}^{-}))}{m_{ sat}+m_{d}} \tag{5}\]

where \(\mathbf{v}_{sat}(t_{c}^{-})\)indicates velocity of satellite, just before the collision and \(\mathbf{v}_{sat}(t_{c}^{+})\) indicates velocity of satellite just after the collision and \(\epsilon\) is the Coefficient of Restitution[25] of the inelastic collision. Then change in velocity of the active satellite can be written as

\[\Delta\mathbf{v}_{sat}(t_{c}^{+})=\mathbf{v}_{sat}(t_{c}^{+})-\mathbf{v}_{sat}(t_{c}^{-}) \tag{6}\]

Define the measurement \(\mathbf{Z}\) as

\[\mathbf{Z}=\begin{bmatrix}\mathbf{r}_{sat}(t_{c}^{+})\\ \Delta\mathbf{v}_{sat}(t_{c}^{+})\end{bmatrix}+\omega(t_{c}^{+})=\begin{bmatrix} \mathbf{r}_{sat}(t_{c}^{+})\\ \mathbf{v}_{sat}(t_{c}^{+})-\mathbf{v}_{sat}(t_{c}^{-})\end{bmatrix}+\omega(t_{c}^{+}) \tag{7}\]

where \(\omega(t_{c}^{+})\) is noise in measurement. As the active satellite is being tracked, \(\mathbf{v}_{sat}(t_{c}^{+})\) and \(\mathbf{v}_{sat}(t_{c}^{-})\) can be measured, and therefore measurement \(\mathbf{Z}\) is available. Using equations (2), (4) and (5), \(\mathbf{Z}\) can be written as

\[\mathbf{Z}=f(\mathbf{r}_{d}(t_{c}^{-}),\mathbf{v}_{d}(t_{c}^{-}),m_{d},\epsilon)+\omega(t_ {c}^{+}) \tag{8}\]

Obtaining \(\mathbf{r}_{d}(t_{c}^{-}),\mathbf{v}_{d}(t_{c}^{-}),m_{d}\) and \(\epsilon\) from \(\mathbf{Z}\) is essentially an estimation problem. Using the estimated parameters \(\mathbf{v}_{d}(t_{c}^{+})\) can be obtained and thus the trajectory of the space debris after the collision can be computed.

Figure 1: Collision
### Classical approach

One can estimate \(\mathbf{r}_{d}(t_{c}^{-}),\mathbf{v}_{d}(t_{c}^{-}),m_{d}\) and \(\epsilon\) from observation \(\mathbf{Z}\) using the Least Squared Estimation (LSE):

\[\widehat{\mathbf{r}}_{d}(t_{c}^{-}),\widehat{\mathbf{v}}_{d}(t_{c}^{-}),\widehat{m}_{d}, \widehat{\epsilon}=\operatorname*{arg\,min}_{\mathbf{r}_{d}(t_{c}^{-}),\mathbf{v}_{d} (t_{c}^{-}),m_{d},\mathbf{\epsilon}}(\mathbf{\Lambda}\mathbf{Z}^{T}\mathbf{\Lambda}\mathbf{Z}) \tag{9}\]

where \(\widehat{\mathbf{r}}_{d}(t_{c}^{-}),\widehat{\mathbf{v}}_{d}(t_{c}^{-}),\widehat{m}_{d}, \widehat{\epsilon}\) are the estimate of the desired quantities and

\[\mathbf{\Lambda}\mathbf{Z}=\mathbf{Z}-f(\widehat{\mathbf{r}}_{d}(t_{c}^{-}),\widehat{\mathbf{v}}_{ d}(t_{c}^{-}),\widehat{m}_{d},\widehat{\epsilon}) \tag{10}\]

Note that the above estimation problem formulation does not consider any mass constraints for the debris. However, we are particularly interested in the untracked debris, which is of small size and essentially has a very small mass. The inclusion of the mass constraint transforms the original problem into a constrained optimization problem, which can be solved using the Lagrange Multiplier approach [9]. The Lagrangian \(\mathcal{L}\) for this constrained optimization problem can be defined as

\[\mathcal{L}=\mathbf{\Lambda}\mathbf{Z}^{T}\mathbf{\Lambda}\mathbf{Z}+\lambda_{1}(1-m_{d})+ \lambda_{2}(m_{d}-\delta m) \tag{11}\]

which includes the mass constraint \(\delta m<m_{d}<1\ kg\). Here, \(\lambda_{1}\) and \(\lambda_{2}\) are Lagrange multipliers. The position, velocity, mass and coefficient of restitution of the debris can be estimated by solving

\[\begin{bmatrix}\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{d}}& \frac{\partial\mathcal{L}}{\partial\mathbf{r}_{d}}&\frac{\partial\mathcal{L}}{ \partial m_{d}}&\frac{\partial\mathcal{L}}{\partial\epsilon}\\ &&\frac{\partial\mathcal{L}}{\partial\lambda_{1}}&\frac{\partial\mathcal{L}}{ \partial\lambda_{2}}\\ \end{bmatrix} =\mathbf{0} \tag{12}\] \[\begin{bmatrix}\frac{\partial\mathcal{L}}{\partial\lambda_{1}}& \frac{\partial\mathcal{L}}{\partial\lambda_{2}}\\ \end{bmatrix} =\mathbf{0} \tag{13}\]

There are various techniques available to solve this minimisation problem. We have used the gradient descent approach to solve the problem numerically.

### Machine Learning-based approaches

An alternate approach for estimating the position, velocity, mass and coefficient of restitution of the debris after a collision can be based on a DNN model. Using (8) one can write

\[E[\mathbf{r}_{d}(t_{c}^{-}),\mathbf{v}_{d}(t_{c}^{-}),m_{d},\mathbf{\epsilon}]=f^{-1}(\mathbf{Z}) \tag{14}\]

considering \(\mathbf{\omega}(t_{c}^{+})\) as a zero mean noise vector. Essentially the trained DNN model approximates the inverse function \(f^{-1}(\cdot)\). The DNN model can be trained using simulated collision data. We consider \(\mathbf{Z}\) as the input to the DNN and \(\mathbf{\Lambda}\mathbf{r},\mathbf{\Delta}\mathbf{v},m_{d},\mathbf{\epsilon}\) as the output of the DNN, where

\[\mathbf{\Delta}\mathbf{r}=\begin{bmatrix}\Delta x_{d}\\ \Delta y_{d}\\ \Delta z_{d}\end{bmatrix}=\mathbf{r}_{d}(t_{c}^{-})-\mathbf{r}_{sat}(t_{c}^{-}) \tag{15}\]

and

\[\mathbf{\Delta}\mathbf{v}=\begin{bmatrix}\Delta v_{xd}\\ \Delta v_{yd}\\ \Delta v_{zd}\end{bmatrix}=\mathbf{v}_{d}(t_{c}^{-})-\mathbf{v}_{sat}(t_{c}^{-}) \tag{16}\]

We consider two distinct DNN architectures for this problem. The first architecture is the usual DNN with an input layer with 6 inputs, multiple hidden layers and an output layer with 8 outputs as shown in Fig. 2. The second architecture comprises of 8 parallel DNNs with 6 inputs, multiple hidden layers and 1 output. Essentially, each of these parallel DNNs learns one of the 8 outputs. The outputs of these 8 parallel DNNs are connected to another DNN with 8 inputs and 8 outputs. The last DNN block ensures learning the dependency of each output variable with other output variables. The architecture is shown in Fig. 3.

We refer to this architecture as the stacked Deep Neural Network (StackDNN). For both architectures, Mean Squared Error (MSE) loss function is used for training the models. We define the MSE loss for this problem as

\[\begin{split}&\mathcal{L}_{MSE}(\widehat{\mathbf{r}}_{d},\widehat{ \mathbf{v}}_{d},\hat{m}_{d},\mathbf{\epsilon})\\ &=\frac{1}{N}\sum_{i=1}^{N}\left[(\mathbf{r}_{d_{i}}-\widehat{\mathbf{r}}_{d_ {i}})^{T}(\mathbf{r}_{d_{i}}-\widehat{\mathbf{r}}_{d_{i}})+(\mathbf{v}_{d_{i}}-\widehat{ \mathbf{v}}_{d_{i}})^{T}(\mathbf{v}_{d_{i}}-\widehat{\mathbf{v}}_{d_{i}})+(m_{d}-\hat{m}_{ d})^{2}+(\mathbf{\epsilon}-\mathbf{\epsilon})^{2}\right]\end{split} \tag{17}\]

where \((\cdot)_{d_{i}}\) and \((\widehat{\cdot})_{d_{i}}\) denote the output data for the \(i^{th}\) collision event and the corresponding predicted output data using the DNN. \(N\) is the number of collision events used to generate the training data.
### Physics Informed Neural Network

Having discussed the DNN-based method for estimating the position, velocity, mass and coefficient of restitution of debris after the collision, it is essential to note that the estimated position and velocity must be the solution to the differential equation that governs the motion of the debris. Based on Newton's Law of Gravitation, the debris dynamics can be written as

\[\begin{bmatrix}\dot{\mathbf{r}}_{d}\\ \dot{\mathbf{v}}_{d}\end{bmatrix}=\begin{bmatrix}\mathbf{v}_{d}\\ -\frac{\mu}{r_{d}^{3}}\mathbf{r}_{d}\end{bmatrix} \tag{18}\]

However, DNN-based solutions do not guarantee the above constraint. The Physics Informed Neural Network [21] preserves the laws of physics in the solution by incorporating a physics loss in the loss function used to train the DNN. Using (18), we defined the physics loss for the PINN as

\[\begin{split}&\mathcal{L}_{phy}(\mathbf{\widehat{r}}_{d},\mathbf{ \widehat{r}}_{d})\\ &=\frac{1}{N}\sum_{i=1}^{N}\left[\left(\frac{\mu}{r_{d_{i}}^{3}} \mathbf{r}_{d_{i}}-\frac{\mu}{\hat{r}_{d_{i}}^{3}}\mathbf{\widehat{r}}_{d_{i}}\right) ^{T}\left(\frac{\mu}{r_{d_{i}}^{3}}\mathbf{r}_{d_{i}}-\frac{\mu}{\hat{r}_{d_{i}}^{ 3}}\mathbf{\widehat{r}}_{d_{i}}\right)+\left(\mathbf{v}_{d_{i}}-\mathbf{\widehat{r}}_{d_{i }}\right)^{T}\left(\mathbf{v}_{d_{i}}-\mathbf{\widehat{r}}_{d_{i}}\right)\right]\end{split} \tag{19}\]

for the problem under consideration. The significance of this loss function is that the minimisation of \(\mathcal{L}_{phy}(\mathbf{\widehat{r}}_{d},\mathbf{\widehat{r}}_{d})\) minimises the difference between the acceleration computed using the output of the training data and the PINN output, and the difference between the velocity computed using the output of the training data and the PINN output. As a result, this loss function minimises the distance between the derivative of the vector \(\left[\mathbf{r}_{d}\ \ \mathbf{v}_{d}\right]^{T}\) computed using the training data, and the derivative of the same computed using the PINN predicted output. Since the velocity loss is already included in \(\mathcal{L}_{phy}(\mathbf{\widehat{r}}_{d},\mathbf{\widehat{r}}_{d})\) we define the MSE loss for the PINN as

\[\mathcal{L}_{MSEP}(\mathbf{\widehat{r}}_{d},\mathbf{\widehat{m}}_{d},\mathbf{\hat{\epsilon }})=\frac{1}{N}\sum_{i=1}^{N}\left[(\mathbf{r}_{d_{i}}-\mathbf{\widehat{r}}_{d_{i}})^ {T}(\mathbf{r}_{d_{i}}-\mathbf{\widehat{r}}_{d_{i}})+(m_{d}-\mathbf{\widehat{m}}_{d})^{2} +(\epsilon-\mathbf{\hat{\epsilon}})^{2}\right] \tag{20}\]

Figure 3: Stacked Deep Neural Network (StackDNN) architecture

Figure 2: Deep Neural Network (DNN) architecture
The complete loss function for the PINN is

\[\mathcal{L}_{PINN}=\mathcal{L}_{phy}\left(\widehat{\mathbf{r}}_{d},\widehat{\mathbf{r}}_{ d}\right)+\mathcal{L}_{MSEP}\left(\widehat{\mathbf{r}}_{d},\hat{m}_{d},\hat{\mathbf{e}}\right) \tag{21}\]

We use this loss function in both DNN and StackDNN for performance comparison.

## Collision data generation

Having discussed the DNN, StackDNN, PINN, and StackPINN approaches for tracking the debris after an inelastic collision event, let us now turn to the preparation of the training data. Publicly available data on inelastic and nondestructive collision events are very few and not sufficient to train the neural networks. Hence simulated collision data are required for training the machine learning models.

The position and velocity information of active satellites can be obtained from publicly available catalogues. However, the position and velocity of the debris which collides with an active satellite at a given epoch must be simulated. The simulated debris must be in a complete orbit, i.e. the orbit eccentricity must be less than 1. Additionally, to form a stable orbit, the periapsis of the debris cannot be at a very low altitude.

The debris position is given by (2) at the time of the collision. From (2), we can approximate that the debris position as nearly equal to the active satellite position because \(r_{sat}\gg\rho_{sat},\rho_{d}\). Corresponding to this position, infinitely many velocity vectors are possible for which the debris will be on a complete orbit. However, the magnitude and the direction of the velocity vector can not be any arbitrary magnitude and direction. In the next subsection, we describe a method of sampling the velocity of debris for a given position vector in space and a given minimum periapsis altitude.

### Velocity sampling

Consider \(r\) is the norm of the active satellite as well as the debris position vector \(\mathbf{r}\) at the time of the collision, \(h\) is the norm of the angular momentum vector \(\mathbf{h}\) of the colliding debris, \(e\) is the eccentricity of the debris orbit, \(\phi\) is the true anomaly of the debris during the collision, \(R_{e}\) is the Radius of the earth and \(a_{min}\) minimum allowable periapsis altitude. If the semi-major axis for the debris orbit is \(a\), then the periapsis \(r_{p}=a(1-e)\). Also

\[\frac{h^{2}}{\mu}=r_{p}(1+e) \tag{22}\]

and

\[r=\frac{h^{2}}{\mu}\frac{1}{1+e\cos\phi} \tag{23}\]

Considering \(R_{min}=R_{e}+a_{min}\) as the minimum allowable distance from the earth centre

\[r_{p}=\frac{r(1+e\cos\phi)}{1+e}>R_{min} \tag{24}\]

then

\[0<e<\frac{r-R_{min}}{R_{min}-r\cos\phi}<1 \tag{25}\]

and

\[\cos^{-1}\frac{2R_{min}-r}{r}<\phi<2\pi-\cos^{-1}\frac{2R_{min}-r}{r} \tag{26}\]

Fig. 4 shows the variation of maximum possible eccentricity with the true anomaly. Now, the position of the debris at the time of collision is

\[\mathbf{r}_{d}=\mathbf{r}=r\begin{bmatrix}r_{ux}\\ r_{uy}\\ r_{uz}\end{bmatrix} \tag{27}\]

and velocity of the debris is

\[\mathbf{v}_{d}=v\begin{bmatrix}v_{ux}\\ v_{uy}\\ v_{uz}\end{bmatrix} \tag{28}\]

where \(\begin{bmatrix}r_{ux}&r_{uy}&r_{uz}\end{bmatrix}^{T}\) and \(\begin{bmatrix}v_{ux}&v_{uy}&v_{uz}\end{bmatrix}^{T}\) are the position and velocity unit vectors respectively, and

\[v=\sqrt{\frac{\mu}{r}\left(2-\frac{1-e^{2}}{1+e\cos\phi}\right)} \tag{29}\]
The angle between \(\mathbf{r}\) and \(\mathbf{v}\) is \(\theta^{26}\)

\[\theta=\sin^{-1}\frac{1+e\cos\phi}{\sqrt{1+2e\cos\phi+e^{2}}} \tag{30}\]

The angular momentum vector of the debris is

\[\mathbf{h} =\mathbf{r}\times\mathbf{v}\] \[=rv\begin{bmatrix}r_{\alpha y}v_{uz}-r_{uz}v_{uy}\\ r_{uz}v_{ux}-r_{ux}v_{uz}\\ r_{ux}v_{uy}-r_{uy}v_{ux}\end{bmatrix} \tag{31}\]

and

\[h=rv\sin\theta \tag{32}\]

then the cosine of the inclination \(i\) of the orbit

\[\cos i=\frac{rv(r_{ux}v_{uy}-r_{uy}v_{ux})}{rv\sin\theta} \tag{33}\]

Consequently, the velocity unit vector needs to satisfy the following:

\[r_{ux}v_{ux}+r_{uy}v_{uy}+r_{uz}v_{uz}=\cos\theta \tag{34}\] \[-r_{uy}v_{ux}+r_{ux}v_{uy}+0v_{uz}=\sin\theta\cos i\] (35) \[\sqrt{v_{ux}^{2}+v_{uy}^{2}+v_{uz}^{2}}=1 \tag{36}\]

Geometrically, the condition for forming a complete orbit is - the line of intersection of planes (34) and (35) must lie within the unit sphere defined by (36). We can write (35) in the normal form:

\[-\frac{r_{uy}}{\sqrt{r_{ux}^{2}+r_{uy}^{2}}}v_{ux}+\frac{r_{ux}}{\sqrt{r_{ux }^{2}+r_{uy}^{2}}}v_{uy}+0.v_{uz}=\frac{\sin\theta\cos i}{\sqrt{r_{ux}^{2}+r_ {uy}^{2}}} \tag{37}\]

Since \(\sqrt{r_{ux}^{2}+r_{uy}^{2}+r_{uz}^{2}}=1\) (34) is already in the normal form. Observing (34) and (37), one can deduce that the angle between the two planes is \(90^{\circ}\). Then from the cross-sectional geometry shown in Fig. 5, the criteria for the intersection of the two planes within the unit sphere is

\[\cos^{2}\theta+\frac{\sin^{2}\theta\cos^{2}i}{r_{ux}^{2}+r_{uy}^{2}}\leq 1 \tag{38}\]
Then

\[-\sqrt{\frac{r_{ux}^{2}+r_{uy}^{2}}{\sin^{2}\theta(1-\cos^{2}\theta)}}\leq\cos i \leq\sqrt{\frac{r_{ux}^{2}+r_{uy}^{2}}{\sin^{2}\theta(1-\cos^{2}\theta)}} \tag{39}\]

Using the ranges of \(\phi\), \(e\) and \(\cos i\) given by (26), (25) and (39) one can independently sample these quantities and subsequently solve for \(v_{ux}\), \(v_{uy}\) and \(v_{uz}\) from (34), (35) and (36). Note that the solution will result in two unit vectors for velocity. Finally, the velocity magnitude can be calculated using (29). Algorithm 1 describes the velocity sampling steps.

```
Input :\(\mathbf{r}\), minimum altitude, Number of samples Output :Velocity samples
1fori=1:Number of samplesdo
2 Given the allowed minimum altitude, find the minimum and maximum true anomaly \(\phi_{min}\) and \(\phi_{max}\) for the given position vector using (26);
3 Sample \(\phi\sim U(\phi_{min},\phi_{max})\);
4 For a sampled \(\phi\) calculate maximum eccentricity \(e_{max}\) using (25);
5 Sample \(e\sim U(0,e_{max})\);
6 Calculate \(\theta\);
7 Compute the lower bound \(\cos i_{min}\) and upper bound \(\cos i_{max}\) for \(\cos i\) using (39);
8 Sample \(\cos i\sim U(\cos i_{min},\cos i_{max})\);
9 Solve for \(v_{ux},v_{uy},v_{uz}\) using (34), (35) and (36);
10 Calculate \(v\) using (29)
11 end for
```

**Algorithm 1**Velocity sampling

Using this algorithm, we have sampled 2000 velocity vectors for an active satellite position at

\[\mathbf{r}=\left[3664.75536654\quad 3893.60049258\quad 5062.87536598\right]^{T}km \tag{40}\]

and a minimum perigee altitude of 300 km. Fig. 6 shows that the generated samples have perigees at altitudes higher than 300 km. Fig. 7 shows the relative velocities for the samples with respect to the active satellite. It can be deduced from Fig. 7 that collision cannot occur from any arbitrary direction for a given position vector in space.

### Simulation of the collision

We have generated inelastic collision observations considering the collision model described in the Problem Formulation section. We have used Python 3.7 with astropy[27] and polastro[28] libraries for orbit simulations of the active satellites as well as the debris.

We have considered 5 collision epochs \(t_{c}\) each at 1200 UTC, starting from April 29, 2023, to May 3, 2023. We consider 1647 satellite orbits from the Starlink constellation and 66 satellite orbits from the LEMUR constellation for generating collision simulations. The Two Line Element (TLE) data for each constellation were obtained from www.spacetrack.org on April
27, 2023. Using TLE data, we have calculated the time difference \(\Delta t\) between the collision time \(t_{c}\) and the the TLE time, and propagated the active satellites. The orbits are propagated considering J2 and J3 zonal harmonics and the exponential atmospheric drag model.

For the debris simulations, we sampled the mass of each untracked debris from a normal distribution with 0.5 kg mean and 0.001 kg of standard deviation. For simplicity, satellites and debris are considered spherical objects of radius 0.5 m and 0.1 m, respectively, and the mass of the satellite is 100 times that of sampled mass of debris. The position of the debris \(\mathbf{r_{d}}\) is considered to be the same as \(\mathbf{r_{sat}}\), whereas the velocity of the debris \(\mathbf{v_{d}}\) is obtained using algorithm 1. The minimum perigee altitude for the debris corresponding to the Starlink satellites was set at 250 km, and for LEMUR, it was 200 km. The Gabbard plots for the debris generated for the Starlink and LEMUR constellations are shown in Fig. 8.

Co-efficient of Restitution (\(\epsilon\)) values for each of the collisions have been sampled from a Gaussian mixture distribution, as shown in Fig. 9. It consists of 4 Gaussian distributions with 0.015 standard deviations sampled around uniformly distributed means ranging from 0.5 to 1 where the mean has been assumed based on the \(\epsilon\) values of the satellite materials (6061-T6 Aluminium, Stainless Steel 304, Nitronic 60A and Titanium)[24].

Inelastic collision simulations are generated using the position, velocity, mass and co-efficient of restitution values for the respective satellite and debris pairs in equations (4) and (5). Post-collision at time \(t_{c}^{+}\), the observations are generated using equation (7). Here, the noise \(\omega(t_{c}^{+})\) signifies the tracking uncertainty of the active satellite. We consider 100 m uncertainty for position and 50 m/sec uncertainty for velocity in each axis.

A total of 8235 collision events were simulated for the orbits of the Starlink constellation using Param Pravega Super Computer. 330 collision events were simulated for the LEMUR constellation using a local workstation. For each collision, we

Figure 8: Gabbard plot of space debris

Figure 9: Co-efficient of Restitution (\(\epsilon\)) distribution using Gaussian mixture
have recorded \(\mathbf{Z}\) and \(\mathbf{r}_{d}\), \(\mathbf{v}_{d}\), \(m_{d}\) and \(\epsilon\). We consider the collision data corresponding to the orbits of the Starlink constellation as the training data and the same corresponding to the orbits of the LEMUR constellation as the test data.

## Results and discussion

We used the training data described in the previous section for training the DNN, PINN, StackDNN and StackPINN. We trained four different models with 2, 4, 6 and 8 hidden layers for each of the ML models, with 5 nodes in each hidden layer. We have selected Rectified Linear Unit (ReLU) as the activation function for each node. All the models are trained for 100 training epochs. The training is performed in a workstation with 8 Intel Xeon processor cores and 62GB of RAM. The training time required for each model is provided in Table 1.

The change in training loss with training epochs for all the ML models is shown in Fig. 10. Fig. 10 shows the training loss for DNN and PINN models. It is observed that for models with an increased number of hidden layers, the training epochs decrease for the training loss to converge. A similar trend can be observed for the StackDNN and StackPINN models in Fig. 10.

Using the observation \(\mathbf{Z}\) for all 330 test collision events in the trained ML models, \(\mathbf{r}_{rso}(t_{c}^{-})\), \(\mathbf{v}_{rso}(t_{c}^{-})\), \(m_{rso}\) and \(\epsilon\) were estimated. Additionally, the same variables were estimated using the Lagrange multiplier method for the same set of observations.

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline \multirow{2}{*}{**Methods**} & \multicolumn{4}{c|}{**Time Taken to train models (in seconds)**} \\ \cline{2-5}  & **2 Layers** & **4 Layers** & **6 Layers** & **8 Layers** \\ \hline
**DNN** & 51.983 & 55.048 & 58.169 & 60.993 \\
**PINN** & 652.623 & 655.689 & 663.668 & 659.781 \\
**StackDNN** & 420.531 & 447.919 & 488.353 & 512.025 \\
**StackPINN** & 1015.206 & 1049.786 & 1076.624 & 1117.294 \\ \hline \end{tabular}
\end{table}
Table 1: Time taken to train various models

Figure 10: Training loss
We recorded the error in estimation for performance comparison.

Fig. 11 shows the box plots for the position estimation errors for all the methods in the x, y and z axes in the ECI frame. It can be observed from Fig. (a)a and Fig.(b)b that the DNN-based methods are more accurate than the Lagrange multiplier method. At the same time, Fig. (a)a shows that there is a significant performance difference between DNN and PINN for 2 hidden layers. The position error in each axis using PINN has significantly less mean and median than the DNN in this case. Additionally, the interquartile range for the PINN error distribution is smaller than the DNN model. However, for 4, 6 and 8 hidden layers, the performance of the DNN and PINN models are similar and does not improve compared to the PINN with 2 hidden layers.

It can be observed from Fig. (b)b that the position estimation performance of the StackDNN and StackPINN are similar for 2, 4, 6 and 8 hidden layers and comparable with the PINN model with 2 hidden layers.

Fig. 12 shows the box plots of the error in velocity estimation along the x, y, and z directions for all the methods. It can be observed that there is no significant difference in the velocity estimation performance among the DNN, PINN, StackDNN, StackPINN and Lagrange multiplier methods. In addition, it should be noted that with the increase in hidden layers, the performance does not vary significantly.

Now, turning to the mass and the coefficient of restitution estimation, it can be observed from Fig. 13 that for 4 hidden layers and above, the performance difference of the DNN, PINN, StackDNN, and StackPINN are not significant. However, Fig. (a)a shows that with 2 hidden layers, the PINN performance is significantly better than the DNN for mass and the coefficient of restitution estimation. On the other hand, although the Lagrange multiplier method is precise in estimating the mass and the Coefficient of restitution, it introduces significant bias in the solution.

From the results shown above, it can be discerned that, for 2 hidden layers, the PINN improves the space debris position, mass and coefficient of restitution estimation performance significantly compared to the Lagrange multiplier-based and DNN-based methods. The performance of the DNN becomes similar to the PINN for higher numbers of hidden layers. However, the results show that there is no advantage to increasing the number of layers for the problem under consideration. In addition, the StackDNN and StackPINN performances are similar for the different numbers of hidden layers and are comparable to the performance of the PINN model with 2 hidden layers. However, the stacked architecture is essentially 9 different Neural Networks, increasing the complexity of implementation. From this discussion, we can conclude that the PINN with 2 hidden layers provides better performance to resource requirement trade-off for tracking untracked space debris from a non-destructive and inelastic collision with an active satellite than DNN, Stacked neural networks and classical methods.

Figure 11: Position error comparison
## Conclusion

We have formulated a problem of space debris position, velocity, mass and coefficient of restitution estimation from a non-destructive inelastic collision event under the assumption that the occurrence of the collision event is known. We have shown that this problem can be posed as a function approximation problem, and hence, a neural network can be trained to approximate the function. We proposed a cost function based on Newton's law of gravitation to design a PINN to solve this

Figure 12: Velocity error comparison

Figure 13: Mass and CoR error comparison
problem. We have also presented an algorithm for selecting the velocity of space debris for collision simulation. Using simulated collision data, we have trained DNN, PINN, StackDNN and StackPINN with varied numbers of hidden layers. We have also compared the estimation performance of these neural network-based methods with the Lagrange multiplier-based optimisation technique. The results demonstrate that the PINN with 2 hidden layers provides better estimation performance than the DNN with 2 hidden layers and the classical method. Additionally, the DNN performance becomes similar to the PINN when the number of hidden layers is increased and when Stacked neural network architecture is used. However, this increase in the model complexity does not improve the estimation performance compared to the PINN with 2 hidden layers. It is anticipated that with further research, the application of the PINN can be further extended to trajectory extraction of the debris after a break-up event. In addition, the velocity sampling formulation can be used to derive a velocity probability density function for a given position in space, which can be used for designing new filters for selecting candidate space debris for collision assessment as well as enhancing the collision probability computation.
