The Cube approach

The Cube approach#

During the first decade of the new millenium computational tools became refined enough for scientists to conceive numerical simulations of the long term evolution of an orbital environment:

“With the help of modern computers, it is now possible to perform numerical simulations of the orbital evolution of an N-body system. Therefore, there is a need for a collision model that can work with an orbital evolution simulation code to allow for source and sink terms of the objects involved as well as to utilize their updated orbital elements as they evolve in time. The “Cube” approach is designed to accomplish just that.” - Liou, J-C. “Collision activities in the future orbital debris environment.” Advances in Space Research 38.9 (2006): 2102-2106.

The Cube approach mentioned in the above quote is based on the idea that if \(P_{ij}\) is a collision probability between orbiting object \(i\) and object \(j\), then the total number of collisions \(N_{tot}\) between the two objects in a time interval where \(P_{ij}\) can be considered constant will be:

\( N_{tot} = \int_{t_s}^{t_f} P_{ij} ds = [t_f - t_s] P_{ij} \)

The formula above is formally exact and shifts the problem to that of finding such an averaged \(P_{ij}\) capturing the collisional dynamics. In the Cube approach this is obtained dividing the entire Cartesian space in cubes of dimension \(L\) and introducing the following two hypothesis:

  • Two objects \(i\) and \(j\) have a non zero impact probability in \([t_f - t_s]\) if and only if they occupy the same cube at \(t_s\).

  • When non-zero, the collision probability between the two object \(i\) and \(j\) is derived from the kinetic theory of gas applied to the cube.

In the following we describe python code implementing this approach for a large orbital population. As usual, we start with some imports:

# Core imports
import pykep as pk
import numpy as np
import pickle as pkl
from copy import deepcopy
import time
from collections import defaultdict
import itertools

# sgp4 imports
import sgp4
from sgp4.api import Satrec, SatrecArray

# Plotting
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
%matplotlib inline

Data imports#

We will use, as satellite population, the current tracked objects in Low Earth Orbit as prepared by the code described in The current LEO population. At the end of that code a pickled file is created containing initial positions and velocities of the objects together with the satcat dictionary to identify them. Let us load that data.

# r is in km and v in km/s
with open("data/leo_population.pk", "rb") as file:
    r,v,to_satcat_index,satcat = pkl.load(file)

The object described by the entry satcat[to_satcat_index[j]] has initial position r[j] and initial velocity v[j]. Let us inspect one entry

print("Dictionary entry: ", satcat[to_satcat_index[3685]])
print("Position (km): ", r[3685])
print("Velocity (km/s): ", v[3685])

Dictionary entry:  {'CCSDS_OMM_VERS': '2.0', 'COMMENT': 'GENERATED VIA SPACE-TRACK.ORG API', 'CREATION_DATE': '2022-02-03T04:23:25', 'ORIGINATOR': '18 SPCS', 'OBJECT_NAME': 'SL-14 DEB', 'OBJECT_ID': '1978-100F', 'CENTER_NAME': 'EARTH', 'REF_FRAME': 'TEME', 'TIME_SYSTEM': 'UTC', 'MEAN_ELEMENT_THEORY': 'SGP4', 'EPOCH': '2022-02-02T22:42:05.237280', 'MEAN_MOTION': '12.56769207', 'ECCENTRICITY': '0.00286110', 'INCLINATION': '82.3526', 'RA_OF_ASC_NODE': '268.9057', 'ARG_OF_PERICENTER': '297.9015', 'MEAN_ANOMALY': '75.7071', 'EPHEMERIS_TYPE': '0', 'CLASSIFICATION_TYPE': 'U', 'NORAD_CAT_ID': '19133', 'ELEMENT_SET_NO': '999', 'REV_AT_EPOCH': '93528', 'BSTAR': '0.02113600000000', 'MEAN_MOTION_DOT': '0.00003527', 'MEAN_MOTION_DDOT': '0.0000000000000', 'SEMIMAJOR_AXIS': '7814.445', 'PERIOD': '114.580', 'APOAPSIS': '1458.668', 'PERIAPSIS': '1413.952', 'OBJECT_TYPE': 'DEBRIS', 'RCS_SIZE': 'MEDIUM', 'COUNTRY_CODE': 'CIS', 'LAUNCH_DATE': '1978-10-26', 'SITE': 'PKMTR', 'DECAY_DATE': None, 'FILE': '3290513', 'GP_ID': '195174860', 'TLE_LINE0': '0 SL-14 DEB', 'TLE_LINE1': '1 19133U 78100F   22033.94589395  .00003527  00000-0  21136-1 0  9996', 'TLE_LINE2': '2 19133  82.3526 268.9057 0028611 297.9015  75.7071 12.56769207935282', 'RADIUS': 0.2631205051777122}
Position (km):  [2129.88342848 3695.71335128 6552.37467945]
Velocity (km/s):  [ 1.38399824  5.88779107 -3.78530589]

