What is the probability of impact?

Question

Space Track is a website dedicated to passing along information about satellites to consumers, primarily satellite operators. The information comes from the US military. In addition to providing the Two Line Elements that inform operators of exactly where their satellites are, they also produce conjunction reports, of predicted near misses to other objects in space.

I am trying to use the information they provide in order to predict with a medium fidelity the probability of impact. Some of the information they provide is:

• The In Track, Radial, and Cross Track miss distances.
• A 3x3 correlation matrix with the above values
• Radar Cross Sectional data for each object.
• Other miscellaneous data.

I'm trying to figure out the following:

1. How can I determine a minimum distance apart? It's clear if the two objects are 1 m diameter spheres, if they are less than 1 meter from each other, it is a hit, but is there a good rule of thumb for an arbitrary object, and a known object (Our satellite) without extensive modeling applied?
2. Given that I have a concept on how much of a miss distance is considered a hit, how do I determine that miss distance?

Answer 1

I... like this question. First things first, TLE sets are nearly useless for precise probability of collision estimation, because they do not contain an estimate of the uncertainty in the orbit, which is ultimately what affects your probability of collision.

So, assuming you get a full message with a high-precision state vector and an associated covariance matrix (full 6x6 state is better, but 3x3 can work) for each satellite, there are a few things you can do. There are quite a few different methods for estimating this, look up work by Alfano and Patera for a couple examples.

What the problem mostly boils down to is estimating the overlap of the uncertainty at the time of closest approach. Imagine "uncertainty" (in 3 dimensions) as an ellipsoid - if you're familiar with concepts of probability, it's your multivariate probability density. Now, the probability of collision can be thought of as the amount of overlap between those two ellipsoids. In your specific situation, you do not have full state and uncertainty information for both objects, but similar concepts apply here (they essentially make some assumptions and combine the uncertainty of the objects).

To try to answer your questions more directly:

1. I assume you're asking what minimum distance you (as an operator) should tolerate. This is hard to put a specific value on. Normally, you get a probability of collision estimate and assign a threshold in that space, such as "anything higher than a 1e-4 probability of collision and we will maneuver". As far as what distance constitutes a "hit", what's usually done is assuming both objects are spherical. This can be hard if you don't know much about the other object.

2. See above for the basic concepts, but implementation obviously requires knowledge of probabilistic concepts, orbit dynamics, and (preferably) state estimation techniques.

Finally, this is very much an active topic of research, and there are often papers describing a "new" way to do this. I should also mention that while what I explained is often done in practice, there are several reasons to believe that it does not, in fact, provide an accurate probability of collision estimate.

Unfortunately, open literature on this topic is sparse, but searching for the authors I mentioned above, and in fact just the terms "satellite probability of collision" will get you some very useful hits if you have the right access (AIAA journal papers, for instance).

---

However, you can estimate this numerically instead of analytically using Monte Carlo methods. Here's a quick and dirty way to do this:

1. Sample from your "primary" satellite (the one you got the conjunction report for). To do this, generate a 3x1 normally-distributed random vector with zero mean and unit variance (the standard normal. Now, take your 3x3 covariance for this object, find its Cholesky decomposition (you should be able to find a routine to do this in your software package - 'chol' in Matlab), and multiply it by your random vector. You've now sampled your probability space.

2. Repeat #1 for the conjuncting object. However, once you sample, you need to translate each point according to its position, which is given to you in RIC coordinates. Just add your position vector to the random vector.

3. Calculate the distance between your two samples. Record said distance.

4. Repeat steps 1-3 for some very large number of samples (think tens of thousands - generally, the more, the merrier).

5. Now, you have a list of "miss distances". There are a couple things you can do with this, but one thing that's particularly helpful is to construct an empirical cumulative distribution function. You can look up how to do this, but you should be able to find a function to do this in most software packages (in Matlab, it's 'cdfplot'). This will give you a plot that looks like this (only with x>0):

Now, one way you can read this is by choosing a "miss distance" you're comfortable with (say, 20m), and finding the value of $F(x)$ at this value of $x$. That's the probability that you will come within 20m or less of the offending satellite. If it's greater than some threshold you've decided upon with your management, you begin the fun task of planning an avoidance maneuver. If not, you wait and watch. Alternatively, you can specify a probability you're comfortable with, and read off a miss distance associated with it.

