1 Introduction

The debris proliferation problem has been brought on by decades of unfettered space activities and the emerging tragedy of the space commons (Radtke et al. 2017; Boley and Byers 2021; Thiele and Boley 2022; Byers and Boley 2023). Robust space situational awareness (SSA) requires the rigorous and comprehensive detection, tracking, identification, characterization, and classification of resident-space objects (RSOs) and orbital events. The dominate source of space debris are fragmentation events due to on-orbit explosions, anti-satellite weapon (ASAT) tests, and collisions (Johnson et al. 2008; Byers and Boley 2023). Research is needed to better identify and characterize such breakup events and to correlate debris fragments to their parent objects (Itaya et al. 2018; Tetreault et al. 2018; Dimare et al. 2019; Wu and Rosengren 2020).

The classical representation of space-debris clouds generated from a breakup event is through a Gabbard diagram; an invention dating back to the early history of satellite-fragmentation characterization (Gabbard 1981; Reynolds et al. 2018). It plots the apogee and perigee heights of debris from a satellite’s fragmentation against their orbital periods. This tool is meant to investigate and show “the nature of the orbit, the location of the fragmentation point, the directionality and intensity of the fragments’ spread, ...(Reynolds et al. 2018)” in a concise manner. The grand opus, (Johnson et al. 2008) contains a tremendous collection of historical breakup events in this form. This classical representation, however, suffers from the time-varying nature of the orbits as both the apogee and perigee of a single object, in the perturbed circumterrestrial environment, is known to change throughout time. It is thus highly desirable to represent space-debris cloud with the latest advances in astrodynamics. Breakup-event detection and fragment correlation can be better represented by revealing their intrinsically connected dynamical nature.

Recently, Celletti et al. (2021, 2022) have proposed the use of RSO proper elements for the classification of debris into families and have worked out the theory in the circumterrestrial context using a detailed Hamiltonian formulation. Proper elements are quasi-integrals of motion (i.e., constants) resulting from the elimination of short- and long-periodic perturbations from their osculating counterparts. Objects forming different clusters in proper-elements space are associated with different breakup events (Celletti et al. 2022). The ingenuity comes from a long line of celestial-dynamics heritage. Sharing the same dynamical properties quantified by the closeness of proper elements, asteroids are also clustered into families. Supported by research in asteroid dynamics, statistical techniques, and spectroscopic composition (Zappalà et al. 1990; Milani et al. 1993; Nesvorný et al. 2015; Cellino et al. 2002), space science take the existence of asteroid families as de facto and correlates family members using asteroid proper elements; the set of analytical and empirical quasi-constant dynamical features that are qualitatively similar to RSO proper elements.

In fact, the practice of finding asteroid families is traditionally an interdisciplinary topic. The calculation of asteroid proper elements has evolved from analytical and semi-analytical approaches to the current state-of-the-art numerical techniques as computation power increased (Hirayama 1918; Beaugá and Roigb 2001; Knežević et al. 2002; Knežević 2017). The modern numerical method is made possible with both the theoretical connection and computational advances in the fast-Fourier transform (FFT) applied to celestial dynamics. Building upon on the concept of asteroid proper elements, Zappala and colleagues led the research on asteroid-family characterization through the applications of hierarchical clustering and wavelet analysis (Zappalà et al. 1990, 1994, 1995). At the intersection of statistics, data mining, and asteroid dynamics, their pioneering practice has been extended by an ever increasing catalog of asteroid proper elements (Milani et al. 2014; Nesvorný et al. 2015).

Despite the similarity in clustering proper elements to define families, asteroids and space debris have distinct differences in their origins. Asteroids are natural celestial bodies that have existed for billions of years, while space debris is the product of human activities in space, dating to the late 1950s with Sputnik. This leads to the fundamental difference that while asteroid have been present throughout human history without a single official “birth-certificate”, some space debris do possess traceable information about their provenance. The information on space debris allows novel integration with emerging deep-learning techniques, which have shared success in space engineering with smart design and reformulation, such as Martin et al. (2020) and Furfaro et al. (2019). These latest works have enabled nonlinear solutions and criteria by applying neural networks in astrodynamics problems; meanwhile in the legacy research of asteroid dynamics, most criteria are still linear such as the cutoff threshold in hierarchical clustering for asteroid-family identification.

The focus of space-debris research needs to shift from objects with known origins to those with unknown origins, as this latter kind of debris is no longer going to be scarce within the circumterrestrial domain. Current research has been focusing mainly on debris with known origins or in a controlled dynamical environment, such as in Celletti et al. (2021) and Celletti et al. (2022). However, the proliferation of space-tracking data will result in a surge in the number of space debris with unknown origins. Given current understanding on RSO proper elements, debris-family association, and the unique information on debris with known origins, many challenging questions are still waiting to be answered: Is there a way to leverage existing knowledge and practice to better deal with these uncharacterized objects? How can we innovate existing schemes with emerging technology such as neural networks? In this paper, we introduce our latest investigation and practice on space debris with unknown origins. First, we briefly review the concept of RSO proper elements in the context of Earth-orbital dynamics. Second, based on the consensus of associating breakup events with proper elements, we explore the huge data collection of space debris with known origins. Then, we introduce the neural-networks-augmented DBSCAN after applying classical DBSCAN for the first time in clustering RSO proper elements. Afterwards, we showcase our investigation on space debris with unknown origins. Lastly, we conclude our research findings.