A few useful functions we will use#

The following function computes the orbital period via the vis-viva equation.

def period(r,v, mu):
    """Computes the orbital period from the vis-viva equation

    Args:
        r (float): The radius (in L).
        v (float): The velocity (in L/T).
        mu (float): The gravitational parameter in L^3/T^2

    Returns:
        The orbital period (in T)
    """
    En = v**2/2 - mu / r
    a = -mu / En / 2
    if a<0:
        raise ValueError("Hyperbola!!!")
    return np.sqrt(a**3/mu)*2*np.pi

The following function is a central part of the Cube approach. It takes \(N\) points and a linear dimension and returns all the points that occupy the same cube of side L. In our implementation the cubes start from the origin. In other words, along the three dimensions, the cube vertices have coordinates \([...,-L,0,L,...]\).

def cubes(points, cube_dimension):
    """Divides the space in cubes of size cube_dimension and returns points within the same cube

    Args:
        def cubes(points, cube_dimension):
 (Nx3 np.array): The cartesian position of the satellites (in L).
        cube_dimension (float): The cube dimentsion (in L).

    Returns:
        a list containing lists of satelites idx occupying the same cube
    """
    # init
    in_same_cube = []
    pairs = []
    cubes = defaultdict(list)

    # We compute the floored Cartesian coordinates identifying the bins.
    pos = points
    pos = pos / cube_dimension
    pos = np.floor(pos).astype(int)
    # We fill the bins
    for i, xyz in enumerate(pos):
        cubes[tuple(xyz)].append(i)
    # We find bins with more than one object
    for key in cubes:
        if len(cubes[key]) > 1:
            in_same_cube.append(cubes[key])
    # We compute all the combinations between the pairs
    for group in in_same_cube:
        pairs.extend(((list(itertools.combinations(group, 2)))))
    return pairs

Note that the great advantage of the above preliminar collision detection is that it has a linear complexity in the number of particles. This makes the Cube approach scale remarkably well as an n-body collision simulation. Lets time its performances in our specific case:

%timeit cubes(r, 10)

28.1 ms ± 4.92 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

The following function computes the SGP4 propagation of selected catalogued objects over a fixed time grid. This numerical propagation is here made using the vectorized version of the SGP4 propagator achieving great efficiency.

def simulate_sgp4(indexes, sim_time=20,time_interval=5, jd = pk.epoch_from_iso_string("20220301T000000").jd):
    """Computes all satellites ephemerides on a time grid

    Args:
        indexes (list of pk.planets): The indexes in the satcat dictionary of the objects to propagate.
        sim_time (float): The total propagation time (in years).
        time_interval(float): The time resolution (in days).
        jd (float): the starting epoch in julian date.

    Returns:
        a list containing lists of idx identifying the object occupying the same cube
    """
    # This list will contain all the sgp4 Satrec objects
    satellite_l = []
    for idx in indexes:
        l1 = satcat[idx]["TLE_LINE1"]
        l2 = satcat[idx]["TLE_LINE2"]
        satellite_l.append(Satrec.twoline2rv(l1, l2))
    # Here we build the vectorized version allowing for speed
    satellites = SatrecArray(satellite_l)
    # The Julian dates are from jd0 to 20 years after
    jds = jd + np.arange(0,sim_time*365.25/time_interval)*time_interval
    frs = jds * 0
    return satellites.sgp4(jds, frs)