One thing you'll notice is that you're really focused on the tail of your plot... and you'll see that to get a smooth representation at the tail, you need a whole hell of a lot of samples. So like I said, the more, the merrier.

One last thing: if you've followed these steps, the probability number you have is $1\sigma$. For, say, a $3\sigma$ estimate, just multiply your two covariance matrices by your factor (in this case 3).

Answer 2

The equation for the distance between the two objects is simply the magnitude of the vector between their positions. Numerical root-finding is made easier by knowing the analytic derivative, which is $$\newcommand{\r}{\mathbf{\vec{r}}} \newcommand{\v}{\mathbf{\vec{v}}} \frac{d}{dt} \sqrt{\r\cdot\r} = \frac{\v\cdot\r + \r\cdot\v}{2\sqrt{\r\cdot\r}} = \frac{\r\cdot\v}{||r||}$$ where $\r$ is the relative position and $\v$ is the relative velocity. To write that as a function of time, which you can minimize to find the time of closest approach (TCA), you need the two state vectors, a propagator to give you a set of points in time for each, and an interpolator to help the minimizer iterate to convergence between the propagated points. There are lots of choices for interpolation. My favorite is Hermite's method, since the propagator already gives us the velocity and position together. I like first derivative on four points (two before and two after, giving a 7th-degree polynomial) but your mileage may vary.

The result is how close the two objects would pass by each other if you knew everything exactly, which of course you don't, so you have to include uncertainty in the model. The way we describe uncertainty in astrodynamics is with covariance matrices. Using them, we express the variables in the equations as standard normal probability distributions. Integration over those distributions, as in the one equation quoted below, calculates the probability that the distance of closest approach (DCA) is less than some user-specified value of concern. Every covariance describes a probability distribution over the values of the associated parameters. The computation can be done in a variety of ways, including the monte carlo method described in another answer to this same question. Position and velocity are the most commonly seen, but anything that might be unknown has its own covariance (which includes its correlations with the other variables), and should be integrated against a probability distribution broadly representative of its associated uncertainty. In this particular problem, if one or both of the objects are not spheres, then attitude estimation and its uncertainty can play a big role if you wish, or you may assume them away for convenience.

At that level of detail, all algorithms are the same. The differences come from exactly which choices each author makes to simplify the calculation into something that can be done relatively quickly, without requiring too much data collection or sacrificing too much accuracy. The method described in the document "Probability of Collision in the Joint Space Operations Center", dated 24 June 2016 (hereafter PCJ) can be used to illustrate one way this might be done, but there are many other possibilities.

Page 3 of PCJ says covariance is distributed in the format described in the Consultative Committee for Space Data Systems (CCSDS)'s Recommended Standard for Conjunction Data Messages (CCSDS 508.0-B-1, hereafter CDM). The best place to start reading it is at the end, with Annex E (pages 66-70 of 72), since that defines the terms used throughout the document. The radial, in-track, and cross-track values in the conjunction message are the components of the vector having DCA as its magnitude when time is TCA. The way those terms are used in both CDM and PCJ, the radial (R) direction is exactly radial, so the in-track (T) direction is not exactly aligned with velocity, and cross-track is called N (normal to the orbit plane).<1>

CDM does not specify how collision probability shall be calculated. Instead, it gives a bibliography of five different methods, and a procedure for registering new options for the method keyword. As of now, the full list at the Space Assigned Numbers Authority (SANA) has 14 options. Of these, PCJ's own bibliography lists Foster 1992 (Foster, J.L., and Estes, H.S., A Parametric Analysis of Orbital Debris Collision Probability and Maneuver Rate for Space Vehicles. NASA/JSC-25898. Houston, Texas: NASA Johnson Space Flight Center, August 1992) and Chan 1997 (Chan, K. Collision Probability Analyses for Earth Orbiting Satellites. In Space Cooperation into the 21st Century: 7th AAS/JRS/CSA Symposium, International Space Conference of Pacific-Basin Societies (ISCOPS; formerly PISSTA) (July 15-18, 1997, Nagasaki, Japan), edited by Peter M. Bainum, et al., 1033-1048. Advances in the Astronautical Sciences Series 96. San Diego, California: Univelt, 1997). The newer Alfano, Alfriend, McKinley, and Patera methods are not mentioned by PCJ, but some are in CDM.