2 RSO dynamics and proper elements

A wide variety of dynamical models are employed to approximate the diversity of trajectories around the Earth (Capderou 2014; Chao and Hoots 2018; Vallado 2022). The dynamics of RSOs (i.e., circumterrestrial satellites and debris), adopting a Newtonian formulation, can be stated as

$$\begin{aligned} \ddot{{\varvec{r}}} = -\frac{\mu }{r^3} {{\varvec{r}}} + {{\varvec{a}}}_\text {grav} + {{\varvec{a}}}_\text {non-grav}, \end{aligned}$$
(1)

where \({{\varvec{r}}}\) and \(\dot{{\varvec{r}}} = {{\varvec{v}}}\) represent the object’s state (relative position and velocity vectors), the first term is the dominate Keplerian acceleration \(\mu = G M\) (Earth’s gravitational parameter), and the disturbing accelerations \({{\varvec{a}}}_d\) are represented by the combined gravitational (grav) and non-gravitational effects. As all of the disturbing accelerations are small compared to that of the two-body acceleration term (see, e.g., Fig. 1), the perturbed two-body approach has used for over a half century for the detection and tracking of objects near Earth (Hoots et al. 2004; Setty et al. 2016; Vallado 2022). The main perturbations are caused by Shapiro (1963):

  1. 1.

    The higher harmonics of the Earth’s gravitational field;

  2. 2.

    Atmospheric drag (both neutral and charged);

  3. 3.

    The lunar and solar gravitational fields;

  4. 4.

    Direct solar radiation pressure (SRP);

  5. 5.

    Tidal effects (land, oceans); and

  6. 6.

    Solar radiation reflected from the Earth (albedo).

