Forecasting, and averting, major asteroid impacts: from science fiction to reality (full version)
Armageddon, Deep Impact, Don’t Look Up, Greenland. Fossil records confirm that major asteroid collisions with Earth happen routinely, and will certainly happen again. Can we forecast such collisions with enough time to actually do something about it? UC San Diego researchers say yes.
1. Introduction. A foundational result in statistics, called the central limit theorem, is that uncertainty distributions generally tend towards “bell”-shaped “Gaussian” (aka “Normal”) distributions. [“Power laws” and systems with “self-organized criticality” are also common in certain settings.] The tendency towards Gaussianity is often visualized via a machine, called a Galton board (see above) that incrementally releases a large set of small spheres near the top/center of an array of identical pegs arranged in a vertical plane. Each sphere bounces its way randomly past each peg that it encounters in the array, going left or right with essentially equal probability at each bounce; after thousands of spheres are released, nearly identically, the resulting distribution of where the spheres wind up, discretized into a finite number of bins at the bottom of the machine, is almost always, quite accurately, shaped something like the cross-section of a “bell”. Mind-blowing, when you first see it in real life.
The assumption that the distribution of unmodeled “disturbances” in nature and finance are usually bell-shaped is ubiquitous, and fits well in many practical circumstances, such as wind speed variations, the demand for airline tickets from LAX to NYC next Tuesday, etc, even in higher-dimensional problems (see above). Uncertain “disturbances” affecting the evolution of most systems of interest are thus often modeled simply in terms of the “center” and the “width” (in each “principle direction”) of a corresponding Gaussian probability distribution representing the range of values that these disturbance variables might have. The customary Kalman filter then models the resulting evolution in time of the “state” of a system of interest, and the associated uncertainty thereof, assuming that these uncertain state variables, driven by Gaussian disturbances, are themselves also endowed with Gaussian uncertainty. This assumption is valid, specifically, when the small perturbations to a system are well modeled as behaving linearly. In short:
Gaussian uncertainty of disturbances + linear perturbation dynamics →
Gaussian uncertainty of the optimal state estimate.
2. Non-Gaussian uncertainty. Our research focuses on cases for which the assumption of Gaussianity of the uncertain states breaks down. This breakdown generally happens in nonlinear systems when uncertainties are large enough that nonlinearities play a significant role in the time evolution of system perturbations.
2.1 Key idea demonstrated on the (3D) Lorenz system. Our “legacy” work, in [1], illustrated graphically the problem of non-Gaussian uncertainty propagation in the estimation of a chaotic nonlinear problem in three variables, governed by the celebrated Lorenz equation. [By “legacy”, I mean the polite way that the lead grad student on the project now, Ben Hanson, refers to the work that his advisor did on the project when Ben was 12 years old…]. The dynamics of this model problem are well understood and, being only three dimensional, easily visualized.
Figure 1. Visualization of the Lorenz attractor, the prototypical “chaotic” system, in which small perturbations of the system’s initial state lead to perturbed trajectories that eventually diverge exponentially. The black curve shown represents a single trajectory of the (simple, 3-state) Lorenz system from t=0 to t=1 over the 3-dimensional (3D) phase space visualizing the evolution in time of the state of the Lorenz system. The 3D “cloud” around t=0 represents Gaussian uncertainty of the initial state (that is, the initial state should be somewhere within this ellipsoidal cloud). This uncertainty cloud evolves in time, as shown in light blue at successive time intervals up to t=1, by which time the uncertainty cloud is strongly non-Gaussian (that is, no longer shaped like an ellipsoid).
The Lorenz system discussed above is governed by a simple 3-state nonlinear ordinary differential equation (ODE), and has a 3-dimensional phase space, which is tractable to analyze using existing computational methods on current-generation laptop/desktop computers. Although the Lorenz system is an ODE, the evolution of the uncertainty “cloud” (representing the possible values of the state over phase space) subject to the Lorenz system is governed by a partial differential equation (PDE) known as the Fokker-Plank equation (FPE). Because the FPE is generally nonintegrable, the actual uncertainty cloud cannot be modeled with a finite set of parameters. Consequently, state estimation and uncertainty propagation methods aim to approximate and evolve this uncertainty using as few parameters as possible, as efficiently as possible.
2.2 The bigger picture in state estimation. One way of calculating the evolution of an uncertainty cloud, called a Monte Carlo (MC) simulation, is to populate the initial cloud with a distribution of N>>1 candidate initial conditions at t=0, each assigned the same initial probability (1/N), and then to march each independent simulation out to t=1. This works well until new measurements require accounting for using Bayes rule.
State estimation techniques, known as filters, assimilate the information from measurement corrections with their prior uncertainty predictions in a tractable manner. Filters can generally be grouped by methodology into three broad categories:
Kalman estimation methods track the mean and second-order statistics (covariance) of an uncertain state. These methods are particularly well suited to problems with Gaussian state uncertainty, and simply propagate the mean and the covariance of the uncertain state in time.
Lagrangian estimation methods march many candidate trajectories though phase space, and infer the mean and uncertainty of the estimated state at any time by examining the mean and distribution of this “cloud” of candidate trajectories as it evolves in time. Particle Filters and Ensemble Filters are the two main classes of Lagrangian estimation methods (there are many distinct variants in both classes):
Particle Filters change the weights on each candidate trajectory at measurement updates, whereas
Ensemble Filters keeps the weights on each candidate trajectory constant, and nudge the position of each candidate trajectory in phase space at measurement updates.
Eulerian estimation methods develop a grid over phase space, or a subset thereof, and model the flux of probability from each grid cell to its neighbors at each timestep.
Note that, for any Lagrangian or Eulerian estimation approach to be valid, it must be consistent with the celebrated Kalman filter in the limit of linear dynamics, high resolution, and Gaussian state disturbances and measurement noise and initial uncertainty.
Ensemble Filters are generally used in very high dimensional problems (e.g., weather forecasting), in which the sampling of the uncertain state is grossly under-resolved. Measurement updates are made based on very-low-order approximations of the second-order statistics of the system’s d>>1 dimensional state using only N=O(100)<<d candidate trajectories.
Particle Filters (PFs), on the other hand, are often used in low dimensional (e.g., astrodynamical) problems, often when the uncertainty is non-Gaussian. An inherent flaw associated with standard PF variants, known as particle degeneracy, occurs when (after many measurement corrections) the probability assigned to some candidate trajectories diminishes to nearly zero, and the probability assigned to other candidate trajectories increases. Eventually, those candidate trajectories with nearly zero probability must be pruned (eliminated) from the simulation, and those candidate trajectories with much larger probability must, to maintain resolution, be somehow “split” into smaller pieces (aka “resampled”). It is this “splitting” process of the PF that is theoretically unsatisfactory (that is, “ad hoc”). A Kalman filter variant, the Unscented Kalman Filter (UKF), is a special PF implemented using only 2d+1 deterministically-sampled candidate trajectories, where d is the dimension of the problem considered (with d=6 in astrodynamics problems), implemented under the Kalman assumption to approximate the mean and covariance of the uncertain state.
Eulerian methods numerically solve the FPE by discretizing phase space, illustrated for example in Figure 1, into a bunch of small “cells” (cubes). They then propagate the probability associated with each cell in time, by computing and applying a “flux” of probability from one cell to its neighbors at each timestep. This principle is straightforward using computational methods developed for hyperbolic systems in the field of fluid mechanics. Instead of discretizing the entire 3D volume illustrated in Figure 1 (or, the d-dimensional volume in harder problems), the key observation motivating our early work is that, at any instant, the probability associated with almost all of the cells covering the domain of interest (see, e.g., Figure 1) is nearly zero. Thus, the central idea implemented in [1] is to make a list of only those “active cells” with nonzero probability at any given time, and their neighbors, and to only calculate and apply the flux of probability between neighboring cells in this list. This general idea is called Grid-based Bayesian Estimation Exploiting Sparsity (GBEES), which is an Eulerian estimation method (aka filter) designed for efficient non-Gaussian uncertainty propagation. For more examples of where the current landscape of contemporary filters fall in this categorization, see Figure 3 of [3].
2.3 Extension of the key idea to higher dimensions. Non-Gaussian estimation problems of practical interest generally have more than d=3 states, and thus a d-dimensional phase space; the (difficult) challenge of extending our original work in 3D to such higher-dimensional problem is often called the curse of dimensionality. Eulerian uncertainty propagation methods are prone to this curse, as a computational grid grows exponentially with each added state dimension. The key challenges in the extension of GBEES to problems in d dimensions, where d>3, are threefold:
how to efficiently manage the list of active cells and their neighbors,
how to quickly eliminate all cells that aren’t actually needed from this list, and
how to accelerate the execution of the resulting algorithm on modern massively-parallel (GPU-based) computational architectures.
3. Application of GBEES to the (6D) motion of spacecraft.
The focus problems considered in our recent work relate to the evolution of the (3D) position and (3D) velocity of bodies moving around in our solar system under Newton’s law of universal gravitation; this problem thus has a 6-dimensional phase space. Our four recent published papers successfully addressed all of the key challenges mentioned in §2.3, facilitating the effective extension of GBEES, demonstrated in 3D in [1], to problems of practical interest in 6D, as illustrated in Figures 2 and 3.
Figure 2. Simulation of an inclined orbit around Enceladus in a rotating (synodic) frame, where the rotation of the frame matches the speed at which Saturn and Enceladus rotate about their barycenter, approximating both orbits as circular.
Figure 3. At time t=0, the orbiter in Figure 2 is at apoapsis, the point in the orbit farthest from Enceladus; we here initialize a Monte Carlo (MC) simulation and an Unscented Kalman Filter (UKF) with Gaussian uncertainty (above left). The dynamic model of the orbiter includes the gravity of Saturn and Enceladus. Both the MC and the UKF are propagated using this force model for 5.95 hours until periapsis, the point closest to Enceladus. As seen (above right), the state uncertainty is now highly non-Gaussian, and the UKF does not accurately represent it. It is in precisely such settings, when new measurements are accounted for, that more sophisticated estimation techniques, like GBEES, become necessary.
In [2], the first attempt to extend GBEES to higher-dimensional, astrodynamic problems was made. To address the computational complexity of the legacy algorithm, we changed the underlying data structure in GBEES from nested lists to a binary search tree (BST). Additionally, we employed directional growing and pruning of the grid, a procedure that substantially enhanced computational efficiency by incorporating knowledge of the dynamics into the direction of the grid evolution. These two changes, amongst others, allowed for feasible uncertainty propagation in 4D phase space. One astrodynamic system with this dimensionality is the Planar Circular Restricted Three-Body Problem (PCR3BP), which includes the direct effect of gravity from two major bodies (e.g., the Earth and Moon) on a third (small) body, approximating the orbits of the two major bodies as circular about their shared barycenter. The “planar” restriction shrinks the phase space from 6D to 4D. We applied GBEES to a Jupiter-Europa Lyapunov orbit with PCR3BP dynamics, and achieved modest accuracy and efficiency.
In [3], we compared the performance of GBEES to other leading non-Gaussian estimation filters on a planar Saturn-Enceladus trajectory with uncertainty that becomes highly non-Gaussian. UKF applied to this problem quickly diverged, whereas GBEES captured the non-Gaussianity accurately when compared to a large MC. The two PF variations also tested, the Bootstrap PF of Candy (2007) and the Ensemble Gaussian Mixture Filter of Anderson et al. (1999), succumbed to particle degeneracy, as defined previously, whereas the GBEES framework ensured the probability density maintained good resolution through the measurement corrections. Overall, GBEES achieved higher accuracy than the other filters evaluated, albeit at a significantly increased computational cost. The work in [3] is formalized in [4]. At this point, the computational efficiency of GBEES lagged behind that of the contemporary state-of-the-art filters mentioned above; significantly, however, the structure of the GBEES algorithm is “embarrassingly parallel”, facilitating its efficient implementation on GPUs, as explored in [5] and discussed below.
In a close collaboration with researchers at the University of León, we addressed in [5] the computational efficiency of GBEES, by optimizing the CPU implementation, and translating the algorithm to CUDA for GPU execution. To begin, we generalized the PDE marching theory from [1] to d-dimensional systems for d>3. We then further changed the GBEES data structure from a BST to a hashtable to enable execution on a GPU. Through careful consideration of the dynamic grid using heaps that store the the free/used hash indices, these modifications preserved the algorithm’s behavior on the GPU, making it function identically to the CPU implementation. To test the enhancements made, we applied the legacy, the CPU-optimized, and the GPU versions of GBEES to the uncertainty propagation of the Lorenz system from [1], and a 6D variation of the Lorenz 96 system.
The CPU-optimized version achieved a speedup of ~14x on the 3D system as compared to the legacy version. The GPU version, tested on four different NVIDIA GPUs, achieved a maximum speedup of ~133x on the 6D system relative to the CPU-optimized version, implying a speedup of over 1000x relative to the legacy version. These results demonstrate effective parallel speedup, and enable efficient uncertainty propagation for 6D systems on relatively modest modern-day computers.
To facilitate the future use of GBEES by the uncertainty propagation community, we have made freely available on GitHub both the CPU-optimized and GPU versions of GBEES, which are the result of a focused development effort. These codes are implemented in an easily extensible manner, and can be modified to implement the user’s preferred dynamics and measurement models. Included in the repositories are thorough explanations of their installation and usage.
4. Application to the (6D) motion of large asteroids. A related problem of particular interest to life on Earth is the motion of large asteroids (known in this context as near-Earth objects) through space, which if they collide with Earth can lead to massive destruction and, possibly, another mass extinction event. This class of problems focuses specifically on “unlikely” events, out in the “tail” of the associated probability distributions. The odds that any given large asteroid will hit Earth, based on what we know many years in advance of such a collision, is generally quite small, but there are lots of such asteroids, and the potential consequences of a single such collision with Earth is quite large indeed, as seen in fossil records. With a sufficiently early and accurate prediction of the probability of such events, it is possible to better prepare, and potentially to redirect the trajectory of an asteroid of concern with a small nudge, averting potential disaster, as demonstrated in the Dart Asteroid Redirection Test. Accurate prediction of unlikely events out in the (non-Gaussian) tails of probability distributions, years in advance of any possible such collision, is essential for this to be done properly; the present work addresses this need.
One near-Earth object of particular interest, dubbed 2024 YR4, currently has (according to our calculations) a 4% chance of colliding with Earth’s moon on Dec 22, 2032. YR4 currently has a negligible probability of hitting Earth on that pass (an older estimate, made in Feb 2025, put this chance at as high as 3.1%). Given the proximity of the moon to the Earth, the reader can well appreciate the importance of getting such predictions right; if YR4 were to hit Earth (in some future pass near Earth), the explosion would be equivalent, roughly, to that of 500 Hiroshima bombs.
YR4 will also (safely) pass, relatively close (within 5 million miles) to Earth on Dec 17, 2028; if as a species humans deemed YR4’s trajectory as an imminent threat to life on Earth (at this time, it is NOT), and if we were to attempt a redirection (“nudge”) of YR4’s orbit around the sun, that would be the ideal time to do it (that is, 4 years in advance of any potential collision). Our study of YR4’s potential collision with our moon on Dec 22, 2032, and what it would take to redirect it to avoid such a collision, better prepares us to refine our prediction (and, our asteroid redirection) infrastructure for the (many) future highly unlikely (yet, each potentially devastating) events of this type.
4.1. Hybrid methods. This brings us to our most recent work: the hybrid estimation of a potential YR4 lunar impact in 2032. In [6], we propose a novel hybrid strategy to uncertainty propagation that leverages the efficiency of a method G (generally some Kalman filter variant, like the UKF) when the uncertainty is nearly Gaussian, with the accuracy of a method NG (generally a Lagrangian or Eulerian filter, like the PF) when the uncertainty is highly non-Gaussian. The key to this hybrid approach is inferring higher-order information of the true uncertainty from the Kalman filter abstraction in the G method, switching to the NG method automatically when this information indicates divergence from Gaussianity, and subsequently switching back to the G method automatically when the information it provides indicates a return Gaussianity.
Previous hybrid filtering techniques [e.g., Raihan and Chakravorty (2018), Frei and Künsch (2013)] rely on measurement corrections to trigger this transition. However, when measurements are sparse (as is the case for YR4 and other distant objects, out of view of Earth-based telescopes) these strategies are inapplicable. Our hybrid filter, dubbed the UKF/PF filter (that is, pairing the UKF for the G method with the PF for NG method), uses a measure of skew derived from the weighted sigma points propagated by the UKF to indicate when the uncertainty in the G calculation is diverging from Gaussianity, and a measure of normality applied to the large sample propagated by the PF to indicate when the uncertainty in the NG calculation is converging back to Gaussianity.
In the remarkable video above, our hybrid UKF/PF method is demonstrated for the estimation of YR4. As YR4 approaches its first close approach with the Earth-Moon system in Dec 2028, the UKF/PF method autonomously switches from UKF mode to PF mode, effectively capturing the evolution of the non-Gaussian uncertainty in this portion of YR4’s trajectory. After the close approach, the UKF/PF method autonomously switches from PF mode back to the (more efficient) UKF mode for the quiescent, slower-moving portion of the trajectory. Later in the video, we demonstrate that the UKF/PF method accurately estimates the potential lunar impact corridor of YR4 as compared to a large (more expensive) MC, achieving a speedup of nearly 5x. This work successfully demonstrates the capability of our hybrid strategy to leverage, automatically, the respective strengths of the (efficient) G method and the (accurate) NG method that it incorporates.
Future work will replace the G and NG methods used by the proposed hybrid approach to further improve both performance and accuracy. Specifically, we plan on replacing the NG method with GBEES for the reasons discussed previously. Additionally, the G method may also be replaced by an extension of the UKF that uses a Conjugate Unscented Transform to match the third and fourth central moments of the random variable. This will allow the hybrid filter to use the more efficient G mode for longer while still maintaining accuracy.
Although the UKF/PF method in [6] was applied to the uncertainty propagation of a natural satellite, this general (hybrid) method also demonstrates considerable potential for use in spacecraft navigation. To date, Lagrangian and Eulerian filters have been underutilized for spacecraft state estimation due to their computational expense in 6-dimensional problems. Kalman filters have mostly sufficed thus far, requiring frequent measurement corrections from ground stations. However, there are many trajectories of interest (particularly, when exploring the moons around other planets in our solar system, in the continuing search for extraterrestrial life) in which state uncertainty may become highly non-Gaussian between measurement updates. In such cases, a hybrid filter allows for an automatic and accurate transition from G to NG methods only when necessary, saving significant computational resources while still accurately capturing the portions of such trajectories with non-Gaussian uncertainty.
[1] Bewley & Sharma (2012) Efficient grid-based Bayesian estimation of nonlinear low-dimensional systems with sparse non-Gaussian PDFs. Automatica 48, 1286.
[2] Hanson, Rosengren, & Bewley (2024) State Estimation of Chaotic Trajectories: A Higher-Dimensional, Grid-Based, Bayesian Approach to Uncertainty Propagation, AIAA 24-426.
[3] Hanson, Rosengren, Bewley, & Ely (2025) Non-Gaussian recursive Bayesian filtering for outer planetary orbilander navigation, AIAA 25-194.
[4] Hanson, Ely, Bewley, & Rosengren (2025) Bayesian Benchmarking of GBEES Applied to Outer Planet Orbiter Estimation, Journal of Guidance, Control, and Dynamics, DOI:10.2514/1.G009146.
[5] Hanson, Rubio, García-Gutiérrez, & Bewley (2025) GBEES-GPU: An efficient parallel GPU algorithm for high-dimensional nonlinear uncertainty propagation, Computer Physics Communications Computer Physics Communications 317, 109819.
[6] Hanson, Carton, Bewley, Ely, Rosengren (preprint) Hybrid, Ephemeris-Quality, Measurement-Free Estimation of the Potential 2024 YR4 Lunar Impact.
5. Postscript If you made it this far in this post, kudos (even if, perhaps, skimming some of the technical stuff - apologies…). As discussed here, science like this is difficult to get right (much more difficult than suggested in the sensational movies mentioned previously). It takes attention to detail and years of focused development. When you do get it right, its real-world impact can be substantial, from successfully navigating the moons of the outer planets of our solar system in order to better understand physics and, possibly, the origins of life itself, to predicting and possibly averting major life-threatening asteroid collisions here on Earth.
Science may, in the end, be the very thing that saves humanity from another major asteroid collision, or perhaps even from ourselves.