Quotations from PCJ which directly answer your specific questions include:

determine the longest distance the two satellite’s centers of mass can be apart and still have the two satellites touch. This defines the "specified distance apart" that is used in the $P_c$ calculation. Note that if the two satellites are not spheres, then a simple change in orientation means that they may not touch and no collision would occur.

circumscribe the primary and secondary objects each by a sphere, add the two sphere radii to generate a supervening sphere that can contain both circumscribing spheres, and project this supervening sphere into the conjunction plane as a circle.

Pre-assigned default values for payloads and platforms (5 meters), rocket bodies and unknown objects (3 meters) and debris (1 meter) were determined through a study of sizes of objects in the space object catalog and are normally used.

The covariance matrix for each object is referenced to its own RTN coordinate frame. For each object matrix components are computed using 5-point Lagrange interpolation of the covariance in the ephemeris file produced by the JSpOC, [assuming] Primary and secondary errors are independent allowing "combined" covariance to be the simple sum of the individual covariances (in a common frame)

Computation of $P_c$ takes place in the collision plane... perpendicular to the relative velocity vector at TCA. This reduces the mathematics from 3D to 2D. the equation used to compute $P_c$ is:

$$\frac{1}{2\pi\sqrt{|C|}} \iint_{x^2+y^2 \leq d^2} \exp \left( -\frac{1}{2} (\r-\r_{SP})^T C^{-1} (\r-\r_{SP}) \right) dx\, dy$$

where $C$ is the 2X2 projection of the combined 3X3 covariance at TCA onto the collision plane, $|C|$ is the determinant of $C$, $C^{-1}$ is the inverse of $C$, $d$ is the sum of the two object sizes, $\r = (x, y)^T$ is any point in the collision plane such that $x^2 + y^2 \leq d^2$, and $\r_{SP}$ is the position of the secondary relative to the primary along the x-axis in the collision plane.

The JSpOC uses error functions (ERF) for computing the double integration in the $P_c$ equation. In addition the JSpOC performs integration over a square circumscribing the circle of radius $d$. This square is aligned with the axes of the combined 2D probability density function in the collision plane. This simplifies the computation of $P_c$ but gives a very slightly larger value.

Every assumption made by PCJ intentionally overestimates the size of the objects, and therefore intentionally overestimates the risk of collision, to err on the side of sending too many warnings rather than too few. If you actually know what the size, shape, and orientation of the two objects are, you could add that to your calculation, as Georges Krier does in Satellite collision probability for long-term encounters and arbitrary primary satellite shape (2017), but to use shape information effectively, you need to have a good estimate for both the attitude states and the covariances in them.

Footnote <1>: CDM names two different frames the words "radial" and "in-track" might mean, but only one of them is actually used in the message format. The one they use they call "RTN", for "Radial, Transverse, Normal". Normal means unit vector parallel to the angular momentum (position cross velocity), Radial means parallel to the vector pointing from the central body to the orbiting object (equivalently, from the object away from the central body), and Transverse means the unit vector that completes the right-handed system, which points in the plane of the orbit somewhere close to but not exactly coinciding with the object's velocity, except at apogee and perigee, or if eccentricity is zero. The name they give to the other frame they don't want you to use is "TVN", meaning "Transverse, Velocity, Normal", where Normal is the same, Velocity exactly coincides with the actual instantaneous velocity direction, and Transverse still means complete the right-handed system, but that means it points in the plane of the orbit somewhere close to but not exactly coinciding with the object's outward radial position vector, except at apogee and perigee, or if eccentricity is zero. Naturally, these names differ between which author you read! What CDM calls RTN (used) is called UVW by PCJ (footnote 13, page 3) and RSW by Vallado (page 157 of the 4th edition), while what CDM calls TVN (not used) is called PTW by PCJ and NTW by Vallado.

Answer 3

The document https://www.space-track.org/documents/How_the_JSpOC_Calculates_Probability_of_Collision.pdf was written and posted in 2016. It discusses the process for the computation of probability of collision in the Combined Space Operations Center (CSpOC) at Vandenberg Space Force Base in California.