Of these perturbing accelerations, all but the gravitational ones depend importantly on the shape and orientation of the orbiting object. Other forces (such as the post-Newtonian relativistic correction (Huang et al. 1990), the gravitational attraction by other planets, the Lorentz force (Paul and Frueh 2017), Poynting–Robertson and solar-wind drag (Lhotka et al. 2016), the Yarkovsky effect (Murawiecka and Lemaitre 2018) are usually quite small (but not necessarily smaller than some listed) and may be incorporated as necessary for high-precision orbit prediction.

Fig. 1
figure 1

Graphical representation of several relevant forces acting on RSOs as a function of the Earth-satellite distance, shown on a log-log scale (left), and a schematic showing the various gravitational and non-gravitational interactions (right); the latter of which depend critically on the physics of the Earth-Moon-Sun dynamical environment (e.g., space weather)

Full size image

The Joint Space Operations Center (JSpOC) maintains an up-to-date SOC of all trackable Earth-orbiting objects, made publicly available in the form of two-line element (TLE) sets. There are currently \(\sim 25{,}000\) RSOs in the SOC with the majority (\(\sim 82\)%) residing in low-Earth orbit (LEO). Though widely used for decades, there remain misunderstandings of the TLE usage, its accuracy, and limitations, as well as the limitations of the simplified general perturbation theory used to produce them (SGP4) (Hoots et al. 2004). While no covariance information is provided, the estimated TLE positional accuracy at epoch is several hundred meters up to several km (1-sigma) for LEO, a few km for MEO, up to tens of km for GEO, and dozens of km for highly eccentric orbits (Flohrer et al. 2008; Früh and Schildknecht 2012). The TLEs provided Brouwer–Lyddane mean elements can be converted to osculating elements, and hence a state vector, by reconstituting the removed short-periodic terms with SGP4 (Levit and Marshall 2011).

Both the Hamiltonian formulation and Fourier-spectrum approach for proper-element computation require delving into a very long list of cumbersome technicalities. It is important enough for our purposes to understand that proper elements are simply quasi-integrals of motion obtained from the osculating orbital elements by removing all periodic oscillations produced by perturbations. Presently, method of choice to compute more “accurate”, in the sense of stable, unchanging, proper elements is to integrate the equations of celestial dynamics, and extract constants of motion directly from a numerical analysis of the predicted positions (Knežević et al. 2002; Nesvorný et al. 2015). Here, we use an ad hoc numerical scheme, akin to the state-of-the-art Fourier-based synthetic method for asteroidal proper elements, to compute proper elements for LEO debris fragments. The details of our method are described in Wu and Rosengren (2019), Wu and Rosengren (2020), Wu and Rosengren (2021) and Wu (2022) and have been submitted for publication elsewhere. We note that any proper-element computation will suffice for our purposes, so long as the obtained elements share the feature of quasi-constancy.

Our dynamical model used herein includes 50\(\times \)50 Earth harmonics, lunisolar perturbation in an ephemeris model, and SRP without shadowing effects. Osculating elements are directly calculated from the TLEs by converting the SGP4-based mean elements in the TEME (true-equator-mean-equinox) frame into osculating elements in the J2000 (mean-equator-mean-equinox) frame. Both the mean and proper elements under this more encompassing force model are then converted from their osculating counterparts in the J2000 frame. An example is illustrated for analyst object 81219, with its TLE given in Table 1.

Table 1 Two-line element set for analyst object 81219
Full size table
Table 2 Osculating, mean, and proper elements for analyst object 81219 at MJD 59013.99744441
Full size table

The TEME-frame TLE is converted to J2000 osculating elements using SGP4 and then into J2000 mean elements using a FFT algorithm; both sets are listed in Table 2. We propagate the evolution of this RSO using THALASSA (Amato et al. 2019) under a 50\(\times \)50 Earth harmonics, lunisolar perturbation, and solar radiation pressure for 1000 days. The generally assumed reflectivity coefficient of 1.2 and drag coefficient of 2.1 are used. The area of reflection, as well as the area of drag, in a cannonball assumption, are estimated to be 0.0009115 m\(^2\) with a mass of 1 kg. The resulting osculating evolution is converted into the mean counterpart for proper-elements calculation. With the specific methodology described in Wu (2022), Wu and Rosengren (2019, 2020, 2021), we extract the final proper elements and list them in Table 2.

Different from the Lie-series Hamiltonian method proposed in Celletti et al. (2021, 2022), our approach (Wu 2022; Wu and Rosengren 2019, 2021, 2020) uses a numerical scheme, which is more accessible and extendable to different models and propagators. In fact, we recommend any formulation that would adhere with the proper elements’ quasi-constant feature for comparison. Figure 2  shows the quasi-constant evolution of proper elements with respect to the changing osculating and mean counterparts, and leads to their application in finding clusters that are also quasi-invariant.

Fig. 2
figure 2

An analyst object’s osculating, mean, and proper-elements evolution. While mean and osculating elements change with respect to time, the proper semimajor axis a, proper eccentricity e and sine of proper inclination \(\sin i\) stay quasi-constant

Full size image

3 Known space-debris families represented in proper-elements space

Unlike their osculating or mean counterparts, proper elements are quasi-constant. Thus, computing RSO proper elements avoids dealing with a time-changing state representation and leads to a mathematically founded static state representation for a single object.

There are many debris with known origins according to Space-Track. Although the exact correlation method is not publicly available, we take their results as true and argue that even admitting some mis-tagged debris would not affect the general picture. We choose breakup events that have more than 100, but less than 500, associated debris fragments. The list of parent space objects and the number of valid debris are shown in Table 3. We note that there are other breakup events in the SOC with more than 500 associated debris as well as fragmentation events with less than 100. We choose this limited range mainly because of the balance between variety, how many different breakup events are included, and quantity, how many debris objects in total are considered.

Table 3 The parent space objects and the number of debris in the illustration of known space-debris families
Full size table

We calculated proper elements of around 3000 debris in Table 3. A collection of debris from the same parent space object or breakup event is defined as a debris family. The conclusion from Celletti et al. (2022) on connecting “the normal forms and proper elements in classifying space debris, especially when associated to break-up events” is linked with our definition of debris families in their dependence on proper elements.

In two separate two-dimensional plots, the distribution is more intuitively presented. Separating proper \((a,\,e,\,si)\) space into \((a,\,e)\) space and \((a,\,si)\) space is commonly practiced in asteroid-families research (Knežević and Milani 2000). Here, in Fig. 3, we plot 5 debris families: CZ-4, DMSP 5D-2 F13, NOAA 16, OPS 4682, THORAD AGENA D and the ranges of their proper elements are listed in Table 4. From the distributions, OPS 4682 is clearly out of the reach of the large chunk in both \((a,\,e)\) and \((a,\,si)\) domains. We may also see that THORAD AGENA D tends to be a little higher in a at the same e and a little lower in e at the same a with a general lower si. The separation inside the chunk composed of CZ-4, DMSP 5D-2 F13, and NOAA 16 seems less intuitive in these plots as we are assuming equal weights to \((a,\,e)\) or \((a,\,si)\) when we are looking at these plots, but not when orbital distance is introduced as a separation criteria.

Fig. 3
figure 3

Illustration of the debris families CZ-4, DMSP 5D-2 F13, NOAA 16, OPS 4682, THORAD AGENA D shown in a separate \((a,\,e)\) plot (left) and \((a,\,si)\) plot (right). Their ranges in proper semimajor axis, eccentricity, and sine of inclination are listed in Table 4

Full size image
Table 4 The range of proper elements for known debris families CZ-4, DMSP 5D-2 F13, NOAA 16, OPS 4682, THORAD AGENA D
Full size table

While there could be some intuitive explanation for the separation of the debris families, it is rather hand-wavy if there is no criteria to follow; let alone the existence of strongly overlapping distributions among families. These plots provide a glimpse of the existing clusters and distributions of debris families in the proper \((a,\,e,\,si)\) space.

Furthermore, a potential limitation with man-made RSOs is that non- gravitational affects tend to be larger than the small-body problem, and such radiation and drag perturbations can behave as de-correlators of common space events. Thus, akin to the Yarkovsky drift in asteroid proper a, the space-debris distributions could suffer slow, but continuous changes. This leads to the question: Which TLE data for a single object is more meaningful in creating a proper elements that could provide more information on defining clusters? In the creation of around 3000 debris’ proper elements, we choose to use the earliest TLE from each debris, largely in consideration of minimizing the effects of drag on proper elements, the quantification of which is still an open research question.

4 Orbital distance

Orbital distance is the concept of quantifying how close two orbiting objects would become. In order to find clusters in proper-elements space, we use the orbital-distance expression in Eq. (2),

$$\begin{aligned} d_{Z} = n_{\text {ave}} a_{\text {ave}}\sqrt{k_1(\delta a/a_{\text {ave}})^2+k_2(\delta e)^2+k_3(\delta i)^2}; \end{aligned}$$
(2)

first introduced in Zappalà et al. (1990). This mathematical quantity uses the parameters \(k_1=5/4,\,k_2=2,\,k_3=2\). Essentially, Eq. (2) takes a scaled first-order difference between two sets of proper orbital elements. After scaling a euclidean-distance calculation from the weighted difference, Eq. (2) gives a quantified orbital distance in proper-elements space.

For a pair of proper elements \((a_{p,A},\,e_{p,A},\,si_{p,A})\) and \((a_{p,B},\,e_{p,B},\,si_{p,B})\), the quantity \(a_{\text {ave}}\) in Eq. (2) is the average of the proper semimajor axis \(a_{\text {ave}}=(a_{p,A}+a_{p,B})/2\). The corresponding angular velocity is given by \(n_{\text {ave}} = \sqrt{\mu /a_{\text {ave}}^3}\), with Earth’s gravitational parameter \(\mu \). The difference in proper semimajor axis is denoted by \(\delta a = a_{p,A}-a_{p,B}\). Similarly, \(\delta e\) is the difference in proper eccentricity \(\delta e = e_{p,A}-e_{p,B}\) and \(\delta i\) is the difference in proper inclination \(\delta i = i_{p,A}-i_{p,B}\) with transformation \(i_{p} = \sin ^{-1} (si_p)\).

We notice that Zappalà et al. (1990) does mention that the three weights \(k_1,\,k_2,\,k_3\) are empirical, but these values have been largely accepted and used in asteroid-families research. This orbital distance is linked with the velocity increment causing orbits separation. The used parameters result in an overestimation of the real separation velocity as well as unbalanced weights on different velocity directions as explained in the asteroids study (Zappalà et al. 1990).

5 Neural-networks-augmented DBSCAN

The current state-of-the-art asteroid families clustering method is described by Zappalà et al. (1990). Despite the details in the two introduced clustering methods, hierarchical clustering and wavelet analysis, there is a general procedure that these two methods both follow: first, calculate the orbital distance in proper-elements space; second, find a threshold through trial and error that would generate reasonable clusters. Although this threshold could vary from debris family to debris family, this cut-through value is deterministic within one. The objects with calculated orbital distance below this threshold will be grouped into one cluster/asteroid family. Still, these thresholds for different families are generally empirical.

This general procedure is valid especially from a model perspective. A breakup would generate a random velocity and this velocity would change the dispersed objects’ proper elements in one family. There should be an upper limit for this breakup velocity in reality as the number of dispersed objects are limited. This breakup velocity could be represented by orbital distance. As long as we could find this upper limit and threshold, we could find all the objects from the same breakup.

The situation grows complicated when there are more than one breakup event in the same space. Moreover, this threshold from each dispersion could vary dramatically. Thus, it could be the case that two or more objects from different families have relatively small breakup velocity but high breakup velocities with respect to their original families objects, and these objects would not be associated with their families but linked with each other as they will “appear” like a family. Let alone the tedious and hand-wavy threshold-setting procedure conducted case-by-case under pure intuition. Wavelet analysis is more complicated than hierarchical clustering in that it could consider “density” in the orbital-elements domain. On the other hand, it would filter non-essential “noisy” asteroids out of the influence of threshold settings. Yet, as this wavelet analysis involves Fourier analysis and fundamentally requires large amounts of data, it is not as popular as hierarchical clustering. Last but not least, one could also challenge the empirical orbital-distance definition since the parameters in the orbital-distance expression \(d_Z\), as mentioned in their defining paper (Zappalà et al. 1990), are not fix values.

Inspired by the wavelet analysis’ concept of filtering out “noise,” we adopt density-based spatial clustering of applications with noise (DBSCAN) to treat searching through the proper-elements space. The classical DBSCAN application still requires the empirical cut-through threshold and empirical orbital distance as these two aspects are tangled with each other. Namely, a changing definition (the weights parameters) of orbital distance will change the associated threshold. Furthermore, we introduce a deep-learning procedure to replace the threshold so that a trained neural network would decide if two objects belong to the same debris family based on information from all the known debris families. This neural network is a first- and second-order difference surrogate model, which could include higher order differences.

5.1 DBSCAN

DBSCAN is an algorithm that is designed to discover clusters and noise in data. The three kinds of points that could be identified are: core points, border points, and noise points (Ester et al. 1996). Both core points and border points would be considered as part of a cluster while noise points are not. We refer readers to either DBSCAN’s original paper (Ester et al. 1996) or the MATLAB DBSCAN function websiteFootnote 1 for more details about the algorithm.

Besides a function defining distance, the two parameters \(\epsilon \) and minpts are used to define the radius of a search area and the minimum-neighbor numbers in this search area. These parameters are used to determine whether a point belongs to a cluster or a point is a noisy data. By choosing a larger \(\epsilon \), if we keep the same minpts, the number of core points is going to increase correspondingly; a smaller \(\epsilon \), the number will decrease. A rule of thumb for choosing an \(\epsilon \) value is on finding the kneeFootnote 2 for an ascending two-point distance plot. As for the choice of minpts, it would rely on trial and error. Yet, however the combination of \(\epsilon \) and minpts changes, the criteria is nevertheless a deterministic one. This deterministic feature is not ideal for realistic cases as there are many uncertainties and this could be improved by incorporating nonlinear surrogate functions by including neural network.

5.2 Learning-based-threshold criteria and neural networks

The universal-approximation theorem tells us that a neural network can approximately approach any function. This means that neural network has the ability to represent a highly nonlinear function. Instead of building a deterministic threshold decided by \(\epsilon \) homogeneously in the proper \((a,\,e,\,si)\) space, a nonlinear criteria is built by neural networks. This criteria is a replacement of the deterministic threshold that is used to decide if a candidate point has enough neighbors to be considered as part of a cluster. While keeping the general structure of DBSCAN, we call this neural-networks-augmented DBSCAN.

Figure 4 shows the change of the threshold between a classical DBSCAN (left) and a neural-networks-augmented DBSCAN (right). In this schematic plot, the target point to be determined is labeled in orange; the search area determined by a threshold \(\epsilon \) is marked by a dashed orange circle; the nonlinear threshold and search area learned by neural networks is shown by the dasher purple line; neighbor objects from both method is shown in green and noise objects in grey.

In the classical DBSCAN, any point with a distance less than the predefined threshold is considered a neighbor of the target point. This distance is calculated by a predefined distance function such as the orbital distance in Eq. (2). This deterministic criteria relies heavily on a homogeneous distribution assumption. No matter of the relative orientation of two points or a single point absolute position, as long as their relative distance is the same, the relationship between these two points do not change. A neural networks based threshold challenges this homogeneous assumption. It learns a threshold based on given data thus if the given data’s “threshold” is homogeneous it would learn the nonlinear homogeneous threshold. Theoretically, it would support a much more complex threshold.

Fig. 4
figure 4

The criteria boundary change from a deterministic \(\epsilon \) threshold (left) to a NN-based nonlinear threshold (right). The figure on the left is a graphical sketch for the classical DBSCAN neighbor determination process by threshold \(\epsilon \). The figure on the right is for neural networks augmented DBSCAN with an idealized nonlinear threshold marked by a purple dashed line. DBSCAN’s initial searching circle (orange dashed line) is maintained to keep the same searching process in the neural networks augmented DBSCAN

Full size image

For two sets of proper elements, orbital distance from Eq. (2) takes the form of a scaled 2-norm calculation due to the existence of \((\delta a)^2,\,(\delta e)^2,\) and \((\delta i)^2 \) (which could be retrieved by converting the usual proper si term by taking the inverse sine function \(i=sin^{-1}(si)\)). \(\delta \) sign means taking the difference between two values (Wu and Rosengren 2021; Zappalà et al. 1990).

In building the neural networks’ input layer, we extend Eq. (2)’s input form to include not only the scaled 2-norm form but also the scaled 1-norm form. The input layer then takes the scaled difference and the scaled square of difference between two proper elements sets. So it follows the form show in Eq. (3). For two proper elements sets \((a_1,\,e_1,\,si_1)\) and \((a_2,\,e_2,\,si_2)\), the parameters expressions are \(a_{\text {ave}} = (a_1 + a_2) / 2\) and \(|\delta i |= |sin^{-1}(si_1) - sin^{-1}(si_2)|\) in radian with \(|\cdot |\) indicating taking absolute value. In our neural networks model, after two hidden layers, the output is a binary classification in Eq. (4). The classifier determines whether these two proper elements sets are from either the same debris family or not the same. The neural networks parameter details could be found in Fig. 5 and its captions.

$$\begin{aligned} \text {Input layer: }&=\left\{ (\delta a / a_{\text {ave}})^2,\,(\delta e)^2,\,(\delta i)^2,\, |\delta a / a_{\text {ave}}|,\,|\delta e |,\, |\delta i |\right\} \end{aligned}$$
(3)
$$\begin{aligned} \text {Output layer: }&= {\left\{ \begin{array}{ll} 1, \quad \text {same family}\\ 0, \quad \text {not same family} \end{array}\right. } \end{aligned}$$
(4)
Fig. 5
figure 5

Neural networks structure for debris family classification. Input layer has 6 neurons. First hidden layer has 100 neurons with ReLU activation functions. Second hidden layer has 50 neurons with ReLU activation functions. Output layer takes 2 neurons with softmax activation functions

Full size image

We notice that the scaling factor \(a_{\text {ave}}\) used in Zappalà et al. (1990) is a good choice to balance the input features and help with neural networks training convergence. Moreover, through our experiments, it is not a good choice to train a neural networks based on raw proper elements sets \((a_1,\,e_1,\,si_1)\) and \((a_2,\,e_2,\,si_2)\) without taking at least the scaled 1-norm preprocessing. The structure of this neural networks does not need to be unique as long as it works as a surrogate orbital-distance function in the form of Eq. (5). With a preprocessing \({\mathcal {P}}\) generating features over two sets of proper elements, a neural network is built upon these input features and is employed to make a classification as the final result.

$$\begin{aligned} d_{DNN} = \mathcal{N}\mathcal{N}\left( {\mathcal {P}}\left( a_1,\,e_1,\,si_1,\,a_2,\,e_2,\,si_2 \right) \right) . \end{aligned}$$
(5)

This surrogate orbital-distance function \(d_{DNN}\) is different from other orbital-distance functions like \(d_Z\) in Eq. (2) in that it combines the problems of defining a distance criteria and setting a threshold into one neural networks learning process. We call this \(d_{DNN}\) orbital distance decision from neural networks. This orbital distance decision from neural networks turns the requirement of a deterministic orbital-distance criteria to a nonlinear classification decision. It mimics the historic process of empirically finding scaled parameters for orbital distance functions by mining some known information (Southworth and Hawkins 1963). More importantly, this neural network provides a way to work around the homogeneous threshold constrain.

We generate the training data from known space-debris families. A consideration is put into the data selection due to the existence of drag. For any object in creating the training data, we extract its proper elements based on the earliest TLE from Space-Track. The proper elements obtained in this way is influenced the least from drag among all the TLEs of this object. First, we randomly select two debris objects from all space debris no matter of their origins. Then, we generate the first- and second-order differences. We also note down the label as whether these two debris are from the same debris family or not. There are 20, 000 training data generated in this way. To keep the balance of the training data, we manually add another 20,000 training data specifically from the same but randomly selected debris family. At last, we have a 40,000 training data with an almost even distribution of being from the same family and not the same family. With optimizer adam and a 70% and 30% cross-validation, our network training ends with about 90% accuracy. We deem this is enough considering there may be wrong classification in the original Space-Track data.

6 Investigation on analyst objects of unknown origins

A “well-tracked analyst object" is an object in orbit with unknown origin. They have been consistently tracked by the U.S. Space Surveillance Network (SSN) but can not be associated with a specific launch or a parent object. These objects of unknown origin are not entered into the satellite catalog, but are maintained using satellite numbers between 80000 and 89999. Historically the United States Space Command (USSPACECOM) has not published these in the Satellite Catalog (SATCAT) due to the absence of critical data points such as launch date or launch country. Their Two Line Elements (TLEs) within any last 30 days could be traced from Recent ELSETs page of Space-TrackFootnote 3 and single object could be traced by Space-Track’s ELSET Search page.Footnote 4 Well-tracked objects are generally objects that have been consistently tracked by the SSN for longer than 6 months that do not frequently cross-tag with other objects. Their TLE quality should be at the same level as other objects with known origins so we take them as they are.

Then, based on proper elements and orbital distance in proper-elements space defined by Eq. (2), we implement classical DBSCAN to find clusters and noise for the identification of debris families. The distribution of analyst objects with semimajor axis between 7100 km and 9000 km are clustered and illustrated in two separate spaces: semimajor axis vs eccentricity (a vs e) space on the left and semimajor axis vs sine of inclination (a vs si) space on the right. A total of 300 analyst objects are analyzed. We show the result with parameters \(\epsilon =0.3,\,\text {minpts}=6\) in Fig. 6. Changing those parameters, including orbital distance, could effectively change the clustering result. Figure 6 shows 5 clusters for analyst objects. Objects in the clusters have their proper elements highlighted by different colors according to their debris families. All the rest objects are recognized as “noisy” unrelated objects.

Fig. 6
figure 6

The DBSCAN clustering result by setting the parameters \(\epsilon =0.3,\,minpts=6\). The clusters are marked with different colors while noise objects out of those clusters are displayed almost in transparent

Full size image

The parameters are chosen based on a visual comparison with results under other parameters settings. Here we show the result of choosing \(\epsilon =0.5,\,\text {minpts}=6\) in Fig. 7. The “noisy” objects are similar in Figs. 6 and 7. This is a preferred result as the similar results work as reinforcing confirmation over those objects’ classification. We also notice that the purple cluster (cluster No. 5 in Fig. 6 and cluster No. 3 in Fig. 7) stays the same. Despite the similar noise objects classification, there is a major difference clearly standing out, and many subtle differences too. First, the total number of clusters reduces from 5 to 3. Clusters No. 3 and No. 4, in Fig. 6, merges into cluster No. 5 because of the increase in \(\epsilon \) from 0.3 to 0.5. This is the most distinct change. For subtle changes, there are multiple new objects included in both cluster No. 2 and No. 1 in Fig. 7.

Fig. 7
figure 7

The DBSCAN clustering result by setting the parameters \(\epsilon =0.5,\,minpts=6\). The clusters are marked with different colors while noise objects out of those clusters are displayed almost in transparent. Instead of 5 clusters in Fig. 6, the simple \(\epsilon \) change reduces the number of clusters to 3. Two clusters in Fig. 6, cluster No. 3 and cluster No. 4, are absorbed into cluster No. 1. Besides, DBSCAN’s robust noise classification ability is shown by the almost similar noise classification results

Full size image

We have shown the different performance under different parameters settings. Generally, it would be up to personal experience and analysis to choose the “right” parameters. In fact, this is also a learning process. Here, we implement the neural networks augmented DBSCAN. The neural networks learn the decision criteria based on known debris families. Moreover, the debris families data come from the same Space-Track database, the same source with analyst objects’ TLEs. They should share a similar generating routine and similar unknown “error” boundaries.

Focusing on clusters, we avoid showing the noise in the following analysis, Fig. 8 shows the different performance of classical DBSCAN clustering algorithm and neural networks augmented DBSCAN clustering algorithm. With the same parameters setting \(\epsilon =0.5,\,\text {minpts}=6\), light squares are the clustering result by classical DBSCAN and dark circles are the clustering result by neural networks augmented DBSCAN. Classical DBSCAN has only three clusters in total and has mis-combined the blue dots cluster and green dots cluster (in the middle part of \((a,\,si)\) plot) with the red dots cluster (at lower part of \((a,\,si)\) plot). Neural networks augmented DBSCAN could distinguish the different clusters.

Fig. 8
figure 8

The comparison of clustering difference from the classical DBSCAN and the neural networks augmented DBSCAN under the same parameters setting \(\epsilon =0.5,\,minpts=6\). Clustering results are marked in different colored squares for the classical DBSCAN and in circles for the neural networks augmented DBSCAN. The classical DBSCAN, the same as shown in Fig. 7, finds 3 clusters (red, yellow, purple squares) in total. Meanwhile, the neural networks augmented DBSCAN finds 5 clusters (red, yellow, blue, green, purple circles). These clusters are the same as in Fig. 6’s setting of a lower \(\epsilon \), equivalently a more strict threshold

Full size image

Besides the difference in total numbers of clusters, there is some subtle difference in terms of whether to include a specific object or not in a cluster. The object, at around a at 7600 km, e at 0.07, si at 0.88, is considered part of the cluster by classical DBSCAN, yet it is not considered part of the cluster by neural networks so it is excluded from the clusters by neural networks augmented DBSCAN. So is another point at around a at 7800 km, e at 0, si at 0.96.

The common fact that these two points locate at the edge of their corresponding clusters lead us to believe this is due to the different “thresholding” criteria. While classical DBSCAN uses an ubiquitous cutoff threshold heavily dependent on the parameters setting, neural networks classification provides more flexibility over the threshold based on learning known clusters. In fact, the assumption for an ubiquitous threshold is a homogeneous eject velocity field. Neural networks do not make this assumption but just get its nonlinear criteria based on what it has been fed.

Moreover, neural networks augmented DBSCAN is robust to parameters change. We show the following plots for parameters settings \(\epsilon =0.6,\,minpts=6\) for Fig. 9 and \(\epsilon =0.4,\,minpts=6\) for Fig. 10. Clusters from neural networks augmented DBSCAN are shown in circles with different colors. Under the same color schemes but lighter, clusters from classical DBSCAN are shown in squares. In Fig. 9, as \(\epsilon \) increases, comparing with classical DBSCAN clusters in Fig. 8, cluster No. 2 in Fig. 8 merges into cluster No. 1. This is clear on the \((a,\,e)\) plane as the squares in it have almost no other colors than red. This phenomenon is expected as \(\epsilon \) gives a single threshold to determine potential neighbors and neighbors are determined solely on this threshold despite of its physical meaning. So we should expect seeing larger clusters with an increasing \(\epsilon \) statistically even though the clusters deviate from their physically meaningful representation. In Fig. 10, with a smaller \(\epsilon \), we see the birth of cluster No. 2 (yellow square) and cluster No. 4 (blue square) in classical DBSCAN’s clustering result. As \(\epsilon \) gets close to the “good” value we found in Fig. 6, the clustering results under these different parameters are becoming much more similar. The separation of cluster No. 2 (yellow square) and cluster No. 4 (blue square) from cluster No. 1 (red square) is cheerful as they support the clustering result in Fig. 6 with an evidence of consistency.

However, no matter how parameters affect classical DBSCAN, neural networks augmented DBSCAN yields an almost invariant clustering result. The five clusters marked by circles in red, yellow, green, blue, and purple do not change no matter how corresponding classical DBSCAN clusters change due to different parameters. The clustering consistency is ideal as they serve as supporting evidences to each clusters. Because we do not use empirical orbital distance function and the clusters are consistent in large, the parameter choice for neural networks augmented DBSCAN is just for the searching process, comparing to its importance in classical DBSCAN’s clustering decision.

Besides, the physical meaning supporting these clusters is a “learned” similarity from existing debris families clusters. Instead of setting homogeneous threshold, the neural networks augmented DBSCAN does not reject the possibility of nonhomogenous threshold. Although choosing parameters based on “personal experience and sophistication” is not totally strange to neural network’s learning process, it is ideal to discard the empirical orbital distance, which is hard to do in classical DBSCAN due to the complex intertwine between choosing thresholds and designing orbital distance. There are tiny member changes within the large clusters. This is still inherent to DBSCAN’s searching process. We keep DBSCAN’s searching process so that we could efficiently search through the proper-element space with its noise identification and filtration feature.

Fig. 9
figure 9

The comparison of clustering difference from the classical DBSCAN and the neural networks augmented DBSCAN under the same parameters setting \(\epsilon =0.6,\,minpts=6\). Clustering results are marked in different colored squares for the classical DBSCAN and in circles for the neural networks augmented DBSCAN. The slight increase in \(\epsilon \), comparing with Fig. 8, changes the classical DBSCAN clustering result yet has minimal influence over the neural networks augmented DBSCAN clustering result. The classical DBSCAN, under this parameters setting, generates only 2 clusters (red and purple squares); the yellow clusters from Figs. 7 and 8 are absorbed into the red cluster because of the existence of a bigger \(\epsilon \). At the same time, the clusters from the neural networks augmented DBSCAN do not change in number. The difference is only at a single-debris level that whether or not specific debris should be involved in or excluded from a cluster

Full size image
Fig. 10
figure 10

The comparison of clustering difference from the classical DBSCAN and the neural networks augmented DBSCAN under the same parameters setting \(\epsilon =0.4,\,minpts=6\). Clustering results are marked in different colored squares for the classical DBSCAN and in circles for the neural networks augmented DBSCAN. The slight decrease in \(\epsilon \), comparing with Figs. 8 and 9, still changes the classical DBSCAN clustering result yet again has minimal influence over the neural networks augmented DBSCAN clustering result. The classical DBSCAN, under this parameters setting, generates now 4 clusters (red, yellow, blue, purple); the blue cluster appears as a stricter \(\epsilon \) is chosen comparing with in Fig. 7. But the green cluster in Fig. 6 is still classified as part of the big red cluster even though there is only a small difference in \(\epsilon \) between these two figures. Once again, clusters from the neural networks augmented DBSCAN do not change in numbers. The difference is still at a single-debris level that is the involving or excluding of some specific debris

Full size image

7 Conclusion

We have shown the first large-scale proper elements calculation in the circumterrestrial space and the first debris-families analysis. With a short introduction to the concept of proper elements and literature reviews on its application in associating debris-breakup events, we create a proper elements database that incorporates 3000 debris of known origins based on Space-Track TLE data. The examples of major debris families including CZ-4, DMSP 5D-2 F13, NOAA 16, OPS 4682, THORAD AGENA D are illustrated. Their overlaying in the proper-element space is discussed. In addition, the first debris-families distribution over well-tracked unknown origin objects/analyst objects from Space-Track is presented.

In the debris and debris families association process, we introduce DBSCAN (density-based spatial clustering of applications with noise) for finding clusters in proper-element space. We present popular orbital-distance functions and experiment with parameters settings to create a reasonable clustering for analyst objects. Then, we propose neural networks augmented DBSCAN to overcome the difficulties not only in classical DBSCAN but also in any similar criteria-based clustering schemes. We break the limits of empirical distance-function design and homogeneous-threshold choice by combining the problems of defining a distance criteria and setting a threshold into one neural networks learning process. Furthermore, our proposed surrogate orbital-distance function changes the problem from defining a function to directly creating decision from neural networks. Our proposed neural networks augmented DBSCAN provides one possible answer to the classical question raised in Zappalà et al. (1990): sensitivity study over parameter choices in cluster algorithm.

This work is a long-standing and cooperative research program that will continue in the future. We plan: (i) to further extend sets of proper elements from known debris and debris families, in such a way that we could facilitate more researches in space debris and we could increase the “learning experiences” for any similar neural networks based decision; (ii) to explore other techniques such as support vector machine (SVM) to replace the neural networks in the proposed neural networks augmented DBSCAN; (iii) to carry out more systematic study over the debris families found in analyst objects, so that the origins of some of them could be found and their breakup events could be retrieved. Debris families will hopefully make us understand better the nature of debris and their evolution.