The full simulation of 20 years of LEO evolution#

We are now ready to set up a simulation of the evolution of the LEO environment for the next 20 yars. In particular we will be only detecting collisions. We start precomputing all the positions and velocities of the orbiting objects in 20 years at a 5 day resolution.

start = time.time()
err,r_sgp4,v_sgp4 = simulate_sgp4(to_satcat_index, sim_time=20,time_interval=5, jd = pk.epoch_from_iso_string("20220301T000000").jd+np.random.random()*5)
end = time.time()
print("Time to perform all numerical integrations: ", end - start)

Time to perform all numerical integrations:  7.937011003494263

We now perfrom the simulation and use within each time interval, the kinetic theory of gasses to determine whether a collision has happened or not.

start = time.time()
undecayed = set(np.arange(r_sgp4.shape[0]))
n_collisions=0
N_tot = 0
for i in range(r_sgp4.shape[1]):
    Lcube = 10. #km
    time_grid  = 5 #days
    # If signalled from the sgp4, we remove the indices of the decayed satellites
    decayed = set(np.where(err[:,i]>0)[0])
    undecayed = undecayed - decayed
    undecayed_l = np.array([j for j in undecayed])
    # We detect all satellites couples in the same cube of Lcube km size
    pairs = cubes(r_sgp4[undecayed_l,i,:], cube_dimension = Lcube)
    #kdt = KDTree(r[undecayed_l,i,:])
    #collision = list(kdt.query_pairs(Lcube))
    #print(collision)
    for pair in pairs:
        # we get the indexes in r,v
        idx1 = undecayed_l[pair[0]]
        idx2 = undecayed_l[pair[1]]
        # we store positions and velocities from r,v
        r1 = r_sgp4[idx1,i,:]
        r2 = r_sgp4[idx2,i,:]
        v1 = v_sgp4[idx1,i,:]
        v2 = v_sgp4[idx2,i,:]
        # we get the collision radiu from debris (indexed differently hence to_satcat is used)
        radius1 = satcat[to_satcat_index[idx1]]["RADIUS"]
        radius2 = satcat[to_satcat_index[idx2]]["RADIUS"]
        # Relative velocity 
        Vrel = np.linalg.norm(v1-v2)
        # Collisional area of the couple (in km^2)
        sigma = np.pi*((radius1+radius2)/1000)**2 
        # Volume of the cube (km^3)
        U = (Lcube)**3
        # We compute the spatial densities
        # densities (from "Assessing collision algorithms for the newspace era" )
        s1 = 1./U
        s2 = 1./U
        # collision probability
        Pij = s1*s2*Vrel*sigma*U*time_grid*pk.DAY2SEC
        N_tot += Pij
        # Store
        if Pij > np.random.random():
            print(f"Collision! pair: {pair}, years: {i*5/365.25}")
            n_collisions+=1
end = time.time()
print("Time elasped: ", end - start)
print("Decayed objects: ", r.shape[0] - len(undecayed))
print("Number of collisions: ", n_collisions)
print("Expected number of collisions: ", N_tot)

Collision! pair: (15977, 16424), years: 0.8213552361396304
Collision! pair: (15935, 16388), years: 1.273100616016427
Collision! pair: (1041, 6320), years: 1.4921286789869952
Collision! pair: (13484, 13591), years: 3.1074606433949348
Collision! pair: (13689, 13950), years: 6.365503080082136
Collision! pair: (14107, 15279), years: 10.704996577686517
Collision! pair: (3897, 15176), years: 11.567419575633128
Collision! pair: (10109, 11757), years: 12.320328542094456
Collision! pair: (12350, 13019), years: 13.566050650239562
Collision! pair: (4910, 12250), years: 14.099931553730322
Collision! pair: (13633, 14012), years: 15.824777549623546
Collision! pair: (2869, 3364), years: 16.892539356605067
Time elasped:  51.43294858932495
Decayed objects:  4992
Number of collisions:  12
Expected number of collisions:  13.336892648579166

Repeating the above code will give different results as we added noise to the initial conditions and the collision probabilities are checked rolling a dice.

Contents